Rubén Ibáñez

Emmanuelle Abisset-Chavanne

Amine Ammar

David González

El-As Cueto
email: ecueto@unizar.es

Antonio Huerta

Jean Louis Duval

Francisco Chinesta

A Multidimensional Data-Driven Sparse Identification Technique: The Sparse Proper Generalized Decomposition

Sparse model identification by means of data is especially cumbersome if the sought dynamics live in a high dimensional space. This usually involves the need for large amount of data, unfeasible in such a high dimensional settings. This well-known phenomenon, coined as the curse of dimensionality, is here overcome by means of the use of separate representations. We present a technique based on the same principles of the Proper Generalized Decomposition that enables the identification of complex laws in the lowdata limit. We provide examples on the performance of the technique in up to ten dimensions.

Introduction

In recent years there has been a growing interest in incorporating data-driven techniques into the field of mechanics. While almost classical in other domains of science like economics, sociology, etc., big data has arrived with important delay to the field of computational mechanics. It is worth noting that, in our field, the amount of data available is very often no so big, and therefore we speak of data-driven techniques instead of big-data techniques.

Among the first in incorporating data-driven technologies to the field of computational mechanics one can cite the works of Kirchdoerfer et al. [START_REF] Kirchdoerfer | Data-driven computational mechanics[END_REF][START_REF] Kirchdoerfer | Data driven computing with noisy material data sets[END_REF], or the ones by Brunton et al. [START_REF] Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF][START_REF] Brunton | Sparse identification of nonlinear dynamics with control (sindyc)[END_REF][START_REF] Quade | Sparse identification of nonlinear dynamics for rapid model recovery[END_REF]. Previous attempts exist; however, to construct data-driven identification algorithms, see, for instance [START_REF] Gonzalez | Proper generalized decomposition based dynamic data driven inverse identification[END_REF][START_REF] Peherstorfer | Data-driven operator inference for nonintrusive projection-based model reduction[END_REF]. More recently, the issue of compliance with general laws like the ones of thermodynamics has been also achieved, which is a distinct feature of data-driven mechanics [START_REF] González | Thermodynamically consistent data-driven computational mechanics[END_REF]. Other applications include the identification of biological systems [START_REF] Mangan | Inferring Biological Networks by Sparse Identification of Nonlinear Dynamics[END_REF] or financial trading [START_REF] Mann | Dynamic mode decomposition for financial trading strategies[END_REF], to name but a few.

The problem with high dimensional systems is that data in these systems is often sparse (due precisely to the high dimensional nature of the phase space) while the system has, on the contrary, low dimensional features-at least very frequently. Based on this, a distinction should be made between methods that require an a priori structure of the sampling points and others which do not require such a regularity.

Regarding the methods that need a rigid structure in the sampling points, the Nonintrusive Sparse Subspace Learning (SSL) method is a novel technique which has proven to be very effective [START_REF] Borzacchiello | Non-intrusive Sparse Subspace Learning for Parametrized Problems[END_REF]. The basic ingredient behind such a technique is that the parametric space is explored in a hierarchical manner, where sampling points are collocated at the Gauss-Lobato-Chebychev integration points. Also, using a hierarchical base allows improving the accuracy adding more hierarchical levels without perturbing the previous ones. To achieve such hierarchical property, just the difference at a given point between the real function minus the estimated value, using the precedent hierarchical levels, is propagated. For more details about the method, the reader is referred to [START_REF] Borzacchiello | Non-intrusive Sparse Subspace Learning for Parametrized Problems[END_REF]. However, in the high-dimensional case, this technique shows severe limitations, as will be detailed hereafter.

On the other hand, nonstructured data-driven techniques are commonly based on Delaunay triangularization techniques, providing an irregular mesh whose nodes coincides with the sampling points. Afterwards, depending on the degree of approximation inside each one of the Delaunay triangles, it gives rise to different interpolation techniques; i.e., linear, nearest, cubic, and natural, among other techniques, are commonly used. Apart from techniques that depend on a given triangularization, it is worth mentioning Kriging interpolants as an appealing technique to provide response surfaces from nonstructured data points. The key ingredient behind such technique is that each sampling point is considered as a realization of a random process. Therefore, defining a spatial correlation function allows to infer the position of unknown points just like providing confidence intervals based on the distance to the measured points. Nevertheless, the calibration of the correlation matrix has an important impact in the performance of the method itself.

Kriging also possesses a very interesting property: it is able to efficiently filter noise and outliers. Therefore, it is expected that it also could help us in problems with noise in the data.

However, in high dimensional settings, all of the just mentioned techniques fail to identify the nature of the system due precisely to the curse of dimensionality. A recent alternative for such a system could be Topological Data Analysis (TDA), which is based on the use of algebraic topology and the concept of persistent homology [START_REF] Epstein | Topological data analysis[END_REF]. A sparse version of this technique also exists [START_REF] Guo | Sparse-TDA: Sparse Realization of Topological Data Analysis for Multi-Way Classification[END_REF].

Hence, if a competitive data-driven identification technique is desired, such a technique should meet the following requirements:

(i) Nonstructured data set: this characteristic provides versatility to the method. Indeed, when evaluating the response surface requiring a lot of computational effort, recycling previous evaluations of the response surface, which do not coincide with a given structure of the data, may be very useful. In addition, the SSL technique establishes sampling points at locations in the phase space with no physical meaning in an industrial setting.

(ii) Robustness with respect to high dimensionality: triangularization-based techniques suffer when dealing with multidimensional data just because a high dimensional mesh has to be generated. Nevertheless, the separation of variables could be an appealing technique to circumvent the problem of generating such a high dimensional mesh.

(iii) Curse of dimensionality: all previous techniques suffer when dealing with high dimensional data. For instance, the SSL needs 2 𝐷 sampling points just to reach the first level of approximation. Thus, when dealing with high dimensional data (𝐷 > 10 uncorrelated dimensions) plenty of sampling points are required to properly capture a given response surface.

In what follows we present a method based on the concept of separate representations to overcome the curse of dimensionality. Such separate representation has previously been employed by the authors to construct a priori reducedorder modeling techniques, coined as Proper Generalized Decompositions [START_REF] Chinesta | Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models[END_REF][START_REF] Chinesta | e proper generalized decomposition for advanced numerical simulations[END_REF][START_REF]Separated representations and PGD-based model reduction[END_REF][START_REF] Gonzalez | Recent advances on the use of separated representations[END_REF][START_REF] Badas | Local proper generalized decomposition[END_REF][START_REF] Cueto | Proper generalized decompositions[END_REF][START_REF] González | Model order reduction for real-time data assimilation through extended Kalman filters[END_REF]. This will give rise to a sparse Proper Generalized Decomposition (s-PGD in what follows) approach to the problem. We then analyze the just developed technique through a series of numerical experiments in Section 4, showing the performance of the method. Examples in up to ten dimensions are shown. The paper is completed with some discussions.

A Sparse PGD (s-PGD) Methodology

. . Basics of the Technique. In this section we develop a novel methodology for sparse identification in high dimensional settings. For the ease of the exposition and, above all, representation, but without loss of generality, let us begin by assuming that the unknown objective function 𝑓(𝑥, 𝑦) lives in R 2 and that is to be recovered from sparse data. As in previous references, see, for instance [START_REF] Loiseau | Constrained sparse Galerkin regression[END_REF]; we have chosen to begin with a Galerkin projection, in the form

∫ Ω 𝑤 * (𝑥, 𝑦) (𝑢 (𝑥, 𝑦) -𝑓 (𝑥, 𝑦)) 𝑑𝑥 𝑑𝑦 = 0, (1)
where Ω ⊂ R 2 stands for the-here, still twodimensional-domain in which the identification is performed and 𝑤 * (𝑥, 𝑦) ∈ C 0 (Ω) is an arbitrary test function. Finally, 𝑢(𝑥, 𝑦) will be the obtained approximation to 𝑓(𝑥, 𝑦), still to be constructed. In previous works of the authors [START_REF] González | Thermodynamically consistent data-driven computational mechanics[END_REF] as well as in other approaches to the problem (e.g., [START_REF] Loiseau | Constrained sparse Galerkin regression[END_REF]), this projection is subject to additional constraints of thermodynamic nature. In this work no particular assumption is made in this regard, although additional constraints could be imposed to the minimization problem.

Following the same rationale behind the Proper Generalized Decomposition (PGD), the next step is to express the approximated function 𝑢 𝑀 (𝑥, 𝑦) ≈ 𝑢(𝑥, 𝑦) as a set of separate one-dimensional functions,

𝑢 𝑀 (𝑥, 𝑦) = 𝑀 ∑ 𝑘=1 𝑋 𝑘 (𝑥) 𝑌 𝑘 (𝑦) . (2)
The determination of the precise form of functional pairs 𝑋 𝑘 (𝑥)𝑌 𝑘 (𝑦), 𝑘 = 1, . . . , 𝑀, is done by first projecting them on a finite element basis and by employing a greedy algorithm such that once the approximation up to order 𝑀-1 is known, the new 𝑀th order term

𝑢 𝑀 (𝑥, 𝑦) = 𝑢 𝑀-1 (𝑥, 𝑦) + 𝑋 𝑀 (𝑥) 𝑌 𝑀 (𝑦) = 𝑀-1 ∑ 𝑘=1 𝑋 𝑘 (𝑥) 𝑌 𝑘 (𝑦) + 𝑋 𝑀 (𝑥) 𝑌 𝑀 (𝑦) , (3)
is found by any nonlinear solver (Picard, Newton,. . .).

It is well-known that this approach produces optimal results for elliptic operators (here, note that we have in fact an identity operator acting on 𝑢) in two dimensions, see [START_REF] Chinesta | Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models[END_REF] and references therein. There is no proof, however, that this separate representation will produce optimal results (in other words, will obtain parsimonious models) in dimensions higher than two. In two dimensions and with 𝑤 * = 𝑢 * it provides the singular value decomposition of 𝑓(𝑥, 𝑦) [START_REF] Chinesta | e proper generalized decomposition for advanced numerical simulations[END_REF]. Our experience, nevertheless, is that it produces almost optimal results in the vast majority of the problems tested so far.

It is worth noting that the product of the test function 𝑤 * (𝑥, 𝑦) times the objective function 𝑓(𝑥, 𝑦) is only evaluated at few locations (the ones corresponding to the experimental measurements) and that, in a general high dimensional setting, we will be in the low-data limit necessarily. Several options can be adopted in this scenario. For instance, the objective function can be first interpolated in the high dimensional space (still 2D in this introductory example) and then integrated together with the test function. Indeed, this will be the so-called PGD in approximation [START_REF] Chinesta | e proper generalized decomposition for advanced numerical simulations[END_REF], commonly used when either 𝑓(𝑥, 𝑦) is known everywhere and a separated representation is sought or if 𝑓(𝑥, 𝑦) is known in a separated format but a few pairs 𝑀 are needed for any reason. Under this rationale the converged solution 𝑢(𝑥, 𝑦) tries to capture the already interpolated solution in the high dimensional space but in a more compact format. As a consequence, the error due to interpolation of experimental measurements on the high dimensional space will persist in the final separate identified function.

In order to overcome such difficulties, we envisage a projection followed by interpolation method. However since information is just known at 𝑃 sampling points (𝑥 𝑖 , 𝑦 𝑖 ), 𝑖 = 1, . . . , 𝑃, it seems reasonable to express the test function not in a finite element context, but to express it as a set of Dirac delta functions collocated at the sampling points,

𝑤 * (𝑥, 𝑦) = 𝑢 * (𝑥, 𝑦) 𝑃 ∑ 𝑖=1 𝛿 (𝑥 𝑖 , 𝑦 𝑖 ) = (𝑋 * (𝑥) 𝑌 𝑀 (𝑦) + 𝑋 𝑀 (𝑥) 𝑌 * (𝑦)) 𝑃 ∑ 𝑖=1 𝛿 (𝑥 𝑖 , 𝑦 𝑖 ) , (4)
giving rise to

∫ Ω 𝑤 * (𝑥, 𝑦) (𝑢 (𝑥, 𝑦) -𝑓 (𝑥, 𝑦)) 𝑑𝑥 𝑑𝑦 = ∫ Ω 𝑢 * (𝑥, 𝑦) 𝑃 ∑ 𝑖=1 𝛿 (𝑥 𝑖 , 𝑦 𝑖 ) (𝑢 (𝑥, 𝑦) -𝑓 (𝑥, 𝑦)) 𝑑𝑥 𝑑𝑦 = 0, (5)
The choice of the test function 𝑤 * (𝑥, 𝑦) in the form dictated by ( 4) is motivated by the desire of employing a collocation approach while maintaining the symmetry of standard Bubnov-Galerkin projection operation.

. . Matrix Form. Let us detail now the finite element projection of the one-dimensional functions 𝑋 𝑘 (𝑥) and 𝑌 𝑘 (𝑦), 𝑘 = 1, . . . , 𝑀, (often referred to as modes) appearing in [START_REF] Kirchdoerfer | Data driven computing with noisy material data sets[END_REF]. Several options can be adopted, ranging from standard piecewise linear shape functions, global nonlinear shape functions, maximum entropy interpolants, splines, kriging, etc. Regarding the kind of interpolant to use, an analysis will be performed in the sequel. Nevertheless, no matter which precise interpolant is employed, it can be expressed in matrix form as

𝑋 𝑘 (𝑥) = 𝑁 ∑ 𝑗=1 𝑁 𝑘 𝑗 (𝑥) 𝛼 𝑘 𝑗 = [𝑁 𝑘 1 (𝑥) . . . 𝑁 𝑘 𝑁 (𝑥)] [ [ [ [ [ 𝛼 𝑘 1 . . . 𝛼 𝑘 𝑁 ] ] ] ] ] = (N 𝑘 𝑥 ) 𝑇 a 𝑘 , (6)
𝑌 𝑘 (𝑦) = 𝑁 ∑ 𝑗=1 𝑁 𝑘 𝑗 (𝑦) 𝛽 𝑘 𝑗 = [𝑁 𝑘 1 (𝑦) . . . 𝑁 𝑘 𝑁 (𝑦)] [ [ [ [ [ 𝛽 𝑘 1 . . . 𝛽 𝑘 𝑁 ] ] ] ] ] = (N 𝑘 𝑦 ) 𝑇 b 𝑘 , (7)
where 𝛼 𝑘 𝑗 and 𝛽 𝑘 𝑗 , 𝑗 = 1, . . . , 𝑁, represent the degrees of freedom of the chosen approximation. We employ N k as the most usual nomenclature for the shape function vector. It is important to remark that the approximation basis could even change from mode to mode (i.e., for each 𝑖). For the sake of simplicity we take the same number of terms for both 𝑋 𝑘 (𝑥) and 𝑌 𝑘 (𝑦), namely, 𝑁.

By combining (1), ( 2), ( 4), [START_REF] Gonzalez | Proper generalized decomposition based dynamic data driven inverse identification[END_REF], and (7) a nonlinear system of equations is derived, due to products of terms in both spatial directions. An alternate direction scheme is here preferred to linearize the problem, which is also a typical choice in the PGD literature. Note that, when computing modes 𝑋 𝑀 (𝑥), the variation in the other spatial direction vanishes, 𝑌 * (𝑦) = 0, and vice versa.

In order to fully detail the matrix form of the resulting problem, we first employ the notation "⊗" as the standard tensorial product (i.e., b ⊗ c = 𝑏 𝑖 𝑐 𝑗 ) and define the following matrices

A 𝑘ℓ 𝑥 = N 𝑘 𝑥 ⊗ N ℓ 𝑥 , A 𝑘ℓ 𝑦 = N 𝑘 𝑦 ⊗ N ℓ 𝑦 , C 𝑘ℓ 𝑥𝑦 = N 𝑘 𝑥 ⊗ N ℓ 𝑦 . (8)
For the sake of simplicity but without loss of generality, evaluations of the former operators at point (𝑥 𝑖 , 𝑦 𝑖 ) are denoted as

A 𝑘ℓ 𝑥 𝑖 = N 𝑘 𝑥 (𝑥 𝑖 ) ⊗ N ℓ 𝑥 (𝑥 𝑖 ) , A 𝑘ℓ 𝑦 𝑖 = N 𝑘 𝑦 (𝑦 𝑖 ) ⊗ N ℓ 𝑦 (𝑦 𝑖 ) , C 𝑘ℓ 𝑥 𝑖 𝑦 𝑖 = N 𝑘 𝑥 (𝑥 𝑖 ) ⊗ N 𝑗 𝑦 (𝑦 𝑖 ) . (9)
Equations ( 10)-( 11) below show the discretized version of the terms appearing in the weak form, (1), when computing modes in the 𝑥 direction. Again, 𝑀 stands for the number of modes in the solution 𝑢(𝑥, 𝑦) while 𝑃 denotes the number of sampling points.

∫ Ω 𝑢 * (𝑥, 𝑦) 𝑃 ∑ 𝑖=1 𝛿 (𝑥 𝑖 , 𝑦 𝑖 ) 𝑢 (𝑥, 𝑦) 𝑑𝑥 𝑑𝑦 = 𝑀 ∑ 𝑘=1 𝑃 ∑ 𝑖=1 ((b 𝑀 ) 𝑇 A 𝑀𝑘 𝑦 𝑖 b 𝑘 ) ((a * ) 𝑇 A 𝑀𝑘 𝑥 𝑖 a 𝑘 ) , (10)
∫ Ω 𝑢 * (𝑥, 𝑦) 𝑃 ∑ 𝑖=1 𝛿 (𝑥 𝑖 , 𝑦 𝑖 ) 𝑓 (𝑥, 𝑦) 𝑑𝑥 𝑑𝑦 = 𝑃 ∑ 𝑖=1 𝑓 (𝑥 𝑖 , 𝑦 𝑖 ) ((a * ) 𝑇 C 𝑀𝑀 𝑥 𝑖 𝑦 𝑖 b 𝑀 ) . (11)
Hence, by defining

M 𝑥 = 𝑃 ∑ 𝑖=1 ((b 𝑀 ) 𝑇 A 𝑀𝑀 𝑦 𝑖 b 𝑀 ) A 𝑀𝑀 𝑥 𝑖 , m 𝑥 = 𝑀-1 ∑ 𝑘=1 𝑃 ∑ 𝑖=1 ((b 𝑀 ) 𝑇 A 𝑀𝑘 𝑦 𝑖 b 𝑘 ) A 𝑀𝑘 𝑥 𝑖 a 𝑘 , f 𝑥 = 𝑃 ∑ 𝑖=1 𝑓 (𝑥 𝑖 , 𝑦 𝑖 ) C 𝑀𝑀 𝑥 𝑖 𝑦 𝑖 b 𝑀 , (12)
allows to write a system of algebraic equations

M 𝑥 a 𝑀 = f 𝑥 -m 𝑥 . ( 13
)
Exactly the same procedure is followed to obtain an algebraic system of equations for b 𝑀 . This allows performing an alternating directions scheme to extract a new couple of 𝑋 𝑀 (𝑥) and 𝑌 𝑀 (𝑦) modes.

This formulation has several aspects that deserve to be highlighted: [START_REF] Kirchdoerfer | Data-driven computational mechanics[END_REF] No assumption about 𝑓(𝑥, 𝑦) has been made other than assuming known its value at sampling points. Indeed, both problems of either interpolating or making a triangulation in a high dimensional space are circumvented due to the separation of variables.

(2) The operator M 𝑥 is composed of 𝑃 rank-one updates, meaning that the rank of such operator is at most 𝑃. Furthermore, if a subset of measured points share the same coordinate 𝑥 𝑖 , the entire subset will increase the rank of the operator in one unity.

(3) The position of the sampling points will constrain the rank of the PGD operators. That is the reason why even if the possibility of having a random sampling of points is available, it is always convenient to perform a smart sampling technique such that the rank in each direction tends to be maximized. Indeed, the higher the rank of the PGD operator is, the more cardinality of a and b can be demanded without degenerating into an underdetermined system of equations.

There are plenty of strategies to smartly select the position of the sampling points. They are based on either knowing an a priori error indicator or having a reasonable estimation of the sought response surface. Certainly, an adaptive strategy based on the gradient of the precomputed modes could be envisaged. However, the position of the new sampling points will depend on the response surface calculated using the previous sampling points, making parallelization difficult. That is the reason why Latin hypercube is chosen in the present work. Particularly, Latin hypercube tries to collocate 𝑃 sampling points in such a way that the projection of those points into 𝑥 and 𝑦 axis are as far as possible.

. . Choice of the D Basis. In the previous section, nothing has been specified about the basis in which each one of the one-dimensional modes was expressed. In this subsection, we will use an interpolant based on Kriging techniques. Simple Kriging has been used throughout history in order to get relatively smooth solutions, avoiding spurious oscillations characteristic of high order polynomial interpolation. This phenomena is called Runge's phenomenon. It appears due to the fact that the sampling point locations are not chosen properly; i.e., they will not be collocated, in general, at the Gauss-Lobato-Chebychev quadrature points. Kriging interpolants consider each point as a realization of a Gaussian process, so that high oscillations are considered as unlikely events.

Hence, by defining a spatial correlation function based on the relative distance between two points, D(𝑥 𝑖 -𝑥 𝑗 ) = D 𝑖𝑗 , an interpolant is created over the separated 1D domain,

𝑋 𝑘 (𝑥) = 𝑁 ∑ 𝑖=1 𝛼 𝑘 𝑖 𝑁 + 𝑁 ∑ 𝑗=1 𝜆 (𝑥 -𝑥 𝑗 ) (𝛼 𝑘 𝑗 - 𝑁 ∑ 𝑙=1 𝛼 𝑘 𝑙 𝑁 ) , (14)
where 𝜆(𝑥 -𝑥 𝑗 ) is a weighting function which strongly depends on the definition of the correlation function and the 𝛼 𝑖 coefficients are the nodal values associated to the 𝑥 𝑖 Kriging control points. Note that these control points are not the sampling points. We have chosen this strategy so as to allow us to accomplish an adaptivity strategy that will be described next. In the present work, these control points are uniformly distributed along the 1D domain. Although several definitions of the correlation function exist, a Gaussian distribution is chosen as

D 𝑖𝑗 = D (𝑥 𝑖 -𝑥 𝑗 ) = 1 𝜎 √ 2𝜋 𝑒 -(𝑥 𝑖 -𝑥 𝑗 ) 2 /2𝜎 2 , ( 15
)
where 𝜎 is the variance of the Gaussian distribution. Several a priori choices can be adopted to select the value of the variance based on the distance between two consecutive control points, e.g., 𝜎 = ℎ √ (𝑥 𝑖+1 -𝑥 𝑖 ) 2 . The magnitude of ℎ should be adapted depending on the desired global character of the support. To ensure the positivity of the variance, ℎ should be in the interval ]0, +∞[. Let us define now a set of 𝐶 control points

x 𝑐𝑝 = [𝑥 𝑐𝑝 1 , 𝑥 𝑐𝑝 2 , . . . , 𝑥 𝑐𝑝 𝐶 ] (16)
and the 𝑃 sampling points

x 𝑠𝑝 = [𝑥 𝑠𝑝 1 , 𝑥 𝑠𝑝 2 , . . . , 𝑥 𝑠𝑝 𝑃 ] . (17)
Let us define in turn a correlation matrix between all control points and a correlation matrix between the control points and the sampling points as

C 𝑐𝑝-𝑐𝑝 𝑖𝑗 = D (𝑥 𝑐𝑝 𝑖 -𝑥 𝑐𝑝 𝑗 ) , C 𝑐𝑝-𝑠𝑝 𝑖𝑗 = D (𝑥 𝑐𝑝 𝑖 -𝑥 𝑠𝑝 𝑗 ) . (18)
Under these settings, we define a weighting function for each control point and for each sampling point as

Λ = (C 𝑐𝑝-𝑐𝑝 ) -1 C 𝑐𝑝-𝑠𝑝 , (19)
where 𝜆(𝑥

𝑐𝑝 𝑖 -𝑥 𝑠𝑝 𝑗 ) = Λ 𝑖𝑗 .
If we reorganize the terms in the same way that we did in the previous section to have a compact and close format of the shape function N 𝑘 𝑥 , we arrive to

𝑋 𝑘 (𝑥 𝑠𝑝 𝑗 ) = 𝑁 ∑ 𝑖=1 N 𝑘 𝑖 (𝑥 𝑠𝑝 𝑗 ) 𝛼 𝑘 𝑖 = [𝑁 𝑘 1 (𝑥 𝑠𝑝 𝑗 ) . . . 𝑁 𝑘 𝑁 (𝑥 𝑠𝑝 𝑗 )] [ [ [ [ [ 𝛼 𝑘 1 . . . 𝛼 𝑘 𝑁 ] ] ] ] ] = (N 𝑘 𝑥 𝑠𝑝 𝑗 ) 𝑇 a 𝑘 , (20)
where each shape function is given by

𝑁 𝑘 𝑖 (𝑥 𝑠𝑝 𝑗 ) = 1 -∑ 𝑁 𝑗=1 Λ 𝑖𝑗 𝑁 + Λ 𝑖𝑗 . (21)
Figures 1 and2 depict the appearance of the simple Kriging interpolants using 7 control points uniformly distributed along the domain, for ℎ = 1 and ℎ = 1/3, respectively. It can be highlighted that both the Kronecker delta (i.e., strict interpolation) and the partition of unity properties are satisfied for any value of ℎ. Moreover, it is worth noting that the higher the variance the correlation function has, the more global the shape functions are. Furthermore, it is known that 99 per cent of the probability of a Gaussian distribution is comprised within an interval of [𝑚 -3𝜎, 𝑚 + 3𝜎], being 𝑚 the mean value of the distribution. This issue explains perfectly well why the support of each Gaussian distribution takes 2 elements for the case, where ℎ = 1/3. Indeed, the shape of the interpolants is quite similar to standard finite element shape functions, but with a Gaussian profile. The remaining 1 per cent of probability is comprised in the small ridges happening in the middle of the elements.

In light of these results, a family of interpolants based on Kriging can be easily created just selecting the value of the variance within the correlation function. Therefore, globality of the support can be easily adjusted always under the framework of the partition of unity. . . Modal Adaptivity Strategy. In a standard PGD framework, the final solution is approximated as a sum of 𝑀 modes or functional products; see [START_REF] Kirchdoerfer | Data driven computing with noisy material data sets[END_REF]. Each one of the separated modes must be projected onto a chosen basis to render the problem finite dimensional. A standard choice is to select the same basis for each one of the modes:

N 1 = N 2 = ⋅ ⋅ ⋅ = N 𝑀 . (22)
Despite of the fact that this choice seems reasonable, when dealing with nonstructured sparse data, it may not be such. In the previous section we proved that the rank of the separated system strongly depends on the distribution of the data sampling. Therefore, the cardinality of the interpolation basis must not exceed the maximum rank provided by the data sampling. Indeed, this constraint, which provides an upper bound to build the interpolation basis, only guarantees that the minimization is satisfied at the sampling points, without saying anything out of the measured points. Hence, if sampling points are not abundant, in the limit of low-data regime, high oscillations may appear out of these measured points. These oscillations are not desirable since the resulting prediction properties of the proposed method could be potentially decimated.

In order to tackle this problem, we take advantage of the residual-based nature of the PGD. Indeed, the greedy PGD algorithm tries to enrich a solution composed by 𝑀 modes,

𝑢 𝑀 (𝑥, 𝑦) = 𝑀 ∑ 𝑘=1 𝑋 𝑘 (𝑥) 𝑌 𝑘 (𝑦) , (23)
just by looking at the residual that accounts for the contribution of the previous modes, as shown in 8).

Therefore, an appealing strategy to minimize spurious oscillations out of the sampling points is to start the PGD algorithm looking for modes with relatively smooth basis (for instance, Kriging interpolants with a few control points). Therefore, an indicator in order to make an online modal adaptive strategy is required. In the present work, we use the norm of the PGD residual,

R 𝑀 P = 1 √ 𝑃 √ ∑ 𝑖∈P (𝑓 (𝑥 𝑖 , 𝑦 𝑖 ) -𝑢 𝑀 (𝑥 𝑖 , 𝑦 𝑖 )) 2 , ( 24
)
where P is the set of 𝑃 measured points and 𝑓(𝑥, 𝑦) is the function to be captured.

In essence, when the residual norm stagnates, a new control mesh is defined, composed by one more control point and always uniformly distributed, following

ΔR 𝑀 P = R 𝑀 P -R 𝑀-1 P < 𝜖 𝑟 . ( 25
)
By doing this, oscillations are reduced, since higher-order basis will try to capture only what remains in the residual. Here, 𝜖 𝑟 is a tolerance defining the resilience of the sPGD to increase the cardinality of the interpolation basis. The lower 𝜖 𝑟 is, the more resilient method is to increase the cardinality.

To better understand the method, we will quantify the error for two set of points: the first set is associated with the sampling points, P,

E P = 1 #P ∑ 𝑠∈P √ (𝑓 (𝑥 𝑠 , 𝑦 𝑠 ) -𝑢 𝑀 (𝑥 𝑠 , 𝑦 𝑠 )) 2 𝑓 (𝑥 𝑠 , 𝑦 𝑠 ) 2 , (26)
where 𝑓(𝑥 𝑠 , 𝑦 𝑠 ) is assumed not to vanish and where L also includes points other than the sampling points. This is done in order to validate the algorithm, by evaluating the reference solution-which is a priori unknown in a general setting-at points different to the sampling ones,

E L = 1 #L ∑ s∈L √ (𝑓 (𝑥 𝑠 , 𝑦 𝑠 ) -𝑢 𝑀 (𝑥 𝑠 , 𝑦 𝑠 )) 2 𝑓 (𝑥 𝑠 , 𝑦 𝑠 ) 2 . ( 27
)
Since the s-PGD algorithm minimizes the error only at the sampling points P it is reasonable to expect that E P ≤ E L .

. . A Preliminary Example.

To test the convergence of the just presented algorithm, we consider

𝑓 1 (𝑥, 𝑦) = (cos (3𝜋𝑥) + sin (3𝜋𝑦)) 𝑦 2 + 4, (28)
which presents a quite oscillating behavior along the 𝑥 direction, whereas the 𝑦 direction is quadratic. We are interested in capturing such a function in the domain

Ω 𝑦 = Ω 𝑥 = [-1, 1].
Figures 3 and4 show the errors E P and E L in identifying the function 𝑓 1 (𝑥, 𝑦). In this case, we consider two distinct possibilities: no modal adaptivity at all and a modal adaptivity based on the residual, respectively. Several aspects can be highlighted. The first one is that E P (asterisks) decreases much faster when there is no modal adaptivity. This is expected, since we are minimizing with a richer basis since the very beginning, instead of starting with smooth functions like in the residual based approach. However, even if the minimization in the sampling points is well achieved, when no modal adaptivity is considered, the error out of the sampling points may increase as the solution is enriched with new modes. Nevertheless, the residual-based modal adaptivity alleviates this problem. As it can be noticed, starting with relatively smooth functions drives the solution out of the sampling points to be smooth as well, avoiding the problem of high oscillations appearing out of the sampling points.

A Local Approach to s-PGD

It is well-known that, as in POD, reduced basis or, in general, any other linear model reduction technique, PGD gives poor results-in the form of a nonparsimonious prediction-when the solution of the problem lives in a highly nonlinear manifold. Previous approaches to this difficulty included the employ of nonlinear dimensionality reduction techniques such as Locally Linear Embeddings [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF], kernel Principal Component Analysis [START_REF] Schölkopf | Nonlinear component analysis as a kernel eigenvalue problem[END_REF][START_REF] Scholkopf | Kernel principal component analysis[END_REF] or isomap techniques [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF]. Another, distinct, possibility, is to employ a local version of PGD [START_REF] Badas | Local proper generalized decomposition[END_REF], in which the domain is sliced so that at every subregion PGD provides optimal or nearly optimal results. We explore this last option here for the purpose of sparse regression, although a bit modified, as will be detailed hereafter.

The approach followed herein is based on the employ of the partition of unity property [START_REF] Melenk | The partition of unity finite element method: basic theory and applications[END_REF][START_REF] Babuška | The partition of unity method[END_REF]. In essence, it is well-known that any approximating function (like finite element shape functions, for instance) that forms a partition of unity can be enriched with an arbitrary function such that the resulting approximation inherits the smoothness of the partition of unity and the approximation properties of the enriching function.

With this philosophy in mind, we proposed to enrich a finite element mesh with an s-PGD approximation. The resulting approximation will be local, due to the compact support of finite element approximation, while inheriting the good approximation properties, already demonstrated, of s-PGD. In essence, what we propose is to construct an approximation of the type

𝑢 (𝑥, 𝑦) ≈ ∑ 𝑖∈I 𝑁 𝑖 (𝑥, 𝑦) 𝑢 𝑖 + ∑ 𝑝∈I enr ∑ 𝑒∈I 𝑝 𝑁 𝑒 (𝑥, 𝑦) 𝑀 ∑ 𝑘=1 𝑋 𝑘 𝑝 (𝑥) 𝑌 𝑘 𝑝 (𝑦) ⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟ 𝑓 enr 𝑝 (𝑥,𝑦) , (29)
where I represents the node set in the finite element mesh, I enr the set of enriched nodes, 𝑢 𝑖 are the nodal degrees of freedom of the mesh, I 𝑝 is the number of finite elements covered by node 𝑝 shape function's support and 𝑋 𝑘 𝑝 (𝑥) and 𝑌 𝑘 𝑝 (𝑥) functions are the 𝑘-th one-dimensional PGD modes enriching node 𝑝, that in fact constitute an enriching function 𝑓 enr (𝑥, 𝑦).

Of course, as already introduced in Eqs. ( 6) and ( 7), every PGD mode is in turn approximated by Galerkin projection on a judiciously chosen basis. In other words,

𝑢 (𝑥, 𝑦) ≈ ∑ 𝑖∈I 𝑁 𝑖 (𝑥, 𝑦) 𝑢 𝑖 + ∑ 𝑝∈I enr ∑ 𝑒∈I 𝑝 𝑁 𝑒 (𝑥, 𝑦) 𝑀 ∑ 𝑘=1 (N 𝑘 𝑥 ) 𝑇 a 𝑘 𝑝 (N 𝑘 𝑦 ) 𝑇 b 𝑘 𝑝 , (30)
with a 𝑘 𝑝 and b 𝑘 𝑝 the nodal values describing each onedimensional PGD mode.

In this framework, the definition of a suitable test function can be done in several ways. As a matter of fact, the test function can be expressed as the sum of a finite element and a PGD contribution,

𝑢 * = 𝑢 * 𝐹𝐸𝑀 + 𝑢 * 𝑃𝐺𝐷 , (31)
so that an approach similar to that of Eq. ( 4) can be accomplished.

An example of the performance of this approach is included in Section 4.4.

Numerical Results

The aim of this section is to compare the ability of sparse model identification for different interpolation techniques. On one hand, the performance of standard techniques based on Delaunay triangulation such as linear, nearest neighbor or cubic interpolation is compared. Even though these techniques are simple, they allow to have a nonstructured sampling point set since they rely on a Delaunay triangulation. On the other hand, the results are compared to the Sparse Subspace Learning (SSL) [START_REF] Borzacchiello | Non-intrusive Sparse Subspace Learning for Parametrized Problems[END_REF]. The convergence and robustness of this method is proven to be very effective since the points are collocated at the Gauss-Lobato-Chebychev points. However, two main drawbacks appear considering this method. The first one is that there is a high concentration of points in the boundary of the domain, so that this quadrature is meant for functions that vary mainly along the boundary. Indeed, if the variation of the function appears in the middle of the domain, many sampling points will be required to converge to the exact function. The second one is that the sampling points have to be located at specific points in the domain. The s-PGD method using simple Kriging interpolants will be compared as well.

The numerical results are structured as follows: first two synthetic 2D functions are analyzed; secondly, two 2D response surfaces coming from a thermal problem and a Plastic Yield function are reconstructed; finally, a 10D synthetic function is reconstructed by means of the s-PGD algorithm.

. . D Synthetic Functions. The first considered function is 𝑓 1 (𝑥, 𝑦), as introduced in the previous section. Figure 5 shows the reconstruction error (E L ) of 𝑓 1 (𝑥, 𝑦) for different sampling points. As it can be noticed, the s-PGD algorithm performs well for a wide range of sampling points. Nevertheless, the SSL method is the one presenting the lower error level when there are more than 150 sampling points.

A second synthetic function is defined as 𝑓 2 (𝑥, 𝑦) = cos (3𝑥𝑦) + log (𝑥 + 𝑦 + 2.05) + 5.

(

) 32
This function is intended to be reconstructed in the domain relatively smooth in the center of the domain, whereas the main variation is located along the boundary of the domain. Indeed, this function is meant to show the potential of the SSL technique.

Ω 𝑥 = Ω 𝑦 = [-
Figure 6 shows the reconstruction error of the 𝑓 2 (𝑥, 𝑦) function for different interpolation techniques. As it can be noticed, both SSL and s-PGD methods are the ones that present the best convergence properties. If the number of points is increased even more, the SSL method is the one that presents the lowest interpolation error. They are followed by linear and natural neighbor interpolations. Finally, the nearest neighbor method is the one presenting the worst error for this particular case.

. . D Response Surfaces Coming from Physical Problems.

Once the convergence of the methods have been unveiled for synthetic functions, it is very interesting to analyze the power of the former methods by trying to identify functions that are coming from either simulations or models popular in the computational mechanics community. Indeed, two functions will be analyzed: the first one is an anisotropic Plastic Yield function, whereas the second one is a solution coming from a quasistatic thermal problem with varying source term and conductivity.

Figure 7 shows the Yld2004-18p anisotropic plastic yield function, defined by Barlat et al. in [START_REF] Yoon | Prediction of six or eight ears in a drawn cup based on a new anisotropic yield function[END_REF]. Under plane stress hypothesis, this plastic yield function is a convex and closed surface defined in a three-dimensional space. Therefore, the position vector of an arbitrary point in the surface can be easily parameterized in cylindrical coordinates as 𝑅(𝜃, 𝜎 𝑥𝑦 ). The 𝑅(𝜃, 𝜎 𝑥𝑦 ) function for the Yld2004-18p is shown in Figure 8, where anisotropies can be easily seen. Otherwise, the radius function will be constant for a given 𝜎 𝑥𝑦 .

Figure 9 shows the error in the identification of Barlat's plastic yield function Yld2004-18p. As it can be noticed, the s-PGD technique outperforms the rest of techniques. Indeed, the s-PGD is exploiting the fact that the response surface is highly separable.

As mentioned above, the second problem is the sparse identification of the solution of a quasistatic thermal problem modeled by ∇ ⋅ (𝜂 (𝑥, 𝑡) ∇ (𝑢 (𝑥, 𝑡))) = 𝑓 (𝑡) ,

in Ω 𝑥 × Ω 𝑡 = [-1, 1] × [-1, 1] , (33)
where conductivity varies in space-time as 𝜂 (𝑥, 𝑡) = (1 + 10 abs (𝑥) + 10𝑥 2 ) log (𝑡 + 2.5)

𝑢 (1, 𝑡) = 2 (34) 𝑓 (𝑥, 𝑡) = 10 cos (3𝜋𝑡) 𝑢 (-1, 𝑡) = 2, (35)
and the source term varies in time. Homogeneous Dirichlet boundary conditions are imposed at both spatial boundaries and no initial conditions are required due to quasistationarity assumptions.

Figure 10 shows the evolution of the temperature field as a function of space time for the set of (33)-(35). It can be noticed how the variation of the temperature throughout time is caused mainly due to the source term. However, conductivity modifies locally the curvature of the temperature along the spatial axis. Symmetry with respect the 𝑥 = 0 axis is preserved due to the fact that the conductivity presents a symmetry along the same axis.

Figure 11 shows the performance of each one of the techniques when trying to reconstruct the temperature field from certain sampling points. As can be noticed, the s-PGD in conjunction with Kriging interpolants is the one that presents the fastest convergence rate to the actual function, which is considered unknown. It is followed by linear and natural interpolations. The SSL method presents a slow convergence rate in this case, due to the fact that the main variation of the function 𝑢(𝑥, 𝑡) is happening in the center of the domain and not in the boundary. . . A D Multivariate Case. In this subsection, we would like to show the scalability that s-PGD presents when dealing with relatively high-dimensional spaces. Since our solution is expressed in a separated format, an 𝑁 dimensional problem (𝑁𝐷) is solved as a sequence of 𝑁 1D problems, which are solved using a fixed-point algorithm in order to circumvent the nonlinearity of the separation of variables.

The objective function that we have used to analyze the properties of the s-PGD is defined as

𝑓 3 (𝑥 1 , 𝑥 2 , . . . , 𝑥 𝑁 ) = 2 + 1 8 𝑁 ∑ 𝑖=1 𝑥 𝑖 + 𝑁 ∏ 𝑖=1 𝑥 𝑖 + 𝑁 ∏ 𝑖=1 𝑥 2 𝑖 , (36)
with 𝑁 = 10 in this case.

Figure 12 shows the error convergence in both sampling points (E P , asterisks) and points out of the sampling (E L , filled points). The L data set was composed by 3000 points, and the P data subset for the s-PGD algorithm was composed by 500 points. The number of points required to properly capture the hypersurface has increased with respect to the 2D examples due to the high dimensionality of the problem. Special attention has to be paid when increasing the cardinality of the interpolant basis without many sampling points, because the problem of high oscillations outside the control points may be accentuated.

.

. An Example of the Performance of the Local s-PGD.

The last example corresponds to the sparse regression of an intricate surface that has been created by mixing three different Gaussian surfaces so as to generate a surface with no easy separate representation (a nonparsimonious model, if we employ a different vocabulary). The appearance of this synthetic surface is shown in Figure 13.

The sought surface is defined in the domain Ω = [0, 1] 2 , which has been split into the finite element mesh shown in Figure 14. Every element in the mesh has been colored according to the number of enriching PGD functions, ranging from a single one to four. The convergence plot of this example as a function of the number of PGD modes added to the approximating space is included in Figure 15.

Conclusions

In this paper we have developed a data-based sparse reducedorder regression technique under the Proper Generalized Decomposition framework. This algorithm combines the robustness typical of the separation of variables together with properties of collocation methods in order to provide with parsimonious models for the data at hand. The performance of simple Kriging interpolation has proven to be effective when the sought model presents some regularity. Furthermore, a modal adaptivity technique has been proposed in order to avoid high oscillations out of the sampling points, characteristic of high-order interpolation methods when data is sparse. For problems in which the result lives in a highly nonlinear manifold, a local version of the technique, which makes use of the partition of unity property, has also been developed. This local version outperforms the standard one for very intricate responses.

The s-PGD method has been compared advantageously versus other existing methods for different example functions. Finally, the convergence of s-PGD method for a high dimensional function has been demonstrated as well.

Although the sparsity of the obtained solution could not seem evident for the reader, we must highlight the fact that the very nature of the PGD strategy a priori selects those terms on the basis that it plays a relevant role in the approximation. So to speak, PGD algorithms automatically discard those terms that in other circumstances will be weighted by zero. Sparsity, in this sense, is equivalent in this context to the number of sums in the PGD separated approximation. If only a few terms are enough to reconstruct the data-as is almost always the case-then sparsity is guaranteed in practice.

Sampling strategies other than the Latin hypercube method could be examined as well. This and the coupling with error indicators to establish good stopping criteria constitute our effort of research at this moment. In fact, the use of reliable error estimators could even allow for the obtention of adaptive samplings in which the cardinality of the basis could be different along different directions.

1 XFigure 1 :XFigure 2 :

 112 Figure 1: Kriging shape functions using 𝜎 = √ (𝑥 𝑖+1 -𝑥 𝑖 ) 2 for 7 control points uniformly distributed along the 1D domain.

Figure 3 :

 3 Figure 3: E L (points) and E P (asterisk) versus the number of modes for 𝑓 1 (𝑥, 𝑦), #P = 100, #L = 1000. No modal adaptivity.

Figure 4 :

 4 Figure 4: E L (points) and E P (asterisk) versus the number of modes for 𝑓 1 (𝑥, 𝑦), #P = 100, #L = 1000. Modal adaptivity based on the residual, 𝜖 𝑟 = 1e-2.

Figure 5 :

 5 Figure 5: E L of 𝑓 1 (𝑥, 𝑦) varying #P for different identification techniques. #L = 1000.

Figure 6 :

 6 Figure 6: E L of 𝑓 2 (𝑥, 𝑦) varying #P for different identification techniques. #L = 1000.

Figure 7 :

 7 Figure 7: Barlat's Yld2004-18p function under plane stress hypothesis.

Figure 8 :

 8 Figure 8: 𝑅(𝜃, 𝜎 𝑥𝑦 ) function for Barlat's Yld2004-18p yield function.

Figure 9 :Figure 10 :

 910 Figure 9: E L of 𝑅(𝜃, 𝜎 𝑥𝑦 ) varying #P for different sparse identification techniques. #L = 1000. 𝜖 𝑟 = 5 ⋅ 10 -4 .

Figure 11 :

 11 Figure 11: E L of 𝑢(𝑥, 𝑡) varying #P for different identification techniques. #L = 1000. 𝜖 𝑟 = 2.5 ⋅ 10 -3 .

Figure 12 :

 12 Figure 12: E L (points) and E P (asterisk) versus the number of modes for 𝑓 3 (𝑥 1 , 𝑥 2 , . . . , 𝑥 𝑁 ), #P = 500, #L = 3000. Modal adaptivity based on the residual, 𝜖 𝑟 = 1𝑒 -3.

Figure 13 :Figure 14 :Figure 15 :

 131415 Figure 13: A synthetic surface generated by superposition of three different Gaussians, that is to be approximated by local s-PGD techniques.

Acknowledgments

This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie Grant agreement no. 675919, also the Spanish Ministry of Economy and Competitiveness through Grants nos. DPI2017-85139-C2-1-R and DPI2015-72365-EXP, and the Regional Government of Aragon and the European Social Fund, research group T24 17R.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


Ph.D Cyril Atkinson-Clement
email: c.atkinson-clement@icm-institute.org

PhD Alban Letanneux

MD Guillaume Baille

MSc Marie-Charlotte Cuartero

MSc Lauriane Véron-Delor

PhD Camille Robieux

MSc Manon Berthelot

MD, PhD Danièle Robert

MD, PhD Jean-Philippe Azulay

MD, PhD Luc Defebvre

MD, PhD Joaquim Ferreira

MD, PhD Alexandre Eusebio

MD, PhD Caroline Moreau

PhD Serge Pinto

L A Dysarthrie Parmi

Psychosocial Impact of dysarthria: the patient-reported outcome as part of the clinical management

Keywords: Dysarthria, Psychosocial impact, Voice, Speech, Self-assessment 1

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

Introduction

Dysarthria is common in neurological movement disorders and represents 53% of acquired communication impairments [1]. Six kinds of dysarthria were initially described by Darley et al. [2] (Fig. 1). Hypokinetic and hyperkinetic dysarthrias result from subcortical dysfunctions respectively associated to reduced and increased movements; ataxic dysarthria is associated with cerebellar dysfunctions. Bilateral upper motoneuron alterations can induce spastic dysarthria, while the dysfunction of the final motor pathway can lead to flaccid dysarthria. A mixed dysarthria results from multiple brain alterations. Duffy [1] introduced two additional dysarthria types resulting from a unilateral upper motor deficit and from undetermined etiology. Dysarthria is generally assessed from a pathophysiological perspective, in that a specific disease is viewed to be the cause of a particular speech impairment. During a routine consultation, neurologists frequently rely on rating speech intelligibility, and if required, further clinical assessments and acoustic analyses can provide additional information, such as objective measures on laryngeal and supra-laryngeal articulations.

Fig.1. The different types of dysarthria according to the location of pathophysiological dysfunctions (based on previous dysarthrias classification).

Hypokinetic and hyperkinetic dysarthrias are associated to basal ganglia dysfunctions, while ataxic dysarthria results from cerebellar dysfunction. Unilateral Upper Motor Neuron (UUMN), spastic and flaccid dysarthrias are the sequel of pyramidal tract dysfunctions (upper motor neuron and unilateral for UUMN, upper motor neuron and bilateral for spastic dysarthria, and lower motor neuron for flaccid dysarthria). Not mentioned in this representation: mixed dysarthrias (resulting from multiple subcortical/cortical lesions) and dysarthria of undetermined etiology. Dysarthria leads to communication deficits, which frequently induce social isolation [3], a decrease in quality of life [4], and alterations in daily-living activities [5]. Such deficits also increase the risk of developing psycho-emotional disorders such as depression [START_REF] Miller | How do I sound to me? Perceived changes in communication in Parkinson's disease[END_REF], or even cognitive deterioration [START_REF] Rektorova | Speech prosody impairment predicts cognitive decline in Parkinson's disease[END_REF]. Yet, the majority of assessments dedicated to speech/voice impairments are based on the symptoms themselves [START_REF] Jacobson | The Voice Handicap Index (VHI): Development and Validation[END_REF][START_REF] Kasper | Voice handicap index and voice-related quality of life in small laryngeal carcinoma[END_REF], and not on the symptoms' consequences. A "biomedical model of illness is inadequate in understanding the full impact of communication disorders" and "psychosocial factors should be considered" [START_REF] Wray | Longer-term needs of stroke survivors with communication difficulties living in the community: a systematic review and thematic synthesis of qualitative studies[END_REF]. This practice leaves unanswered the symptoms' effects on quality of life, despite quality of life being dependent on those consequences that affect physical, psychological and social components [START_REF]WHO | International Classification of Functioning, Disability and Health (ICF[END_REF]. A method of bridging the gap between symptom consequences on quality of life is to monitor patient-reported outcomes (PROs). PROs appear as particularly important contributions to clinical practice, especially because they report functional facts that are not systematically in agreement with the degree of severity of the symptom [START_REF] Wiklund | The relationship between exercise tolerance and quality of life in angina pectoris[END_REF]. As stated by the PRO Harmonization Group (www.eriqa-project.com), PROs allow clinicians: 1) to better acquaint themselves with symptom consequences; 2) to assess the relevance of a treatment; 3) to guide medical choices; 4) and to improve the interpretation of clinical outcomes [START_REF] Acquadro | Incorporating the patient's perspective into drug development and communication: an ad hoc task force report of the Patient-Reported Outcomes (PRO) Harmonization Group meeting at the Food and Drug Administration, February 16, 2001[END_REF].

Among the assessments of PROs is the Dysarthria Impact Profile (DIP), which was developed and validated in English, with patients presenting with different kinds of dysarthria [START_REF] Walshe | Dysarthria impact profile: development of a scale to measure psychosocial effects[END_REF]. Its adaptations in French [START_REF] Letanneux | The Dysarthria Impact Profile: A Preliminary French Experience with Parkinson's Disease[END_REF] and European Portuguese [START_REF] Cardoso | Psychosocial impact of Parkinson's disease-associated dysarthria: Cross-cultural adaptation and validation of the Dysarthria Impact Profile into European Portuguese: EP: DIP validation for PD[END_REF] have been performed, but only for hypokinetic dysarthria in PD patients. We propose here a revised French version for which all previously acknowledged limitations [START_REF] Letanneux | The Dysarthria Impact Profile: A Preliminary French Experience with Parkinson's Disease[END_REF] were overcome. The objectives of the present study were to confirm the psychometric properties of this revised version of the French DIP with patients presenting with different kinds of dysarthria, and to further demonstrate the relevance of using psychosocial indicators as part of the global assessment and management of dysarthria.

Materials and Methods

Participants

We recruited 105 patients pertaining to seven pathological groups ( For intelligibility, DIP and SF-36: the higher the score, the lower the degree of perceived handicap. For VHI: the higher the score, the higher the degree of perceived handicap.

Acoustic measures & clinical assessments

Speech intelligibility was assessed by a speech pathologist on the basis of speech records provided by the intelligibility section of the Frenchay Dysarthria Assessment (FDA-2). The participants read ten words and ten sentences, randomly and blindly chosen from a list of 116 words and 50 sentences. Intelligibility scores were reported as percentages of correctly transcribed items [START_REF] Enderby | FDA-2 : Frenchay Dysarthria Assessment : Examiner's Manual[END_REF].

Among acoustic parameters, means and standard deviations of fundamental frequency (F0) and speech loudness (SPL, sound pressure level) were extracted from the sentences read by the patients. These measures were focused on in our study since monopitch and monoloudness are particular deficits often displayed in dysarthric speech.

The Montreal Cognitive Assessment (MoCA [START_REF] Nasreddine | The Montreal Cognitive Assessment, MoCA: a brief screening tool for mild cognitive impairment[END_REF]) was performed to detect any cognitive impairment.

The 36-Item Short Form Health Survey (SF-36 [START_REF] Ware | The MOS 36-Item Short-Form Health Survey (SF-36): I. Conceptual Framework and Item Selection[END_REF]) is a non-specific quality of life self-questionnaire, which examines four physical and four mental subdomains. Each subdomain is scored from 0 (maximal disturbance) to 100 (normal).

The Voice Handicap Index (VHI [START_REF] Jacobson | The Voice Handicap Index (VHI): Development and Validation[END_REF]) is often considered as the gold standard for voice self-perception [START_REF] Kasper | Voice handicap index and voice-related quality of life in small laryngeal carcinoma[END_REF]. It includes 30 items split into three domains: physical, functional, emotional. Each item is scored from 0 to 4: the higher the score, the higher the degree of perceived handicap.

The DIP [START_REF] Walshe | Dysarthria impact profile: development of a scale to measure psychosocial effects[END_REF][START_REF] Letanneux | The Dysarthria Impact Profile: A Preliminary French Experience with Parkinson's Disease[END_REF][START_REF] Cardoso | Psychosocial impact of Parkinson's disease-associated dysarthria: Cross-cultural adaptation and validation of the Dysarthria Impact Profile into European Portuguese: EP: DIP validation for PD[END_REF] is composed of four numerical subscales (A: 'the effect of dysarthria on me as a speaker'; B: 'accepting my dysarthria'; C: 'how I feel others react to my speech'; D: 'how dysarthria affects my communication with others'). Each item is scored from 1 (major impact) to 5 (no impact): the lower the score, the higher the degree of psychosocial handicap. The last subscale (part E) was not included in this analysis since it is a qualitative section that does not provide any numerical score.

Statistical analyses

All statistical analyses were performed with R software. We calculated Cronbach's α to assess the internal consistency of the DIP (an adequate consistency would correspond to a coefficient of at least .70). We performed Pearson's correlations (after Holm's correction) between the scales in order to estimate the concurrent validity of the DIP, and Receiver Operating Characteristic (ROC) analyses between: 1) HC and all patients to determine discriminant validity (sensitivity, specificity, accuracy [area under the curve], power and cut-off) of the DIP, and 2) for each items of the DIP, between HC and patients to determine discriminant items and between each group of patients and the other patients to determine any specificity of each disease. Also, Hedges' g effect sizes were calculated to compare HC to patients. Considering the number of analyses performed, and the rather small number of participants by group, we considered as meaningful only large effect sizes (g≥.8). Finally, we performed one linear regression to determine the influences of both intelligibility and the pathological group on the DIP total score.

Results

Psychometric properties of the DIP

For the total DIP score, internal consistency was confirmed for all groups (α>.7; Table.2.A); it was also confirmed for the four subscales when we considered all participants together: A (α=.81), B (α=.71), C (α=0.79) and D (α=.85). The correlations we generated demonstrated that the DIP had good construct validity (Table .2.B): the total DIP score was highly correlated to the total VHI (r=-.77), SF-36 (r=.53) and speech intelligibility (r=.43). No correlation was found between the DIP and any acoustic measure. The highest correlations were found between the emotional part of the VHI and the subscales A, B and C of the DIP (respectively, r=-.65, r=-.62, r=-.63), and between the functional part of the VHI and the subscale D of the DIP (r=-.78). Discriminant validity was confirmed for the total DIP score (accuracy=.95; cut-off=180/240) and for each of the subscales: A (accuracy=.91; cut-off=44/60), B (accuracy=.87; cut-off=38/50), C (accuracy=.87; cut-off=48/70) and D (accuracy=.96; cut-off=48/60; Table .2.C). The cutoff of the total DIP score allowed us to determine the proportion of patients who had a communication complaint: 66.6% of PD patients, 73.3% of ataxic patients, 80% of dystonic patients, 86.6% of HD patients, 93.3% of MSA patients; 100% of PSP patients and 100% of ALS patients (Fig. 2).

Fig.2. DIP score distributions by pathological group according to the calculated cut-off.

The vertical line represents the cut-off (180/240) for the DIP total score obtained with the ROC analysis: below 180, the DIP reveals statistically a psychosocial impact of dysarthria; above 180, the DIP score does not report an impact with functional meaning. ALS: Amyotrophic Lateral Sclerosis; DIP: Dysarthria Impact Profile; HC: Healthy control; HD: Huntington's disease; MSA: Multiple System Atrophy; PD: Parkinson's disease; PSP: Progressive Supranuclear Palsy.

Patient phenotypes according to the communication complaint

Patient phenotypes were ranked according to the total DIP score, from the smallest to the highest complaint: PD (x=169.1), Ataxia (𝑥=155.2), HD (𝑥=152.5), Dystonia (x=147.9), PSP (𝑥=142.5), MSA (x=140.3) and ALS (x=137.7). Hedges' g effect sizes between patients and HC for the DIP and VHI subscales demonstrated that all patients reported a dysarthria complaint (g>.8). However, PD, ataxic and PSP patients did not show any decrease of intelligibility (g<.8). Most of the items were discriminant between each group of patients and the controls (Fig. 3, upper part of the right panels). No items reached the threshold of accuracy between PD patients and other patients, while only one item was significant for ataxic patients (C13; p=.008) and dystonic patients (D04; p=.027). For HD, PSP and MSA patients, between 4 and 5 items were significantly discriminant from other patients (Fig. 3). Lastly, for ALS patients, eleven items, ranging from A06 to D12 were significant (p<.05). We also observed that both the pathological group (F(7;104)=12.97; p<.001) and the intelligibility score (F(1;104)=6.87; p=.01) had an effect on the total DIP score, but without any significant interaction (F(7;104)=0.97; p=.45). Left panel: the profile diagrams are based on Hedges' g effect sizes between healthy controls and each group of patients. Right panel: accuracy scores for each DIP items (from A01 to D12) between the pathological group and the healthy control group (top panel), and between the pathological group and the other pathological groups (bottom panel). The colored bars correspond to meaningful items. ALS: Amyotrophic Lateral Sclerosis; DIP: Dysarthria Impact Profile (A: "the effect of dysarthria on me as a speaker"; B: "accepting my dysarthria"; C: "how I feel others react to my speech"; D: "how dysarthria affects my communication with others"); HD: Huntington's disease; MSA: Multiple System Atrophy; ns: not significant; PD: Parkinson's disease; PSP: Progressive Supranuclear Palsy; VHI: Voice Handicap Index (F: functional; P: physical; E: emotional).

Discussion

All patients reported a communication complaint, regardless of the presence of speech impairment clinically estimated via speech intelligibility. Interestingly, PD, ataxic and PSP patients did not show any deficit of speech intelligibility. Our study suggests that it would be relevant to take into account psychosocial indicators of dysarthria in patients with neurodegenerative disorders.

Psychometric confirmation of the DIP

Our study confirmed the psychometric properties of the revised French version of the DIP, thus validating the use of this self-questionnaire for estimating the impact of dysarthria in patients with several neurodegenerative disorders. The original DIP was administered to 31 people with different types of dysarthria [START_REF] Walshe | Dysarthria impact profile: development of a scale to measure psychosocial effects[END_REF], and consecutive adaptations in other languages confirmed its relevance for the management of PD [START_REF] Letanneux | The Dysarthria Impact Profile: A Preliminary French Experience with Parkinson's Disease[END_REF][START_REF] Cardoso | Psychosocial impact of Parkinson's disease-associated dysarthria: Cross-cultural adaptation and validation of the Dysarthria Impact Profile into European Portuguese: EP: DIP validation for PD[END_REF]. Our results go one step beyond the use of the DIP as part of the clinical assessment of dysarthria, providing such evidence from a protocol that overcame the limitations of the previous version. Notably, we included a larger number of participants with different diseases and different kinds of dysarthria, we rephrased unclear statements and we perfected the formal presentation of the questionnaire (Supplementary material). In addition, the new formal organization of the DIP makes it easy to quickly determine whether the questionnaire has been properly completed: answering the questionnaire requires alternating between the responses "agree" and "disagree". As suggested by the absence of correlation between the DIP and the MoCA, the DIP score is either not or only weakly driven by the cognitive status of the patients. The DIP can thus be administered to patients with mild cognitive impairment, but we recommend that patients with higher cognitive impairments fill out the DIP with the help of a clinician or a caregiver.

Speech intelligibility is not functional communication

Our findings suggest that intelligibility loss is crucially correlated to the psychosocial impact of dysarthria. This means that, for a similar level of intelligibility impairment, the DIP total score is similar regardless of the pathological group. Nevertheless, the intelligibility score remains insufficient to address the psychosocial impact of dysarthria for at least two reasons: 1) the intelligibility score explained only 18.5% of the total DIP score (determination coefficient; r 2 ); 2) two patient groups (PD and Ataxia) displayed a significant DIP score (in comparison to the controls), but a level of speech impairment not so significant to impair intelligibility. Our findings suggest that the psychosocial impact measured by the DIP could be partially independent from the severity of dysarthria (indirectly addressed here via speech intelligibility): the DIP was able to detect patients without any intelligibility impairment, but with a psychosocial impact. This might be related to the observed independence between the perspective of the patients and those of the caregivers on communication effectiveness [START_REF] Mcauliffe | Speech intelligibility and perceptions of communication effectiveness by speakers with dysarthria following traumatic brain injury and their communication partners[END_REF]. Thus, the assessment of speech intelligibility, which is central in the neurological investigation of dysarthria, remains insufficient to address the patient's complaints and could not represent per se the sole indicator of speech evaluation and management in clinical routine [START_REF] Barnish | Speech and communication in Parkinson's disease: a cross-sectional exploratory study in the UK[END_REF]. Indeed, one study showed that the majority of patients with Parkinsonian syndromes reported a communication impairment, even if some of them remained intelligible [START_REF] Wiblin | The Importance of Connection to Others in QoL in MSA and PSP[END_REF]. It is necessary to consider that communication ability depends on intelligibility, but also on how patients feel about the perception of others [START_REF] Wiblin | The Importance of Connection to Others in QoL in MSA and PSP[END_REF]. For example, based on a literature review, only 5% of scientific research in 2007 focused on psychosocial aspects of dysarthria [START_REF] Yorkston | The degenerative dysarthrias: a window into critical clinical and research issues[END_REF]. The ROC analysis of DIP items demonstrated (especially for PD, ataxic and PSP patients who had no significant intelligibility deficit) that psychosocial consequences of dysarthria could be due to the reactions of others ("I never feel that others laugh or make fun of my speech" [C13, significant for ataxic patients]). The analysis also showed that psychosocial consequences could lower patients' self-perception ("I am as confident now as I was before I had a speech problem" [A05, significant for PSP patients]) and contribute to the isolation of patients ("Although I have difficulty speaking I do not avoid communicating with people I know" [D01, significant for PSP patients]; "My social life has not changed as a result of my dysarthria" [D02, significant for PSP patients]). Thus, our results are in agreement with a previous qualitative report that focused on acquired dysarthria and for which several themes were frequently mentioned by patients: "communication has changed", "people treat me differently", "dysarthria resulting in negative emotions", "barriers to communicate", "life is different now" [START_REF] Walshe | Living with acquired dysarthria: the speaker's perspective[END_REF]. So far, this statement remains topical and underlies the fact that dysarthria is a functional symptom which is still undervalued.

Enhancing clinical assessment and management of dysarthria with the PRO

On the exclusive basis of the neurological assessment, patients with dysarthria are generally referred to a speech therapist/pathologist only when intelligibility is significantly impaired. This might explain the gap between the number of patients with dysarthria complaints and the number of patients who have received speech therapy. For example, it has been reported that 70 to 90% of PD patients complain about speech impairment; 29% of the patients reported that dysarthria had a significant impact in their life; but only 3% accessed speech therapy [START_REF] Hartelius | Speech and swallowing symptoms associated with Parkinson's disease and multiple sclerosis: a survey[END_REF]; though these numbers have recently improved [START_REF] Schalling | Speech and Communication Changes Reported by People with Parkinson's Disease[END_REF]. Speech therapy can improve communication abilities and reduce social isolation [START_REF] Fogg-Rogers | Choral singing therapy following stroke or Parkinson's disease: an exploration of participants' experiences[END_REF]. Thus, patients should be seen earlier than they are currently to help reduce the psychosocial impact. Indeed, a recent study on PD reported that a longer disease duration is associated with a lower speech therapy efficacy [START_REF] Boutsen | Prosodic Improvement in Persons with Parkinson Disease Receiving SPEAK OUT!® Voice Therapy[END_REF]. The ROC analyses we conducted provided a cut-off of the DIP, below which the complaint of the patients should be further discussed with them and possibly, the patients could be referred to a speech/language therapist for additional counselling: clinicians could use the DIP as a tool to precisely estimate the psychosocial consequences of dysarthria [START_REF] Pinto | Dysarthria in individuals with Parkinson's disease: a protocol for a binational, cross-sectional, case-controlled study in French and European Portuguese (FraLusoPark)[END_REF]. Furthermore, to better manage dysarthria, it is crucial to include patient-reported difficulties in functional communication [START_REF] Schalling | Speech and Communication Changes Reported by People with Parkinson's Disease[END_REF] as indicators of both dysarthria's psychosocial impact [START_REF] Pinto | Dysarthria in individuals with Parkinson's disease: a protocol for a binational, cross-sectional, case-controlled study in French and European Portuguese (FraLusoPark)[END_REF] and speech therapy management [START_REF] Fogg-Rogers | Choral singing therapy following stroke or Parkinson's disease: an exploration of participants' experiences[END_REF][START_REF] Ramig | Speech treatment in Parkinson's disease: Randomized controlled trial (RCT): SPEECH TREATMENT IN PARKINSON'S DISEASE: RCT[END_REF][START_REF] Mackenzie | Behavioural intervention effects in dysarthria following stroke: communication effectiveness, intelligibility and dysarthria impact[END_REF]. Thus, using the DIP as a psychosocial indicator can be useful to correctly adjusting a speech rehabilitation program [START_REF] Atkinson-Clement | Behavioral treatments for speech in Parkinson's disease: meta-analyses and review of the literature[END_REF][START_REF] Enderby | Disorders of communication[END_REF].

Limitations of the study

Our study included some limitations. First, our patient groups included an amount of intra-group heterogeneity in terms of symptom expression, etiology, disease severity, disease duration, etc. However, our objective was to highlight the relevance of using the DIP with patients suffering from different neurodegenerative diseases, with no specific criteria of exclusion (except other neurological or psychiatric history). Second, we did not have access to information related to social environment, professional context and even cultural information. These data could be of importance to better identify complaints of communication, and have to be considered in clinical practice on a patient-bypatient individual basis. Finally, we did not evaluate depression, apathy or anhedonia, which could contribute to altering functional communication. Additional and specific studies are still required to overcome these limitations.

Conclusion

Dysarthria is a frequent and debilitating symptom of neurological disorders, and is from the patient perspective, a multidimensional symptom, impacting also social and emotional daily-living components [START_REF] Walshe | Living with acquired dysarthria: the speaker's perspective[END_REF]. However, the proportion of patients who access speech therapy is significantly lower than the proportion of patients who complain about communication. Integrating an additional tool such as the DIP, which provides the patient's point-of-view, is important to better identify patients with a communication complaint, especially when these patients remain intelligible. Thus, communication ought to be assessed from a holistic perspective, including both professional expertise and the patient's feelings.

Author Contributions

Fig. 3 .

 3 Fig.3. Dysarthria phenotypes according to the pathological groups.Left panel: the profile diagrams are based on Hedges' g effect sizes between healthy controls and each group of patients. Right panel: accuracy scores for each DIP items (from A01 to D12) between the pathological group and the healthy control group (top panel), and between the pathological group and the other pathological groups (bottom panel). The colored bars correspond to meaningful items. ALS: Amyotrophic Lateral Sclerosis; DIP: Dysarthria Impact Profile (A: "the effect of dysarthria on me as a speaker"; B: "accepting my dysarthria"; C: "how I feel others react to my speech"; D: "how dysarthria affects my communication with others"); HD: Huntington's disease; MSA: Multiple System Atrophy; ns: not significant; PD: Parkinson's disease; PSP: Progressive Supranuclear Palsy; VHI: Voice Handicap Index (F: functional; P: physical; E: emotional).

Table 1 : Demographic data and clinical outcomes of participants.

 1

	15 per group) and associated with

ALS: Amyotrophic Lateral Sclerosis; DIP: Dysarthria Impact Profile; F: Female; HC: Healthy control; HD: Huntington's disease; M: Male; MoCA: Montreal Cognitive Assessment; MSA: Multiple System Atrophy; N: number of participants; PD: Parkinson's disease; PSP: Progressive Supranuclear Palsy; SD: Standard Deviation; SF-36: 36-Item Short Form Health Survey; VHI: Voice Handicap Index.

Table 2 : Psychometric properties of the DIP: A) internal consistency (Cronbach's α), B) construct validity (Pearson's correlations) and C) discriminant validity (ROC).

 2

	DIP

  1. Research Project: A. Conception, B. Organization, C. Execution; 2. Statistical Analysis: A. Design, B. Execution, C. Review and Critique; 3. Manuscript Preparation: A. Writing the First Draft, B. Review and Critique.

	CAC: 1A;1B;1C;2A;2B;3A
	AL: 1A;1B;1C;3B
	GB: 1C;2C;3B
	MCC: 1C;2C;3B
	LVD: 1C;2C;3B
	CR: 1C;3B
	MB: 1A;1C;3B
	DR: 1C;3B
	JPA: 1C;3B
	LD: 1C;3B
	JF: 3B
	AE: 1C;3B
	CM: 1C;3B
	SP: 1A;1B;2C;3A

Funding Sources

This work was supported by the PACA Regional Council and Orthomalin [Ph.D grant]; and by a France/Portugal bilateral funding [ANR-13-ISH2-0001-01; FCT-ANR/NEU-SCC/0005/2013].

Statements

Acknowledgement

The authors would like to thank the neurologists, the ENTs and the clinical research associates who contributed to the patient recruitment: in Lille (Nicolas Carriere, Adrian Degardin, David Devos, Alexandre Kreisler, Clémence Simonin) and in Marseille (Frédérique Fluchère, Laura Mundler, Manel Nouira, Elsa Robin, Sophie Soulayrol, Tatiana Witjas). The authors would also like to thank Marie-Hélène Mignard, Oriana Collins and Karl Neergaard for helpful English revisions of the article. Also, the authors wish to thank Margaret Walshe, Richard K. Peach and Nick Miller for their pioneering work on the Dysarthria Impact Profile.

Statement of Ethics

Subjects have given their informed consent.

The study protocol was approved by the local research institute's committee on human research.

Disclosure Statement

The authors have no conflicts of interest to declare.