Elodie Duraffourg
email: elodie.duraffourg@onera.fr

Laurent Burlion
email: laurent.burlion@onera.fr

Tarek Ahmed-Ali
email: tarek.ahmed-ali@greyc.ensicaen.fr

Françoise Lamnabhi-Lagarrigue

Adaptive Finite-time Observer for a Nonlinear and Flexible Space Launch Vehicle

Keywords: Hybrid observer, Nonlinear control, Parameter estimation, Flexible mode NOMENCLATURE

In this paper we propose a new finite-time state observer for a nonlinear space launch vehicle with flexible dynamics and uncertain parameters. Indeed, flexible states are required to ensure nonlinear control objectives of both reference path tracking and bending mode damping. Our main contribution is to show that it is enough to observe a sublinear uncertain system to ensure a finite-time convergence of the whole state. For that purpose, a Luenberger observer is mixed with a parameter and initial state estimator, based on algebraic estimation tools. Closed loop simulations show the effectiveness of the observer in combination with a backstepping control design (extended to flexible systems).

The problem of finite-time observation for linear and nonlinear systems has been widely investigated over the past decades. Two classes of observers emerge in the series of methods that achieve finite-time convergence. The first one, based on the use of delays has deserved a lot of attention [Menold (2003)], [Engel (2002)]. Recently in [Karafyllis (2011)], a novel hybrid dead-beat observer which uses delays has been proposed. The history of the output is used in order to estimate the state of the system. Sliding mode observers that contain large study in the literature make the second class. (see [Shtessel (2010b)], [Ahmed-Ali (1999)] for instance). More recently, homogeneous finite-time observers have been developed for a specific class of nonlinear systems [Perruquetti (2008)]. Most of these approaches make the assumption that the system structure and parameters are known.

In this paper we propose to design a finite-time observer for a space launch vehicle which belongs to the class of uncertain nonlinear systems. Due to mass constraints, space vehicles tend to have lightweight and flexible structures with low natural frequencies, distorting sensors measurement and adding stability problems during flight. Researchers have recently investigated this subject in the field of nonlinear control. Several solutions have been proposed. Some of them, on the one hand, addressed the problem of unknown parameters and uncertainties using direct-adaptive [Fiorentini (2009)], or time varying controller [Hervas (2012)]. Nevertheless, only rigid states, (that means measured states), are used in these proposed methods. On the other hand, using sliding mode state observer, the authors of [Shtessel (2010a)] reconstruct the flexible states in order to remove the undesirable dynamics from the measurements. This approach requires unfortunately a strong knowledge of the mathematical model of the system, in particular flexible modes parameters. However, to the best of our knowledge, the design of a finite-time observer has not been achieved on uncertain nonlinear aerospace models.

As far as we are concerned, we recently designed a Lyapunov-based nonlinear controller, which uses the flexible states, to ensure control objectives of both reference path tracking and bending mode damping for a class of nonlinear and flexible system [Duraffourg (2013a)], [Duraffourg (2013b)], [Burlion (2013)]. Assuming that the whole state is available, this control law has been applied to the rotational dynamics of a space launch vehicle. [Duraffourg (2013c)]. Such assumptions do not hold in practical applications since flexible states are generally not measured. Consequently we need to reconstruct the flexible states. Besides, noting that flexible parameters are subject to uncertainties or variation during flight, this paper proposes to extend existing theory by proposing an indirect adaptive hybrid observer that no longer requires system parameter knowledge. The proposed approach consists in estimating flexible parameters (natural frequency and damping), and state initial conditions using algebraic tools [Fliess (2003)]. The first ones improve the accuracy and the robustness of the observer through indirect adaptive feature. The second ones are used to regularly update the estimated state and so guarantee a finite-time convergence. This paper is organised as follows. Section 2 describes the space launch vehicle and gives the problem statement. Section 3 develops a parameter estimator and a state observer which are then mixed to design a hybrid adaptive finite-time observer. In section 4, estimated state is used in a nonlinear control law and a closed-loop simulation is presented. Finally Section 5 contains our conclusions and future research directions. The rotational dynamics of a nonlinear flexible launch vehicle is extracted from [Duraffourg (2013c)], where Lagrange's formalism was developed to get the fullmathematical model. The equations of motion are given by:

PROBLEM STATEMENT

Launcher mathematical model

     ψ = q q = - l aero I L L(ψ) + T I L (L T r -h)η + T L T I L β η = -(ω 2 -hrT )η -2ξω η + hT β (1)
where ψ and η are real variables and the lift L is a nonlinear function of the attitude, given by:

L(ψ) = qS C 1 L ψ -C 2 L ψ 2
The launch vehicle is equipped with an inertial unit and a rategyro that give attitude and pitch rate informations. As the first bending mode distorts the sensors measurement, the available outputs are:

y 1 = ψ + r ci η y 2 = q + r gy η (2)
Measurement corruption terms make the control law design more difficult.

Nonlinear control law

A nonlinear control law, denoted Flexible Backstepping, that achieves the control objectives of both reference path tracking and bending mode damping has been developed in [Duraffourg (2013c)]. This controller limits the interaction of the rigid-dynamics on the transient of the flexible dynamics, and so, improves the damping of the bending mode. However, it is based on the deep knowledge of the flexible states and parameters.

Using the notation y 3 = y 2 + (r cir gy ) η, system (1) can be reformulated as follow:

     ẏ1 = y 3 ẏ3 = ḡ(y 1 , N )y 1 + KN + hT C β β + O(η 2 ) η = -ω 2 η -2ξω η + hT β (3)
where:

ω2 = ω 2 -hrT C β = I L h L T + I L r ci h ḡ(y 1 , N ) = - l aero qS I L C 1 L -C 2 L y 1 + 2C 2 L r ci η K = l aero qSC 1 L r ci I L + T I L (L T r -h) -r ci ω2 -2ξωr ci N = (η η) T
Then, the following change of coordinates is applied to the flexible-dynamics:

z = N -C β y 1 y 3 (4) System (3) becomes:      ẏ1 = y 3 ẏ3 = ḡ(y 1 , N )y 1 + KN + hT C β β + O(η 2 ) ż = A z z + F (y 1 , N )y 1 + Gy 3 + O(η 2 ) (5)
with:

A z = 0 1 -ω 2 -2ξω - 0 C β K G = A z 0 C β - C β 0 F (y 1 , N ) = A z C β 0 - 0 C β ḡ(y 1 , N )
A z being Hurwitz, system (5) belongs to the nonlinear class of system on which a flexible backstepping controller can be developed (see [Duraffourg (2013a)], [Duraffourg (2013b)]). The design of this flexible backstepping control law requires the whole state. Hence, outputs (y 1 and y 2 ) and the flexible states (η and η) are required.

Remark 1. The exact expression of the y 3 -dynamics is:

ẏ3 = ḡ(y 1 , η)y 1 + KN + hT C β β + l aero qSC 2 L r 2 ci I L η 2
Inertial flexible rotation r ci takes very low values (in the order of 10 -4 ). Thus the last term can be seen as an error term and is represented by a O(η 2 ).

PARAMETERS AND STATE ESTIMATION

Since the bending mode is not measured and the flexible parameters are generally distorted, this section proposes a way to estimate the flexible states and parameters and so to make the use of this control law possible via output feedback.

State observer design with known parameter (ideal case)

The flexible states η and η must be estimated. They are described by linear differential equations that result from system (1). Since the outputs (2) involve both the rigid states and the flexible ones, the idea is to consider the augmented state X = [ψ η η] T instead of just the required flexible states. It is important to note that this system is linear, contrary to the original nonlinear one (1).

Working with this augmented state gives the possibility to design a linear observer, by focusing on:

Ẋ = AX + B y y 2 + B β β y 1 = CX (6)
where y 2 acts as an input and:

A =   0 0 -r gy 0 0 1 0 -ω 2 -2ξω   B y = 1 0 0 B β = 0 0 hT C = (1 r ci 0)
Observability conditions holds for the pair (A, C) and the estimated state X = [ ψ η η] T is given by the classical Luenberger observer:

Ẋ = A X + B y y 2 + B β β + L(y 1 -C X) (7) where L ∈ R 3×1 is chosen such that A -LC is stable.

Parameters estimation

Bending mode natural damping and pulsation are generally subject to uncertainties and variations during the flight. Besides, since our flight control law intends to attenuate the oscillations of the bending mode, that is, to add damping, it is important to know these parameters accurately.

Moreover the accuracy and convergence time of the state observer can be improved by a better knowledge of the state initial conditions.

In this way we choose to estimate the following parameters:

θ =    θ 1 θ 2 θ 3 θ 4    =    -2ξω -ω 2 y 3 (0) ẏ3 (0) -ẏ2 (0) + 2ξωy 3 (0)    (8)
The first two parameters (θ 1 and θ 2 ) give the natural damping and pulsation of the bending mode. The flexible state initial conditions depend on the four parameters θ.

η(0) = 1 θ 2 θ 4 + θ 1 y 2 (0) r ci -r gy -hT β(0) η(0) = θ 3 -y 2 (0) r ci -r gy
r ci being small, we use the first initial output to approximate the state initial condition ψ(0):

ψ(0) = y 1 (0) -r ci η(0) ≃ y 1 (0) (9)
The flexible state initial conditions must be estimated through the outputs and their time derivatives. Noting that η = ẏ1y 2 r cir gy , it comes the following equation that links the outputs and their time derivatives to the input derivative, and consists in the basic equation of the parameter estimation: y 10) An algebraic methodology for parameter identification is described in [Fliess (2003)] . This approach was used here to estimate θ. The main steps are presented here:

(3) 1 = θ 2 ( ẏ1 -y 2 ) + θ 1 (ÿ 1 -ẏ2 ) + hT (r ci -r gy ) β + ÿ2 (
(1) Take the Laplace transformation of (10) to reveal the four parameters to be estimated. It comes:

θ 1 -s 2 Y 1 +s Y 2 +y 1 (0)s -y 2 (0) +θ 2 -sY 1 +Y 2 +y 1 (0) -sθ 3 -θ 4 = hT (r ci -r gy ) sB -β(0) -s 3 Y 1 + s 2 Y 2 + y 1 (0) -sy 2 (0) (11)
where s represents the Laplace variable and Y i (resp. B) is the Laplace transformation of signal y i (i ∈ {1, 2}) (resp. β).

(2) Take derivatives with respect to s, (three times) to get as equations as unknown parameters.

(3) Multiply by s -3 both sides to avoid time derivations. (4) Come back to time domain using inverse transformations. (5) Define time dependent matrices P ∈ R 4×4 and Q ∈ R 4×1 such that:

P (t)θ = Q(t) (12) Assumption 2. Persistence of excitation condition: ∃ τ > 0 such that: ∀ u ≥ τ, τ u-τ P (s) T P (s)ds > 0 (13)
Under Assumption 2, θ is finally given by:

θ = τ u-τ P (s) T Q(s)ds τ u-τ P (s) T P (s)ds ( 14
)
Since τ is unknown, we compute estimated parameters θ and then  and X0 at discrete time instants τ k where {τ k } +∞ k=0 is a partition of R + , using:

θ(τ k ) = τ k 0 P (s) T Q(s)ds τ k 0 P (s) T P (s)ds (15) Â(τ k ) = Âk =   0 0 -r gy 0 0 1 0 θ2 (τ k ) θ1 (τ k )   (16) X0 (τ k ) =       y 1 (0) 1 θ2 (τ k ) θ4 (τ k ) + θ1 (τ k )y 2 (0) r ci -r gy -hT β(0) θ3 (τ k ) -y 2 (0) r ci -r gy       (17)
Remark 3. It is important to note that  and X0 are not completely described by estimated parameters. Parameters r gy , r ci , h and T are supposed to be known. This may be assumed for flexible displacements and rotations since they take low values, and cannot vary significantly during the flight.

Figure 2 shows in blue (resp. in red) the evolution of the flexible estimated parameters ( θ1 , θ2 ) and state initial conditions X0 = ψ(0), η(0), η(0) when the sensors give ideal (resp. noisy) measurements. The red curve was obtained applying a zero-mean periodic noise on signals y 1 and y 2 . Remark 4. With noisy measurements, estimated parameter η(0) is biased and the estimation error η = η(0) -η( 0) is quite important (about 1.5). This problem is detailed in section 3.3.

Robustness improvement of the observer with respect to noisy measurements

An estimation bias appears in simulation on estimated parameter η(0) when measurements are noisy (solid red curve on the last plot of figure 1). It is given by:

η(0) = η(0) -η(0) = θ 3 -θ3 r ci -r gy (18)
This bias term is due to the low values of r cir gy which accentuates the (low) difference between θ 3 and θ3 . This section proposes a way to identify this estimation bias η(0) and thus to correct X0 (τ k ).

From system (6), it comes:

y 1 (t) = CX(t) = Ce At X 0 (t) (19) + C t 0 e A(t-u) B y y 2 (u) + B β β(u) du
Similarly with the estimated state:

X(t) = e Âk t X0 (t) + t 0 e Âk (t-u) B y y 2 (u) + B β β(u) du (20)
Under persistence of excitation condition, algebraic parameter estimation still converges in finite-time. Thus, there exists k * > 0 such that for all k ≥ k * Âk = A and X0 = ψ(0) η( 0) η(0) -η(0) T (21) This last equation underlines the fact that only estimated initial condition η is biased.

Then, ∀t > t k * y 1 (t) -C X(t) = Ce At e 3 η(0) (22)
with e 3 = (0 0 1) T . That can be written as:

t 0 (Ce Au e 3 ) 2 du η(0) = t 0
Ce Au e 3 y 1 (u) -C X(u) du (23) Supposing that t 0 (Ce Au e 3 ) 2 du = 0, estimation error is given by:

η(0) = t 0 Ce Au e 3 y 1 (u) -C X(u) du t 0 (Ce Au e 3 ) 2 du (24)
This term is added on the last simulation. The result is represented by the dashed red line on figure 1.

Adaptive Finite-time observer

Mixing the results of the last subsections, it is now possible to design an observer when the parameters are unknown. Estimated parameters and initial conditions are used to improve the accuracy of the observer. In particular, θ1 and θ2 are used in the design such that the observer no longer depends on the natural damping and pulsation of the bending mode, that are subject to variations. Proposition 5. Under Assumption 2, the following hybrid observer converges in finite-time.

             ∀t ∈ [τ k , τ k+1 [ Ẋ = Âk X + B y y 2 + B β β + L(y 1 -C X) when t = τ k X(τ k ) = e Âk τ k X0 (τ k ) + τ k 0 e Âk (τ k -u) B y y 2 (u) + B β β(u) du (25)
At each τ k , the state is updated and the dynamics ( Âk ) as well.

= A -Â, it comes: Ẋ = AX -Â X -LC X = (A -LC) X -Ã X (26)
Because of the persistence of excitation condition, algebraic parameter estimation converges in finite-time. Thus, there exists k * such that τ k * > τ > 0 and

∀ k > k * , θ(τ k ) = θ (27)
Consequently ∀ k > k * , Âk = A and the estimation error satisfies:

∀ t ≥ τ * k , Ẋ = (A -LC) X (28)
A -LC being Hurwitz, the estimation error X converges asymptotically to zero.

Besides, the estimated state X(t) is updated on X(τ k ) at each

τ k verifying τ k > τ k * . ∀t ≥ τ * k , X0 (t) = X 0 (t) ie X(τ k ) = X(τ k ) ( 29
)
Estimation error is then given by: ∀t ≥ τ * k , X(t) = X(τ k )e (A-LC)(t-τ k ) = 0 (30) Finally the estimation error vanishes in finite-time.

Figure 3 compares this finite-time observer that updates the state from the knowledge of the initial condition (IC) (in blue) with the same Luenberger observer without any update of the estimated state (in green). In the two cases, the observer parameter L is the same. Moreover, initial conditions of the classical Luenberger observer have been chosen very close to the estimated parameters, so that the comparison is fair. This figure shows how the observer convergence time is improved. The dramatic slope on the blue curve corresponds to the moment where the estimated state is updated.

CLOSED-LOOP SIMULATION

Finite-time observer was designed in open-loop. In this section, the loop is closed and estimated flexible states and parameters are directly used in the controller.

As explained in Section 2.2, a flexible backstepping control law can be designed for system (5). It is extracted from [Duraffourg (2013c)], [Duraffourg (2013b)] and involved the following Lyapunov function:

V = c f 1 2 y 2 1 + c f 2 2 Z T P z Z + c f 3 2 ỹ2 3 (31)
where c f i are positive constants (i ∈ {1, 2, 3}), P z ∈ R 2×2 is the positive and symmetric matrix that verifies:

A T z P z + P z A z = -2Q z (32) with Q z a positive and symmetric matrix of R 2×2 , and

A z = 0 1 θ 2 θ 1 - 0 C β K K = l aero qSC 1 L r ci I L + T I L (L T r -h) + r ci θ 2 r ci θ 1 Z = z -Gy 1 , G = A z 0 C β - C β 0 F (y 1 , N ) = A z C β 0 - 0 C β ḡ(y 1 , N ) ỹ3 = y 3 -y 3 cmd y 3 cmd = - 1 c f 1 λ y1 y 1 + c f 2 F (y 1 , N ) T P z Z F (y 1 , N ) = A z G + F (y 1 , N
) Choosing the control law as follow:

β f (y 1 , y 3 , N ) = C β hT ẏ3 cmd -ḡ(y 1 , N )y 1 -KN - 1 c f 3 (c f 1 y 1 + λ y3 ỹ3 ) λ y1 > 0 λ y3 > 0 (33)
Lyapunov function time derivative is given by:

V = -λ y1 y 2 1 -c f 2 Z T Q z Z -λ y3 ỹ2
3 (34) This nominal control law is then blended with the adaptive observer, using estimated flexible states N = η η . The tested control law is thus β f (y 1 , y 3 , N ), reminding that flexible states converge in finite time.

This output-feedback flexible backstepping control law is compared with a classical backstepping one, applied on the sole rigid dynamics:

β c = I l T L T qcmd + L aero I L L(ψ) - 1 c c 2 (c c 1 ψ -λ q q) (35)
where q = qq cmd q cmd = -λ ψ ψ (36) and c c 1 , c c 2 , λ ψ , λ q are positive constants. The classical backstepping law also consists of an output-feedback controller, assuming that state ψ (resp. q) approximately corresponds to output y 1 (resp. y 2 ). This approximation makes sense thanks to the low values of the sensors flexible rotations r ci and r gy in equation (2). Simulation is realised in this way.

It is worth noting that both flexible (33) and classical (35) backstepping laws converge asymptotically when the flexible mode is collocated.

Figure 4 presents a simulation result of the state temporal evolution with the two different backstepping laws. Gains have been tuned so that the time delay response is the same. Initial conditions are also the same. Regardless of the control law, attitude ψ has nearly the same behaviour. As expected, the difference appears on the flexible states: the bending mode oscillations are better damped by the flexible backstepping controller. However it is worth noting that peaks appear at the beginning of the simulation (on the blue curve and particularly on q and η). This phenomenon is caused by state initial condition estimation and estimated state update. Indeed, at time equals to 0.2 s state initial conditions are obtained and estimated state is thus updated at that time. Once this process has been done, states get a smooth behaviour. The complete proof of the stability of the closed-loop will be addressed in the future.

Finally, although we have not given a formal proof, simulations have been performed with noisy measurements. More precisely, a zero-mean periodic noise has been applied on the outputs. In this case, as illustrated in red on figure 4 for the closed loop system, our finite-time convergent estimation algorithm well performs.

CONCLUSION

An adaptive finite-time observer has been designed in this paper. Based on a linear state observer and a parameter estimator, it identifies the unmeasured flexible states, without requiring parameter knowledge. Flexible states are then used to realise an output-feedback controller via flexible backstepping control law. Methodology and simulation results have been developed on the rotational dynamics of a space launch vehicle with only one flexible mode. Although this first flexible mode is the most relevant to consider in the design, more flexible modes will be considered in the future. In this study, flexible parameters (natural pulsation and damping) are estimated and then used in the control law. Uncertainties also affect other parameters such as aerodynamic coefficients. Theses uncertainties will be considered in our future work. Finally we will focus on the stability proof of the closed-loop system.

  r ci Inertial unit flexible rotation, rad r gy Rategyro flexible rotation, rad T Thrust, kg.m/s 2 L Lift, kg.m/s 2 G L Launcher center of mass C T Gimbal joint F L Aerodynamic center I L Launcher body inertia, kg.m 2 L T Algebraic distance from G L to C T , m l aero Algebraic distance from F L to G L , m q Dynamic pressure, Pa S Reference area of the vehicle, m 2 ω Natural frequency of the first bending mode, rad/s ξ Natural damping of the first bending mode R n Frame (G L , x n , y n ) linked to the reference trajectory R L Frame (G L , x L , y L ) linked to the launcher 1. INTRODUCTION

Fig. 1 .

 1 Fig. 1. Schematic model for a flexible launcher

Fig

  Fig. 2. Estimated parameters

  Fig. 3. Observation errorsProof. Using the notations X = X -X and à = A -Â, it comes:Ẋ = AX - X -LC X = (A -LC) X -à X(26)

Fig. 4 .

 4 Fig. 4. Evolution of attitude, pitch rate and flexible states

P.H. Menold, R. Findeisen and F. Allgöwer (2003), Finite time convergent observers for nonlinear systems, in Proc. of the 42nd IEEE Conference on Decision and Control, Hawaï, USA R. Engel and G. Kreisselmeier (2002), A Continuous-Time Observer Which Converges in Finite Time, IEEE Transactions on Automatic Control, Vol.47, No. 7, pp. 1202-1204 I. Karafyllis and Z.P. Jiang (2011), Hybrid dead-beat observers for a class of nonlinear systems, Systems and Control Letters, vol.60, pp. 608-617 Y.B. Shtessel, S. Baev, C. Edwards and S. Spurgeon (2010b), HOSM Observer for a Class of Non-Minimum Phase Causal Nonlinear MIMO Systems, IEEE Transactions on Automatic Control, Vol.55, No.2, pp. 543-548 T. Ahmed-Ali and F. Lamnabhi Lagarrigue (1999), Sliding Observer-Controller Design for Uncertain Triangular Nonlinear Systems, IEEE Transactions on Automatic Control, Vol. 44, No.6, pp. 1244-1249 W. Perruquetti, T. Floquet and E. Moulay (2008), Finite-Time Observers: Application to Secure Communication, IEEE Transactions on Automatic Control, Vol. 53, No. 1, pp. 356-360 L. Fiorentini, A. Serrani, M.A. Bolender, D.B. Doman (2009), Nonlinear Robust Adaptive Control of Flexible Air-Breathing Hypersonic Vehicles, Journal of Guidance, Control and Dynamics, Vol. 32, No. 2, pp. 401-416 J.R. Hervas and M. Reyhanoglu (2012), Control of a Spacecraft with Time-Varying Propellant Slosh Parameters, in Proc. of the 12th International Conference on Control Automation and Systems, Jeju Island, Korea Y.B. Shtessel, C.E. Hall, S. Baev and J.S. Orr (2010a), Flexible Modes Control Using Sliding Mode Observers: Application to Ares I, In Proc. of AIAA Guidance, Navigation and Control Conference, Canada E. Duraffourg, L. Burlion, T. Ahmed-Ali and F. Lamnabhi Lagarrigue (2013a), Non-linear Control of the Longitudinal Rotational Dynamics of a Flexible Aircraft, in Proc. of the 12th European Control Conference, Zürich, Switzerland, 2013 E. Duraffourg, L. Burlion, T. Ahmed-Ali and F. Lamnabhi Lagarrigue (2013b), Non-linear full-state control of a flexible hypersonic vehicle, in Proc. of 11th IEEE International Workshop on Electronics Control Measurement Signals and their application to Mechatronics, Toulouse, France E. Duraffourg, L. Burlion and T. Ahmed-Ali (2013c), Longitudinal Modeling and Preliminary Control of a Nonlinear Flexible Launch Vehicle, in Proc. of 11th IFAC International Workshop on Adaptation and Learning in Control and Signal Processing, Caen, France L. Burlion, E. Duraffourg, T. Ahmed-Ali and F. Lamnabhi Lagarrigue (2013), Global asymptotic stabilization for some nonlinear models of flexible aerospace vehicles, in Proc. of 52nd IEEE Conference on Decision and Control, pp. 4230-4235, Florence, Italy M. Fliess and H. Sira-Ramirez (2003), An Algebraic Framework for Linear Identification, ESAIM Control, Optimisation and Calculus of Variations, Vol.9, pp. 151-168


Yves Verhertbruggen

Camille Alvarado

Axelle Bouder

Xavier Falourd

Loïc Foucat

Audrey Geairon

Fabienne Guillon

Sophie Le Gall

Jacqueline Vigouroux

Mark D Wilkinson

Fabian Stritt

Markus Pauly

Mi Yeon Lee

Jenny C Mortimer

Henrik V Scheller

Rowan A C Mitchell

Catalin Voiniciuc

Luc Saulnier

Anne-Laure Chateigner-Boutin

ELUCIDATION OF THE STRUCTURE AND BIOSYNTHESIS OF WHEAT GRAIN MANNAN

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

making, pasta, brewing...). In wheat grain, the two major cell wall polysaccharides, arabinoxylans and β-glucans, have been extensively studied. By contrast, mannan was poorly documented.

Nevertheless, this hemicellulosic polysaccharide might have a determinant role in wheat grain development since, in Arabidopsis thaliana, mutants with a reduced amount of mannan show an altered seed development. Our project aimed at deciphering the biochemical structure and the biological function of mannan.

To study the biological function of mannan, we have studied its occurrence during wheat grain development. The detection of mannan using monoclonal antibodies showed that mannan is mainly located in the grain endosperm and that its deposition starts at an intermediate step of endosperm cellularization. We have purified mannan from wheat endosperm to determine its fine structure.

Former studies have speculated that the class of mannan present in wheat was glucomannan. Our data suggest that this class of hemicelluloses is only represented by short unsubstituted chains of 1,4 linked D-mannose residues that are weakly acetylated [1].

We identified candidate mannan biosynthetic genes expressed in the endosperm. After unsuccessful attempts to generate wheat lines depleted in mannan, the heterologous expression of our best candidate genes was carried out in Pichia pastoris and Arabidopsis triple mutant depleted in glucomannan. Introducing TaCslA12 gene in yeast and Arabidopsis led to the production of wheat-like linear mannan lacking glucose residues and with moderate acetylation [2]. Therefore, this gene encodes a mannan synthase and is responsible for the synthesis or wheat endosperm mannan.


Jérôme J Mariette

Madalina Olteanu
email: madalina.olteanu@univ-paris1.fr

Nathalie Vialaneix

Kernel and dissimilarity methods for exploratory analysis in a social context

come

Introduction

While most of statistical methods for prediction or data mining have been built for data made of independent observations of a common set of p numerical variables, many real-world applications do not fit in this framework. Typical such examples include categorical variables, relations between entities (e.g., a graph or network) or even more complex frameworks such as categorical time series. A particularly useful simplification of these more general situations is the case where a relevant similarity or dissimilarity can be computed between the observations, providing a summary of their relations to each other. In addition, when this similarity has some mild additional properties, it is called a kernel and provides a strong mathematical framework [START_REF] Berlinet | Reproducing Kernel Hilbert Spaces in Probability and Statistics[END_REF] for extending most of standard statistical supervised and unsupervised methods to any type of data for which a relevant such kernel can be obtained [START_REF] Cristianini | An Introduction to support vector machines[END_REF][START_REF] Shawe-Taylor | Kernel methods for pattern analysis[END_REF]. This approach has already proven useful in computational biology [START_REF] Schölkopf | Kernel Methods in Computational Biology[END_REF] or in social sciences and humanities [START_REF] Boulet | Batch kernel SOM and related Laplacian methods for social network analysis[END_REF][START_REF] Massoni | Which distance use when extracting typologies in sequence analysis? An application to school to work transitions[END_REF].

Nevertheless, the choice of a relevant kernel is still an open problem. Some authors have proposed to combine all candidate kernels into a "meta kernel" which is an "optimal" linear or convex combination of the individual kernels. This approach is known as the "multiple kernel learning problem" and has been widely studied in the supervised framework [START_REF] Gönen | Multiple kernel learning algorithms[END_REF]. The present chapter aims at presenting the less addressed unsupervised framework. More precisely, after a brief introduction to kernels and their relation with the more general similarity/dissimilarity settings (Section 2), we describe how statistical methods can be extended to the kernel framework by using the so-called "kernel trick" (Section 3). Section 4 focuses more precisely on the extension of an exploratory method, called Self-Organizing Maps [START_REF] Kohonen | Self-Organizing Maps, 3rd Edition[END_REF], to the kernel framework and discusses the issue of complexity and how it can be solved in this particular setting. Section 5 explains how kernels can be combined in an unsupervised setting, as a processing prior the unsupervised methods presented before. The overall approach is illustrated in Section 6 on categorical time series in a social science context: originally developed in bioinformatics, sequence analysis is indeed increasingly used in social sciences for the study of life-course processes. In this section, we discuss how the choice of a given type of dissimilarity or group of dissimilarities influences the output of the exploratory analysis and allows to extract relevant patterns from this particular kind of data.

2 Kernels and more general proximity data

Kernels and RKHS

Kernel methods consider the case where data are described by a kernel obtained from a Reproducing Kernel Hilbert Space (RKHS; [START_REF] Berlinet | Reproducing Kernel Hilbert Spaces in Probability and Statistics[END_REF]). Usually, the sample of interest takes values in an arbitrary space, X that encompasses a variety of data types. This sample is then described by a function K : X × X → R, which is symmetric (∀ x, x ∈ X, K(x, x ) = K(x , x)) and positive (∀ N ∈ N, ∀(α i ) i=1,..., N ⊂ R and ∀ (x i ) i=1,..., N ⊂ X, N i,i =1 α i α i K(x i , x i ) ≥ 0) and is called the kernel. Indeed, in this case, it is known [START_REF] Aronszajn | Theory of reproducing kernels[END_REF][START_REF] Berlinet | Reproducing Kernel Hilbert Spaces in Probability and Statistics[END_REF] that there exists a unique Hilbert space (H, ., . H ) and a unique application φ : X → H , such that ∀ x, x ∈ X, φ(x), φ(x ) H = K(x, x ).

(H, ., . H ) is the RKHS of K and is also often called feature space; φ is the feature map of K.

In statistics and machine learning, this framework is often used to deal with observations that are not just multidimensional vectors (e.g., categorical time series or graphs, among others) or to incorporate expert knowledge in the analysis (see Examples 1 and 2 below with standard examples of kernels often used in practice). The sample (x i ) i=1,...,n is then described by pairwise relations between observations, as measured by the kernel. This leads to the computation of the kernel matrix K = (k ii ) i,i =1,...,n , with k ii = K(x i , x i ), which is symmetric and semi-definite positive, by definition of the kernel K.

The idea of kernel methods is to perform standard linear statistical analyses in the feature space (H, ., . H ). Since the only operations involved in these analyses are related to the computation of dot products and norms, the Hilbert space H and the feature map φ are usually not explicitly given but used implicitly through the kernel K instead. This principle, which we illustrate below, is called the kernel trick.

Example 1 Some useful kernels Kernels in R p . Kernel methods are often used for standard multidimensional data to provide more flexibility and non linearity in the analyses. In these spaces, a trivial kernel is given by using the standard dot product of R p : K(x, x ) = (x ) x, which leads to the trivial feature map φ = Id. The feature space is then unchanged as compared to the original space (X = H = R p ) and the performed statistical analysis is thus still linear. Among more interesting kernels for R p , one of the most popular is the Gaussian kernel (also called Radial Basis Function -RBF -kernel) K γ (x, x ) = e -γ x-x 2 , which shape is controlled by a hyper-parameter γ > 0. This kernel is of special importance since it is continuous and universal for every compact set of X, C (meaning that the set of all functions induced by x ∈ C → K(x, .) is dense in the set of all continuous functions C ⊂ X → R). This property allowed Steinwart [START_REF] Steinwart | On the influence of the kernel on the consistency of support vector machines[END_REF][START_REF] Steinwart | Support vector machines are universally consistent[END_REF] to demonstrate the consistency of kernel classification and regression methods in the statistical sense (when the sample size grows to infinity and in terms of convergence of the error loss to its optimum). The polynomial kernel (K(x, x ) = (1 -(x ) x) γ for γ > 0) and the exponential kernel (K(x, x ) = e (x ) x ) are also universal kernels.

Kernels on graphs. In many application fields including social sciences and biology, graphs (also called networks) are widely used to represent pairwise relations between entities (friendship, professional contacts, regulation between genes, ...). A number of kernels for graphs have been proposed to provide a similarity measure between nodes based on the graph structure. Most of them are derived from regularized version of the Laplacian of the graph [START_REF] Kondor | Diffusion kernels on graphs and other discrete structures[END_REF][START_REF] Smola | Kernels and regularization on graphs[END_REF] and have been used in prediction or exploratory analyses in biology (e.g., for introducing known relations between genes [START_REF] Vert | Extracting active pathways from gene expression data[END_REF]39]) or in social sciences (e.g., to extract information from a medieval social network [START_REF] Boulet | Batch kernel SOM and related Laplacian methods for social network analysis[END_REF]).

From general similarities to kernels

In practice, data are often described by similarities (or dissimilarities) that are not necessarily definite positive (see Example 2). This situation is addressed either by generalizing kernel methods to the "pseudo-Euclidean" framework [START_REF] Goldfarb | A unified approach to pattern recognition[END_REF][START_REF] Ong | Learning with non-positive kernels[END_REF], by embedding the sample directly into a Euclidean space whose dot product ressembles the original similarity (Multidimensional Scaling -MDS-is one of these approaches [START_REF] Cox | Multidimensional Scaling. Chapman and Hall/CRC[END_REF]), or by using a proper definite kernel instead of the original indefinite similarity.

In the latter case, the chosen kernel is often obtained by a simple transformation of its spectrum meant to obtain only positive eigenvalues. [START_REF] Chen | Similarity-based classification: concepts and algorithm[END_REF][START_REF] Schleif | Indefinite proximity learning: a review[END_REF] are two reviews describing the topic of general similarity learning and its relation with kernel methods.

Example 2 Some more general similarities and dissimilarities Categorical sequences or time series. Categorical sequences are naturally used in biology to represent the DNA sequences or proteins (with the categories being the amino acids). Among the many proposals for quantifying the similarities between two sequences, edit distances (also known as "Levenshtein distances" or "optimal matching dissimilarities") [START_REF] Needleman | A general method applicable to the search for similarities in the amino acid sequence of two proteins[END_REF] are one of the most famous. Their main idea is to quantify the minimum number of transformations needed to obtain a sequence from another one. A cost is associated to insertion, deletion and substitution transformations to allow a flexible customization of these dissimilarities. These measures have been increasingly used in social sciences as well, for studying life-course processes [START_REF] Abbott | Sequence analysis and optimal matching methods in sociology[END_REF][START_REF] Massoni | Which distance use when extracting typologies in sequence analysis? An application to school to work transitions[END_REF]. Section 6 describe in further details those dissimilarities.

Dissimilarities based on phylogeny. As already mentioned, kernels and dissimilarities can also embed prior expert information in their computation.

A typical example is the case where variables are the abundances of different species for which a phylogeny information (a parental information between those species) is given. Such frameworks are met when studying the biodiversity of different places or in metagenomics for instance. In these applications, data are described by vectors of counts that represent the number of times given species or Operational Taxonomic Units (OTUs) have been found for a given individual. For these data, computing a measure of proximity between observations that accounts for the distances between the species has been shown to provide a more relevant information than the simple Euclidean distance between counts [START_REF] Lozupone | UniFrac: a new phylogenetic method for comparing microbial communities[END_REF][START_REF] Lozupone | Quantitative and qualitative β eiversity measures lead to different insights into factors that structure microbial communities[END_REF]. Such distances include the (weighted) UniFrac distance or the generalized UniFrac distance [START_REF] Chen | Associating microbiome composition with environmental covariates using generalized UniFrac distances[END_REF].

3 Basics of statistical learning with kernels

Supervised setting

A simple example of a supervised learning method that has been extended to kernels is the ridge regression. More precisely, when given a training sample {(x i , y i ) i=1,...,n } for which x i ∈ R p and y i is a real number, the ridge regression finds the best linear predictor for (y i ) i based on (x i ) i that minimizes the squared loss plus a regularization term based on the 2 norm:

β * = argmin β ∈R p n i=1 y i -β x i 2 + λ β 2 (1)
for a given λ > 0, called regularization parameter, which is usually tuned by a cross validation approach. The solution of Equation ( 1) is given by:

β * = n i=1 x i x i + λI p -1 n i=1 y i x i ,
that can also be written as

β * = XX + λI n -1 X y,
with y = (y 1 , . . . , y n ) and X being the (n × p)-matrix with rows containing the x i so that the matrix XX is the (n × n)-matrix with entries the pairwise dot products x i x i for all i, i = 1, . . . , n. In summary, the solution writes

β * = n i=1 α * i x i with α * = XX + λI n -1 y. (2)
The extension of this approach to samples (x i ) i taking values in an arbitrary space X through the use of kernels is called kernel ridge regression [START_REF] Saunders | Ridge regression learning algorithm in dual variables[END_REF]. The idea is simply to search for a linear predictor in the feature space induced by the kernel, H , which transforms the optimization criterion of Equation (1) into:

w * = argmin w ∈H 1 n n i=1 (y i -w, φ(x i ) H ) 2 + λ w 2 H . (3)
The best linear predictor in H is thus given similarly as the solution of Equation ( 2) but in the feature space, replacing x i by φ(x i ) and the R p dot product by ., . H . In particular, this means that the matrix XX is replaced by a matrix with entries equal to φ(x i ), φ(x i ) H , which, by the so-called kernel tricks, turns out to simply be equal to K. We thus have that

w * = n i=1 α * i φ(x i ) with α * = (K + λI n ) -1 y.
This result can also be found as a consequence of the Representer Theorem [START_REF] Kimeldorf | A correspondence between Bayesian estimation on stochastic processes and smoothing by splines[END_REF][START_REF] Schölkopf | A generalized representer theorem[END_REF] or directly solving the dual of Equation ( 3):

α * = argmin α ∈R n 1 n n i=1 y i -α K i 2 + λ α 2 K , (4)
in which K i is the i-th row of the kernel matrix K and . K is the 2 norm induced by this matrix in R n : α 2 K = α Kα. Variants of this framework include Support Vector Machines (SVM, [START_REF] Boser | A training algorithm for optimal margin classifiers[END_REF]), for the classification case, or -SVM, for the regression case. In both cases, the main difference with the kernel ridge regression lies in the loss function but the main principle of the approach remains identical: the Representer Theorem allows to express the solution as a linear combination of the images by φ of the observations and the solution is obtained by solving a dual optimization problem obtained thanks to the use of the kernel trick.

Unsupervised setting

Kernel methods have also been developed for the unsupervised setting. Among the most direct of these extensions, the generalization of PCA [START_REF] Schölkopf | Nonlinear component analysis as a kernel eigenvalue problem[END_REF] and that of k-means [START_REF] Dhillon | Kernel k-means, spectral clustering and normalized cuts[END_REF] are probably the most known and used. They both use approaches similar to the supervised case described in the previous section, and more precisely:

• computations related to the original method (i.e., standard PCA and k-means) are performed in the feature space; • to do so, the kernel trick is used instead of the standard computation of dot products or norms.

Kernel PCA.

Standard PCA is often presented as the eigendecomposition of the variance/covariance matrix associated to the (n × p)-matrix of sample measures, X.

Assuming without loss of generality that X is centered, this eigendecomposition is equivalent to the dual eigendecomposition of XX , that provides the coordinates (or scores) of the projection of X on the different principal components. More precisely, if T is the (n × k) column matrix with the first k eigenvectors of XX , orthogonal and with a norm equal to 1 √ λ j

, then the (p × k) column matrix of the principal components (orthogonal and with a norm equal to 1) is X T.

Kernel PCA uses a similar approach taking advantage of the analogy between XX and K and between the i-th row of X and φ(x i ). More precisely, 1. assuming that K is centered in the feature space1, the eigendecomposition of K is obtained. It gives (λ j ) j=1,...,k , the first k eigenvalues of K, and (t j ) j=1,...,k , the associated first k orthogonal eigenvectors with a norm equal to 1 √ λ j ; 2. the first k (orthogonal) principal components are thus obtained as

w j = n i=1 t ji φ(x i ),
and have a norm equal to 1 in H . The coordinate of φ(x i ) on the j-th principal component is thus w j , φ(x i ) H = λ j t ji .

Kernel k-means.

Similarly, kernel k-means performs a standard k-means algorithm in the feature space H . To do so, in addition to computing dot products and norms using the kernel trick, it is necessary to obtain a representation of the cluster barycenters. More precisely, if (x i ) i ∈C are the observations assigned to a given cluster C, then, the barycenter is given by:

xC = 1 |C| i ∈C φ(x i )
and its (squared) distance to any other observation, x i , in the sample is obtained by:

1 if K is not centered, the centering operation is simply K -1 n 1 n K -1 n K1 n + 1 n 2 1 n K1 n , in which 1 n is an n × n matrix with all entries equal to 1. φ(x i ) -xC 2 H = φ(x i ) - 1 |C| i ∈C φ(x i ) 2 H = φ(x i ) 2 H - 2 |C| i ∈C φ(x i ), φ(x i ) H + 1 |C| 2 i i ∈C φ(x i ), φ(x i ) H = k ii - 2 |C| i ∈C k ii + 1 |C| 2 i ,i ∈C k i i .
In both situations (kernel PCA and kernel k-means), the adaptation of the algorithms to kernel data is made by their direct rewriting in the feature space. New data points, that were not previously in the feature space (principal components or barycenters) are represented by linear combinations of the images by the feature map, φ, of observations. In addition, distances to these new elements can be expressed in function of the kernel, using the kernel trick. These adaptations are thus very similar to the supervised case situations.

Kernel self-organizing maps and complexity reduction

In this section, we present an extension of kernel k-means to a more general method, which simultaneously performs clustering and dimensionality reduction for visualization, namely the self-organizing map (SOM) algorithm. Originally designed for unsupervised exploration of standard numerical datasets [START_REF] Kohonen | Self-Organizing Maps, 3rd Edition[END_REF], the method has been extended to handle non numeric data by using approaches based on Multiple Correspondence Analysis [START_REF] Cottrell | How to use the Kohonen algorithm to simultaneously analyse individuals in a survey[END_REF] or by relying on an algorithm that represents all the clusters by a prototype chosen among the data (median SOM, [START_REF] Kohonen | Self-organizing maps of symbol strings[END_REF]). Even if very general, the latter approach is very restrictive and generates representation issues, with associated biases in the obtained maps. In the present section, we present the extension of SOM to kernels that has been introduced by several authors for batch and online versions [START_REF] Donald | The kernel self organising map[END_REF][START_REF] Boulet | Batch kernel SOM and related Laplacian methods for social network analysis[END_REF] and has also been generalized to data represented by general dissimilarities (rather than kernels) [START_REF] Hammer | Topographic mapping of large dissimilarity data sets[END_REF]. The second part of this section will discuss associated complexity issues when the sample size is large, and review the different strategies that can be implemented to overcome them.

Kernel self-organizing maps

For the standard case of a dataset (x i ) i=1,...,n of multidimensional observations x i ∈ R p , SOM algorithm is close to k-means algorithm, except that the clusters are organized on a map equipped with a distance, d. More precisely, a map (also sometimes called a grid) is a set of U clusters (also sometimes called units or neurons) associated to physical locations in a low dimensional space. The clusters are frequently positioned in R 2 at coordinates (a, b) a=1,..., A, b=1,...,B with AB = U. Clus-ters are related to each other using pairwise distances, that can be, for instance, the Euclidean distances between their coordinates in R 2 . In addition, every cluster, u, is summarized by a prototype, p u that takes its values in the input space R p .

When fixing the number of neurons U, one should take into account that SOM is mainly intended as a method for nonlinear mapping and dimensionality reduction -in the sense of vector quantization-than as a method for clustering the data into a small number of clusters. Some authors [START_REF] Ultsch | Kohonen's self organizing feature maps for exploratory data analysis[END_REF] suggest to build very large SOMs, with a number of units U larger than the sample size, n. In this context, SOM essentially reduces to nonlinear mapping and to mining the underlying distribution of the data. A second option, which is more commonly used in practice, consists in building medium size maps that are smaller than the sample size, but still large enough to have a few input observations representing each unit [START_REF] Kohonen | Self-Organizing Maps, 3rd Edition[END_REF]. This strategy is a good trade-off between mapping and clustering, and a heuristic suggests to set U close to n/10 [START_REF] Villa-Vialaneix | Stochastic self-organizing map variants with the R package SOMbrero[END_REF].

The method aims at assigning every observation in the dataset to one of the clusters, while minimizing the distortion of the topology between the original space (here, R p ) and the map (as seen through the distance d). The prototypes are thus expected to be representative of the observations assigned to their cluster, as the barycenter is representative of its cluster in kernel k-means. To do so, the stochastic version of the method iterates over two steps:

• an assignment step in which an observation, x i is randomly chosen and assigned to the unit with the closest prototype:

f (x i ) = argmin u=1,...,U x i -p u 2
where f (x i ) is the cluster to which observation x i is assigned; • a representation step in which all prototypes are updated according to the new assignment:

∀ u = 1, . . . , U, p u ← p u + µH(d( f (x i ), u))(x i -p u )
where µ > 0 is chosen so as to vanish during the training process and H is a decreasing function, generally chosen such that H(0) = 1 and lim z→+∞ H(z) = 0.

The method is usually initialized with a random choice of the prototypes in R p . Different heuristics can be used to choose T: either it is defined proportionally to the sample size (typically, T = 5n) or it is not fixed in advance and the algorithm stops when the iterations no longer modify the solution.

The extension of SOM to kernel data is based on the same key tools as the ones described for kernel PCA and kernel k-means, and which allow to re-write the algorithm in the feature space H :

• the prototypes are expressed as convex combinations of the images by φ of the observations, and the assignment step is written in terms of coefficients related to each image φ(x i ); • the representation step is expressed with K by means of the kernel trick.

The full version of the method is provided in Algorithm 1.

Algorithm 1 Stochastic kernel SOM

1: ∀u = 1, . . . , U and ∀i = 1, . . . , n, random initialization of the prototypes:

p 1 u = n i=1 β 1 ui φ(x i ) with β 1 ui ∈ [0, 1] and i β 1 ui = 1 2: for t = 1 to T do 3:
Select randomly one observation i ∈ {1, . . . , n} Assignment step

f t +1 (x i ) = argmin u=1, ...,U φ(x i ) -p t u 2 H = argmin u=1, ...,U k ii -2 n l=1 β t ul k il + n l, l =1 β t ul β t ul k ll 4:
For all u = 1, . . . , U, Representation step

p t +1 u = p t u + µ t H t (d( f t +1 (x i ), u))(φ(x i ) -p t u ) ⇔ β t +1 u = β t u + µ t H t (d( f t +1 (x i ), u)) 1 n i -β t u ,
where 1 n i is a vector of length n with all entries equal to 0 except for the ith, which is equal to 1. 5: end for 6: return (p T +1 u ) u (prototypes) and

( f T +1 (x i )) i (clustering)

Complexity of kernel SOM

Kernel methods are generally considered efficient to deal with large dimensional data (when the original space is a standard multidimensional space, R p , with p large) but often encounter scalability issues when the sample size, n, becomes large. As noted by [START_REF] Rossi | How many dissimilarity/kernel self organizing map variants do we need?[END_REF], the complexity of kernel SOM is O(n 2 UT) whereas the complexity of the numeric SOM in R p is O(npUT). When p n and n is large, this cost can be very prohibitive since, typically, T is of order O(n). Different strategies have been developed to overcome this difficulty, among which approximations with dimensionality reduction or sparse representations [START_REF] Hofmann | Efficient approximations of robust soft learning vector quantization for non-vectorial data[END_REF][START_REF] Mariette | Efficient interpretable variants of online SOM for large dissimilarity data[END_REF] or exact approaches using a storage of intermediate results [START_REF] Mariette | Accelerating stochastic kernel som[END_REF]. Other alternatives to reduce the complexity of kernel methods directly use accelerated computations of the kernel dot products (used in most kernel methods) with tiled reduction schemes on GPU, without even storing the kernel itself (see KeOps, https://www.kernel-operations.io/keops/index.html). They would be a practicable approach to accelerate the assignment step of kernel SOM when the full kernel matrix itself does not fit in memory but we will restrict to the simpler case where the kernel is already computed and stored, in the remaining of this section.

Low rank and sparse approximations.

The first type of approaches relies on a simpler representation of the prototypes, with a reduced number of (non zero) coefficients. [START_REF] Mariette | Efficient interpretable variants of online SOM for large dissimilarity data[END_REF] proposes two types of solutions. The first one is very similar to the strategies developed in [START_REF] Hofmann | Efficient approximations of robust soft learning vector quantization for non-vectorial data[END_REF] and uses a direct sparse approach in which an additional step is added to each iterations, aiming at thresholding the smallest coefficients (β ui ) i for every prototype p u . The second method relies on a prior step (a kernel PCA) to provide inputs that are the coordinates of the original observations on the first k principal components of the kernel PCA. The SOM algorithm then used is a simple numeric SOM with a complexity of O(nkTU) with k n. However, this prior step has a high computational cost itself: the full eigendecomposition of K has a computational cost of O(n 3 ) but it can be reduced with the Nyström approximation [START_REF] Williams | Using the Nyström method to speed up kernel machines[END_REF]. This method allows to obtain an approximation of the eigendecomposition of K using an eigendecomposition of a submatrix K (m) based on m observations chosen at random in the original sample. The eigendecomposition approximation is even exact if the rank of K is smaller than m. The complexity of the approach is reduced to O(nm 2 ) where m is usually chosen n.

Exact approaches.

Most of the complexity of the kernel SOM comes from the assignment step, which is O(n 2 U). Re-formulating this step and transforming it into the update of stored results, we reduced it to O(U). More precisely, the assignment step writes:

f t+1 (x i ) = argmin u=1,...,U
A t u -2B t ui with A t u = n j, j =1 β t u j β t u j k j j and B t ui = n j=1 β t u j k i j . Storing these quantities in memory (a vector of size U and a (U × n)-matrix), the representation step reduces to an update of A t and B t :

A t+1 u = (1 -λ u (t)) 2 A t u + λ u (t) 2 k ii + 2λ u (t)(1 -λ u (t))B t ui B t+1 ui = (1 -λ u (t))B t ui + λ u (t)k i i , with λ u (t) = µ t H t (d( f t+1 (x i ), u))
. The representation step thus has a complexity of O(nU) (update of B t ) and the total complexity of the approach is reduced to O(nUT). This reduction of the computational time is thus obtained at the cost of storing operations with a memory cost of O(U) and O(nU) for A t and B t respectively.

Combining kernels

Kernel methods have proven to be particularly efficient when data are described by multi-source and multi-type information obtained on the same n observations. In this case, each source of data, of a given particular type (numerical, graph data, factors, ...), can be passed through a kernel, K m (m = 1, . . . , M): this kernel provides the similarity information between observations, seen from the point of view of the source m. The advantage of such an approach is that it provides a common representation of the different sources that can be easily combined. A similar framework is the one where multiple kernels can be obtained from a single dataset, each capturing a specific feature. Combining these kernels avoids having to choose between them, and also benefits of the information coming from different aspects of the data. Among the combination approaches [START_REF] Gönen | Multiple kernel learning algorithms[END_REF], one that has been widely developed is the computation of a convex combination of the M kernels into a single meta-kernel:

K γ = M m=1 γ m K m , st γ m ≥ 0, ∀ m = 1, . . . , M M m=1 γ m = 1 .
In the context of supervised methods, the choice of (γ m ) m is usually done by solving a global optimization problem that aims at minimizing a prediction loss, with respect to the parameter of a given method (SVM for instance) and to the value of (γ m ) m [START_REF] Zhao | Multiple kernel clustering[END_REF][START_REF] Yu | Optimized data fusion for kernel k-means clustering[END_REF][START_REF] Huang | Multiple kernel fuzzy clustering[END_REF][START_REF] Gönen | Localized data fusion for kernel k-means clustering with application to cancer biology[END_REF]. In the unsupervised setting, choosing relevant (γ m ) m is harder since the objective function might not be as easily designed or because, as it is the case for kernel PCA, its joint optimization to estimate the principal components and the (γ m ) m is degenerate [START_REF] Speicher | Towards multiple kernel principal component analysis for integrative analysis of tumor samples[END_REF]. Several propositions have thus been made [START_REF] Lin | Multiple kernel learning for dimensionality reduction[END_REF][START_REF] Zhuang | Unsupervised multiple kernel clustering[END_REF][START_REF] Speicher | Integrating different data types by regularized unsupervised multiple kernel learning with application to cancer subtype discovery[END_REF][START_REF] Wang | Visualization and analysis of single-cell RNA-seq data by kernel-based similarity learning[END_REF][START_REF] Mariette | Unsupervised multiple kernel learning for heterogeneous data integration[END_REF] to overcome this issue and, in the latter, we proposed two solutions that can cope with non numerical observations, contrary to the others. The first method, named STATIS-UMKL, is based on the STATIS method [START_REF] Des Plantes | Structuration des tableaux à trois indices de la statistique[END_REF][START_REF] Lavit | The ACT (STATIS method)[END_REF] and aims at searching for a consensual metakernel. More precisely, the method searches for the kernel that is the most similar, on average, to all the kernels to be combined, (K m ) m=1,..., M :

max v M m=1 K v , K m K m F F for K v = M m=1 v m K m ,
and v ∈ R M such that v 2 R M = 1, where ., . and . F stand for the Frobenius dot product and norm. It is easy to show that the solution is given by the spectral decomposition of a M × M-matrix, C, such

that C mm = K m ,K m F K m F K m F
and γ is thus chosen as v m v m . The second method first creates a proxy of the local geometry induced by each kernel K m using a k nearest neighbor graph and the global adjacency matrix of these M graphs, W is then used in a global criterion. This criterion is designed to preserve at best the local geometry measured by W in the feature space induced by the meta-kernel K γ :

argmin γ ∈R M n i,i =1 W ii C i (γ) -C i (γ) 2 R n , st γ m ≥ 0 and M m=1 γ m = 1, with C i (γ) = φ γ (x i ), φ γ (x 1 ) . . . φ γ (x n ) H γ = K γ (x i , x 1 ) . . . K γ (x i , x n )
.

This problem has a sparse solution that performs a selection of the kernels (some of the entries of (γ m ) m are forced toward 0) because of the convexity constraint M m=1 γ m = 1 but this can be relaxed by replacing the 1 constraint with a constraint of the 2 norm instead (the two versions are called sparse-UMKL and full-UMKL where "UMKL" stands for Unsupervised Multiple Kernel Learning).

Once the kernel is obtained, it can be used as input to kernel based algorithm, like kernel PCA, kernel k-means or kernel SOM for exploratory purpose.

Application

The combined-kernel SOM algorithm is illustrated on data related to school-to-work transitions, extracted from the survey "Generation 98"2. The dataset contains information on 16,040 young people having graduated in 1998 and monitored during 94 months after having left school. The labor-market status has nine categories, labeled as follows: permanent contract, fixed-term contract, apprenticeship, public temporary contract, on-call contract, unemployed, inactive, military service, education. The following stylized facts are highlighted by a first descriptive analysis of the data, as shown in Figure 13: • permanent contracts represent more than 20% of all status after one year and their ratio continues to increase up to 50% after three years and almost 75% after seven years; • the ratio of fixed-term contracts is more than 20% after one year on the labor market, but it is decreasing to 15% after three years and then seems to converge to 8%; • almost 30% of the young graduates are unemployed after one year. This ratio is decreasing and becomes constant, 10%, after the fourth year.

Fig. 1 Chronogram of the labor market structure as illustrated by the "Génération 98" dataset.

x-axis is time (in month) and y-axis is the proportion of each type of contract.

Some trajectories were duplicated (some people had exactly the same job trajectories) so, to reduce redundancy and computational time, we used only the 12,471 unique trajectories.

Three optimal combined kernels were computed from three sets of dissimilarities with different features, and each of these kernels was then used as input to the kernel SOM. The three sets of dissimilarities are described with more details in the next section. For each set of dissimilarities, the optimal combined kernel was obtained as follows:

• first, each dissimilarity matrix, D, was transformed into a (centered) similarity matrix by computing

∀ i, j = 1, . . . , n, s(x i , x j ) = - 1 2 D 2 i j - 1 n n l=1 D 2 il - 1 n n l=1 D 2 jl + 1 n 2 n l,l =1 D 2 ll .
Each of these resulting similarity matrices were used as kernels, even though a small part of their spectra was non positive; • second, the resulting kernel matrices were optimally combined using STATIS-UMKL, as described in Section 5.

All dissimilarities were computed using the R package TraMineR [START_REF] Gabadinho | Analyzing and visualizing state sequences in R with TraMineR[END_REF] and the combination of kernels was obtained using the R package mixKernel. Each combined kernel was processed through a kernel self-organizing map using the implementation provided in the R package SOMbrero. For each map to be trained, a 10 × 10 configuration was selected and default values of the package were chosen for the initialization step, the topology of the map (choice of H t in Algorithm 1) and the decreasing value µ t (also as in Algorithm 1). We decided to use 10×10 maps as a trade-off between having a meaningful visualization and a reduced number of meaningful typical trajectories. Since many trajectories are redundant -the permanent contracts are overrepresented-a map with a number of units equal to about a tenth of the number of inputs was a reasonable choice. Final results were represented as chronograms: more precisely, each unit of the map was featured by a rectangle containing the chronogram of the subsample of trajectories assigned to this unit.

Three sets of dissimilarities and relations between them

There is currently a vast literature devoted to measuring similarities for longitudinal data, and the community agrees that the differences between the various criteria focus on three different aspects of sociological importance: the sequencing or the order in which the states appear, the timing, and the duration of the states. A recent and detailed review of these methods, focusing on these different aspects, and also introducing some new and versatile criteria, is available in [START_REF] Studer | What matters in differences between life trajectories: a comparative review of sequence dissimilarity measures[END_REF]. Starting from these considerations, three sets of distances were selected for the present study: the first aimed at focusing on the sequencing and possibly the duration, the second on the timing, and the third on the duration only.

The first group of dissimilarities contains one criterion based on the number of matching subsequences, and two based on generalizations of the classical OM metric [START_REF] Needleman | A general method applicable to the search for similarities in the amino acid sequence of two proteins[END_REF][START_REF] Abbott | Optimal matching methods for historical sequences[END_REF]:

• The SVRspell distance, proposed by [START_REF] Elzinga | Spell sequences, state proximities, and distance metrics[END_REF], is computed using the number of matching subsequences within the distinct sequences of states, where the durations of the spells are weighted by some parameter b. This method is built to be sensitive to sequencing, and, depending on b, it may be also made sensitive to durations. In the following, we set b = 0, so that the sensitivity to the order only is put forward. • The OMspell distance, introduced in [START_REF] Studer | What matters in differences between life trajectories: a comparative review of sequence dissimilarity measures[END_REF], generalizes the OM distance to sequences of spells. The increasing size of the alphabet and of the costs to be specified are controlled through a linear function depending on a parameter δ, representing the cost of extending or compressing a spell by one unit of time. For small values of δ, the method favors the expansion or the compression of existing spells. For δ = 0, OMspell reduces to the usual OM distance between the distinct sequences of states. Let us remark here that the smaller δ, the less sensitive the criterion is to the duration of the spells and the more it is to the sequencing. • The OMstran distance, also introduced in [START_REF] Studer | What matters in differences between life trajectories: a comparative review of sequence dissimilarity measures[END_REF], adapts the OM distance to sequences of transitions. Here also, the alphabet and the number of costs to be specified are much larger than in the usual setting, but the dimension of the parameters is reduced by considering a convex combination between the costs of the spells and the transition costs, controlled by some w. In the following, the value of w was fixed so as to favor a criterion sensitive to differences in sequencing.

The second set of dissimilarities was focused toward highlighting timing. Four distances were tested: the Hamming distance (HAM), based on the number of non-matching states, the Euclidean distance (EUCLID), which, in this case, is the squared root of the Hamming, the χ 2 (CHI2), which is similar to the two previous dissimilarities except that it gives more weights to the infrequent states, and the Dynamic Hamming distance (DHD). For both EUCLID and CHI2 distances, the sensitivity between duration and timing is controlled through a parameter, L. When L = 94 (i.e., in our case the trajectory length), scores are similar to those of the Hamming family regarding timing, when L = 1, dissimilarity measures are more sensitive to duration. DHD, proposed by [START_REF] Lesnard | Setting cost in optimal matching to uncover contemporaneous socio-temporal patterns[END_REF], generalizes the Hamming distance by considering an OM dissimilarity without insertions or deletions, and with the substitution costs defined at each temporal instant from the corresponding transition matrices. While taking the position in the sequence and time into account, DHD has often been criticized for its risk of over-parameterization. In the following, distance parameters used are specified between parentheses.

The third group of dissimilarities was defined to be sensitive to duration of states. Again, four distances were selected: Euclidean and χ 2 distances between state distributions in the whole trajectories, OMspell with δ = 1, which is sensitive both to sequencing and duration, and a distance based on the length of the longest common subsequence of two trajectories, LCS, as described in [START_REF] Bergroth | A survey of longest common subsequence algorithms[END_REF]. Whereas the Euclidean distance is more sensitive to differences between states with a high duration, the χ 2 gives more importance to rare states.

Figure S1 of Supplementary material illustrates the relations between the different distances on a cosine matrix (computed from the Frobenius dot product) for the school-to-work transition dataset. Within the first group of distances, SVRspell with b = 0 appears to behave very specifically, and is even very different from the other distances of its own group, except for the OMspell with δ = 0. This is explained by the fact that SVRspell is sensitive to the sequencing only, and not at all to timing or duration. Furthermore, since it is based on the number of common subsequences only, its principle and computation are quite different from its OM-based counterparts.

The only distance highly correlated to SVRspell is the OMspell with δ = 0, which is also sensitive to sequencing only. The OMspell with δ = 0 is also very different from the distances in its own group, except for OMspell with δ = 0.1 and OMstran. These two distances introduce some sensitivity with respect to the duration of the spells, while still mainly favoring sequencing. These results are consistent with the conclusions in [START_REF] Studer | What matters in differences between life trajectories: a comparative review of sequence dissimilarity measures[END_REF], based on simulated data.

In the second group, the four distances are all very similar, and more particularly HAM, Euclid and DHD. The χ 2 -distance stands out because of its particular weighting. We can also note that all four distances in Group 2 are also very similar to OMstran and OMspell with δ = 0.1, which favor sequencing and duration, and very different from SVRspell and OMspell with δ = 0, which favor sequencing only.

Within the third group of distances, all distances are also very similar, even if the value of the cosine is a bit lower than within distances in the second group. The distances in the third group are also similar to the ones in the second group, to OMspell with δ = 0.1 and to OMstran in the first group.

In conclusion, SVRspell and OMspell with δ = 0 are very different from all other distances, which are globally similar. This indicates that, for this dataset, criteria favoring sequencing are quite opposite to those favoring timing or duration.

Results of the clusterings

The map obtained with the first four distances is provided in Figure 2. For the sake of conciseness, the other two maps are available in Figures S2 andS3 and in Section S2 of the Supplementary material.

As one may easily see, most of the clusters show smooth sigmoidal transitions between states, and in some cases "sandwich"-like representations. These patterns are inherent to the fact that the first set of distances was built to point out similarities in terms of sequencing and not in terms of timing or duration. A chronogram representation is not a well suited representation in this case because it allows for temporal shifts. Despite that, some clusters of particular interest can be identified, such as the lower-left corner of the map, showing the outcomes of public fixed-term contracts. It gives more weight to the OM-inspired distances, which appear to be more discriminant than the SVRspell. As it has already been stressed in the comparison between distances (Figure S1 of the Supplementary material), the SVRspell distance behaves very differently from all the other distances and, according to [START_REF] Studer | What matters in differences between life trajectories: a comparative review of sequence dissimilarity measures[END_REF], it is very sensitive to sequencing and small random perturbations. The latter may be a reason for which the other dissimilarities appear to be fitter for clustering trajectories based on their sequencing properties.

Concluding remarks

This chapter has presented a global approach to perform exploratory analysis in the presence of multiple sources of data or of multiple kernels describing different features of the data. The approach has been illustrated on a dataset of categorical time series representing labor-market status of recently graduated people. We have shown the effectiveness of the approach to identify a relevant typography of the dataset, with contrasting results depending on which features the focus is put on. Different dissimilarities led to highlight different characteristics of the trajectories, some less suited to chronogram representations than others. To fully exploit that diversity, alternative representations would be needed, which could be able to better represent similar duration of states or similar global distributions of the trajectories and thus to highlight the distance features.

Fig. 2

 2 Fig. 2 Final map obtained with the first group of distances.

Available thanks to Génération 1998 à 7 ans -2005, [producer] CEREQ, [diffusion] Centre Maurice Halbwachs (CMH).

The graphical illustrations were carried out using the TraMineR package[START_REF] Gabadinho | Analyzing and visualizing state sequences in R with TraMineR[END_REF].

Acknowledgements. Nathalie Vialaneix would like to sincerely thank Professor Christine Thomas-Agnan for so many interesting discussions and for making her participate to the organization of "Journées de Statistique de la SFdS 2013" (an unforgettable adventure!). From a scientific point of view, Christine's book (written with Alain Berlinet) has been her main starting point for her work on kernel methods and use of spline based regularizations in functional data analysis. For all this, I am more than grateful!