The Impact of Regularization on High-dimensional Logistic Regression

Introduction

Logistic regression is the most commonly used statistical model for predicting dichotomous outcomes . It has been extensively employed in many areas of engineering and applied sciences, such as in the medical and social sciences . As an example, in medical studies logistic regression can be used to predict the risk of developing a certain disease (e.g. diabetes) based on a set of observed characteristics from the patient (age, gender, weight, etc.) Linear regression is a very useful tool for predicting a quantitive response. However, in many situations the response variable is qualitative (or categorical) and linear regression is no longer appropriate . This is mainly due to the fact that least-squares often succeeds under the assumption that the error components are independent with normal distribution. In categorical predictions, however, the error components are neither inependent nor normally distributed . In logistic regression we model the probability that the label, $Y$, belongs to a certain category. When no prior knowledge is available regarding the structure of the parameters, maximum likelihood is often used for fitting the model. Maximum likelihood estimation (MLE) is a special case of maximum a posteriori estimation (MAP) that assumes a uniform prior distribution on the parameters. In many applications in statistics, machine learning, signal processing, etc., the underlying parameter obeys some sort of structure (sparse, group-sparse, low-rank, finite-alphabet, etc.). For instance, in modern applications where the number of features far exceeds the number of observations, one typically enforces the solution to contain only a few non-zero entries. To exploit such structural information, inspired by the Lasso algorithm for linear models, researchers have studied regularization methods for generalized linear models . From a statistical viewpoint, adding a regularization term provides a MAP estimate with a non-uniform prior distribution that has higher densities in the set of structured solutions.

Classical results in logistic regression mainly concern the regime where the sample size, $n$, is overwhelmingly larger than the feature dimension $p$. It can be shown that in the limit of large samples when $p$ is fixed and $n\rightarrow\infty$, the maximum likelihood estimator provides an efficient estimate of the underlying parameter, i.e., an unbiased estimate with covariance matrix approaching the inverse of the the Fisher information . However, in most modern applications in data science, the datasets often have a huge number of features, and therefore, the assumption $\frac{n}{p}\gg 1$ is not valid. Sur and Candes have recently studied the performance of the maximum likeliood estimator for logistic regression in the regime where $n$ is proportional to $p$. Their findings challenge the conventional wisdom, as they have shown that in the linear asymptotic regime the maximum likelikehood estimate is not even unbiased. Their analysis provides the precise performance of the maximum likelihood estimator. There have been many studies in the literature on the performance of regularized (penalized) logistic regression, where a regularizer is added to the negative log-likelihood function (a partial list includes ). These studies often require the underlying parameter to be heavily structured. For example, if the parameters are sparse the sparsity is taken to be $o(p)$. Furthermore, they provide orderwise bounds on the performance but do not give a precise characterization of the quality of the resulting estimate. A major advantage of adding a regularization term is that it allows for recovery of the parameter vector even in regimes where the maximum likelihood estimate does not exist (due to an insufficient number of observations.)

2 Summary of contributions

Preliminaries

and the proximal operator is the solution to this optimization, i.e.,

2 Mathematical Setup

where $\rho^{\prime}(t):=\frac{e^{t}}{1+e^{t}}$ is the standard logistic function. The goal is to compute an estimate for $\bm{\beta}^{*}$ from the training data $\mathcal{D}$. The maximum likelihood estimator, $\hat{\bm{\beta}}_{ML}$, is defined as,

Main Results

In this section, we present the main result of the paper, that is the characterization of the asymptotic performance of regularized logistic regression (RLR). When the estimation performance is measured via a locally- Lipschitz function (e.g. mean-squared error), Theorem 1 precisely predicts the asymptotic behavior of the error. The derived expression captures the role of the regularizer, $f(\cdot)$, and the particular distribution of $\bm{\beta}^{*}$, through a set of scalars derived by solving a system of nonlinear equations. In Section 3.1 we present this system of nonlinear equations along with some insights on how to numerically compute its solution. After formally stating our result in Section 3.2, we use that to predict the general behavior of $\hat{\bm{\beta}}$. In particular, in Section 3.3 we compute its correlation with the true signal as well as its mean-squared error.

As we will see in Theorem 1, given the signal strength $\kappa$, and the ratio $\delta$, the asymptotic performance of RLR is characterized by the solution to the following system of nonlinear equations with six unknowns $(\alpha,\sigma,\gamma,\theta,\tau,r)$.

Here $Z,Z_{1},Z_{2}$ are standard normal variables, and $\beta\sim\Pi$, where $\Pi$ denotes the distribution on the entries of $\bm{\beta}^{*}$. The following remarks provide some insights on solving the nonlinear system.

2 Asymptotic performance of regularized logistic regression

We are now able to present our main result. Theorem 1 below describes the average behavior of the entries of $\hat{\bm{\beta}}$, the solution of the RLR. The derived expression is in terms of the solution of the nonlinear system (6), denoted by $(\bar{\alpha},\bar{\sigma},\bar{\gamma},\bar{\theta},\bar{\tau},\bar{r})$. An informal statement of our result is that as $n\rightarrow\infty$, the entries of $\hat{\bm{\beta}}$ converge as follows,

In other words, the RLR solution has the same behavior as applying the proximal operator on the "perturbed signal", i.e., the true signal added with a Gaussian noise.

where $Z\sim\mathcal{N}(0,1)$, $\beta\sim\Pi$ is independent of $Z$, and the function $\Gamma(\cdot,\cdot)$ is defined in (8).

We defer the detailed proof to the Appendix. In short, to show this result we first represnt the optimization as a bilinear form $\mathbf{u}^{T}\mathbf{X}\mathbf{v}$, where $\mathbf{X}$ is the measurement matrix. Applying the CGMT to derive an equivalent optimization, we then simplify this optimization to obtain an unconstrained optimization with six scalar variables. The nonlinear system (6) represents the first-order optimality condition of the resulting scalar optimization. Before stating the consequences of this result, a few remarks are in order.

The assumptions in Theorem 1 are chosen in a conservative manner. In particular, we could relax the separability condition on $f(\cdot)$, to some milder condition in terms of asymptotic convergence of its proximal operator. Furthermore, one can relax the assumption on the entries of $\bm{\beta}^{*}$ being i.i.d. to a weaker assumption on the empirical distribution of its entries. However, for the applications of this paper, the theorem in its current form is adequate.

The performance measure in Theorem 1 is computed in terms of evaluation of a locally-Lipschitz function, $\Psi(\cdot,\cdot)$ . As an example, $\Psi(u,v)=(u-v)^{2}$ can be used to compute the mean-squared error. Later on, we will appeal to this theorem with various choices of $\Psi$ to evaluate different performance measures on $\hat{\bm{\beta}}$.

3 Correlation and variance of the RLR estimate

As the first application of Theorem 1 we compute common descriptive statistics of the estimate $\hat{\bm{\beta}}$. In the following corollaries, we establish that the parametrs $\bar{\alpha}$, and $\bar{\sigma}$ in (6) correspond to the correlation and the mean-squared error of the resulting estimate.

As $p\rightarrow\infty$, $\frac{1}{||\bm{\beta}^{*}||^{2}}~{}\hat{\bm{\beta}}^{T}\bm{\beta}^{*}\overset{P}{\longrightarrow}\bar{\alpha}$ .

Recall that $||\bm{\beta}^{*}||^{2}=p\kappa^{2}$. Applying Theorem 1 with $\Psi(u,v)=uv$ gives,

where the last equality is derived from the first equation in the nonlinear system (6), along with the fact that $(\bar{\alpha},\bar{\sigma},\bar{\gamma},\bar{\theta},\bar{\tau},\bar{r})$ is a solution to this system. ∎

We appeal to Theorem 1 with $\Psi(u,v)=(u-\bar{\alpha}v)^{2}$,

where the last equality is derived from the third equation in the nonlinear system (6) together with the result of Corollary 1. ∎

This indicates that the proximal operator in this case is just a simple rescaling. Substituting (13) in the nonlinear system (6), we can rewrite the first three equations as follows,

Consider the optimization (12) with parameters $\kappa$, $\delta$, and $\gamma$, and the same assumptions as in Theorem 1. As $p\rightarrow\infty$, for any locally-Lipschitz function $\Psi(\cdot,\cdot)$, the following convergence holds,

where $Z$ is standard norma, $\beta\sim\Pi$, and $\bar{\alpha}$,$\bar{\sigma}$ are the unique solution to the following nonlinear system of equations,

The proof is deferred to the Appendix. Theorem 2 states that upon centering the estimate $\hat{\bm{\beta}}$, it becomes decorrelated from $\bm{\beta}^{*}$ and the distribution of the entries approach a zero-mean Gaussian distribution with variance ${\bar{\sigma}}^{2}$. Figure 1 depicts the performance of the regularized estimate for different values of $\lambda$. As observed in the figure, increasing the value of $\lambda$ reduces the correlation factor $\bar{\alpha}$ (Figure 1(a)) and the variance ${\bar{\sigma}}^{2}$ (Figure 1(b)). Figure 1(c) shows the mean-squared-error of the estimate as a function of $\lambda$ . It indicates that for different values of $\delta$ there exist an optimal value $\lambda_{\text{opt}}$ that achieves the minimum mean-squared error.

When $\lambda=0$ in (12), we obtain the optimization with no regularization, i.e., the maximum likelihood estimate. When we set $\lambda$ to zero in (16), Theorem 2 gives the same result as Sur and Candes reported in . In their analysis, they have also provided an interesting interpretation of $\bar{\gamma}$ in terms of the likelihood ratio statistics. Studying the likelihood ratio test is beyond the scope of this paper.

Sparse Logistic Regression

In Section 5.1, we explicitly describe the expectations in the nonlinear system (6) using two $q$-functionsThe $q$-function is the tail distribution of the standard normal r.v. defined as, $Q(t):=\int_{t}^{\infty}\frac{e^{-x^{2}/2}}{\sqrt{2\pi}}dx$ .. In Section 5.2, we analyze the support recovery in the resulting estimate and show that the two $q$-functions represent the probability of on and off support recovery.

For our analysis in this section, we assume each entry $\bm{\beta}^{*}_{i}$, for $i=1,\ldots,p$, is sampled i.i.d. from a distribution,

In the next section, we provide an interpretation for $t_{1}$ and $t_{2}$. In particular, we will show that $Q(\bar{t_{1}})$, and $Q(\bar{t_{2}})$ are related to the probabilities of on and off support recovery. We can rewrite the first three equations in (6) as follows,

Appending the three equations in (20) to the last three equations in (6) gives the nonlinear system for sparse LR. Upon solving these equations, we can use the result of Theorem 1 to compute various performance measure on the estimate $\hat{\bm{\beta}}$.

Figure 2 shows the performance of our estimate as a function of $\lambda$. It can be seen that the bound derived from our theoretical result matches the empirical simulations. Also, it can be inferrred from Figure 2(c) that the optimal value of $\lambda$ ($\lambda_{\text{opt}}$ that achieves the minimum mean-squared error) is a decreasing function of $\delta$.

2 Support recovery

In this section, we study the support recovery in sparse LR. As mentioned earlier, sparse LR is often used when the underlying paramter has few non-zero entries. We define the support of $\bm{\beta}^{*}$ as $\Omega:=\{j|1\leq j\leq p,\bm{\beta}^{*}_{j}\neq 0\}$. Here, we would like to compute the probability of success in recovery of the support of $\bm{\beta}^{*}$. Let $\hat{\bm{\beta}}$ denote the solution of the optimization (17). We fix the value $\epsilon>0$ as a hard-threshold based on which we decide whether an entry is on the support or not. In other words, we form the following set as our estimate of the support given $\hat{\bm{\beta}}$,

In order to evaluate the success in support recovery, we define the following two error measures,

In our estimation, $E_{1}$ represents the probability of false alarm, and $E_{2}$ is the probability of misdetection of an entry of the support. The following lemma indicates the asymptotic behavior of both errors as $\epsilon$ approcahes zero .

Let $\hat{\bm{\beta}}$ be the solution to the optimization (17), and the entries of $\bm{\beta}^{*}$ have distribution $\Pi$ defined in (18). Assume $\lambda$ is chosen such that the nonlinear system (6) has a unique solution $(\bar{\alpha},\bar{\sigma},\bar{\gamma},\bar{\theta},\bar{\tau},\bar{r})$. As $p\rightarrow\infty$ we have,

Conclusion and Future Directions

References

Appendix

Appendix A Convex Gaussian Min-max Theorem (CGMT)

Our analysis is based on the convex gaussian min-max theorem (CGMT). Here, we formally state this theorem. The CGMT associates with a Primary Optimization (PO) problem an Auxiliary Optimization (AO) problem from which we can investigate various properties of the primary optimization, such as the phase transition. In particular, the (PO) and the (AO) problems are defined respectively as follows:

In (24), let $\mathcal{S}_{\mathbf{w}}$, $\mathcal{S}_{\mathbf{u}}$, be convex and compact sets, and assume $\psi(\cdot,\cdot)$ is convex-concave on $\mathcal{S}_{\mathbf{w}}\times\mathcal{S}_{\mathbf{u}}$. Also assume that $\mathbf{G},~{}\mathbf{g},$ and $\mathbf{h}$ all have entries i.i.d. standard normal. The following statements are true,

Let $\mathcal{S}$ be an arbitrary open subset of $\mathcal{S}_{\mathbf{w}}$ and $\mathcal{S}^{c}:=\mathcal{S}_{\mathbf{w}}/\mathcal{S}$. Denote $\Phi_{\mathcal{S}^{c}}(\mathbf{G})$ and $\phi_{\mathcal{S}^{c}}(\mathbf{g},\mathbf{h})$ be the optimal costs of the optimizations in (24a), and (24a), respectively, when the minimization over $\mathbf{w}$ is now constrained over $\mathbf{w}\in\mathcal{S}^{c}$. If there exists constants $\bar{\phi}$, ${\bar{\phi}}_{\mathcal{S}^{c}}$, and $\eta>0$ such that,

${\bar{\phi}}_{\mathcal{S}^{c}}\geq\bar{\phi}+3\eta$ ,

$\phi(\mathbf{g},\mathbf{h})<\bar{\phi}+\eta$, with probability at least $1-p$ ,

$\phi_{\mathcal{S}^{c}}(\mathbf{g},\mathbf{h})>{\bar{\phi}}_{\mathcal{S}^{c}}-\eta$, with probability at least $1-p$ ,

The probabilities are taken with respect to the randomness in $\mathbf{G}$, $\mathbf{g}$, and $\mathbf{h}$.

We also use the following corollary that is true in the asymptotic regime,

using the same notations and assumptions as in Theorem 3, suppose there exists constants $\bar{\phi}<{\bar{\phi}}_{\mathcal{S}^{c}}$ such that $\phi(\mathbf{g},\mathbf{h})\overset{p}{\longrightarrow}\bar{\phi}$, and $\phi_{\mathcal{S}^{c}}(\mathbf{g},\mathbf{h})\longrightarrow{\bar{\phi}}_{\mathcal{S}^{c}}$. Then,

We refer the interested reader to for furder reading on the subject, its premises and applications.

Appendix B Useful Mathematical Tools

We gathered here some useful lemmas that are used in the proof of our main results. The first lemma provides the partial derivatives of the Moreau envelope function.

and the proximal operator is the solution to this optimization, i.e.,

The derivative of the Moreau envelope function can be computed as follows,

We refer the interested reader to for the proof as well as a detailed study of the properties of the Moreau envelope.

The next two lemmas present some properties of the proximal operator for the function $\rho(z)=\log(1+e^{z})$.

Let $\rho(z)=\log(1+e^{z})$, then the following identity holds,

Since the function $\rho(\cdot)$ is differentiable the proximal operator satisfies the following equation,

The derivative of the proximal operator of the function $\rho(\cdot)$ can be computed as follows,

Taking derivative with respect to $x$ of (31),

Appendix C Proof of Theorem 1

We present the proof of our main result that is a precise characterization on the performance of the optimization program (5) in the limit where $p,n\rightarrow\infty$ at a fixed ratio $\delta:=\frac{n}{p}$. We assume the data points are drawn independently from Gaussian distribution, $\mathbf{x}_{i}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,\frac{1}{p}\mathbf{I}_{p})$. We first rewrite (5) as follows,

Note that the matrix $\mathbf{H}$ is defined in such a way that its entries have i.i.d. standard normal distribution. We use the CGMT framework for our analysis. The proof strategy consists of three main steps:

Finding the auxiliary optimization: In order to apply the result of Theorem 3, we need to rewrite the optimization as a bilinear form and find its associated auxiliary optimization.

Analyzing the auxiliary optimization: The goal of this step is to simplify the auxiliary optimization in such a way that its performance can be characterized via a scalar optimization.

Finding the optimality condition on the scalar optimization: We investigate the solution to the resulting scalar optimization. Specifically, by writing the first-order optimality conditions, we will derive the nonlinear system of equations (6).

We explain each of the three steps in more details in the following subsections.

In order to apply the CGMT, we need to have a min-max optimization. Introducing a new variable $\mathbf{u}$, we have the following optimization,

Next, we use the Lagrange multiplier $\mathbf{v}$ to rewrite (39) as a min-max optimization,

Since $\mathbf{y}$ depends on $\mathbf{H}$ we can not directly apply CGMT to the bilinear form $\mathbf{v}^{T}\mathbf{H}\bm{\beta}$. To solve this issue, we first introduce, $\mathbf{P}:=\frac{1}{{||\bm{\beta}^{*}||}_{2}^{2}}\bm{\beta}^{*}{\bm{\beta}^{*}}^{T}$, and $\mathbf{P}^{\perp}:=\mathbf{I}_{p}-\mathbf{P}$, the projection matrices on the direction of $\bm{\beta}^{*}$ and its orthogonal complement, respectively. We use these projections to decompose the matrix $\mathbf{H}$ as, $\mathbf{H}=\mathbf{H}_{1}+\mathbf{H}_{2}$, with $\mathbf{H}_{1}:=\mathbf{H}\times\mathbf{P}$, and $\mathbf{H}_{2}:=\mathbf{H}\times\mathbf{P}^{\perp}$. Rewriting (40) with the decomposition of $\mathbf{H}$ would give,

It is worth noting that after performing this decomposition, the label vector ($\mathbf{y}$) would be independent of $\mathbf{H}_{2}$ since,

where we used $\mathbf{P}\bm{\beta}^{*}=\bm{\beta}^{*}$. Exploiting this fact, one can check that all the additive terms in the objective function of (41) except the last one are independent of $\mathbf{H}_{2}$. Also, the objective function is convex with respect to $\bm{\beta}$ and $\mathbf{u}$ and concave with respect to $\mathbf{v}$. In order to apply the CGMT framework, we only need an extra condition which is restricting the feasible sets of $\bm{\beta},\mathbf{u}$, and $\mathbf{v}$ to be compact and convex. We can introduce some artificial convex and bounded sets $\mathcal{S}_{\mathbf{u}}$, $S_{\mathbf{v}}$, and $\mathcal{S}_{\bm{\beta}}$, and perform the optimization over these sets. Note that these sets can be chosen large enough such that they do not affect the optimization itself. For simplicity, in our arguments here we ignore the condition on the compactness of the fesible sets and apply the CGMT whenever our feasible sets are convex.

The optimization program (41) is suitable to be analyzed via the CGMT as the conditions are all satisfied. Having identified (41) as the (PO) in our optimization, it is straightforward to write its corresponding (AO) as in (24). Therefore, the Auxiliary Optimization (AO) can be written as follows,

C.2 Analyzing the auxiliary optimization

In this section, we analyze the auxiliary optimization (C.1). Ideally, we would like to solve the optimizations with respect to the direction of the vectors, in order to finally get a scalar-valued optimization over the magnitude of the variables. Proceeding onwards, we first perform the maximization with respect to the direction of $\mathbf{v}$. We can write the following maximization with respect to $\mathbf{v}$,

In order to maximize the objective function, $\mathbf{v}$ chooses its direction to be the same as the vector it is multiplied to. Define $r:=||\mathbf{v}||/\sqrt{n}$, then maximizing over the direction of $\mathbf{v}$ would give,

where $\delta:=\frac{n}{p}$ is the oversampling ratio. Next, we use a trick adopted from where by introducing two new scalar variables, namely $\upsilon$ and $\tau$, we can change $||\cdot||$ to $||\cdot||^{2}$ which simplifies the next steps of our analysis. The new optimization would be,

Next, in order to compute the optimal $\mathbf{w}$, we use the following completion of squares,

where we also used the following equality:

Consequently, by flipping the order of $\min$ and $\max$, we first compute the minimization with respect to $\bm{\mu}$. Hence, the optimal $\bm{\mu}$ would be the solution to the following optimization:

Using the Lagrange multiplier $\theta$ we can rewrite this optimization as,

Applying the completion of squares we have,

where we omit the term $\frac{1}{p}{\mathbf{g}}^{T}\bm{\beta}^{*}=\mathcal{O}(\frac{1}{\sqrt{p}})$ as its negligible compare to the other terms (which are of constant orders). We are able to represent the solution of (53) in terms of the Moreau envelope of the function $f(\cdot)$ as follows,

Substituting (56) in (51), we have the following optimization:

We now focus on the optimization with respect to $\mathbf{u}$. Recall that \mathbf{y}=Ber\big{(}\rho^{\prime}(\frac{1}{\sqrt{p}}\mathbf{H}\bm{\beta}^{*})\big{)}=Ber\big{(}\rho^{\prime}(\kappa\mathbf{q})\big{)}. We are interested in solving the following optimization:

Similar to the previous steps, we first do a completion of squares as follows,

Next, we use the distribution of $\mathbf{y}$ to simplify the expressions in the right-hand side of (59). We can write,

where $Z\sim\mathcal{N}(0,1)$. Also note that we can ignore the term $\frac{\sigma}{n}\mathbf{y}^{T}\mathbf{h}$ since it is of order $\frac{1}{\sqrt{n}}$. Hence, we are able to rewrite the optimization (58) with respect to $\mathbf{u}$ in the following form:

We can rewrite the equation (62) in terms of the Moreau envelope, $M_{\rho(\cdot)}$, as follows,

As the last step, we want to analyze the convergence properties of (AO). Recall that $f(\cdot)$ is a separable function. Therefore, using the result of Lemma 6, we have:

Using the strong law of large numbers, we have,

where $Z$ is a standard normal random variable and $\beta\sim\Pi$ is independent of $Z$. Similarly, we can write,

We appeal to Lemma 9 in Appendix A of to analyze the convergence properties of (AO). Due to the convergence we are getting from the LLN, applying this lemma enables us to replace the Moreau envelopes with their expected value. Hence, We need to analyze the following optimization,

C.3 Finding the optimality condition of the scalar optimization

In this section, we conclude the proof of the main result of the paper. For this, we need to show that the optimizer of the optimization (C.2) can be found by solving the nonlinear system of equations (6). Let $C(\alpha,\sigma.r,\tau,\upsilon,\theta)$ denote the objective function in (C.2). We want to find the optimer of $C(\cdot)$, i.e., the point $(\alpha^{\star},\sigma^{\star},r^{\star},\tau^{\star},\upsilon^{\star},\theta^{\star})$. Since the objective function is smooth, when the optimal values are all non-zero, they should satisfy the first order optimality condition, i.e.,

We will show that the (68) would simplify to our system of nonlinear equations. We start by putting the derivative w.r.t. $\theta$ equal to zero. We have the following,

where we used Lemma 2 for taking the derivative of the Moreau envelope, $M_{\lambda f(\cdot)}$. We can simplify (69) and reqrite it as follows,

Next, we take derivative of the objective function $C(\cdot)$ w.r.t. $r$ and $\upsilon$ and put that equal to zero. We state the following lemma which will be exploited in taking the derivatives.

, then the derivative of $F(\cdot)$ would be as follows:

In order to compute the last derivative we exploit Lemma 2. We have,

Next, we use the result of Lemma 7 to find the optimality conditions with respect to $r$ and $\upsilon$. We have,

In order to simplify the equations, we define a new variable $\gamma:=\frac{1}{r\upsilon}$. We can rewrite the equations (76) as follows,

So, far we have shown that three of the optimality conditions are the same as the nonlinear equations $1$,$2$, and $5$ in (6). Next, we take the derivative w.r.t. $\tau$. We have,

The derivative of the expected Moreau envelope can be computed as follows,

which is the third equation in the nonlinear system (6). Next, putting the derivative w.r.t. $\sigma$ equal zero gives the following,

We can compute the partial derivative of the expected Moreau envelopes as follows,

To derive the last equality, we used Lemma 4 and Lemma 5. Replacing (82), and (83) in (81) gives,

As the last step, we take the derivative with respect to $\alpha$ in order to derive the fourth equation in the nonlinear system (6). We have,

Using Stein’s lemma, we can rewrite the RHS as,

where we exploit (84) to derive the last equation. Substituting in (87) would give,

Therefore, we have shown that the nonlinear system (6) is equivalent to the optimality condition in (C.2).

Recall in the process of simplifying (AO) in Section C.2, we introduced the Moreau envelope of $f(\cdot)$ in (56). The optimizer of that Moreau envelope gives the solution of the Auxiliary optimization. Let $(\bar{\alpha},\bar{\sigma},\bar{\gamma},\bar{\theta},\bar{\tau},\bar{r})$ be the unique solution of the nonlinear system. Hence, we can present the solution of the (AO) in terms of the proximal operator as follows,

As the last step we want to show the convergence of the locally-Lipschitz function $\Psi(\cdot,\cdot)$. Recall in Section C.1, in order to apply the CGMT, we have introduced some artificial bounded sets on the optimization variables and state that we can perform the optimization over these sets. Considering the variables belong to those bounded sets, we can state the function $\Psi(\cdot,\cdot)$ is Lipschitz, as constraining a locally-Lipschitz function to a bounded set gives a Lipschitz function. Next, using the strong law of large numbers along with the fact that the entries of $\bm{\beta}^{*}$ are i.i.d. and drawn from distribution $\Pi$, we have,

where $Z$ is a standard normal random variable and $\beta\sim\Pi$ is independent of $Z$.

Exploiting the assymptotic convergence of CGMT (Corollary 26), we can introduce the set $\mathcal{S}$ as follows,

The convergence in (91) would establish that as $p\rightarrow\infty$, ${\hat{\bm{\beta}}}^{AO}\in\mathcal{S}$ with probability approaching $1$. Therefore, as the result of Corollary 26, $\hat{\bm{\beta}}={\hat{\bm{\beta}}}^{PO}\in\mathcal{S}$ with probability approaching $1$. This concludes the proof.

Appendix D Proof of Theorem 2

This result can be derived using the result of Theorem 1. We just need to show that the parameters $\theta$, $r$, and $\tau$ can be explicitely computed from the first three equations in the nonlinear system (6). Recall that we characterize the performance of the RLR in terms of the solution of the following nonlinear equation,

Replacing in the first equation of (93) gives,

and finally from the thrid equation in (93) we can compute,

We can rewrite the equations (95), (96), and (97) as follows,

Replacing the derived expressions in (98) for $\theta,r$, and $\tau$ in the last three equations of (93) would gives the following system of three equations with three unknowns,