A partial graphical model with a structural prior on the direct links between predictors and responses
Eunice Okome Obiang Pascal Jézéquel Frédéric Proïa
math.ST stat.AP stat.TH


A partial graphical model with prior on the direct links]
A partial graphical model with a structural prior on the direct links between predictors and responses


E. Okome Obiang]Eunice Okome Obiang
Univ Angers, CNRS, LAREMA, SFR MATHSTIC, F-49000 Angers, France.
okome@math.univ-angers.fr
P. Jézéquel]Pascal Jézéquel
1 Unité de Bioinfomique, Institut de Cancérologie de l'Ouest, Bd Jacques Monod, 44805 Saint Herblain Cedex, France.

2 SIRIC ILIAD, Nantes, Angers, France. 

3 CRCINA, INSERM, CNRS, Université de Nantes, Université d'Angers, Institut de Recherche en Santé-Université de Nantes, 8 Quai Moncousu - BP 70721, 44007, Nantes Cedex 1, France.
pascal.jezequel@ico.unicancer.fr
F. Proïa]Frédéric Proïa
Univ Angers, CNRS, LAREMA, SFR MATHSTIC, F-49000 Angers, France.
frederic.proia@univ-angers.fr






This paper is devoted to the estimation of a partial graphical model with a structural Bayesian penalization. Precisely, we are interested in the linear regression setting where the estimation is made through the direct links between potentially high-dimensional predictors and multiple responses, since it is known that Gaussian graphical models enable to exhibit direct links only, whereas coefficients in linear regressions contain both direct and indirect relations (due e.g. to strong correlations among the variables). A smooth penalty reflecting a generalized Gaussian Bayesian prior on the covariates is added, either enforcing patterns (like row structures) in the direct links or regulating the joint influence of predictors. We give a theoretical guarantee for our method, taking the form of an upper bound on the estimation error arising with high probability, provided that the model is suitably regularized. Empirical studies on synthetic data and a real dataset are conducted.


[
    [
    
=====








AMS 2020 subject classifications: Primary 62A09, 62F30; Secondary 62J05.







§ INTRODUCTION



We are interested in the recovery and estimation of direct links between high-dimensional predictors and a set of responses. Whereas the graphical models seem a natural way to go, we propose to take account of a prior knowledge on the predictors, when possible. This is typically the case when dealing with genetic markers whose joint influence may be anticipated thanks to some kind of genetic distance, or when the predictors are supposed to represent a continuous phenomenon so that consecutive covariates probably act together. In this regard, while taking up the graphical approach, we introduce some Bayesian information in a structural regularization of the estimation procedure, although the inference remains frequentist, thereby following the idea of Chiquet et al. <cit.>. This strategy also enables to affect the amount of shrinkage by playing with some hyperparametrization in the prior, while sparsity may be obtained via usual penalty-based patterns. Regarding the mathematical formalization of the graphical models that we will just briefly discuss in this introduction, we refer the reader to the very complete handbook recently edited by Maathuis et al. <cit.>. We also refer the reader to the book of Hastie et al. <cit.> and to the one of Giraud <cit.>, both related to the standard high-dimensional statistical methods. Before introducing the model and the organization of this work, let us describe the notation used throughout the paper.



 §.§ Notation


For any matrix A, | A |_* = ‖(⃗A) ‖_* is the elementwise ℓ_* norm of A and | A |_*^- is | A |_* deprived of the diagonal terms of A. We also note ‖ A ‖_F = | A |_2 the Frobenius norm of A and ‖ A ‖_2 the spectral norm of A. The Frobenius inner product between any matrices A and B of same dimensions is A, B  = ⟨(⃗A), (⃗B) ⟩ = (A^t  B) whereas ⟨ u, v ⟩ = u^t  v is the inner product of the Euclidean real space. For any vector u, | u |_0 is the number of non-zero values in u. For a matrix A, [A]_C is to be understood as the matrix A whose elements outside of the set of coordinates C are set to zero and (⃗A) is the vectorization of A into a column vector. The eigenvalues of a square matrix A of size d with spectrum sp(A) are λ_i(A) taken in decreasing order (from λ_1(A) = (A) to λ_d(A) = (A)). The cones of symmetric positive semi-definite and definite matrices of dimension d are _+^  d and _++^  d respectively.



 §.§ The partial graphical model


In the classic Gaussian graphical model (GGM) setting, we aim at estimating the precision matrix Ω = Σ^-1 of jointly normally distributed random vectors Y ∈^q and X ∈^p with zero mean and covariance Σ. The point is that it induces a graphical structure among the variables and the support of Ω is closely related to the conditional interdependences between them. Let us consider, now and in all the study, the sample covariances of n independent observations (Y_i, X_i), denoted by

    = 1/n∑_i=1^n Y_i  Y_i^t,  = 1/n∑_i=1^n Y_i  X_i^tand = 1/n∑_i=1^n X_i  X_i^t.

Maximizing the penalized likelihood of a GGM boils down to finding Ω∈_++^  p+q that minimizes the convex objective

    L_n(Ω) = -ln(Ω) +  S^  (n), Ω + λ pen(Ω)

where S^  (n) is the full sample covariance built from the blocks (<ref>). The penalty function pen(Ω) is usually |Ω|_1 or even |Ω|_1^-. Efficient algorithms exist to get solutions for (<ref>), see e.g. Banerjee et al. <cit.>, Yuan and Lin <cit.>, Lu <cit.> or the graphical Lasso of Friedman et al. <cit.>. The reader may also look at the theoretical guarantees of Ravikumar et al. <cit.>. However, thinking at X_i as a predictor of size p associated with a response Y_i of size q, the partial Gaussian graphical model (PGGM), developped e.g. by Sohn and Kim <cit.> or Yuan and Zhang <cit.>, appears as a powerful tool to exhibit direct relationships between the predictors and the responses. To understand this, consider the decomposition into blocks

    Ω = [  ;   ]andΣ = [  ;   ]

where ∈_++^  q, ∈^q × p and ∈_++^  p and where the same goes for ,  and . The precision matrix Ω = Σ^-1 satisfies, by blockwise inversion,

    =  -   and = -( -   )^-1  .

The conditional distribution peculiar to Gaussian vectors

    Y_i |  X_i ∼ (-   X_i, )

gives a new light on the multiple-output regression Y_i = B^  t  X_i + E_i with Gaussian noise E_i∼(0, R), through the reparametrization B = - and R =. Whereas B contains direct and indirect links between the predictors and the responses (due e.g. to strong correlations among the variables),  only contains direct links, as it is shown by the graphical models theory. In other words, the direct links are closely related to the concept of partial correlations between X and Y (see Meinshausen and Bülmann <cit.> or Peng et al. <cit.>, for the univariate case). For example, the direct link between predictor k and response ℓ may be evaluated through the partial correlation (Y_ℓ, X_k |  Y_≠ ℓ, X_≠  k) contained, apart from a multiplicative coefficient, in the ℓ-th row and k-th column of  (see e.g. Cor. A.6 in <cit.>) with the particularly interesting consequence that the support of  is sufficient to identify direct relationships between X and Y. Hence, in the partial setting, the objective reduces to the estimation of the direct links  together with the conditional precision matrix of the responses . Maximizing the penalized conditional log-likelihood of the model now comes down to minimizing the new convex objective

    L_n(, )   =   -ln() + ,  + 2 , 
     + ,    + λ pen() + μ pen()

over (, ) ∈_++^  q×^q × p for some usual penalty functions. It is worth noting that pen() often plays a crucial role in modern statistics dealing with high-dimensional predictors (and the natural choice is ||_1 to get sparsity) while we may choose λ=0, for the responses are hardly numerous. In the seminal papers <cit.> and <cit.>, the authors consider ||_1 and ||_1^- for pen(), respectively. Yuan and Zhang <cit.> also point out that no estimation of  is needed anymore. In a graphical model, the estimation of  and  depends on the accuracy of the estimation of Ω which, in turn, is strongly affected by the one of , especially in a high-dimensional setting. The partial model overrides this issue, the focus is on  and  while  has disappeared from the objective function (<ref>). The latter is obtained either by considering the multiple-output Gaussian regression scheme, or, as it is done in <cit.>, by eliminating  thanks to a first optimization step in (<ref>). In this paper, we will consider the penalties

    pen() = ||_1^-andpen() = ||_1

which correspond to the PGGM (Gm) of <cit.>. The Spring (Spr) of <cit.> can also be seen as a PGGM but with no penalty on  (replaced with an additional structuring one on , we will come back to this point thereafter), so for Spr we may consider λ=0. The generalized procedure (GenGm) at the heart of the study relies on a combination between these two approaches. We will see in due time that we keep both the penalties of Gm and the structuring one of Spr on . Finally, the intermediate solution consisting in estimating  and B through the conditional distribution Y_i |  X_i∼(B^  t  X_i, ) with penalizations both on B and , presented and analyzed by Rothman et al. <cit.> and by Lee and Liu <cit.>, is better known as a multivariate regression with covariance estimation (MRCE). However, it has been shown that the objective function suffers from a lack of convexity and that the optimization procedure may be debatable, in addition to the less convenient setup for statistical interpretation (B contains both direct and indirect influences) compared to PGGM. Without  claiming  to  be  exhaustive, let us conclude this quick introduction by citing some related works, like the structural generalization of the Elastic-Net of Slawski et al. <cit.>, the Dantzig approach of Cai et al. <cit.> put in practice on genomic data <cit.>, the greedy research of the non-zero pattern in Ω of Johnson et al. <cit.>, the approach of Fan et al. <cit.> using a non-convex SCAD penalty to reduce the bias of the Lasso in the estimation of Ω, the eQTL data analysis of Yin and Li <cit.> which makes use of a sparse conditional GGM, and so on. All the references inside will complete this concise list.



 §.§ Organization of the paper


To sum up, we have two goals in this paper :


  * Give some theoretical guarantees to the (slightly modified) model introduced in Chiquet et al. <cit.>.

  * Generalize the result of Yuan and Zhang <cit.> to the case where a structural penalization is added in the estimation step.

In Section <ref>, we introduce the model, consisting in putting a generalized Gaussian prior on the direct links before the procedure of estimation of  and , and we detail the new convex objective. Then we provide some error bounds for our estimates, useful as theoretical guarantees of performance. Section <ref> is devoted to empirical considerations. We explain how we deal with the minimization of the new objective and we test the method on simulations first, and next on a real dataset (a Canadian average annual weather cycle, see e.g. <cit.>). After a short conclusion in Section <ref>, we finally prove our results in Section <ref>. The numerous constants appearing in the results and the proofs are gathered in the Appendix, for the sake of readability.



§ A GENERALIZED GAUSSIAN PRIOR ON THE DIRECT LINKS



We use the definition given in formulas (1)-(2) of <cit.> for the so-called d-dimensional multivariate generalized Gaussian (0, 1, V, β) distribution with mean 0, scale 1, scatter parameter V ∈_++^  d and shape parameter β > 0. According to the authors, the density takes the form of

    ∀  z ∈^d,  f_V, β(z) = β Γ(d/2)/π^d/2 Γ(d/2 β)  2^d/2 β√((V)) exp(-⟨ z, V^-1 z ⟩^β/2)

where Γ is the Euler Gamma function.



We clearly recognize the Gaussian (0, V) setting for β=1. Moreover, for β=1/2, it can be seen as a multivariate Laplace distribution whereas it is known to converge to some uniform distribution as β→ +∞. The marginal shapes (d=1 and V=1) of the distribution are represented on Figure <ref>, depending on whether β < 1, β=1 or β > 1. Our results hold for all β≥ 1 but, as will be explained in due course, we shall not theoretically deviate too much from the Gaussianity in the prior (even if we will allow ourselves some exceptions in the practical works). The usual Bayesian approach for multiple-output Gaussian regression having B as matrix of coefficients and R as noise variance consists in a conjugate prior (⃗B) ∼(b, R ⊗ L^-1) for some information matrix L ∈_++^  p and a centering value b (see e.g. Sec. 2.8.5 of <cit.>). In the PGGM reformulation, we have R = and B = - as explained in Section <ref>, and of course we shall choose b=0 to meet our purposes. Thus, 

    (⃗) = -(⊗ I_p) (⃗B) ∼ (0, ⊗ L^-1)

is a natural prior for the direct links (this is in particular the choice of the authors of <cit.>). Following the same logic, let us choose ⊗ L^-1 for scatter parameter and suppose that

    (⃗) ∼ (0, 1, ⊗ L^-1, β).

In this way, we can play on the intensity of the constraint we want to bring on , from a non-informative prior to quasi-boundedness through Laplace and Gaussian distributions. This prior entails an additional smooth term acting as a structural penalization in the objective (<ref>) that becomes

    L_n(, )   =   -ln() + ,  + 2 , 
     + ,    + η  L, ^  t  ^β + λ ||_1^- + μ ||_1

with three regularization parameters (λ, μ, η). The smooth penalization lends weight to the prior on  and thereby plays on the extent of shrinkage and structuring through β, whereas ||_1 and ||_1^- are designed to induce sparsity. One can note that this is closely related to the log-likelihood of a hierarchical model of the form

    {[ Y_i |  X_i,  ∼ (-   X_i, );    (⃗) ∼ (0, 1, ⊗ L^-1, β) ].

where the emphasis is on  in the prior and  remains a fixed parameter, although it is important to see that, in this work, the estimation step does not rely on a posterior distribution. The following proposition is related to the existence of a global minimum for our objective (<ref>) with respect to (, ) as soon as β≥ 1.


Assume that β≥ 1. Then, L_n(, ) defined in (<ref>) is jointly convex with respect to (, ).


See Section <ref>.


Now and throughout the rest of the paper, denote by θ = (, ) ∈Θ = _++^  q×^q × p the (q × (q+p))-matrix of parameters of the model, with true value = (, ). As it is usually done in studies implying sparsity, we will also consider S of cardinality | S |, the true active set of  defined as S = { (i,j), _i,j≠ 0 }, and its complement S̅. Our results also depends on some basic assumptions related to the true covariances of the Gaussian observations, and we will assume that the following holds.

    H_1∈^  p_++, ∈^  q_++,  B ≠ 0 (that is, ≠ 0)and  L ∈_++^  q.

This is a natural hypothesis in our framework, in particular we suppose that there is at least a link between X and Y.



Even if it is of less interest, our study does not exclude the case where =0. Indeed, we might as well consider that =0 and get the same results, but some constants should be refined. On the other hand, ∈^  p_++ and ∈^  q_++ are crucial.


Under (<ref>), the random matrices

    A_n = ( - ) -    ( - )  with h_a = | A_n|_∞

and

    B_n = 2  (( - ) +    ( - )) with h_b = | B_n|_∞

are going to play a fundamental role, especially  h_a and h_b. Let us now provide some theoretical guarantees for the estimation of θ in our model, provided that the regularization parameters are located in a particular area (λ, μ, η) ∈Λ. Consider the penalized likelihood ℓ_λ, μ, η(θ) given in (<ref>), and estimate θ by the global minimum

    = min_Θ ℓ_λ, μ, η(θ)

obtained for β≥ 1. To facilitate reading, we postpone the precise definition of the numerous constants to the Appendix. We recall that p is the number of predictors, q is the number of responses and | S | is the size of the true active set.


Fix d_λ > c_λ > 1, d_μ > c_μ > 1, e_λ > 0 and e_μ > 0, and assume that the regularization parameters satisfy (λ, μ, η) ∈Λ = [c_λ  h_a,  d_λ  h_a] × [c_μ  h_b,  d_μ  h_b] × [0, η], where

    η = min{(c_λ-1) λ/c_λ ℓ_a, (c_μ-1) μ/c_μ ℓ_b, e_λ  h_a/ℓ_a, e_μ  h_b/ℓ_b}/β  s_L^ β-1

for some non-random constants s_L, ℓ_a and ℓ_b defined in (<ref>) and (<ref>), and the random constants h_a and h_b given above. Then, under (<ref>), there exists absolute constants b_1 > 0 and b_2 > 0 such that, for any 0 < b_3 < 1 and as soon as n > n_0, with probability no less that 1 - ^-b_2 n-b_3, the estimator (<ref>) satisfies

    ‖ - ‖_F ≤ 16  m^*  c_λ, μ√(| S |)/γ_r,η,β,p √(ln(10 (p+q)^2) - ln(b_3)/n)

where γ_r,η,β,p, c_λ, μ and m^* are technical constants defined in (<ref>), (<ref>) and (<ref>), respectively, and where the minimal number of observations is given by

    n_0     =    max{(ln(10 (p+q)^2) - ln(b_3))  c_λ, μ^2 | S |  (16  m^*)^2/r^*  2 γ_r,η,β,p^2, 
           b_1  (q + ⌈ s_α⌉ln(p+q)), ln(10 (p+q)^2) - ln(b_3) }

with s_α defined in (<ref>) and r^* in (<ref>).


See Section <ref>.


Among all these constants, we can note that s_L, ℓ_a, ℓ_b, h_a and h_b are useful to properly describe and restrict Λ, the domain of validity of (λ, μ, η) for the theorem to hold. Once Λ is fixed, the other constants take part in the upper bound of the estimation error. However, as it stands, the theorem is very difficult to interpret. The next two remarks seem essential to have an overview of the orders of magnitude involved for the number of observations, for p and q, for the estimation error and for the regularization parameters. 



Of course the degree of sparsity | S | is crucial in the estimation error, but it also plays an indirect role in the probability associated with the theorem and in the numerous constants. In virtue of Lemma <ref>, we can hope that λ and μ have a wide validity band, by playing on c_λ, c_μ, d_λ and d_μ. In turn, η also has a non-negligible area of validity, provided of course that ℓ_a, ℓ_b and s_L, all depending on combinations between ,  and L, are small enough. Accordingly, it would be to our advantage if L was both sparse and not chosen with too large elements. As it always appears together with η, we may as well take a normalized version of L (e.g. | L |_∞≤ 1).




Even if the result holds for any β≥ 1, the terms ∝ p^ β-1 appearing in some upper bounds of the proof clearly argue in favor of a moderate choice β∈ [1, 1+ϵ] for a small ϵ > 0, depending on p. In other words, we cannot deviate too much from the Gaussianity in the prior on the direct links. For example in a very high-dimensional setting (p ∼ 10^7), choosing ϵ = 0.1 leads to p^ β-1≈ 5 whereas we may try larger values of ϵ for the more common high-dimensional settings p ∼ 10^3 or p ∼ 10^4. By contrast, we can see that n_0 must (at least) grow like q for the theorem to hold, so high-dimensional responses are excluded. However in multiple-output regressions, even when p is extremely large, q generally remains small. According to all these considerations, we may roughly say that, in a high-dimensional setting with respect to p,

    ‖ - ‖_F ≲ √(| S | ln p/n)

with a large probability, under a suitable regularization of the model. We recognize the usual terms appearing in the error bounds of regressions with high-dimensional covariates, like the ℓ_2 error of the Lasso (see e.g. Chap. 11 of <cit.>). This is the same bound as in <cit.>, but our additional structural penalty restricts Λ.




§ SIMULATIONS AND REAL DATASET



The minimization problem (<ref>) is solved using a coordinate descent procedure, alternating between the computations of

    = min__++^  q ℓ_λ, μ, η(, ) and = min_^q × p ℓ_λ, μ, η(, ).

Each step is done by an Orthant-Wise Limited-Memory Quasi-Newton (OWL-QN) algorithm (see e.g. <cit.>). The first subproblem is performed through half-vectorization () to ensure symmetry and we set the objective to +∞ on _++^  q to ensure positive definiteness of the solution. The coordinate descent is stopped when

    ‖^  (t) - ^  (t-1)‖_2 ≤ ϵmax(1, ‖^  (t-1)‖_2) and‖^  (t) - ^  (t-1)‖_2 ≤ ϵmax(1, ‖^  (t-1)‖_2)

following two consecutive iterations t-1 and t, where ϵ > 0 is a small threshold depending on the desired precision. We are now going to try our method on synthetic data first, and then on a real dataset. We will pay attention to the role played by β, in particular we will see that it can be useful as well as counterproductive, depending on the situations.



 §.§ Simulations

 For each scenario, we first generate i.i.d. standard Gaussian vectors X_i∈^p, then Y_i∈^q is simulated according to the setting and we estimate  and . From the relations detailed in Section <ref>, we recall that Y_i = B^  t  X_i + E_i with E_i∼(0, R) is an equivalent formulation, provided that B = - and R =. In a compact form, we may also write

    Y = X B + E or(⃗Y) = (I_q⊗ X) (⃗B) + (⃗E)

where the i-th row of Y is Y_i^  t and the i-th row of X is X_i^  t. Thus, we can estimate B using the Lasso (Las) and the Group-Lasso (GLas) in the vectorized form, to provide a basis for comparison between our method and the usual penalized methods. The Lasso penalty is obviously ‖(⃗B) ‖_1 to promote coordinate sparsity while, for the Group-Lasso, we use the penalty ‖ B_1‖_2 +  + ‖ B_p‖_2 where B_i is the i-th row of B, to promote row sparsity and exclude altogether some predictors from the model. We also implement some variants of our generalized graphical model (GenGm). The case where = R^-1 is known and does not need to be estimated is the Oracle (Or) and the case where η=0 so that β has no influence is the classic PGGM (Gm). The case where λ=0 and β=1 is called the Spring (Spr) by the authors of <cit.>. We will focus on structured scenarios. With no structure in , there is no reason why our method should outperform the usual PGGM. In a completely random setting, we have observed that all PGGM procedures perform identically. In fact, a slight gain can be obtained compared to Spr and Gm simply due to the flexibility induced by the additional parameter (Spr and Gm are particular cases of GenGm). However, that clearly cannot counterbalance the extended computational times, and GenGm should not be used for such situations. The calibration of the regularization parameters is made using a cross-validation on a training set of size n_t = 150 and the accuracy is evaluated thanks to the mean squared prediction error (MSPE) on a validation set of size n_v = 1000,

    MSPE = ‖ Y + X ^  t ^  -1‖_F^2/q  n_v.

Due to the large amount of treatments, the grids for cross-validation are not very sharp here but they will be carefully refined for the real datasets of the next section. The covariance between the outputs is R = (r^ | i-j |)_1 ≤  i,  j ≤  q for r = 1/2 and we work with p=100. Each scenario is repeated N = 500 times and GenGm is evaluated with numerous values of β, from 0.25 to 2 with a step of 0.25. The results of the following scenarios are summarized on Figures <ref>, <ref> and <ref> below, respectively.


  * Scenario 1 (q=1). We draw ω_i = ±1/2 for i=1, , 10 and we fill 10 randomly selected sections of size 3 in  with ω_i. The remaining part of  is 0.

  * Scenario 2 (q=2). We draw ω = ±1/2 and one randomly selected row of  is filled with ω while the other is identically 0.

  * Scenario 3 (q=3). We draw ω_i = ±1/2 and we fill a randomly selected section of size 30 on the i-th row of  with ω_i, for i=1,2,3. The remaining part of  is 0.


The row structure is promoted by a normalized first finite difference operator

    L = 1/2[   1 - 1   0       0;  -1   2   ⋱   ⋱   ⋮;   0   ⋱   ⋱   ⋱   0;   ⋮   ⋱   ⋱   2  -1;   0       0  -1   1 ]

which, through L, tends to penalize the difference between two consecutive values on a same row (as does Fused-Lasso with ℓ_1 penalty). Yet, the Fused-Lasso is not a suitable alternative to GLas and Las in this precise context because B = - is not supposed to have a row structure even if  has one. For this choice of L, one can note that, in the particular case where R = (σ^  2_1, , σ^  2_q),

    L, ^  t  ^β = ( ∑_i=1^qσ^  2_i∑_j=2^p (ω_i,j - ω_i,j-1)^2 )^β≥∑_i=1^qσ_i^  2 β∑_j=2^p|ω_i,j - ω_i,j-1|^  2 β

where ω_i,j is the (i,j)-th element of , so we may fairly expect that β≥ 1 is going to strengthen the smoothness of the estimation and to enforce all the more the structuring.



We could as well add a small diagonal element in the matrix L defined above, positive semi-definite but not invertible. The resulting effect would be a negligible ridge-like penalization on the elements of . This is not required for the estimation procedure but useful for Theorem <ref> to hold (see e.g. (<ref>)). Likewise, it seemed interesting to test some settings with β < 1 even if the theory developped in the paper does not give any guarantee for them, as a basis for comparison.








First of all, one can observe that Las and GLas are left behind in all our simulations. This is not surprising since the covariance between the outputs cannot be recovered with the standard Lasso, at least for q ≥ 2. Generally, GLas remains more robust compared to Las, probably due to the high level of sparsity in  approximately passed to B (provided that the covariances in R are small enough), and exploited by the grouping effect. In the weakly structured setting (Scenario 1), we also observe that, as expected, all PGGM procedures perform almost identically, with obviously an advantage for Or (although small, illustrating the accuracy of the estimation). In the strongly structured settings (Scenarios 2 and 3), Gm gives results below the expected level, because it is not designed to promote such layouts. On the contrary, thanks to this choice of L showing here geat efficiency, GenGm and Spr are doing pretty well. Note that, in this context, GenGm with β=1 is almost the same as Spr since, q being small, λ does not play a crucial role. However, some empirical facts draw our attention: the prediction error decreases with β to some extent, but the most interesting fact seems to be the simultaneous decrease of its variance. It is likely that the increasing pressure exerted by β on the estimation procedure leads to a higher homogeneity in the numerical results, despite the repetitions of random experiments under random settings. In other words, the structuring seems to be strengthened and we also observe an acceleration in the procedure of estimation that logically follows from the latter remarks (especially clear when we compare β=0.25 and β=2). On the other hand, for the opposite reason, we notice that the predictions are hardly better than Gm (even worse in some cases), both on average and in terms of variability, for β < 1, and these simulations tend to undermine such values of the hyperparameter. On the whole, GenGm with β > 1 might be a sound approach for practitioners who place a high priority on structuring the estimations, even if Remark <ref> below should probably temper this statement. To conclude, let us consider the strongly structured scenarios with L = I_p (without structuring) in the Oracle setting with β=2, and let us compare the results with those of Figures <ref> and <ref>, obtained with the correct version of L given in (<ref>). The results are displayed on Figure <ref> where we can see that the benefit of structuring is manifest. Unsurprisingly, the results without structuring are close to those of Gm since L=I_p only strengthens the shrinkage effect with ridge-like additional penalties.





To estimate (, ) in the model Spr, the authors of <cit.> use a very judicious and efficient method relying, in each step of the coordinate descent procedure, on a direct computation of the estimation of  together with an Elastic-Net estimation of . This is possible for λ=0 and β=1, but unfortunately cannot be implemented in the general setting. As a result, computational times remain an issue that should be paid attention to.




The mean value of the estimation errors ‖ - ‖_F^2 leads to the same kind of observations for the models being compared in the simulations. But the minimal prediction error does not always coincide with an optimal support recovery due to the shrinkage effect on the estimation of . The so-called F-score is given by

    F = 2  p_r  r_e/p_r+r_ewhere p_r = TP/TP + FPand r_e = TP/TP + FN

are the precision and the recall, respectively, and where T/F and P/N stand for true/false and positive/negative. In the strongly structured scenarios, F is generally located between 0.60 and 0.65, and a deeper analysis shows that a proportion of more than 0.99 of true non-zero values are recovered (that is, the part of the true active set S related to ). If the models are not calibrated to reach the best prediction error but the best F-score, F regularly exceeds 0.90, at least for the structured procedures.


Nevertheless, Scenarios 2 and 3 are very strongly structured, more than one would expect from an unknown underlying generating process, and the real dataset of the next section is going to highlight the fact that the improvement may be hardly noticeable with respect to β, even in terms of prediction error variance. But we will see that β can still be useful for variable selection.



 §.§ A real dataset



The dataset available as  in the  package  contains daily temperature and precipitation at 35 different locations in Canada, averaged over annual reports starting in 1960 and ending in 1994 (see e.g. <cit.>). We intend to look at the direct links between the minimal and maximal rainfall (on the log_10 scale) and the temperature pattern in the 35 weather stations, so as to identify the times of the year that have a strong effect on rainfall (positive as well as negative). In this context, n=35, q=2 and p=365. Figure <ref> shows temperature and log-precipitation measured over a year in Montreal, chosen as an example, together with the empirical distribution of the minimal and maximal log-precipitation for the 35 weather stations. We can note that, since the data are averaged over numerous years, outliers are unlikely even for the extremes (min and max).



Some authors (see e.g. <cit.>) have already highlighted the pertinence of using the matrix L defined in (<ref>) in this dataset, because the predictors are ordered temporally so that the selection of isolated days instead of relevant sequences of days seems an unreliable procedure for statistical interpretation. To assess the models, we repeat N=100 times the following experiment: n_t = 25 observations are randomly selected for calibration (via 2-fold cross-validation) and estimation, the remaining n_v=10 observations are used to compute the MSPE (<ref>) related to the prediction of the minimum (min_p) and maximum (max_p) precipitation. We can see on Figure <ref> that all structured PGGM perform almost identically, with the phenomenon described in the previous section still visible but to a lesser extent. We can even notice that structuring is hardly beneficial for this dataset, from a purely numerical point of view. This conclusion can also be found in <cit.>, where the author compares the structured Elastic-Net with unstructured alternatives to predict the 0.25-, 0.50- and 0.75-quantiles of the log-precipitation. But we will see that, in terms of variable selection and statistical interpretation, L and β still have a substancial role to play.



The point is that the best prediction error does not usually coincide with a sparse solution (see Remark <ref> above). On the basis of the MSPE, most of the time we must retain μ≪ 10^-2 and only a few direct links are set to zero. To look for sequences of days directly related to min_p and max_p, let us constraint μ≥ 10^-2 and focus on variable selection. The active set of  is evaluated on the basis of n_t = 25 randomly chosen observations. The experiment is repeated N=100 times, and the locations having a frequency of occurrence that exceeds 0.5 are retained (or, equivalently, those whose estimates have a non-zero median). This can be seen as a measure of variable importance. The results are given on Figures <ref> and <ref> below for min_p and max_p, respectively, with a fixed set of regularization parameters and increasing values of β. The objective is to show the influence of the latter, all other things being equal. Note that, since we retain λ=0 in these experiments, GenGm for β=1 coincides with Spr.





We observe that the increasing pressure exerted by β on the estimation procedure tends to refine the selection by giving priority to the most important variables and by dropping the others much more easily, at the cost of prediction results: 
we are undoubtedly in a selection process. The sequence of inclusions

    S_β_2⊂S_β_1for β_1 < β_2

that we observe for the estimated active sets is clearly a guarantee of quality for the selected variables. The median values of the estimated direct links between the temperature of the days and the pair (min_p, max_p) are represented on Figure <ref> together with the estimated regression coefficients, for β=2. We detect sequences of influent days in November, December, January and February, especially related to min_p, positively at the end of the year and negatively at the beggining. This is broadly consistent with the analysis of <cit.> – even if the responses are not extremes but quantiles in this reference – with however two differences: the regression coefficients associated with max_p are much lower, and an activity is also detected between July and August. The main explanation, at least for the first of them, probably lies in the use of graphical models that take into account the correlation between responses. Indeed, as can be seen on Figure <ref> which gives an overview of the estimation of R obtained from the repeated experiments, a non-zero correlation is detected between the responses (≈ 0.32). The influence of November and December on all quantiles and that of January and February on the 0.75-quantile in <cit.> might actually be an artificial effect of the correlation with the 0.25-quantile. This is what our study suggests by highlighting min_p compared to max_p. From this point of view, the interest of graphical models is particularly obvious.





Let us also mention that, interestingly enough, we notice that the role of η tends to depreciate for the large values of β. For example, for the same regularization parameters (λ,μ)=(0, 0.05) and β=2, the difference between the estimated active sets for η=0.1 and η=1 is almost negligible (depending on the experiments, between 1 and 3 days are concerned, on average). Based on these studies and observations, we might conclude that β is insignificant when we are interested in the best prediction error on a validation set (even counterproductive with respect to computational times, e.g. compared to Spr), whereas it seems to have a substancial role to play when focusing on selection, by accelerating the discrimination of variables. In the first case, η has to be carefully adjusted while in the second case, β will quickly help to reach the desired sparsity.



For the simulations and the real dataset, we have used the popular first finite difference operator given in (<ref>). Other examples can be found in the literature, like the promotion of a genetic distance for genomic selection in Brassica napus <cit.> or the bidimensional discretization of the Laplacian to work on handwritten digit recognition <cit.>. More generally, L can be used in a classic Bayesian prior supposed to promote some covariance structure on the direct links, with no `physical' structuring in mind (like temporal, spatial or genetic proximity).




§ CONCLUSION



In conclusion, our work is a generalization of <cit.>, using the same technical tools to establish an upper bound on the estimation error when a prior on the direct links generates an additional structural penalty in the objective, provided that the model is suitably regularized. Our work is also an improvement of <cit.> since, while being inspired by the methodology of the authors, we generalize the prior and give some theoretical guarantees. The empirical study shows that the hyperparametrization in the prior, although more expensive in adjusting the parameters, is likely to refine the selection results but clearly, this does not appear as a crucial improvement compared to the two previous points. Let us conclude the paper by highlighting two weaknesses that might be trails for future studies. On the one hand, the Laplace distribution is often used as a prior in the Bayesian Lasso (see e.g. Sec. 6.1 of <cit.>). However, our reasonings do not allow β=1/2, which may correspond to a multivariate Laplace distribution on the direct links. Combined with the first finite difference operator L, the choice β=1/2 could generate a Fused-Lasso-type penalty. In this regard, it would be challenging and interesting to obtain some theoretical guarantees for β≥ 1/2 and not only for β≥ 1, even if our probably too brief simulation study does not encourage the choice of β < 1. On the other hand, λ=0 is a natural choice when q is small (this is in particular the configuration of <cit.>), not to mention that it is computionally faster. But, the proof of our theorem needs λ > c_λ  h_a > 0 to hold. We think that a reasoning enabling to deal with λ = 0 should also be beneficial to the study. More generally, it would be instructive to consider a very high-dimensional setting (p ≫ n and not only p ∼ 10^2 although always larger than n, as in our experiments). Such studies should follow with omic data.



Acknowledgements and Fundings. The authors warmly thank the two anonymous reviewers for the careful reading and for making numerous useful corrections to improve the paper. We thank ALM (Angers Loire Métropole) and the ICO (Institut de Cancérologie de l'Ouest) for the financial support. This work is partially financed through the ALM grant and the “Programme opérationnel régional FEDER-FSE Pays de la Loire 2014-2020" noPL0015129 (EPICURE). The authors also thank Mario Campone (project leader and director of the ICO), Mathilde Colombié (scientific coordinator of EPICURE clinical trial) and Fadwa Ben Azzouz, biomathematician in Bioinfomics, for the initiation, the coordination and the smooth running of the project.



§ TECHNICAL PROOFS



We start in a first part by some useful linear algebra lemmas that will be repeatedly used subsequently, well-known for most of them. In a second part, we prove the joint convexity of the objective and our main result.



 §.§ Linear algebra





Let A ∈_+^  d and U ∈^d ×ℓ. Then, U^  t A U ∈_+^ ℓ.


Since A is symmetric with non-negative eigenvalues, there is an orthogonal matrix P such that A = P D P^  t with D = (sp(A))∈_+^  d. Thus, for all v ∈^ℓ, it follows that ⟨ v, U^  t A  U  v ⟩ = ‖ D^1/2  P^  t  U  v ‖^2 ≥ 0. 




Let A ∈_++^  d and B ∈_+^  d. Then for all i, λ_i(AB) ≥ 0.


The equality A B = A^1/2  (A^1/2  B  A^1/2)  A^-1/2 shows that A B and A^1/2  B  A^1/2 are similar, so they must share the same eigenvalues. From Lemma <ref>, λ_i(A^1/2  B  A^1/2) ≥ 0 .




Let A ∈_+^  d and B ∈_+^  d. Then,

    (A) (B) ≤ (A B) ≤ (A) (B).



Since A - (A) I_d∈_+^  d and B ∈_+^  d,

    ((A - (A) I_d)  B) = (B^1/2  (A - (A) I_d)  B^1/2) ≥ 0

from Lemma <ref>, thus (A B) ≥(A) (B). The other inequality is obtained through (A) I_d - A ∈_+^  d.




Let A ∈_++^  d and B ∈_+^  d. Then,

    (A) (B) ≤ (A B) and(A B) ≤ (A) (B).



On the one hand, (A B) ≤‖ A B ‖_2 ≤‖ A ‖_2 ‖ B ‖_2 = (A) (B), since A and B are symmetric and since, from Lemma <ref> and by hypothesis, all eigenvalues appearing in the relation are non-negative. Suppose now that B is invertible so that both A^-1 and B^-1 belong to _++^  d. Then, ( (A B)^-1 ) ≤(A^-1) (B^-1) and this immediately gives (A B) ≥(A) (B). If B is not invertible, the relation trivially holds since we still have (A B) ≥ 0 from Lemma <ref>.




Let A ∈_+^  d and U ∈^d ×ℓ. Then,

    (A) ‖ U ‖_F^2 ≤ (U^  t A U) ≤ (A) ‖ U ‖_F^2.



Denote by u_i the i-th column of U. It is not hard to see that the i-th diagonal element of U^  t A U is u_i^  t A  u_i≥(A) ‖ u_i‖^2 ≥ 0. Thus,

    (U^  t A U) = ∑_i=1^ℓ u_i^  t A  u_i≥(A) ∑_i=1^ℓ‖ u_i‖^2 = (A) ‖ U ‖_F^2.

The upper bound stems from 0 ≤ u_i^  t A  u_i≤(A) ‖ u_i‖^2.





Let A and B be symmetric matrices of same dimensions. Then,

    (A) + (B) ≤ (A+B) and(A+B) ≤ (A) + (B).



These are just two special cases of Weyl inequalities. We refer the reader to Thm. 4.3.1 of <cit.>, for example.




 §.§ Convexity of the objective



We know from Prop. 1 of <cit.> and the convexity of the elementwise ℓ_1 norm that L_n(, ) - η  L,   ^β is itself convex, but it remains to show that this is still the case with the additional smooth penalty.



  §.§.§ Proof of Proposition <ref>

Recall that Θ = _++^  q×^q × p and consider the mapping Φ : Θ→_+^  p defined as

    ∀  (A,B) ∈Θ, Φ(A, B) = B^  t A^-1 B.

We can already note from Lemma <ref> that (Φ(A, B)) ≥ 0. Moreover, for all 0 ≤ h ≤ 1 and all Z_i = (A_i,B_i) ∈Θ, i=1,2, it is easy to see that

    S_h(Z_1,Z_2) = h Φ(Z_1) + (1-h) Φ(Z_2) - Φ(h Z_1 + (1-h) Z_2)

is the Schur complement of h A_1 + (1-h) A_2 in the matrix

    M_h(Z_1, Z_2) = h [                 A_1                 B_1;             B_1^  t B_1^  t  A_1^-1 B_1 ] + (1-h) [                 A_2                 B_2;             B_2^  t B_2^  t  A_2^-1 B_2 ].

But the decomposition

    [     A^1/2 A^-1/2  B;         0         0 ]^ t[     A^1/2 A^-1/2  B;         0         0 ] = [             A             B;         B^  t B^  t  A^-1 B ]

directly shows that M_h(Z_1, Z_2) in (<ref>) is symmetric and positive semi-definite. It is well-known (see e.g. Appendix A.5.5 of <cit.>) that in that case, the Schur complement (<ref>) must also be positive semi-definite. Consequently, for Z_i = (Ω_i,yy, Ω_i,yx L^1/2), i=1,2, taking the trace of S_h(Z_1,Z_2) and considering β≥ 1,

    L, P_h^  t  Q_h^-1 P_h^β    =     ((Φ(h Z_1 + (1-h) Z_2)))^β
       ≤    (h (Φ(Z_1)) + (1-h) (Φ(Z_2)))^β
        =     (h  L, Ω_1,yx^  t Ω_1,yy^-1 Ω_1,yx + (1-h)  L, Ω_2,yx^  t Ω_2,yy^-1 Ω_2,yx)^β
       ≤    h  L, Ω_1,yx^  t Ω_1,yy^-1 Ω_1,yx^β + (1-h)  L, Ω_2,yx^  t Ω_2,yy^-1 Ω_2,yx^β

where P_h = h Ω_1,yx + (1-h) Ω_2,yx and Q_h = h Ω_1,yy + (1-h) Ω_2,yy. This convexity inequality concludes the proof. 



 §.§ Theoretical guarantees





  §.§.§ Proof of Theorem <ref>

Let R_n(θ) be the the smooth part of the objective (<ref>),

    R_n(θ)     =     -ln() + ,  + 2 , 
           + ,    + η  L, ^  t  ^β.

For any θ∈Θ and t ∈, by a Taylor expansion,

    R_n( + t  (θ-)) = R_n() + t ∇ R_n(), θ- + e_t(θ, )

for some second-order error term e_t(θ, ). Consider the reparametrization

    ϕ(t) = R_n( + t  (θ-))

so that ϕ^'(0) = ∇ R_n(), θ-. Let =  - and =  -, let also = θ - in a compact form. The estimation error is denoted

    =  -  = ( - ,  - ) = (, ).

Before we start the actual proof, some additional lemmas are needed. They constitute a local study in a sort of r^*-neighborhood of  that we define as

    N_r, α() = {θ∈Θ, ‖‖_F≤ r^* and | []_S̅|_1 ≤α| []_S|_1 }.

Our strategy can be summarized as follows:


  * (Lemma <ref>) Show that there exists a configuration for the regularization parameters (λ, μ, η) so that the estimation error satisfies | []_S̅|_1 ≤α| []_S|_1 for some α > 0.

  * (Lemma <ref>) Find some r^* > 0 and γ_r,η,β,p > 0 such that e_1(θ, ) > γ_r,η,β,p‖‖_F^2 as soon as θ∈ N_r, α().

  * (Lemma <ref>) Exploit this result to show that the estimation error must also satisfy ‖‖_F≤ r^* provided that max{h_a, h_b} is small enough.

  * (Lemma <ref>) Conclude that the theorem holds with high probability, provided that n is large enough.

For the sake of readability, we refer the reader to the Appendix for the numerous constants that are about to appear in the following lemmas and proofs. Thereafter, N_r, α() will always refer to α in (<ref>) and r^* in (<ref>), while the second hypothesis (<ref>) given below is to be assumed with the smallest integer greater than s_α in (<ref>). This is a random hypothesis, which will be controlled with a probability, related to the proximity between the empirical covariance and the true covariance of the predictors, since we recall that S^  (n) has no reason to be an excellent approximation of Σ^* when p ≫ n. This is also assumed by the authors of <cit.>, it is a kind of restricted isometry propertie (RIP), well-known in high-dimensional studies. In particular, we will see through Lemma <ref> that it is satisfied with high probability provided that n is large enough.

    H_2∀  u ≠ 0 such that | u |_0 ≤⌈ s_α⌉, 1/2  u^  t   u ≤  u^  t   u ≤ 3/2  u^  t   u.


    In addition,(  ) ≤ 7/5 (  ).

The next two lemmas give some bounds for expressions that will appear repeatedly.



Under (<ref>) and (<ref>), for all θ∈ N_r, α(), we have the bound

    (   ) ≤

where  is given in (<ref>). In addition,

    (  ) ≥ ()/10 ‖‖_F^2.



Similar reasonings may be found in the proofs of Lem. 1-2 of <cit.>. We simply reworked the constants to make them stick to our study.




Under (<ref>), for all θ∈ N_r, α(), we have the bounds

    (   L ) ≥ and(   L ) ≤

where  and  are given in (<ref>). As a corollary,

    p  ≤  L,    ≤  p .



From Lemmas <ref> and <ref>,

    2 (  L )    ≥    2 ( (  L ) + (  L  +   L ) ) 
       ≥    2 ( (  L ) - ‖  L  +   L ‖_2) 
       ≥    2 ( (  L ) - 2 ‖‖_F ‖ L ‖_2)  ≥ (  L )

as soon as ‖‖_F≤ r^*. From Lemma <ref>, we get

    (   L ) ≥(  L )/()≥(  L )/4 ()

where the inequality in the denominator comes from () ≤() + (), via Lemma <ref>, and the fact that () ≤‖‖_F≤ r^*≤(). For the upper bound, a similar logic gives, with Lemma <ref>,

    √((  L ))   ≤   √((  L )) + √((  L ))
       ≤   √((  L )) + ‖‖_F √((L)) ≤ √(2 (  L ))

for ‖‖_F≤ r^*. It follows from Lemma <ref> that

    (   L ) ≤(  L )/()≤4 (  L )/()

where the inequality in the denominator comes from () ≥() + (), via Lemma <ref>, and the fact that 2 () ≥ -2 ‖‖_F≥ -2  r^*≥ -(). The corollary that concludes the lemma is now immediate.




Assume that λ, μ and η are chosen according to the configuration of the theorem. Then, under (<ref>), the estimation error satisfies

    | [ ]_S̅|_1 ≤α | [  ]_S|_1

where α > 0 is given in (<ref>).


Taking t=1 in the Taylor expansion (<ref>) with θ = and considering the definition of ϕ in (<ref>), by convexity,

    R_n() - R_n() ≥ϕ^'(0).

The first derivative of ϕ will be explicitely computed in (<ref>). For t=0, we find

    ϕ^'(0)     =     -,   + ,  + 2 , 
           +  2 ,    - ,     
           + ηβ  s_L^ β-1[ 2  L,    -  L,     ] 
        =     A_n + ηβ  s_L^ β-1 C_A,  +  B_n + ηβ  s_L^ β-1 C_B,

where s_L is given in (<ref>), where, through the blockwise relations (<ref>), we recognize the random matrices A_n (with max norm h_a) and B_n (with max norm h_b) defined in (<ref>) and (<ref>), and where, coming from the structural regularization term,

    C_A = -     L  and C_B = 2    L.

Whence it follows from the well-known relation |(M_1 M_2) |≤| M_1 |_∞ | M_2 |_1, where M_1 and M_2 are compatible matrices, that

    ϕ^'(0) ≥ -λ/c_λ ||_1 - ηβ  s_L^ β-1ℓ_a ||_1 - μ/c_μ ||_1 - ηβ  s_L^ β-1ℓ_b ||_1

making use of the constants (<ref>), λ≥ c_λ  h_a and μ≥ c_μ  h_b. For the sake of clarity, let

    Δ_n(θ, ) = R_n(θ) + λ ||_1^- + μ ||_1 - R_n() - λ ||_1^- - μ ||_1.

For all θ∈Θ,

    ||_1^- - ||_1^-    =    | [ + ]_S|_1^- + | []_S̅|_1^- - | []_S|_1^-
       ≥   || []_S|_1^- - | []_S|_1^-| + | []_S̅|_1^- - | []_S|_1^-
       ≥   | []_S̅|_1 - | []_S|_1

from the triangle inequality and the fact that, as  is positive definite, the diagonal must belong to S, i.e. (j,j) ∈ S for all 1 ≤ j ≤ q so that any square matrix M of size q is such that [M]_S̅ has diagonal elements all equal to zero. A similar bound obviously holds for ||_1 - ||_1. Now, a straightforward calculation shows that

    Δ_n(, ) ≥c ( | []_S̅|_1 + | []_S̅|_1 ) - c ( | []_S|_1 + | []_S|_1 )

where

    c = max{(c_λ+1) λ/c_λ + ηβ  s_L^ β-1ℓ_a, (c_μ+1) μ/c_μ + ηβ  s_L^ β-1ℓ_b}

and

    c = min{(c_λ-1) λ/c_λ - ηβ  s_L^ β-1ℓ_a, (c_μ-1) μ/c_μ - ηβ  s_L^ β-1ℓ_b}.

Thus, provided that c > 0, which is stated in the configuration of the theorem, it only remains to note that, necessarily,

    Δ_n(, ) ≤ 0

since  is the global minimizer of θ↦ R_n(θ) + λ ||_1^- + μ ||_1. The identification of α given in (<ref>) easily follows.




Under (<ref>) and (<ref>), the second-order error term of (<ref>) satisfies, for t=1 and all θ∈ N_r, α(),

    e_1(θ, ) > γ_r,η,β,p ‖‖_F^2

where γ_r,η,β,p > 0 is given in (<ref>).


From the definition of ϕ in (<ref>) and the fact that ϕ^'(0) = ∇ R_n(), θ-, there exists h ∈  ]0,1[ satisfying

    e_1(θ, ) = 1/2 ϕ^''(h).

To simplify the calculations, let

    u_L =  L,   .

We are going to study the behavior of R_n(, ) in the directions =  + t and =  + t through ϕ(t), where we recall that =  - and =  -. One can see that ϕ(t) moves from R_n(, ) to R_n(, ) as t decreases from 1 to 0. The first derivative is

    ϕ^'(t)     =     -,  + ,  + 2 , 
           +  2 ,    - ,     
           + ηβ  u_L^ β-1[ 2  L,    -  L,     ].

The second derivative is tedious to write but straightforward to establish,

    ϕ^''(t)     =    ,    + 2 [ ,    - 2 ,     
           + ,       ] 
           +  2 ηβ  u_L^ β-1[  L,    - 2  L,     
           +  L,       ] 
           + ηβ (β-1)  u_L^ β-2[ 2  L,    -  L,     ]^2.

First, from the combination of Lemmas <ref> and <ref>, we clearly have u_L≥ 0. We also note that 0 ≤‖2/c M_1 - c M_2 ‖^2_F = 4/c^2 ‖ M_1 ‖^2_F - 4  M_1, M_2  + c^2 ‖ M_2 ‖^2_F for any c ≠ 0 and any matrices M_1 and M_2 of same dimensions. It follows, after some reorganizations, that for any c ≠ 0 and d ≠ 0,

    ϕ^''(t)    ≥   ,   
           +  c_1 ,    + c_2 ,       
           + ηβ  u_L^ β-1[ d_1 ,   L  + d_2  L,       ]

where c_1 = 2 - 4/c^2, c_2 = 2 - c^2, d_1 = 2 - 4/d^2 and d_2 = 2 - d^2. Here we exploited the previous inequality twice, u_L≥ 0 and β≥ 1. From Lemmas <ref>, <ref>, <ref> and <ref>, using sp(M_1  M_2) = sp(M_2  M_1) for square matrices M_1 and M_2, we obtain

    L,       ≤ ,

where  is defined in (<ref>). Replacing L by  and  by , a similar bound obviously holds. Suppose that c and d are chosen so that c_1 > 0, d_1 > 0, c_2 < 0 and d_2 < 0. Then,

    ϕ^''(t)    ≥   ,    [ 1 - | c_2 |  - ηβ  u_L^ β-1| d_2 | ] 
           +  c_1 ,    + ηβ  u_L^ β-1 d_1 ,   L 
       ≥   ,    [ 1 - | c_2 |  - ηβ  (p )^β-1| d_2 | ] 
           +  c_1 ,    + ηβ  (p )^β-1 d_1 ,   L .

Now choose ϵ_S > 0 and ϵ_L > 0 small enough so that ϵ_S  + ηβ  p^ β-1 ^ β ϵ_L < 1 and fix c = √(2 + ϵ_S) and d = √(2 + ϵ_L). We finally obtain

    ϕ^''(t) ≥ a_1 ,    + a_2 ,    + a_3 ,   L

where these positive constants are respectively given by

    a_1 = 1 - ϵ_S  - ηβ  p^ β-1 ^ β ϵ_L,  a_2 = 2 ϵ_S/2+ϵ_Sand a_3 = ηβ  (p )^β-12 ϵ_L/2+ϵ_L.

The combination of Lemmas <ref>, <ref> and <ref> gives, uniformly in t ∈ [0, 1],

    ,   ≥() (  ) ≥‖‖_F^2/4 ^2()

where the inequality in the denominator comes from () ≤ 2 () already established in the proof of Lemma <ref>. Similarly,

    ,   L ≥() (  L ) ≥(L) ‖‖_F^2/2 ().

Lemma <ref> directly enables to bound the last term,

    ,   ≥() (  ) ≥() ‖‖_F^2/20 ().

In conclusion, combining (<ref>), (<ref>) and the upper bounds above,

    e_1(θ, )    ≥   a_1 ‖‖_F^2/8 ^2() + a_2 (L) ‖‖_F^2/4 () + a_3 () ‖‖_F^2/40 ()
       ≥   min{a_1/8 ^2(), a_2 (L)/4 () + a_3 ()/40 ()}‖‖_F^2

and we clearly identify γ_r,η,β,p > 0.




Assume that λ, μ and η are chosen according to the configuration of the theorem. Suppose also that h_a in (<ref>) and h_b in (<ref>) satisfy

    max{h_a, h_b} < r^* γ_r,η,β,p/c_λ, μ√(| S |)

where r^* is given in (<ref>), γ_r,η,β,p in (<ref>) and c_λ, μ in (<ref>). Then, under (<ref>) and (<ref>), the estimation error satisfies ‖‖_F≤ r^*.


By convexity of the objective and optimality of , each move from  in the direction t for t ∈ [0, 1] must lead to a decrease of the objective, i.e.

    R_n( + t ) + λ | + t |_1^- + μ | + t |_1 - R_n() - λ ||_1^- - μ ||_1 ≤ 0.

Taking the notation of (<ref>), this can be rewritten as Δ_n( + t , ) ≤ 0. If ‖‖_F≤ r^* then choose t=1, otherwise calibrate 0 < t < 1 such that ‖ t ‖_F = r^*. Then, from Lemma <ref>, it clearly follows that + t ∈ N_r, α(). Hence, the reasoning preceding (<ref>) still holds and, together with Lemma <ref>, we obtain

    0    ≥   c ( | [t ]_S̅|_1 + | [t ]_S̅|_1 ) - c ( | [t ]_S|_1 + | [t ]_S|_1 ) + γ_r,η,β,p ‖ t ‖_F^2 
       ≥    -c | [t ]_S|_1 + γ_r,η,β,p ‖ t ‖_F^2 
       ≥    -c_λ, μmax{h_a, h_b} √(| S |) ‖ t ‖_F + γ_r,η,β,p ‖ t ‖_F^2

where we used c > 0 and Cauchy-Schwarz inequality to get | [ · ]_S|_1^2 ≤| S | ‖ [ · ]_S‖_F^2. The constant c_λ, μ may be explicitely computed from the configuration of (λ, μ, η) and is given in (<ref>). Note that in the proof of Lemma <ref>, it was sufficient to see that R_n(θ) - R_n() ≥ϕ^'(0) whereas here, we must consider R_n(θ) - R_n() = ϕ^'(0) + e_1(θ, ) to meet our purposes. That explains the presence of γ_r,η,β,p ‖ t ‖_F^2 in the inequality. We deduce that the error must satisfy

    ‖ t ‖_F≤c_λ, μ√(| S |)max{h_a, h_b}/γ_r,η,β,p.

As a corollary, it holds that ‖‖_F > r^* ⇒  c_λ, μ√(| S |)max{h_a, h_b}≥ r^* γ_r,η,β,p or, conversely written, c_λ, μ√(| S |)max{h_a, h_b} < r^* γ_r,η,β,p ⇒ ‖‖_F≤ r^*.




Assume that λ, μ and η are chosen according to the configuration of the theorem. Then, under (<ref>), there exists absolute constants b_1 > 0 and b_2 > 0 such that, for any b_3 ∈  ]0,1[ and as soon as

    n ≥max{ b_1  (q + ⌈ s_α⌉ln(p+q)), ln(10 (p+q)^2) - ln(b_3) },

with probability no less that 1 - ^-b_2 n-b_3 both the random hypothesis (<ref>) is satisfied and the upper bound

    max{h_a, h_b}≤ 16  m^*√(ln(10 (p+q)^2) - ln(b_3)/n)

holds, where h_a and h_b are given in (<ref>) and (<ref>), s_α is defined in (<ref>) and m^* in (<ref>). Hence, one can find a minimal number of observations n_0 such that the theorem holds with high probability as soon as n > n_0.


All the ingredients of the proof are established in <cit.>. The authors start by recalling that there exists absolute constants b_1 > 0 and b_2 > 0 such that hypothesis (<ref>) is satisfied with probability no less than 1 - ^-b_2 n as soon as n ≥ b_1  (q + ⌈ s_α⌉ln(p+q)). We also refer the reader to Lem. 5.1 and Thm. 5.2 of <cit.>, or to Lem. 7.4 of <cit.> for the random bounds of the restricted isometry constants. Afterwards, they prove (see Prop. 4) that, as soon as n ≥ln(10 (p+q)^2) - ln(b_3) for some b_3 > 0, with probability 1-b_3,

    max{h_a, h_b}≤ 16  m^*√(ln(10 (p+q)^2) - ln(b_3)/n).

To find the minimal number of observations, we just need to make sure that the above bound is itself smaller than the one of Lemma <ref>. It is then not hard to see that we may retain the minimal size n_0 given in (<ref>).




§ SOME CONSTANTS


This appendix is entirely dedicated to the constants appearing in the theoretical guarantees. Indeed, a centralization seemed necessary to clarify the rest of the paper, especially the understanding of the main theorem. First, we need to define some constants related to L and to the true values of the model. The bounds

    = (  L )/4 (),  = 4 (  L )/(),  = 4 (  )/().

are useful to control the eigenvalues of some recurrent expressions (Lemmas <ref> and <ref>), uniformly in a neighborhood of = (, ). The true value of the term at the heart of the structural regularization is

    s_L =  L,   .

It plays a role in the proof of Lemma <ref> and, as a consequence, in the definition of the area of validity Λ. This important lemma also  requires to define

    ℓ_a = |   L  |_∞andℓ_b = 2 |   L |_∞

and, in the context of the theorem,

    α = max{(c_λ+1) λ/c_λ + ηβ  s_L^ β-1ℓ_a, (c_μ+1) μ/c_μ + ηβ  s_L^ β-1ℓ_b}/min{(c_λ-1) λ/c_λ - ηβ  s_L^ β-1ℓ_a, (c_μ-1) μ/c_μ - ηβ  s_L^ β-1ℓ_b}.

From α and the cardinality of the true active set | S |, let

    s_α = | S | [ 1 + 12 α^2 ()/()]

which serves as an upper bound in the random hypothesis (<ref>). Similarly, let

    r^* = min{ r_1^*, r_2^*, r_3^*, r_4^*}

where

    r_1^* = ()/2,  r_2^* = √(10) - √(7)/√(5)√((  ))/3 √(3)/2 √(2)√(()),  r_3^* = (  L )/4 ‖ L ‖_2

and

    r_4^* =  (√(2)-1) √((  L ))/√((L)).

Together with α given above, r^* is necessary to build the so-called neighborhood N_r, α() defined in (<ref>), which plays a fundamental role in all our reasonings. It is important to note that, under the configuration of the theorem and hypothesis (<ref>), α > 0 and r^* > 0. Then, Lemma <ref> highlights a new constant, characterizing a strong local convexity of the smooth part of the objective in the neighborhood N_r, α(),

    γ_r,η,β,p = min{a_1/8 ^2(), a_2 (L)/4 () + a_3 ()/40 ()}

where, as it is detailed in the proof of the lemma in question,

    a_1 = 1 - ϵ_S  - ηβ  p^ β-1 ^ β ϵ_L,  a_2 = 2 ϵ_S/2+ϵ_Sand a_3 = ηβ  (p )^β-12 ϵ_L/2+ϵ_L

for some well-chosen ϵ_S > 0 and ϵ_L > 0. Here again, we make sure that γ_r,η,β,p > 0. In the same way, in the context of the theorem,

    c_λ, μ = max{(c_λ+1)  d_λ/c_λ + e_λ, (c_μ+1)  d_μ/c_μ + e_μ}

is needed through Lemma <ref>. Finally, independently of the structure matrix L,

    m^* = |() |_∞ + |(    ) |_∞

is going to play a significative role in the upper bound of the theorem.

*

acm


