Asymptotic Errors for Teacher-Student Convex Generalized Linear Models (or : How to Prove Kabashima's Replica Formula)

Introduction

In the modern era of statistics and machine learning, data analysis often requires solving high-dimensional estimation problems with a very large number of parameters. Developing algorithms for this task and understanding their limitations has become a major challenge. In this paper, we consider this question in the framework of supervised learning under the teacher-student scenario: (i) the data is synthetic and labels are generated by a “teacher”rule and (ii) training is done with a convex Generalized Linear Model (GLM) . Such problems are ubiquitous in machine learning, statistics, communications, and signal processing.

The study of asymptotic (i.e. large-dimensional) reconstruction performance of generalized linear estimation in the teacher-student setting has been the subject of a significant body of work over the past few decades [SST92, WRB93, EVdB01, BM11b, EKBB+13, DM16, ZK16], and is currently witnessing a renewal of interest, especially for the case of identically and independently distributed (i.i.d.) standard normal data matrices, see e.g. [SCC19, HMRT22, MM22]. The aim of this paper is to provide a general analytical formula describing the reconstruction performance of such convex generalized linear models, but for a broader class of more adaptable matrices.

2 Main contributions

We provide a set of equations characterizing the asymptotic statistical properties of the estimator defined by problem (2) with data generated by (1) in the asymptotic setup, for separable, convex losses and penalties (including for instance Logistic, Hinge, LASSO and Elastic net), for rotationally invariant sequences of matrices $\mathbf{F}$. For sufficiently strongly convex problems (in the sense of Lemma 3), our assumptions are classical with respect to earlier work. To extend the result to convex problems however, we require a concentration assumption that we discuss further in section 3.

By doing so, we give, under the aforementioned set of assumptions, a mathematically rigorous proof, of a replica formula obtained heuristically through statistical physics for this problem, notably by Y. Kabashima[Kab08]. This is a significant step beyond the setting of most rigorous work on replica results, which assume matrices to be i.i.d. random Gaussian ones.

Our proof method builds on a detailed mapping between alternating directions descent methods [BPC+11] from convex optimization and a set of algorithms called multi-layer vector approximate message-passing algorithms [MKMZ17, SRF16]. This enables us to use convergence results from convex analysis and dynamical systems to study the trajectories of vector approximate message-passing algorithms.

Beyond the high-dimensional result on the estimator defined by the GLM, our convergence analysis provides a generic condition for the convergence of 2-layer MLVAMP, regardless of the randomness of the design matrix and of the dimensions of the problem, for sufficiently strongly convex problems.

3 Related work

The simplest case of the present question, when both $f$ and $g$ are quadratic functions, can be mapped to a random matrix theory problem and solved rigorously, as in e.g. [HMRT22]. Handling non-linearity is, however, more challenging. A long history of research tackles this difficulty in the high-dimensional limit, especially in the statistical physics literature where this setup is common. The usual analytical approach in statistical physics of learning [SST92, WRB93, EVdB01] is a heuristic, non-rigorous but very adaptable technique called the replica method [MPV87, MM09]. In particular, it has been applied on many variations of the present problem, and laid the foundation of a large number of deep, non-trivial results in machine learning, signal processing and statistics, e.g. [GD89, OKKN90, OK96, Bie03, KWT09, GS10, AG16, Mit19, ESAP+20]. Among them, a generic formula for the present problem has been conjectured by Y. Kabashima, providing sharp asymptotics for the reconstruction performance of the signal ${\bf x_{0}}$ [Kab08].

Proving the validity of a replica prediction is a difficult task altogether. There has been recent progress in the particular case of Gaussian data, where the matrix $\mathbf{F}$ is made of i.i.d. standard Gaussian coefficients. In this case, the asymptotic performance of the LASSO was rigorously derived in [BM11a], and the existence of the logistic estimator discussed in [SCC19]. A set of papers managed to extend this study to a large set of convex losses $g$, using the so-called Gordon comparison theorem [TAH18]. We broaden those results here by proving the Kabashima formula, valid for the set of rotationally invariant matrices introduced above and any convex, separable loss $g$ and sufficiently strongly convex regularizer $f$ under classical conditions. We extend this result to any convex, separable $g$ and $f$ under stronger assumptions.

Our proof strategy is based on the use of approximate-message-passing [DMM09, Ran11], as pioneered in [BM11b], and is similar to a recent work [GAK20] on a simpler setting. This family of algorithms is a statistical physics-inspired variant of belief propagation [Méz89, Kab03, KU04] where local beliefs are approximated by Gaussian distributions. A key feature of these algorithms is the existence of the state evolution equations, a scalar equivalent model which allows to track the asymptotic statistical properties of the iterates at every time step. A series of groundbreaking papers initiated with [BM11a] proved that these equations are exact in the large system limit, and extended the method to treat nonlinear problems [Ran11] and handle rotationally invariant matrices [RSF19, TK22]. We shall use a variant of these algorithms called multi-layer vector approximate message-passing (MLVAMP) [SRF16, FRS18]. The key technical point in our approach is an analysis of the convergence of MLVAMP. This is achieved by phrasing the algorithm as a dynamical system, and then determining sufficient conditions for convergence with linear rate. Our analysis guarantees converging trajectories above a threshold value of the strong convexity parameter of the problem, which is sufficient to complete the proof in that region. We use an analytic continuation to extend the result to convex problems, at the cost of an additional condition discussed after stating our main set of assumption.

Background on MLVAMP

In this section, we present background on the multilayer vector approximate message-passing algorithm developed in [FRS18]. In doing so, we will introduce the key quantities involved in our main theorem. MLVAMP was initially designed as a probabilistic inference algorithm in multilayer architectures. Here, we only focus on the 2-layer version for inference in GLMs, and use the notations of [TK22]. The algorithm can be derived in several ways, notably from expectation-consistent variational inference frameworks such as expectation propagation [Min01], where the target posterior distribution is approximated by a simpler one with moment matching constraints. In the maximum a posteriori setting (MAP), the frequentist optimization framework is recovered, with additional parameter prescriptions due to the probabilistic models, as we will see below. The derivation of the algorithm is, however, not our point of interest. We focus on providing a self-contained interpretation from the convex optimization point of view, in particular in terms of variable splitting.

A common procedure to tackle nonlinear optimization problems involving several functions is variable splitting, so that each non-linearity may be treated independently. Augmenting the Lagrangian with a square penalty on the slack variable equality constraint leads to the family of alternating direction methods of multipliers (ADMM) [BPC+11], where the objective is iteratively minimized in the direction of each initial variable and slack variable. The descent steps then take the form of proximal operators of the non-linearities. For example, on problem (2), adding a slack variable $\mathbf{z}=\mathbf{F}\mathbf{x}$ would lead to the augmented Lagrangian:

where $\alpha>0$ is a free parameter that can enforce strong convexity of the objective if large enough and $\theta$ is a Lagrange multiplier. Updating $\mathbf{x}$ from an update on $\mathbf{z}$ amounts to a linear estimation problem, which can be solved by least squares. This is implemented, for example, in linearized ADMM [BPC+11], where the proximal descent steps are coupled to least-square ones. MLVAMP solves problem (2) by introducing the same splitting as in (3) with an additional trivial splitting for each variable: $\mathbf{x}_{1},\mathbf{x}_{2},\mathbf{z}_{1},\mathbf{z}_{2}$ such that $\mathbf{x}_{1}=\mathbf{x}_{2},\thickspace\mathbf{z}_{1}=\mathbf{F}\mathbf{x}_{1},\thickspace\mathbf{z}_{2}=\mathbf{F}\mathbf{x}_{2}$. In the convex optimization framework, parameters like gradient step sizes, or proximal parameters need to be chosen. In the expectation propagation framework, they are prescribed by expectation-consistency constraints, which leads to additional steps in the algorithm. MLVAMP thus consists in four descent steps on $\mathbf{x}_{1},\mathbf{x}_{2},\mathbf{z}_{1},\mathbf{z}_{2}$, and the updates on the parameters of the functions corresponding to those descent steps. This is shown in the MLVAMP iterations (see (1) further), where $\mathbf{x}_{1},\mathbf{z}_{1}$ are updated using the proximal operators of the loss and regularizer, while $\mathbf{z}_{2}$ and $\mathbf{x}_{2}$ are obtained through least-squares. As mentioned above, the parameters of proximal operators (or denoisers in the signal processing literature) and least-squares are set by probabilistic inference rules (here moment-matching of marginal distributions). It is shown in [FSARS16] that, in the MAP setting, these updates amount to adapting the parameters to the local curvature of the cost function.

2 2-layer MLVAMP and its state evolution

The denoising functions $g_{1x}$ and $g_{1z}$ can be written as proximal operators in the MAP setting:

The LMMSE denoisers $g_{2z}$ and $g_{2x}$ in the MAP setting read (see [SRF16]):

where we defined the matrices $\mathbf{M}_{1}^{(t)}=(\hat{Q}_{2z}^{(t)}\mathbf{F}^{T}\mathbf{F}+\hat{Q}_{2x}^{(t)}\mbox{Id})^{-1}$, and $\mathbf{M}_{2}^{(t)}=(\hat{Q}_{2z}^{(t+1)}\mathbf{F}^{T}\mathbf{F}+\hat{Q}_{2x}^{(t)}\mbox{Id})^{-1}$. As mentioned in the previous section, MLVAMP returns at each iteration two sets of estimators $(\mathbf{\hat{x}}_{1}^{(t)},\mathbf{\hat{x}}_{2}^{(t)})$ and $(\mathbf{\hat{z}}_{1}^{(t)},\mathbf{\hat{z}}_{2}^{(t)})$ which respectively aim at reconstructing the minimizer $\mathbf{\hat{x}}$ and $\mathbf{\hat{z}}=\mathbf{F}\mathbf{\hat{x}}$. At the fixed point, we have $\mathbf{\hat{x}}_{1}^{(t)}=\mathbf{\hat{x}}_{2}^{(t)}$ and $\mathbf{\hat{z}}_{1}^{(t)}=\mathbf{\hat{z}}_{2}^{(t)}$, as proven in [PSAR+20]. The intermediate vectors $\mathbf{h}_{1x}^{(t)}$, $\mathbf{h}_{2x}^{(t)}$, $\mathbf{h}_{1z}^{(t)}$ and $\mathbf{h}_{2z}^{(t)}$ have the key feature that they behave asymptotically as Gaussian centered around $\mathbf{x}_{0}$ and $\mathbf{z}_{0}=\mathbf{F}\mathbf{x}_{0}$, under the set of assumptions given in appendix E.2. More precisely, at each iteration, they converge empirically with second order moment (PL2) towards Gaussian variables:

where $\xi_{1x}^{(t)},\xi_{2x}^{(t)},\xi_{1z}^{(t)},\xi_{2z}^{(t)}$ are i.i.d standard normal random variables independent of all other quantities. The definition of PL(2) convergence is reminded in Appendix A, and we use the notation $\stackrel{{\scriptstyle PL(2)}}{{=}}$ following [RSF19, FRS18]. We can roughly say that the $\hat{Q},\hat{m},\hat{\chi}$’s parameters characterize the distributions of the $\mathbf{h}$’s. Using the representation (10) in the iterations of MLVAMP results in a scalar recursion that tracks the evolution of the parameters of the aforementioned Gaussian distributions. This recursion provides the so-called state evolution equations. The existence of state evolution equations is the reason why we use 2-layer MLVAMP in our proof. Indeed, they allow the construction of iterate paths that lead to the solution of problem (1), while knowing their statistical properties.

Main result

Our main result characterizes the asymptotic empirical distribution of the estimator $\mathbf{\hat{x}}$ defined in (2) with data generated by (1), and of $\mathbf{\hat{z}}=\mathbf{F\hat{x}}$. We start by stating the necessary assumptions.

the functions $f$ and $g$ are proper, closed, convex and separable functions.

the cost function $g(\mathbf{F}.,\mathbf{y})+f(.)$ is coercive, i.e. $\lim_{\lVert\mathbf{x}\rVert\to\infty}g(\mathbf{F}\mathbf{x},\mathbf{y})+f(\mathbf{x})=+\infty$.

for any $\mathbf{x}\in\mbox{dom}(f)$ and any $\mathbf{x}^{\prime}\in\partial f(\mathbf{x})$, there exists a constant $C$ such that $\lVert\mathbf{x}^{\prime}\rVert_{2}\leqslant C(1+\lVert\mathbf{x}\rVert_{2})$. The same holds for $g$ on its domain.

there exist sequences of real analytic functions $g_{\epsilon},f_{\epsilon}$ such that for any $x$, $\lim_{\epsilon\to 0}g_{\epsilon}(x)=g(x)$, $\lim_{\epsilon\to 0}f_{\epsilon}(x)=f(x)$, and for all $\epsilon>0$, $g^{\prime\prime}_{\epsilon}$ and $f^{\prime\prime}_{\epsilon}$ belong to the Schwartz space.

the empirical distributions of the underlying truth $\mathbf{x}_{0}$, eigenvalues of $\mathbf{F}^{T}\mathbf{F}$, and noise vector $w_{0}$, respectively converge empirically with second order moments, as defined in appendix A, to independent scalar random variables $x_{0},w_{0},\lambda$ with distributions $p_{x_{0}}$, $p_{\lambda}$, $p_{w_{0}}$. We assume that the distribution $p_{\lambda}$ is not all-zero and has compact support.

the solution to the set of fixed point equations (13) exists and is unique, for any convex $g$ and $f$ verifying the assumptions above

finally assume that $M,N\to\infty$ with fixed ratio $\alpha=M/N$.

where $H_{x}=\frac{\hat{m}_{1x}^{*}x_{0}+\sqrt{\hat{\chi}_{1x}^{*}}\xi_{1x}}{\hat{Q}^{*}_{1x}}$, $H_{z}=\frac{\hat{m}_{1z}^{*}z_{0}+\sqrt{\hat{\chi}_{1z}^{*}}\xi_{1z}}{\hat{Q}^{*}_{1z}}$ and expectations are taken with respect to the random variables $x_{0}\sim p_{x_{0}}$, $z_{0}\sim\mathcal{N}(0,\sqrt{\rho_{z}})$, $\xi_{1x},\xi_{1z}\sim\mathcal{N}(0,1)$. The parameters $\hat{Q}_{1x}^{*},\hat{Q}_{1z}^{*},\hat{m}_{1x}^{*},\hat{m}_{1z}^{*},\hat{\chi}_{1x}^{*},\hat{\chi}_{1z}^{*}$ are determined by the fixed point of the system:

where $\chi_{x}=(\hat{Q}_{1x}+\hat{Q}_{2x})^{-1}$, $\chi_{z}=(\hat{Q}_{1z}+\hat{Q}_{2z})^{-1}$, and expectations are taken with respect to the random variables $x_{0}\sim p_{x_{0}}$, $z_{0}\sim\mathcal{N}(0,\sqrt{\rho_{z}})$, $y\sim\varphi(z_{0},\omega_{0})$, $\xi_{1x},\xi_{1z}\sim\mathcal{N}(0,1)$, and eigenvalues $\lambda\sim p_{\lambda}$. $\eta$ is a shorthand for the scalar proximal operator:

The set of fixed point equations from Theorem 1 naturally stems from the "replica-symmetric" free energy commonly used in the statistical physics community [MPV87, MM09]. The free energy depends on a set of parameters, and extremizing it with respect to all parameters, i.e. writing the zero gradient condition for each parameter, provides the set of equations (13). We state this correspondence in the following corollary to Theorem 1 :

The fixed point equations from theorem 1 can equivalently be rewritten as the solution of the extreme value problem (15) defined by the replica free energy from [TK22].

$\beta$ is a parameter that corresponds in the physics approach to an inverse temperature. In the $\beta\to\infty$ limit (the so-called zero temperature limit), the integrals defining $\phi_{x}$ and $\phi_{z}$ concentrate on their extremal value. Note that they are closely related to the Moreau envelopes $\mathcal{M}$ [PB+14, BC+11] of $f$ and $g$, which represent a smoothed form of the objective function with the same minimizers:

As immediate corollaries to Theorem 1, we can determine the asymptotic errors of the GLM and the optimal value of the loss function. To characterize the asymptotic reconstruction errors and angles, we can define the norms of the estimators and their overlaps with the ground-truth vectors as the limits

Under the set of Assumptions 1, the squared norms $m_{x}^{*},m_{z}^{*}$ of estimator $\mathbf{\hat{x}}$ defined by (2) and $\mathbf{\hat{z}}=\mathbf{F\hat{x}}$, and their overlaps $q_{x}^{*},q_{z}^{*}$ with ground-truth vectors are almost surely given by:

with $H_{x}$ and $H_{z}$ defined as in Theorem 1.

With the knowledge of the asymptotic overlap $m_{x}^{*}$, and squared norms $q_{x}^{*}$, $\rho_{x}$, most quantities of interest can be determined. For instance, the quadratic reconstruction error is obtained from its definition as $\mbox{E}=\rho_{x}+q_{x}^{*}-2m_{x}^{*}$, while the angle between the ground-truth vector and the estimator is $\theta=\mbox{arccos}(m_{x}^{*}/(\sqrt{\rho_{x}q_{x}^{*}}))$. One can also evaluate the generalization error for new random Gaussian samples, as advocated in [EVdB01], or compute similar errors for the denoising of ${\bf z_{0}}$.

Numerical results

Obtaining a stable implementation of the fixed point equations can be challenging. We provide simulation details in appendix F along with a link to the script we used to produce the figures. Theoretical predictions (full lines) are compared with numerical experiments (points) conducted using standard convex optimization solvers from [PVG+11]. The comparison with finite size ($N\equiv$ a few hundreds) numerical experiments shows that, despite being asymptotic in nature, the predictions are accurate even at moderate system sizes. All experimental points were done with $N=200$ and averaged one hundred times.

We start with a simple verification of the replica prediction in Figure1, on a classification problem where data is generated as $\mathbf{y}=\operatorname{sign}(\mathbf{Fx_{0}})$. We consider two types of singular value distributions for $\mathbf{F}$ and three types of losses: a square loss, a linear support vector classification (SVC) loss and a logistic loss. Technical details and expressions are given in appendix F. We use ridge regularization with penalty $f=\frac{\lambda_{2}}{2}\lVert\cdot\rVert_{2}^{2}$. We plot the reconstruction angle $\theta$ as a function of the aspect ratio of the problem $\alpha$ in Figure 1. A first plot is done with a Marchenko-Pastur eigenvalue distribution for $\mathbf{F}^{T}\mathbf{F}$ corresponding to $\mathbf{F}$ being i.i.d Gaussian. We then move out of the Gaussian setting and change the eigenvalue distribution for (137), which has a qualitatively similar behaviour: it has bounded support, and includes vanishing singular values at a given value $\alpha=1$ of the aspect ratio. We recover a result close to the i.i.d. Gaussian one, including the error peak for the square loss when $\alpha=1$. In both cases, the SVC and the logistic regression perform similarly. Note that error peaks can also be obtained for the max-margin solution as shown in [GLK+20], using a more elaborate teacher.

2 Sparse logistic regression

We now use the replica prediction to study sparse logistic regression with i.i.d Gaussian and row-orthogonal data, the latter being ubiquitous in signal processing. Row-orthogonal data gives rise to a discrete eigenvalue distribution for $\mathbf{F}^{T}\mathbf{F}$ of zeroes and ones:

and is often found to outperform Gaussian sensing matrices for recovery tasks, see e.g. [KWT09] or [GAK20]. In what follows, we define the sparsity $\rho$ of the ground truth vector as the fraction of non-zero components which are sampled from a standard normal distribution. Labels are generated with $\mathbf{y}=\operatorname{sign}(\mathbf{Fx_{0}})$ as for Figure 1.

2.2 Varying the regularization parameter at constant sparsity

2.3 Comparing case

2.4 Discussion

Sketch of proof of Theorem 1

Our proof follows an approach pioneered in [BM11b] where the LASSO risk for i.i.d. Gaussian matrices is determined. The idea is to build a sequence of iterates that provably converges towards the estimator $\mathbf{\hat{x}}$, while also knowing the statistical properties of those iterates through a set of equations. We must therefore concern ourselves with three fundamental aspects:

construct a sequence of iterates with a rigorous statistical characterization that matches their equations of Theorem 1 at the fixed point,

verify that the sequence’s fixed point corresponds to the estimator $\hat{\mathbf{x}}$,

check that this sequence is provably convergent, otherwise the iterates might drift off on a diverging trajectory, and the fixed point would never be reached. We thus make sure the statistical characterization indeed applies to the point of interest $\hat{\mathbf{x}}$.

The following lemma establishes the link between the state evolution equations and our main theorem.

(Fixed point of 2-layer MLVAMP state evolution equations) The state evolution equations of 2-layer MLVAMP from [FRS18], reminded in appendix E, match the equations of Theorem 1 at their fixed point.

This confirms that 2-layer MLVAMP is a good choice to design the sequences that we seek. We know that the iterates of 2-layer MLVAMP can be characterized by state evolution equations which correspond, at their fixed point, to the equations of Theorem 1 by virtue of Lemma 1. The necessary assumptions for the state evolution equations to hold are verified in appendix E.2. We must now show that the estimator of interest defined by (1) and (2) can be reached using 2-layer MLVAMP. We thus continue with point (ii).

(Fixed point of 2-layer MLVAMP) The fixed point of algorithm (1) matches the optimality condition of the unconstrained convex problem Eq.(2)

(Linear convergence of 2-layer MLVAMP for strongly convex problems) Assume $f$ and $g$ are twice differentiable. Define the constrained problem

the functions $f$ and $g$ are proper, closed, convex and separable functions.

the cost function $g(\mathbf{F}.,\mathbf{y})+f(.)$ is coercive, i.e. $\lim_{\lVert\mathbf{x}\rVert\to\infty}g(\mathbf{F}\mathbf{x},\mathbf{y})+f(\mathbf{x})=+\infty$.

there exists a constant $B_{1}$ such that $\frac{1}{N}\lVert\hat{\mathbf{x}}\rVert^{2}_{2}\leqslant B_{1}$ almost surely as $N\to\infty$.

for any $\mathbf{x}\in\mbox{dom}(f)$ and any $\mathbf{x}^{\prime}\in\partial f(\mathbf{x})$, there exists a constant $C$ such that $\lVert\mathbf{x}^{\prime}\rVert_{2}\leqslant C(1+\lVert\mathbf{x}\rVert_{2})$. The same holds for $g$ on its domain.

the empirical distributions of the underlying truth $\mathbf{x}_{0}$, eigenvalues of $\mathbf{F}^{T}\mathbf{F}$, and noise vector $w_{0}$, respectively converge empirically with second order moments, as defined in appendix A, to independent scalar random variables $x_{0},w_{0},\lambda$ with distributions $p_{x_{0}}$, $p_{\lambda}$, $p_{w_{0}}$. We assume that the distribution $p_{\lambda}$ is not all-zero and has compact support.

the solution to the set of fixed point equations (13) exists and is unique for any convex functions $f,g$ verifying the

finally assume that $M,N\to\infty$ with fixed ratio $\alpha=M/N$.

where the scalars $\hat{Q}_{1x},\hat{Q}_{1z}$ and the random variables $H_{x},H_{z}$ are defined as in Theorem 1.

Convergence analysis of 2-layer MLVAMP

where $\mathbf{I}_{d_{v}},\mathbf{I}_{d_{u}}$ are the identity matrices with dimensions of $\mathbf{v},\mathbf{u}$, i.e. $M$ or $N$ in our case. Any co-coercivity property (verified by proximal operators) can be rewritten in matrix form but yields non block diagonal constraint matrices. We will thus directly use the Lipschitz constants for our proof, as they lead to simpler derivations and suffice to prove the required result. The main theorem from [LRP16], adapted to ADMM in [NLR+15], then establishes a sufficient condition for convergence with a linear matrix inequality, involving the matrices defining the linear recast of the algorithm and the constraints. Let us now detail how this approach can be used on 2-layer MLVAMP.

We start by rewriting 2-layer MLVAMP in a more compact form:

For the linear recast, we then define the variables:

This leads to the following linear dynamical system recast of (6.1)-(34):

where $\otimes$ denotes the Kronecker product. We then use a time dependent form of Theorem 4 from [LRP16] in the appropriate form for 2-layer MLVAMP, as was done in [NLR+15] for ADMM.

If, at each time step, (55) is feasible for some $\mathbf{P}\succ 0$ and $\beta_{1}^{(t)},\beta_{2}^{(t)}\geqslant 0$, then for any initialization $\mathbf{h}^{(0)}$, $\mathbf{h}^{(t)}$ converges to $\mathbf{h^{*}}$, the fixed point of (46)-(48):

where $\kappa(\mathbf{P})$ is the condition number of $\mathbf{P}$ and we defined $\tau^{*}=\sup_{t}\tau_{(t)}$.

where $\mathcal{H}_{f},\mathcal{H}_{g}$ are the Hessian of the loss and regularization functions taken at the fixed point. These bounds are obtained from the definitions of $\chi_{x},\chi_{z}$ in the state evolution equations (or equivalently in Theorem 1), and the fact that the derivative of a proximal operator reads, for a twice differentiable function:

2 Numerical experiments for Lemma 3

Acknowledgments

The authors would like to thank Andrea Montanari, Benjamin Aubin, Yoshiyuki Kabashima and Lenka Zdeborová for discussions. This work is supported by the French Agence Nationale de la Recherche under grant ANR-17-CE23-0023-01 PAIL and ANR-19-P3IA-0001 PRAIRIE. Additional funding is acknowledged from “Chaire de recherche sur les modèles et sciences des données”, Fondation CFM pour la Recherche-ENS.

References

Appendix A Convergence of vector sequences

This section is a brief summary of the framework originally introduced in [BM11a] and used in [FRS18, RSF19]. We review the key definitions and verify that they apply in our setting. We remind the full set of state evolution equations from [FRS18] at (118), when applied to learning a GLM, in appendix E, along with the required assumptions for them to hold in appendix E.2. The main building blocks are the notions of vector sequence and pseudo-Lipschitz function, which allow to define the empirical convergence with p-th order moment. Consider a vector of the form

Note that defining an empirically converging singular value distribution implicitly defines a sequence of matrices $\mathbf{F}(N)$ using the definition of rotational invariance from the introduction. This naturally brings us back to the original definitions from [BM11a]. An important point is that the almost sure convergence of the second condition holds for random vector sequences, such as the ones we consider in the introduction. Note that the noise vector $\mathbf{\omega}_{0}$ must also satisfy these conditions, and naturally does when it is an i.i.d. Gaussian one. We also remind the definition of uniform Lipschitz continuity.

for all $\mathbf{x}_{1},\mathbf{x}_{2}\in\mathcal{X}$ and $A\in U$; and

for all $\mathbf{x}\in\mathcal{X}$ and $A_{1},A_{2}\in U$.

We discuss the required assumptions for the state evolution equations to hold in detail, and why they are verified in our setting, in appendix E.2.

Appendix B Convex analysis and properties of proximal operators

(Strong convexity) A proper closed function is $\sigma$-strongly convex with $\sigma>0$ if ${f-\frac{\sigma}{2}\lVert.\rVert^{2}}$ is convex. If f is differentiable, the definition is equivalent to

(Smoothness for convex functions) A proper closed function $f$ is $\beta$-smooth with $\beta>0$ if $\frac{\beta}{2}\lVert.\rVert^{2}-f$ is convex. If f is differentiable, the definition is equivalent to

An immediate consequence of those definitions is the following second order condition: for twice differentiable functions, $f$ is $\sigma$-strongly convex and $\beta$-smooth if and only if:

for all $\mathbf{x},\mathbf{y}\in\mathcal{X}$.

Proximal operators are 1 co-coercive or equivalently firmly-nonexpansive.

(Remark 4.24 [BC+11]) A mapping $T:\mathcal{X}\to\mathcal{X}$ is $\beta$-cocoercive if and only if $\beta$T is half-averaged. This means that T can be expressed as:

(Resolvent of the sub-differential [BC+11]) The proximal mapping of a convex function $f$ is the resolvent of the sub-differential $\partial f$ of $f$:

The following proposition is due to [GB16], and is useful to determine upper bounds on the Lipschitz constant of update functions involving proximal operators.

(Proposition 2 from [GB16]) Assume that $f$ is $\sigma$-strongly convex and $\beta$-smooth and that $\gamma\in]0,\infty[$. Then $\mbox{Prox}_{\gamma f}-\frac{1}{1+\gamma\beta}\rm{Id}$ is $\frac{1}{\frac{1}{1+\gamma\beta}-\frac{1}{1+\gamma\sigma}}$-cocoercive if $\beta>\sigma$ and 0-Lipschitz if $\beta=\sigma$. If $f$ has no smoothness constant, the same holds by taking $\beta=+\infty$.

For any convex and differentiable function $f$, we have:

For a twice differentiable $f$, applying the chain rule then yields:

where $\mathcal{D}$ is the Jacobian matrix and $\mathcal{H}$ the Hessian. Since f is a convex function, its Hessian is positive semi-definite, and, knowing that $\gamma$ is strictly positive, the matrix (\rm{Id}+\gamma\mathcal{H}_{f}(\mbox{\mbox{Prox}_{\gamma f}})) is invertible. We thus have:

Proximal of ridge regularized functions Since we consider only separable functions, we can work with scalar version of the proximal operators. The scalar proximal of a given function with an added ridge regularization can be written:

where the second equality is true only for differentiable $f$. If $f$ is real analytic, we can apply the analytic inverse function theorem [KP02] and verify analyticity in $\lambda_{2}$ of the proximal.

Finally, we remind a result from [BC+11] describing the limiting behavior of regularized estimators for vanishing regularization.

(Theorem 26.20 from [BC+11]) Let f and h be proper, lower semi-continuous, convex functions defined on $\mathcal{X}$. Suppose that $\operatorname*{arg\,min}f\cap\mbox{dom}(h)\neq\emptyset$ and that $h$ is coercive and strictly convex. Then $h$ admits a unique minimizer $\mathbf{x}_{0}$ over $\operatorname*{arg\,min}f$ and , for every $\epsilon\in]0,1[$, the regularized problem

admits a unique solution $\mathbf{x}_{\epsilon}$. If we assume further that $h$ is uniformly convex on any closed ball of the input space, then $\lim_{\epsilon\to 0}\mathbf{x}_{\epsilon}=\mathbf{x}_{0}$.

Appendix C From replica potentials to Moreau envelopes

We can now match the replica potentials with the Moreau envelope. We start from the definition of said potentials, to which we apply Laplace’s approximation:

This is an unconstraint convex optimization problem, thus its optimality condition is enough to characterize its set of minimizers:

Replacing this in the replica potential and completing the square, we get:

where we used the shorthand $X=\frac{\hat{m}_{1x}x_{0}+\sqrt{\hat{\chi}_{1x}}\xi_{1x}}{\hat{Q}_{1x}}$.

Appendix D Fixed point of multilayer vector approximate message passing

Here we show that the fixed point of 2-layer MLVAMP coincides with the optimality condition of the convex problem 2, proving Lemma 2. Writing the fixed point of the scalar parameters of algorithm (1), we get the following prescriptions on the scalar quantities:

and the following ones on the estimates, as proved in [PSAR+20] section III:

We would like the fixed point of MLVAMP to satisfy the following first-order optimality condition

which characterizes the unique minimizer of the unconstraint convex problem (2). Replacing $\mathbf{h}_{1x}$’s expression inside $\mathbf{h}_{2x}$ reads

and using (92) we get $\mathbf{\hat{x}}_{1}=\mathbf{\hat{x}}_{2}$, and a similar reasoning gives $\mathbf{\hat{z}}_{2}=\mathbf{\hat{z}}_{1}$. From (8) and (9), we clearly find $\mathbf{\hat{z}}_{2}=\mathbf{F\hat{x}_{2}}$. Inverting the proximal operators in (5) and (7) yields

Starting from the MLVAMP equation on $\mathbf{h}_{1x}$, we write

which is equal to the left-hand term in (100). Using this equality, as well as $\mathbf{\hat{z}}_{1}=\mathbf{F\hat{x}_{2}}$ and relations (92) and (94) yields

Hence, the fixed point of MLVAMP satisfies the optimality condition (97) and is indeed the desired estimator: $\mathbf{\hat{x}}_{1}=\mathbf{\hat{x}}_{2}=\mathbf{\hat{x}}$.

Appendix E State evolution equations

This appendix is intended mainly for completeness, to show that the fixed point equations from Theorem 1, stemming from the heuristic state evolution written in [TK22] are indeed made rigorous by the results presented in [FRS18].

The state evolution equations track the evolution of MLVAMP (1) and provide statistical properties of its iterates. They are derived in [TK22] taking the heuristic assumption that $\mathbf{h_{1x}},\mathbf{h_{1z}},\mathbf{h_{2x}},\mathbf{h_{2z}}$ behave as Gaussian estimates, which comes from the physics cavity approach:

where $\stackrel{{\scriptstyle PL2}}{{=}}$ denotes $PL2$ convergence. $\mathbf{U}$ and $\mathbf{V}$ come from the singular value decomposition $\mathbf{F}=\mathbf{UDV}^{T}$ and are Haar-sampled; $\mathbf{\xi}_{1x}^{(t)},\mathbf{\xi}_{2x}^{(t)},\mathbf{\xi}_{1z}^{(t)},\mathbf{\xi}_{2z}^{(t)}$ are normal Gaussian vectors, independent from $\mathbf{x_{0}},\mathbf{z_{0}},\mathbf{V}^{T}\mathbf{x_{0}}$ and $\mathbf{U}^{T}\mathbf{z_{0}}$. Parameters $\hat{Q}_{1x}^{(t)},\hat{Q}_{1z}^{(t)}$, $\hat{Q}_{2x}^{(t)},\hat{Q}_{2z}^{(t)}$ are defined through MLVAMP’s iterations (1); while parameters $\hat{m}_{1x}^{(t)},\hat{m}_{1z}^{(t)},\hat{m}_{2x}^{(t)},\hat{m}_{2z}^{(t)}$ and $\hat{\chi}_{1x}^{(t)},\hat{\chi}_{1z}^{(t)},\hat{\chi}_{2x}^{(t)},\hat{\chi}_{2z}^{(t)}$ are prescribed through SE equations. Other useful variables are the overlaps and squared norms of estimators, for $k\in\{1,2\}$:

Starting from assumptions (107), and following the derivation of [TK22] adapted to the iteration order from (1), the heuristic state evolution equations read:

We are interested in the fixed point of these state evolution equations, where $\chi_{1x}^{(t)}=\chi_{2x}^{(t)}=\chi_{x}$, $q_{1x}^{(t)}=q_{2x}^{(t)}=q_{x}$, $m_{1x}^{(t)}=m_{2x}^{(t)}=m_{x}$, $\chi_{1z}^{(t)}=\chi_{2z}^{(t)}=\chi_{z}$, $q_{1z}^{(t)}=q_{2z}^{(t)}=q_{z}$, and $m_{1z}^{(t)}=m_{2z}^{(t)}=m_{z}$ are achieved. From there we easily recover eq. (13). However, these equations are not rigorous since the starting assumptions are not proven. Therefore, we will turn to a rigorous formalism to consolidate those results.

E.2 Necessary assumptions for the rigorous state evolution equations

Here we remind the main assumptions needed for the rigorous state evolution equations to hold, as they are listed for Theorem 1 of [FRS18], and show they are verified in our setting.

the empirical distributions of the underlying truth $\mathbf{x}_{0}$, eigenvalues of $\mathbf{F}^{T}\mathbf{F}$, and noise vector $w_{0}$, respectively converge with second order moments, as defined in appendix A, to independent scalar random variables $x_{0},w_{0},\lambda$ with distributions $p_{x_{0}}$, $p_{\lambda}$, $p_{w_{0}}$. We assume that the distribution $p_{\lambda}$ is not all-zero and has compact support.

assume that $M,N\to\infty$ with fixed ratio $\alpha=M/N$ independent of $M,N$.

the activation function $\phi(.,\mathbf{w}_{0})$ from Eq.(1) is pseudo-Lipschitz of order 2.

the constants $\left\langle\partial_{\mathbf{h}_{1x}^{(t)}}g_{1x}(\mathbf{h}_{1x}^{(t)},\hat{Q}_{1x}^{(t)})\right\rangle,\left\langle\partial_{\mathbf{h}_{1z}^{(t)}}g_{1z}(\mathbf{h}_{1z}^{(t)},\hat{Q}_{1z}^{(t)})\right\rangle,\left\langle\partial_{\mathbf{h}_{2x}^{(t)}}g_{2x}(\mathbf{h}_{2x}^{(t)},\mathbf{h}_{2z}^{(t+1)},\hat{Q}_{2x}^{(t)},\hat{Q}_{2z}^{(t+1)})\right\rangle$ $\left\langle\partial_{\mathbf{h}_{2x}^{(t)}}g_{2z}(\mathbf{h}_{2x}^{(t)},\mathbf{h}_{2z}^{(t)},\hat{Q}_{2x}^{(t)},\hat{Q}_{2z}^{(t)})\right\rangle$ from algorithm (1) are all in $$.

the component estimation functions $g_{1x}(\mathbf{h}_{1x}^{(t)},\hat{Q}_{1x}^{(t)}),g_{1z}(\mathbf{h}_{1z}^{(t)},\hat{Q}_{1z}^{(t)}),g_{2x}(\mathbf{h}_{2x}^{(t)},\mathbf{h}_{2z}^{(t+1)},\hat{Q}_{2x}^{(t)},\hat{Q}_{2z}^{(t+1)}),$ $g_{2z}(\mathbf{h}_{2x}^{(t)},\mathbf{h}_{2z}^{(t)}\hat{Q}_{2x}^{(t)},\hat{Q}_{2z}^{(t)})$ from algorithm (1) are uniformly Lipschitz continuous, at all time steps $t$, respectively in $\mathbf{h}_{1x}^{(t)}$ at $\hat{Q}_{1x}^{(t)}$, in $\mathbf{h}_{1z}^{(t)}$ at $\hat{Q}_{1z}^{(t)}$, $\mathbf{h}_{2x}^{(t)}$ at $\hat{Q}_{2x}^{(t)}$ and in $\mathbf{h}_{2z}^{(t)}$ at $\hat{Q}_{2z}^{(t)}$.

The first four points are included in the set of assumptions 1 and are therefore verified. We need to check the last two points, starting with the function $g_{1x}(\mathbf{h}_{1x}^{(t)},\hat{Q}_{1x}^{(t)})=\mbox{Prox}_{f/\hat{Q}_{1x}^{(t)}}(\mathbf{h}_{1x}^{(t)})$. Since proximal operators are firmly nonexpansive, they are 1-Lipschitz and we thus have, using the separability of the function $f$:

where the last line is obtained using the scaling conditions on the subdifferential of $f$ from assumption 1. Then, for any $\hat{Q}_{1x}^{(t)},\hat{Q}_{1x}^{(t^{\prime})}$, $\lVert\mbox{Prox}_{f/\hat{Q}_{1x}^{(t)}}-\mbox{Prox}_{f/\hat{Q}_{1x}^{(t^{\prime})}}\rVert_{2}\leqslant C(1+\lVert\mathbf{h}_{1x}^{(t)}\rVert_{2})\lvert\hat{Q}_{1x}^{(t)}-\hat{Q}_{1x}^{(t^{\prime})}\rvert$ and $g_{1x}(\mathbf{h}_{1x}^{(t)},\hat{Q}_{1x}^{(t)})$ is uniformly Lipschitz in $\mathbf{h}_{1x}^{(t)}$ at $\hat{Q}_{1x}^{(t)}$, at any time index $t$. The argument is identical for $g_{1z}(\mathbf{h}_{1z}^{(t)},\hat{Q}_{1z}^{(t)})=\mbox{Prox}_{f/\hat{Q}_{1z}^{(t)}}(\mathbf{h}_{1z}^{(t)})$. The functions $g_{2x}(\mathbf{h}_{2x}^{(t)},\mathbf{h}_{2z}^{(t+1)},\hat{Q}_{2x}^{(t)},\hat{Q}_{2z}^{(t+1)}),g_{2z}(\mathbf{h}_{2x}^{(t)},\mathbf{h}_{2z}^{(t)},\hat{Q}_{2x}^{(t)},\hat{Q}_{2z}^{(t)})$ have explicit expressions and it is straightforward to check the last two points using linear algebra and the assumptions on the spectrum of $\mathbf{F}^{\top}\mathbf{F}$.

E.3 Rigorous state evolution formalism

We now look into the state evolution equations derived for MLVAMP in [SRF16]. Those equations are proven to be exact in the asymptotic limit, and follow the same algorithm as (1). In particular, they provide statistical properties of vectors $\mathbf{h}_{1x},\mathbf{h}_{2x},\mathbf{h}_{1z},\mathbf{h}_{2z}$. We can read relations from [FRS18] using the following dictionary between our notations and theirs, valid at each iteration of the algorithm:

Placing ourselves in the asymptotic limit, [FRS18] shows the following equalities:

where $\mathbf{Q}_{0}^{-}\sim\mathcal{N}(0,\tau_{0}^{-})^{N}$ and $\mathbf{Q}_{1}^{-}\sim\mathcal{N}(0,\tau_{1}^{-})^{N}$ are i.i.d. Gaussian vectors. $\mathbf{Q}_{0}^{+}$, $\mathbf{Q}_{1}^{+}$ have the following norms and non-zero correlations with ground-truth vectors $\mathbf{Q}_{0}^{0},\mathbf{Q}_{1}^{0}$:

With simple manipulations, we can rewrite (112) as:

Simple bookkeeping transforms equations (115) into a rigorous statement of starting assumptions (112) from [TK22]. Since those assumptions are now rigorously established in the asymptotic limit, the remaining derivation of state evolution equations (108) holds and provides a mathematically exact statement.

E.4 Scalar equivalent model of state evolution

By virtue of Lemma 5 from [RSF19], the six previous vectors have elements that converge empirically to a Gaussian variable. Hence, all defined vectors have an element-wise separable distribution, and we can write the state evolution as a scalar model on random variables sampled from those distributions. To do so, we will simply write the variables without the bold font: for instance $Z_{0}^{0}\sim p_{x_{0}}$, $s_{\nu}\sim p_{\nu}$, and $Q_{0}^{-}$ refers to the random variable distributed according to the element-wise distribution of vector $\mathbf{Q}_{0}^{-}$. The scalar random variable state evolution from [FRS18] now reads:

E.5 Direct matching of the state evolution fixed point equations

To be consistent, we should be able to show that equations (118) allow us to recover equations (108) at their fixed point. Although somewhat tedious, this task is facilitated using dictionaries (111) and (117). We shall give here an overview of this matching through a few examples.

Let us start from the rigorous scalar state evolution, in particular equation (118h) that defines variable $Q_{0}^{+}$. We get rid of time indices here since we focus on the fixed point. We first compute the correlation

hence we perfectly recover equations (108e) at the fixed point.

We start again from (118h) and square it:

Since $Q_{0}^{-}$ is a Gaussian variable, independent from $Q_{0}^{0}$, we can use Stein’s lemma and use equation (118f) to get

Replacing (126) and (128) into (125), we reach

and $s_{\nu}^{2}\sim p_{\lambda}$. Applying this to $h(s)=\dfrac{\gamma_{1}^{-}s^{2}}{\gamma_{1}^{-}s^{2}+\gamma_{0}^{+}}$ and starting from (118m), we rewrite

with $\lambda\sim p_{\lambda}$, which translates into equation (108t):

In a similar fashion, we can recover all equations (108) by writing variances and correlations between scalar random variables defined in (118), and using the independence properties established in [FRS18]; thus directly showing the matching between the two state evolution formalisms at their fixed point.

Appendix F Numerical implementation details

The plots were generated using the toolbox available at https://github.com/cgerbelo/Replica_GLM_orth.inv.git Here we give a few derivation details for implementation of the equations presented in Theorem 1. We provide the Python script used to produce the figures in the main body of the paper as an example. The experimental points were obtained using the convex optimization tools of [PVG+11], with a data matrix of dimension $N=200,M=\alpha N$, for $\alpha\in[0.1,3]$. Each point is averaged 100 times to get smoother curves. The theoretical prediction was simply obtained by iterating the equations from Theorem 1. This can lead to unstable numerical schemes, and we include a few comments about stability in the code provided with this version of the paper. For Gaussian data, the design matrices were simply obtained by sampling a normal distribution $\mathcal{N}(0,\sqrt{1/M})$, effectively yielding the Marchenko-Pastur distribution [TV04] for averaging on the eigenvalues of $\mathbf{F}^{T}\mathbf{F}$ in the state evolution equations :

where $a=\sqrt{1-\left(\frac{1}{\alpha}\right)^{2}},b=\sqrt{1+\left(\frac{1}{\alpha}\right)^{2}}$, and $(0,x)^{+}=\max(0,x)$. For the example of orthogonally invariant matrix with arbitrary spectrum, we chose to sample the singular values of $\mathbf{F}$ from the uniform distribution $\mathcal{U}(\left[(1-\alpha)^{2},(1+\alpha)^{2}\right])$. This leads to the following distribution for the eigenvalues of $\mathbf{F}^{T}\mathbf{F}$:

For the elastic net regularization, we can obtain an exact expression, avoiding any numerical integration. The proximal of the elastic net reads:

where $s\left(.,\frac{\lambda_{1}}{\hat{Q}_{1x}}\right)$ is the soft-thresholding function:

We assume that the ground-truth $x_{0}$ is pulled from a Gauss-Bernoulli law of the form:

Note that we did our plots with $\rho=1$, but this form can be used to study the effect of sparsity in the model. Writing $X=\frac{\hat{m}_{1x}x_{0}+\sqrt{\hat{\chi}_{1x}}\xi_{1x}}{\hat{Q}_{1x}}$, and remembering that $\xi_{1x}\sim\mathcal{N}(0,1)$, some calculus then shows that:

F.2 Loss functions

The loss functions sometimes have no closed form, as is the case for the logistic loss. In that case, numerical integration cannot be avoided, and we recommend marginalizing all the possible variables that can be averaged out. In the present model, if the teacher $y$ is chosen as a sign, one-dimensional integrals can be reached, leading to stable and reasonably fast implementation (a few minutes to generate a curve comparable to those of Figure 1 for the non-linear models, the ridge regression being very fast). The interested reader can find the corresponding marginalized prefactors in the code jointly provided with this paper.

its proximal and partial derivative then read:

Assuming $y\in\{-1,+1\}$, its proximal and partial derivative then read:

Its proximal (at point p) is the solution to the fixed point problem:

and its derivative, given that the logistic loss is twice differentiable, reads:

Appendix G Proof of Lemma 3: Convergence analysis of 2-layer MLVAMP

In this section, we give the detail of the convergence proof of 2-layer MLVAMP.

This proof is quite straightforward and close to the one of Theorem 4 from [LRP16]. Multiplying Eq.(55) on the left and right by $[(\mathbf{h}^{(t)}-\mathbf{h}^{(t-1)})^{\top}\quad(\mathbf{u}^{(t)}-\mathbf{u}^{(t-1)})^{\top}]$ and its transpose respectively, we get

Using the definition of the iteration (46)-(48), this simplifies to

Letting $\tau^{*}=\sup_{t}\tau_{(t)}$, an immediate induction concludes the proof.

The matrices on the r.h.s. of the previous equation are all diagonalisable in the same basis. Then each eigenvalue has the form

G.3 Operator norms and Lipschitz constants

The norms of the linear operators $\mathbf{W_{1}}^{(t)},\mathbf{W_{2}}^{(t)},\mathbf{W_{3}}^{(t)},\mathbf{W_{4}}^{(t)}$ can be computed or bounded with respect to the singular values of the matrix $\mathbf{F}$. The derivations are straightforward and do not require any specific mathematical result. Denoting $\lVert\mathbf{W}\rVert$ the operator norm of a given matrix $\mathbf{W}$, we have the following:

Proposition 3 gives the following expression:

In this case, we have from Proposition 3:

The upper bound on the Lipschitz constant is therefore:

This setting is not necessary for our proof, because we only handle penalty functions which have a strictly positive strong convexity constant, by adding a ridge term. However, we list it for completeness. In this case, the only information we have is the firm nonexpansiveness of the proximal operator, which leads us to the same derivation as the previous one up to (188), where the first term in the sum can be positive or negative. This yields the Lipschitz constant:

G.4 Dynamical system convergence analysis

where all the matrices constituting the blocks have been defined in section 6. This gives the following form for the constraint matrix:

We are then trying to find the conditions for the following problem to be feasible with $0<\tau<1$:

Schur’s lemma then gives that the strict version of (210), which we will consider, is equivalent [HJ12] to:

We start with $-(\mathbf{B}^{T}\mathbf{P}\mathbf{B}+\mathbf{H}_{3})$.

where $\mathbf{K}_{1}=(\beta_{1}\mathbf{I}_{N\times N}-\beta_{2}\omega_{2}^{2}\mathbf{W}_{4}^{T}\mathbf{W}_{4}-\mathbf{W}_{1}^{T}\mathbf{P}_{2}\mathbf{W}_{1})^{-1}$. We start with (213). A sufficient condition for it to hold true is:

Using the bounds from appendix G.3, we have:

where $b_{1},b_{2}$ are constants independent of $\lambda_{2},\hat{Q}_{2x},\hat{Q}_{1z}$. We now turn to (214). A sufficient condition for it to hold is:

Note that condition (213) ensures that the denominator in (G.4.1) is non-zero. We then have:

This quantity can be bounded above by a constant independent of $\lambda_{2},\hat{Q}_{2x},\hat{Q}_{1z}$ for arbitrarily large $\hat{Q}_{2x}$. Let $b_{3}$ be such a constant . Then a sufficient condition for condition (214) to hold is:

we see that $\beta_{1}$ must scale linearly with $\hat{Q}_{1z}$ which is one of the parameters that is made arbitrarily large. Then $\beta_{1}$ also needs to become arbitrarily large for the conditions to hold. We choose $\beta_{1}=2\beta_{2}\hat{Q}_{1z}b_{2}+\lambda_{max}(\mathbf{P}_{2})b_{1}$ for the rest of the proof. Condition (220) is then verified, and $\beta_{2}$ needs to be chosen according to condition (G.4.1), which becomes:

This has a bounded solution for large values of $\hat{Q}_{1z}$. We now turn to the second part of (211).

We thus only need to characterize the lower right block of $\mathbf{Z}$. It is easy to see that conditions (213) and (214) also enforce that both the Schur complements associated with the upper left and lower right blocks of $-(\mathbf{B}^{T}\mathbf{P}\mathbf{B}+\mathbf{H}_{3})$ are invertible, thus giving the following form for $\mathbf{Z}_{3}$ using the block matrix inversion lemma [HJ12]:

where $\mathbf{K}_{2}=(\beta_{1}\mathbf{I}_{M}-\mathbf{P}_{1}-\mathbf{W}_{2}^{T}\mathbf{P}_{2}\mathbf{W}_{2})^{-1}$. We thus have the following upper bound on the largest eigenvalue of $\mathbf{Z}_{3}$:

where $k=\frac{b_{1}b_{3}\lambda_{max}^{2}(\mathbf{P}_{2})}{\beta_{2}-\lambda_{max}(\mathbf{P}_{1})-b_{2}\lambda_{max}(\mathbf{P}_{2})}$. Using the prescription $\beta_{1}=2\beta_{2}\hat{Q}_{1z}b_{2}+\lambda_{max}(\mathbf{P}_{1})b_{1}$, we get:

where $b_{4}$ is a constant independent of the arbitrarily large parameters $\lambda_{2},\hat{Q}_{2x},\hat{Q}_{1z}$. Thus $\lambda_{max}(\mathbf{Z}_{3})$ can be made arbitrarily small by making $\lambda_{2}$ arbitrarily large. We now want to find conditions for $\tau^{2}\mathbf{P}-\mathbf{H}_{1}+\mathbf{H}_{2}(\mathbf{B}^{T}\mathbf{P}\mathbf{B}+\mathbf{H}_{3})^{-1}\mathbf{H}_{2}^{T}\succ 0$ which is equivalent to:

We start with the upper matrix inequality, for which a sufficient condition is:

Using the bounds from appendix G.3, we have:

Thus there exists a constant $b_{5}$, independent of $\lambda_{2},\hat{Q}_{1z},\hat{Q}_{2x}$ such that, for sufficiently large $\hat{Q}_{1z}$:

which gives the following sufficient condition for the upper left block in (G.4.2):

A sufficient condition for the lower right block in (G.4.2) then reads:

We remind the reader that $\hat{Q}_{1z},\hat{Q}_{2x}$ grow linearly with $\lambda_{2}$. Thus the dominant scaling at large $\lambda_{2}$ is (exchanging $\hat{Q}_{2x}$ with $\hat{Q}_{1z}$ up to a constant):

where $b_{6}$ is a constant independent of the arbitrarily large quantities. The final condition becomes:

Appendix H Analytic continuation

where $H_{x}=\frac{\hat{m}_{1x}^{*}x_{0}+\sqrt{\hat{\chi}_{1x}^{*}}\xi_{1x}}{\hat{Q}_{1x}}$ is defined in Theorem 1. We would like to show that this equality still holds for any $\lambda_{2}>0$. To do so we will show that, for a real analytic approximation of problem Eq.(2), both sides of Eq.(247) are real analytic in $\lambda_{2}$. We may then use the real analytic continuation theorem, as given in [KP02] to extend to any $\lambda_{2}>0$. We will treat the case $\lambda_{2}=0$ separately. In what follows, we will write the dependency in $\lambda_{2}$ of the estimator explicitly, i.e., $\hat{\mathbf{x}}=\hat{\mathbf{x}}(\lambda_{2})$.

We remind a useful characterization of real analytic functions from [KP02]:

Let $f\in\mathcal{C}^{\infty}(I)$ for some open interval I. The function f is in fact real analytic on I if and only if, for each $\alpha\in I$, there are an open interval J, with $\alpha\in J\subset I$, and finite constants $C>0$ and $R>0$ such that the derivatives of f satisfy :

We also remind the formula for the higher order derivatives of a composition of two infinitely differentiable functions:

where $k=k_{1}+k_{2}+...+k_{n}$ and the sum is taken over all $k_{1},k_{2},...,k_{n}$ for which $k_{1}+2k_{2}+...+nk_{n}=n$.

The following lemma establishes bounds on the higher order derivatives of $\hat{\mathbf{x}}(\lambda_{2})$ with respect to $\lambda_{2}$.

$\hat{\mathbf{x}}(\lambda_{2})$ is infinitely differentiable w.r.t. $\lambda_{2}$ and, for any integer $p$, there exists a constant $K^{\prime}$ such that its elementwise p-th derivative, denoted $D_{\lambda_{2}}^{(p)}\hat{\mathbf{x}}(\lambda_{2})$ verifies, almost surely

Furthermore, $D^{(p)}_{\lambda_{2}}\hat{\mathbf{x}}(\lambda_{2})$ is a Lipschitz function of $\hat{\mathbf{x}}(\lambda_{2})$.

Recall the strongly convex problem, for any finite N,

Owing to assumption 1, we have almost surely

and the identity is a Lipshchitz function of $\hat{\mathbf{x}}(\lambda_{2})$ The function of $\lambda_{2}$ defined by :

is always zero valued from the definition of $\mathbf{\hat{x}}(\lambda_{2})$, thus all its derivatives are zero. Taking the first derivative with respect to $\lambda_{2}$ yields:

where $D^{p}$ is the $(N\times 1)$ dimensional element-wise p-th differential of $\hat{\mathbf{x}}(\lambda_{2})$. We then define the operator

We obtain a simple expression for $D\mathbf{\hat{x}}(\lambda_{2})$

Since $f$ and $g$ are convex, the operator norm of $\mathcal{O}^{-1}(\lambda_{2})$ is bounded with probability one, and $D\mathbf{\hat{x}}(\lambda_{2})$ is a Lipschitz function of $\mathbf{\hat{x}}(\lambda_{2})$ where $\frac{1}{N}\lVert D\mathbf{\hat{x}}(\lambda_{2})\rVert_{2}^{2}$ is almost surely bounded.

Assume the property is verified up to $p-1$. For higher order derivatives, applying Leibniz’s rule on Eq.(H.1) gives, denoting $\mathcal{O}^{(i)}(\lambda_{2})$ the i-th derivative of $\mathcal{O}(\lambda_{2})$, for the (p-1)-th derivative of (H.1) :

We obtain the recursion on the differentials of $\mathbf{\hat{x}}(\lambda_{2})$ :

Under assumption 1, the function $\psi(\lambda_{2})$ defined as

so there exists a constant $C_{\psi}^{\prime}$ such that

which is almost surely bounded. We have also proved in the previous lemma that $D\hat{\mathbf{x}}(\lambda_{2})$ is a Lipschitz function of $\lambda_{2}$, thus $D\psi(\lambda_{2})$ is a PL2 function of $\hat{\mathbf{x}}(\lambda_{2})$ and its limit exists according to Assumption 1 (c). For the higher order derivatives, we use proposition 6 to obtain, for any coordinate $1\leqslant i\leqslant n$ :

The assumption on the higher order derivatives of $\phi$ from Theorem 1 and Lemma 8 implies that the term $\phi^{(k)}\left(\hat{x}_{i}(\lambda_{2})\right)\left(\frac{\hat{x}_{i}^{(1)}(\lambda_{2})}{1!}\right)^{k_{1}}\left(\frac{\hat{x}_{i}^{(2)}(\lambda_{2})}{2!}\right)^{k_{2}}...\left(\frac{\hat{x}_{i}^{(p)}(\lambda_{2})}{p!}\right)^{k_{p}}$ has bounded absolute value with probability one, for all coordinates $i$. Using the characterization of real analytic functions and assumption 1 (c) from proposition 5, this concludes the proof. ∎

H.3 Real analytic approximation of strongly convex problems

where $g_{\epsilon},f_{\epsilon}$ are real analytic approximations of the loss $g$ and regularizer $f$ verifying assumption 1(e). To relax the analytic approximation, we need to prove the following equality.

Under assumption 1 (c) and owing to the definition of PL2 functions, it is sufficient to prove

Since minimizers of convex functions are fixed points of the corresponding proximity operators, it holds that

The results from appendix G.3.2 show that proximity operators of strongly convex functions are contractions, thus their exists a positive constant $L_{\lambda_{2}}<1$ such that for any realisation of $\mathbf{F},\mathbf{x}^{0},\boldsymbol{\omega}_{0}$

Furthermore, the function $\mbox{prox}_{\mathcal{C}_{\epsilon}(.)+\frac{\lambda_{2}}{2}\lVert.\rVert_{2}^{2}}(.)$ converges uniformly to $\mbox{prox}_{\mathcal{C}(.)+\frac{\lambda_{2}}{2}\lVert.\rVert_{2}^{2}}(.)$ when $\epsilon\to 0$, and thus

for functions $f,g$ that may not be strictly convex. To do so we will use Theorem 26.20 from [BC+11], which is reminded in appendix B, proposition 4. Under assumption 1 and since the $l_{2}$ norm is strongly convex thus uniformly convex, we have, denoting $\hat{\mathbf{x}}_{0}$ the unique least $l_{2}$ norm element in $\operatorname*{arg\,min}_{\mathbf{x}\in\mathcal{X}}g(\mathbf{F}\mathbf{x},\mathbf{y})+f(\mathbf{x})$,

We can therefore uniquely define the continuous extension of any continuous observable $\phi$ of $\hat{\mathbf{x}}(\lambda_{2})$ such that $\phi(\lambda_{2}=0)=\phi(\hat{\mathbf{x}}_{0})$. Then this observable and the corresponding function implicitly defined by the set of fixed point equations are continuous on $[0,+\infty[$ and equal for any $\lambda_{2}\in]0,+\infty[$, and thus also equal at $\lambda_{2}=0$ using the definition of continuity and the fact that $]0,+\infty[$ is dense in $[0,+\infty[$.

H.6 Real analytic approximation of usual cost functions with fast decaying higher-order derivatives

thus $\mathcal{W}_{\epsilon}\left[\lvert.\rvert)\right]$ is strongly convex and its higher order derivatives all decay faster than any finite order polynomial. A similar computation shows that, for the hinge loss,