Daniel Martin Xavier
email: daniel.martinxavier@ens-paris-saclay.fr

Ludovic Chamoin
email: ludovic.chamoin@ens-paris-saclay.fr

Laurent Fribourg
email: laurent.fribourg@ens-paris-saclay.fr

Data-driven MPC applied to non-linear systems for real-time applications *

Model Predictive Control (MPC) is a traditional technique widely employed on the control of constrained non-linear systems. In light of the increasing popularity of neural networks, data-driven MPC has emerged as an alternative to alleviate the computation burden of traditional control strategies. This paper aims to replace the constrained optimization problem by training a feed-forward neural network with data collected from an offline MPC simulation. The network's performance is then compared to that of traditional MPC using the Van der Pol oscillator as toy example. Finally, the convergence of the network training error is analytically proven by extending the analysis of neural tangent kernels (NTK) to underparameterized networks.

Introduction

The most prevalent type of control strategy in industrial applications is the feedback controller, which involves the comparison of a reference signal r with a measured variable y to determine an appropriate value for the manipulated variable u using the error e = ry. This broad class of control strategies can be divided into three categories: classical controllers such as Proportional-Derivative-Integral (PID), predictive controllers such as Model Predictive Control (MPC) and Linear Quadratic Regulator (LQR), and repetitive controllers [START_REF] Schwenzer | Review on model predictive control: an engineering perspective[END_REF].

Within the domain of classical control strategies, the PID controller stands out as the most renowned, being extensively employed in industrial applications. Over the years, several setup rules have been proposed to adjust the controller's parameters, however finding an optimal parameterization can often be a challenging task. This is especially true for nonlinear and time-variant systems, where even higher control methods, such as back-stepping controllers, struggle to yield a suitable performance.

Predictive controllers, on the other hand, use a model of the system (digital twin) to predict its future behavior, being able to anticipate deviations from the reference. Within this class of controllers, MPC assumes a prominent role, relying on a repeated real-time optimization of a mathematical system model. Using the predictions of future states, MPC is able to determine the optimal trajectory of the manipulated variable u, which is then adjusted at each time step by subsequent optimizations.

Model Predictive Control offers several advantages over classical control techniques, as it can anticipate the system's behavior and naturally consider hard constraints on the optimization problem. Furthermore, the parameterization of the controller is often more intuitive compared to other methods: the performance is usually improved at the cost of a higher computational effort, which is defined by a weighting coefficient.

In certain real-world applications, however, the mathematical model of the system may not be available due to the complexity of modeling or the lack of a priori knowledge regarding the phenomenon of interest. In such cases, data-driven models have proven to be effective to perform predictive control, using data obtained from sensors to identify the underlying dynamics of the system through auto-regressive, machine learning or physics-guided models [START_REF] Kaheman | Learning discrepancy models from experimental data[END_REF].

Nonetheless, there has been limited progress in mitigating the computational burden of the MPC controller on real-time applications. That's why the present paper investigates the use of data-driven MPC to achieve fast computation of the control law while providing an analytical guarantee of the neural networks (NN). The goal is to replace the traditional constrained optimization problem in MPC by a NN to control a nonlinear system.

In recent years, the interest in neural networks has increased due to its remarkable performance across a wide range of applications. However, the comprehension of why such networks converge when trained via gradient descent, and generalize well over unseen data remained an open question for a while. Some studies have attempted to characterize the geometric landscape of objective functions in order to explain training convergence [START_REF] Nguyen | The loss surface of deep and wide neural networks[END_REF][START_REF] Freeman | Topology and geometry of half-rectified network optimization[END_REF].

Another approach consisted in investigating the trajectory of optimization and demonstrating convergence toward a global minimum [START_REF] Li | Convergence analysis of two-layer neural networks with relu activation[END_REF][START_REF] Du | Gradient descent learns one-hidden-layer cnn: Don't be afraid of spurious local minima[END_REF]. In this context, Zhang et al. [START_REF] Zhang | Understanding deep learning requires rethinking generalization[END_REF] were among the first to study the generalization of neural networks, showing that they can attain zero training error on sufficiently large nets. This class of NNs are called overparameterized neural networks, characterized by a number of parameters m much larger than the number of data n (m ≫ n).

Recently, the notion of Neural Tangent Kernel (NTK) matrix was introduced, providing an elegant explanation of the linear convergence towards 0 of the training error in overparameterized NNs with smooth activation functions [START_REF] Jacot | Neural tangent kernel: Convergence and generalization in neural networks[END_REF]. The explanation essentially lies in the fact that, during the Gradient Descent (GD) procedure, the n × n NTK matrix stays close to its infinite limit, with all eigenvalues remaining positive.

Thereafter, Du et al. [START_REF] Du | Gradient descent provably optimizes over-parameterized neural networks[END_REF][START_REF] Du | Gradient descent finds global minima of deep neural networks[END_REF] expanded this analysis to prove convergence of neural networks with non-smooth and non-convex activation functions. Arora et al. [START_REF] Arora | Fine-grained analysis of optimization and generalization for overparameterized two-layer neural networks[END_REF] also employed this framework to analyze convergence rates when using random and true labels for gradient descent.

However, most of the focus of these papers lies in the analysis of overparameterized neural networks, which is not the case of a considerable number of real-world applications.

In this paper, we address the problem of underparameterized neural networks (m ≪ n), where the NTK matrix is positive semi-definite. Thus, its nullspace N is no longer reduced to 0, but it is a space of dimension larger than 1. Moreover, the NTK matrix evolves over time and its coefficients do not remain close to those of initialization. Nevertheless, under certain assumptions, this paper demonstrates that the training error converges linearly to a constant which can be accurately evaluated at the beginning of the GD procedure.

The remainder of the paper proceeds follows this structure: Section 2 presents the traditional MPC strategy and data-driven approaches with preliminary results using the Van der Pol oscillator as a toy example. Section 3 outlines the theoretical guarantees of the network convergence and the results are illustrated using the same toy example. Finally, Section 4 contains the conclusions and prospects for this work.

Model Predictive Control

Model Predictive Control (MPC) is a receding horizon strategy that uses a model of the system being controlled to predict its future states and compute suitable commands through a constrained optimization problem [START_REF] Schwenzer | Review on model predictive control: an engineering perspective[END_REF]. This strategy was first applied in chemical plants and oil refineries in order to deal with constraints in nonlinear systems [START_REF] Qin | A survey of industrial model predictive control technology[END_REF]. Following its successful industrial implementation in the 1980s, the research community introduced theoretical tools to mathematically assure the stability and optimality of MPC controllers [START_REF] Mayne | Constrained model predictive control: Stability and optimality[END_REF].

The advantage of MPC over other predictive strategies lies in its ability to continuously recompute a new trajectory based on data measured at each time step. Consequently, the controller is robust by design, capable of accounting for unexpected disturbances by dynamically recalculating paths to steer the system back to the desired state.

Data-driven MPC

Current research focuses on replacing the mathematical model on traditional MPC by datadriven models, which do not require a priori information about the system dynamics. Draeger et al. [START_REF] Draeger | Model predictive control using neural networks[END_REF], for instance, employed a MPC based on the predictions generated by a neural network to control the pH levels in a neutralization reactor. Furthermore, Kaiser et al. [START_REF] Kaiser | Sparse identification of nonlinear dynamics for model predictive control in the low-data limit[END_REF] proposed the Sparse Identification of Nonlinear Dynamics (SINDY) to determine the prediction model of the system through a collection of predefined functions.

Alternative approaches concentrate on the derivation of the feasibility set or the cost function by data-driven methods in order to improve the controller performance. In this context, Berberich et al. [START_REF] Berberich | On the design of terminal ingredients for data-driven mpc[END_REF] proposed a methodology to compute the terminal cost and terminal set for data-driven MPC using measured input-output data. Furthermore, Collet et al. [START_REF] Collet | Data-driven fatigueoriented mpc applied to wind turbines individual pitch control[END_REF] presented a fatigue-oriented MPC approach that employs a data-driven cost function to optimize the fatigue trade-off on wind turbines.

This paper aims to replace the constrained optimization problem in traditional MPC by a neural-network to alleviate its computational burden. The proposed strategy involves trading optimality for a fast computation of the control law. The neural network is trained using data generated by running MPC simulations offline. As a result, the data-driven MPC controller proposed in this paper is more suitable for real-time applications.

The majority of works in this topic focus on using neural networks to approximate the performance of an Explicit MPC controller. Chen et al. [START_REF] Chen | Approximating explicit model predictive control using constrained neural networks[END_REF], for instance, proposed a modified reinforcement learning policy for controlling a linear system, ensuring feasibility by projecting the control inputs onto polytopes. Furthermore, Winqvist et al. [START_REF] Winqvist | On training and evaluation of neural network approaches for model predictive control[END_REF] developed a 2-layer neural network with constrained ReLU based activation function, incorporating a projection layer to guarantee recursive feasibility and asymptotic stability.

An Explicit MPC computes the control law offline in the form of a piecewise affine function for each control region of the state space. These regions are convex polytopes for linear MPC, making this strategy unsuited for highly nonlinear systems as the number of regions grows exponential with the model complexity. This work stands out from the existing literature by proposing the use of neural networks to replace the optimization problem on the control of non-linear systems.

To validate this approach, the proposed strategy is tested on a toy example of the Van der Pol oscillator, which is presented in the following subsection. In the remainder of this paper, the MPC problem is implemented using the MPC Toolbox in Python [START_REF] Lucia | Rapid development of modular and sustainable nonlinear model predictive control solutions[END_REF], which solves the differential algebraic equations of the nonlinear model through the Sundials CVODE, and the optimization problem using the Interior-Point Method [START_REF] Yu | Interior-point methods: An old and new approach to nonlinear programming[END_REF].

Toy Example -Van der Pol Oscillator

The Van der Pol oscillator is a classical proof-of-concept dynamical system for optimal control applications. It was originally introduced by Balthazar Van der Pol for modeling the triode oscillations in electrical circuits [START_REF] Aislan Antonelo | Physics-informed neural nets for control of dynamical systems[END_REF]. The system possesses 2 states x = [x 1 , x 2 ], a control action u, and a damping coefficient µ = 1. Its dynamics are described by the following equations:

ẋ1 = x 2 ẋ2 = µ(1 -x 2 1 )x 2 -x 1 + u (1)
The goal is to use MPC to make the system converge to a desired position x ref 1 , and then get this data to train offline a feed-forward neural network. The MPC problem is formulated using a time step of T s = 0.5s, an initial condition of x 0 = [1, 0] and a prediction horizon of N = 5. The command applied to the oscillator is constrained to the interval u ∈ [-1, 1] and the cost function associated to the optimization problem is defined as:

J(x, u) = N -1 i=0 ∥x ref 1 (k + i) -x 1 (k + i)∥ 2 +α∥∆u(k + i)∥ 2 (2) s. t. u lb ≤ u(k + i) ≤ u ub , ∀i ∈ [0, ..., N -1] (3)
where u lb and u up represent the lower and upper bounds for the command, respectively, ∆u(k) = u(k) -u(k -1) the command variation, and α = 0.1 a weighting coefficient that regulates the trade-off between reducing the error and increasing the command. Furthermore, the current time step is denoted by k, and the error is represented by ϵ

(k) = x 1 -x ref 1 .
The MPC simulations are conducted offline through the definition of different references that are randomly created between [-1,1], each reference being applied for 10s. The nonlinear model of the system is implemented using CasADi [START_REF] Joel A E Andersson | CasADi -A software framework for nonlinear optimization and optimal control[END_REF], being defined using the same time step as the MPC controller (T s = 0.5s). The MPC simulation is 1800s long, resulting in a dataset of S = 3600 data points, which is divided into 3 sets: S train = 2160 for training, S val = 720 for validation, and S test = 720 for testing.

Once the data is collected, the neural network is designed seeking a good balance between model complexity and performance. Its architecture is depicted in Figure 1, being a one-stepahead predictive controller with one hidden layer. The network has 4 entries in the input layer, 10 neurons in the hidden layer, and 1 output in the last layer. It was implemented using PyTorch [START_REF] Paszke | Pytorch: An imperative style, highperformance deep learning library[END_REF] with a ReLU activation function applied to the hidden layer, and a Hardtanh to the output layer to constrain the command within the interval [-1, 1].

Input Layer ¡ ¢

Hidden Layer ¡ ¢¹¤ Output Layer ¡ ¢¹ The neural network was trained using the AdamW [START_REF] Loshchilov | Decoupled weight decay regularization[END_REF] optimization algorithm during 10 4 epochs with a learning rate of 10 -3 , and a weight decay of 10 -6 . The performance of the trained network is assessed using two metrics: the root mean squared error (RMSE) in (4), and the coefficient of determination (R2 score) in [START_REF] Chen | Approximating explicit model predictive control using constrained neural networks[END_REF]. In these equations, u MPC i represents the command computed by MPC, u NN i the command predicted by the neural network, and ūMPC an average of MPC commands.

RMSE = 1 S test Stest i=1 u MPC i -u NN i 2 (4) R 2 = 1 - Stest i=1 u MPC i -u NN i 2 Stest i=1 ūMPC -u NN i 2 (5)
The trained network is capable of simulating the MPC behavior over the test dataset with an error of RMSE = 0.0314, indicating a satisfactory level of generalization over unseen data. Furthermore, the NN obtains a high coefficient of determination R 2 = 0.9975 ∼ 1, showing its efficiency in the regression task. Figure 2 illustrates the command computed by MPC compared to the one predicted by the neural network using the test set, highlighting its capacity in learning the behavior of an optimization problem. Finally, to evaluate the performance of the data-driven MPC, the NN is used to control the Van der Pol oscillator online. The simulation is initialized at x 0 = [-0.3, 0], following the same reference trajectory as in the testing dataset. Figure 3 presents the simulation results, illustrating the capability of the neural network to effectively steer the system to the desired reference. In comparison to traditional MPC, the proposed strategy yields a trajectory with more oscillations, which is empirically verified by the RMSE as shown in Table 1. Remark 1. These closed-loop oscillations suggest that the traditional MPC was not finely tuned to achieve zero overshoot. As the NN is trained to reproduce the MPC performance, the network will learn to oscillate if this is the trajectory yield by MPC. In order to avoid such oscillations, the traditional MPC should be carefully tuned by adjusting the prediction horizon N and the weighting coefficient α. Another option is to enforce precision by penalizing the deviations from the reference in the function being minimized by the GD algorithm.

However, the primary goal of the proposed strategy is to prioritize fast computation of the control law at the expense of optimality, which was satisfactorily achieved using neural networks. The average run-time of traditional MPC is 2.4201 ms, over 50 times slower than data-driven MPC, primarily due to the time-consuming nature of the constrained optimization problem. This difference would only become more pronounced with a larger prediction horizon N > 5, confirming the suitability of data-driven over traditional MPC for real-time applications.

Convergence of underparameterized neural networks

In recent years, the interest in neural networks has increased due to its outstanding performance across a wide range of applications. Nevertheless, they lacked a formal formulation of theoretical properties regarding their guarantees of convergence for years. Within the existing literature, Jacot et al. [START_REF] Jacot | Neural tangent kernel: Convergence and generalization in neural networks[END_REF] proposed that the convergence of gradient descent could be characterized by a kernel in the case of overparameterized neural networks with smooth activation functions. Subsequently, Du et al. [START_REF] Du | Gradient descent provably optimizes over-parameterized neural networks[END_REF] expanded this analysis to prove the convergence of neural networks with non-smooth activation functions. In this paper, the analysis of convergence using a kernel is further expanded to encompass underparameterized networks, in which the number of trainable parameters are smaller compared to available data (m ≪ n). After some preliminaries (Subsection 3.2), the proof of the error convergence on underparameterized neural networks is presented (Subsection 3.3), followed by some numerical results on the Van der Pol oscillator (Subsection 3.4).

Notation

In this paper, we denote by R and N the set of real and natural numbers, respectively. These symbols are annotated with subscripts to restrict them in the usual way, e.g., R >0 denotes the positive real numbers. We also denote by R n an n-dimensional Euclidean space, and by R n×m a space of real matrices with n rows and m columns.

We use bold letters for vectors and bold capital letters for matrices. Given a matrix A, let A i,j be its (i, j)-th entry, λ min (A) its minimal eigenvalue, and A ⊤ its transpose. The Euclidean norm is denoted by ∥ • ∥, the Frobenius norm by ∥ • ∥ F , and the inner product by ⟨•, •⟩. Let I n be the n × n identity matrix, [n] the set {1, . . . , n}, and σ(•) the ReLU function σ(z) = max{z, 0}.

We denote by I{E} the indicator function for an event E, by B ⊎ C the disjoint union of sets B and C, and by V 1 ⊕V 2 the direct sum of vector spaces V 1 and V 2 . We also use the abbreviation i.i.d. to indicate that a collection of random variables is independent and identically distributed. Finally, the time-discrete version of a time-continuous object ξ is denoted as ξ.

Preliminaries

We recall from Du et al. [START_REF] Du | Gradient descent provably optimizes over-parameterized neural networks[END_REF] some definitions regarding the application of the gradient descent algorithm to neural networks. We consider a NN of the form: This work focuses on the empirical risk minimization problem with the quadratic loss defined in [START_REF] Draeger | Model predictive control using neural networks[END_REF]. Indeed, given a training data set S = {(x i , y i )} n i=1 , the objective of the optimization algorithm is to minimize the following loss function:

f (θ, x) = 1 √ m m r=1 a r σ(w ⊤ r x) (6
L(θ) = 1 2 n i=1 (f (θ, x i ) -y i ) 2 (7)
In order to do so, Du et al. [START_REF] Du | Gradient descent provably optimizes over-parameterized neural networks[END_REF] fixed the parameters of the second layer of the network and applied the gradient descent algorithm to the weights of the first layer. In this paper, however, both layers are optimized simultaneously during training through the following update rule:

θ k+1 = θ k -h ∂L(θ k ) ∂θ k (8)
The gradient formula for each layer is, then, given by:

∂L(θ) ∂w r = 1 √ m n i=1 (f (θ, x i ) -y i )a r x i I{w ⊤ r x i ≥ 0} (9) ∂L(θ) ∂a r = 1 √ m n i=1 (f (θ, x i ) -y i )σ(w ⊤ r x i ) ( 10
)
Given the discontinuous nature of the ReLU activation function, one may consider ∂L(θ) ∂θ as a convenient notation for the right hand side of ( 9) and [START_REF] Du | Gradient descent provably optimizes over-parameterized neural networks[END_REF]. Lastly, the discrete equation ( 8) corresponds to the Euler discretization, with a step size h, of the set of ordinary differential equations defined by:

dθ(t) dt = - ∂L(θ(t)) ∂θ(t) . ( 11
)

Gradient Descent for Neural Networks

When training a neural network, the convergence of the GD algorithm to a globally optimal solution comes down to showing the convergence of the error (i.e. difference between the prediction of the NN and the ground truth) to zero [START_REF] Du | Gradient descent provably optimizes over-parameterized neural networks[END_REF]. For this reason, the rest of this paper will focus on the analysis of the error dynamics v : R ≥0 → R n , defined by:

v = p -y. ( 12
)
where p = (p 1 , . . . , p n ) ⊤ ∈ R n is a vector with all n predictions p i = f (θ, x i ) and y = (y 1 , . . . , y n ) ⊤ ∈ R n . As demonstrated by Du et al. [START_REF] Du | Gradient descent provably optimizes over-parameterized neural networks[END_REF], the continuous dynamics of the error v can be written in a compact way:

d dt v(t) = -H[θ(t)]v(t), v(0) = v 0 ( 13
)
where H[θ] : R ≥0 → R n×n is the NTK matrix from a kernel associated with the ReLU function, being symmetric positive semi-definite as follows:

H ij [θ] = ∂f (θ, x i ) ∂θ , ∂f (θ, x j ) ∂θ
The discrete version resulting from the Euler discretization of ( 13), corresponding to the GD algorithm of ṽ, reads:

ṽk+1 -ṽk = -hH[ θk ]ṽ k ( 14
)
where h is the step size, and ṽ(0) = ṽ0 . We are now ready to formalize the problem:

Problem 1. Given the discrete time system in [START_REF] Kaheman | Learning discrepancy models from experimental data[END_REF], provide conditions over the matrix H[θ] : R ≥0 → R n×n , the loss function L : R d → R ≥0 and the step size h to ensure the convergence of ṽk to zero, together with an explicit bound on its convergence rate.

In the literature, Jerray et al. [START_REF] Jerray | Using euler's method to prove the convergence of neural networks[END_REF] provided a solution to Problem 1 under the following assumptions:

Assumption 1. The gradient descent algorithm used for updating the weights of the neural network in (8) converges to a local minima θ * , i.e. ∂L( θk ) ∂θ converges to 0 as k goes to the infinity.

Assumption 2. There exist λ * > 0 and t 0 ≥ 0 such that, for all t ≥ t 0 : λ min (H[θ(t)]) ≥ λ * , where the time-varying matrix H[θ](t) is given in [START_REF] Jerray | Using euler's method to prove the convergence of neural networks[END_REF].

Assumption 3. L is locally strongly convex around every local minimizer θ * of L, that is: for every local minimizer θ * , there is a neighborhood around θ * on which L is strongly convex. 1 These assumptions are natural in the context of overparameterized NNs, as discussed in [START_REF] Jerray | Using euler's method to prove the convergence of neural networks[END_REF], leading to: Theorem 1. Under Assumptions 1 and 3, if the step size h satisfies h < 2 L , where L is the Lipschitz constant of the loss function L, then the sequence ∥ θkθ k ∥, k ∈ N converges to 0, where θk is defined by (8) and θ k = [w ⊤ 1 (kh), . . . , w ⊤ m (kh), a 1 (kh), . . . , a m (kh)] with θ(t) defined by [START_REF] Freeman | Topology and geometry of half-rectified network optimization[END_REF].

Besides the above theorem, Jerray et al. [START_REF] Jerray | Using euler's method to prove the convergence of neural networks[END_REF] used Assumption 2 to derive the following theorem:

1 Remark: The local strong convexity property is equivalent to the fact that -∂L( θk ) ∂θ belongs to a contractive region, i.e, there exists k 0 ∈ N, a convex region D and a positive real γ > 0 such that θk ∈ D for all k ≥ k 0 , θ ∈ D and θ ′ ∈ D the following relation is true:

∂L(θ) ∂θ - ∂L(θ ′ ) ∂θ , θ -θ ′ ≥ γ∥θ -θ ′ ∥ 2
This is also equivalent to the positive definiteness of the Hessian of L in the neighborhood of each local minimizer θ * . Moreover, in view of [START_REF] Wensing | Beyond convexity-contraction and global convergence of gradient descent[END_REF], the strong convexity corresponds to the special case of contracting gradient descent in the identity metric.

Theorem 2. Under Assumptions 1, 2 and 3, if the step size h satisfies h < 2 L , where L is the Lipschitz constant of the loss function L, then there exist µ ∈ (0, 1) and k 0 ∈ N such that, for all k ≥ k 0 : ∥ṽ k ∥ ≤ ∥ṽ k0 ∥µ k-k0 .

In this paper, the Assumption 2 does not hold, and without loss of generality, we suppose that the space

V(t) of eigenvectors of H[θ(t)] decomposes as V 1 (t) ⊕ V 2 (t), where V 1 (resp. V 2 )
is the cluster of eigenvectors of eigenvalues λ i (t) greater than (resp. less than or equal to) some α > 0. Formally Assumption 2 is replaced by: Assumption 4. There exist λ * ∈ R >0 , t 0 ≥ 0, I 1 = {0, . . . , n 1 -1} and I 2 = {n 1 , . . . , n 1 + n 2 -1} with n 1 + n 2 = n such that the eigenvalues can be divided into:

• λ i (t) ≥ λ * for all i ∈ I 1 , t ≥ t 0 , • λ i (t) = 0 for all i ∈ I 2 , t ≥ t 0 .
where {λ i (t)} i∈I1⊎I2 is the eigenvalues set of the NTK matrix H[θ(t)], and I 1 (resp. I 2 ) the index set of eigenvectors spanning V 1 (t) (resp. V 2 (t)).

Under Assumption 4 instead of Assumption 2, the result of Theorem 2 becomes: Theorem 3. Under Assumptions 1, 3 and 4, if the step size h satisfies h < 2 L and h < 2 λ * , where L is the Lipschitz constant of the loss function L, then there exists k 0 such that for all k ≥ k 0 :

∥ṽ k ∥ ≤ (1 - 1 2 λ * h) 2(k-k0) ∥ṽ 1 k0 ∥ 2 + ∥ṽ 2 k0 ∥ 2 (15)
where ṽi k (i = 1, 2) is the projection of the vector ṽk on the span of the n i eigenvectors of H[ θ(k)]. The right-hand side of (15) thus converges to the constant b = ∥ṽ 2 k0 ∥.

Proof. By Theorem 1, which relies on Assumptions 1 and 3, but not Assumption 2, we know that for h < 2 L , ∥ θkθ k ∥ converges to 0 as k → ∞. Let λi (kh) (i ∈ [n]) denote the eigenvalues of the matrix H[ θk ] obtained from GD, which is discrete in time. Recall that Assumption 4 specifies that the eigenvalues λ i (t) of H[θ(t)] with i ∈ [n 1 ] are in [λ * , ∞), and those with i ∈ [n 2 ] are null. It follows by continuity, using the fact that H[ θk ] is positive semi-definite, that there exists k 0 such that for all k ≥ k 0 :

λi (kh) ∈ [ λ * 2 , ∞) for all i ∈ I 1 , (16) λi
(kh) ∈ [0, λ * 2 ) for all i ∈ I 2 . ( 17
)
We know that for k ≥ k 0 , V 1 k and V 2 k are orthogonal (since all the eigenvalues λi (kh) with i ∈ I 1 are different from the eigenvalues λj (kh) with j ∈ I 2 ). The discretized dynamics of the error in [START_REF] Kaheman | Learning discrepancy models from experimental data[END_REF] then writes: ṽk+1 = (I n -hH[ θk ])ṽ k .

We also know that the NTK matrix H[ θk ] can be decomposed as:

H[ θk ] = P k D k P ⊤ k ( 19
)
where P k is the transition matrix whose columns are the eigenvectors of H[ θk ], and D k is the diagonal eigenvalue matrix with first n 1 elements in [ λ * 2 , ∞), and last n 2 elements in [0, λ * 2 ). From [START_REF] Lucia | Rapid development of modular and sustainable nonlinear model predictive control solutions[END_REF], it follows:

P ⊤ k ṽk+1 = (I n -hD k )P ⊤ k ṽk (20)
The

k+1 = (I n1 -hD 1 k )ṽ 1 k (21) ṽ2 k+1 = (I n2 -hD 2 k )ṽ 2 k ( 22
)
where D 1 k is the n 1 × n 1 top left submatrix of D k , and D 2 k the n 2 × n 2 bottom right submatrix. For k ≥ k 0 and h < 2 λ * , equations ( 21) and ( 22) lead to:

∥ṽ 1 k+1 ∥ ≤ ∥I n1 -hD 1 k ∥∥ṽ 1 k ∥ ≤ (1 - 1 2 λ * h)∥ṽ 1 k ∥ (23) ∥ṽ 2 k+1 ∥ ≤ ∥I n2 -hD 2 k ∥∥ṽ 2 k ∥ ≤ ∥ṽ 2 k ∥ (24)
It follows from ( 23) and ( 24) that there exists k 0 such that for all k ≥ k 0 and h < 2 λ * :

∥ṽ 1 k ∥ ≤ (1 - 1 2 λ * h) k-k0 ∥ṽ 1 k0 ∥ ∥ṽ 2 k ∥ ≤ ∥ṽ 2 k0 ∥
Finally, using the fact that V 1 k and V 2 k are orthogonal:

∥ṽ k ∥ 2 = ∥ṽ 1 k ∥ 2 + ∥ṽ 2 k ∥ 2 ≤ (1 - 1 2 λ * h) 2(k-k0) ∥ṽ 1 k0 ∥ 2 + ∥ṽ 2 k0 ∥ 2 (25)
Remark 2. In practice k 0 is small, which means that ∥ṽ 2 k ∥ is rapidly constant from t 0 = k 0 h. The training error ∥ṽ k ∥ is thus asymptotically equal to b = ∥ṽ 2 k0 ∥. The constant b depends on the initial value θ(0) of the NN parameters. Note that, apart from the choice of initial weight θ(0), the GD procedure considered here is deterministic.

Numerical Results

In order to validate the above analysis, a toy example using the Van der Pol oscillator is used. The same NN architecture presented in Figure 1 is used, but this time only 1 batch of 20 data points are employed for training using the SGD (Stochastic Gradient Descent) optimization algorithm with a learning rate of 10 -2 . The goal is to demonstrate that the error can be decomposed into two subsets (V 1 and V 2 ), the projection of the error onto V 1 converges to 0 as t → ∞, and the projection of the error onto V 2 remains nearly constant during training.

Figure 4 illustrates the evolution of the eigenvalues during gradient descent, in which λ * = 10 -3 is the threshold separating both space of eigenvalues: λ i > λ * . Furthermore, it can be noticed that the eigenvalues evolve in time as the network is underparameterized, which makes the derivation of convergence bounds a difficult task. However, from the projection of the error on V 2 , a lower bound of convergence can be derived as depicted in Figure 5.

From the evolution of ||v 2 ||, it is possible to verify that the lower bound is not zero, remaining nearly constant (b = 0.04) during training. We can also verify some small oscillations of ||v 2 || due to the evolution of eigenvectors V 2 during training as their associated eigenvalues. Despite the fact that the network error v is unable to converge to zero, the lower bound b is small, leading to good performances.

Conclusion

In this paper we were interested in reducing the computational burden of MPC for the realtime synthesis of a control law. The method consists in training a neural network on the data generated by a MPC controller, then using this neural network to control the system in real time. The theoretical guarantees of convergence are derived by expanding the analysis of neural tangent kernels to underparameterized networks. From the results in Section 2, we conclude that data-driven MPC is a viable solution for simulating the behavior of a traditional MPc controller. It shifts the computational burden of MPC to offline training, speeding up the computation of a control law. Regarding performance, the traditional MPC is capable of following the reference more accurately than the data-driven approach, however this gap can be reduced by increasing the neural network complexity.

One of the main focuses of future research is the integration of hard constraints to the neural network, which is not evident as just adding a penalization factor to the cost function as it is done for soft constraints. Furthermore, other applications will be studied as the forging process, which is highly non-linear, and, finally, a proof-of-concept will be carried out in the context of Structural Health Monitoring (SHM) applications.

Regarding the theoretical aspects of network convergence, future work consists in deriving an upper bound of convergence. Moreover, this analysis can be extended to derive guarantees of convergence using Federative Learning for underparameterized neural networks. The interest in this particular type of learning is that training can be performed separately over different datasets, opening the opportunity for parallelization to improve computational time.

Figure 1 :

 1 Figure 1: Architecture of the neural network used to mimic the behavior of an optimization problem.

Figure 2 :

 2 Figure 2: Comparison between the true command (u MPC ) and the command predicted by the neural network (u NN ) using the test dataset.

Figure 3 :

 3 Figure 3: Comparison between traditional and data-driven MPC on following the reference x ref 1

  ) where x ∈ R d represents the input, W = [w 1 , . . . , w m ] ∈ R d×m the weight matrix of the first layer composed by vectors w r ∈ R d for r ∈ [m], and a r ∈ R the output weight. Additionally, θ = [w 1 ⊤ , . . . , w m ⊤ , a 1 ⊤ , . . . , a m ⊤ ] denotes a vector containing the parameters from both layers.

10 Figure 4 :

 104 Figure 4: Evolution of the 11 positive eigenvalues ( λ0 , . . . , λ10 ) of the NTK matrix H[ θk ](the 9 other eigenvalues being approximately 0).

Figure 5 :

 5 Figure5: Log-scale evolution of the training error ∥ṽ k ∥ (blue) and it's projection over V 2 resulting in ∥ṽ 2 k ∥ (orange), which is a lower bound for the error.

Table 1 :

 1 Metrics for the closed-loop simulation.

	Strategy RMSE Computation Time (ms)
	MPC	0.252	2.420
	NN	0.312	0.045

  n×n matrix P ⊤ k is made of an upper n 1 ×n matrix P ⊤ k,1 and a lower n 2 ×n matrix P ⊤ k,2 . The rows of the matrix P ⊤ k,1 are the n 1 eigenvectors of H[ θk ] associated with eigenvalues in [ λ * 2 , ∞], and those of P ⊤ k,2 are the n 2 other eigenvectors. Hence, for a vector v, P ⊤ k,i v (i = 1, 2) corresponds to the projection ṽi k of v on the span V i k of the n i eigenvectors of H[ θk ]. Equation (20) thus decomposes as: ṽ1

Acknowledgment

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program under grant agreement No. 101002857.

* This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program under grant agreement No. 101002857.


Juhyun Park
email: juhyun.park@ensiie.fr

Jeongyoun Ahn

Sparse Functional Linear Discriminant Analysis

Keywords: Sparse functional linear discriminant analysis Domain selection, Functional classification, Functional sparsity, Interpretability, L 1 penalty, Linear discriminant analysis

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

Introduction

We consider a classification problem when data are curves. Denote the functional predictor by X ∈ L 2 (I), observable on a real interval I, and the class label by Y . Assuming that X is a member of two possible groups (Y = 0 or Y = 1), we seek a classification rule that depends on a linear map, with unknown function β,

F β (X) = I X(t)β(t) dt.
This defines a functional linear classification problem, where we wish to determine β(•) in such a way that the linear map yields a good class separation. As functional data models are inherently infinite-dimensional, dimension reduction techniques are essential in constructing a solution.

The standard logistic regression framework for classification can be extended with functional variables [START_REF] Müller | Functional modelling and classification of longitudinal data[END_REF] as

log pr(Y = 1 | X) pr(Y = 0 | X) = a 0 + I X(t)β(t) dt .
Alternatively, one can try to directly extend the Bayes classifier. If X were finitedimensional, the Bayes rule is defined as maximizing the conditional probability given by

pr(Y = 1 | X) = π 1 f 1 (X) π 0 f 0 (X) + π 1 f 1 (X)
,

where i = 0, 1, π i is the probability of X coming from group i and f i (X) is the marginal density of X if it belongs to group i. An optimal classification can then be achieved when one classifies new variable X to the class 1 if π 1 f 1 (X) > π 0 f 0 (X). Consequently, the main effort is devoted to estimating f 1 (X) and f 0 (X). If X is Gaussian, further simplifications can be made and such an approach is broadly known as linear discriminant analysis. The main difficulty with extending the Bayes rule to functional data is linked to the fact that the marginal densities do not exist for functional data [START_REF] Delaigle | Defining probability density for a distribution of random functions[END_REF]. Nevertheless, approaches attempting to directly approximate the marginal densities with dimension reduction techniques have proven useful, even in the absence of well-defined target densities. For example, [START_REF] James | Functional linear discriminant analysis for irregularly sampled curves[END_REF] use a Gaussian framework to develop regularization methods, while [START_REF] Bongiorno | Classification methods for Hilbert data based on surrogate density[END_REF] and [START_REF] Dai | Optimal Bayes classifiers for functional data and density ratios[END_REF] develop nonparametric approaches without Gaussian assumptions.

It is well known that functional classification can achieve a perfect classification, if the infinite-dimensionality is well exploited. This implies that for the purpose of classification, it is not necessarily advantageous to have a well-defined finite-dimensional representation. [START_REF] Delaigle | Achieving near perfect classification for functional data[END_REF] demonstrate such a phenomenon with a simple linear centroid classifier using asymptotic analysis and suggest a practical representation using components obtained from functional principal component analysis and partial least squares. [START_REF] Kraus | Classification of functional fragments by regularized linear classifiers with domain selection[END_REF] propose L 2 regularization methods to obtain the representation. [START_REF] Berrendero | On the use of reproducing kernel Hilbert spaces in functional classification[END_REF] further clarify this phenomenon under a reproducing kernel Hilbert space framework for Gaussian processes, suggesting an alternative finite dimensional approximation.

While optimal performance is an important criterion to consider, the increasing impact of statistical analysis on modern scientific investigations has created the need to carefully consider interpretability of the outcomes of the analysis. Some attempts have been made to address interpretability in functional data, based on the idea that a simpler form of function is easier to interpret, and thus more useful in practice. As the function is an infinite-dimensional object, the formulation is often given in terms of a basis function representation. Under this setting, three different approaches have been proposed to construct a simpler form of functions. The first one is to impose a constraint on the coefficients directly with an 1 norm (e.g., [START_REF] Zhou | Wavelet-based LASSO in functional linear regression[END_REF]. Assuming that β can be well approximated by a finite number of basis functions, say K, this can be expressed as

β(t) = K k=1 α k B k (t) subject to K k=1 |α j | ≤ C
and thus encourages the coefficients to be zero. The second approach is to limit the class of the functions {B k } in terms of their shape such as constant or strictly linear functions only (e.g., [START_REF] Tian | Interpretable dimension reduction for classifying functional data[END_REF]. A difficulty with a standard sparse regularization for a function is that the penalty that encourages a sparse representation of the function does not necessarily have a control over domain selection: i.e., even if α k = 0 for some k, β(t) = 0 for t ∈ I k (t) = {t ∈ I : B k (t) = 0}. The third one is to limit the support of the function β to include zero regions [START_REF] James | Functional linear regression that's interpretable[END_REF][START_REF] Zhou | Functional linear model with zero-value coefficient function at sub-regions[END_REF][START_REF] Martin-Barragan | Interpretable support vetor machines for functional data[END_REF][START_REF] Lin | Interpretable functional principal component analysis[END_REF][START_REF] Picheny | Interpretable sparse sir for functional data[END_REF]. We wish to incorporate interpretability in the latter notion of obtaining zero regions in the solution. However, as noted by [START_REF] Kneip | Functional linear regression with points of impact[END_REF] and [START_REF] Roche | Variable selection and estimation in multivariate functional linear regression via the lasso[END_REF], the theoretical framework appropriate to deal with a discrete notion of sparsity in high dimensions does not necessarily offer an insight into a problem in an infinite dimensional setting.

A functional formulation on sparsity is relatively scarce. [START_REF] Wang | Functional sparsity: Global versus local[END_REF] introduce the notion of functional sparsity, distinguishing global sparsity, which relates to functional variable selection, from local sparsity, which relates to domain selection with zero regions. [START_REF] Tu | Estimation of functional sparsity in nonparametric varying coefficient models for longitudinal data analysis[END_REF] develop a regularization method to achieve simultaneous estimation of both types of sparsity in a time varying functional regression setting and [START_REF] Lin | Locally sparse estimator for functional linear regression models[END_REF] propose an alternative regularization to achieve local sparsity in functional linear regression. Both approaches rely on a clever grouping of the sparse coefficients in the basis function representation. [START_REF] Hall | Truncated linear models for functional data[END_REF] study the issue of identifiability of domain selection problem in functional linear regression and suggest a domain search strategy. [START_REF] Kraus | Classification of functional fragments by regularized linear classifiers with domain selection[END_REF] follow a similar line.

In this work, we seek an alternative approach to functional linear classification with a direct estimation method that addresses dimensionality, optimality and interpretability. With the need to embed additional constraints on the form of the solution, we reformulate the functional linear discriminant analysis as a regularization problem with an appropriate choice of penalty functions. So far, the penalty-based regularization methods have been used mostly for either smoothness in regression (e.g., [START_REF] Cardot | Spline estimators for the functional linear model[END_REF][START_REF] Crambes | Smoothing spline estimators for functional linear regression[END_REF] or invertibility in classification [START_REF] Kraus | Classification of functional fragments by regularized linear classifiers with domain selection[END_REF], based on an L 2 -type penalty. The idea of sparsity as variable selection with an 1 -type penalty is actively developed in the high-dimensional setting with scalar variables, but much less for the infinite-dimensional functional setting. We develop a new regularization method for functional linear discriminant analysis with an L 1 -type penalty in the functional setting for inducing zero regions, as opposed to zero coefficients. We study the underlying optimization problem in detail, showing that it has a well defined solution that contains zero regions, achieving a functional sparsity in the sense of domain selection. Furthermore, we show that our regularized classifier asymptotically behaves similarly to the optimal classifier. We do not presume, even in the case where the optimal classifier is well defined, that the underlying projection function β has a local sparse property with true zero regions in the domain. The role of the L 1 penalty is to provide a sparse estimator that consistently estimates the true function whether or not the latter is sparse, enhancing interpretability without compromising prediction performance. Unlike the approaches based on direct dimension reduction techniques such as principal component analysis or partial least squares, our proposed regularization method does not require assumptions on the eigenvalue sequences of the covariance function or any other unknown quantities. Numerical comparisons with existing methods demonstrate the properties of the proposed estimator in finite samples with both simulated and real data examples.

Methodology

Functional data framework

We use the following standard notation: for f ∈ L p (I), 1 ≤ p < ∞, the norm is defined by

f p = ( I |f (t)| p dt) 1/p . When p = ∞, f ∞ = inf{C : |f (t)| ≤ C a.e. on I}.
Assume that X ∈ L 2 (I) with the mean µ(t) and covariance function cov{X(s), X(t)} = γ(s, t), s, t, ∈ I. Define the inner product in L 2 (I) by f, g = I f g for f, g ∈ L 2 (I). Assuming I I γ(s, t) 2 dsdt < ∞, define the covariance operator Γ : L 2 (I) → L 2 (I) as

Γ(β)(t) = I γ(s, t)β(s) ds , t ∈ I.
It is known that Γ is a compact operator, and is Hilbert-Schmidt, that admits a spectral representation given by

Γ(β) = ∞ j=1 θ j β, φ j φ j , θ j ≥ 0 , θ j → 0 as j → ∞ ,
where θ j and φ j (•) correspond to eigenvalues and eigenfunctions of the covariance operator Γ, respectively. One of the important properties of a compact operator is that it is not invertible unless it has only finitely many distinct eigenvalues. At the same time, since θ j → 0 as j → ∞, a finite rank approximation to Γ is also well understood, which has been the basis of many regularization and dimension reduction techniques for functional data, notably functional principal component analysis and its applications (e.g., [START_REF] Hall | A functional data-analytic approach to signal discrimination[END_REF]. More details on theoretical foundations of functional data are found in [START_REF] Hsing | Theoretical Foundations of Functional Data Analysis, with an Introduction to Linear Operators[END_REF].

For our classification problem, we denote the variable in each group by X 0 and X 1 and we make the following standard assumptions.

Assumption 1. For k = 0, 1, X k is square integrable function on a compact interval I with µ k = E(X k ) and µ 0 = µ 1 with common covariance function γ. The mean functions and the covariance function are continuous.

Assumption 2. For k = 0, 1, E( X k 4 2 ) < ∞.
Without loss of generality, we assume that µ 0 = 0 and µ 1 = δ.

An optimal linear classifier

The optimality of the linear methods in the homoskedastic Gaussian scenario is well established by [START_REF] Delaigle | Achieving near perfect classification for functional data[END_REF] based on an asymptotic centroid-based classifier. We briefly review their framework to motivate our regularization method, which is introduced in the following section.

For given X, the linear classifier with β can be defined in multiple ways. We assume that β ∈ L 2 (I) is continuous. The class assignment based on β is done via

T 0 (X) = X, β -δ, β 2 -X, β 2 , which assigns to X the group label Y = 1 if T 0 (X) < 0 and Y = 0 if T 0 (X) > 0. Then, the misclassification error is π 0 pr{T 0 (X) < 0 | Y = 0} + π 1 pr{T 0 (X) > 0 | Y = 1}. If X
is Gaussian, the misclassification probability can be expressed as

err X (β) = 1 -Φ | δ, β | 2 β, Γβ 1/2 , (1)
with its minimal error given by 1 -Φ( Γ -1/2 δ 2 /2) [START_REF] Delaigle | Achieving near perfect classification for functional data[END_REF][START_REF] Kraus | Classification of functional fragments by regularized linear classifiers with domain selection[END_REF]. It can be seen that an optimal function β that minimizes the misclassification error is

max β =0 δ, β 2 Γβ, β = max β =0 δ, β 2 var( X, β ) .
This is equivalent to the extended criterion for Fisher's discrimination analysis for functional data [START_REF] Shin | An extension of Fisher's discriminant analysis for stochastic processes[END_REF] defined by max

β =0 var[E{F β (X) | Y }] E[var{F β (X) | Y }] = max β =0 π 0 π 1 δ, β 2 Γβ, β . (2)
Attempting to directly solve (2) would lead to a form of generalized eigenvalue problem.

Equivalently, this can be formulated as

max β δ, β 2 subject to Γβ, β -1 = 0. ( 3
)
The solution can be derived from the Lagrangian formulation of (3)

J ρ (β) = δ, β 2 -ρ{ Γβ, β -1} ,
for some constant ρ, which leads to

Γβ = δ . (4)
Hence, we see that the optimal linear classifier can be defined as the solution to (4). However, since a bounded inverse of Γ does not exist in the infinite dimensional setting (e.g., [START_REF] Cardot | CLT in functional linear regression models[END_REF], this equation clearly illustrates that linear discrimination analysis for functional data is an ill-posed inverse problem. When the covariance operator further satisfies Γ -1 δ 2 < ∞, a unique classifier exists in (Ker(Γ)) ⊥ with the optimal error 1 -Φ( Γ -1/2 δ 2 /2). When Γ -1 δ 2 = ∞, an optimal solution does not exist, but optimal classification can be achieved asymptotically along a non-convergent path and perfect classification may be possible if Γ -1/2 δ 2 = ∞. Therefore, unlike regression, it is not necessary to assume the existence of a unique solution to obtain an approximate solution, whose performance can mimic the optimal classifier asymptotically [START_REF] Delaigle | Achieving near perfect classification for functional data[END_REF][START_REF] Kraus | Classification of functional fragments by regularized linear classifiers with domain selection[END_REF]. In the following, we will denote the solution to (4) by β 0 when it exists and is finite.

A regularized solution to discrimination

The primary purpose of regularization is to solve an ill-posed inverse problem. We introduce our regularization method to solve the ill-posed inverse problem (4) above. A sensible strategy would be to impose some constraints on the solution set or a penalty term to an objective function. In order to impose the functional equation ( 4), we introduce a corresponding minimization problem. Specifically, viewing the discriminant equation in (4) as a type of score equation leads to defining an objective function as

(β) = 1 2 Γβ, β -δ, β . (5)
Hence the functional derivative of ( 5) with respect to β yields (4). In practice, where Γ and δ are not available, replacing (Γ, δ) by their empirical counterpart (Γ n , δ n ) gives

n (β) = 1 2 Γ n β, β -δ n , β .
For example, the standard sample covariance operator and the sample mean could be used for Γ n and δ n . The standard sample estimators often lack smoothness (e.g., [START_REF] Cardot | Functional linear model[END_REF] so it is often desirable to use their smoothed version by taking into account the underlying design schemes [START_REF] Zhang | From sparse to dense functional data and beyond[END_REF]. Our formulation does not rely on the specific choice of estimators as long as they are consistent. To incorporate additional constraints, we consider the following optimization problem

min β∈H { n (β) + P (β)} , ( 6
)
where H is an appropriate function space and P (β) is a penalty corresponding to the constraint. The type of constraint requires some further consideration. As the underlying problem is defined in L 2 (I), it seems natural to add a penalty that depends on the L 2 norm of β such as β 2 2 as in [START_REF] Kraus | Classification of functional fragments by regularized linear classifiers with domain selection[END_REF], while a smoothness constraint on β would lead to the L 2 penalty on the derivatives. Nevertheless, L 2 regularization in general does not produce a sparse solution (e.g., [START_REF] Lin | Locally sparse estimator for functional linear regression models[END_REF][START_REF] Tu | Estimation of functional sparsity in nonparametric varying coefficient models for longitudinal data analysis[END_REF]. This phenomenon is well investigated in multivariate data, and the popularity of the 1 or 0 penalty demonstrates the effectiveness of a sparse solution in a broader context.

In order to introduce a type of sparsity in the discriminant function, we propose to directly impose the functional norm constraints β 2 2 ≤ C 1 and β 1 ≤ C 2 . Although other constraints are possible, it turns out that this combination of constraints allows direct control over the L 2 norm. One might wonder whether the derivative penalty is necessary in our sparse regularization. Our experiences suggest that a combined norm give more stable results. We have assumed compact support to simplify the statement in Assumption 1. Under the assumption of compact support for β, the role of the derivative penalty does not seem so significant with respect to the classification performance. In general, since we only have β 1 ≤ β 2 , the bound on the L 1 norm alone is not sufficient to ensure the L 2 properties of β and the techniques we have developed below do not necessarily work. On the other hand, we have (e.g., [START_REF] Gabushin | Inequalities for the norms of a function and its derivatives in metric l p[END_REF][START_REF] Li | L 1 regularization for compact support[END_REF] that

β 2 ≤ β 1 + β 2 ,
which suggests that the proposed method is able to control the L 2 norm of the function more effectively by regularizing both the derivative and the L 1 norm of the function.

Taking into account the functional constraints in (6) leads to the following objective function

J n (β) = 1 2 Γ n β, β -δ n , β + P (β) , (7)
where P (β) = λ β 1 + (η/2) β 2 2 and λ and η are tuning parameters, chosen from data. In this formulation, the derivative penalty is a standard smoothness penalty. The non-standard part is the L 1 penalty, which is an infinite-dimensional counterpart of the 1 sparsity penalty. As is demonstrated below, the L 1 norm can be linked to a sparse solution in function space in terms of domain selection [START_REF] Li | L 1 regularization for compact support[END_REF], which is a special case of functional sparsity [START_REF] Wang | Functional sparsity: Global versus local[END_REF]. Arguably, the result would be easier to interpret as the localized effects are automatically identified.

We note that it is possible to use the ridge penalty ( β 2 2 ) instead of the derivative norm. In fact, this choice makes the analysis of the optimization problem much easier as it directly controls the L 2 norm. In this case, the fact that

Γ n β, β + τ β 2 2 ≥ τ β 2 2
for all τ > 0 is sufficient to guarantee the existence of a solution. Nevertheless, we have chosen the derivative norm that has the added consideration of smoothness of the solution, since a smoother solution is easier to interpret.

Properties of the regularized classifier in the population model

In this section we study the characteristics of the proposed regularized solution to the infinite-dimensional convex optimization problem. In order to simplify the discussion, we first consider the following optimization problem that replaces the sample quantities by the population counterparts. Let

J(β) = 1 2 Γβ, β -δ, β + λ β 1 + η 2 β 2 2 . ( 8
)
The smoothness constraint with β reduces the feasible set of β from L 2 (I) to a differentiable subspace [START_REF] Brezis | Functional analysis, Sobolev spaces and partial differential equations[END_REF]. Since L 2 (I) ⊂ L 1 (I), we consider H = {β ∈ H 1 (I) : β 1 + β 2 < ∞} as the space of feasible solutions, equipped with the norm given by β = β 1 + β 2 . Then it is clear that H is a convex set.

H 1 (I) = {β ∈ L 2 (I) : β is absolutely continuous, β ∈ L 2 (I)}, which is the Sobolev space H 1 = W 1,2 in standard notation
Under this setting, it can be shown that J is convex, so that the existence of a solution is guaranteed. Moreover, we can show that J is strictly convex so the solution is unique. We summarize the result in the following proposition.

Proposition 1. The function J : H → R defined in (8) has a global minimizer β ∈ H such that J( β) ≤ J(β)
for all β ∈ H. Furthermore, the solution is unique.

In order to understand the property of the solution β, additional characterizations of the solution beyond its existence are necessary. Note that if the objective function is differentiable, the first derivative (along with the second derivative) condition sufficiently characterizes the solution and a Newton-type algorithm can be used to find it. Nondifferentiable functions require an alternative form to replace the derivative condition. Following similar arguments in Reyes (2015, ch 6) and Glowinski (1984, p. 70-71), we derive a set of conditions that β satisfies in the following proposition, which could be used for developing an optimization algorithm.

Proposition 2. If β is the minimizer of J(β) in ( 8), there exists a function α ∈ L 2 (I) with |α(t)| ≤ λ a.e. for t ∈ I such that β satisfies the following relation:

Γ β, β + η β , β -δ, β + α, β = 0 , for all β ∈ H , (9)
and

   α(t) = λ on {t ∈ I : β(t) > 0} , |α(t)| ≤ λ on {t ∈ I : β(t) = 0} , α(t) = -λ on {t ∈ I : β(t) < 0} . ( 10
)
The equations ( 9) and ( 10) are necessary conditions for the solution to satisfy. Equation (9) means that the subgradient of the objective at the minimizer β contains zero. Furthermore, Proposition 2 demonstrates the domain selection property of the solution.

In order to make an analogy to the finite high dimensional case, consider a simplistic example where Γ = I, the identity operator, and η = 0 in (8). Note that, since the identity operator is not compact, this is not a realistic example from the functional data perspective but serves as a toy example. Then, the optimal solution β can be given explicitly as

β(t) =    δ(t) -λ if δ(t) > λ , 0 if δ(t) ∈ [-λ, λ] , δ(t) + λ if δ(t) < -λ .
Since δ is continuous, the set {t ∈ I : β(t) = 0} will be a union of intervals and then β joins the zero intervals continuously at the boundary. Hence the solution exhibits a thresholding behaviour, similar to the sparse estimators in the finite dimensional case, but on the continuous domain, thus justifying the notion of functional sparsity in the sense of domain selection. When η = 0, the optimization occurs in L 2 (I), so the boundary points of the zero intervals are not necessarily differentiable. For the general case, it is difficult to express the solution in an analytical form but the same argument holds and when η > 0 as in ( 9) and (10), the function values at zero intervals and non-zero intervals are joined smoothly.

Moreover, suppose that β c ∈ H 1 0 (I) = {u ∈ H 1 (I) : u is zero at the boundary of I}.

Then, (9) leads to

I (Γ β -δ + α)β c + I (η β )β c = 0 , for all β c ∈ H 1 0 (I) . (11)
This implies (e.g., [START_REF] Brezis | Functional analysis, Sobolev spaces and partial differential equations[END_REF]) that β , the derivative of the solution, is also absolutely continuous and satisfies

η β = Γ β -δ + α , a.e. , (12)
which, together with (10), gives an additional necessary condition for the solution.

The following observation helps us to determine the upper bound of the tuning parameter λ. Suppose that δ ∈ L ∞ (I). Then, we have δ, u ≤ δ ∞ u 1 . In addition, it can be shown that the solution β satisfies

Γ β, β + η β , β + λ β 1 = δ, β ≤ δ ∞ β 1 . It follows that Γ β, β + η β , β + (λ -δ ∞ ) β 1 ≤ 0 .
Since the first two terms are non-negative, if δ ∞ ≤ λ, then β = 0. This gives a range of λ values for a non-zero solution and this depends on the magnitude of δ. Figure 1 illustrates the effect of the tuning parameters η and λ on β, under two different scenarios for β 0 : sparse and non-sparse. The curves in the top and bottom panels respectively correspond to the simulation settings 3 and 5 in Section 3. We observe that for fixed λ, larger η gives smoother estimates and for fixed η, larger λ gives a more sparse solution.

Suppose that Ker(Γ) = {0}. If Γ -1 δ 2 < ∞, then β 0 = Γ -1 δ ∈ L 2 (I) is the unique solution to (4). Denote by β the minimizer of J in (8). The following proposition shows that the regularized solution β approximates the underlying model solution to (4) whether the latter is finite or not. We write λ = (λ, η), and λ → 0 means that λ → 0 and η → 0.

Proposition 3. Fix any u ∈ L 2 (I). If Γ -1 δ 2 < ∞, then β, u → β 0 , u as λ → 0. If Γ -1 δ 2 = ∞, then β → ∞, as λ → 0.
Using Proposition 3 together with (1), we can study the misclassification error of the regularized classifier. We show below that our regularized classifier mimics the behaviour of the optimal classifier in terms of the misclassification rate.

Proposition 4. Assume that X is Gaussian with covariance operator Γ. The misclassification probability of the regularized classifier β, denoted by err X ( β), converges to 1 -Φ( Γ -1/2 δ 2 /2) as λ → 0.

Remark 1. It may be suspected that the two classes have different covariance functions, say Γ 0 and Γ 1 , in which case a full generalization of our linear approach is not trivial. For example, a direct generalization of (2) replaces the denominator by Γβ, β where Γ = π 0 Γ 0 +π 1 Γ 1 , the pooled covariance function. Then the same formulation as (8), with Γ replaced by Γ, follows. Alternatively, we can consider the two-dimensional projection approach for quadratic discriminant analysis proposed by [START_REF] Gaynanova | Sparse quadratic classification rules via linear dimension reduction[END_REF] in the multivariate setting. In our functional framework, we would optimize the following two objective functions

J i (β) = 1 2 (Γ i β, β -δ, β + λ i β 1 + η i 2 β 2 2 ,
for i = 0, 1. The classification will be then made based on the two linear maps with the respective solutions. An advantage of this method over standard quadratic discriminant analysis is that the computational complexity is essentially the same as the homoscedastic case being considered here. In either case, the theoretical guarantee on the domain selection and classification accuracy need to be carefully investigated in the functional context.

Properties of the regularized classifier in the sample model

The development in the previous section justifies the use of the L 1 penalty in our formulation, by ensuring that we have a well-defined regularization problem (Proposition 1) with the desired sparsity in the solution (Proposition 2). We can show that these properties also hold for the empirical estimator, denoted by β, the minimizer of J n defined in ( 7).

In the sample case, our formulation of J n is based on two estimators δ n and Γ n . Provided that the sample covariance operator Γ n is non-negative definite and δ n ∈ L 2 (I), it can be shown that Propositions 1 and 2 hold for J n and β. In fact, the proofs are essentially the same, thus we omit the formal statements of the Corollaries here.

In order to study the limit of the empirical estimator β and the corresponding misclassification probability, we make the following assumptions. Assumption 3 says that some consistent estimators are available for δ n and Γ n (or γ n ) and Assumption 4 assumes that the tuning parameters decrease not too rapidly as the sample size gets large.

Assumption 3. Let (δ n , γ n ) be consistent estimators such that

γ n -γ ∞ = O p (a n ), δ n -δ ∞ = O p (b n ) ,
for some sequences a n , b n converging to 0.

Assumption 4. Let λ = λ n → 0, η = η n → 0 with λ n /η n = O(1), a n /η n = o(1) and b n /η n = o(1).
For sufficiently dense data as we consider here, the standard sample mean and sample covariance function can be used with a root-n convergence. Alternatively, any preferred smoothing methods could be employed, then the convergence rates are given as a function of smoothing parameters used, for example, see [START_REF] Yao | Functional data analysis for sparse longitudinal data[END_REF]; [START_REF] Li | Uniform convergence rates for nonparametric regression and principal component analysis in functionallongitudinal data[END_REF] for local linear smoothers. Normally the error for covariance estimation is larger than that for the mean estimation (a n ≥ b n ). [START_REF] Zhang | From sparse to dense functional data and beyond[END_REF] generalize those results under a general design scheme and show that a root-n rate can be achieved for dense and ultra-dense data, in which case the uniform convergence rates can be

a n = b n = (log n/n) 1/2 . Proposition 5. Fix any u ∈ L 2 (I). Suppose that Assumptions 1-4 hold. If Γ -1 δ 2 < ∞, then β, u → β 0 , u in probability as n → ∞. On the other hand, if Γ -1 δ 2 = ∞, then β → ∞ in probability as n → ∞.
Proposition 6. Assume that X is Gaussian with covariance operator Γ. Under Assumptions 1-4, the misclassification probability of the regularized classifier β, denoted by err X ( β), converges to 1 -Φ( Γ -1 δ 2 /2) in probability as n → ∞.

Numerical algorithms

A fundamental strategy to solve infinite-dimensional optimization problems could be said to be either to discretize-and-optimize or to optimize-and-discretize [START_REF] Reyes | Numerical PDE-constrained optimization[END_REF]. For example, [START_REF] Glowinski | Numerical methods for nonlinear variational problems[END_REF] analysed numerical methods based on the former approaches using piecewise linear and quadratic approximations on H 1 (I), while [START_REF] Reyes | Numerical PDE-constrained optimization[END_REF] includes examples of the latter. Generally, the former is easier, while the latter may give a more elegant solution. According to Proposition 2, the second approach requires solving equations ( 9), or (12), and (10) simultaneously. Based on these equations, Newton-like update steps for both β and α could be derived. However, as this involves further development of the testable conditions at the iteration to be adapted to a proper function space [START_REF] Reyes | Numerical PDE-constrained optimization[END_REF], we leave it for future work.

Here we have used piecewise linear approximations in our implementation with equidistant grid points. We use estimates from the method by [START_REF] Yao | Functional data analysis for sparse longitudinal data[END_REF] for Γ n and δ n and evaluate the derivative by the finite-difference approximation. This strategy, combined with the rather simple form of the optimization problem, leads to a lasso-type problem, for which many efficient optimization algorithms are available. We implement the coordinate descent algorithm by [START_REF] Fu | Penalized regressions: The bridge versus the lasso[END_REF] to solve our discretized problem.

Implications in functional linear regression

Although we started with a classification problem, our regularization method with an L 1 -type penalty can be motivated from a regression point of view. Consider a functional regression problem with scalar response Y and functional predictor X ∈ H = L 2 (I). For simplicity, assume that E(Y ) = 0 and E(X) = 0 and consider a linear regression given by

Y = I β(t)X(t) dt + ,
where E( ) = 0 and var( ) = σ 2 . As before, the covariance operator of X is denoted by Γ. In addition, let Λ : L 2 (I) → R denote the covariance operator between X and Y , defined as

Λ(u) = I E{X(s)Y }u(s) ds .
Then, the population least squares criterion can be expressed as

E{(Y -β, X ) 2 } = Γβ, β -2Λ(β) + const ,
from which it follows that if β is the minimizer if and only if it satisfies

Γ β, u = Λ(u) for all u ∈ H . (13)
To be consistent with the earlier notation, write Λ(u) = ∆, u where ∆(

•) = E{X(•)Y }.
Then, (13) can be expressed as

Γ β, u = ∆, u for all u ∈ H .
Note that when Y is binary with values in {-1, 1}, ∆ is equal to δ in (4), the difference in the group means in X between two classes. The existence of the solution requires the so-called Picard condition [START_REF] Cardot | Spline estimators for the functional linear model[END_REF][START_REF] Hsing | Theoretical Foundations of Functional Data Analysis, with an Introduction to Linear Operators[END_REF]

) ∞ j=1 E{X(•)Y }, φ j 2 θ 2 j = ∞ j=1 ∆, φ j 2 θ 2 j < ∞ , which is equivalent to Γ -1 δ 2 < ∞ in the case of classification problems (c.f., Sec- tion 2.2).
The above discussion shows that even though the motivation differs, the underlying problem for functional linear regression is equivalent to that of linear classification. Hence, our methodology is equally applicable in a regression setting and offers a regularization method to estimate the coefficient function β. Given a sample of size n of (Y i , X i ), i = 1, . . . , n, denote the sample version of the operators of Γ and Λ by Γ n and Λ n , respectively, defined as

Γ n (u)(t) = 1 n n i=1 X i , u X i (t) , Λ n (u) = 1 n n i=1 X i , u Y i = ∆ n , u u ∈ L 2 (I) .
where

∆ n (•) = (1/n) n i=1 X i (•)Y i .
Assuming that the true coefficient function β is a smooth function, [START_REF] Cardot | Spline estimators for the functional linear model[END_REF] proposed a generalization of a ridge regression with an L 2 derivative penalty as

1 2 Γ n β K , β K -∆ n , β K + η 2 β (m) K 2 2 ,
where

β K = K j=1 c k B k with B k is a B-spline basis function and β (m)
K , m ≥ 1 is the mth derivative of β K . Our analysis suggests that adding an L 1 penalty would lead to a sparse estimator that consistently estimates the true function, whether or not the latter is sparse.

Remark 2. Although the theoretical framework differs, our main consideration of interpretability is not limited to this type of problem. For example, semi-parametric problems of single-index or multiple-index regression or classification can be expressed as

Y = F k X, β k +
with a known or unknown link function F . Then an L 1 -type functional regularization for β k can be developed with an appropriate loss function. We leave such a generalization to our future work.

3 Numerical studies

Simulation models

We use simulated curves under six different underlying population scenarios to assess the performance of the proposed method and compare it with existing approaches. Settings 1 -5 represent scenarios where the underlying structure follows model (4) in Section 2.2 and Setting 6 is borrowed from [START_REF] Delaigle | Achieving near perfect classification for functional data[END_REF]. The proposed approach, which we name SFLDA (Sparse Functional Linear Discriminant Analysis), is compared with eight different approaches: the non-sparse version of the proposed approach (FLDA) that sets the tuning parameter λ = 0 in (7); ridge FLDA (RFLDA) proposed by [START_REF] Kraus | Classification of functional fragments by regularized linear classifiers with domain selection[END_REF]; partial least squares (PLS) by [START_REF] Delaigle | Achieving near perfect classification for functional data[END_REF]; the three Bayes methods (B-Gauss, B-NPD, B-NPR) in [START_REF] Dai | Optimal Bayes classifiers for functional data and density ratios[END_REF]; the quadratic discriminant analysis (QDA) by [START_REF] Galeano | The Mahalanobis distance for functional data with applications to classification[END_REF]; functional logistic regression (Logistic) [START_REF] Müller | Functional modelling and classification of longitudinal data[END_REF].

The first four settings have a sparse underlying discriminant function. The sparse discriminant functions, β 1 (t) and β 2 (t), are generated with 33 cubic B-spline basis functions B i,4 (i = 1, . . . , 33) with the knots ξ j = j/30, j = 0, . . . , 30, where ξ 0 and ξ 30 are the boundary knots. We set β 1 (t) = 0.2B 5,4 (t) + 0.2B 28,4 (t) for Settings 1 and 2, and β 2 (t) = 0.2B 5,4 (t) -0.2B 24,4 (t) for Settings 3 and 4. For the non-sparse β 3 (t) for Setting 5, we use 8 cubic B-spline basis and set β 3 (t) = 0.1B 1,4 (t) -0.3B 3,4 (t) -0.2B 5,4 (t) + 0.2B 7,4 (t) -0.1B 8,4 (t).

For the covariance function Γ, we used the Matérn covariance function:

γ(s, t) = σ 2 2 1-ν G(ν) (2ν) 1/2 |s -t| ρ ν K ν (2ν) 1/2 |s -t| ρ ,
where G is the gamma function, K ν is the modified Bessel function of the second kind, for which the R-function besselK was used, and the parameters were set as σ = 1, ρ = 0.2, and ν = 3. We set the mean of the first group as µ 0 (t) = 5t + 5 j=1 (c j /j)φ(t), where φ(t) = 2 1/2 sin(πjt), (c 1 , . . . , c 5 ) = (2.19, -0.18, -0.19, -2.51, -0.56), and the second group mean as µ 1 (t) = µ 0 (t) + δ(t). Here the mean difference function δ(t) is determined by δ(t) = I Γ(t, s)β k (s)ds, k = 1, 2, 3.

The data were assumed to be available on a fine grid. We used T = 100 equispaced grid points on [0, 1]. For grid points t j , j = 1, . . . , T , let Γ be the T × T matrix with (i, j) entry Γ(t i , t j ), and let m 0 and m 1 be the T × 1 vectors with ith entry µ 0 (t i ) and µ 1 (t i ), respectively. Data were generated by x 0 = m 0 + Γ1/2 e 0 and x 1 = m 1 + Γ1/2 e 1 where e 0 and e 1 are the vectors of independent noise for each class, generated from either a standard normal for Settings 1, 3, and 5 or a t-distribution with 5 degrees of freedom for Settings 2 and 4. As mentioned earlier, Setting 6 is borrowed from [START_REF] Delaigle | Achieving near perfect classification for functional data[END_REF]. Each curve from the ith group is generated as

x i = 40 j=1 (j -1 Z ij + µ ij )φ j (t)
, where Z ij are centred exponential variables, (µ 01 , . . . , µ 06 ) = (0, -0.5, 1, -0.5, 1, -0.5), and (µ 11 , . . . , µ 16 ) = (0, -0.75, 0.75, -0.15, 1.4, 0.1).

We generated n = 100 curves for each group to train each classifier and independently generated 300 curves for each group to evaluate the misclassification error. The tuning parameters in each method are chosen via 5-fold cross-validation within the training data based on the misclassification error.

Results

Table 1 lists test errors based on 100 repetitions. The proposed approach is generally competitive in all settings, even though the differences among the methods do not stand out. A main benefit of the proposed sparse estimation can be seen in its ability to estimate zero regions in the discriminant function. Table 2 displays the norm differences of the estimated discriminant function β(t) by the proposed method with sparse functional penalty, its non-sparse version, ridge regularization, and partial least squares, in Settings 1, 3, and 5. The proposed method has the lowest error with respect to both norms in the two sparse settings (1 and 3), while its non-sparse version and partial least squares are better in the non-sparse setting (5). Also shown in Figure 2 are the estimated curves for each of the four methods for Settings 1, 3, and 5, along with the median in solid black and the true discriminant curve in a dot-dashed black line. The median curve of the estimators from the proposed method resembles the true function the most in Settings 1 and 3, successfully identifying the non-signal regions. The estimators with a ridge penalty show larger variability than other methods, which suggests that a ridge-type regularization in functional classification does not necessarily reduce variance, unlike in ridge regression.

Real data examples

We tested the same set of the functional classification methods as in the previous section for real data examples. The Tecator dataset [START_REF] Delaigle | Achieving near perfect classification for functional data[END_REF] has 240 curves of near infrared absorbance spectra (850 -1050 nm) of finely chopped meat, using a Tecator Infratec Food & Feed Analyser. Here the groups are defined according to the fat content. The Wine spectra dataset [START_REF] Dai | Optimal Bayes classifiers for functional data and density ratios[END_REF] contains 123 samples of mid infrared spectra (4000 -400 cm -1 ), which are divided into two groups based on the alcohol content level. The Growth dataset [START_REF] Gasser | Nonparametric regression analysis of growth curves[END_REF][START_REF] Sheehy | An analysis of variance of the pubertal and midgrowth spurts for length and width[END_REF] has height growth curves of 112 boys and 120 girls from births to the 18th year. We also analysed the velocity of the Growth curves. The Wheat data [START_REF] Kalivas | Two data sets of near infrared spectra[END_REF] have near infrared spectra of 100 wheat samples, divided into two groups according to the protein content. The Phoneme data are log-periodograms constructed from digitized speech of two different sounds "aa" and "ao", as described in [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF]. Each panel of Figure 3 respectively shows 40 randomly selected curves from each data set, with each class shown with different colors.

To estimate test errors, we split into 2/3 training and 1/3 testing data, and tune each method in the same way as the simulation study, and repeat this process 30 times. Mean test errors are displayed in Table 3. While there are no meaningful differences among the test errors of these methods, we can see that the proposed method achieves compatible accuracies with a much better ability to find localized discriminant regions, as can be seen in Figure 4. In particular, we have found that the selected regions of Tecator, Wine, and Phoneme are consistent with findings from the spectroscopy literature. Specifically, according to the body composition study [START_REF] Conway | A new approach for the estimation of body composition: infrared interactance[END_REF], fat has high absorbance around 930nm, the absorbance region for ethanol is concentrated in 1200 -850 cm -1 [START_REF] Debebe | Non-destructive determination of ethanol levels in fermented alcoholic beverages using fourier transform mid-infrared spectroscopy[END_REF], and it has been found by [START_REF] Hastie | Penalized discriminant analysis[END_REF] that the discriminating feature in phonemes are in the low frequencies about 500 -1000Hz, corresponding to the frequencies 16-32. All these known regions are consistent with the estimates from the proposed method in the figure.

Figure 1 :

 1 Figure 1: Illustration of the roles of the tuning parameters η and λ. Each panel plots the true β 0 (t), as well as β(t) with a varying tuning parameter.

Figure 2 :

 2 Figure 2: Discriminant functions from 100 repetitions of the simulation study estimated by SFLDA, FLDA, RFLDA, and PLS, for Settings 1, 3, and 5. The thicker black line in each plot represents the median of the estimated curves, while the dot-dashed curve is the true discriminant function.

Figure 3 :

 3 Figure 3: For each real data example, 20 sample curves from each class are shown along with average curves. Solid and dotted curves represent curves from different classes.

Figure 4 :

 4 Figure 4: Estimated discriminant curves for real data examples from 30 repetitions. The median curve is shown as a thicker black solid line.

Table 1 :

 1 Test error rates (%) for simulated data based on 100 repetitions. Standard errors are in the parentheses.

	Setting SFLDA FLDA RFLDA	PLS	B-Gauss B-NPD B-NPR Logistic QDA
	1	33.3 (0.22)	33.4 (0.22)	33.6 (0.24)	33.7 (0.23)	33.4 (0.24)	34.6 (0.25)	35.1 (0.25)	34.7 (0.27)	34.2 (0.28)
	2	33.5 (0.21)	33.5 (0.20)	33.8 (0.20)	33.7 (0.20)	33.6 (0.20)	34.7 (0.22)	35.4 (0.24)	34.9 (0.25)	34.2 (0.22)
	3	37.2 (0.23)	36.9 (0.20)	37.5 (0.20)	37.1 (0.21)	37.2 (0.22)	37.8 (0.25)	39.0 (0.29)	38.5 (0.27)	38.0 (0.25)
	4	37.7 (0.24)	37.3 (0.22)	38.1 (0.24)	37.6 (0.23)	37.5 (0.24)	39.1 (0.27)	40.3 (0.30)	39.2 (0.28)	38.2 (0.26)
	5	3.2 (0.08)	3.1 (0.07)	3.1 (0.07)	3.1 (0.07)	3.2 (0.08)	3.5 (0.10)	3.7 (0.10)	3.6 (0.14)	4.2 (0.08)
	6	3.3 (0.08)	3.2 (0.07)	3.3 (0.07)	3.1 (0.07)	4.0 (0.08)	3.4 (0.08)	2.5 (0.10)	2.5 (0.10)	3.6 (0.14)

Table 2 :

 2 Average norm difference ( β -β 0 j ) between the estimated discriminant function and the true function. Numbers are multiplied by 100.

	Setting j	SFLDA	FLDA	RFLDA	PLS
	1	1 5.37 (0.19) 2 9.64 (0.22) 10.07 (0.11) 11.29 (0.18) 10.14 (0.10) 7.97 (0.10) 8.81 (0.13) 8.20 (0.09)
	3	1 5.84 (0.18) 2 10.39 (0.22) 10.40 (0.14) 11.69 (0.18) 10.26 (0.11) 8.10 (0.11) 9.00 (0.12) 8.30 (0.11)
	5	1 6.47 (0.14) 2 8.21 (0.18)	4.15 (0.11) 5.25 (0.13)	5.78 (0.17) 7.23 (0.20)	4.50 (0.10) 5.64 (0.11)

Table 3 :

 3 Mean test errors (%) from real data examples, based on 30 repetitions. Standard errors are in the parentheses.

	Data	SFLDA FLDA RFLDA	PLS	B-Gauss B-NPD B-NPR Logistic QDA
	Tecator	5.0 (0.44)	3.4 (0.41)	7.4 (0.68)	4.8 (0.46)	6.3 (0.48)	4.2 (0.42)	4.8 (0.59)	1.8 (0.27)	3.3 (0.42)
	Wine	10.3 (0.81)	10.0 (0.74)	9.2 (0.72)	9.2 (0.81)	11.1 (0.91)	10.4 (0.89)	9.9 (0.57)	8.7 (0.65)	11.6 (0.80)
	Growth	6.0 (0.39)	5.7 (0.39)	8.3 (0.43)	6.1 (0.32)	6.2 (0.44)	7.2 (0.36)	7.1 (0.40)	5.9 (0.48)	6.7 (0.40)
	dGrowth	7.0 (0.57)	6.4 (0.47)	7.4 (0.52)	6.3 (0.52)	7.4 (0.48)	8.1 (0.49)	7.6 (0.55)	7.1 (0.56)	9.5 (0.60)
	Wheat	0.1 (0.09)	0.0 (0.00)	0.1 (0.33)	0.0 (0.00)	0.1 (0.09)	0.0 (0.00)	0.2 (0.13)	0.2 (0.17)	0.2 (0.13)
	Phoneme	19.7 (0.53)	21.6 (0.59)	20.9 (0.53)	21.4 (0.58)	23.7 (0.60)	23.8 (0.57)	24.2 (0.60)	22.1 (0.62)	24.4 (0.65)

Acknowledgment

We are grateful to Aurore Delaigle and Xiongtao Dai for sharing the Matlab codes for compared methods. We are also grateful to two referees, an associate editor and the editor for careful reading and helpful suggestions. The third author's work was supported by Basic Science Research Program of the National Research Foundation of Korea (NRF-2019R1A2C1005979) funded by the Korean government.

Supplementary material

Supplementary material includes additional details on relevant background material and proofs, as well as an example of matlab codes for implementation.