Pierre Courrieu
email: courrieu@free.fr

On the convergence of rank-one multi-target linear regression

Keywords: multi-target linear regression, multi-output linear regression, reduced rank multivariate regression, regression coefficient estimator consistency

This paper presents a new method for solving rank-one multivariate regression problems, providing a solution that maximizes the sum of squared correlations of the one-dimensional fitted pattern with the target variates. The suitability of the method and the consistency of the estimator are formally proved and experimentally tested. In particular, it is shown that the estimate converges not only as a function of the number of items, but also as a function of the number of target variates.

An equivalent conventional reduced-rank regression case is identified, and it inherits the convergence properties of the new approach. Application programs in Matlab/Octave code are provided and numerical examples using artificial data as well as real data are presented.

Introduction

In this paper, we consider the following maximization problem. Given a random sample of m items, a set of n explicative variables (typically item attributes), and a set of d target variates (typically observation data), let Ξ ∈ R m×n be the matrix of regressors, and let Θ ∈ R m×d be the target matrix. Consider the matrix X ∈ R m×n whose columns are those of Ξ centred on their mean, and the matrix T ∈ R m×d whose columns are those of Θ centred on their mean and normalized to 1. Then find w ∈ R n such that: w = arg maxv∈R n Σk=1..d r 2 (Xv, Tk), [START_REF] Anderson | Estimating linear restrictions on regression coefficients for multivariate normal distributions[END_REF] where Tk is the kth column of T, and r(x, y) is the Pearson correlation coefficient between vectors x and y. Note that one can as well use Ξ and Θ in (1) instead of X and T, without changing the result since squared correlations are invariant to linear transforms of the variables. This is a multiple linear regression problem if n>1, and it is a multitarget linear regression problem if d>1. Note that in the multi-target case, only one linear combination of regressors (Xw) is fitted to all the target variates, which is called a "rank-one multi-target linear regression", and in the present case is also "onedimensional" since the solution is a simple vector. This is different from the common practice where one fits regressors to each target variate independently of other ones.

Thus, the question here is not simply to know in what measure the considered regressors can account for each target variate, but rather to know what the target variates have in common that the regressors can account for. This type of approach is known as "multi-target linear regression", "multi-output linear regression", or "multivariate linear regression", and the "rank-one" qualifier refers to the fact that this corresponds to a particular case of the more general "reduced-rank regression"

approach [START_REF] Anderson | Estimating linear restrictions on regression coefficients for multivariate normal distributions[END_REF][START_REF] Izenman | Reduced-rank regression for the multivariate linear model[END_REF][START_REF] Izenman | Modern Multivariate Statistical Techniques[END_REF], where one estimates matrices of regression coefficients of various ranks.

The rank-one, one-dimensional case is of special interest in a number of situations. For instance, one can need to estimate a unique vector of regression coefficients allowing combining a given set of predictive variables in order to predict new target variates and possibly using new items. Such a generalization process is common in machine-learning, however, it can also be useful to build reusable regressors in the context of data analysis. Another important application field in neurosciences is the regression analysis of electrophysiological data such as electroencephalograms (EEG), in particular event related potentials (ERP), where one uses stimulus and/or participant characteristics to explain a sequence of electrophysiological responses at various cortical (scalp) locations [START_REF] Hauk | The time course of visual word recognition as revealed by linear regression analysis of ERP data[END_REF][START_REF] Rey | Accounting for variance in single-word ERPs[END_REF]. Using a classical multiple regression technique, one can observe significant fits at various locations and time delays in the EEG, while the electrophysiological response patterns are actually very different at these various locations and latencies. Using a rank-one, one-dimensional multi-target regression avoids such a situation by fitting only what is common in the various target response patterns, possibly with different strengths and signs.

Non-linear multi-target regression methods are common in machine-learning studies [START_REF] Borchani | A survey on multi-output regression[END_REF], while linear multi-target regression methods have been more particularly investigated in the framework of multivariate data analysis [START_REF] Anderson | Estimating linear restrictions on regression coefficients for multivariate normal distributions[END_REF][START_REF] Breiman | Predicting multivariate responses in multiple linear regression[END_REF][START_REF] Brillinger | Time Series: Data Analysis and Theory[END_REF][START_REF] Izenman | Reduced-rank regression for the multivariate linear model[END_REF][START_REF] Izenman | Modern Multivariate Statistical Techniques[END_REF][START_REF] Reinsel | Multivariate Reduced-Rank Regression: Theory and Applications[END_REF][START_REF] Similä | Input selection and shrinkage in multiresponse linear regression[END_REF]. As mentioned above, the problem (1) closely relates to the rank-one case of the well-known reduced-rank regression (commonly abbreviated "RRR") approach [START_REF] Izenman | Reduced-rank regression for the multivariate linear model[END_REF][START_REF] Izenman | Modern Multivariate Statistical Techniques[END_REF]. Several variants of the reduced-rank regression have been proposed [START_REF] Reinsel | Multivariate Reduced-Rank Regression: Theory and Applications[END_REF], and recent studies often focus on the determination of an optimal rank for the matrix of regression coefficients [START_REF] Buena | Optimal selection of reduced rank estimators of high-dimensional matrices[END_REF][START_REF] Giraud | Low rank multivariate regression[END_REF][START_REF] Kargin | On estimation in the reduced-rank regression with a large number of responses and predictors[END_REF].

In the present study, however, we a priori fix this rank to one, and we concentrate on certain convergence properties that are relevant for the targeted applications, in particular the convergence of the regression coefficients estimate as a function of the number of target variates. The convergence as a function of the number of items is known [26, section 2.5], however the convergence as a function of the number of targets, independently of the number of items, has not been studied (to our knowledge), and it is commonly assumed that n+d ≤ m, see [26, p. 3]. We will not retain this assumption, letting d grow freely, and we will define a particular rank-one multitarget regression method that is quite simple and allows us to easily derive the desired convergence properties. We also provide indications for the implementation of this tool, and a ready to use Matlab/Octave code program for applications. Finally, the method and the convergence properties are tested on numerical examples.

T' denotes the transpose of T, and X (1,2,3) denotes any {1,2,3}-inverse of X, such as the unique {1,2,3,4}-inverse of Moore-Penrose [START_REF] Ben-Israel | Generalized Inverses: Theory and Applications[END_REF], or possibly some fast computation {1,2,3}-inverse (Theorem 5 from [START_REF] Courrieu | Fast solving of Weighted Pairing Least-Squares systems[END_REF]), for instance. It can be useful to remember the four Penrose equations for the generalized inverse matrices. Let M (k) denote a {k}-inverse of the matrix M, then:

MM (1) M = M, M (2) MM (2) = M (2) , (MM (3) )'= MM (3) , (M (4) M)'= M (4) M.

We will also use the Kronecker (or tensor) product, denoted A⊗B, of matrices A and B, and the column-wise vectorization, denoted vec(A), of a matrix. The notation [xi]i=1..n corresponds to a column vector of n components (the xi's).

Problem solution

One-dimensional solution

Theorem 1. Problem ( 1) is equivalent to the maximal eigen-element problem:

XX (1,2,3) TT'y = λmax y, with ||y|| = 1, and the solution of problem (1) is:

w = X (1,2,3) y, with Σk=1..d r 2 (Xw, Tk) = λmax.
Proof. Given a vector v∈R n , the kth component of the vector T'Xv is equal to:

< Xv, Tk > = || Xv || r(Xv, Tk), thus v'X'TT'Xv = || Xv || 2 Σk=1..d r 2 (Xv, Tk). ( 2
)
As a consequence, problem (1) is equivalent to:

w = arg maxv∈R n v'X'TT'Xv, subject to v'X'Xv = 1, (3)
which, at first glance, is a quadratically constrained quadratic programming problem.

Setting z = Xv, one can reformulate problem (3) as: y = arg maxz∈R m z'TT'z, subject to ||z|| = 1, and XX (1,2,3)

z = z, (4)
where the last constraint expresses that z must belong to the range of X, and thus there exists v such that Xv = z, that is v = X (1,2,3) z. The matrix XX (1,2,3) is an orthogonal projector, it is idempotent (equal to its square) and symmetric [START_REF] Ben-Israel | Generalized Inverses: Theory and Applications[END_REF].

We can now solve problem (4) using the solution of the maximal eigen-element problem of Theorem 1:

with ||y|| = 1,
XX (1,2,3) TT'y = λmax y ⇒ λmax XX (1,2,3) y = XX (1,2,3) XX (1,2,3) TT'y = XX (1,2,3) TT'y ⇒ y = XX (1,2,3) y (5)

⇒ y'TT'y = y'XX (1,2,3) TT'y = λmax y'y = λmax

On the other hand, using [START_REF] Brillinger | Time Series: Data Analysis and Theory[END_REF]: (1,2,3) y ⇒ w = X (1,2,3)

y = XX
y ⇒ y = Xw ⇒ w'X'TT'Xw = Σk=1..d r 2 (Xw, Tk) = λmax ,
(using ( 2) and ( 6))

which completes the proof. ☐

Practical considerations

In practice, there are various ways of solving the maximal eigen-element problem of Theorem 1. However, given that we need only the dominant eigenvalue and the corresponding eigenvector, while the problem m×m matrix, say A = XX (1,2,3) TT', can be huge if the number m of items is very large, it is preferable to use a suitable method based on Krylov subspaces and Arnoldi iterations [START_REF] Lehoucq | ARPACK Users Guide: Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods[END_REF][START_REF] Sorensen | Numerical methods for large eigenvalue problems[END_REF]. However, even in this case, the computation of the m×m matrix A can be too heavy. Fortunately, the computational complexity can be lowered exploiting the fact that the number n of regressors is usually not very large.

First, we note that in several {1,2,3}-inverse matrix computation methods, X (1,2,3) can be written as X (1,2,3) =HX', where H is a symmetric n×n real matrix. For instance one can find an H matrix for computing the Moore-Penrose inverse of X in [START_REF] Courrieu | Fast computation of Moore-Penrose inverse matrices[END_REF], while Theorem 5 from [START_REF] Courrieu | Fast solving of Weighted Pairing Least-Squares systems[END_REF] provides an H matrix for computing another type of {1,2,3}-inverse.

The first thing to do is to compute the matrix H for the chosen type of {1,2,3}-inverse of X, which is usually a reasonable cost operation. Eigenvalue methods based on Krylov subspaces use only matrix-vector products, avoiding matrix-matrix products that are very expensive in the case of huge matrices. This principle can be extended to the computation of the problem matrix itself, given that this matrix is used only in matrixvector products. Let z denote the vector in a matrix-vector product to be computed: (1,2,3) TT'z = X(H(X'(T(T'z)))). [START_REF] Cohen | Applied Multiple Regression/Correlation Analysis for the Behavioral Sciences[END_REF] This way, the computation includes only matrix-vector products, and the possibly huge matrix A is in fact never computed.

Az = XX
One can also note that the sign of the solution vector w depends on the sign of the eigenvector y, and thus it is arbitrary. This determines arbitrarily the signs of the correlation coefficients of Xw with the columns of T, thus the sign of w can possibly be modified in order to satisfy various criteria for the application. For instance, one can choose the sign of the largest correlation, or the sign of the subset of correlations having the same sign and the greatest sum of squares.

Associated least-squares solution

Well-known methods such as the reduced-rank regressions do not use a correlation criterion to be maximized as in [START_REF] Anderson | Estimating linear restrictions on regression coefficients for multivariate normal distributions[END_REF], but various weighted least-squares criteria to be minimized [START_REF] Izenman | Reduced-rank regression for the multivariate linear model[END_REF][START_REF] Izenman | Modern Multivariate Statistical Techniques[END_REF]. Contrarily to the correlation, the quadratic error is not invariant to linear transformations of the variates, thus in order to minimize it in the rank-one case, the one-dimensional solution of Theorem 1 is not sufficient, and it must be weighted with a specific coefficient while a specific bias coefficient must be append for each target variate. Fortunately, given the one-dimensional solution w of Theorem 1, the required coefficients are very easy to compute. Let µ Ξ ∈ R n denote the row vector of Ξ column means, let µ Θ ∈ R d denote the row vector of Θ column means, and let Y be the matrix Θ whose columns have been centred (but not normalized), then we must first solve:

b' = arg minv∈R d || Xwv' -Y || 2 ,
where || . || denotes the Frobenius norm. The solution is simply the vector:

b' = (Xw)'Y/|| Xw || 2 .
The row vector µ ∈ R d of bias coefficients is given by:

µ = µ Θ -µ Ξ (wb').
Finally, the matrix of least-squares regression coefficient B ∈ R (n+1) × d is given by:

B = [µ', (wb')']',
where [. , .] here denotes the concatenation operator. The approximation of Θ is:

Θ * = [1m , Ξ ]B,
and the mean square error (MSE) that will be used in certain numerical test is:

MSE = || Θ * -Θ || 2 / (md).
Note that the coefficients in b not only change the norm of the regression coefficient vectors, but they can also change their sign independently for each variate. As a result, if one computes the correlations between the columns of Θ * and those of Θ, all these correlations have the same positive sign, while possible opposite behaviors of different target variates are hidden. This is harmful for applications using correlation statistics, thus in this case the correlations must be computed between the onedimensional approximation Xw and the columns of Θ (or equivalently of T), which preserves the magnitude of the correlations as well as the correlation sign variations.

Note also that a similar problem occurs in usual reduced-rank regression methods, where the regression coefficients are always in matrix form, even in the rank-one case. This requires applying a sign correction procedure when one needs to use correlation statistics.

Relation with the usual R 2 statistic

A squared Pearson correlation coefficient whose one of the arguments is a linear combination of regressors, as in (1), is similar to a squared multiple correlation coefficient, and it is equal to the well-known R 2 statistic if, and only if, the linear combination of regressors is optimal in the least-squares sense:

uk = arg minv∈R n || Xv -Tk || 2 , 1≤k≤d , (8)
which has solutions of the form uk = X (1,2,3) Tk. In this case, the set of optimal regression coefficients is an n×d matrix U = X (1,2,3) T, with one column of regression coefficients per target variate. This provides the maximum possible value to each R 2 (X, Tk), k = 1..d.

However, in the case of a multi-target regression (1), there is only one vector w of regression coefficients, and we must clarify the relation of w with the columns of U. With w, λmax and y as in Theorem 1, we have the following result: 3) T, and c = (λmax) -1 T'y , then:

Lemma 1. Set U = X (1,2 ,
If λmax > 0, then w = Uc,
Proof. Using the maximal eigen-element of Theorem 1, we have:

XX (1,2,3) TT'y = λmax y ⇒ XX (1,2,3) Tc = y ⇒ X (1,2,3) XX (1,2,3) Tc = X (1,2,3) y = w.
On the other hand, since X (1,2,3) is a {2}-inverse, we have:

X (1,2,3) XX (1,2,3) Tc = X (1,2,3) Tc = Uc,
and thus Uc = w , which proves Lemma 1. ☐

Corollary 1. R 2 (X, Tk) = r 2 (Xuk , Tk) ≥ r 2 (Xw, Tk), k = 1..d .
Proof. The inequality in Corollary 1 results from the fact that the linear combination Uc of the columns of U in Lemma 1 is not necessarily optimal in the least squares sense for the target variates considered individually. ☐

A consequence of this is that the distribution of the squared correlations in the multitarget case (d>1) does not reduce to the distribution of the well-known R 2 statistic, and thus the F-test associated with the R 2 statistic [START_REF] Cohen | Applied Multiple Regression/Correlation Analysis for the Behavioral Sciences[END_REF] is not suitable to test the null hypothesis of the r 2 statistics in this case.

Lemma 2. ck = r(Xw, Tk) / (Σj=1..d r 2 (Xw, Tj)), k = 1..d.
Proof. This is straightforward:

c = (λmax) -1 T'y = (Σj=1..d r 2 (Xw, Tj)) -1 T'Xw = [r(Xw, Tk)]k = 1..d / (Σj=1..d r 2 (Xw, Tj)). ☐

Statistical significance of the multi-target regression correlation coefficients

Computing p-values

As we have seen in Section 3 (Corollary 1), one cannot use significance tests associated with usual multiple correlation coefficients to test the significance of multi-target correlation coefficients. Given that usual parametric tests are not suitable, we will turn to robust distribution-free methods such as permutation tests.

One knows that, in order to build a valid permutation test, it is necessary and sufficient to have exchangeable observations under the null hypothesis [START_REF] Good | Extensions of the concept of exchangeability and their applications[END_REF][START_REF] Good | Permutation, Parametric, and Bootstrap Tests of Hypotheses[END_REF][START_REF] Legendre | Numerical ecology, 2nd English edition[END_REF]. It is usually the case of the items in regression problems, and in particular in problem [START_REF] Anderson | Estimating linear restrictions on regression coefficients for multivariate normal distributions[END_REF]. As we shall see in section 5.1, the used statistical model allows the kth target variate to be written as:

Tk = αkXω + ξk , k=1..d
where αkXω is the regression part (αk ∈ R, ω ∈ R n ), and ξk ∈ R m is the residual whose components are assumed to be independent and identically distributed (with mean 0 and variance σ 2 ), which implies that these residue components are exchangeable. On the other hand, the null hypothesis for the kth target variate is:

H0: αk = 0.
Thus, under the null hypothesis, the linear relation of the target with the regressors vanishes, and it remains only the exchangeable residue components (one per item).

Exploiting this, one can use the data to perform repeated Monte-Carlo simulations to generate a large sample from the distribution of multi-target r 2 statistics under the null hypothesis. As we have seen, if the null hypothesis is true, all possible pairings of the regressors values with the target variates values are equally likely to occur. Thus it suffices to repeatedly randomly permute the rows (items) of the matrix X of regressors, or equivalently (but not simultaneously), to repeatedly randomly permute the rows of the target matrix T, then to compute and store the resulting multi-target r 2 statistics. Note, however, that the r 2 distributions associated to the different target variates are not necessarily the same, and we must consider a specific distribution for each target variate. The estimated p-value for the kth r 2 statistic is the proportion of r 2 statistics in the kth permutation test distribution that are greater or equal to the kth r 2 statistic obtained without permutation.

In order to solve problem (1) for each permutation, let P be a random permutation matrix of order m, then [START_REF] Cohen | Applied Multiple Regression/Correlation Analysis for the Behavioral Sciences[END_REF] becomes:

Apz = PXX (1,2,3) P'TT'z = Xp(H((Xp)'(T(T'z)))), with Xp = PX (9a)
or alternatively:

Apz = XX (1,2,3) PTT'P'z = X(H(X'(Tp((Tp)'z)))), with Tp = PT (9b)
Note, however, that in practice one simply permutes the rows of X or of T using an index permutation operator (what is denoted Xp or Tp), not products with permutation matrices that would be computationally heavy. This procedure (with the option (9b)) is implemented in Matlab/Octave code as the function "MTRegPV" listed in the Appendix.

In addition, the listed function "MTRegLS" computes the associated least-squares solution if necessary.

Controlling for the Family Wise Error Rate

When the number d of target variates is large, there is a potential problem of test inflation that requires a control procedure. A priori, the target variates are related in some way, and the dependences of tests can be of any type. So we must use a method not requiring restrictive hypotheses on the dependences of the tests. This is the case of methods derived from Bonferroni's correction, in particular the Family Wise Error Rate (FWER) control of Holm-Bonferroni [START_REF] Holm | A Simple sequentially rejective multiple test procedure[END_REF]. The FWER can be controlled using, for instance, the p-values computed by permutation tests as described in Section 4.1. The FWER control procedure is implemented in Matlab/Octave code as the function "Holm" listed in the Appendix, and it is optionally called by the function "MTRegPV" if the input argument "fwerc" of "MTRegPV" is set to a strictly positive value.

Statistical model and consistency of the estimator

Statistical model

As usually described in mathematical statistics handbooks (e.g. [START_REF] Montfort | Cours de Statistique Mathématique[END_REF] pp. 237-257), the multiple linear regression model could be summarized as:

Tk = Xυk + εk , k=1..d ( 10
)
where υk ∈ R n is the true vector of regression coefficients for the kth target variate, and εk ∈ R m is a vector of residues.

One usually retains the following hypotheses:

H1: E(εk) = 0,

H2: E(εk εk') = σ 2 Im,
where E is the expected value operator, and Im the identity matrix of order m.

Proving the consistency of the estimators also requires the additional hypothesis:

H3: rank(X) = n, which supposes that m>n, and the n regressors are linearly independent.

In what concerns the multi-target linear regression model, it is exactly the same as above, except that the regression coefficients must be decomposed as:

υk = αkω + δk, k=1..d ( 11
)
where ω ∈ R n is the true vector of regression coefficients common to the d target variates, αk = ω' υk/(ω'ω) is the projection coefficient of υk on ω, and δk ∈ R n is the residual rejection vector of the kth target variate.

Substituting [START_REF] Courrieu | Validated intraclass correlation statistics to test item performance models[END_REF] in [START_REF] Courrieu | Density Codes, Shape Spaces, and Reading. ERMITES 2012 : Representations and Decisions in Cognitive Vision[END_REF], one obtains:

Tk = X(αkω + δk) + εk = αkXω + (Xδk + εk) = αkXω + ξk , k=1..d (12)
where the quantity Xδk is transferred from the regression to the residue.

Concerning the multi-target case, while H3 remains unchanged, we must modify H1 and H2 in the following way. Define the global residue ξ ∈ R md as:

ξ = vec(ξ1, ξ2, …, ξd) .
This is the residue associated to the target vector t = vec(T). Modify the hypotheses as:

H1: E(ξ) = 0 , H2: E(ξ ξ') = σ 2 Imd .
The new versions of the hypotheses imply the old ones. In addition, the new version of H2 supposes that the target variates are linearly independent at the population level.

Note that ω is not defined if the distribution of αk's concentrates on zero (δ(0) Dirac distribution). In this case, all multi-target correlations are random, but not necessarily zero (except if all target variates are orthogonal to all regressors). In such a case, we say that the regression problem (1) is "insubstantial".

Definition. We say that the multi-target regression problem (1) is "not insubstantial" if there exists a population parameter ω ∈ R n and a sample of d coefficients {αk ; k=1..d}, as defined in [START_REF] Courrieu | Validated intraclass correlation statistics to test item performance models[END_REF], belonging to a population whose probability distribution does not reduce to the Dirac distribution on zero δ(0). This implies that a non-zero proportion of the coefficient population has a non-zero value.

Intuitively, this means that the considered regressors are not completely irrelevant for the considered population of target variates.

Consistency of the estimator

Theorem 2. Under the hypotheses H1-H3, consider w ∈ R n as defined in Theorem 1, αk and ω as defined in [START_REF] Courrieu | Validated intraclass correlation statistics to test item performance models[END_REF], and assume that the problem (1) is not insubstantial. Then:

(i) E(w) = ω, (ii) limd →∞ E( || w -ω || 2 ) = 0.
If in addition the regressor variates have finite moments of order 1 and 2, then:

(iii) limm →∞ E( || w -ω || 2 ) = 0.
Proof. After H3, X is of full column rank, thus all {1,2,3}-inverses of X are equivalent, and one can simply consider X † = (X'X) -1 X'.

Let c ∈ R d be defined as in Lemma 1, t = vec(T) ∈ R dm , and define M ∈ R dm×n as:

M = (c') † ⊗ X.
Then consider the following linear least squares system:

Mv ≈ t ⇔ v = M † t = (c' ⊗ X † ) vec(T) = X † Tc = w.
Thus, if the vector c is given, then we can transform (1) into an equivalent linear least squares system of matrix M and single target t.

H1 implies that E(t) = Mω, and one has:

E(w) = (M'M) -1 M'E(t) = (M'M) -1 M'Mω = ω,
which proves (i).

Moreover, let Σ(w) denote the covariance matrix of w, then, given H2:

Σ(w) = (M'M) -1 M'(σ 2 I)[ (M'M) -1 M']' = σ 2 (M'M) -1 M'M(M'M) -1 = σ 2 (M'M) -1 ,
for some variance σ 2 of the residues. Now, since E(w) = ω, we have:

E( || w -ω || 2 ) = Trace(Σ(w)) = σ 2 Trace((M'M) -1 ).
Using Lemma 2, one easily verify that (c') † = [r(Xw, Tk)]k = 1..d , and thus:

(M'M) -1 = (((c') † ⊗ X)'((c') † ⊗ X)) -1 = ((Σj=1..d r 2 (Xw, Tj))X'X) -1 = (Σj=1..d r 2 (Xw, Tj)) -1 (X'X) -1 .
This implies that:

E( || w -ω || 2 ) = σ 2 Trace((X'X) -1 )/ (Σj=1..d r 2 (Xw, Tj)), ( 13
)
where the numerator does not depend on d, and the denominator (= λmax) tends to infinity as d tends to infinity, because the problem (1) is not insubstantial, by hypothesis, thus the number of non-zero coefficients αk's tends to increase proportionally to d and the resulting sum of squared correlations tends to infinity. This proves (ii).

Given that the regressors are centred, X'X/m is a convergent estimator of the true covariance matrix V of the regressor variates, which is regular after H3. Thus:

limm →∞ X'X/m = V ⇒ limm →∞ (X'X/m) -1 = V -1 ⇒ limm →∞ m(X'X) -1 = V -1 ⇒ limm →∞ Trace((X'X) -1 ) = limm →∞ Trace(V -1 )/m = 0.
Substituting this in [START_REF] Courrieu | General or idiosyncratic item effects: what is the good target for models[END_REF], one obtains (iii), which completes the proof of Theorem 2. ☐

Remark. As noted in section 2.2, the sign of the eigenvector y in Theorem 1 is arbitrary, that is, if y is a solution, then -y is also a solution. As a result, the sign of w, the sign of c (Lemma 1), and the sign of the correlations (Lemma 2) are also arbitrary. Now, considering the unknown vector ω as the asymptotic solution of problem (1), it is clear that the sign of ω is also arbitrary. This indeterminacy is lifted in the above proof because the vector c is used in the calculation of the matrix M, thus the sign of ω is implicitly the sign corresponding to that of c and to that of w.

Note that [START_REF] Courrieu | General or idiosyncratic item effects: what is the good target for models[END_REF] shows how and why increasing the number of correlated linearly independent target variates (d) enhances the convergence of the estimator (w). Now, if rank(X) < n , for instance because m ≤ n, then H3 is not met and Theorem 2 does not apply. However, Theorem 1 remains valid and the method still works, but the estimator is generally biased since E(w) = M (1,2,3) E(t) = M (1,2,3) Mω , which is usually different from ω when M is not of full column rank.

Relation with reduced-rank regression methods

Let X ∈ R m×n and Y ∈ R m×d be the centred columns versions of the regressor and target matrices, respectively, and let T ∈ R m×d be the normalized columns version of Y, as previously.

Definition 1. Set sk = a|| Yk || , k=1..d, for any fixed real a>0. Then define S ∈ R d×d as the diagonal matrix whose diagonal elements are the sk 's, k=1..d.

Theorem 3.

The following rank-one reduced-rank regression problem with S -2 weighting:

find v ∈ R n , and q ∈ R d that minimize Trace{(Xvq' -Y)S -2 (Xvq' -Y)'} can be solved using Theorem 1 since it is equivalent to :

(i) find v ∈ R n that maximizes v'X'TT'Xv, subject to v'X'Xv = 1, (ii) then set q' = (Xv)'Y .

Conversely, let Q ∈ R n×d be the matrix of regression coefficients provided by the rankone RRR with S -2 weighting, and let Qk be its kth column, 1≤k≤d, then:

(iii) v = ± Qk/|| XQk || .
Proof. First, suppose that one knows v, then the solution of:

arg minq ∈ R d Trace{(Xvq' -Y)S -2 (Xvq' -Y)'} = arg minq ∈ R d || (Xvq' -Y)S -1 || 2
is simply given by:

q' = (Xv)'Y/|| Xv || 2 = (Xv)'Y,
which proves (ii).

Next, one must determine v. One has: minv,q Trace{(Xvq' -Y)S -2 (Xvq' -Y)'} = minv,q Trace{ S -1 (Xvq' -Y)'(Xvq' -Y) S -1 } = minv,q Trace{ S -1 qv'X'Xvq'S -1 -S -1 Y'Xvq'S -1 -S -1 qv'X'YS

-1 + S -1 Y'YS -1 } = minv Trace{ S -1 || Xv || -2 Y'Xv || Xv || 2 v'X'Y || Xv || -2 S -1 -S -1 Y'Xvv'X'Y || Xv || -2 S -1 -S -1 || Xv || -2 Y'Xvv'X'YS -1 + S -1 Y'YS -1 } = minv Trace{S -1 || Xv || -2 Y'Xvv'X'Y S -1 -2 S -1 || Xv || -2 Y'Xvv'X'YS -1 + S -1 Y'YS -1 } = minv Trace{ S -1 Y'YS -1 -|| Xv || -2 S -1 Y'Xvv'X'YS -1 } = minv Trace{ a -2 T'T -a -2 || Xv || -2 T'Xvv'X'T} = minv Trace{ a -2 TT' -a -2 || Xv || -2 v'X'TT'Xv} ,
whose solution is given by:

arg maxv ∈ R n (v'X'TT'Xv), subject to || Xv || 2 = 1,
which is similar to (i) and to equation (3), and thus can be solved using Theorem 1.

Finally, if Q = vq', then all columns of Q are collinear to v, while the constraint || Xv || = 1 implies (iii). Note, however, that the solution sign is always arbitrary. ☐ Corollary 2. The consistency properties stated in Theorem 2 for the method of Theorem 1 also apply to the rank-one reduced-rank regression with a weight matrix of type S -2 .

Proof. This immediately follows from Theorem 3 since the two considered methods are equivalent. ☐

Computational tests

Generating artificial data

Using artificial data allows us to know the exact solution of problems and thus to compare the results provided by various methods with the exact ones. In the present case, we use the following model to generate each column of the target matrix.

Tk = αkXω + βkb + γkgk, k=1..d ( 14
)
where X ∈ R m×n is the matrix of regressors, ω ∈ R n is the vector of exact regression coefficients, b ∈ R m with b'X=0 is a shared factor not accountable for by the regression, and gk ∈ R m is a random Gaussian component of mean 0. The column vectors Xω, b, and gk have unit norms and the real coefficients αk, βk, and γk are such that αk 2 + βk 2 + γk 2 =1.

Varying αk allows one to approximately control the kth exact correlation value r(Xω, Tk), but this one must be precisely computed afterwards.

Typical observations

We used [START_REF] Giraud | Low rank multivariate regression[END_REF] to generate data with m=300, n=30, d=30, the αk's were uniformly sampled in various intervals, and the βk's were set to zero, the γk's providing the complements to the αk's. We studied the multi-target regression r-values, computed with the "MTRegPV" function listed in the Appendix, and the classical multiple regression R values, as functions of the exact r values. The statistical significance was established at the risk α = 0.05, using classical parametric F-tests for the exact r's and the multiple regression R's, while we used permutation tests as described in section 4.1

for the multi-target regression r-values. The FWER control was not used for the basic studies, but it was controlled afterwards. substantially overestimated the exact positive correlations, as well as the absolute value of negative ones (remember that the R coefficient is not signed). There were 25 significant exact correlation coefficients, 22 significant multi-target correlation coefficients, and 18 significant multiple regression R coefficients, without FWER control.

After FWER control, it remained respectively 22, 16, and 16 significant coefficients.

The results presented in Figure 2 were obtained with the αk's sampled in the interval [0, 0], that is, all coefficients were random. In this case, significant coefficients are type I errors. Without FWER control, 2 exact correlation coefficients were significant (6.7%) at the 0.05 risk, only 1 multi-target correlation coefficient was significant (3.3%), and 3 multiple regression R's were significant (10%). With FWER control, there was no longer significant coefficient. coefficients as functions of the exact r coefficients. The data were generated using the model ( 14) with αk = 0, βk = 0, k=1..30. The α-risk was set to 0.05 and there was no FWER control.

Effect of the number of target variates

In order to see how the number of target variates and the number of items influence the regression accuracy in the framework of the data-generating model ( 14 Multi-Target Regression r significant MTR r Multiple Regression R significant MR R significant exact r and analysed 100 independent random data sets, comparing the obtained regression coefficients (w) with the exact ones (ω), and comparing the multi-target r coefficient values together with their significance to those of the exact r coefficients.

Instead of being arbitrarily chosen, the sign of w was determined in such a way that w'ω ≥ 0. Then the accuracy of w was measured by the error norm || wω ||.

The accuracy of the estimated r coefficients was measured using their "discrepancy" (D) from the exact r-values:

D(rMT, rexact) = maxk=1..d |rMT(k) -rexact(k)| .
The obtained statistical significance (at the .05 level) was also compared to that of the exact coefficients without FWER control. We computed the "significance agreement" defined as the proportion of target variates for which both the estimate and the exact coefficient lead to the same conclusion (significant or non-significant).

Table 1 shows the averages over 100 tests of these accuracy measures as functions of the number of target variates and of the number of items. The average sums of squared correlations (λmax) and the average computation times (in seconds) for the multi-target regressions are also reported (for Matlab 9.4 on Mac OS X 10.11.6). For each performance measure, a one-factor standard analysis of variance was performed to test the effect of the number of target variates (d). The obtained significances are reported in the last column of Table 1. As one can see in Table 1, regardless of m, both the w error and the rMT discrepancy significantly decreased, converging towards 0, while the significance agreement converged towards 1, as the number of target variates increased.

In the same time both λmax and the computation time increased. As expected, the accuracy of the estimates was better with 320 items than with 160 items.

Table 1. Variation of the averages (over 100 tests) of the accuracy measures for the multi-target linear regression, as functions of the number of target variates (d) and the number of items (m). The data were generated using the model ( 14 Table 2 shows the averages over 100 tests of the accuracy measures as functions of the number of target variates and of the number of regressors. Large values, greater than the number of items, have been selected for the number of targets. One of the two numbers of regressors is greater than the number of items, resulting in a violation of the hypothesis H3. As one can see in Table 2, regardless of n, both the w error and the rMT discrepancy significantly decreased as the number of target variates increased. However, we a priori know that in the case where n>m the estimator is biased, resulting in a non-zero error limit. We also observe in this case that the significance agreement is quite bad and does not improve as d increases.

Table 2. Variation of the averages (over 100 tests) of the accuracy measures for the multi-target linear regression, as functions of the number of target variates (d) and the number of regressors (n). The data were generated using the model ( 14 Globally, all these observations perfectly illustrate Theorem 2. In summary, one can say that in order to obtain a good rank-one estimate, one must use as many items and target variates as possible, while always having a number of items greater than the number of regressors and using linearly independent regressors. We also know that the target variates must be correlated but linearly independent at the population level. For instance, it is intuitively clear that if one artificially increases the number of targets by repeating many times the same small set of variates, this provides no more information, and there is no reason that this improves the convergence. Note however that, at the data sample level, the target vectors (T columns) cannot be linearly Table 3. Numerical test on artificial data with m=615, n=24, d=100, and the actual rank of the true regression coefficient matrix is 1 or 3. The compared methods are the MTR (Theorem 1), the rank-one RRR or the 5-folds, cross validated rank (CV-rank) RRR, both with the S -2 weight matrix, the identity weight matrix (I), and the cov(Y) -1 weight matrix (Canonical). The table reports the sum of squared correlations ( Σr 2 ), the mean squared error (MSE), the obtained solution rank, and the Bayesian Information Criterion (BIC).

Comparison with reduced-rank regressions on real data

In this section, we compare the performance of the same methods as in the previous section, but on a real data regression problem whose actual rank is a priori unknown. The data are z-scores of speeded naming latencies of 615 printed French words, provided by 100 human participants. These data were previously studied in [START_REF] Courrieu | Validated intraclass correlation statistics to test item performance models[END_REF][START_REF] Courrieu | Missing data imputation and corrected statistics for large-scale behavioral databases[END_REF]. The purpose here is to study the direct effect of the orthographic form of the words on the naming latencies. In order to do this, we need a suitable numerical representation of the character strings that could be used as a multidimensional regressor. A possible numerical coding of character strings has previously been proposed in [START_REF] Courrieu | Density Codes, Shape Spaces, and Reading. ERMITES 2012 : Representations and Decisions in Cognitive Vision[END_REF] as the output layer of artificial neural networks for handwriting recognition, and was also used in [START_REF] Rey | Accounting for variance in single-word ERPs[END_REF] as a multidimensional orthographic regressor to analyse cerebral event related potentials (EEG/ERP) in a printed word naming task. The coding principle is as follows.

Consider an alphabet of N symbols {s1, s2, ..., sN}, for instance the 26 lower-case letters of the Roman alphabet. The coding associates to each symbol of the alphabet one component of a real vector (c1, c2, ..., cN). Let χ be a symbol string of L characters, one first determines the "symbol position bits" as bk,i = 1 if the symbol si appears at rank k in χ, else one has bk,i = 0. Then the components of the orthographic code are given by:

ci(χ) = (Σk=1..L bk,i 2 -k ) p , i=1..N, 0 < p ≤1, ( 15
)
where p is a free parameter (we use p=1/3). Such a numerical code unequivocally represents the corresponding character string and it can always be decoded back into this string.

The necessary alphabet for describing the 615 experimental words had only 24 letters since "k" and "w" are very rare in French, and these two letters never appeared in the corpus. Moreover, the used words did not include letters with diacritic marks. Thus the regressor matrix was a 615×24 real matrix of numerical string codes [START_REF] Good | Extensions of the concept of exchangeability and their applications[END_REF].

One can observe in Table 4 that the rank-one RRR with identity matrix weighting provided the same performance as those of the MTR and the rank-one S -2 weighted RRR. This is because the data are z-scores, thus all columns of the target matrix have the same norm and the identity matrix, in this case, is also a S -2 matrix. Optimized rank RRRs detected a rank-two case, however, the BIC suggests that the best choice is in fact the MTR model and its equivalent rank-one RRRs. In other words, the improvement of the MSE does not compensate the increase of model complexity with rank-two solutions, and these observations finally suggest that a rank-one model is the most suitable. In fact, one knows that this type of data is essentially described by a general rank-one model, but with an additional small amount of non-random idiosyncratic effect that could explain the detected rank-two solution [START_REF] Courrieu | General or idiosyncratic item effects: what is the good target for models[END_REF]. A common practice is to calculate the average response vector over all participants, and then to fit this vector with a linear combination of regressors, using an ordinary multiple regression. Doing this, one implicitly assumes an underlying model where the αk's in [START_REF] Courrieu | Missing data imputation and corrected statistics for large-scale behavioral databases[END_REF] are equal for all participants, which is possibly unrealistic and can lead to problematic results. Using the multi-target regression avoids this drawback and can provide a much better solution.

Observing the estimator convergence with real data

With real data, such as those used in the previous section, one does not know the exact solution. However one can exploit the fact that if two independent estimates converge on the same solution, then they also converge on each other. The idea here is that if the used real data fulfil the assumptions of Theorem 2, then one should observe the predicted convergence as a function of d and as a function of m.

Using the regressors and the behavioural data of section 7.5, we repeatedly sampled pairs of independent random samples of d = 10, 20, 30, 40, or 50 participants, we computed the regression solution (Theorem 1) for each of the two samples (w1 and w2), and the Euclidean distance between these two solutions. This was repeated 100 times for each modality of d, and a standard analysis of variance was performed on the solution distances. In addition, the same was done not varying d but the number m of items, with m = 61, 122, 183, 244, or 305. The results are reported in Table 5, where one can clearly observe the convergence as a function of d and as a function of m.

Table 5. Average distance of two independent estimates of the regression coefficient vector w (from Theorem 1), as a function of the number d of target variates (participants), and of the number m of items (words). The regressors are the 24 dimensions of a numerical orthographic code [START_REF] Good | Extensions of the concept of exchangeability and their applications[END_REF], and the target data are z-scores of printed word naming times. Each test was repeated 100 times, and a standard analysis of variance was performed for each factor (d and m).

Conclusion

We have defined a new method for solving multivariate linear regression problems, with only one shared vector of regression coefficients. This way, only what is common to the target variates is fitted while the sum of squared correlations of the fitted pattern with the target variates is maximized. This is equivalent to a rank-one variant of Izenman's reduced rank regression method [START_REF] Izenman | Reduced-rank regression for the multivariate linear model[END_REF][START_REF] Izenman | Modern Multivariate Statistical Techniques[END_REF], with a weighting proportional to a diagonal inverse-variance matrix. However, the formulation of the new method allowed us to derive certain new convergence properties of interest for practical applications. In particular, one clearly shows how and why the estimation of the regression coefficients converges as a function of the number of target variates, as well as a function of the number of items. Of course, these consistency properties also concern the equivalent rank-one, reduced-rank regression method, however, we don't know whether or not the convergence as a function of the number of target variates concerns other RRR approaches. A ready to use Matlab/Octave code program is provided for practical applications, and numerical tests on artificial data as well as on real data fully corroborated the formal statements.

Figure 1

 1 Figure 1 shows a typical result with the αk's sampled in the interval [-0.4, 0.7]. As

Figure 1 .Figure 2 .

 12 Figure 1. Multi-target linear regression r coefficients and classical multiple regression R coefficients as functions of the exact r coefficients. The data were generated using the model (14) with αk ∈ [-0.4, 0.7], βk = 0, k=1..30. The α-risk was set to 0.05 and there was no FWER control.

  ), we varied the parameter d from 8 to 128 (in powers of 2), the parameter m ∈ {160, 320}, while the other parameters were set to n=24, αk ∈ [-0.4, 0.7], βk = 0, k=1..d, and the α-risk was set to 0.05. For each number of target variates and of items, we generated

  ) with n=24, αk ∈ [-0.4, 0.7], βk = 0, k=1..d. The α-risk was set to 0.05 and there was no FWER control. Number of target variates: d = 8 d = 16 d = 32 d = 64 d =128 d effect significance m = 160, n = 24 λmax || wω || rMT discrepancy from exact r Significance agreement Computation time (seconds) m = 320, n = 24 λmax || wω || rMT discrepancy from exact r

  ) with m=64, αk ∈ [-0.4, 0.7], βk = 0, k=1..d. The α-risk was set to 0.05 and there was no FWER control. Number of target variates: d=64 d=128 d=256 d=512 d=1024 d effect significance m = 64, n = 24 λmax || wω || rMT discrepancy from exact r Significance agreement Computation time (seconds) m =64, n =96 (H3 violation) λmax || wω || rMT discrepancy from exact r

Table 4 .

 4 Comparison of 7 multivariate regression methods on the analysis of real zscores of naming latencies of m=615 printed words, by d=100 human participants, regressed on numerical orthographic codes of the words with n=24 components. The 100 correlation coefficients provided by MTR (Theorem 1) and equivalent methods are in the range [0.1306, 0.5044] and significant at the 0.05 level, with a mean r of 0.3225.

	Method	Σr 2	MSE	rank	BIC
	MTR	10.9564	0.8890	1	-0.0956
	Rank-1 RRR:
	S -2 weight	10.9564	0.8890	1	-0.0956
	I weight	10.9564	0.8890	1	-0.0956
	Canonical	10.1679	0.8969	1	-0.0868
	CV-rank RRR:
	S -2 weight	11.7183	0.8814	2	-0.0825
	I weight	11.7183	0.8814	2	-0.0825
	Canonical	11.4674	0.8839	2	-0.0797

Funding: This work, carried out within the Labex BLRI (ANR-11-LABX-0036) and the Institut Convergence ILCB (ANR-16-CONV-0002), has benefited from support from the French government, managed by the French National Agency for Research (ANR) and the Excellence Initiative of Aix-Marseille University (A*MIDEX)

Data Analysis, 52, 406-422. doi:10.1016/j.csda.2007.01.025

independent if the number of targets is greater than the number of items, but this does not prevent the convergence, as illustrated in Table 2.

Comparison with reduced-rank regressions on controlled rank artificial data

In this section, we compare the behaviour of the proposed rank-one regression method with that of conventional reduced-rank regression methods, on rank-one and rank-three artificial data. The data were generated using a random 615×24 matrix of regressors and a random 24×100 coefficient matrix whose rank was controlled (using the singular value decomposition) to be 1 or 3. Target variates were obtained multiplying the regressors by the coefficient matrix and adding Gaussian noise. An available Matlab implementation of Izenman's reduced-rank regression (RRR) was used to test alternative methods [START_REF] Mccomb | Reduced Rank Regression[END_REF]. We used the rank-one RRR and the 5-folds crossvalidation optimized rank RRR, both with the S -2 weight matrix, the identity weight matrix (I), and the cov(Y) -1 weight matrix ("canonical"). Table 3 reports, for each case, the sum of squared correlations ( Σr 2 ), the mean squared error (MSE), the obtained solution rank, and the Bayesian Information Criterion (BIC). The BIC is commonly used as a model selection criterion [START_REF] Reinsel | Multivariate Reduced-Rank Regression: Theory and Applications[END_REF][START_REF] Schwarz | Estimating the dimension of a model[END_REF], while the model having the lowest BIC must be preferred.

As one can see in Table 3, the new method (noted "MTR") and the rank-one RRR with S -2 weighting have exactly the same performance (as expected from Theorem 3).

For rank-one data, these methods provided the maximum sum of squared correlations, while the minimum MSE was obtained with the identity weighting, which also resulted in the minimum BIC. For rank-three data, the optimized rank RRRs detected the appropriate rank, leading these methods to the best performance, including the BIC.

Appendix. Matlab/Octave code of useful functions (for academic use only)

function [sr2,psr2,r,p,s,w] = MTRegPV(X,T,fwerc,a) % Multi-Target Multiple Linear Regression % For m items, n regressors, and d target variates % p-values are computed by permutation test % Input arguments: % X: m x n matrix of regressors % T: m x d matrix of target variates % fwerc: control the FWER [START_REF] Anderson | Estimating linear restrictions on regression coefficients for multivariate normal distributions[END_REF] or not (0) % a: chosen alpha risk (default: a=0.05) % Output arguments: % sr2: sum of squared correlations % psr2: p-value of sr2 (by permutation test) % r: d x 1 vector of r estimates % p: d x 1 vector of p-values (by permutation test) % s: d x 1 significance vector (1=signif., 0=n.s.) % w: n x 1 vector of regression coefficients % - ----------------------------------------------- -----------------------------------------------------[