Nicolas Oudart

Gommers, R P Virtanen

Peterson, P., Weckesser, W E Burovski

Bright, J S J Van Der Walt

M Brett

PyBWE: Python tools for Bandwidth Extrapolation of planetary radar signals

PyBWE is a Python library containing radar super-resolution methods known as "Bandwidth Extrapolation" (BWE).

Range resolution enhancement is one of the main challenges in radar signal processing. It is driven by the time resolution of radar soundings, and the speed of electromagnetic waves in the sounded material. Time resolution being limited by the frequency bandwidth of the instrument, the same applies to range resolution: the larger the bandwidth, the better the resolution.

. Based on a signal model, assuming deterministic properties of the signal, the output of such techniques yields a better resolution. However, parametric techniques are less robust than classic Fourier transform ones, in particular in the presence of noise or distortions.

a State-Space modelling approach.

47

Each package contains an integrated solution to directly apply the complete BWE process to 48 a given radar spectrum, as well as all the individual functions for modelling and extrapolation. 49

The Bandwidth Extrapolation technique (BWE) is a compromise between a classic Fourier transform and a parametric spectral estimation technique [START_REF] D R A F T Cuomo | A bandwidth extrapolation technique for improved range resolution of coherent radar data[END_REF]. A parametric model is fitted to the measured signal's spectrum, this model is then used to extrapolate this spectrum forward and backward, and the spectrum is eventually Fourier transformed using an IFFT. The extrapolation factor is equal to the resolution enhancement, and can be up to 3 in practical cases.

The extrapolation of a radar spectrum is indeed possible, as it can be modelled by a sum of complex sine-waves with different amplitudes, frequencies and phases. Each of these complex sine-waves corresponds to a target echo in time-domain. In practical cases, this deterministic signal is corrupted by various sources of noise, distorsions, and can be dampened by losses in the sounded material. Hence the limitation of a radar spectrum extrapolation.

An example of application on a synthetic planetary radar spectrum, inspired by the WISDOM (ExoMars, ESA) planetary radar (Ciarletti et al., 2017), is shown in Figure 1.

In the regular BWE, the signal is modelled by an autoregressive (AR) model, using the Burg algorithm. Several improvements to the BWE have been proposed: the Polarimetric BWE (PBWE) using the correlation between several polarimetric radar channels for an improved extrapolation (K. [START_REF] Suwa | A bandwidth extrapolation technique of polarimetric radar data and a recursive method of polarimetric linear prediction coefficient estimation[END_REF][START_REF] Suwa | A two-dimensional bandwidth extrapolation technique for polarimetric synthetic aperture radar images[END_REF], or the State-Space BWE (SSBWE) using a State-Space model accounting for noise and exponential distortions in the signal [START_REF] Piou | A state-space technique for ultrawide-bandwidth coherent processing[END_REF]. Also, the BWE can be used to fill a gap between two spectra of multiband radars [START_REF] Moore | Enhanced imagery using spectral-estimation-based techniques[END_REF]. This process is known as Bandwidth Interpolation (BWI).

Figure 1 :

 1 Figure 1:An example of BWE application to a synthetic planetary radar spectrum, with 2 targets in free-space separated by a distance of 7 cm, below the instrument's resolution limit with IFFT and Hamming windowing. This spectrum is generated with 2 complex sine-waves of amplitude 1, corrupted by a white-noise of standard deviation 0.1 (SNR = 20 dB), for frequencies between 0.5 and 3 GHz. (a) Extrapolation of this spectrum with a parametric model. (b) Corresponding time-domain signal obtained with classic FFT or BWE.

Oudart. (2024). PyBWE: Python tools for Bandwidth Extrapolation of planetary radar signals. Journal of Open Source Software, 0(0), 6319. https://doi.org/10.xxxxxx/draft.

Acknowledgments

81

The author would like to thank Valérie Ciarletti, Alice Le Gall and Emile Brighi for their helpful 82 comments and suggestions.


David Dalmazzo

Ken Déguernel
email: ken.deguernel@cnrs.fr

Bob L T Sturm
email: bobs@kth.se*

The Chordinator: Modeling Music Harmony By Implementing Transformer Networks and Token Strategies ⋆

Keywords: Chord progressions, Transformer Neural Networks, Music Generation

This paper compares two tokenization strategies for modeling chord progressions using the encoder transformer architecture trained with a large dataset of chord progressions in a variety of styles. The first strategy includes a tokenization method treating all different chords as unique elements, which results in a vocabulary of 5202 independent tokens. The second strategy expresses the chords as a dynamic tuple describing root, nature (e.g., major, minor, diminished, etc.), and extensions (e.g., additions or alterations), producing a specific vocabulary of 59 tokens related to chords and 75 tokens for style, bars, form, and format. In the second approach, MIDI embeddings are added into the positional embedding layer of the transformer architecture, with an array of eight values related to the notes forming the chords. We propose a trigram analysis addition to the dataset to compare the generated chord progressions with the training dataset, which reveals common progressions and the extent to which a sequence is duplicated. We analyze progressions generated by the models comparing HITS@k metrics and human evaluation of 10 participants, rating the plausibility of the progressions as potential music compositions from a musical perspective. The second model reported lower validation loss, better metrics, and more musical consistency in the suggested progressions.

Introduction

The domain of music computing has seen a variety of applications of state-ofthe-art Natural Language Processing (NLP) models or neural networks, such as music generated from symbolic music notation [START_REF] Huang | Music transformer: Generating music with long-term structure[END_REF][START_REF] Briot | Deep Learning Techniques for Music Generation[END_REF], sound synthesis [START_REF] Van Den Oord | WaveNet: A Generative Model for Raw Audio[END_REF][START_REF] Caillon | Rave: A variational autoencoder for fast and high-quality neural audio synthesis[END_REF], or text-to-sound models such as MusicLM [START_REF] Agostinelli | Musiclm: Generating music from text[END_REF]. Nevertheless, compared to melodic sequences, the modeling of harmonic progressions has seen comparatively little research. Hence, we are interested in modeling chords from symbolic music notation by presenting an approach to encode chord progressions (CP). We designed a system to suggest a possible next chord after a CP, a complete CP for a specific bar length, or a parallel plausible CP, where the human user is the leading composer.

When encoding harmony information and chord vocabulary, there is no unified method in the literature. In Western modern popular music, the basic cell of the chord is the triad as a stack of thirds, containing the scale's 1st, 3rd, and 5th notes. The triad has four possible natures: major, minor, augmented, or diminished. A chord can also be expressed in tetrads, adding the 7th scale degree (e.g., Cmaj7: C, E, G, B), having a richer range of possible natures (e.g., major 7th, minor 7th, half-diminished, diminished, dominant, minor-major, among others). Non-triad-based chords are also in the vocabulary, for instance, a stack of fourths or fifths or suspended chords where the 3rd is substituted by the 2nd or 4th. Chords can also be extended further by stacking even more thirds, adding scale degrees 9th, 11th, or 13th (e.g. Dmaj7♯11: D, F♯, A, C♯, G). Moreover, chords can have alterations that modify one or several notes and subtractions that eliminate one or several notes. If we consider chord inversions, the range of vocabulary is considerably extended. Hence, modeling symbolic music notation by implementing current NLP architectures to encode form, structure, and harmonic information without losing vocabulary richness is one of the current challenges in the field.

The following three points give the contribution of the proposed system:

• A dataset with 70,812 CP from several popular music styles.

• A tokenization method to encode chord symbols without reducing its vocabulary. • A positional encoding strategy with domain-specific embeddings of the transformer network to better suit spatial music information.

The following section reviews work applying machine learning to modeling and generating chord progressions. Section 3 presents our application of the transformer architecture to this problem. Section 4 presents the results of our models and analyzes them from a few different perspectives. We conclude with a look at the future development of our system and its integration into a pipeline for music creation. The supplementary material and related code to the publication can be found on GitHub. 1 notation. ChordRipple [START_REF] Huang | Chordripple: Recommending chords to help novice composers go beyond the ordinary[END_REF] presented a Chord2Vec encoder by adapting Word2Vec using datasets derived from the Bach chorales and The Rolling Stone Top 200 songs corpus. The chord vocabulary includes triads, tetrads, and inversions. The authors evaluated the application with music students in a creative exercise and found benefits, but more precise stylistic development was still needed.

A text-based machine-learning model for chord progressions is proposed by Choi et al. [START_REF] Choi | Text-based lstm networks for automatic music composition[END_REF]. The authors train two Long Short-Term Memory (LSTM) networks using textual representations of chord progression, e.g., "C:7 E:min". One LSTM has a vocabulary of 39 symbols corresponding to the characters which combined create chords, e.g., "E:min" consists of five characters in this vocabulary. The other LSTM has a vocabulary of 1,259 chord symbols. In this case, "E:min" is one vocabulary element. The authors train both systems by using 2,486 scores from The Realbooks and The Fakebooks. They find both models generate progressions that are plausible within jazz styles.

Transformers [START_REF] Vaswani | Attention is all you need[END_REF] have shown an improvement over the previous neural networks, enhancing the long-term dependencies by implementing the attentional mechanism as a parallelizable architecture. Most research in the context is focused on adapting or expanding it to specific tasks.

In 2020, Wu and Yang [START_REF] Wu | The jazz transformer on the front line: Exploring the shortcomings of ai-composed music through quantitative measures[END_REF] introduced Jazz Transformer, a generative composition system including chords, melodies, and form. The model is trained using the Weimar Jazz Database [START_REF]Inside the Jazzomat -New Perspectives for Jazz Research[END_REF], which contains 456 jazz standards, including styles such as Swing, Bebop, Cool, Hard Bob, and Fusion. The data format is MIDI. The authors suggested a set of objective measurements included in the MusDr2 framework that reveals the deficiencies of machine-generated music, such as unpredictable pitch class usage, inconsistent grooving patterns, and chord progressions, or the lack of recurring structures. These metrics demonstrate the limitations of the transformer architecture as a music generator and provide practical quantitative standards for evaluating the efficacy of future automated music composition efforts. Currently, model assessment largely depends on human assessment, thus concluding that besides the training loss and evaluation rating numbers, the music generated still needs to be revised by experts.

The transformer, by default, is not fully suitable as a neural network to encode chord progressions by only providing symbolic music notation as it misses the configuration of voicings and chord construction; thus, in 2021, different publications suggest updating the positional embedding layer of the transformer architecture in order to provide the network with more precise information. Chen et al. [START_REF] Chen | Attend to chords: Improving harmonic analysis of symbolic music using transformer-based models[END_REF] propose an encoder-based transformer model to classify chords and provide information on harmonic progressions in the classical music domain. The dataset used is the BPS-FH dataset and the Bach Preludes. They improve the capacity of transformers to encode harmonic features by providing extra information by adding a relative positional encoding, adding a layer of contextual information with the relative function of the chord in the chord sequence.

MusicBERT, by Zeng et al. [START_REF] Zeng | Musicbert: Symbolic music understanding with large-scale pre-training[END_REF] is a symbolic music understanding model based on the BERT: Pre-training of Deep Bidirectional Transformers for Lan-guage Understanding architecture, to provide melody completion, accompaniment suggestion, and style classification and recall. They use the Million MIDI Dataset (MMD), which contains 436.631 samples3 , from different styles such as Rock, Rap, Latin, Classic, or Jazz. The data format as a symbolic music notation is shaped as bar-level masking using MIDI. The authors propose adding musical note embeddings to the positional embeddings, including time signature, tempo, bar, position, instrument, pitch, duration, and velocity.

MrBert, by Li and Sung [START_REF] Li | Mrbert: Pre-training of melody and rhythm for automatic music generation[END_REF] is an automatic music generator that creates melodies and rhythm and then creates chord progressions based on the melody. The masked language model is trained with the OpenEWLD dataset (Mu-sicXML), which is translated to music events using the Python library music21.

The system is based on parallel transformer encoders: one encoding melodies, and the other its rhythmic pattern. Authors argue that chords are generated after the melodies, using a Seq2Seq generation that changes chords when a new symbol outside the chord pattern appears. HITS@k metrics are used to evaluate the compositions. The chord vocabulary is restricted to triads.

Li and Sung [START_REF] Li | Transformer-based seq2seq model for chord progression generation[END_REF] propose an encoder-decoder (Seq2Seq) architecture to generate chord progressions given a melody. The data format is music XML. The authors compare three architectures: a) the bidirectional long short-term memory (BLSTM), b) the bidirectional encoder representation from transformers (BERT), and c) the generative pre-trained transformer (GPT2). They evaluated the models using HITS@k [START_REF] Ali | Bringing light into the dark: A large-scale evaluation of knowledge graph embedding models under a unified framework[END_REF]. The proposed methodology uses a pre-trained encoder that takes song melodies and couples their relations with the chord progression. The decoder receives the same melodic material and produces suggested chord matchings. The authors argue that the model does not need music theory inference as it relies on the patterns reflected in the symbolic music notation; however, it is not proven the model learns harmonic principles. The transformer is trained using the OpenEWLD dataset and Enhanced Wikifonia Leadsheet Dataset (EWLD), which contain 502 samples in musicXML format, including melodies and chords. 382 chord types were extracted as a vocabulary. Examples of generated music or chord progressions are not available.

The Chordinator

We compare two strategies to model chord progressions, maintaining the full range of vocabulary available in the music corpus. All approaches are based on the transformer encoder, changing the token strategy. We also propose complementary embeddings to the transformer to encode chord information by adding a new level of relative embeddings using an array of defined MIDI values added to the positional and tokens embeddings layer. In the second tokenization strategy, we include a style token reference, providing the transformer with contextual information about the style and form.

Database

The database was formed by exporting 4,300 CP included in the iReal Pro application (as songs). The data provides metadata (composer, song name, key, etc.) and information regarding the form, tempo, and chords. The database includes styles such as Jazz, Samba, Ballad, etc. (see the full list in Figure 1). We filter the CP, keeping only those in 4/4 time signature (around 90% of the dataset), to keep the sequences consistent. We apply data augmentation by transposing all CP into the 18 enharmonic keys, creating 70,812 in total. Furthermore, we prepare the data stream to train a GPT-2 4 transformer architecture with a fixed-length sequence of 512-time steps in the first model and 1024 in the second model. The < pad > token is used to fulfill the defined length.

Model Strategies

Our proposal is based on two main strategies to train a GPT-2-based transformer encoder architecture with chord progressions.

• Model 1: It has the entire token vocabulary. All chords are unique elements, producing a vocabulary of size 5,202 unique chords/tokens. The chords are quantized in time, creating a chord per quarter note. When the suggestions are generated, four consecutive chords are equivalent to an entire 4/4 bar. • Model 2: We split the chord units into subsections, simplifying the tokenization range without reducing the vocabulary in the corpus. Starting the chord element with an identifier token '.' (dot), followed by the root, nature, and extensions as shown in Table 1. For instance, 'Dm7 add 9', is split into '.', 'D', 'm7', 'add 9'. This creates a 1D array of dynamic lengths depending on how much information the chord has. Slash chords are also split, following the same structure. For example, 'Cmaj7 / Bb' results in a sequence '.', 'C', 'maj7', '/', 'Bb'. The second root in slash chords becomes the actual root reference. This strategy reduces the vocabulary to a more specific format with 59 Tokens: Song formatting (¡start¿, ¡and¿, ¡pad¿, '.'), roots ('A' ''A♯' 'Ab' 'B' 'Bb' 'C' 'C♯' 'Cb' 'D' 'D♯' 'Db' 'E' 'Eb' 'F' 'F♯' 'G' 'G♯' 'Gb'), natures ('dim', 'm', 'm6', 'maj7', 'o7', 'sus', 'power'), extensions (additions, subtractions, and alterations), and sharp ('/'). In the cases where the chord is notated with only one note, that chord is considered a major triad. • Extended Model 2: We extend Model 2 in two ways. First, we add a context token at the beginning of the samples with a style specification, meaning we start with a token (¡style¿) followed by the actual style (such as jazz, pop, blues, samba, etc.). Second, we update the transformer model, adding a MIDI Embedding relative position array into the positional embedding layer. Based on the second strategy, when we have chords divided into root, nature, and extension, the MIDI array is composed of an array of eight numbers where the first is the MIDI root reference, the next three values are used to describe the nature, and extensions could use the next value. The longest sequence of the MIDI array for a chord is seven MIDI values; hence, the eighth MIDI value in the array is an identifier (number 127) to express that there is a slash token.

MIDI Embeddings

MIDI Embeddings (ME) are extracted from the symbolic chord representation using the Python libraries music21 (symbol to notes translation) and librosa (notes to MIDI translation). We decided to maintain the chord pitches in only a two-octave range. During training and generation, ME is paired with token inputs, providing only information about the current time-step in the sequence; therefore, when only one note is presented in a specific chord, by default, the ME reads it as a Major chord (triad). When a second token appears next to the original note, the first token is read as a root note (only one MIDI note is then placed), and the nature is updated, adding its related notes to the MIDI array.

The same logic is applied to extension and sharp chord, where each step in the sequence updates the ME information, having the complete chord MIDI array only in the last time step of the chord tuple, as shown in Figure 2. We avoid information leaking following this format. The formula to update the ME into the positional embedding in the transformer network is defined as follows:

midi emb = Linear(Emb(midi vocab, dm)(m).view(m.size(0), m.size(1), -1), n embd )

Positional embeddings are added to the input embeddings before the transformer encoder is fed. The positional embedding layer enables the model to distribute information spatially in order to learn specific patterns. The ME takes information from MIDI and adjusts it into an embedding tensor to align with the required architecture format. This, then, is summed with the standard positional encoding layer. It is composed as follows:

1. Embedding MIDI Values:

• 'Emb(midi vocab, dm)(m)' : This defines an embedding layer where 'midi vocab' is the total number of unique MIDI events and 'dm' is the dimension of the embeddings. '(m)' is the time-step in the chord tuple sequence. • '.view(m.size(0), m.size(1), -1)' : The tensor matrix is reshaped to adjust its format to the standard positional embedding. Here, 'm.size(0)' represents the batch size, while 'm.size(1)' represents the sequence length. '-1' means that the size in the last dimension is computed so that the total size remains constant. 2. Linear Transformation: The reshaped data is passed through a Linear layer. The size of the output data is adapted to match a tensor of shape (batch size, sequence length, n embd ). It is then summed with the positional encoding layer.

Training The Models

The models are trained on 4 GPUS NVIDIA GeForce RTX 3090 (24576MiB). While we tested many configurations, this manuscript only reports the best versions of both strategies, defined as Model 1 and Model 2.

Model 1 is trained with 120 Epochs, 6 attention heads, 6 layers, 192 embeddings sizes, 64 batch sizes, 192 embeddings, number of workers 6, and a learning rate of 3e-5. Model 1 is prone to overfit after 120 epochs; hence, these configurations are chosen to maintain a small network that adapts to the size of the vocabulary and dataset. Extra attention heads tend to increase the overfitting tendency and do not necessarily lead to better performance.

Model 2 is trained with 250 epochs without showing overfiting tendencies (it could be trained with more epochs). Thus, it has an increment in 8 attention

Metrics

HITS@k (k = 1, 3, and 5) [START_REF] Zeng | Musicbert: Symbolic music understanding with large-scale pre-training[END_REF] is used to evaluate the predictions in both models. HITS@k calculates the proportion of the correct answer given by the candidates, denoted by the letter k. It was calculated as follows:

HIT S@k = 1 n n i=1 I(rank i ≤ k) (1)
where,

• HIT S@k: This represents the HITS at k metric.

• n: Total number of instances or items being evaluated.

• n i=1 : This denotes the sum over all instances. • I(rank i ≤ k): An indicator function that evaluates to 1 if the rank of instance i is less than or equal to k, and 0 otherwise.

Generated Chord Progressions

We generated 20 CP per model in one unique iteration loop where no selection of the CP was performed, and they were shown to the evaluators as they were; no voicing, arrangement, or editing was performed. All CP were constrained with a length of 8 to 24 chords, producing a collection of 40 randomized CP. The Model behind the generated CPs are not identifiable by the evaluators. Ten participants, with an average age of 36 (s.d. +/-8), reported using a 1 to 10 scale the plausibility of the chord progressions as a starting point of a conceivable music piece: 1 implies no plausibility as a sequence, and 10 is very plausible to start a potential composition. The overall median score received was Model 1: 5.7; Model 2: 6.7, as shown in Figure 3. All suggestions begin with a root token, which can be any scale note; nevertheless, Model 2 includes a random number defining the style from the first six most common styles in the dataset (Jazz, Samba, Ballad, Pop, Bossa, Rock). To generate a sound sample, we implemented the Python library librosa and music21 to generate the chords in MIDI format using a standard Piano. Figure 4 shows examples of generated CP with form tokens in Model 2.

Trigrams

A supplementary dataset comprises three chord cells, including information on the song and composer. Those cells are called trigrams, which provide information on the suggested progressions to check for original versions and sources.

Trigrams are formed by creating groups of subsequent three chords, moving the windowed observation one time step forward. A sequence with Em7, A7, Dmaj7, Ab7 add ♯11, Gm7, are then packed as (Em7, A7, Dmaj7), (A7, Dmaj7, Ab7 add ♯11) and (Dmaj7, Ab7 add ♯11, Gm7). That means four chords are equal if two consecutive trigrams from a suggested progression are found in a particular song. Hence, we can define a threshold of matching trigrams to determine duplication. The trigram analysis can be executed with a list of suggested CP or immediately after a generation. In this case, analyzing 40 CP, the report is too extensive to include in the manuscript. A list of the most 20 trigrams cells found is shown in Figure 8. The report includes a tag for Sequence, which is the CP analyzed; Location, which reports the chord position in the sequence; Trigram, which is the progressions cell; and a list of Composer and Song where it was found. Both Models show similar trigram repetitions found in the dataset. However, it does not mean the whole progression is copied; it only provides information about the original trigram cell location. However, it is possible to determine if a whole progression is, in fact, replicated from the training data. No generated CP reported more than two trigram cells from the same song. In Figure 7, some of the most common trigrams found are shown, and Model 1 is more prone to duplicate cells. Based on trigrams analysis, we observed that sequences based on ii-V-I progressions and their different versions in major and minor are prominent in the dataset (see Figure 5). However, the ii-V-I is not the most used progression, as shown in Figure 6. The most used progressions are combinations of IV-V-I motions and some permutations, commonly audible in the Blues, Rock, Pop, or Cinematic Music.

Discussion and Conclusion

We have presented two tokenization strategies with an extension of the positional embeddings (MIDI Embeddings) for Model 2. We treat the chord vocabulary with a specific methodology differentiating the root, nature, and extensions, thereby reducing tokens without diminishing the chord vocabulary and its richness. We trained each model on a large database of chord progressions and used it to generate new suggested progressions or chord completion. Our database of 70,812 songs has diverse music styles in popular music. Our evaluation of sequences generated by these models shows that reducing the vocabulary and adding contextual information improves the encoding process of the chord vocabulary. Thus, adding MIDI notes for chord constructions to the Positional Embeddings layer has helped the transformer network stabilize the learning process, allowing more training epochs without overfitting; it also helped to report better HITS@K and validation metrics (see Figure 3.

Our system contributes to a growing literature on modeling harmonic progressions as a system that can be used as a composer assistant. Future work will explore the usefulness of this system with users performing musical tasks to address further development and usability. We will also explore ways of generating melodies from chord progressions and vice versa. This will be implemented as an accessible web application for others to study.

As expressed in [START_REF] Wu | The jazz transformer on the front line: Exploring the shortcomings of ai-composed music through quantitative measures[END_REF], and also based on the outputs generated by Model 1, the vanilla encoder transformer architecture (GPT-2) is not an ideal network for encoding music information from only providing symbolic music notation. Even though it can generate plausible textual tokens, the musicality is still unclear. From the perspective that chord construction and voicing formats are hidden from symbolic music notation and that decision is the responsibility of the musician's knowledge, it is still necessary to propose new versions of related neural networks and test architectural modifications to find better strategies to encodes all musical information related to chord progressions. Therefore, our next model is based on multi-hot encoders (2 octaves) exemplifying the notes in the piano from an expert performer playing chord progressions and an encoder-decoder architecture to suggest, on the one hand, the symbolic music notation but, on the other, the voicings for more musical outputs.

The human evaluation reported a median score of 5.7 for Model 1 and 6.7 for Model 2. Even though there is a difference in appreciation of Model 2 quality, the low score could be directly affected by the test format shown to the participants. Raw MIDI piano sounds were chosen to evaluate the chord progressions, and no voicing was added to avoid human completion from expert performers. The CP were played at 120BPM; in that term, the suggested CPs are expressed in a raw MIDI version that can be highly improved by musicians' knowledge in terms of performance and interpretation and also in sound quality. However, HIT@K metrics reported a difference where Model 2 is much more accurate, particularly in the validation loss. The interpretation of its musicality based on NLP metrics is still open for discussion.

As a further development, to enhance Model 2, when translating symbolic music notation into MIDI chords, we will add an array of plausible voicings repertoire per chord and a distance calculation of note proximities to suggest more natural voicing. By default, the Python library music21 is unsuitable for this purpose as it reported errors, particularly with slash chord formatting and notes related to suspended chords.

Model 2 reported a tendency over Model 1 to generate more consistent progression in musical terms. The explanation is given by the modeling of the contextual token, which provides a general view of the music style; therefore, Model 1 mixes styles and formats, while Model 2 is more prone to maintain some degree of consistency. For instance, when Blues is given, the form and chords are common to the style. Also, Model 1 chooses chord tonalities wrongly formatted in terms of notation; in other words, it mixes enharmonic, e.g., when Ebmaj7 makes sense in the context, it suggests D♯maj7. Further work is also needed to adjust this duality in names as it is relevant to be fixed in a musical context.

Model 2 also suggests form and structures; further development will incorporate that suggestion into the web application.

For future work, we will add a melody generator that is conditioned on a CP and vice versa, where a melody conditions the generation of a CP. A web application as an interface to play with the model and generate suggestions is under development; it will open its usability to a broader audience. We will perform further human evaluations by expert musicians to test the impact and usability of the system.

Fig. 1 .

 1 Fig. 1. List of musical styles in the iReal Pro dataset. Values are presented in the log10 ratio.

Fig. 3 .

 3 Fig. 3. Boxplot of human evaluation of Model 1 (orange) and Model 2 (blue).

[

  '<start>' '<style>' 'Soul' '|' 'Form_A' '.' 'C' 'm' '.' 'C' 'm' 'add #7' '|' '.' 'C' 'm7' '.' 'C' 'm6' '|' '.' 'F' 'm' '.' 'F' 'm' 'add #7' '|' '.' 'F' 'm7' '.' 'E' 'm7' 'alter b5' '.' 'A' '7' '|' '.' 'D' '7sus' '.' 'D' '7' '|' '.' 'G' '7' '|' '.' 'C' 'm' '.' 'C' 'm' 'add #7' '|' '.' 'C' 'm7' '.' 'C' 'm6' '|' '.' 'F' 'm' '.' 'E' 'm7' 'alter b5' '.' 'A' '7' '|' '.' 'D' '7' '|' '.' 'G' '7' '|:' '.' 'C' 'm7'] ['<start>' '<style>' 'Pop' '|' 'Form_A' '.' 'C' 'm7' '|' '.' 'Db' 'maj7' '|' '.' 'C' 'm7' '|' '.' 'D' 'sus' '.' 'D' '|' '.' 'D' 'sus' '.' 'D' '|' '.' 'C' 'm7' '|' '.' 'Db' 'maj7' '|' '.' 'C' 'm7' '|' '.' 'D' 'sus' '.' 'D' '|' '.' '|' '.' 'Db' 'maj7' '.' 'Db' '7' '|' '.' 'C' 'm7' '|' 'F#' 'maj7' '#11' '|' '.' 'F' 'm7' '|'] Tokens Array

Fig. 4 .

 4 Fig. 4. Two samples from Model 2 using Pop and Soul context token. The tokens array provides information about chords and form, using the token '|' as a bar identification.

Fig. 5 .

 5 Fig. 5. Most used trigrams found in the dataset. Most progressions are sequences of II-V-I in major and minor functional formation.

Fig. 6 .

 6 Fig. 6. Most common harmonic functions found in the dataset.

Fig. 7 .

 7 Fig. 7. Model 1 (orange) and Model 2 (blue) report trigrams cells also found in the dataset. However, Model 1 is more prone to replicate common trigram cells found in the dataset

Fig. 8 .

 8 Fig. 8. Trigram analysis provides information about trigram cells found in the dataset. It contains the sequence number; we have 40 CP suggested in this case. The Sequence tab will point to that suggested CP origin. The Location tab provides information on the chord position within the progression; the Trigram tag is the chord cell, and the Count tag expresses how many times it was found in the dataset. The Figure contains 20 of the most trigram cells found

Table 1 .

 1 Chord Tokens Format

	Chord	Start Root Nature Ext. Slash Root 2
	C	.	C
	Dm7 add 9	.	D	m7	add 9
	Bbmaj7 ♯11	.	Bb	maj7	♯11
	C♯o7	.	C♯	o7
	Ebmaj7/D	.	Eb	maj7	/	D

  and 8 layers. The other features are maintained as they are: 192 embeddings sizes, 64 batch sizes, 192 embeddings, number of workers 8, a learning rate of 3e-5, and a MIDI vocabulary size of 128, which refers to the size of the possible MIDI values. (see table 2).

		Pred.	.	Bb	7	add #11
					Positional Encoding
	Input	.		D	m7	add 9
	Midi Emb.	M.	MD		Mm7	Madd9
	Pos. Emb.	P.	PD		Pm7	Padd9
	Token Emb.	W.	WD		Wm7	Wadd9

MIDI

Fig.

2

. MIDI Embeddings: Formatted as an array of eight values for each time-step of the chord tuple. In Model 2, the chord information is split into tuples, and the input layer processes it as a dynamic sequence. The input sequence is translated into tokens embeddings, positional embeddings, and midi embeddings; all those layers share the same tensor matrix shape, and they are summed inside the positional encoding layer, which is passed to the multi-head attention. heads

Table 2 .

 2 Transformer Format

	Configuration	M1: Full Tokens	M2: MIDI Emb
	Epochs	120	250
	Heads	6	8
	Layers	6	8
	Embedding	192	192
	Batch size	64	64
	Learning rate	3e-5	3e-5
	Workers	8	8
	MIDI vocab	-	128

Table 3

 3

compares the metrics of Model 1 and 2.

Table 3 .

 3 Metrics

		Train Loss Val. Loss HITS@3 HITS@5
	Model 1	0.3699	0.4108	0.7383	0.8294	0.8556
	Model 2	0.04982	0.03645	0.9138	0.9707	0.984
	4 Results

https://github.com/slSeanWU/MusDr

https://github.com/jeffreyjohnens/MetaMIDIDataset

https://github.com/karpathy/minGPT

⋆ This paper is an outcome of: MUSAiC, a project that has received funding from the European Research Council under the European Union's Horizon 2020 research and innovation program (Grant agreement No. 864189); and a Margarita Salas Grant, UPF, Barcelona,


Catherine Y Tu
email: catherine.tu@gmail.com

Juhyun Park
email: juhyun.park@lancaster.ac.uk

Haonan Wang
email: wanghn@stat.colostate.edu

ESTIMATION OF FUNCTIONAL SPARSITY IN NONPARAMETRIC VARYING COEFFICIENT MODELS FOR LONGITUDINAL DATA ANALYSIS

Keywords: Functional sparsity, group bridge, longitudinal data, model selection, nonparametric regression

We study the simultaneous domain selection problem for varying coefficient models as a functional regression model for longitudinal data with many covariates. The domain selection problem in a functional regression mostly appears within a functional linear regression with a scalar response; however, there is no direct correspondence to functional response models with many covariates. We reformulate the problem as a nonparametric function estimation problem under the notion of functional sparsity. Sparsity encapsulates interpretability in a regression with multiple inputs, and the problem of sparse estimation is well understood in the context of variable selection in a parametric setting. For nonparametric models, interpretability not only concerns the number of covariates involved, but also the zero regions in the functional form. Thus, the sparsity consideration is much more complex. To distinguish the types of sparsity in nonparametric models, we refer to the former as global sparsity and to the latter as local sparsity, both of which constitute functional sparsity. Most existing methods focus on directly extending the framework of parametric sparsity for linear models to nonparametric models to address one type of sparsity, but not both. We develop a penalized estimation procedure that simultaneously addresses both types of sparsity in a unified framework. We establish the asymptotic properties of estimation consistency and sparsistency of the proposed method. Our method is illustrated by means of a simulation study and real-data analysis, and is shown to outperform existing methods in terms of identifying both local and global sparsity.

Introduction

We study the simultaneous domain selection problem for varying coefficient models as a functional regression model for longitudinal data, where the response variable changes over time, recorded for multiple subjects with multiple predictors. The varying coefficient models [START_REF] Hastie | Varying-coefficient models[END_REF]; [START_REF] Hoover | Nonparametric smoothing estimates of time-varying coefficient models with longitudinal data[END_REF]) are defined as y(t) = x T (t)β(t) + (t), (1.1)

where y(t) is the response at time t, x(t) = (x 1 (t), . . . , x p (t)) T is a vector of predictors at time t, (t) is an error process independent of x(t), and β(t) = (β 1 (t), . . . , β p (t)) T is a vector of time-varying regression coefficient functions. This model assumes a linear relationship between the response and predictors at each observation time point, but allows the coefficients to vary over time, thus greatly enhancing the utility of the standard linear model formulation. For generality, we consider the predictors to be functions. However, note that the varying coefficient models are equally applicable when the predictors take scalar values.

The domain selection problem in functional regression is known to be intrinsically difficult [START_REF] Müller | Peter hall, functional data analysis and random objects[END_REF]). Most prior studies examine the problem as a functional linear regression with a scalar response and a single functional covariate. [START_REF] Hall | Truncated linear models for functional data[END_REF] formulated the problem as a truncated regression model with a single unknown domain, in order to study the identifiability issues in nonparametric function estimation. James, [START_REF] James | Functional linear regression that's interpretable[END_REF] approached the problem from the viewpoint of sparsity estimation as interpretable solutions. Using a grid approximation, they imposed parametric sparsity constraints on the derivatives of the underlying function at a large number of grid points, which produces an estimate that distinguishes between zero and non-zero regions. However, as Zhou, Wang and Wang (2013) have noted, the overlapping contribution of each coefficient to the neighboring regions means the independent shrinkage of the coefficients does not necessarily induce zero values in the coefficient function, in general. Thus, the procedure tends to over-penalize. As a remedy, Zhou, Wang and Wang (2013) suggested a two-step estimation procedure. [START_REF] Wang | Functional sparsity: Global versus local[END_REF] studied a similar problem under a standard nonparametric regression, suggesting the need to distinguish the functional features from the parametric variable selection.

We consider the regression problem under a functional response variable with varying coefficient models, involving multiple domain selection under the general setting, where the true number of covariates is unknown. Although the views and approaches adopted in prior studies are quite different, the domain selection problem can be motivated as a means to enhance the interpretability of the model selection in nonparametric models. In this regard, we share the view that considering some form of sparsity could be useful. For nonparametric models, however, interpretability not only concerns the number of covariates involved [START_REF] Wang | Variable selection in nonparametric varying-coefficient models for analysis of repeated measurements[END_REF]; [START_REF] Noh | Sparse varying coefficient models for longitudinal data[END_REF]; [START_REF] Wei | Variable selection and estimation in high-dimensional varying-coefficient models[END_REF]; [START_REF] Xue | Variable selection in high-dimensional varying-coefficient models with global optimality[END_REF]), but also the zero regions in the functional form [START_REF] James | Functional linear regression that's interpretable[END_REF]; Zhou, Wang and Wang (2013)). To distinguish between the types of sparsity in nonparametric models, we refer to the former as global sparsity and to the latter as local sparsity, both of which constitute functional sparsity [START_REF] Tu | Functional model selection for sparse binary time series with multiple input[END_REF]; [START_REF] Wang | Functional sparsity: Global versus local[END_REF]). More formally, a function has global sparsity if it is zero over the entire domain. This indicates that the corresponding covariate is irrelevant to the response variable. A function has local sparsity if it is nonzero, but remains zero for a set of intervals. Thus, this identifies an inactive period for the corresponding covariate. These notions of interpretability are used informally in a separate context of the analysis. Thus, the significance of local sparsity estimation has not been well recognized.

We reformulate the domain selection problem as a nonparametric function estimation problem under the unified theme of functional sparsity. Then, we propose a one-step penalized estimation procedure that automatically determines the type of functional sparsity (i.e., local or global). Although we distinguish between the two types of sparsity on a conceptual level, our unified formulation does not require this distinction for the implementation. We directly exploit the fact that global sparsity is a special case of local sparsity from the viewpoint of domain selection, but not the other way around. Furthermore, the consistency of the coefficient function estimation does not necessarily provide us with information on local sparsity. This feature distinguishes our approach from the majority of existing methods that target global sparsity, such as [START_REF] Wei | Variable selection and estimation in high-dimensional varying-coefficient models[END_REF]). Note that there is a fundamental difference in the underlying assumption on sparsity between parametric and nonparametric models, because we focus on function estimation with dependent variables over an unknown domain. This difference was also recognized by [START_REF] Kneip | Functional linear regression with points of impact[END_REF]. Moreover, for parametric sparsity, an underlying sparse vector is specified. In contrast, the true sparse representation for functional sparsity may not be well defined in the function approximation. These differences pose different conceptual challenges in the development. Our proposed penalized procedure resembles a type of parametric sparsity estimation. However, our analysis is not comparable with those that focus on high-dimensional parametric sparsity estimation (e.g., [START_REF] James | Functional linear regression that's interpretable[END_REF]).

We provide a theoretical analysis of the proposed method. In particular, we show that the local sparsity can be recovered consistently, and even diluted with the problem of global sparsity estimation. We study the properties of our proposed method under the standard assumptions on nonparametric smooth function estimation, and exploit the functional property in a more natural manner. In this way, we contribute to bridging the gap between parametric variable selection and nonparametric functional sparsity in a coherent manner.

Our formulation is given in Section 2. Our approach is a one-step procedure that allows us to directly control the functional sparsity through the coefficient functions themselves, rather than using a pointwise evaluation. In Section 3, we study the large-sample properties of the proposed method and establish the consistency and sparsistency of the function estimates. Section 4 describes our simulation studies under different scenarios, and a real-data analysis is provided in Section 5, demonstrating the utility of functional sparsity in relation to the interpretability of the results. All technical assumptions and proofs are provided in the online Supplementary Material.

Methodology

Suppose that, for n randomly selected subjects, observations of the kth subject are obtained at {t kl , l = 1, . . . , n k }, and the measurements satisfy the varying coefficient linear model relationship in (1.1):

y k (t kl ) = x T k (t kl )β(t kl ) + k (t kl ), (2.1)
where x k (t kl ) = (x 1 (t kl ), . . . , x p (t kl )) T and y k (t kl ) is the response of the kth subject at t kl . We assume that β i (t), for i = 1, . . . , p, are smooth coefficient functions with bounded second derivatives for t ∈ T . We use spline approximations to represent β(t) and formulate a constrained optimization problem for the parameter estimation.

Least squares estimation under a b-spline approximation

B-spline approximations are widely used to estimate smooth nonparametric functions. For a detailed discussion on B-splines, see de [START_REF] De Boor | A Practical Guide to Splines[END_REF] and [START_REF] Schumaker | Spline Functions: Basic Theory[END_REF]. Specifically, for a smooth function β(t), t ∈ [0, 1], its approximant can be written as

β(t) = J j=1 α j B j (t), (2.2)
where {B j (•), j = 1, . . . , J} is a group of B-spline basis functions of degree d ≥ 1 and knots 0 = η 0 < η 1 < • • • < η K < η K+1 = 1. Note that K is the number of interior knots and J = K + d + 1. Here, we adopt the definition of a B-spline as stated in Definition 4.12 of [START_REF] Schumaker | Spline Functions: Basic Theory[END_REF]. In general, the performance of Bspline approximations has been well studied. For instance, under mild conditions, there exists a function β(t) of the form (2.2), such that the approximation error goes to zero. See Theorem 6.27 of [START_REF] Schumaker | Spline Functions: Basic Theory[END_REF] for further detail.

We write the B-spline approximation for each smooth nonparametric coefficient function as

β i (t) = Ji j=1 α ij B ij (t) = B i (t) T α i , t ∈ [0, 1], i = 1, . . . , p, (2.3)
where B i (t) = (B i1 (t), . . . , B iJi (t)) T , α i = (α i1 , . . . , α iji ) T , and

J i = K i + d + 1.
Here, K i is the number of interior knots for β i (t), which may vary over i. For simplicity, we assume that the knots are evenly distributed over [0, 1]. Define a block diagonal matrix B(t) as

B(t) = diag{B T 1 (t), . . . , B T p (t)}. Using (2.3) in the varying coefficient model (2.1) leads to y k (t kl ) ≈ x T k (t kl )B(t kl )α + k (t kl ) = U k (t kl )α + k (t kl ) , where U k (t kl ) = x T k (t kl )B(t kl ) and α = (α T 1 , . . . , α T p ) T .
The least squares criterion of α [START_REF] Huang | Varying-coefficient models and basis function approximations for the analysis of repeated measurements[END_REF]) is defined as

(α) = n k=1 ω k y k -U k α 2 2 ,
where y k = (y k (t k1 ), . . . , y k (t knk )) T and U k = (U T k (t k1 ), . . . , U T k (t knk )) T . The weights ω k , for k = 1, . . . , n, are usually chosen as ω k ≡ 1 or ω k ≡ 1/n k [START_REF] Huang | Polynomial spline estimation and inference for varying coefficient models with longitudinal data[END_REF]). In this study, for simplicity, we set equal weights to every subject (i.e., ω k ≡ 1). Setting U = (U T 1 , . . . , U T n ) T and y = (y T 1 , . . . , y T n ) T , the least squares criterion l(α) can be written in matrix form; that is, l(α) = y -U α 2 2 . Huang, [START_REF] Huang | Polynomial spline estimation and inference for varying coefficient models with longitudinal data[END_REF] proved that, under certain assumptions, the matrix U T U is invertible for fixed p. Consequently, l(α) has a unique minimizer,

α LSE = (U T U ) -1 U T y,
which is the least squares estimator (LSE) of α. Thus, the LSEs of the coefficient functions are

β LSE i (t) = Ji j=1 α LSE ij B ij (t), i = 1, . . . , p,
where α LSE ij denotes an entry of α LSE . Here, we take a marginal approach -Figure 1. Top: A graphical display of a smooth function (solid thick line) and two approximating functions from a family of cubic B-spline basis functions with nine equally spaced interior knots. Bottom: A graphical display of the set of B-spline functions used in the approximation. [START_REF] Wu | Nonparametric Regression Methods for Longitudinal Data Analysis: Mixed-effects Modelling Approaches[END_REF]) to construct the LSE criterion, without accounting for within-subject correlation. Proper modeling of the covariance structure would require further parametric assumptions [START_REF] Diggle | Analysis of Longitudinal Data[END_REF]) or nonparametric smoothing techniques [START_REF] Wu | Nonparametric Regression Methods for Longitudinal Data Analysis: Mixed-effects Modelling Approaches[END_REF]), which are not the focus of this study.

B-spline approximation and sparsity

From B-spline approximation theory, there exists a function of the form given in (2.2) that is very close to the true underlying function. However, this function is not capable of characterizing the functional sparsity of the true function. Here, the term "functional sparsity" is a generalization of the "parameter sparsity" in regression models; see [START_REF] Wang | Functional sparsity: Global versus local[END_REF] for further detail.

For better illustration, we consider a toy example in Figure 1. Here, in the top panel, a smooth function β(t) (thick line) with two spline approximants (dashed, dotted) is depicted. In the bottom panel, a family of cubic B-spline basis functions with nine interior knots is shown. The "best" fitted function from the L 2 criterion is shown as the dashed line in the upper panel, which signifies good performance of the approximation. Furthermore, β(t) is zero on [0, 0.1] and [0.9, 1], but its approximation is not zero, except for some singletons. From this aspect, the approximation does not capture the sparsity of the true underlying function. In contrast, the dotted curve depicted in the upper panel, also a linear combination of the B-spline basis functions, automatically corrects the function to reflect local sparsity, with almost indistinguishable performance.

The other extreme case arises when the function is close to zero, for part of or the whole of the interval. Our goal is to pursue a sparse solution, up to a function approximation error, within the linear space spanned by the B-spline basis functions. From a nonparametric estimation viewpoint, such a solution preserves the statistical accuracy and enhances interpretability; in fact, it is indistinguishable from the true underlying function.

Inspired by the above observations on functional sparsity, we develop a new procedure that equips the least squares criterion with a regularization term. Usually, the regularization on parameters is expressed in terms of a penalty function. Below, we introduce a composite penalty based on the B-spline approximation of the coefficient functions.

Penalized least squares estimation with a composite penalty

It is not too difficult to see that global sparsity corresponds to the group variable selection of α i , as a whole. To achieve local sparsity, these estimates need to be adjusted so that some of the estimates can be exactly zero. As demonstrated in Section 2.2, for the B-spline approximation, when α j = 0 for j = l, . . . , l+d, the approximation β(t) = 0 on the interval [η l-1 , η l ). In particular, when α j = 0 for all j, β(t) = 0 over the entire domain of [0, M ]. This suggests local sparsity needs to be imposed at the level of a group of neighboring coefficients. To incorporate global sparsity in the varying coefficient model, we require another layer for the group structure. These considerations lead to a composite penalty, defined as follows:

L γ 1 (α) = p i=1 Ki+1 m=1   m+d j=m |α ij |   γ ,
which can be written simply as

L γ 1 (α) = p i=1 Gi g=1 α Aig γ 1 , (2.4)
where α Aig = (α ig , . . . , α i(g+d) ) , for i = 1, . . . , p, g = 1, . . . , G i . The number of groups for the ith coefficient function is

G i = K i + 1.
Equipping the least squares criterion with the penalty defined in (2.4), we obtain the following penalized least squares (PLS) criterion:

pl(α) = y -U α 2 2 + λ p i=1 Gi g=1 α Aig γ 1 , (2.5)
where λ > 0 and 0 < γ < 1 are tuning parameters. The proposed penalized LSE (PLSE) α = α(λ, γ) is defined as the minimizer of pl(α). Consequently, the functional estimate of β i (t) is given by

β i (t) = B i (t) T α i , where α i is a subvector of α.
Note that, for γ ∈ (0, 1), the penalized criterion pl(α) is not a convex function of α. Thus, we implement the following iterative algorithm, proposed by [START_REF] Huang | A group bridge approach for variable selection[END_REF], to minimize (2.5):

Step 1. Obtain an initial value α (0) .

Step 2. For a given tuning parameter λ n , and for l = 1, 2, . . . , compute

θ (l) ig = 1 -γ τ n γ γ α (l-1) Aig γ 1 , for i = 1, . . . , p, g = 1, . . . , G i , where τ n = (λ n ) 1/(1-γ) γ γ/(1-γ) (1 -γ).
Step 3. Compute

α (l) = arg min α y -U α 2 2 + p i=1 Gi g=1 (θ (l) ig ) 1-1/γ α Aig 1 .
Step 4. Repeat steps 2 and 3 until convergence.

Note that, unlike the standard LASSO method, step 3 requires an overlapping LASSO. Because the grouping does not change at each iteration, this can be solved easily using a simple linear transformation with a grouping indicator matrix for α.

The motivation for this algorithm is a reparametrization of the nonconvex optimization problem into a complex optimization problem in terms of (θ, τ ), which reaches an equivalent solution. In essence, the suggested algorithm performs an iteratively reweighted LASSO until convergence. Thus, steps 2 and 3 can be expressed in more compact form. Given (λ n , γ),

Step 1. Obtain an initial value α (0) .

Step 2. For l = 1, 2, . . . , define ν

(l) ig = γ α (l-1) Aig γ-1
1 , for i = 1, . . . , p; g = 1, . . . , p.

Step 3. Solve

α (l) = arg min α y -U α 2 2 + λ n p i=1 Gi g=1 ν (l) ig α Aig 1 .
Step 4. Repeat steps 2 and 3 until convergence.

Variance estimation

In this section, we consider the problem of finding the asymptotic variance of our proposed estimator of the coefficient functions. Let α S denote the nonzero estimators of the coefficients α ij . Then from Step 3 in the aforementioned algorithm and the Karush-Kuhn-Tucker condition, we have

α S = U T S U S + 1 2 Θ S -1 U T S y,
where U S is a sub-matrix of U , with each column corresponding to the selected α ij , and

Θ S is a diagonal matrix, diag    g:Aig j θ 1-1/γ ig | α ij | , for α ij = 0    .
In the absence of covariance modeling of y, we further approximate the variance of y by σ 2 I, where σ 2 can be estimated by σ 2 = y -U α 2 2 /n. Thus, similar to [START_REF] Wang | Variable selection in nonparametric varying-coefficient models for analysis of repeated measurements[END_REF], the asymptotic variance of α S may be expressed as

avar( α S ) = U T S U S + 1 2 Θ S -1 U T S U S U T S U S + 1 2 Θ S -1 σ 2 .
Let B i (t) be the i-th row of the basis matrix B(t). Thus, the functional estimate of β i (t) can be written as

β i (t) = B i (t) α. Correspondingly, the asymptotic variance of β i (t) is avar( β i (t)) = B iS (t)avar( α S )B T iS (t), (2.6)
where B iS (t) is a sub-vector of B i (t), with each element corresponding to the selected α ij . Note that the estimator of α depends on the choice of λ; thus, asymptotic variances of α S and β i (t) are also tuning-parameter dependent. Although this is a naive estimator, as shown in our numerical studies, its approximation is nevertheless found to be effective in capturing the level of variability. Alternatively, we can estimate the full covariance function nonparametrically. However, owing to the additional complexity in the implementation with irregular design points, this is not very practical. The literature takes a more pragmatic approach of using a random-effects formulation (e.g., [START_REF] Wu | Nonparametric Regression Methods for Longitudinal Data Analysis: Mixed-effects Modelling Approaches[END_REF]). However, the difficulty of selecting the covariates in the random-effects terms under the current context of sparse function estimation outweighs the potential benefits; thus, we do not pursue this. Instead, we investigated using a fully nonparametric approach to estimate the covariance surface by means of a functional principal component analysis [START_REF] Yao | Functional data analysis for sparse longitudinal data[END_REF]). However, our numerical study did not identify a clear advantage to this approach. Further investigation is left for future work.

Choice of tuning parameters

To fit the model using a finite sample, we consider how to calibrate the tuning parameters. The tuning parameter λ > 0 balances the trade-off between the goodness-of-fit and the model complexity. When λ is large, we have strong penalization, and thus are more likely to obtain a sparse solution with poor model fitting. With a small λ, we would select more variables and obtain better estimation results, but lose control of the functional sparsity. In the classical nonparametric approaches, criteria such as the AIC, BIC, and GCV [START_REF] Wahba | Spline Models for Observational Data[END_REF]) are commonly used for model selection. Previous analyses noted that the AIC and GCV tend to select more variables, and thus are better suited to predictions. We use a BIC-type criterion in our analysis, reported in Section 4. To account for the increasing number of parameters when comparing models with varying dimensions, we use the extended BIC (EBIC) [START_REF] Huang | Variable selection in nonparametric additive models[END_REF]), which also penalizes the size of the full model. The EBIC is given by

EBIC(λ) = log y -U α(λ) 2 2 N + K(λ) log(N ) N + νK(λ) log( p i=1 J i ) N ,
where N = n k=1 n k , α(λ) is the penalized estimator of α, given λ, and K(λ) is the total number of nonzero estimates in α(λ). Then, p i=1 J i = p i=1 (K i +d+1) is the total number of parameters in the full model. Note that when ν = 0, the EBIC is the same as the BIC, but when ν > 0, the EBIC imposes greater penalty on overfitting. We use ν = 0.5, as suggested in [START_REF] Huang | Variable selection in nonparametric additive models[END_REF]).

Note that the tuning parameter γ influences the performance of the group selection. A value of γ that is too small or too large could lead to inefficient group variable selection. When γ is close to 1, (2.4) is close to the L 1 penalty. Consequently, the minimizer of (2.5) may not achieve functional sparsity in its solution. Unlike λ, however, 0 < γ < 1 is more often viewed as a higher-level model parameter (often set equal to 0.5, [START_REF] Huang | A group bridge approach for variable selection[END_REF])), similarly to how the LASSO (γ = 1) estimator may be chosen over the Ridge (γ = 2) estimator in advance. Our theoretical results suggest that γ is closely related to λ, in the asymptotic sense, similarly to [START_REF] Huang | A group bridge approach for variable selection[END_REF]; [START_REF] Knight | Asymptotics for lasso-type estimators[END_REF]). Thus an adaptive selection of λ in a finite sample is expected to reflect this relation automatically. This is confirmed numerically (see Section 4), suggesting a value of γ = 0.5 as a rule of thumb.

In addition, because the parametric model formulation arises as an approximation to the nonparametric model, the parameter space to explore is not fixed, and is potentially very large. Even with known covariates, the spline approximation with the fully adaptive choice of the degree, knot locations and number of knots is impractical. Following a similar strategy to that in the literature (e.g., [START_REF] Huang | Polynomial spline estimation and inference for varying coefficient models with longitudinal data[END_REF]; Wang, Li and Huang ( 2008)), we use equally spaced knots with cubic splines and select the number of knots K adaptively. We attempted to simultaneously optimize the parameter K inside the model selection criterion. However, we found that the penalty was not effective in controlling the systematic increase in the parameter space and that the criterion favored the smallest possible K in the majority of cases. Instead, we select the number of knots K adaptively to the sample using 10-fold cross-validation without a penalty, leaving the potentially adaptive choice of sparsity to be controlled solely by the other tuning parameters.

Large-Sample Properties

We study the large-sample properties of our proposed PLSE β i (t), for i = 1, . . . , p, when the number of sampled subjects n goes to infinity. In the proofs, we assume that the number of observations for each subject n k is bounded. However, a similar argument can be applied to the case when n k increases to infinity with n (Huang, [START_REF] Huang | Polynomial spline estimation and inference for varying coefficient models with longitudinal data[END_REF]). The number of interior knots increases with n; therefore, we write K i = K in for each i = 1, . . . , p, and denote K n = max 0≤i≤p K i . The standard regularity conditions for varying coefficient linear models [START_REF] Huang | Polynomial spline estimation and inference for varying coefficient models with longitudinal data[END_REF]; [START_REF] Wang | Variable selection in nonparametric varying-coefficient models for analysis of repeated measurements[END_REF]) are provided in the online Supplementary Material.

For mathematical convenience, we classify all group indices {1, . . . , G i } for the coefficient function β i (t) into two groups, defined as

A i1 = g : max t∈[ηg-1,ηg) | β i (t) |> C i K -2 n , A i2 = g : 0 ≤ max t∈[ηg-1,ηg) | β i (t) |≤ C i K -2 n ,
for some positive constant C i . For sufficiently large C i , the zero region {t :

β i (t) = 0} is a subset of ∪ g∈Ai2 [η g-1 , η g ).
Theorem 6.27 of [START_REF] Schumaker | Spline Functions: Basic Theory[END_REF] shows that any smooth coefficient function β i (t) with a bounded second derivative has a B-spline approximant β i (t) of the form given in (2.3) and that the approximation error is of order O(K -2 in ). We denote the sparse modification introduced in Wang and Kai (2015) by β 0 i (t), with coefficients α 0 .

Note that for a vector-valued square integrable function

A(t) = (a 1 (t), . . . , a m (t)) T , with t ∈ [0, M ], A 2 denotes the L 2 -norm defined by A 2 = ( m l=1 a l 2
2 ) 1/2 , where a l 2 is the usual L 2 -norm in the function space. Next, we establish the consistency of our proposed penalized estimator.

Theorem 1 (Consistency). Suppose that assumptions (A1)-(A6) in the online Supplementary Material are satisfied. For some 0 < γ < 1 and K n = O(n 1/5 ), we have the following assumption:

(S1) For α 0 defined above,

λ n (d + 1) 1/2   p i=1 g∈Ai1 α 0 Aig 2(γ-1) 1   1/2 = O(n 1/2 ) .
If (S1) holds, then we have ββ 2 = O p (n -2/5 ), where β = (β 1 , . . . , β p ) T . Assumption (S1) provides a bound on the rate of λ n growing with n. The convergence rate established in Theorem 1 is essentially optimal [START_REF] Stone | Optimal global rates of convergence for nonparametric regression[END_REF]). In fact, the result remains valid for more general classes of functions, for example, the collection of functions the derivatives of which satisfy the Hölder condition. Next, Theorem 2 states that our proposed penalized method is consistent in detecting functional sparsity. That is, if β i (t) = 0 for t ∈ [η l-1 , η l ), then the proposed estimator produces α Ail = 0 to identify local sparsity with probability converging to one. In addition, β i (t) = 0 for all t, then the proposed method yields α Ail = 0, for all l = 1, . . . , K i + 1, with probability converging to one.

Theorem 2 (Sparsistency). Consider the following assumption:

(S2) λ n K γ-1 n n -γ/2 -→ ∞ . If (S2
) and the assumptions in Theorem 1 are satisfied, then we have for every i, i = 1, . . . , p, ( α Aig : g ∈ A i2 ) = 0, with probability converging to one as n goes to ∞.

It is not surprising that our proposed method may yield a slightly more sparse functional estimate. This is because, for all intervals with indices belonging to A i2 , the value of β i (t) is quite small (the same order as the optimal rate) and is indistinguishable from zeros. Moreover, such intervals can be further partitioned into two groups, including the intervals on which the function is zero and the intervals on which the function is not always zero. However, the total length of the latter converges to zero as n increases.

The above discussion is related to the notion of selection consistency, an important and well-studied problem of variable selection under parametric settings; for instance, see (Zhao and Yu (2006)). However, for nonparametric models, particularly when local sparsity exists, selection consistency has not been widely studied. For the convenience of our discussion, we begin with some notation. For a coefficient function β(t), let N (β) and S(β) denote the zero region and nonzero region, respectively. The (closed) support of β, denoted by C(β), is defined as the closure of the nonzero region S(β). Assume that N (β) has finitely many singletons (as zero crossing), and that C(β) can be expressed as a finite union of closed intervals.

If β(t 0 ) = 0, for some t 0 , the consistency property in Theorem 1 and the smoothness constraint of the function and its estimate ensure that β(t 0 ) = 0 for sufficiently large n. However, such a result may not be of great interest, given that β(t) lies in an infinite, not necessarily countable dimensional space. Next, consider a simple case of an interval [a, b] ⊂ C(β) and β(t) = 0, for all t ∈ [a, b]. Thus, β(t) is bounded away from zero over [a, b]. Similarly, as a consequence of Theorem 1, β(t) is also bounded away from zero over [a, b] for sufficiently large n. A more challenging case arises when β(a) = 0 and β(t) = 0 over (a, b]. We further assume that there is a sequence of knots such that

η k ≤ a < η k+1 • • • < η k < b ≤ η k +1 .
The subinterval formed by two adjacent knots is either in A i1 or A i2 . The total length of the subintervals in A i2 converges to zero as n increases. For those intervals in A i1 , a suitable choice of the constant C i suggests that the estimated function deviates from zero.

Simulation Study

We conducted simulation studies to assess the performance of our proposed method, with the main emphasis on understanding the effects of the tuning parameters and the increasing dimension p on the functional sparsity estimation. We consider three scenarios. In Scenario 1, we choose our tuning parameters (λ, K) as described in Section 2.5, and compare the results under various γvalues. In Scenario 2, we assess the impact of the increasing dimension p, given γ, assuming the number of relevant covariates, p 0 , is set to four. In Scenario 3, we assess the performance with respect to K to study the effect of the adaptive choice of knots on the sparsity estimation. In addition, the relative performance is measured against that of the LSE and LASSO methods. The simulation results are summarized based on 400 replications. In each iteration, subjects are generated randomly according to the following varying coefficient model specification:

y k (t kl ) = p i=1 x ki (t kl )β i (t kl ) + k (t kl ), l = 1, . . . , n k , k = 1, 2, . . . , n,
where x 1 (t) is a constant-equal to one, x i (t), for i = 2, 3, 4, are similar to those considered in [START_REF] Huang | Varying-coefficient models and basis function approximations for the analysis of repeated measurements[END_REF]: x 2 (t) is a uniform random variable over [4t, 4t + 2]; x 3 (t), conditioning on x 2 (t), is a normal random variable with mean zero and variance (1 + x 2 (t))/(2 + x 2 (t)); and x 4 (t), independent of x 2 (t) and x 3 (t), is Bernoulli(0.6). The number of measurements available varies across subjects. For each subject, a sequence of 40 possible observation time points {(i -0.5)/40 : i = 1, . . . , 40} is considered, but each time point has a chance of 0.4 being selected. We further added a random perturbation from U (-0.5/40, 0.5/40) to each observation time. The random errors k (t kl ) are independent of the predictors, but include serial correlation and a measurement error of k (t) =

(1)

k (t) + (2)
k (t). The serial correlation component (1) (t) is generated from a Gaussian process with mean zero and covariance function cov(

(1) k (t),

(1) k (s)) = exp(-10|t -s|) for the same subject k, and is uncorrelated for different subjects. Then,

(2) k (t) follows a normal distribution with mean zero and variance one, and is independent and identically distributed (i.i.d.).

The nonzero coefficient functions used in all scenarios are displayed in Figure 2. The coefficient functions do not belong to the B-spline function space.

In Scenario 1, we add a redundant variable x 5 (t) from a normal distribution with mean zero and variance 0.1 exp(t) to illustrate global sparsity. In Scenario 2, with increasing p, the extra predictors with zero coefficient functions are defined as x i (t) = Z i (t) + 3/20 5 l=1 x l (t), for i = 6, . . . , p, where Z i (t) is i.i.d. from a standard normal distribution.

The overall performance is measured in terms of bias and the mean integrated squared error (MISE), based on R = 400 repetitions, computed as

Bias i (u) = 1 R R r=1 β i (r) (u) -β i (u), i = 1, . . . , p, u ∈ [0, 1],
Figure 2. A graphical illustration of the coefficient functions β i , for i = 1, . . . , 4 (from left to right, top to bottom).

M ISE i = 1 R R r=1 1 0 ( β (r) i (u) -β i (u)) 2 du, i = 1, . . . , p, where β (r) i
is the estimated coefficient function from the rth repeated study. In addition, we use the following summary measures to compare the functional sparsity: Scenario 1: Effect of γ

Here, we consider the varying coefficient model with p = 5 and two different numbers of subjects, n = 100, 200. In each iteration, our proposed PLSE method is implemented with γ = 0.25, 0.35, 0.5, and 0.75. The MISE values for each coefficient function are summarized in Table 1. In general, as n increases, all methods yield decreasing MISE values. Notably, the results for the PLSE indicate comparable performance across different γ; in fact, the PLSE and LASSO methods show similar performance in terms of the function estimation. In addition, the PLSE with γ = 0.35, 0.50 successfully identifies the global sparsity of β 5 (•) with zero MISE values for both choices of n, as does the PLSE with γ = 0.75 for n = 200. The bias of PLSE 0.5 (PLSE with γ = 0.5) and the LASSO and LSE methods with n = 200 is compared in Figure 3, showing that PLSE 0.5 has zero bias in estimating β 5 (•).

Table 2 describes the performance when identifying local sparsity. The true values of sparsity in terms of C i,0 and I i,0 are given in the last row of true model as a reference. Hence, the closer the values of C i,0 are to those of the true model, the better. In contrast, the value of I i,0 in the true model is the maximum error each method can make; thus, the smaller I i,0 , the better. In general, the LASSO Here, for the true model, C i,0 , for i = 1, . . . , 6, are the lengths of the zero intervals, I i,0 denotes the length of a nonezero interval, C 0 is the number of zero-coefficient functions, and I 0 is the number of nonzero-coefficient functions.

Method

C 1,0 I 1,0 C 2,0 I 2,0 C 3,0 I 3,0 C 4,0 I 4,0 C 5,0 I 5,0 C 0 I 0 n = 100 LSE 0 0 0 0 0 0 0 0 0 0 0 0 LASSO 0 0 0.0219 0 0 0 0.0616 0.0009 0.8799 0 0.5675 0 PLSE 0.25 0 0 0.1468 0.0003 0 0 0.1626 0.0004 0.8005 0 0.4975 0 PLSE 0.35 0 0 0.3332 0.0048 0 0 0.3723 0.0062 1.0000 0 1 0 PLSE 0.5 0 0 0.3360 0.0040 0 0 0.3809 0.0072 1.0000 0 1 0 PLSE 0.75 0 0 0.2559 0.0005 0 0 0.3453 0.0042 0.9990 0 0.9975 0 n = 200 LSE 0 0 0 0 0 0 0 0 0 0 0 0 LASSO 0 0 0.0166 0 0 0 0.0736 0.0004 0.9087 0 0.6850 0 PLSE 0.25 0 0 0.1299 0.0001 0 0 0.1433 0.0003 0.7510 0 0.4175 0 PLSE 0.35 0 0 0.3178 0.0022 0 0 0.3512 0.0030 1.0000 0 1 0 PLSE 0.5 0 0 0.3343 0.0025 0 0 0.3735 0.0047 1.0000 0 1 0 PLSE 0.75 0 0 0.2696 0.0005 0 0 0.3462 0.0026 1.0000 0 1 0 true model 0 1 0.4000 0.6000 0 1 0.4286 0.5714 1 0 1 4 and PLSE methods show better performance in terms of functional sparsity. In addition, it can be seen that the PLSE with γ = 0.35, 0.5, 0.75 has an advantage in achieving both global and local sparsity, as compared with the LSE and LASSO methods. The case of γ = 0.5 performs slightly better than the others. For the remaining part, we use γ = 0.5 for the comparisons.

Scenario 2: Effect of dimension p

In this scenario, we study the effect on the performance of increasing p, for a given sample size. In particular, we consider the three choices of p = 5, 20, and 50. Figure 4 and Table 3 show the results for the bias and MISE, respectively. The last column of Tables 3 shows the maximum MISE among the zero-coefficient functions as the selected variables vary between samples. Compared with the LSE and LASSO, the PLSE method shows remarkable stability in terms of performance over the increasing dimension p.

The performance for the sparsity is summarized in Table 4. The additional two columns in C i,0 and I i,0 are added to summarize the performance on all other redundant variables within an interval range of [min i≥6 C i,0 , max i≥6 C i,0 ] and [min i≥6 I i,0 , max i≥6 I i,0 ]. Together with the global sparsity measure in C 0 and I 0 , we conclude that PLSE 0.5 systematically outperforms the other methods for all dimensions. The variation in the selection of the knots is expected to mainly influence the estimation of local sparsity. Increasing the number of knots helps to identify the boundary of local sparsity, but runs the risk of over-fitting nonzero estimates. Fine-tuning this parameter is much more delicate, because all model selection criteria are developed to control the squared error loss (MISE) as a goodnessof-fit, and thus are insensitive to the loss of missing local sparsity. That is, the balance between global and local sparsity is beyond the usual control of the bias-variance trade-off, and developing a new measure is still an open problem. Our knot selection based on cross-validation is essentially tuned toward global sparsity. Here, we assess the performance of our proposed estimator from the point of view of robustness to these variations. For comparison, we include the results for fixed knots (K = 11) across the sample.

Figure 5 and Table 5 summarize the bias and MISE. The sparsity summary is given in Table 6. We conclude that the overall performance is fairly comparable to that in Scenario 2, with no major concern over the sensitivity of the knots selection in the comparison of the results.

In addition, in order to assess the usefulness of the asymptotic formula for the standard errors in (2.6), we calculate both the asymptotic and the empirical standard errors based on 400 repetitions. Then, we compare these standard errors in Figure 6 for an adaptive number of knots, and in Figure 7 for a fixed Table 4. Sparsity summary measures (a)-(d) for Scenario 2. Here, for the true model, I i,0 , for i = 1, . . . , 4, are the lengths of the nonzero intervals and C i,0 are the lengths of the zero intervals. Method number of knots, with the results showing good agreement. Thus, the variation in the number of knots greatly increases the variation in the estimation of the Here, for the true model, I i,0 , for i = 1, . . . , 4, are the lengths of the nonzero intervals and C i,0 are the lengths of the zero intervals. Method In summary, the simulation results demonstrate that our proposed method not only has an advantage in achieving local sparsity compared with the LASSO and LSE methods, but also ensures global sparsity for finite-dimensional models. Moreover, this advantage applies to models with an increasing dimension.

C 1,0 I 1,0 C 2,0 I 2,0 C 3,0 I 3,0 C 4,0 I 4,0 C 5,0 I 5,0 [C (min) i,0 , C (max) i,0 ] [I (min) i,0 , I (max) i,0 ] C 0 I 0 p = 5
C 1,0 I 1,0 C 2,0 I 2,0 C 3,0 I 3,0 C 4,0 I 4,0 C 5,0 I 5,0 [C (min) i,0 , C (max) i,0 ] [I (min) i,0 , I (max) i,0 ] C 0 I 0 p = 5

Real-Data Analysis

We demonstrate our method by analyzing yeast cell cycle gene expression data [START_REF] Spellman | Comprehensive identification of cell cycle-regulated genes of the yeast saccharomyces cerevisiae by microarray hybridization[END_REF]; [START_REF] Lee | Transcriptional regulatory networks in saccharomyces cerevisiae[END_REF]).

In the biological sciences, gene expression data are common. Scientists believe that transcription factors (TFs) might have an effect on a genome's cell cycle regulation. As a result, they have attempted to identify the key TFs in the regulatory network, based on a set of gene expression measurements. In this study, we analyze the relationship between the level of a gene expression and the physical binding of the TFs from chromatin immunoprecipitation (ChIP-chip) data [START_REF] Lee | Transcriptional regulatory networks in saccharomyces cerevisiae[END_REF]). One set of gene expression data comes from an α-factor synchronization experiment of 542 genes, in which mRNA levels are measured every 7 minutes during a period of 119 minutes, resulting in 18 measurements in total [START_REF] Spellman | Comprehensive identification of cell cycle-regulated genes of the yeast saccharomyces cerevisiae by microarray hybridization[END_REF]). For our analysis, the time has been rescaled to [0,1].

The ChIP-chip data contain the binding information of 106 TFs, of which 21 TFs are confirmed to be related to cell cycle regulation, by experiment. [START_REF] Wang | Group scad regression analysis for microarray time course gene expression data[END_REF] demonstrated that a variable selection procedure is able to identify some of the key TFs. It is believed that the effects of TFs vary during the cell cycle. In [START_REF] Chun | Sparse partial least squares regression for simultaneous dimension reduction and variable selection[END_REF], the authors considered a sparse partial least squares regression to study which TFs are important in a gene expression. However, they did not focus on the active periods of TFs. In this study, we apply our method to identify the key TFs, and to estimate the effects of those selected TFs over time. In addition, our approach allows us to investigate whether active and inactive periods during the cycle can be identified for each TF. Let y kt denote the gene expression level for gene k at time t, for k = 1, . . . , 542 and t = 1, . . . , 18, and let x ki denote the binding information of TF i for gene k, for i = 1, . . . , 106. Then, the varying coefficient model can be written as where β i (t) models the effect of the ith TF on a gene expression at time t, and for the kth gene, kt are independent over time.

Similarly to the simulation study, we apply our method together with LSE and LASSO methods, and compare the identification of the active period of each TF within the cell cycle process. Each coefficient function is approximated using quadratic B-splines defined on time interval [0, 1], with seven equally spaced knots. The number of knots is selected using cross-validation. It is not surprising that the LSE selects all TFs. The LASSO method identifies 32 TFs as important, while our proposed method identifies 16 TFs, which are a subset of those identified by the LASSO method. In Figure 8, the estimated coefficient functions for 21 experimentally confirmed TFs are shown. The figure shows that eight are selected by both methods. The LASSO method selects an additional four TFs, namely, SWI4, STB1, FKH1, and REB1, which show very low levels of activities. In [START_REF] Chun | Sparse partial least squares regression for simultaneous dimension reduction and variable selection[END_REF], the authors selected 32 TFs, 10 of which are verified TFs. In addition, our proposed method identifies some inactive periods for selected TFs. For example, STE12 tends to be inactive in the later period, and ACE2 is inactive in the early period.

  (a) C 0 : average number of correctly identified constant zero-coefficient functions (b) I 0 : average number of incorrectly identified constant zero-coefficient functions (c) C i,0 : average length of correctly identified zero intervals for the ith coefficient function (d) I i,0 : average length of incorrectly identified zero intervals for the ith coefficient function. Note that (a) and (b) summarize global sparsity, whereas (c) and (d) summarize local sparsity.

Figure 3 .

 3 Figure 3. Comparison of the bias of the coefficient functions based on the LSE (dotdashed), the LASSO (dashed), and PLSE 0.5 (solid) in Scenario 1, with n = 200. Note that PLSE 0.5 has zero bias in estimating the zero coefficient function β 5 (•).

Figure 4 .

 4 Figure 4. Comparison of the bias of the nonzero coefficient functions β 1 , β 2 , β 3 , and β 4 (from left to right) based on the LSE (dot-dashed), LASSO (dashed), and PLSE 0.5 (solid), for p = 5 (top row), p = 20 (middle), and p = 50 (bottom) in Scenario 2.

Figure 5 .

 5 Figure 5. Comparison of the bias of the nonzero coefficient functions β 1 , β 2 , β 3 , and β 4 (from left to right) based on the LSE (dot-dashed), LASSO (dashed), and PLSE 0.5 (solid), for p = 5 (top row), p = 20 (middle), and p = 50 (bottom) in Scenario 3.

Figure 6 .

 6 Figure 6. Asymptotic standard error (grey solid line) and empirical standard deviation (black solid line) of the coefficient functions with an adaptive number of knots in Scenario 2.

Figure 7 .

 7 Figure 7. Asymptotic standard error (grey solid line) and empirical standard deviation (black solid line) of the coefficient functions with a fixed number of knots in Scenario 3.

Figure 8 .

 8 Figure8. Subplots of estimated coefficient functions for the 21 confirmed TFs using the LSE (dashed), LASSO (dot-dashed), and PLSE 0.5 (solid).

y

  kt = β 0 (t) + 106 i=1 β k (t)x ik + kt ,

Table 1 .

 1 Comparison of MISE for each coefficient function in Scenario 1.

	Method			MISE
		β 1	β 2	β 3	β 4	β 5
				n = 100
	LSE	0.9519 0.0825 0.0365 0.1145 1.3199
	LASSO	2.9156 0.1636 0.0314 0.0591 0.0114
	PLSE 0.25 1.1115 0.0686 0.0281 0.0613 0.1330
	PLSE 0.35 1.2199 0.0633 0.0307 0.0440 0
	PLSE 0.5	1.3156 0.0674 0.0319 0.0459 0
	PLSE 0.75 1.8267 0.0948 0.0317 0.0471 0.0005
				n = 200
	LSE	0.4232 0.0367 0.0165 0.0563 0.5745
	LASSO	1.4561 0.0845 0.0153 0.0299 0.0041
	PLSE 0.25 0.7259 0.0424 0.0138 0.0329 0.0731
	PLSE 0.35 0.6421 0.0351 0.0152 0.0235 0
	PLSE 0.5	0.7193 0.0382 0.0166 0.0250 0
	PLSE 0.75 0.8615 0.0469 0.0157 0.0251 0

Table 2 .

 2 Sparsity summary measures (a)-(d) in Scenario 1.

Table 3 .

 3 Comparison of MISE for each coefficient function with p = 5, 20, and 50 in Scenario 2.

	Method			MISE		max i≥6 MISE i
		β 1	β 2	β 3	β 4	β 5
				p = 5
	LSE	0.4232 0.0367 0.0165 0.0563 0.5745 -
	LASSO	1.4561 0.0845 0.0153 0.0299 0.0041 -
	PLSE 0.5	0.7193 0.0382 0.0166 0.0250 0	-
				p = 20
	LSE	0.5157 0.0434 0.0197 0.0612 0.6694 0.0151
	LASSO	17.8758 0.8520 0.0331 0.0347 0	0.0016
	PLSE 0.5	0.7422 0.0391 0.0166 0.0240 0	2.1886e-05
				p = 50
	LSE	0.8269 0.0724 0.0292 0.0897 1.1484 0.0281
	LASSO	35.1543 1.6475 0.0497 0.0415 0	8.4758e-04
	PLSE 0.5	0.7360 0.0396 0.0149 0.0205 0	2.7551e-05
	Scenario 3: Effect of knots selection

Table 5 .

 5 Comparison of MISE for each coefficient function in Scenario 3 . Here, the number of knots is fixed at 11.

	Method			MISE		max i≥6 MISE i
		β 1	β 2	β 3	β 4	β 5
				p = 5
	LSE	0.2783 0.0253 0.0108 0.0379 0.3888 -
	LASSO	1.2993 0.0753 0.0096 0.0195 0.0029 -
	PLSE 0.5	0.6888 0.0405 0.0107 0.0154 0	-
				p = 20
	LSE	0.3376 0.0292 0.0127 0.0408 0.4429 0.0097
	LASSO	11.6604 0.5420 0.0199 0.0211 0	0.0011
	PLSE 0.5	0.7180 0.0412 0.0114 0.0156 0	2.1444e-05
				p = 50
	LSE	0.4521 0.0387 0.0167 0.0528 0.6608 0.0138
	LASSO	30.2698 1.3683 0.0352 0.0290 0	7.0743e-04
	PLSE 0.5	0.7279 0.0417 0.0114 0.0150 0	1.9744e-05

Table 6 .

 6 Sparsity summary measures (a)-(d) for Scenario 3.

Acknowledgments

The authors are grateful to the three referees and the associate editor for their careful reading and helpful comments. The research of Haonan Wang was partially supported by NSF grants DMS-1106975, DMS-1521746, and DMS-1737795.

Supplementary Material

The online Supplementary Material includes the technical assumptions and proofs of the theoretical properties of our proposed method.


Ciara R Wigham

Kristi Jauregi-Ondarra

Cathryn Bennett

Trainee teachers' engagement in initial synchronous virtual exchange interactions

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Trainee teachers' engagement in initial synchronous virtual exchange interactions

Ciara R. Wigham, Laboratoire ACTé, Université Clermont Auvergne Kristi Jauregi-Ondarra, Utrecht Univeristy Cathryn Bennett, Laboratoire ACTé, Université Clermont Auvergne Engagement has been discussed in relation to the effort invested during task accomplishment (Bygate & Samuda, 2009;Helme & Clarke, 2001;Skinner & Pitzer, 2012). Recent studies have been framed by social constructivist views of task processing (Pekrun & Linnenbrink-Garcia, 2012;Philp & Duchesene, 2016) and have approached task engagement more holistically as the effort invested in accomplishing a task, guided by participants' requirements of task success and involving behavioral, cognitive and attitudinal qualities in interaction with each other and for which the social element permeates the aforementioned three dimensions (Gijsen, 2021).

This study addresses the task performance parameter of partner orientation: the effort participants invest in establishing and deepening a relationship with partners to strengthen interaction (Gijsen, 2021). We adopt multimodal (inter)action analysis (Norris, 2004) to investigate how different semiotic resources contribute to partner orientation.

Our context is a virtual exchange within the E-LIVE project (Engaging Languages in Intercultural Virtual Exchange) involving 48 trainee language teachers from French, Dutch, and Colombian institutions. Our analysis focuses on the first synchronous interaction via videoconferencing in which teaching experiences were discussed. This followed an asynchronous language biography task.

Using multimodal transcriptions from a subset of interactions, we examine three dimensions of partner orientation: The behavioral dimension with performance indicators of providing and accepting support/feedback, giving compliments, using comprehension checks/ production prompts; the attitudinal dimension through comments on working together, task ownership, technological tools used; and the cognitive dimension explored through reflective blogs and questionnaire/focus group data during which participants reflected on the relevance of strengthening the interaction with their partners, the effort they made to reach mutual understanding and address intercultural aspects.

We anticipate study findings will advance our understanding of task engagement and help identify effective strategies to maximise engagement, and thus learning, in virtual exchange projects.