Statistical mechanical analysis of sparse linear regression as a variable selection problem

Direct computation of the entropy is not easy, and a more systematic method of evaluation is available using a Legendre transform of the entropy, which we call free entropy throughout this study [40].

We introduce a partition function $\newcommand{\V}[1]{\boldsymbol{#1}} G(\mu|\rho, \V{y}, A)$ as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\MSEy}{\epsilon_{y}} \newcommand{\rb}{\right)} \newcommand{\lb}{\left(} \newcommand{\V}[1]{\boldsymbol{#1}} \displaystyle G(\mu|\rho,\V{y},A) \equiv \lb \prod_{i=1}^{N}\sum_{c_i=0,1} \rb \delta\lb \sum_{i}c_i-N\rho \rb {\rm e}^{-M\mu \MSEy(\V{c}|\V{y},A)}, \label{G} \nonumber \end{align} \tag{ 7 }$$

where $\delta(\cdot)$ is the delta function. The free entropy is represented by the exponential rate of G. Considering the definition of the entropy and assuming that the saddle-point method is applicable in the large N limit, we can easily see that the following relation holds

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\MSEp}{\epsilon_{y}} \newcommand{\supp}[1]{{\rm supp}({#1})} \newcommand{\rbb}{\right\}} \newcommand{\lbb}{\left\{} \newcommand{\rb}{\right)} \newcommand{\lb}{\left(} \newcommand{\V}[1]{\boldsymbol{#1}} \displaystyle g(\mu|\rho,\V{y},A)\equiv \frac{1}{N}\log G(\mu|\rho,\V{y},A)=\max_{\MSEp \in \supp{s(\MSEp|\rho,\V{y},A)}} \lbb -\alpha \mu \MSEp + s(\MSEp|\rho,\V{y},A) \rbb, \label{g} \nonumber \end{align} \tag{ 8 }$$

where $\newcommand{\supp}[1]{{\rm supp}({#1})} \supp{f(x)}$ denotes the support of the function $f(x)$. Hence g and s are the Legendre transforms of each other. Assuming $\newcommand{\e}{{\rm e}} \newcommand{\MSEp}{\epsilon_{y}} s(\MSEp)$ is concave with respect to $\newcommand{\e}{{\rm e}} \newcommand{\MSEp}{\epsilon_{y}} \MSEp$, then g and s have a one-to-one correspondence and are connected by the control parameter $\mu(\geqslant 0)$. This control parameter plays the role of 'inverse temperature' in physics. The maximiser in equation (8) and the corresponding entropy, $\newcommand{\V}[1]{\boldsymbol{#1}} \newcommand{\e}{{\rm e}} \newcommand{\MSEp}{\epsilon_{y}} \MSEp(\mu|\rho, \V{y}, A)$ and $\newcommand{\V}[1]{\boldsymbol{#1}} \newcommand{\e}{{\rm e}} \newcommand{\MSEp}{\epsilon_{y}} s(\mu|\rho, \V{y}, A)=s(\MSEp(\mu|\rho, \V{y}, A)|\rho, \V{y}, A)$, are thus parameterised as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\MSEp}{\epsilon_{y}} \newcommand{\V}[1]{\boldsymbol{#1}} \newcommand{\Part}[3]{\frac{\partial^{#3} #1}{\partial #2^{#3}}} \displaystyle \MSEp(\mu|\rho,\V{y},A)=-\frac{1}{\alpha}\Part{}{\mu}{}g(\mu|\rho,\V{y},A), \nonumber \end{align} \tag{ 9a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\V}[1]{\boldsymbol{#1}} \newcommand{\Part}[3]{\frac{\partial^{#3} #1}{\partial #2^{#3}}} \displaystyle s(\mu|\rho,\V{y},A)=g(\mu|\rho,\V{y},A)-\mu\Part{}{\mu}{}g(\mu|\rho,\V{y},A). \nonumber \end{align} \tag{ 9b }$$

Employing these relations, we can easily handle the dominant output MSE and determine the corresponding entropy from g.

However, there are two difficulties in the evaluation of g. One is the dependence on $\newcommand{\V}[1]{\boldsymbol{#1}} \V{y}$ and A, and the other is the presence of the least squares problem in the definition of $\newcommand{\V}[1]{\boldsymbol{#1}} \newcommand{\e}{{\rm e}} \newcommand{\MSEy}{\epsilon_{y}} \MSEy(\V{c})$.

The first problem is overcome as follows. The free entropy g is a self-averaging quantity, as is the entropy. Typical values of self-averaging quantities are in accordance with their averaged values in the large N limit. This means that we can calculate the averaged value of g over A and $\newcommand{\V}[1]{\boldsymbol{#1}} \V{y}=A\V{x}_0+\V{\xi}$ instead of directly considering the entropy's dependence on those quantities. We denote the average over $\newcommand{\V}[1]{\boldsymbol{#1}} \V{x}_0, \V{\xi}$ and A by square brackets with appropriate subscripts as $\newcommand{\V}[1]{\boldsymbol{#1}} [\cdots ]_{\V{x}_0, \V{\xi}, A}$. Unfortunately, this average is not easy to calculate. The so-called replica method is a great aid in such a situation, and is symbolised by the following identity

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rsb}{\right]} \newcommand{\lsb}{\left[} \newcommand{\V}[1]{\boldsymbol{#1}} \displaystyle g(\mu|\rho)=\lsb g(\mu|\rho,\V{y},A) \rsb_{\V{x}_0,\V{\xi},A} = \lim_{n\to 0}\frac{1}{Nn}\log \lsb G^n(\mu|\rho,\V{y},A) \rsb_{\V{x}_0,\V{\xi},A}. \label{replica trick} \nonumber \end{align} \tag{ 10 }$$

In addition to this identity, we assume that n is a positive integer. This assumption enables us to compute the average $\newcommand{\V}[1]{\boldsymbol{#1}} [\cdots ]_{\V{x}_0, \V{\xi}, A}$. After computing this average, we take the limit $n\to 0$ by employing the analytical continuation from $n \in \mathbb{N}$ to $n \in \mathbb{R}$ under the so-called replica symmetric (RS) or the replica symmetry breaking (RSB) ansatz, which will be explained later.

The second problem is solved by introducing a variable β and taking a limit as follows

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\MSEy}{\epsilon_{y}} \newcommand{\mc}[1]{\mathcal{#1}} \newcommand{\V}[1]{\boldsymbol{#1}} \displaystyle {\rm e}^{-M\mu \MSEy(\V{c}|\V{y},A)} =\lim_{\beta \to \infty} {\rm e}^{-\mu \mc{H}(\V{c}|\beta,\V{y},A) }, \nonumber \end{align} \tag{ 11 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\mc}[1]{\mathcal{#1}} \newcommand{\rbb}{\right\}} \newcommand{\lbb}{\left\{} \newcommand{\rb}{\right)} \newcommand{\lb}{\left(} \newcommand{\V}[1]{\boldsymbol{#1}} \displaystyle \mc{H}(\V{c}|\beta,\V{y},A) =-\frac{1}{\beta} \log \int{\rm d}\boldsymbol{x} \prod_{i}\lbb (1-c_i)\delta(x_i) + c_i \rbb {\rm e}^{-\frac{\beta}{2}||\boldsymbol{y}-\boldsymbol{A}(\boldsymbol{c}\circ\boldsymbol{x})||_2^2}, \label{Hamiltonian} \nonumber \end{align} \tag{ 12 }$$

where the factor $\newcommand{\lb}{\left(} \newcommand{\rb}{\right)} \newcommand{\lbb}{\left\{} \newcommand{\rbb}{\right\}} \prod_{i}\lbb (1-c_i)\delta(x_i) + c_i \rbb$ is introduced to make the integral well-defined and can be regarded as a prior for $\newcommand{\V}[1]{\boldsymbol{#1}} \V{x}$. In the limit $\beta \to \infty$, only the contribution corresponding to the solution of the least squares problem in equation (3) survives the integration, and $\newcommand{\V}[1]{\boldsymbol{#1}} \newcommand{\mc}[1]{\mathcal{#1}} \newcommand{\e}{{\rm e}} \newcommand{\MSEy}{\epsilon_{y}} \mc{H}(\V{c}|\beta, \V{y}, \V{A}) \to M \MSEy(\V{c}|\V{y}, \V{A})$. We further assume that $\nu=\mu/\beta$ is a positive integer as well as n in equation (10). This enables us to treat in parallel the summation over $\newcommand{\V}[1]{\boldsymbol{#1}} \V{c}$ and the integration over $\newcommand{\V}[1]{\boldsymbol{#1}} \V{x}$, as well as the average $\newcommand{\V}[1]{\boldsymbol{#1}} [\cdots ]_{\V{x}_0, \V{\xi}, A}$. The limit $\beta\to \infty \Leftrightarrow \nu \to 0$ is taken after those operations through the analytic continuation.

These operations can be summarised in a line

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Tr}[1]{\mathop{\rm Tr}_{#1}} \newcommand{\re}{{\rm Re}} \newcommand{\T}[1]{\tilde{#1}} \newcommand{\rsb}{\right]} \newcommand{\lsb}{\left[} \newcommand{\rbb}{\right\}} \newcommand{\lbb}{\left\{} \newcommand{\rb}{\right)} \newcommand{\lb}{\left(} \newcommand{\V}[1]{\boldsymbol{#1}} \displaystyle g(\mu|\rho)= \lim_{n \to 0} \lim_{\nu \to 0} \frac{1}{Nn} \log \lsb \lbb \Tr{\V{c}} \lb \Tr{\V{x}|\V{c}} {\rm e}^{-\frac{1}{2}\frac{\mu}{\nu}||\V{y}-\V{A}(\V{c}\circ \V{x}) ||_2^2 } \rb^{\nu} \rbb^n \rsb_{\V{x}_0,\V{\xi},\V{A}}, \label{phi-replica} \nonumber \end{align} \tag{ 13 }$$

with abbreviations $\newcommand{\V}[1]{\boldsymbol{#1}} \newcommand{\lb}{\left(} \newcommand{\rb}{\right)} \newcommand{\T}[1]{\tilde{#1}} \newcommand{\re}{{\rm Re}} \newcommand{\Tr}[1]{\mathop{\rm Tr}_{#1}} \Tr{\V{c}}=\lb \prod_{i=1}^{N}\sum_{c_i=0, 1} \rb \delta\lb \sum_{i}c_i-N\rho \rb$ and $\newcommand{\V}[1]{\boldsymbol{#1}} \newcommand{\lb}{\left(} \newcommand{\rb}{\right)} \newcommand{\lbb}{\left\{} \newcommand{\rbb}{\right\}} \newcommand{\T}[1]{\tilde{#1}} \newcommand{\re}{{\rm Re}} \newcommand{\Tr}[1]{\mathop{\rm Tr}_{#1}} \Tr{\V{x}|\V{c}}=\int{\rm d}\V{x} \prod_{i}$$\newcommand{\V}[1]{\boldsymbol{#1}} \newcommand{\lb}{\left(} \newcommand{\rb}{\right)} \newcommand{\lbb}{\left\{} \newcommand{\rbb}{\right\}} \newcommand{\T}[1]{\tilde{#1}} \newcommand{\re}{{\rm Re}} \newcommand{\Tr}[1]\left \lbb (1-c_i)\delta(x_i) + c_i \rbb$. Overall, to calculate the entropy, we assess the free entropy g. The averages over $\newcommand{\V}[1]{\boldsymbol{#1}} \V{x}_0, \V{\xi}, A$ are taken through the replica method with a replica number n, and the internal variables $\newcommand{\V}[1]{\boldsymbol{#1}} \V{x}$ are integrated with an additional replica number ν.

The description above is generic, but for technical reasons we need additional assumptions to complete the computation. The most crucial assumption is applied to the distribution of A: each component of A is assumed to be i.i.d. from $\newcommand{\mc}[1]{\mathcal{#1}} \mc{N}(0, 1/N)$. Relaxing this assumption, i.e. introducing correlations between components, makes the analysis much more complicated. Admittedly, this assumption is not necessarily realistic. However, the purpose of this paper is to provide an analytical basis to understand the variable-selection performance, and we consider that the random-matrix assumption can provide sufficiently nontrivial implications for this purpose.

Thus assuming the absence of correlations among components of A, we note that only the average behaviour of the signal components is relevant. According to this observation, without loss of generality, we may assume a factorised prior of the signal vector $\newcommand{\V}[1]{\boldsymbol{#1}} \V{x}_0$ as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rbb}{\right\}} \newcommand{\lbb}{\left\{} \newcommand{\rb}{\right)} \newcommand{\lb}{\left(} \newcommand{\V}[1]{\boldsymbol{#1}} \displaystyle P(\V{x}_0)=\prod_{i=1}^{N}\lbb (1-\rho_0)\delta(x_{0i}) + \rho_0P_{0}(x_{0i})\rbb, \label{factorised} \nonumber \end{align} \tag{ 14 }$$

where $\newcommand{\nnzero}{K} \rho_0=\nnzero_0/N$ is the density of nonzero components and P₀(x) is a prior distribution for the nonzero component. Our theoretical computation can be performed for any prior P₀(x) having no probability mass at x  =  0, and we keep it unspecified for a while.

Usually, the entropy function $\newcommand{\e}{{\rm e}} \newcommand{\MSEp}{\epsilon_{y}} s(\MSEp)$ enjoys some useful properties such as non-negativity, boundedness, bounded support, concavity, and analyticity. We assume these properties, but as shown below, the concavity and analyticity are partially broken in the present problem. This causes problems in evaluating the parametric form (9), but they can be bypassed by some additional considerations when conducting the saddle-point method. The analyticity breaking of the entropy is actually related to algorithmic performances and is one of the central issues discussed in this paper.

Our formulation has a probabilistic meaning which involves two different intrinsic hierarchies in the present problem. We can define the distribution for $\newcommand{\V}[1]{\boldsymbol{#1}} \V{c}$ as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\MSEy}{\epsilon_{y}} \newcommand{\rb}{\right)} \newcommand{\lb}{\left(} \newcommand{\V}[1]{\boldsymbol{#1}} \displaystyle P^{(1)}(\V{c}|\mu,\rho,\V{y},A)=\frac{1}{G(\mu|\rho,\V{y},A)}\delta\lb \sum_{i}c_i-N\rho \rb {\rm e}^{-M\mu \MSEy(\V{c}|\V{y},A)}. \label{P1} \nonumber \end{align} \tag{ 15 }$$

Let us denote the average over P⁽¹⁾ by angular brackets as $\newcommand{\V}[1]{\boldsymbol{#1}} \newcommand{\Ave}[1]{\left\langle {#1} \right\rangle} \Ave{\cdots}_{\V{c}}$. This distribution is clearly conditioned by $\newcommand{\V}[1]{\boldsymbol{#1}} \V{y}$ and A. We also define another distribution for $\newcommand{\V}[1]{\boldsymbol{#1}} \V{x}$ given $\newcommand{\V}[1]{\boldsymbol{#1}} \V{c}$ as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rbb}{\right\}} \newcommand{\lbb}{\left\{} \newcommand{\rb}{\right)} \newcommand{\lb}{\left(} \newcommand{\V}[1]{\boldsymbol{#1}} \displaystyle P^{(2)}(\V{x}|\beta,\V{c},\V{y},A)= \frac{1}{Z} \prod_{i}\lbb (1-c_i)\delta(x_i)+c_i \rbb {\rm e}^{-\frac{\beta}{2} ||\V{y}-A\lb \V{c}\circ\V{x} \rb||_2^2}, \nonumber \end{align} \tag{ 16 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rbb}{\right\}} \newcommand{\lbb}{\left\{} \newcommand{\rb}{\right)} \newcommand{\lb}{\left(} \newcommand{\V}[1]{\boldsymbol{#1}} \displaystyle Z=\int {\rm d} \V{x} \prod_{i}\lbb (1-c_i)\delta(x_i)+c_i \rbb {\rm e}^{-\frac{\beta}{2} ||\V{y}-A\lb \V{c}\circ\V{x} \rb||_2^2} ={\rm e}^{-\beta \mathcal{H}(\V{c}|\beta,\V{y},A)} . \nonumber \end{align} \tag{ 17 }$$

This distribution is conditioned by $\newcommand{\V}[1]{\boldsymbol{#1}} \V{c}$ in addition to $\newcommand{\V}[1]{\boldsymbol{#1}} \V{y}$ and A. Hence, we denote the average over P⁽²⁾ by $\newcommand{\V}[1]{\boldsymbol{#1}} \newcommand{\Ave}[1]{\left\langle {#1} \right\rangle} \Ave{\cdots}_{\V{x}|\V{c}}$. The simultaneous average over both P⁽¹⁾ and P⁽²⁾ is denoted by double angular brackets $\newcommand{\dAve}[1]{\left\langle \left\langle {#1} \right\rangle \right\rangle} \dAve{\cdots}$.

Recalling that $\newcommand{\V}[1]{\boldsymbol{#1}} \V{y}$ and A are also random variables, we note that there are three different hierarchies of random variables: $\newcommand{\V}[1]{\boldsymbol{#1}} \V{x}$ is conditioned by $\newcommand{\V}[1]{\boldsymbol{#1}} \V{c}$ which is conditioned by $\newcommand{\V}[1]{\boldsymbol{#1}} \V{y}$ and A. This discrimination is a natural consequence of the structure of the present problem⁵.

Postponing the details of analysis to the appendix, we present the resultant formulas of the free entropy and related quantities. The expression for g in the RS level is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rbb}{\right\}} \newcommand{\lbb}{\left\{} \newcommand{\rb}{\right)} \newcommand{\lb}{\left(} \newcommand{\Extr}[1]{\mathop{\rm Extr}_{#1}} \displaystyle g_{\rm RS}(\mu|\rho)= &\,\Extr{\Omega_{\rm RS}} \Biggl\{\tilde{\rho} \rho+\frac{1}{2}\tilde{Q}Q-\frac{\tilde{\chi}\chi}{2\mu}+\frac{1}{2}\tilde{q}q-\tilde{m}m \nonumber \\ & +\rho_0\int {\rm d}x_0P_0(x_0) \int Dz \log \lb 1+Y^{{\rm RS}}_{\tilde{m}}\rb +(1-\rho_0)\int Dz \log \lb 1+Y^{{\rm RS}}_{0}\rb \nonumber \\ & +\frac{\alpha}{2}\lbb \log\frac{1+\chi}{D_{{\rm RS}}}-\frac{\mu (V+q)}{D_{{\rm RS}}}\rbb \Biggr\}, \label{g_RS} \nonumber \end{align} \tag{ 27 }$$

where for simplicity of notation we let $\newcommand{\lb}{\left(} \newcommand{\rb}{\right)} \newcommand{\lbb}{\left\{} \newcommand{\rbb}{\right\}} \Omega_{{\rm RS}}=\lbb \chi, Q, q, m, \tilde{\rho}, \tilde{\chi}, \tilde{Q}, \tilde{q}, \tilde{m} \rbb$, $Dz={\rm d}z{\rm e}^{-\frac{1}{2}z^2}/\sqrt{2\pi}$, $\newcommand{\POWx}{\sigma_x^2} V=\rho_0 \POWx+\sigma_{\xi}^2-2m$, and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta_{{\rm RS}}=Q-q, \nonumber \end{align} \tag{ 28 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle D_{{\rm RS}}=1+\chi+\mu\Delta_{{\rm RS}} \nonumber \end{align} \tag{ 29 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle h^{{\rm RS}}_{\tilde{m}}=\tilde{m}x_0+\sqrt{\tilde{q}}z, \nonumber \end{align} \tag{ 30 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rb}{\right)} \newcommand{\lb}{\left(} \displaystyle Y^{{\rm RS}}_{\tilde{m}}=\sqrt{\frac{\tilde{\chi}+\tilde{Q}}{\tilde{Q}+\tilde{q}}} {\rm e}^{-\tilde{\rho}+\frac{1}{2}\frac{1}{\tilde{Q}+\tilde{q}}\lb h^{{\rm RS}}_{\tilde{m}} \rb^{2}}. \label{hY_RS} \nonumber \end{align} \tag{ 31 }$$

$Y^{{\rm RS}}_{0}$ is obtained by substituting $\tilde{m}=0$ into $Y^{{\rm RS}}_{\tilde{m}}$, and $h^{{\rm RS}}_0$ is defined similarly. The symbol $\newcommand{\Extr}[1]{\mathop{\rm Extr}_{#1}} \Extr{\Omega}$ denotes the extremisation condition with respect to Ω coming from the saddle-point method. This extremisation condition yields the following equations of state (EOS):

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rbb}{\right\}} \newcommand{\lbb}{\left\{} \newcommand{\rb}{\right)} \newcommand{\lb}{\left(} \displaystyle \tilde{\chi}=\alpha \lbb \frac{\mu^2\Delta_{{\rm RS}}}{(1+\chi)D_{{\rm RS}}} +\frac{\mu^2(V+q) }{D_{{\rm RS}}^2}\rbb, \label{chitilde_RS} \nonumber \end{align} \tag{ 32a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rbb}{\right\}} \newcommand{\lbb}{\left\{} \newcommand{\rb}{\right)} \newcommand{\lb}{\left(} \displaystyle \tilde{Q}=\alpha \lbb \frac{\mu }{D_{{\rm RS}}} -\frac{\mu^2 (V+q) }{D_{{\rm RS}}^2}\rbb, \label{Qtilde_RS} \nonumber \end{align} \tag{ 32b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{q}=\alpha \frac{\mu^2 (V+q)}{D_{{\rm RS}}^2}, \label{qtilde_RS} \nonumber \end{align} \tag{ 32c }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{m}= \frac{\alpha \mu}{D_{{\rm RS}}}, \label{mtilde_RS} \nonumber \end{align} \tag{ 32d }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho=\rho_0\int {\rm d}x_0 P_0(x_0) \int Dz \frac{Y^{{\rm RS}}_{\tilde{m}}}{1+Y^{{\rm RS}}_{\tilde{m}}}+(1-\rho_0)\int Dz \frac{Y^{{\rm RS}}_{0}}{1+Y^{{\rm RS}}_{0}}, \label{rho_RS} \nonumber \end{align} \tag{ 32e }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rbb}{\right\}} \newcommand{\lbb}{\left\{} \newcommand{\rb}{\right)} \newcommand{\lb}{\left(} \displaystyle \chi=\frac{\mu}{\tilde{\chi}+\tilde{Q}} \lbb \rho_0\int {\rm d}x_0 P_0(x_0)\int Dz \frac{Y^{{\rm RS}}_{\tilde{m}}}{1+Y^{{\rm RS}}_{\tilde{m}}}+(1-\rho_0)\int Dz \frac{Y^{{\rm RS}}_{0}}{1+Y^{{\rm RS}}_{0}} \rbb, \label{chi_RS} \nonumber \end{align} \tag{ 32f }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Brbb}{\Biggr\}} \newcommand{\Blbb}{\Biggl\{} \newcommand{\Brb}{\Biggr)} \newcommand{\Blb}{\Biggl(} \newcommand{\rbb}{\right\}} \newcommand{\lbb}{\left\{} \newcommand{\rb}{\right)} \newcommand{\lb}{\left(} \displaystyle &Q=\frac{\tilde{\chi}-\tilde{q}}{(\tilde{\chi}+\tilde{Q})(\tilde{Q}+\tilde{q})} \lbb \rho_0\int {\rm d}x_0 P_0(x_0)\int Dz \frac{Y^{{\rm RS}}_{\tilde{m}}}{1+Y^{{\rm RS}}_{\tilde{m}}}+(1-\rho_0)\int Dz \frac{Y^{{\rm RS}}_{0}}{1+Y^{{\rm RS}}_{0}} \rbb \nonumber \\ & +\frac{1}{(\tilde{Q}+\tilde{q})^2}\Blbb \rho_0 \int {\rm d}x_0 P_0(x_0)\int Dz~ \frac{\lb h^{{\rm RS} }_{\tilde{m}} \rb^2 Y^{{\rm RS} }_{\tilde{m}} }{1+Y^{{\rm RS}}_{\tilde{m}}} + (1-\rho_0)\int Dz~ \frac{\lb h^{{\rm RS} }_{0} \rb^2 Y^{{\rm RS}}_{0} }{1+Y^{{\rm RS}}_{0}} \Brbb, \label{Q_RS} \nonumber \end{align} \tag{ 32g }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Brbb}{\Biggr\}} \newcommand{\Blbb}{\Biggl\{} \newcommand{\Brb}{\Biggr)} \newcommand{\Blb}{\Biggl(} \newcommand{\rb}{\right)} \newcommand{\lb}{\left(} \displaystyle q= \frac{1}{(\tilde{Q}+\tilde{q})^2} \Blbb \rho_0\int {\rm d}x_0 P_0(x_0) \int Dz~ \lb \frac{h^{{\rm RS} }_{\tilde{m}} Y^{{\rm RS}}_{\tilde{m}}}{1+Y^{{\rm RS}}_{\tilde{m}}} \rb^2 + (1-\rho_0)\int Dz~ \lb \frac{h^{{\rm RS} }_{0} Y^{{\rm RS}}_{0}}{1+Y^{{\rm RS}}_{0}} \rb^2 \Brbb , \label{q_RS} \nonumber \end{align} \tag{ 32h }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rbb}{\right\}} \newcommand{\lbb}{\left\{} \newcommand{\rb}{\right)} \newcommand{\lb}{\left(} \displaystyle m= \frac{1}{\tilde{Q}+\tilde{q}}\lbb \rho_0\int {\rm d}x_0 P_0(x_0) \int Dz \frac{x_0 h^{{\rm RS}}_{\tilde{m}} Y^{{\rm RS}}_{\tilde{m}}}{1+Y^{{\rm RS}}_{\tilde{m}}} \rbb . \label{m_RS} \nonumber \end{align} \tag{ 32i }$$

Note that all tilde variables are conjugates of their respective variables and are introduced to expand delta functions with the Fourier transform as shown in the appendix. From equation (32), we obtain some simple relations

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \chi=\frac{\rho}{\alpha-\rho}, \nonumber \end{align} \tag{ 33a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{\chi}+\tilde{Q}=\mu (\alpha-\rho), \nonumber \end{align} \tag{ 33b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{Q}+\tilde{q}=\frac{\alpha \mu}{D_{{\rm RS}}}=\tilde{m}, \nonumber \end{align} \tag{ 33c }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{\chi}-\tilde{q}=\frac{\alpha \mu^2 \Delta_{{\rm RS}}}{(1+\chi)D_{{\rm RS}}}, \nonumber \end{align} \tag{ 33d }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\MSEp}{\epsilon_{y}} \newcommand{\Part}[3]{\frac{\partial^{#3} #1}{\partial #2^{#3}}} \displaystyle \MSEp(\mu|\rho)=-\frac{1}{\alpha}\Part{g(\mu|\rho)}{\mu}{}=\frac{\tilde{\chi}}{2\alpha \mu^2}. \nonumber \end{align} \tag{ 33e }$$

The RS solution can be inaccurate when the configuration space of $\newcommand{\V}[1]{\boldsymbol{#1}} \V{c}$ exhibits spontaneous breaking into many locally separated components. In such a situation, the RSB ansatz should be adopted. The RSB solution is categorised by the level of the emerging hierarchical structure of the separation. In the simplest case called the 1st step RSB (1RSB) solution, only one level of hierarchy is taken and we examine this in the present paper. The 1RSB solution is actually sufficient to expose the instability of the RS solution and hence is sufficient to achieve the present purpose of obtaining implications to local search algorithms.

We postpone the detailed derivation of the 1RSB solution to the appendix. The explicit 1RSB formula involving g is given as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rbb}{\right\}} \newcommand{\lbb}{\left\{} \newcommand{\rb}{\right)} \newcommand{\lb}{\left(} \newcommand{\Extr}[1]{\mathop{\rm Extr}_{#1}} \displaystyle & g_{{\rm 1RSB}}(\mu,\tau;\rho)= \Extr{\Omega_{{\rm 1RSB}}} \Biggl\{\tilde{\rho} \rho+\frac{1}{2}\tilde{Q}Q-\frac{\tilde{\chi}\chi}{2\mu}-\frac{1}{2}(\tau-1)\tilde{q}_1q_1+\frac{1}{2}\tau\tilde{q}_0q_0-\tilde{m}m \nonumber \\ & +\frac{\rho_0}{\tau}\int {\rm d}x_0P_0(x_0) \int Dz_0 \log \int Dz_1 \lb 1+Y^{{\rm 1RSB}}_{\tilde{m}}\rb^{\tau} +\frac{1-\rho_0}{\tau}\int Dz_0 \log \int Dz_1 \lb 1+Y^{{\rm 1RSB}}_{0}\rb^{\tau} \nonumber \\ & +\frac{\alpha}{2} \lbb \log\frac{1+\chi}{D_1} +\frac{1}{\tau}\log\frac{D_1}{D_0} -\frac{\mu (V+q_0)}{D_0} \rbb \Biggr\}, \label{g_1RSB} \nonumber \end{align} \tag{ 34 }$$

where τ is Parisi's breaking parameter, $\newcommand{\lb}{\left(} \newcommand{\rb}{\right)} \newcommand{\lbb}{\left\{} \newcommand{\rbb}{\right\}} \Omega_{{\rm 1RSB}}=\lbb \chi, Q, q_1, q_0, m, \tilde{\rho}, \tilde{\chi}, \tilde{Q}, \tilde{q}_1, \tilde{q}_0, \tilde{m}, \tau \rbb$, and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta_1=Q-q_1,~\Delta_0=q_1-q_0, \nonumber \end{align} \tag{ 35 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle D_1=1+\chi+\mu\Delta_{1}, \nonumber \end{align} \tag{ 36 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle D_0=1+\chi+\mu\Delta_{1}+\tau\mu\Delta_0, \nonumber \end{align} \tag{ 37 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle h^{{\rm 1RSB}}_{\tilde{m}}=\tilde{m}x_0+\sqrt{\tilde{q}_1-\tilde{q}_0}z_1+\sqrt{\tilde{q}_0}z_0. \nonumber \end{align} \tag{ 38 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rb}{\right)} \newcommand{\lb}{\left(} \displaystyle Y^{{\rm 1RSB}}_{\tilde{m}}=\sqrt{\frac{\tilde{\chi}+\tilde{Q}}{\tilde{Q}+\tilde{q}_1}} {\rm e}^{-\tilde{\rho}+\frac{1}{2}\frac{1}{\tilde{Q}+\tilde{q}_1}\lb h^{{\rm 1RSB}}_{\tilde{m}} \rb^{2}}. \nonumber \end{align} \tag{ 39 }$$

As in the RS case, $h^{{\rm 1RSB}}_0$ and $Y^{{\rm 1RSB}}_{0}$ are given by inserting $\tilde{m}=0$ into $h^{{\rm 1RSB}}_{\tilde{m}}$ and $Y^{{\rm 1RSB}}_{\tilde{m}}$, respectively. The corresponding EOS are involved and are given in the appendix.

By examining the 1RSB solution, we can determine the instability points of the RS solution. Empirically, two types of instabilities are known to appear in a wide range of systems:

• Global instability  
  The RS solution is locally stable but there emerges another solution involving exponentially many metastable states, which induces the so-called random first-order transition (RFOT). We thus call the associated instability RFOT instability in this paper.
• Local instability  
  The local instability of the RS solution, which can be signaled by expanding the 1RSB solution with respect to $\Delta_0$ and observing its coefficient. This is also known as the de Almeida–Thouless (AT) instability.

According to these empirical facts, we derive two instability conditions below.

The RFOT instability is known to emerge at $\tau=1$ and can be detected in an easy manner as follows. For $\tau=1$, we can identify q₀ and $\tilde{q}_0$ with q and $\tilde{q}$ in the RS solution, respectively, because their EOS formally accord with each other. As well, the 1RSB EOS of all other order parameters except for q₁ and $\tilde{q}_1$ become identical to their corresponding RS EOS. Hence, we should compute q₁ and $\tilde{q}_1$ on top of the RS solution and examine whether the nontrivial solution $q_1\neq q_0=q$ exists or not. The equations to be solved are written in terms of $\Delta_0$ and $\newcommand{\e}{{\rm e}} \tilde{\Delta}_0\equiv \tilde{q}_1-\tilde{q}_0=\tilde{q}_1-\tilde{q}$ as follows

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{\Delta}_0=\frac{\alpha \mu^2}{1+\chi+\mu(Q-q)}\frac{\Delta_0}{1+\chi+\mu(Q-q-\Delta_0)}, \label{Delta0tilde_1RSB} \nonumber \end{align} \tag{ 40 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Brbb}{\Biggr\}} \newcommand{\Blbb}{\Biggl\{} \newcommand{\Brb}{\Biggr)} \newcommand{\Blb}{\Biggl(} \newcommand{\rb}{\right)} \newcommand{\lb}{\left(} \displaystyle &\Delta_0=\frac{1}{\lb \tilde{Q}+\tilde{q}+\tilde{\Delta}_0 \rb^2} \Blbb \rho_0\int {\rm d}x_0 P_0(x_0) \int Dz_0~ \frac{\int Dz_1 \lb h^{{\rm 1RSB} }_{\tilde{m}} Y^{{\rm 1RSB}}_{\tilde{m}} \rb^2 (1+Y^{{\rm 1RSB}}_{\tilde{m}})^{-1} }{1+Y^{{\rm RS}}_{\tilde{m}} } \nonumber \\ & + (1-\rho_0)\int Dz_0~ \frac{\int Dz_1 \lb h^{{\rm 1RSB} }_{0} Y^{{\rm 1RSB}}_{0} \rb^2 (1+Y^{{\rm 1RSB}}_{0})^{-1} }{1+Y^{{\rm RS}}_{0} } \Brbb -q. \label{Delta0_1RSB} \nonumber \end{align} \tag{ 41 }$$

In these equations, we should read $\tilde{q}_1=\tilde{q}+\tilde{\Delta}_0$ in $h^{{\rm 1RSB} }$ and $Y^{{\rm 1RSB} }$. The trivial RS solution $\Delta_0=0$ always exists and the question is whether a nontrivial solution $\Delta_0\neq 0$ exists or not. Such a nontrivial solution is absent for the low μ region but is present at sufficiently large values of μ. The lowest value of μ for which the nontrivial solution exists defines the RFOT point $\mu_{\rm RFOT}$.

The AT instability is observed by examining the presence of a nontrivial solution around $\Delta_0=0$. This can be accomplished by expanding the right hand side of equation (41) with respect to $\Delta_0$ up to the first order after inserting equation (40) into $\tilde{\Delta}_0$. If the coefficient of the first-order term is greater than unity, a nontrivial solution emerges. This condition is written using the RS solution only, because the small $\Delta_0$ limit implies that the order parameters q₁ and q₀ of the 1RSB solution can be identified as q in the RS solution, and the corresponding tilde variables can also be identified. The explicit stability condition of the RS solution against the AT instability is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rb}{\right)} \newcommand{\lb}{\left(} \displaystyle & \alpha > \rho_0\int {\rm d}x_0 P_0(x_0) \int Dz \lb \frac{Y^{{\rm RS}}_{\tilde{m}}}{1+Y^{{\rm RS}}_{\tilde{m}}} + \frac{\lb h^{{\rm RS}}_{\tilde{m}}\rb^2}{\tilde{Q}+\tilde{q}} \frac{Y^{{\rm RS}}_{\tilde{m}}}{\lb 1+Y^{{\rm RS}}_{\tilde{m}} \rb^2} \rb^2 \nonumber \\ & +(1-\rho_0)\int Dz \lb \frac{Y^{{\rm RS}}_{0}}{1+Y^{{\rm RS}}_{0}} + \frac{\lb h^{{\rm RS}}_{0} \rb^2}{\tilde{Q}+\tilde{q}} \frac{Y^{{\rm RS}}_{0} }{\lb 1+Y^{{\rm RS}}_{0} \rb^2} \rb^2. \label{AT} \nonumber \end{align} \tag{ 42 }$$

This condition always holds for sufficiently low values of μ. The lowest μ value violating equation (42) indicates the AT transition point $\mu_{\rm AT}$.

These two instabilities are known to affect the performance of local search algorithms. The origin of this degraded performance by the RFOT transition is clear—there emerge exponentially many local minima and thus the search/dynamics is easily trapped in one of those states; the typical trapping state will be the most numerous one and will be far from the true global minimum. Each trapping state in this case is separated by high energy barriers and hence escaping will take an exponentially long time [40]. Meanwhile, the influence of the AT instability is less trivial than the RFOT. According to the standard physical picture, when the AT instability occurs, the structure of the configuration space of $\newcommand{\V}[1]{\boldsymbol{#1}} \V{c}$ has many saddle-point-like structures, which leads to a complicated critical slowing down of the dynamics and thus the performance of local search algorithms will be strongly degraded. However, this degradation will be less serious compared to that caused by the RFOT transition, because in the AT case it is considered that there exist certain directions along which the system can escape from a local saddle point. This may take a long time, as the energy landscape will be very flat along the escape directions owing to the saddle-point nature, but will be shorter than the RFOT case, where an exponentially long time is required. These descriptions have actually been supported by numerical simulations of some metaheuristic algorithms in several optimisation problems [53–58]. We, however, stress that there are still many unclear points about the dynamics around the AT instability and further investigations using concrete algorithms, such as Monte-Carlo methods or message-passing algorithms, are desired.

As will be seen in the next section, we have several characteristic regions depending on the parameters ρ and $\sigma_{\xi}^2$. In some regions, the RSB transitions induced by AT and RFOT instabilities occur, which makes it difficult for the system's dynamics to converge to the equilibrium distribution. In other regions, the RS solution is always stable when μ is changed. However, the RS-stable regions are separated into two small regions, one has no phase transitions and thus the metaheuristic algorithm can work well, and the other has another 1st order transition which prevents the algorithm from approaching the global minimum. These descriptions will be actually confirmed by numerical experiments of a metaheuristic algorithm in section 4.2.

To check our replica results, we summarise some simple solutions obtained at particular limits below.

High temperature solution.

A trivial solution in the high temperature limit, $T=\mu^{-1} \to \infty$, is derived from equation (32). From simple algebra based on equations (32) and (33), we obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{Y^{{\rm RS}}_{\tilde{m}}}{1+Y^{{\rm RS}}_{\tilde{m}}}=\frac{Y^{{\rm RS}}_{0}}{1+Y^{{\rm RS}}_0}=\frac{{\rm e}^{-\tilde{\rho}}}{1+{\rm e}^{-\tilde{\rho}}}=\rho. \label{rho-mu=0} \nonumber \end{align} \tag{ 43 }$$

Accordingly,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\POWx}{\sigma_x^2} \displaystyle m= \rho \rho_0 \POWx, \nonumber \end{align} \tag{ 44a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\POWx}{\sigma_x^2} \newcommand{\rbb}{\right\}} \newcommand{\lbb}{\left\{} \newcommand{\rb}{\right)} \newcommand{\lb}{\left(} \displaystyle q=\frac{\rho^2}{\alpha-\rho^2}\lbb (1+\alpha-2\rho)\rho_0\POWx+\sigma_{\xi}^2 \rbb, \nonumber \end{align} \tag{ 44b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\POWx}{\sigma_x^2} \newcommand{\rbb}{\right\}} \newcommand{\lbb}{\left\{} \newcommand{\rb}{\right)} \newcommand{\lb}{\left(} \displaystyle Q=\frac{\rho}{\alpha-\rho}\lbb (1+\alpha-2\rho)\rho_0\POWx+\sigma_{\xi}^2 \rbb. \nonumber \end{align} \tag{ 44c }$$

Using the relation $\tilde{\rho}=-\log (\rho/(1-\rho))$, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle g(\mu=0,\rho)=-\rho \log \rho -(1-\rho)\log (1-\rho)\equiv H_2(\rho), \nonumber \end{align} \tag{ 45a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\POWx}{\sigma_x^2} \newcommand{\MSEp}{\epsilon_{y}} \newcommand{\rb}{\right)} \newcommand{\lb}{\left(} \displaystyle \MSEp(\mu=0,\rho)=\frac{1}{2}\frac{\alpha-\rho}{\alpha} \lb (1-\rho) \rho_0\POWx +\sigma_{\xi}^2 \rb, \nonumber \end{align} \tag{ 45b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle s(\mu=0,\rho)=H_2(\rho), \nonumber \end{align} \tag{ 45c }$$

where $H_2(\rho)$ is the binary entropy, giving a reasonable result.

The local stability of this solution can be checked by substituting equations (43) and (44) into equation (42), which yields

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha > \rho^2. \nonumber \end{align} \tag{ 46 }$$

This always holds in the meaningful setup of $\rho\leqslant 1$ and $\rho\leqslant \alpha$, and hence the high temperature solution is stable.

Perfect reconstruction in the noiseless limit.

Of particular interest for the noiseless limit $\sigma_{\xi}=0$ is whether we can achieve the perfect reconstruction of $\newcommand{\V}[1]{\boldsymbol{#1}} \V{x}_0$. We call the corresponding solution the perfect reconstruction (PR) solution, which is defined by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle c_i=1,~(\forall{i}~{\rm s.t.}~|x_{0i}|_0=1). \label{PRsupport} \nonumber \end{align} \tag{ 47 }$$

Note that the support components outside the true support may take unity: c_i can be 1 even if $|x_{0i}|_0=0$. This is because the coefficients $\newcommand{\V}[1]{\boldsymbol{#1}} \hat{\V{x}}(\V{c})$ outside the true support become automatically zero as a result of the optimization if all the components of the true support is covered as equation (47), leading to the vanishing MSEs. In other words, the PR solution always exists if the estimated nonzero density is greater than or equal to the true one, $\rho \geqslant \rho_0$. This seemingly indicates that it is safer to overestimate the nonzero density ρ. This is, however, not necessarily true because larger estimates of ρ tend to involve unfavourable phase transitions, as shown below.

The output MSE $\newcommand{\e}{{\rm e}} \newcommand{\MSEp}{\epsilon_{y}} \MSEp$ of the PR solution is zero and the entropy becomes

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\MSEp}{\epsilon_{y}} \displaystyle s(\MSEp=0)&=\frac{1}{N}\log \left(\begin{array}{@{}c@{}} N(1-\rho_0) \nonumber \\N(\rho-\rho_0)\end{array}\right) \nonumber \\ &= (1-\rho_0)\log(1-\rho_0)-(\rho-\rho_0)\log (\rho-\rho_0) -(1-\rho)\log(1-\rho) \equiv s_{\rm PR}. \label{s_PR} \nonumber \end{align} \tag{ 48 }$$

We can actually find this PR solution in the limit $\mu \to \infty$ of our RS formula (32). To derive this, we have to carefully treat the scaling of $V+q$ and $\Delta_{{\rm RS}}$ within that limit. We first assume the following scaling:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle ~ \mu (V+q) \to 0,~ \mu \Delta_{{\rm RS}} \to 0,~ {\rm and }~\tilde{\rho}=O(1),~(\mu \to \infty). \label{PRA} \nonumber \end{align} \tag{ 49 }$$

The consistency of this assumption is confirmed after the computation. Using this and equation (33), we can easily determine the following limits

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rb}{\right)} \newcommand{\lb}{\left(} \displaystyle \frac{\tilde{\chi}+\tilde{q} }{\tilde{Q}+\tilde{q} } \to 1, \frac{\lb h_{\tilde{m}}^{{\rm RS}} \rb^2 }{\tilde{Q}+\tilde{q} } \to (\alpha-\rho)\mu x_0^2 +2\sqrt{(\alpha-\rho)(V+q) }\mu x_0 z ,~ \frac{\lb h_{\tilde{0}}^{{\rm RS}} \rb^2 }{\tilde{Q}+\tilde{q} } \to 0. \label{factors-perfect} \nonumber \end{align} \tag{ 50 }$$

These relations mean that $Y^{{\rm RS}}_{\tilde{m}}$ diverges or vanishes depending on the values of x₀ and z while $Y^{{\rm RS}}_{0}$ converges to ${\rm e}^{-\tilde{\rho}}$. The vanishing region of $Y^{{\rm RS}}_{\tilde{m}}$ is expressed in terms of z as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle z > \frac{1}{2}\sqrt{\frac{\alpha}{(1+\chi)(V+q)}}x_0 \quad (x_0>0), \nonumber \end{align} \tag{ 51 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle z < \frac{1}{2}\sqrt{\frac{\alpha}{(1+\chi)(V+q)}}x_0 \quad (x_0<0). \nonumber \end{align} \tag{ 52 }$$

As we have assumed $(V+q)\to 0~(\mu \to \infty)$, this vanishing region rapidly shrinks and does not provide any meaningful contribution. Hence, we may treat $Y^{{\rm RS}}_{\tilde{m}}$ as a diverging factor in all contributing regions. This consideration yields, from equation (32e):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm e}^{-\tilde{\rho}} \to \frac{\rho-\rho_0}{1-\rho}. \label{rho-perfect} \nonumber \end{align} \tag{ 53 }$$

Additionally, some algebraic operations from equations (32) and (33) lead to

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Brbb}{\Biggr\}} \newcommand{\Blbb}{\Biggl\{} \newcommand{\Brb}{\Biggr)} \newcommand{\Blb}{\Biggl(} \displaystyle & (V+q)= \frac{\rho_0}{\alpha-\rho_2} \Blbb \alpha \int {\rm d}x_0 P_0(x_0)x_0^2 \int Dz (1+Y^{\rm RS}_{\tilde{m}})^{-2} \nonumber \\ & +2\sqrt{\alpha(V+q)} \int {\rm d}x_0 P_0(x_0) \int Dz~x_0 zY_{\tilde{m}}(1+Y^{\rm RS}_{\tilde{m}})^{-2} \Brbb, \nonumber \end{align} \tag{ 54 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rb}{\right)} \newcommand{\lb}{\left(} \displaystyle \rho_2\equiv \rho_0 \int {\rm d}x_0P_0(x_0) \int Dz z^2 \lb \frac{Y^{\rm RS}_{\tilde{m}}}{1+Y^{\rm RS}_{\tilde{m}}}\rb^2 + (1-\rho_0)\int Dz z^2 \lb \frac{Y^{\rm RS}_{0}}{1+Y^{\rm RS}_{0}}\rb^2, \label{rho_2} \nonumber \end{align} \tag{ 55 }$$

which yields a finite contribution in the limit $\mu \to \infty$. The divergence of $Y^{{\rm RS}}_{\tilde{m}}$ implies that

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rb}{\right)} \newcommand{\lb}{\left(} \displaystyle (V+q) \propto \int {\rm d}x_0 P_0(x_0) x_0^2 \int Dz \lb Y^{\rm RS}_{\tilde{m}} \rb^{-2}. \nonumber \end{align} \tag{ 56 }$$

The precise scaling of the right hand side strongly depends on the choice of the prior distribution $P_0(x_0)$. For examples, if it is Gaussian $\newcommand{\POWx}{\sigma_x^2} P_0(x_0)=\mathcal{N}(0, \POWx)$, then $(V+q)=O(\mu^{-\frac{3}{2}})$; if the signal strength is a constant at $P_0(x_0)=\delta(x_0-C)$ with $C\neq 0$, then $(V+q)$ exponentially decays as μ increases. Similar calculations apply to $\Delta_{{\rm RS}}$, which involves the same scaling of $\Delta_{{\rm RS}}$ as $(V+q)$. These scalings are consistent with the assumed ones (49).

Summarising these calculations in conjunction with equations (32) and (33), we determine that

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\MSEp}{\epsilon_{y}} \newcommand{\rb}{\right)} \newcommand{\lb}{\left(} \displaystyle \MSEp(\mu)=O\lb V+q \rb \to 0. \label{eps_scaling} \nonumber \end{align} \tag{ 57 }$$

This implies that the free entropy g coincides with the entropy in the limit $\mu\to \infty$. Substituting the scalings obtained so far into equation (27), we find that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \lim_{\mu\to \infty}s(\mu)=\lim_{\mu\to \infty}g(\mu)=s_{\rm PR}. \nonumber \end{align} \tag{ 58 }$$

The limiting behaviours of $\Delta_{{\rm RS}}$ and $V+q$ cause the input MSE to vanish:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\MSEx}{\epsilon_{x}} \newcommand{\dAve}[1]{\left\langle \left\langle {#1} \right\rangle \right\rangle} \displaystyle V+q\to V+Q=\dAve{\MSEx} \to 0. \nonumber \end{align} \tag{ 59 }$$

Hence, the PR solution is successfully derived.

The local stability (42) of this PR solution should be checked. Inserting equations (50) and (53) into equation (42) results in the stability condition

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle (\alpha-\rho_0)(1-\rho_0)> (\rho-\rho_0)^2. \label{PRstability} \nonumber \end{align} \tag{ 60 }$$

Hence, the PR solution is stable as long as $\rho_0 \leqslant \rho \leqslant \alpha$ where the PR solution exists. Furthermore, this stability is considered to be rather 'robust'. It is physically reasonable that the RS solution is stable even for $\rho < \rho_0$ if ρ is sufficiently close to $\rho_0$ because the inequality (60) is safely satisfied even at $\rho=\rho_0$. This will be demonstrated in the phase diagram shown below. Note that equation (60) is just a necessary condition and not a sufficient one.

The stability of the RS solution at and around the PR solution has an algorithmic implication to the $\newcommand{\e}{{\rm e}} \ell_0$-minimisation approach [59–61]. Namely, local search algorithms such as the message-passing algorithm will not be degraded by the rugged energy landscape if the initial condition is sufficiently close to the PR solution and hence the predictions based on the RS computation will be precise.

We start from the noiseless limit $\sigma_{\xi}=0$. Figure 1 shows the plot of entropy $\newcommand{\e}{{\rm e}} \newcommand{\MSEp}{\epsilon_{y}} s(\MSEp)$ versus $\newcommand{\e}{{\rm e}} \newcommand{\MSEp}{\epsilon_{y}} \MSEp$ for several different values of ρ. Other parameters are common and are set as $\alpha=0.5$ and $\rho_0=0.2$. The drawn curves are based on the RS solution, and the RSB solution is only used to indicate the respective instability point. As seen from figure 1, we observe four different characteristic behaviours of $\newcommand{\e}{{\rm e}} \newcommand{\MSEp}{\epsilon_{y}} s(\MSEp)$. For small ρ, the AT instability occurs at small $\newcommand{\e}{{\rm e}} \newcommand{\MSEp}{\epsilon_{y}} \MSEp$, where the full-step RSB will be needed to correctly describe that region. This RS unstable region vanishes for larger ρ, defining a critical value $\rho_{{\rm AT}}(<\rho_0)$. In the region $\rho_{{\rm AT}} < \rho < \rho_0$, the RS solution is stable for the entire $\newcommand{\e}{{\rm e}} \newcommand{\MSEp}{\epsilon_{y}} \MSEp$ region having the nonnegative entropy $\newcommand{\e}{{\rm e}} \newcommand{\MSEp}{\epsilon_{y}} s(\MSEp)\geqslant 0$. A somewhat surprising fact in this region is that the entropy crisis (EC), $\newcommand{\e}{{\rm e}} \newcommand{\MSEp}{\epsilon_{y}} s(\MSEp)=0$, occurs at a finite critical value of $\mu_{\rm EC}(<\infty)$. As it approaches $\rho_0$, this critical value $\mu_{\rm EC}(\rho)$ diverges and the entropy curve is continuously connected to the PR solution in the region $\rho_0 \leqslant \rho$. For a wide range of $\rho(\geqslant \rho_0)$, the RS solution is again stable for the entire region of $\newcommand{\e}{{\rm e}} \newcommand{\MSEp}{\epsilon_{y}} \MSEp$. However, at larger values of ρ, there emerges a new phase transition, defining another critical value, $\rho_{\rm SP}$. For $\rho_{\rm SP} \leqslant \rho$, the PR solution is detached from the high temperature limit, and there are two different branches for the small $\newcommand{\e}{{\rm e}} \newcommand{\MSEp}{\epsilon_{y}} \MSEp$ or large μ region. Two critical temperature points are accordingly defined—the spinodal point T_SP at which the low-temperature branch connected to the PR solution vanishes, and the first-order transition point T_F at which g values of the two branches coincide. To locate the corresponding critical temperatures at a given $\rho(>\rho_{\rm SP})$, we plot g against $T=1/\mu$ in figure 2. As a reference, the output MSE $\newcommand{\e}{{\rm e}} \newcommand{\MSEp}{\epsilon_{y}} \MSEp$ is also plotted against T around the critical temperatures in the right panel. The presence of the first-order phase transition implies that the simple SA algorithm with a rapid annealing schedule will fail to find the PR solution in $\rho_{\rm SP}\leqslant \rho$. The system's dynamics goes along the branch connected to the high temperature limit and cannot move to the PR solution. Actually, this strongly affects the SA performance, as the values of input MSEs $\newcommand{\e}{{\rm e}} \newcommand{\MSEx}{\epsilon_{x}} \MSEx$ are quite different between the two branches, despite the output MSEs $\newcommand{\e}{{\rm e}} \newcommand{\MSEy}{\epsilon_{y}} \MSEy$ both being small, as shown in the right panel of figure 2. This will be demonstrated in section 4.2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Summarising the above findings, we draw a phase diagram in the ρ-T plane for $\alpha=0.5$ and $\rho_0=0.2$ in figure 3. Note that the entropy-crisis line $T_{\rm EC}(\rho)$ below the AT line $T_{\rm AT}(\rho)$ has no direct physical consequence because the RS solution is unreliable in that region. The exact entropy-crisis line would be derived by the full step RSB solution, but this is beyond the scope of the present paper.

Zoom In Zoom Out Reset image size

Figure 3 implies that finding the ground state is easy in the region $\rho_{\rm AT}< \rho < \rho_{\rm SP}$, while it is difficult in other regions: $\rho\leqslant \rho_{\rm AT}$ and $\rho_{\rm SP}\leqslant \rho$. We focus on how the easy region behaves when changing the external parameters α and $\rho_0$. We examine the parameters and find that the easy region shrinks as $\rho_0$ increases against a fixed α and finally vanishes at a certain critical value of $\rho_0$. As an example, the ρ-T phase diagram for $\alpha=0.5$ and $\rho_0=0.3$ is shown in figure 4. The spinodal and first-order transition lines cover the entire $\rho_0< \rho$ region and hence the easy region disappears, implying that local search algorithms cannot find the PR solution for any ρ in this case. This vanishing of the easy region thus defines the algorithmic limit. By searching all parameter regions of $\rho_0$ and α, we can draw a phase diagram of the algorithmic limit, which is shown in figure 5. The performance of our formulation is clearly better than the $\newcommand{\e}{{\rm e}} \ell_1$ relaxation [14] and is competitive with the Bayesian result [21]. This implies that the present formulation, which can be regarded as a MAP estimation in the Bayesian framework, does not significantly lose its reconstruction performance despite discarding the signal source information. This is encouraging the use of the present formulation in the context of signal recovery and is one of the main results of this paper.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A new behaviour specific to the noisy case is the presence of the RFOT transition for strong noises at middle values of ρ. As an example, the entropy curves for $\sigma_{\xi}^2=10$ are given in figure 6. Two critical ρ values, $\rho_{\rm AT}$ and $\rho_{\rm RFOT}$, accordingly emerge. For $\rho\leqslant \rho_{\rm AT}$, the AT instability first occurs as the temperature decreases. For $\rho_{\rm AT} < \rho \leqslant \rho_{\rm RFOT}$, the RFOT begins to emerge above the AT instability temperature. For larger ρ, the RS solution is accurate for all temperature regions. The RS EC occurs at finite T in that region, which is somewhat similar to the Ising perceptron problem [37, 38]. Note that $\rho_{\rm RFOT}$ does not exist if the noise is sufficiently weak.

Zoom In Zoom Out Reset image size

Summarising the above findings, we show phase diagrams for three noise strengths, $\sigma_{\xi}^2=0.001, ~0.1$ and 10 in figure 7. The first-order transition in the noiseless case is quite fragile and disappears for very weak noise as seen in the left panel of figure 7. We attempted to capture the critical values of $\sigma_{\xi}$ for the disappearance, but it proved numerically difficult and we did not pursue this point. As the noise increases, the phase boundary's shape expand more and more from the noiseless limit and the RSB region, $\rho<\rho_{\rm RFOT}$, seems to grow. We consider whether this RSB region will cover the low-temperature region completely as the noise is very large. To answer this, we tested the no signal case $\rho_0=0$ and $\sigma_{\xi}^{2}=1$ and observed that $\rho_{\rm RFOT}$ takes a value similar to the one in the right panel of figure 7. Hence, an RS region exists even in the strong noise case, which is consistent with our previous analysis in the data compression context [24].

Zoom In Zoom Out Reset image size

The free-entropy values evaluated by numerical simulations and the extrapolation to the infinite size limit of the noiseless case $\sigma_{\xi}^2=0$ are presented in figure 8. The extrapolation lines result from linear regression using an asymptotic form $g\approx a+bN^{-1}+cN^{-1}\log N^{-1}$. The regression is conducted by applying the least squares method as follows:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \min_{a,b,c}\sum_N \left(a+b\frac{1}{N}+c\frac{1}{N}\log\frac{1}{N}-g(N)\right)^2. \label{regression} \nonumber \end{align} \tag{ 65 }$$

Zoom In Zoom Out Reset image size

This asymptotic form is based on Stirling's formula and is exact at $\mu=0$, which motivates us to use the form even when $\mu\not=0$. The same asymptotic form is used for obtaining the extrapolated values of the output MSE $\newcommand{\e}{{\rm e}} \newcommand{\MSEp}{\epsilon_{y}} \MSEp$ and the entropy s. Using these values for the limit $N \to \infty$, we present the curves $g(T)$ and $\newcommand{\e}{{\rm e}} \newcommand{\MSEp}{\epsilon_{y}} s(\MSEp)$ in figure 9. The lines represent the RS analytical results. The circles represent the extrapolated values obtained from the numerical results. The extrapolated values show fairly good agreement with the RS analytical ones, justifying our analytical results. For the case of $\rho=0.1$, the AT instability occurs at $T\approx 0.04$, but even below this temperature, the agreement between the RS analytical result and the numerical one is fairly good, suggesting a weak RSB effect on s and $\newcommand{\e}{{\rm e}} \newcommand{\MSEp}{\epsilon_{y}} \MSEp$. This is, however, not the case for the input MSE $\newcommand{\e}{{\rm e}} \newcommand{\MSEx}{\epsilon_{x}} \MSEx$, as demonstrated in section 4.2 below.

Zoom In Zoom Out Reset image size

A similar analysis for the case of strong noise ($\sigma_{\xi}^2=10$) is performed and the results are shown in figures 10 and 11. These figures again demonstrate good agreement between the RS and numerical results as long as the RSB does not occur. Below the RSB transition point, we observe a deviation between them, as shown in the left upper panel of figure 11. In such a situation, generally speaking, the RS entropy curve can be regarded as an upper bound of the entropy values [64], although a meaningful difference is not observed in the right upper panel of figure 11. Again, the RSB effect on s and $\newcommand{\e}{{\rm e}} \newcommand{\MSEp}{\epsilon_{y}} \MSEp$ appears to be weak.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have a noteworthy remark to make on the EC phenomenon for $\rho=0.3$. We stress that this EC phenomenon can be described in the RS level and occurs at a finite temperature. This may be somewhat surprising for readers familiar with other models exhibiting similar EC phenomena, because in most of such systems the 1RSB treatment is needed to describe the EC phenomenon. We note that the energy levels of the present system around the ground state can be very dense and the energy gap between the ground and excited states can be extremely small, which can be argued by the fact that the PR solution in the noiseless case has numerous degeneracies for $\rho >\rho_0$. This gap is supposed to vanish in the $N\to \infty$ limit, presumably enabling the RS EC phenomenon to appear at a finite μ. The agreement between the RS and numerical results strongly argues in favour of this description. Note that the EC phenomenon at small values of ρ is in a different situation and its RS description is not accurate. This is because the AT instability occurs at higher temperatures in that region and hence the full step RSB treatment is needed. This is in contrast to the large-ρ region in which no instability occurs at higher temperatures than T_EC.

Noiseless case.

Let us begin by showing the results for the noiseless case. The MSEs at $\alpha=0.5$ and $\rho_0=0.2$ are plotted against T in figure 12. The assumed values of ρ are $\rho=0.02$, $\rho=0.3$, and $\rho=0.45$ for the left, middle, and right panels, respectively. The numerical results agree well with the black solid line representing the RS solution connected to the high temperature limit in all cases. The middle panels show the successful region for finding the PR solution, $\rho_{0} < \rho < \rho_{\rm SP}$, and the vanishing $\newcommand{\e}{{\rm e}} \newcommand{\MSEx}{\epsilon_{x}} \MSEx$ means that we actually find the PR solution. The right panels are in $\rho_{\rm SP} < \rho$, meaning that the search is trapped in the metastable state connected to high temperatures. The SA result follows the high-temperature branch and cannot reach the low-temperature one denoted by the blue dashed dotted line. Overall, the SA experiments demonstrate that our theoretical predictions are very precise, and the presence of phase transitions strongly degrades the SA's performance in finding the minimum-MSE configuration.

Zoom In Zoom Out Reset image size

Noisy case.

Next, we show the results of the noisy case. The SA results for the strong noise case $\sigma_{\xi}^2=10$ at $\alpha=0.5$ and $\rho_0=0.2$ are given in figure 13. As seen from the figure, the MSEs show good agreement to the analytical curve (black solid line) up to the transition points for the left and middle panels. An exceptional deviation is observed at low temperatures in the upper right panel, but this is considered to be a finite-size effect because the range of $\newcommand{\e}{{\rm e}} \newcommand{\MSEp}{\epsilon_{y}} \MSEp$ in this region is very small and supposedly unreachable by $N\approx 100$ systems. Hence, these behaviours are very consistent with the analytical predictions that the system's dynamic behaviour is affected by the RSB transitions and ceases to follow the equilibrium state.

Zoom In Zoom Out Reset image size

The effect of the RSB on the reconstruction performance becomes much clearer by examining the achievable-limit values of the MSEs for a moderate noise case. Figure 14 shows the plots against ρ of the limit values obtained at very low temperatures by the rapid SA with $\tau=5$ for $\sigma_{\xi}^2=0.1$. The $\newcommand{\e}{{\rm e}} \newcommand{\MSEx}{\epsilon_{x}} \MSEx$ values at middle ρ values are clearly dominated by the ones at the AT transition points rather than the ones at the EC points, implying that the system's search is trapped by local minima emerging at the transition points. An interesting outcome of this phenomenon is a better reconstruction of the planted signal $\newcommand{\V}[1]{\boldsymbol{#1}} \V{x}_0$. As seen from figure 14, the $\newcommand{\e}{{\rm e}} \newcommand{\MSEx}{\epsilon_{x}} \MSEx$ values at the AT points are lower than the ones at the RS EC points, implying that the reconstruction performance of the local minima induced by the RSB is better than that of the minimum-$\newcommand{\e}{{\rm e}} \newcommand{\MSEp}{\epsilon_{y}} \MSEp$ configuration $\newcommand{\V}[1]{\boldsymbol{#1}} \hat{\V{c}}$ by solving equation (2). This means that the generalisation capability of the rapid SA is no worse than exactly solving equation (2) in this case because the input MSE $\newcommand{\e}{{\rm e}} \newcommand{\MSEx}{\epsilon_{x}} \MSEx$ is proportional to the generalisation error when each component of the design matrix and the noise is i.i.d. from the zero-mean Gaussian. This encourages the use of the presented formulation and algorithm for practical purposes, as reported in [27].

Zoom In Zoom Out Reset image size