L0 regularization-based compressed sensing with quantum-classical hybrid approach

L0-regularization-based compressed sensing (L0-RBCS) has the potential to outperform L1-regularization-based compressed sensing (L1-RBCS), but the optimization in L0-RBCS is difficult because it is a combinatorial optimization problem. To perform optimization in L0-RBCS, we propose a quantum-classical hybrid system consisting of a quantum machine and a classical digital processor. The coherent Ising machine (CIM) is a suitable quantum machine for this system because this optimization problem can only be solved with a densely connected network. To evaluate the performance of the CIM-classical hybrid system theoretically, a truncated Wigner stochastic differential equation (W-SDE) is introduced as a model for the network of degenerate optical parametric oscillators, and macroscopic equations are derived by applying statistical mechanics to the W-SDE. We show that the system performance in principle approaches the theoretical limit of compressed sensing and this hybrid system may exceed the estimation accuracy of L1-RBCS in actual situations, such as in magnetic resonance imaging data analysis.


Introduction
Quantum machines have attracted significant interest because of their potential to overcome the difficulty of solving large-scale combinatorial optimization problems.Many quantum machines, such as the quantum annealers (QA) of D-Wave systems [1], the quantum approximate optimization algorithm (QAOA) [2,3], quantum bifurcation machines [4,5,6] and coherent Ising machine (CIM) [7,8,9,10,11,12], have been proposed in the past decade.Other examples include classical annealers, which have been implemented in electromechanical resonators [13], nanomagnet arrays [14], electronic oscillators [15], silicon photonic weight banks [15], complementary metal?oxide?semiconductor static random access memory circuits [16,17,18], and fieldprogrammable gate arrays (FPGAs) [19,20].Interest has been centered on implementing quantum machines and understanding their behavior, whereas there have been few practical applications [21,22,23].To open the door to practical use of quantum machines, we show that they can be used for implementing compressed sensing (CS).Furthermore, we demonstrate, using non-equilibrium statistical mechanics [24], that the system performance in principle approaches the theoretical limit of CS.
On the other hand, L0-regularization-based CS (L0-RBCS) can be formulated with the following L0 norm instead of the L1 norm [43]: L0-RBCS, as defined in Eq. (2), can be equivalently reformulated as a two-fold optimization problem [43,44]: Here the vector r is the value of the N-dimensional source signal and the element r i in r represents the real-number value of the i-th element in the source signal.The vector σ is called a support vector, which represents the places of the non-zero elements in the Ndimensional source signal.The element σ i in σ takes either 0 or 1 to indicate whether the i-th element in the source signal is zero or non-zero.The symbol • denotes the Hadamard product.From the elementwise representation of Eq. ( 3), the Hamiltonian (or cost function) of L0-RBCS can be written as where A µ i is an element in an M-by-N observation matrix A, and y µ is an element in an M-dimensional observation signal.
The minimization of H with respect to r under the condition that σ is fixed is the same as the problem of solving a system of simultaneous linear equations that gives the minimum point of the quadratic potential for r.On the other hand, the minimization of H with respect to σ under the condition that r is fixed is the same as the problem of finding the ground state of an Ising Hamiltonian, where − M µ=1 A µ i A µ j r i r j can be considered to be the mutual interaction between σ i and σ j .
It has been suggested that L0-RBCS has the potential to outperform L1-RBCS, because L1 regularization imposes a shrinkage on variables over a threshold (soft-thresholding) but L0 regularization does not impose such a shrinkage (hardthresholding) [43].However, the optimization of the support vector is a combinatorial optimization problem, which can be mapped into an Ising model, as mentioned above.In this optimization problem, there are a lot of meta-stable states because the effective interaction − M µ=1 A µ i A µ j r i r j induces frustration in the Ising model in the minimization of H with respect to σ under the condition that r is fixed.Thus, it is difficult to solve this kind of problem.Because of this difficulty, only a few approximation algorithms have been proposed and they only work under special conditions [45,46,47].
In this paper, to overcome the difficulty of optimizing the support vector σ, we focus on quantum machines.We propose a quantum-classical hybrid system composed of a quantum machine and a classical digital processor (CDP) (Fig. 1).This system solves the two-fold optimization problem by alternately performing two minimization processes; (i) the quantum machine optimizes σ to minimize H under the condition that r is fixed, and (ii) the CDP optimizes r to minimize H under the condition that σ is fixed.If the quantum machine can find the ground state of H under the condition that r is fixed, the quantum-classical hybrid system is expected to outperform L1-RBCS.
Several quantum machines can potentially be used for optimizing σ, such as QA [1], QAOA [2,3], CIM [7,8,9,10,11,12], and so on.As defined in Eq. ( 4), the number of non-zero connections from one spin to the others is O(N); thus, it is necessary to form a densely connected network on a quantum machine in order to optimize σ.A comparison of these candidates reveals that a measurement-feedback (MFB) CIM is one of most suitable machines for this purpose.In fact, an MFB-CIM can construct any densely connected network composed of degenerate optical parametric oscillators (OPOs) because it uses a time-division multiplexing scheme and MFB [9,10].In contrast, QA and almost all other machines can only support local graphs, including chimera graphs, and thus, a densely-connected network for optimizing σ has to be embedded in a fixed hardware local graph by using the Lechner-Hauke-Zoller scheme, which requires additional physical spins [48,49].Furthermore, it was reported [50] that an MFB-CIM experimentally outperformed QA on two problem sets, i.e., a fully connected Sherrington-Kirkpatrick model [51] and dense graph MAX-CUT.In contrast to QA having an exponential computation time proportional to exp (O (N 2 )), a CIM has an exponential computational time proportional to exp O √ N , where N is the problem size [50].
Here, we evaluate the performance of a quantum-classical hybrid system composed of an MFB-CIM and CDP (Fig. 2).We derive a truncated Wigner stochastic differential equation (W-SDE) from the master equation for the density operator of a network consisting of OPOs.Then, we develop a statistical mechanics method based on self-consistent signal-to-noise analysis (SCSNA) [52,24,53] and derive a macroscopic equation (ME) for the whole system [54,55,56].We show that the performance of the hybrid system approaches the theoretical limit of L0-minimization-based CS [44] and the hybrid system may exceed the estimation accuracy of L1-RBCS in actual situations, such as MRI data analysis.

Configuration of CIM-CDP hybrid system
The CIM-CDP hybrid system (Fig. 2) executes the L0-RBCS defined as Eqs.( 3) and (4).This system optimizes by alternately performing the following two minimization processes.The CIM optimizes σ to minimize H under the condition that r is fixed and forwards σ to the CDP.The CDP then optimizes r to minimize H under the condition that σ is fixed and then forwards r to the CIM.
At a stationary point r and σ that satisfies ∂H ∂σ i = 0 and ∂H ∂r i = 0, the following equations hold: where h i is the local field and H(X) is the Heaviside step function taking 0 for X ≤ 0 or +1 for X > 0.
In this paper, we assume that M µ=1 (A µ i ) 2 = 1 is satisfied.This assumption does not lose any generality because it is possible to normalize the observation matrix A to satisfy M µ=1 (A µ i ) 2 = 1 for any case.Under this assumption, r i = σ i h i is satisfied in Eq. ( 6), and according to the Maxwell rule [57], a stationary point of σ i = H (σ i h 2 i − λ) obtained by substituting r i = σ i h i into Eq.( 5) can be dertermined as follows, where one of two different functions F + (h) and F ± (h) is chosen depending on the source signal.F χ (h) is the identity function if the source signal is non-negative, and F χ (h) is the absolute value function if the source signal is signed.In the presence of noise, this conversion increases the threshold-to-noise ratio, which allows the low threshold to work as a sparse bias, as shown in the experiment below.
The CIM estimates the support vector σ, i.e. the places of the non-zero elements in the source signal.According to Eq. ( 8), the optical field injecting to the target (i-th) OPO pulse is set as where h i is the local field explained below, K is the gain of the feedback circuit, and η is the threshold.η is related to the regularization parameter λ in Eq. ( 4) by η = √ 2λ, as shown in Eq. ( 8).We use one of two different functions F + (h) and F ± (h) depending on the source signal.F + (h) is the identity function: it is used as a non-negative source signal.F ± (h) is the absolute value function: it is used as a signed source signal.
The local field for the support estimation in the CIM, which corresponds to Eq. ( 7), is set as where X j is the in-phase amplitude (generalized coordinate) of the j-th OPO pulse measured by a homodyne detector, H(X j ) is the binarized in-phase amplitude of the j-th OPO pulse through the Heaviside step function, and r j is a solution for the signal value given by the CDP.During the support estimation on the CIM, all r j are fixed.The first term of Eq. ( 10) is the mutual interaction term, while the second term is the Zeeman term.The CDP estimates r, i.e. the values of the non-zero elements in the source signal.In accordance with the simultaneous equations (6) satisfied by the stationary point that minimizes H with respect to r, the CDP solves the following simultaneous equations: Here, h i in Eq. ( 12) is the local field for the signal estimation in the CDP, which corresponds to Eq. (7).Equation ( 12) for the signal estimation in the CDP is the same as Eq. ( 10) for the support estimation in the CIM.During the signal estimation in the CDP, all H(X j ) are fixed.The solution of the simultaneous equations (Eq.( 11)) is Algorithm 1 is an outline of the alternating minimization processes.In this algorithm, to make the basin of attraction wider, we heuristically introduce a linear threshold reduction whereby the threshold η is linearly lowered from η init to η end as the alternating minimization proceeds.

Derivation of W-SDE for CIM
As shown in Fig. 2, the pump pulses are injected into the main ring cavity through a second harmonic generation (SHG) crystal.A periodically poled lithium niobate (PPLN) waveguide is a highly efficient nonlinear medium for optical parametric oscillation.Suppose that the amplitude of the pump field injected into the main cavity is ǫ and the parametric coupling constant of the PPLN waveguide between the signal field and the pump field is κ.Then, the pumping Hamiltonian is Ĥ1 = i ǫ(â † p − âp ) and the parametric interaction Hamiltonian is Ĥ2 = i κ/2(â †2 s âp − â † p â2 s ).Here, âp and âs are the annihilation operators for the intra-cavity pump and signal fields.If the round-trip time of the ring cavity is correctly adjusted to N times the pump pulse interval, N independent and identical OPO pulses are simultaneously generated inside the cavity.The photon annihilation and creation operators for the j-th OPO signal pulse are denoted by âj and â † j .The intra-cavity pump field and signal field have loss rates γ p and γ s , respectively.If γ p ≫ γ s , the pump field can be eliminated by invoking the slaving principle: the following master equation of the density operator for a solitary j-th OPO signal pulse is obtained by adiabatic elimination of the pump mode [58,59], where S = ǫκ/γ p and B = κ 2 /(2γ p ) are the linear parametric gain coefficient and two photon absorption (or back conversion) rate, respectively.[x, ŷ] denotes the bosonic commutator.
Next, let us examine the measurement-feedback circuit shown in Fig. 2. The circuit is connected to the main cavity by extraction and injection couplers with reflection coefficients R ex = j ex ∆t and R in = j in ∆t, where j ex and j in are coarsegrained out-coupling and in-coupling constants and ∆t is the cavity round trip time.When B/γ s << 1 and vacuum fluctuations are incident on the open ports of the extraction and injection couplers, the measurement-feedback circuit can be described with a Gaussian quantum model [60,61].The master equation consists of a linear loss term, measurement-induced state reduction term, and coherent feedback signal injection term (see Eqs. (12) (13) (14) in ref. [61]).
The Fokker-Planck equation is derived using the Wigner W (α) representation of the density operator ρ in the master equations, and we arrive at the following truncated Wigner stochastic differential equation (W-SDE) by applying Ito's rule [62,61], where j = j ex + j in , α i is the complex Wigner amplitude, and υ i is the c-number noise amplitude satisfying υ i (t) = 0, υ * i (t)υ j (t ′ ) = 2δ ij δ(t − t ′ ).Then, by introducing a saturation parameter A s = 2γ p (γ s + j)/κ 2 and applying the following scale transformation: α i /A s = c i + is i , t(γ s + j) = t, p = S/(γ s + j) and Kj in /A s (γ s + j) = K, we obtain where c i and s i are the in-phase and quadrature-phase normalized amplitudes of the i-th OPO pulse.The second term of the R.H.S. in the upper equation of Eq. ( 15) is the optical injection field, which has only an in-phase component.p is the normalized pump rate.p = 1 corresponds to the oscillation threshold of a solitary OPO without mutual coupling.If p is above the oscillation threshold (p > 1), each of the OPO pulses is either in the 0-phase state or π-phase state.The 0-phase of an OPO pulse is assigned to an Ising-spin up-state, while the π-phase is assigned to the down-state.The last terms of the upper and lower equations express the vacuum fluctuations injected from external reservoirs and the pump fluctuations coupled to the OPO system via gain saturation.W i,1 and W i,2 are independent real Gaussian noise processes satisfying W i,k (t) = 0, W i,k (t)W j,l (t ′ ) = δ ij δ kl δ(t − t ′ ).The saturation parameter A s determines the nonlinear increase (abrupt jump) of the photon number at the OPO threshold.Finally, we should mention that a more general quantum model of MFB-CIM without a Gaussian approximation was derived for both discrete time models [63] and continuous time models [64].

Precondition for applying statistical mechanics
To solve the W-SDE (15) and the simultaneous equations (11) using statistical mechanics methods, we introduce the following observation model in which the values of all variables are randomly chosen, where Here, we introduce a measurement compression rate α, defined as α = M/N.

Each element of [A
T randomly and independently takes 1 or 0 with probability a and 1 − a, respectively.Here, we introduce a sparseness a defined as the probability of non-zero elements in the source signal. Each T is also an independent and identically distributed value generated from some probability distribution g(x).To verify the system performance, we use the following probability density functions for generating the source signal: 3).To verify the invariance of our results relative to the type of probability distribution of the source signals, we used two different probability distributions in each of the non-negative and signed cases.Gaussian and bilateral Gamma distributions were used to generate the signed source signals.On the other hand, half-Gaussian and Gamma distributions were used to generate the nonnegative source signals.The second moments of the half-Gaussian and Gaussian were set to x 2 x = 1.The shape and scale parameters of the Gamma and bilateral Gamma were set to k = 2 and θ = 0.4; thus, the second moment of both distributions was x 2 x = 0.96.The figures in the main text show results for source signals generated from the half-Gaussian and Gaussian, while the supplementary figures show results for source signals from the Gamma and bilateral Gamma.

Outline of derivation of MEs for whole hybrid system
Here, we summarize the derivation of the MEs by solving the W-SDE (15) and the simultaneous equations (11) under the precondition described in section 2.3.A detailed derivation of the MEs is provided in sections 2.5 and 2.6.
As described in section 2.1, the CIM and CDP share the same local field (Eqs.( 10) and ( 12)).Therefore, the CIM and CDP can be unified into a single mean field system in the steady state.By substituting the observation model ( 16), the shared local field can be rewritten as where X j is replaced with c j , which is the real part of the normalized complex Wigner amplitude c j + is j explained in section 2.2.The W-SDE ( 15) of the i-th OPO implies that H(c i ) is a stochastic variable depending on the local field h i and time t in the steady state.Here, we introduce the following formal time-dependent stochastic transfer function from h i to H(c i ) [24,53]: Substituting X(h i , t) into Eq.( 11) and noting that 1 , we can write a formal time-dependent stochastic transfer function from h i to r i as It is difficult to specify such a transfer function concretely, but it is possible to introduce one formally.As a premise that this transfer function holds, we assume that the microscopic memory effect can be neglected in the steady state [65].
The local field can be defined self-consistently through the formal transfer function G, as follows: Following the recipe of the SCSNA [52,24,53], the local field h i can be separated into a pure local field and an Onsager reaction term (ORT) [24,53]: The pure local field hi of the i-th OPO pulse is statistically independent of hj of the j-th OPO pulse when i = j.The ORT is an effective self-feedback via other OPO pulses.By substituting Eq. ( 19) into Eq.( 15), we obtain Eq. ( 28).The W-SDE (28) can be regarded as describing N independent one-body OPO pulses, because each pure local field hi is statistically independent of the other.In other words, the W-SDE (28) can be regarded as N independent equations.
The i-th independent equation in the W-SDE (28) implies that H(c i ) is a stochastic variable depending on the pure local field hi and time t in the steady state.Thus, X(h i , t) and G(h i , t) can be redefined with a pure local field hi as follows: where X and G are formal time-dependent stochastic transfer functions from the pure local field hi to H(c i ) and r i , respectively.It is difficult to specify these transfer functions concretely, but it is possible to calculate the expectations, such as X( hi , t) SDE , G( hi , t)

SDE
and their higher order moments.X( h, t) SDE and G( hi , t) SDE respectively denote the conditional expectations of X( h, t) and G( hi , t) given the pure local field h, and they can be approximately calculated with the W-SDE (29), as described in section 2.5.
The macroscopic parameters sought when separating the local field into the pure local field and the ORT using SCSNA are Q, R and U, which are called the overlap, the mean square magnetization, and the susceptibility, respectively.R, Q ,and U are defined as the site averages of G( hi , t), its square, and its sensitivity to the local field (see Eqs. (40) (41) (42)).Because the pure local fields hi are statistically independent of each other, the site averages in R, Q and U can be replaced with the averages of G( h, t) SDE and G( h, t) 2

SDE
with respect to h.Finally, the following MEs are obtained: where • x,ξ denotes the average with respect to x and ξ, and Ξ c (z, xξ) and Ξ s (z, xξ) can be determined self-consistently from the following equations, A s is the saturation parameter, which diverges in the infinite limit of the amplitude of the injected pump field ǫ → +∞ (see section 2.2).In the limit A s → +∞, we obtain the following simplified MEs, Here, X(h p , h m ) is an effective output function obtained from the Maxwell rule [57]:

Mean-field behavior of OPO pulses and CDP
In this section, we approximately calculate the conditional expectations of X( h, t), G( h, t) and G( h, t) 2 given the pure local field h, which are denoted by X( h, t) Under the premise that the pure local field and ORT in Eq. ( 19) can be separated by following the recipe of the SCSNA [52,24,53], substituting Eq. ( 19) into Eq.( 11) and because 1 M M µ=1 (A µ i ) 2 = 1, we can write r i as Furthermore, substituting Eqs. ( 19) and ( 27) into the W-SDE (15) gives: Equation ( 28) of the i-th OPO pulse only depends on the pure local fields hi , which are statistically independent of each other in the steady state.The W-SDE (28) can be regarded as describing N independent one-body OPO pulses in the steady state.In other words, the W-SDE (28) can be regarded as N independent equations in the steady state, and it is not necessary to solve the W-SDE (28) simultaneously.Since the steady-state solution of Eq. ( 28) depends only on the value of the pure local field, the site index i in Eq. ( 28) can be deleted.It is difficult to solve Eq. ( 28) analytically even after the N-body system has been reduced to a one-body system.To obtain a mathematically tractable form, we replace the state variables in the secondorder coefficient of the Kramers-Moyal expansion [62] (representing the power of the quantum noise) with the average values of these state variables [54]: From Eq. ( 29), we can derive the following equations to determine the approximate value of X( h, t)
where V (c, s) is the potential appearing in the CIM-ferromagnetic and the CIM-finite loading Hopfield models [54].Ξ c and Ξ s are parameters for calculating X( h, t) Ξ c and Ξ s are equal to c 2 and s 2 , and by giving h and Γ, they can be self-consistently determined from the above equation.
Similarly, from Eq. ( 27), G( h, t) SDE and G( h, t) 2 SDE can be obtained as follows: 2.6.Details of SCSNA for the whole hybrid system Under the precondition described in section 2.3, we separate the local field into the pure local field and the ORT (Eq.( 19)) with SCSNA [52,24,53,55,56] and reduce the N-body system composed of N mutually coupled OPO pulses to an effective one-body system.After that, we derive the ME for the whole hybrid system.Let us start by introducing the following parameters.
Below, we assume that ) is satisfied, because, under the precondition, the correlation between A µ j and G(h j , t) − ξ j x j is O(1/ √ N) for any µ if the reconstruction succeeds.
Substituting Eq. (30) into Eq.(18) gives where the first term is the cross-talk noise part, and the second term is introduced to subtract the direct self-coupling term from the local field h i because Eq. ( 18) does not contain the direct self-coupling.
Next, we split the local field into a signal term, independent Gaussian noise, and the ORT.g µ , as defined in Eq. ( 30), recursively contains A µ j g µ in G(h j , t), so it is a factor causing correlation between OPO pulses.Because g µ = O(1/ √ N), we perform the following expansion on Eq. ( 30): where h is the cavity field [66] and U (µ) is a macroscopic parameter called the susceptibility, which are given by The cavity field h (µ) j does not contain A µ j g µ , so G(h ) is uncorrelated with A µ j and U (µ) is also uncorrelated with A µ j .The terms that cause the correlation between OPO pulses are extracted by performing a first-order Taylor expansion around g µ = 0 and these extracted terms form the third term on the right side of Eq. (32).
From Eq. ( 32), we redefine g µ on the basis of the cavity fields h The terms causing the correlation between OPO pulses in Eq. ( 30) are converted into the scale coefficient of g µ defined by the cavity fields in Eq. (35).Substituting Eq. ( 35) into the crosstalk noise in Eq. ( 31), we split up the local field into three terms, as follows: where the first term is the signal term, Z i in the second term is Gaussian random noise defined by Eq. ( 37), and the third term is the self-coupling term.Here, µ) .These three terms are obtained under the conditions A µ i = 0 and A µ i A ν j = δ ij δ µν , and G(h i , t) and U (µ) are uncorrelated with A µ i .Moreover, the third term is obtained under the approximation G(h From the central limit theorem, Z i becomes Gaussian random noise in the thermodynamic limit.The average of Z i and the covariance between Z i and Z j are where R (µ) and Q (µ) are macroscopic parameters that are respectively called the overlap and the mean square magnetization and are given by x j ξ j G(h Because Z i is statistically independent of Z j when i = j, the first and second terms in Eq. ( 36) are statistically independent of those of other sites.The third term is the difference between the self-coupling term in the crosstalk noise rescaled by the extracted correlation factors and the direct self-coupling term in Eq. ( 31), and it represents selffeedback via other OPO pulses.Therefore, the first and second terms are the pure local field and the third term is the ORT.By comparing Eqs. ( 19) and (36), hi and Γ are determined as follows: As explained in section 2.5, substituting Eq. ( 19) into the W-SDE (15) reduces the N-body system to an effective one-body system.The W-SDE (28) can be regarded as N independent equations.The i-th independent equation in the W-SDE (28) implies that H(c i ) is a stochastic variable depending on the pure local field hi and time t in the steady state.Thus, X(h i , t) and G(h i , t) can be redefined as X( hi , t) and G( hi , t), as shown in Eq. (20).Through the manipulations in Eqs.(32) and (35), the pure local field and the ORT are defined on the basis of the cavity field.Therefore, in the thermodynamic limit (N → ∞), the cavity field can be consistently replaced with the pure local field and the ORT, and G(h ) (39) (34) can be safely replaced with G( hi , t).As a result of this replacement, the cavity indexes (µ) of R (µ) , Q (µ) , and U (µ) become negligible, and these macroscopic parameters are redefined as follows: where U expresses the average sensitivity of G( hi , t) to the bare local field h i using the chain rule because of the definition of U (µ) in Eq. ( 34).
Because the pure local fields are independent of each other, the site averages in Eqs.(40)(41)(42) can be replaced with the averages of G( h, t) SDE and G( h, t) 2 SDE with respect to the Gaussian random noise Z and the source signal ξx.The replacement for U can be achieved by integration by parts.Finally, we obtain the MEs (21)(22)(23) for finite A s and the MEs (24)(25)(26) for infinite A s in section 2.4.

Perturbation expansion for ME in the limit
By introducing a new macroscopic parameter W defined by W = Q − 2R, when there is no observation noise, i.e. β = 0, we can rewrite the ME in the limit A s → ∞ as Here, we put 2η/ 1 + 1/F χ 1 + a α U = ζ 2 .In the limit η → +0, i.e. ζ → 0, the ME (43) has a solution W = − x 2 x corresponding to perfect reconstruction.We assume that the above ME has the following solution when ζ ≪ 1.
x + ζ 2 w.Substituting this into the ME (43) and expanding around ζ = 0, we obtain the following relation independent of the probability distribution of x, g(x), if g(0) and g ′ (0) are finite.
This equation suggests that the solution w = 0, i.e. the perfect reconstruction solution, is stable when a < α, neutral when a = α, and unstable when a > α.Thus, when there is no observation noise, in the infinite limit of the amplitude of the injected pump field (i.e.A 2 s → ∞) and in the infinitesimal limit of η, the phase transition line of the perfect reconstruction solution becomes a = α independent of g(x).

ME of LASSO
Under the observation model ( 16), the update rule of the LASSO is given by where T χ,η (h j ) is a soft-thresholding function with threshold η, defined as We use two different functions T +,η (h) and T ±,η (h) depending on the source signal.T +,η (h) is for non-negative source signals, and T ±,η (h) is for signed source signals.
Following the same recipe of the SCSNA, we obtain the following MEs, where the pure local field of the LASSO becomes Tχ,η ( h) is an effective output function into which the ORT is renormalized: .
T+,η ( h) is for non-negative source signals, and T±,η ( h) is for signed source signals.

Root-mean-square error
The numerical experiments used the root-mean-square error (RMSE) as a measure of estimation accuracy.The RMSE of CIM L0-RBCS and LASSO is where R and Q are the overlap and mean square magnetization defined above, a is sparseness, and a x 2 x is the second moment of the source signal.The RMSE is zero if CIM L0-RBCS / LASSO perfectly reconstructs the source signal.

Accuracy of MEs and comparison with LASSO when β = 0
First, to confirm the accuracy of the MEs, we compared solutions to the MEs with those given by Algorithm 1. Figures 4a and 4b show the root-mean-square errors (RMSEs) (defined in section 2.9) of the solutions to the MEs with A 2 s = 250 (Eqs.(21)(22)( 23)) and those in the limit A 2 s → ∞ (Eqs.(24)(25)( 26)) for various values of the threshold η and compression rate α (red and green solid lines).The figures also indicate the RMSEs of solutions obtained using Algorithm 1 with A 2 s = 250 and A 2 s = 10 7 (circles with error bars).Note that A 2 s = 10 7 is on the same order as A 2 s in real experimental CIMs.The results in Fig. 4 are for the case when there is no observation noise (i.e.β = 0) and the source signals are from a half-Gaussian (+) or Gaussian (±).To confirm if Algorithm 1 has solutions corresponding to the near-zero RMSE states obtained by the MEs, r was initialized to the true signal value, i.e. x • ξ, in the alternating minimization process in Algorithm 1.However, even in this situation, the c-amplitude was always initialized to c = 0 in the initial stage of the support estimation.
As results obtained by MEs (24)(25) (26), in the case of the half-Gaussian (+), two macroscopic states with non-zero RMSE (red solid lines) and near-zero RMSE (green solid lines) coexist, as in a CIM-implemented CDMA multiuser detector [67,56].On the other hand, in the case of the Gaussian (±), we only found a single macroscopic state with near-zero RMSE (red solid lines).In both the half-Gaussian case (+) and Gaussian case (±), the near-zero-RMSE states of the MEs (24)(25)(26) (red solid lines in Fig. 4b) matched the numerical results of Algorithm 1 with A 2 s = 10 7 (circles with error bars in Fig. 4b), and the phase transition points given by the MEs (24)(25) (26) coincided with those of Algorithm 1.Under these conditions, as η was lowered to 0.01, RMSE decreased monotonically and the phase transition point a c from the near-zero-RMSE state grew monotonically.On the other hand, the theoretical results obtained from the MEs (21)(22)(23) with A 2 s = 250 (red solid lines on the left of Fig. 4a) were in good agreement with the numerical results of Algorithm 1 with A 2 s = 250 (circles with error bars on the left of Fig. 4a) in the half-Gaussian case (+), whereas the phase transition points given by the MEs (21) (22) (23) became lower than those of Algorithm 1 when η = 0.01 in the Gaussian case (±), as shown on the right of Fig. 4a.
Furthermore, to compare the abilities of CIM L0-RBCS and LASSO, we computed the RMSE profiles of LASSO using the MEs (44)(45) (46) with the same threshold value as CIM L0-RBCS; these profiles are superimposed upon Fig. 4 (blue solid lines).The RMSEs of CIM L0-RBCS in the limit A 2 s → ∞ (red solid lines in Fig. 4b) were lower than those of LASSO (blue solid lines) at the same compression rate α and sparseness a, and the first-order phase transition points of CIM L0-RBCS were higher than those of LASSO.On the other hand, the RMSEs of CIM L0-RBCS with A 2 s = 250 (red solid lines and circles with error bars in Fig. 4a) were lower than those of LASSO (blue solid lines) when η = 0.1 and 0.05, but the theoretical RMSEs of CIM L0-RBCS became higher than those of LASSO when η = 0.01.
We numerically checked that qualitatively the same results were obtained even in the case of source signals from the Gamma (+) and the bilateral Gamma (±) (See Supplementary Fig. 1).

Phase diagrams of CIM L0-RBCS and LASSO when β = 0
We drew phase diagrams of CIM L0-RBCS for various values of the threshold η when there was no observation noise (i.e.β = 0).Figure 5a and Supplementary Fig. 2 show first-order phase transition lines from the near-zero-RMSE state (red lines) in the half-Gaussian case (+) and Gaussian case (±).The phase transition lines in Fig. 5a are for the limit A 2 s → ∞, while the ones in Supplementary Fig. 2 are for A 2 s = 250.As demonstrated in Fig. 5a, in the limit A 2 s → ∞, the phase transition lines from the near-zero-RMSE state in CIM L0-RBCS become asymptotic to the black solid line a = α as η decreases, while the RMSEs of the near-zero-RMSE state in CIM L0-RBCS decrease to zero (the red lines in Fig. 4b).Nakanishi et al. showed using statistical mechanics that L0-minimization-based CS (minimize x 0 s.t.y = Ax) has a stable solution corresponding to reconstructing x with zero error as long as a is below a critical value a c = α [44].Thus, the near-zero-RMSE solutions of CIM L0-RBCS are asymptotic to the perfect reconstruction solution of L0-minimization-based CS as η decreases.
To compare the properties of CIM L0-RBCS with those of LASSO, Fig. 5b shows the phase diagrams of LASSO; the blue lines are the first-order phase transition lines from the near-zero-RMSE state for various η.As η decreases, the phase transition lines from the near-zero-RMSE state in LASSO for the half-Gaussian (+) and Gaussian (±) become asymptotic to the two black dotted lines, while the RMSEs of the nearzero-RMSE state in LASSO decrease to zero (the blue lines in Fig. 4).Donoho et al. showed that L1-minimization-based CS (minimize x 1 s.t.y = Ax) can reconstruct a sparse source signal x with zero error as long as a is below a critical value, shown by the black dotted lines in Fig. 5 [68].Thus, the near-zero-RMSE solutions of LASSO become asymptotic to the perfect reconstruction solution of L1-minimization-based CS as η decreases.
CIM L0-RBCS and LASSO have these asymptotic properties even when the source signals are from the Gamma (+) and bilateral Gamma (±) (see Supplementary Fig. 3b).Note that we have theoretically proved that the asymptotic property of CIM L0-RBCS is invariant to differences in the probability distributions of the source signal by applying a perturbation expansion to the MEs (24)(25)(26) in the limit η → +0 (see section 2.7).Thus, we have confirmed this theoretical result numerically.
On the other hand, when A 2 s = 250, the first-order phase transition lines of CIM L0-RBCS are not asymptotic to the black solid line a = α, as shown in Supplementary Figs. 2 and 3a.Around η = 0.1, the phase transition line is closest to a = α.
The black dotted dashed lines in Fig. 5a shows the lower bounds of the first-order phase transition lines of CIM L0-RBCS in the limit A 2 s → ∞.The lower bound lines are above the critical line (black dotted line) of L1-minimization-based CS when the compression rate α is lower than around 0.5 for the half-Gaussian (+) and 0.7 for the Gaussian (±).The lower boundary property in Fig. 5a is satisfied even in the case of source signals from the Gamma (+) and bilateral Gamma (±) (see Supplementary Fig. 3b).On the other hand, there are no such lower bounds when A 2 s = 250 (Supplementary Figs. 2 and 3a).

Basin of attraction when β = 0
To check the practicality of CIM L0-RBCS, we verified the basin of attraction of Algorithm 1.To make the basin wider, we heuristically introduced a linear threshold attenuation wherein the threshold η was linearly lowered from η init to η end as the minimization process was alternated (see Algorithm 1).First, we carried out numerical experiments to verify the size of the basin of attraction for various initial values η init for fixed η end = 0.01 in the case of no observation noise (i.e.β = 0).As shown in Fig. 6a, the basin of attraction tended to be widened by selecting a higher initial threshold η init than η end .As the compression rate α decreased, this tendency became more marked, especially in the Gaussian case (±).Next, we sought to confirm how well Algorithm 1 converged on the near-zero RMSE state given by the MEs (24)(25)(26) when starting from an initial state r = 0 for various η init (Fig. 6b).As demonstrated in Fig. 6b, when the sparseness a was lower than the lower bound of the first-order phase transition points (the black dotted dashed line in Fig. 5a), Algorithm 1 with η init = 0.6 converged to the solutions (red lines) of the MEs (24)(25) (26), whereas it failed to converge to the solutions for other values of η init .Compared with the RMSE profiles of LASSO in Fig. 6b, Algorithm 1 exceeded LASSO's estimation accuracy under almost all of the conditions in which LASSO has a small error.
The properties for the source signals taken from the Gamma (+) and bilateral Gamma (±) (see Supplementary Fig. 4) are similar to those in Fig. 6.

Performance of CIM L0-RBCS and LASSO when β = 0
Moreover, to check the practicality of CIM L0-RBCS, we verified its accuracy and convergence in the presence of observation noise (i.e.β = 0).We searched for the optimal threshold values that would give the minimum RMSEs of CIM L0-RBCS and those of LASSO (Figs. 7a and 8a) and computed the difference between their minimum RMSEs (Figs. 7b and 8b) under the optimal threshold for each method when β = 0.01, 0.05, and 0.1.The minimum RMSE was obtained by conducting a grid search on the set of solutions to the MEs (24)(25) (26) and the MEs (44)(45) (46) in the range 0.002 ≤ η ≤ 0.5 at each point (a, α).These figures show cases of the half-Gaussian (+) and Gaussian (±) source signals.As indicated in Figs.7a and 8a, as β decreases, the phase transition lines from the-near-zero RMSE state in CIM L0-RBCS under the optimal threshold approaches the critical line (black solid line) of L0-minimizationbased CS, and the RMSEs of CIM L0-RBCS under the optimal threshold decreases.As shown in Figs.7b and 8b, the RMSEs of LASSO are higher than those of CIM L0-RBCS under almost all of the conditions in which LASSO has an error less than 0.2; thus, CIM L0-RBCS exceeds LASSO's estimation accuracy under the optimal threshold for each method.
Next, for the case of observation noise, we determined whether the output of Algorithm 1 with A 2 s = 10 7 converged on solutions to the MEs (24)(25)(26) when starting from the initial state r = 0 and η init = 0.6.As shown in Figs.7c and 8c, near or at the phase transition points, Algorithm 1 converged to the solutions of the MEs (24)(25) (26).
The properties for the source signals from the Gamma (+) and bilateral Gamma (±) (see Supplementary Figs. 5 and 6) are similar to those in Fig. 7 and 8.

Performance of CIM L0-RBCS on realistic data
Finally, we evaluated the performance of CIM L0-RBCS and other methods on realistic data.We used MRI data obtained from the fastMRI datasets [69].A Haar-wavelet transform (HWT) was applied to the data, and 86.6% of the HWT coefficients were set to zero to create a signal spanned by Haar basis functions with a sparseness of 0.134 (left panel of Fig. 9a).The k-space data shown in the middle panel of Fig. 9a was obtained by calculating the discrete Fourier transform (DFT) from the signal of the left panel of Fig. 9a, and 40% of k-space data were undersampled at random red points in the middle panel of Fig. 9a to create an observation signal with a compression rate of 0.4.The right panel of Fig. 9a shows an image with incoherent artifacts obtained by zero-filling Fourier reconstruction from the randomly undersampled k-space data.
To achieve higher reconstruction accuracy from the undersampled signal, we formulated the following implementable optimization problem on CIM with L0 and L2 norms [70]: where x is a source signal, y is a k-space undersampling signal, F is a DFT matrix, S is an undersampling matrix, Ψ is a HWT matrix, ∆ v and ∆ h are the secondderivative matrices for the vertical and horizontal directions respectively, and γ and λ are regularization parameters.Under the variable transformation r = Ψx, the local field vector and the mutual interaction matrix for CIM L0-RBCS can be set as 9b shows images (and RMSEs) reconstructed from CIM L0-RBCS (left panel of Fig. 9b), LASSO [41] (middle panel of Fig. 9b), and L1-minimization-based CS implemented in CVX [71,72] (right panel of Fig. 9b).As indicated in the images surrounded by the red circles in these panels, CIM L0-RBCS gave the most accurate reconstruction.
We evaluated the RMSEs of the three methods as a function of the threshold η.As shown in Fig. 9c, the blue line with error bars is the RMSE of CIM L0-RBCS obtained from ten trials, the red line is the RMSE of LASSO, and the circle is the RMSE of L1 minimization-based CS.There is an optimal value of η to minimize the RMSEs of both CIM L0-RBCS and LASSO because of the trade-off between detecting small nonzero elements and eliminating incoherent artifacts by thresholding.The RMSE of CIM L0-RBCS was lower than those of the other methods in a wide range of η.

Summary and Conclusion
We proposed a quantum-classical hybrid system that performs CIM and CDP alternately to optimize r and σ.To evaluate the performance of CIM L0-RBCS, we derived the W-SDE from the master equation for the density operator of a system consisting of N OPOs and a measurement-feedback circuit.After that, we obtained the MEs for CIM L0-RBCS from the W-SDE (15) and the simultaneous equations (11).
As shown in Figs. 4, 7c and 8c and Supplementary Figs. 1, 5c and 6c, we confirmed that the theoretical results obtained from the MEs were consistent with the numerical results of Algorithm 1 regardless of whether observation noise exists in the observed signal y.In particular, the theoretical results in the limit A 2 s → ∞ were in good agreement with the results of Algorithm 1 with A 2 s = 10 7 .Because A 2 s = 10 7 is on the same order as A 2 s in the experimental CIMs [9,10], we expect that the MEs (24)(25)(26) can be used to evaluate real experimental CIMs.
In the case of no observation noise, we theoretically showed that the performance of CIM L0-RBCS in principle approaches the theoretical limit of L0-minimization-based CS [44] at high pump rates (see Fig. 5a and Supplementary Fig. 3b).From a mathematical perspective, the theoretical limit line a = α is the condition when the rank of a matrix composed of column vectors of an observation matrix corresponding to non-zero elements of the source signal is full.Thus, it is impossible for any system to go beyond this line mathematically.As described above, because the theoretical results in the limit A 2 s → ∞ are in good agreement with those of Algorithm 1 with A 2 s = 10 7 , we expect that the theoretical performance limit of real experimental CIMs is close to this ideal limit.
In the case of observation noise, we theoretically showed that the RMSEs of CIM L0-RBCS are lower than those of LASSO for almost all conditions in which LASSO has an error less than 0.2, and thus CIM L0-RBCS exceeds LASSO's estimation accuracy under the optimal threshold for each method (see Figs. 7a, 7b, 8a and 8b and Supplementary Figs.5a, 5b, 6a and 6b).
However, there is a problem regarding the basin of attraction.As numerically demonstrated in Fig. 6 and Supplementary Figs. 4, when there is no observation noise, Algorithm 1 cannot reach the theoretical performance limit if it starts from the practical initial condition r = 0.However, even in such a situation, Algorithm 1 exceeds LASSO's estimation accuracy until the lower bound of the first phase transition points of CIM L0-RBCS (Fig. 6 and Supplementary Fig. 4).On the other hand, when there is observation noise, under the practical initial condition r = 0, Algorithm 1 gets very close to or achieves the theoretical performance limit given by the ME (see Figs. 7c and  8c and Supplementary Figs.5c and 6c).
Finally, we confirmed using realistic data that CIM L0-RBCS gave the most accurate reconstruction compared with LASSO and L1-minimization-based CS (Fig. 9).
Therefore, we can conclude that the performance of CIM L0-RBCS in principle approaches the theoretical limit of L0-minimization-based CS at high pump rates, exceeds that of LASSO, and moreover in practical situations exceeds LASSO's estimation accuracy.
A detailed interpretation and discussion of these results is given below.

Effectiveness of CIM in support estimation
In Algorithm 1, the c-amplitude is always initialized to c = 0 in the initial stage of the support estimation even when r is initialized to the true signal value, i.e. x • ξ.In this situation, the solutions of Algorithm 1 match the near-zero-RMSE states of the MEs very well, as demonstrated in Fig. 4 and Supplementary Fig. 1; hence, the simulated CIM can reconstruct the support vector up to the theoretical limit.Therefore, the simulated CIM can reach the ground state which minimizes H with respect to σ when r is fixed.

Correctness of assumptions
To derive the MEs (21)(22) (23), we derived an approximate value for X( h, t) SDE of each OPO pulse by replacing the state variables in the second-order coefficient of the power of the quantum noise with average values of the state variables (see Eq. ( 29)).As shown in Figs.4b, 7c, and 8c, the ME derived under this approximation has good accuracy at the values of A 2 s used in the actual equipment of the CIM.However, as shown in Fig. 4a, some solutions of the ME did not match the numerical solutions of Algorithm 1 for smaller values of A 2 s .Thus, this approximation is possible if the mutual injection field is much larger than the noise in the steady state where the c-amplitude has grown.

Basin of attraction and its dependency on threshold
To make the basin of attraction of Algorithm 1 wider, we heuristically introduced a linear threshold attenuation in which the threshold η linearly decreases as the alternating minimization proceeds.We confirmed that the basin of attraction widens as a result of lowering η from a higher initial threshold η init to a lower terminal threshold η end (see Fig. 6 and Supplementary Fig. 4).
According to the definition of the injection field for each OPO pulse in Eq. ( 9), the threshold η acts as an external field to give a negative bias for the OPO pulses to take the down state.By initially giving a large negative external field, almost all of the OPO pulses take the π-phase state, and thus, almost all of the {H(X j )} j=1,•••,N take zero in the initial stage of the alternating minimization process.In the initial stage, the system can easily reach the ground state under a strong negative bias because the phase space, which consists of a small number of up-state OPO pulses, is simple.Then, through the alternating minimization process, the system tracks gradual changes in the ground state due to incremental increases in the number of up-state OPO pulses by gradually sweeping out a negative external field.Finally, the system achieves the ground state at the terminal threshold η end .This is the qualitative interpretation of how the basin of attraction of Algorithm 1 widens by linearly lowering the threshold.
However, as demonstrated in Fig. 6b, when there is no observation noise, the system fails to converge to the near-zero-RMSE solutions beyond the lower bound line of the first-order phase transition points.We suspect that there might be many quasi steady states beyond the lower bound line, as in the spin-glass phase [73]; thus, the system might become trapped in one of the quasi steady states.
On the other hand, when there is observation noise, as demonstrated in Figs.7c and 8c, the system converges to t near-zero-RMSE solutions even nearby the phase transition line when it starts from the practical initial condition r = 0.It was suggested that the symmetries of the system allow for the creation of quasi steady states [74].We conjecture that observation noise could break the symmetries for quasi steady states.

Plan to improve CIM L0-RBCS
To achieve the theoretical performance limit when there is no observation noise, we need to construct a full quantum system in which both the support estimation and the signal estimation are implemented on the CIM.We expect that due to the minimum gain principle [7,50], the full quantum system could overcome the quasi-steady-state problem discussed above.

3:
Minimize H with respect to σ by CIM: σ = CIM support estimation(r, η) # Initialize the c-amplitude as c = 0, and numerically integrate the W-SDE while increasing the normalized pump rate from 0 to 1.5 for five times the photon's lifetime when A 2 s = 10 7 or for two hundred times the photon's lifetime when A 2 s = 250.

4:
Minimize H with respect to r by CDP: Decrement η: η = max(η init (1?t/50), η end ) 6: end for 7: return σ and r for signal estimation.This system is realizes the alternating minimization described in Algorithm 1. Pump pulses are injected into an optical parametric oscillator (OPO) formed in a fiber ring cavity through second harmonic generation (SHG) crystal.A periodically poled lithium niobate (PPLN) waveguide device induces a phase-sensitive degenerate optical parametric amplification of the signal pulses, and each of the OPO pulses take either the 0-phase state (corresponding to the up-spin) or the π-phase state (corresponding to the down-spin) above the oscillation threshold.Part of each pulse is taken from the main cavity by the output coupler, and it is measured by optical homodyne detectors.A field programmable gate array (FPGA) calculates the feedback signal, which is then provided to the intensity modulator (IM) and phase modulator (PM) to produce the injection field described in Eq. ( 9) to each of the OPO pulses through the input coupler.H(X i ) is a binarized value, either 0 or 1, of the in-phase amplitude of the i-th OPO pulse, which is the support estimate to be transferred to the CDP.The CDP solves the linear simultaneous equation (Eq.( 11)), and the solution r i is transferred to the CIM.Gamma(+) Bilateral Gamma(±) Half-Gaussian(+) Gaussian(±) Figure 3. Four probability density functions used for generating the source signal in the numerical experiments.The half-Gaussian (+) and Gamma (+) are defined over a non-negative random variable.The Gaussian (±) and bilateral Gamma (±) are defined over a signed random variable.In a, the red lines show the first-order phase transition point a c from the nearzero-RMSE state as a function of α.The black dotted-dashed in each plot the lower bound of the first-order transition points the state in CIM L0-RBCS.In b, the blue lines the first-order transition point a c of LASSO as a function of α.The black solid line in each is a critical line indicating a boundary whether or not L0-minimization-based CS can perfectly reconstruct the source signal with no error for the non-negative case and signed case [44], while the black dotted line is a critical line indicating whether or not L1-minimization-based CS can perfectly reconstruct the source signal with no error for the non-negative case and signed case [68].

Figure 1 .Figure 2 .
Figure1.Quantum-classical hybrid system for L0-RBCS.To estimate the Ndimensional support vector σ and N -dimensional signal vector r, this system solves the two-fold optimization problem by alternately performing two minimization processes; a the quantum machine optimizes σ to minimize H with the given r, and b the classical digital processor optimizes r to minimize H with given σ.

Figure 4 .
Figure 4. Comparison of solutions of MEs with solutions of Algorithm 1: cases of no observation noise (i.e.β = 0) and half-Gaussian (+) and Gaussian (±) source signals.The RMSEs of the solutions are plotted as a function of sparseness a for various thresholds η and compression rates α. a Comparison of solutions of the MEs (21)(22)(23) and those of Algorithm 1 with A 2 s = 250.b Comparison of solutions of the MEs (24)(25)(26) for A s → ∞ and those of Algorithm 1 with A 2 s = 10 7 .In a and b, the red and green lines respectively indicate RMSEs of the near-zero RMSE state and non-zero RMSE state in CIM L0-RBCS, which were obtained with the MEs (21)(22)(23) with A 2 s = 250 and the MEs (24)(25)(26).The blue lines are the RMSEs of the near-zero RMSE state in LASSO, which were obtained with the MEs (44)(45)(46) with the same threshold value of η indicated above the graphs.The circles and error bars represent the mean values and standard deviations of ten trial solutions numerically obtained by Algorithm 1.To confirm the existence of solutions of Algorithm 1 corresponding to the near-zero RMSE states indicated by the MEs, r was initialized to the true signal value, i.e. x • ξ, and η was kept constant by setting η init = η end to the value of η indicated above the graphs.For all graphs, K = 0.25 and N = 2000. 0

Figure 5 .
Figure 5. Phase diagrams of CIM L0-RBCS in the limit A 2 s → ∞ and LASSO for various η: cases of no observation noise (i.e.β = 0) and half-Gaussian (+) and Gaussian (±) source signals.a Phase diagrams of CIM L0-RBCS.b Phase diagrams of LASSO.In a, the red lines show the first-order phase transition point a c from the nearzero-RMSE state as a function of α.The black dotted-dashed in each plot the lower bound of the first-order transition points the state in CIM L0-RBCS.In b, the blue lines the first-order transition point a c of LASSO as a function of α.The black solid line in each is a critical line indicating a boundary whether or not L0-minimization-based CS can perfectly reconstruct the source signal with no error for the non-negative case and signed case[44], while the black dotted line is a critical line indicating whether or not L1-minimization-based CS can perfectly reconstruct the source signal with no error for the non-negative case and signed case[68].

ηHalfηFigure 6 .Figure 9 .
Figure 6.Basin of attraction of CIM L0-RBCS depending on the initial threshold η init : cases of no observation noise (i.e.β = 0) and half-Gaussian (+) and Gaussian (±) source signals.a Final states of Algorithm 1 when starting from various initial states for various η init .The pairs of points connected by a line indicate the initial and final states for each trial when η end = 0.01 and A 2 s = 10 7 .b Final states of Algorithm 1 when starting from an initial state r = 0 for various η init .The circles and error bars represent the mean values and standard deviations of twenty trial solutions numerically obtained by Algorithm 1 with η end = 0.01 and A 2 s = 10 7 .The red lines show the solutions of the MEs (24)(25)(26) with near-zero RMSE when η = 0.01 and A 2 s → ∞, while the blue lines indicate RMSEs of LASSO when η = 0.01.The black lines are the lower bounds of first-order phase transition points of the CIM L0-RBCS.K = 0.25 and N = 4000.