Mean field analysis of reverse annealing for code-division multiple-access multiuser detection

We evaluate the typical ARA performance of the CDMA multiuser detection by means of statistical mechanics using the replica method. At first, we consider the oracle cases where the initial candidate solution is randomly generated with a fixed fraction of the original signal in the initial state. In the oracle cases, the first-order phase transition can be avoided or mitigated by ARA if we prepare for the proper initial candidate solution. We validate our theoretical analysis with quantum Monte Carlo simulations. The theoretical results to avoid the first-order phase transition are consistent with the numerical results. Next, we consider the practical cases where we prepare for the initial candidate solution obtained by commonly used algorithms. We show that the practical algorithms can exceed the threshold to avoid the first-order phase transition. Finally, we test the performance of ARA with the initial candidate solution obtained by the practical algorithm. In this case, the ARA can not avoid the first-order phase transition even if the initial candidate solution exceeds the threshold to avoid the first-order phase transition.


I. INTRODUCTION
The code-division multiple-access (CDMA) multiuser detection has been used in various communication systems [1].The theoretical performance of CDMA multiuser detection has been analyzed by means of statistical mechanics [2][3][4][5].The CDMA multiuser detection is regarded as a type of signal recovery problems, similar to compressed sensing [6][7][8][9].Statistical-mechanical analyses for signal recovery problems focus on the inference of the original information from the degraded information with noise.The noise can be physically regarded as thermal fluctuations.By tuning the strength of the thermal fluctuations, the original signal can be estimated from the degraded one.
In addition to thermal fluctuations, quantum fluctuations may be used to estimate the signal.Several studies have demonstrated that quantum fluctuations such as the transverse field do not necessarily improve the performance of the inferences for image restoration, Sourlas codes, and CDMA [10][11][12].The optimal decoding performance with quantum fluctuations is inferior to that with thermal fluctuations in Bayes optimal cases.However, in certain non-Bayes optimal cases; for example, where a lower temperature than the true noise level is set, the decoding performance with finite quantum fluctuations and thermal fluctuations is superior to that with only * shunta.arai.s6@dc.tohoku.ac.jp thermal fluctuations.That implies the potential of the combination of quantum and thermal fluctuations for inference problems.
The performance of optimization algorithms with quantum fluctuations, which is known as quantum annealing (QA) [13][14][15][16][17][18] or adiabatic quantum computation (AQC) [19,20], is equal to or better than that of an optimization algorithm with thermal fluctuations [21,22], which is known as simulated annealing (SA) [23].The physical implementation of QA is realized by the quantum annealer [24][25][26][27][28].The quantum annealer has been implemented in numerous applications, such as portfolio optimizations [29,30], traffic optimization [31], item listing for E-commerce [32], automated guided vehicles in factories [33], machine learning [34][35][36][37][38], quantum simulation [39][40][41], material design [42] and decoding algorithm [43].In a closed system, QA begins from the ground state of the transverse field term and the transverse field strength is gradually reduced.Following the Schrdinger equation, the trivial ground state evolves adiabatically into a nontrivial ground state of the target Hamiltonian, which corresponds to the solution of combinatiorial optimization problems.The quantum adiabatic theorem guarantees a theoretically sufficient condition to obtain the ground state in QA [44].The theorem indicates that the total computational time for obtaining the ground state is characterized by the minimum energy gap between the ground state and first exited state.The energy gap is related to the order of the phase transition.In the case of the first-order phase transition, the computational time for search-ing the ground state increases exponentially [45][46][47][48], which is the worst case of QA.
To avoid the first-order phase transition, many methods are proposed for example, QA with a non-stoquastic Hamiltonian [49][50][51][52] , inhomogeneous driving of the transverse field [53,54] , and reverse annealing (RA) [55,56].The RA is a protocol to restart the quantum dynamics starting form the resulting state of the standard procedure of QA.We expect that the RA leads to a closer solution to the ground state as its output.To assess the performance of RA, we carefully classified its implementation into two methods: adiabatic reverse annealing (ARA) and iterated reverse annealing (IRA).The main difference between the ARA and IRA is to incorporate the resulting state.One is to implement the resulting state by modification of the initial Hamiltonian and the other introduce it as the initial condition.
In the ARA, we modify the initial Hamiltonian according the resulting state.We assume that the resulting state is a candidate solution, which is sufficiently close to the ground state of the original problem we wish to solve.We prepare the initial Hamiltonian in the ARA such that its ground state is the candidate solution.The procedure of the ARA is outlined as follows: We start from the ground state of the initial Hamiltonian.Next, we gradually increase the effects of quantum fluctuations and search locally around the candidate solution.Thereafter, we gradually decrease the effects of quantum fluctuations.When the effects of quantum fluctuations disappear, the ground state or lower energy state of the original problem can be obtained.The ARA is rather theoretical approach for understanding the property of the RA.Theoretical analysis has indeed shown that the ARA can avoid the first-order phase transition for the p-spin model [57].However, this protocol has not been implemented in the current quantum annealer.
The procedure of the IRA is slightly different from the ARA.The difference is that the IRA starts from a classical state, which corresponds to the candidate solution without introducing the additional Hamiltonian.A similar protocol to the IRA is feasible in the current quantum annealer.The performance of the IRA can be analyzed from the dynamics and it significantly depends on effect of heat bath.In a closed system, the IRA has not enhanced the performance of the QA [58].In an open system, the IRA has improved the performance of the QA by incorporating the relaxation mechanisms [59].
In this study, we focus on the ARA because it dramatically enhances the performance of the QA, and its performance can be analyzed by the statistical mechanics.To the best of our knowledge, it remains unknown whether or not the ARA is useful for certain practical problems.We apply the ARA to the CDMA multiuser detection, which is a representative example in signal recovery problems.The CDMA model is mainly characterized by the ratio of the number of users to that of the measurements, which is called the pattern ratio.In the low-temperature regions and the intermediate pattern ratio, the CDMA model has two solutions.This phenomenon reveals the existence of the first-order phase transition, which degrades the estimation efficacy.We use the ARA to mitigate or avoid the estimation difficulty.In the ARA process, we set the initial Hamiltonian.The initial Hamiltonian is interpreted as prior information of the original signal in the context of the inference problems.We expect that the prior information of the original signal will mitigate the estimation difficulty We consider the marginal posterior mode (MPM) estimation by the ARA [60].The estimated signal corresponds to the expectation of the signal over the Gibbs-Boltzmann distribution.The MPM estimation can be performed in the current quantum annealer, which provides samples from the density matrix incorporating both thermal and quantum fluctuations in about 20µs [61,62].We analyze the average MPM estimation performance with ARA at a finite temperature using the replica method.The MPM estimation with the ARA is regarded as the MPM estimation with quantum fluctuation incorporating the prior information of the original signal.The typical performance of the MPM estimation with quantum fluctuations for the CDMA multiuser detection has been analyzed in the previous study [12].The connection between the present and previous studies is presented in Sec.II A.
In the ARA, we need to prepare for the initial candidate solution.In the previous study [57,58], they have not cared about how to prepare for the initial candidate solution.We investigate whether or not we can prepare for the proper initial candidate solution to avoid the first-order phase transition with commonly used algorithms.We test the performance of the ARA with the initial candidate solution obtained by these practical algorithms.Although the implementation of the ARA in the current quantum annealer has not yet been realized, our results provide the first theoretical demonstration of the ARA as a practical technique for signal recovery problems.
The remainder of this paper is organized as follows.In Section II, we review the previous study and present the formulation of the CDMA model with quantum fluctuations.In Section III, we extend the formulation for the ARA.We derive the free energy under the replica symmetry (RS) ansatz and the static approximation.In Section IV, we illustrate the phase diagrams in the ARA.At first, we consider oracle cases where the initial candidate solution is randomly generated from the probability distribution given the fraction of the original signal in the initial state.To verify the RS solutions, we perform quantum Monte Carlo simulations.Next, we check whether or not we can prepare for the proper initial candidate solution to avoid the first-order phase transition with commonly used algorithms.Finally, we test the performance of the ARA with the initial candidate solution attained from these practical algorithms.In Section V, we conclude the study and discuss the future research directions.

II. CDMA MODEL WITH QUANTUM FLUCTUATIONS
At first, we review the previous study [12] and show its relationship with the MPM estimation with the ARA in Sec.II A. Next, we formulate the classical CDMA model and move onto the quantum system in Sec.II B.
The previous study [12] analysed the performance of the MPM estimation for the CDMA model with quantum fluctuations under the standard protocol of QA.In particular, they shed light on the difference between quantum and thermal fluctuations.In other words, they compared the performance between SA and QA.In the case by SA, one controls the strength of thermal fluctuation through a parameter of temperature.Depending on the noise in the received signal, they found the optimal strength of the thermal fluctuation known as the Nishimori temperature [63] to retrieve the original signal in the context of the CDMA.In the previous study, they investigated the existence of the optimal strength of the quantum fluctuation similarly to the case with thermal fluctuation.They showed that the MPM estimation with quantum fluctuations could partially improve its performance compared to the case without quantum fluctuation.However, the MPM estimation with quantum fluctuations did not archive the optimal MPM performance found in the case only with thermal fluctuations.In this sense, the thermal fluctuation is superior to the quantum fluctuation in the retrieval of the original signal of the CDMA.Nevertheless, one of the crucial bottlenecks of the protocol in both of SA and QA to retrieved the original signal still exists.There is the first-order phase transition in the case with the intermediate pattern ratio in the low-temperature region.Here the temperature is a control parameter of the MPM estimation.The existence of the first-order phase transition hampers efficient retrieval of the original signal and needs the long computation time of its execution.As shown in the previous study, quantum fluctuation could not avoid nor mitigate the first-order phase transition.We thus investigate the potential of ARA, which is slightly different from the standard protocol of QA, in the present study.In this sense, our study is placed in position as an extension of the MPM estimation with quantum fluctuations by using a different protocol of the standard QA.

B. Formulation
The main concept of the CDMA model is as follows: The digital signal of each user is modulated and transmitted to a base station through fully synchoronous channels.By demodulating the received signal composed of the multiuser signals and noises, we infer the original signal from the provided information.The following formulation is mainly based on the previous study of the CDMA model with quantum fluctuations [12].They add the transverse field to the original CDMA model and compute the partition function following the prescription of the statistical mechanics.They used the Suzuki-Trotter decomposition to deal with quantum fluctuation written in the transverse field and the replica method to compute the averaged free energy over the quenched randomness related to the signals and modulation.In the present study, we employ the same methods to tackle the MPM estimation of CDMA by using the ARA and setting the initial Hamiltonian depending on the initial candidate solution.
We consider that N users communicate via fully synchronous channels.At the base station, the receiver obtains the signal as follows: where ξ i ∈ {±1}, (i = 1, . . ., N) is the original information and is the spreading code for each user i.The length of the spreading codes for each user i is represented by K.The channel noise ǫ µ is added into the received signal.The received signal (1) can be expressed as for which the following notations are used: We assume that the spreading codes and original signal are independently generated from the uniform distribution: We consider the Gaussian channels and ǫ k is independently generated from the Gaussian distribution as follows: where T 0 = β −1 0 is the true noise scale.In the CDMA multiuser detection, we estimate the original signal from the received output signal and the spreading codes that are prepared for each user in advance.Because the output signal fluctuates owing to noise, we formulate this problem as Bayesian inference.Subsequently, we introduce the posterior distribution of the estimated signal σ = (σ 1 , . . ., σ N ) T as We define the likelihood as where β = 1/T is the inverse temperature in statistical mechanics and corresponds to the estimated channel noise scale.
If the true noise level is known, the estimation performance is the best and Bayes optimal.According to Eqs. ( 8) and ( 9), the posterior distribution can be written using the Gibbs-Boltzmann distribution with the Hamiltonian H (σ), as follows: where Z is the partition function and H init (σ) is the initial Hamiltonian, which represents the prior information of the estimated signal.We generally assume that the prior of the estimated signal follows the uniform distribution In this case, we can omit the initial Hamiltonian from Eqs. ( 10) and (11).
To estimate the original signal, we consider the MPM estimation.The estimation performance can be evaluated by the overlap between the original and estimated signal as M = 1/N N i=1 ξ i sgn σ i , where • is the expectation over the posterior distribution P(σ|y) and sgn(•) is the signum function.This quantity is expected to exhibit a "self-averaging" property in the thermodynamics limit N → ∞.This means that the observables, such as the overlap for a quenched realization of the data y, η, and ξ, are equivalent to the expectation of itself over the data distribution P(η)P(ξ)P(y|ξ).In this case, the overlap can be expressed as lim N→∞ M = [ξ i sgn σ i ], where the bracket [•] indicates the expectation over the data distribution.
It is straightforward to extend the above formulation into the quantum mechanical version: where σz i and σx i are the z and x components of the Pauli matrices at site i, respectively.In this case, Ĥ0 consists of the z components of the Pauli matrices and ĤTF is composed of the x components of the Pauli matrices.We parameterize the Hamiltonian (14) with the annealing parameter s for application to the ARA.
As in the classical case, we consider the MPM estimation with quantum fluctuations.The posterior distribution can be written as ρ = exp −β Ĥ /Tr exp −β Ĥ where Tr denotes the summation over all possible spin configurations in the z-basis.The performance of the MPM estimation with quantum fluctuations can be evaluated by M = 1/N N i=1 ξ i sgn Tr σz i ρ .

III. MEAN FIELD ANALYSIS
Following Ref. [57], we extend the CDMA model with quantum fluctuations [12] to the ARA formulation as where λ (0 ≤ λ ≤ 1) is the RA parameter.The initial candidate solution is denoted by τ i = ±1(i = 1, . . ., N) that is expected to be close to the original signal ξ i .We introduce the probability distribution of the initial candidate solution as follows: where we define c 1 = c and c −1 = 1 − c.The number c (0 ≤ c ≤ 1) denotes the fraction of the original signal τ i = ξ i in the initial state as The prior information can be incorporated through Eq.( 19).
The main concept of the MPM estimation with the ARA is to avoid or mitigate the first-order phase transition by controlling the RA parameter and utilizing the prior information.
The typical behaviors of the order parameters such as the overlap can be obtained via the free energy.We calculate the partition function Z = Tr exp −β Ĥ and derive the RS free energy in the limit of N, K → ∞, while maintaining the pattern ratio α ≡ K/N = O(1).The free energy density f can be evaluated as −β f = lim N→∞ (1/N)[ln Z] where [•] denotes the configuration average over the data distribution P(y|ξ)P(η)P(ξ)P(τ).When computing f , we have two difficulties.The first one is the non-commutativity of the spin operator from Eq.( 16).We can not apply the mean-field analysis directly into the partition function of Eq.( 17).The second one is to compute [ln Z].In general, it is difficult to directly evaluate [ln Z].We remove two difficulties by using two techniques.
Firstly, to exclude the non-commutativity of the spin operator, we employ the Suzuki-Trotter (ST) decomposition [64] in the partition function: where in which the symbol t is the index of the Trotter slice and M is the Trotter number.We impose the periodic boundary conditions, σ i (M + 1) = σ i (1) for all i.By using the ST decomposition, we can map the quantum system into the identical classical system.The difficulty from the non-commutativity of the spin operator is removed.Above expressions, we replace σ z i (t) with the classical spin σ i (t) ∈ {−1, +1}.In this case, the symbol Tr represents the trace over the classical spins.The x component of the Pauli matrix yields the last term in Eq. (22).
Secondly, to evaluate [ln Z], we exploit the replica method [65]: /n.The symbol n denotes the number of replicas.By using the replica method, we can take the configuration average for the replicated partition function Z n and the limit of n → 0. When manipulating the configuration average over P(y|ξ)P(η)P(ξ)P(τ), we introduce the order parameters and their conjugate parameters through the delta function and its Fourier integral representation as follows: the magnetization m a (t) = (1/N) N i=1 ξ i σ ia (t), the spin glass order parameter q ab (t, , and the correlation between each Trotter slice R a (t, The conjugate parameters are denoted by ma (t), qab (t, t ′ ) (a b) ,and Ra (t, t ′ ).These conjugate parameters appear in manipulation of several integrals over order parameters to compute the partition function as detailed in Supplemental Material [66] The symbols a and b represent the replica indices.Under the RS ansatz and static approximation: we can finally obtain the RS free energy density: where in which Dz means that the Gaussian measure Dz := dz/ √ 2πe −z 2 /2 and Dy is the same as Dz.Here extr represents the extremization by changing the order parameters m, q and R, and their conjugate parameters as m, q and R. The extremum point is determined by the saddle-point conditions and characterizes the free energy density.The expression in Eq.( 23) for λ = 1 can be reduced to the RS free energy density derived in [12].The detailed derivation of Eq. ( 23) is written in Supplemental Material [66].The saddle-point equations are referred in Appendix A. Below we investigate the phase transition of the order parameters while tuning the external parameters as the strength of the transverse field, the pattern ratio etc.. Then we numerically solve the saddle-point equations for each set of external parameters.

IV. NUMERICAL RESULTS
In this section, we evaluate the typical performance of the ARA based on the results attained in Section III.In Section IV A, we consider the oracle cases where the initial candidate solution is randomly generated from Eq.( 19) given the fraction of the original signal in the initial state.In Section IV B, we consider the practical cases where we prepare for the initial candidate solution with commonly used algorithms.We compare the performance of the ARA with the oracle cases and the practical cases.

A. Analysis of the ARA in oracle cases
We numerically solve the saddle-point equations in Eqs.(A1) to (A6) with the temperature T = 0.1.We set the several RA parameter λ.We start from the ARA with λ = 1 which corresponds to the vanilla QA [12].We show that the CDMA model has the first-order phase transition in the intermediate pattern ratio.Next, we consider the classical case with λ = 0 to validate the RS ansatz without the static approximation.Finally, we move onto the finite λ.We exhibit that the ARA can avoid or mitigate the first-order phase transition if we prepare for the proper initial candidate solution.Let us begin with the ARA with λ = 1.The phase diagrams for the true noise scale T 0 = 0, 0.05 and 0.1 are displayed in Fig. 1.The blue solid curve and orange dash-dotted curve indicate the spinodal curves where the solutions for each initial condition disappears in Figs.1(a) -1(c).Two solutions coexist between the two spinodal curves.From these figures, we can establish the existence of the first-order phase transition in the intermediate pattern ratio and under the weak strength of the transverse field.The green dotted curve denotes the critical point at which the RS free energy takes the same value.In Fig. 1(c), we do not write down the curve because we can not distinguish the critical point from the spinodal points in this scale.Higher noise results in a narrower region in which the two solutions coexist.Although the noise mitigates the first-order phase transition, it decreases the overlap between the original signal and the estimated one.
Next, we consider why the first-order phase transition is troublesome in estimating the original signal.The difficulty of estimating the original signal is related to the free energy landscape.We take Fig. 1(a) as an example.On the right side of spinodal curve 2, it is easy to estimate the original signal because the free energy exhibits a minimum, which is a good estimator.When we set the pattern ratio as α = 0.6, we encounter the first-order phase transition at s ≃ 0.8.The free energy landscape has two valleys.At spinodal curve 2, the free energy landscape is transformed into a simple valley.In this case, it is comparatively easy to estimate the original signal.For α = 0.5, the spinodal curve 2 does not exist.The free energy landscape maintains two valleys.We can not efficiently estimate the original signal because the metastable state remains.For α = 0.4, the critical point does not exist.In this case, we can not obtain the original signal informationtheoretically.The ground state or low energy state does not correspond to the original signal at s = 1.The minima of the free energy do not provide us with an effective estimation.
To verify the RS ansatz and the static approximation, we perform quantum Monte Carlo simulations for the CDMA model.We set the system size as N = 500, the Trotter number as M = 50, the temperature as T = 0.1, and the true noise scale as T 0 = 0. We use a 100000 Monte Carlo step (MCS) average after 50000 MCS equilibrations for each instance.We take the configuration average over the spreading codes and the original signals by randomly generating 50 instances.We plot the behavior of the order parameters with respect to the pattern ratio for the fixed annealing parameter s = 0.9 in Fig. 2 and the annealing parameter for the fixed pattern ratio α = 0.6 in Fig. 3.The error bar is given by the standard deviation.The results obtained by the quantum Monte Carlo simulations are the averages over all of the Trotter slices.Following Ref. [67], we adopt the magnetization to quantify the performance of the MPM estimation.In this study, we refer to the solution representing the "spinodal 1" curve as "branch 1" and to the solution representing the "spinodal 2" curve as "branch 2".According to Fig. 2, the results obtained by the quantum Monte Carlo simulations are consistent with the RS solutions, with the exception of the low pattern ratio.Figure.3shows that the numerical results for the magnetization is consistent with the RS solutions, except for the intermediate values of the annealing parameter.The numerical result of the correlation between the Trotter slices does not follow the RS solutions other than the large annealing parameter, which is close to 1 [68].To investigate the deviations between the numerical results and the RS solutions due to the replica symmetry breaking (RSB), we compute the Almeida-Thouless (AT) condition and the entropy.The details of these formula are written in Appendix A. In these problem settings, the AT condition is not broken, and the entropy is positive.The deviations between the the numerical results and the RS solutions probably result from the breaking of the static approximation.

ARA with λ = 0
To support the RS ansatz without the static approximation, we consider the ARA with λ = 0.In this case, the quantum part in Eq. ( 17) disappears.The experimental settings are the same as those in Fig. 3.We set α = 0.6 and T 0 = 0 in Fig. 4(a), and α = 0.62 and T 0 = 0.05 in Fig. 4(b).We consider three initial conditions: c = 0.7, 0.8, and 0.95.The initial candidate solutions are generated from Eq. ( 19) given a fixed fraction c.The error bar is given by the standard deviation.Each curve represents the RS solutions, and each symbol denotes the numerical results obtained by the Markov-chain Monte Carlo simulations.It can be observed that the numerical results are consistent with the RS solutions.We can see that the deviations between the numerical results and the RS solutions are not the breaking of the RS ansatz to the breaking of the static approximation.For λ = 0 with or without noise, the ARA can avoid the first-order phase transition if we prepare for the proper initial conditions.In the next Section, we analyze the general cases in detail.

ARA with finite λ
We consider the ARA with finite λ.The experimental settings are the same as those in Fig. 1(a).Figure 5 presents the phase diagram of the CDMA model in the ARA for α = 0.6 and 0.5.We consider four initial conditions: c = 0.7, 0.8, 0.9, and 0.95.Each curve represents a point of the first-order phase transition.We can observe from Figs. 5(a) and 5(b) that the first-order phase transition can be avoided if the initial state is close to the original signal.As the information regarding the original signal is increased, the region for avoiding the firstorder phase transition is broadened.In Fig. 5(b), the region in which the first-order phase transition can be avoided is narrower than that in Fig. 5(a).For a lower pattern ratio, further information regarding the original signal is initially required to avoid the first-order phase transition.We also investigate the stability of the RS solutions and find that the RSB does not happen for finite λ To analyze the extent to which the difficulty in obtaining the original signal is mitigated by the ARA, we plot the differences in the magnetization ∆m between the two local minima at the first-order phase transition in the case of α = 0.6 and 0.5 in Fig. 6.As discussed in Ref. [57], the rate of quantum tunneling between two local minima in the free energy landscape is related to ∆m. Figure 6 indicates that ∆m decreases as c increases.For finite λ, ∆m is smaller than that of the vanila QA (λ = 1).Even though the ARA can not eliminate the first-order phase transition, the two local minima of the free energy become closer than those of the original one.The result demonstrates that the ARA enhances the effects of the quantum tunneling for the CDMA model.In the ARA, we add the bias towards the original signal through the initial Hamiltonian.Since the bias removes or softens the free energy barrier, the ARA can avoid or mitigate the first-order phase transition.
We consider the noise effects for the CDMA model in the ARA.The experimental settings are the same as those illustrated in Fig. 1(b).Figure 7 displays the phase diagrams of the CDMA in the ARA for α = 0.62 and 0.57.The qualitative behaviors of the systems are approximately the same as those in the noiseless cases.The regions in which the firstorder phase transition can be avoided are larger than those of the noiseless cases because the first-order phase transition is weakened owing to the noise effects.Figure 8 presents ∆m in the case of α = 0.62 and 0.57.We can see that ∆m is smaller than in the noiseless cases.In the small noisy cases,the free energy barrier is lower than in the noiseless cases.The ARA works better in the small noisy cases than the noiseless cases.
To validate the RS solutions under the static approximation for finite λ, we perform quantum Monte Carlo simulations.The experimental settings are the same as those in Figs. 3  and 4. We set the RA parameter as λ = 0.8, and the initial conditions as c = 0.7 and 0.9.In Fig. 8, the order parameters for c = 0.7 still exhibit a jump.In the case of c = 0.9, it can be observed that the first-order phase transition can be avoided.We can see the deviations between the RS solutions and the numerical results are the same as in Fig. 3.In these problem settings, the RSB does not happen from the results of the saddle-point equations.Although the numerical results do not entirely match the RS solutions, the qualitative behaviors of the numerical results to avoid the first-order phase transition are similar to those of the RS solutions.

B. Analysis of the ARA in practical cases
In Section.IV A, we assume that the initial candidate solution is randomly generated from Eq.( 19) with a fixed c. Practically, we need to prepare for the initial candidate solution by some algorithms.At first, we examine whether or not we can prepare for the proper initial condition to avoid the first-order phase transition with commonly used algorithms.Next, we evaluate the performance of the ARA with the initial candidate solutions obtained by the algorithms.

How to prepare for the initial candidate solution
To prepare for the initial candidate solution, we adopt SA, simulated quantum annealing (SQA), and the approximate message passing (AMP) algorithm [69][70][71].To perform SA and SQA, we take advantage of OpenJij, an open-source library for heuristic optimization problems in Python [72].The implementation of the AMP algorithm is based on the paper [70].We perform three algorithms for different 50 instances.For SA and SQA, we carry out 51 different initial conditions for each instance.We set the system size as N = 8, 16, 32, 64, 128, 256 and 512.We check the dependence of c on N with these algorithms in Fig. 10.We compute c from the relationship c = (1+M)/2.We define the threshold required to avoid the first-order phase transition as c min .We calculate c min from the saddle-point equations.We consider the region where the spinodal curve 2 does not exist, and the critical curve exists, for example 0.45 ≤ α ≤ 0.57 in Fig. 1.In this region, to avoid or mitigate the first-order phase transition is crucial to estimate the original signal efficiently.
At first, we consider the noise less case : α = 0.5 and T 0 = 0.0.In this case, the threshold c min is nearly equal to 0.816.Figure.10(a) shows that the results of the three algorithms almost converge to the RS solutions c replica ≃ 0.864 as we increase N.These practical algorithms can lead to the candidate solutions exceeding c min .Next, we consider the noisy case: α = 0.57 and T 0 = 0.05.The threshold c min is nearly equal to 0.756.Figure .10(b)exhibits that these practical algorithms can accomplish c min in the noisy case.

Performance evaluation
We evaluate the performance of the ARA with the initial candidate solution attained by the practical algorithms.We adopt the AMP algorithm to prepare for the initial candidate solution.The experimental settings are the same as those in Fig. 3.We set the RA parameter as λ = 0.6.At first, we perform the AMP algorithm for each instance.We utilize the final result as the initial candidate solution.We call this setting "AMP init.".For the same instance, we randomly generate the initial candidate solution from Eq.( 19) with a fixed fraction c, which is the same as one obtained by the AMP algorithm.We call this setting "random init.".
We plot the magnetization for two initializations in Fig. 11.The error bar is given by the standard deviation.We consider two cases α = 0.5 and T 0 = 0 in Fig. 11(a), and α = 0.57 and T 0 = 0.05 in Fig. 11(b).The dashed curves represent the RS solutions with c replica .In the "random init."setting, the numerical results are consistent with the RS solutions.In this setting, the first-order phase transition can be avoided by the ARA with and without noise.In the "AMP init."setting, the numerical results do not match the RS solutions.The ARA  can not exclude the first-order phase transition.In Fig. 12, we plot the histogram of the magnetization at s = 0.99 used in Fig. 11 to check the existence of the first-order phase transition in detail.In the "random init."setting, only one peak exists at m ≃ 1.In the "AMP init."setting, there are two peaks around m ≃ 1 and m 1.The ARA can not eliminate the firstorder phase transition even though the fraction c obtained by the AMP algorithm exceeds the threshold c min .For the other practical algorithms, the similar behaviors seem to occur.The deviations between the numerical results of the ARA in the "AMP init."setting and the RS solutions are due to the assumption of the probability distribution of the initial candidate solution in our replica analysis.We assume that the probability distribution of the initial candidate solution follows Eq. (19).Since the initial candidate solution obtained by the AMP algorithm is not generated from Eq.( 19), we can not directly apply our analytical results to the "AMP init."setting.
Finally, we consider why the ARA in the "random init."setting can avoid the first-order phase transition and the ARA in the "AMP init."setting can not.In the "AMP init."setting, the initial candidate solution depends on the original signal, the received signal and the spreading codes.Therefore, the initial candidate solution is correlated with the received signal and the spreading codes.Meanwhile, in the "random init."setting, the initial candidate solution depends only on the original signal and does not depend on the received signal and the spreading codes.Originally, the information about the original signal can be attained through the received signal and the spreading codes.In the "random init."setting, the initial can-didate solution has the information about the original signal not included in the received signal and the spreading codes.The initial candidate solution increases the effective pattern ratio.As we increase the pattern ratio, the free energy barrier gets smaller.Thus, the ARA in the "random init."setting can avoid the first-order phase transition.In the "AMP init."setting, the initial candidate solution only has the same information about the original signal obtained through the received signal and the spreading codes.Since the effective pattern ratio is the same as the original one, the free energy landscape in the "AMP init."setting is the same as the original one.Consequently, the ARA in the "AMP init."setting can not eliminate the first-order phase transition even if the candidate solution obtained by the AMP algorithm exceeds c min which is attained from the saddle-point equations in oracle cases.To analyze the performance of the ARA in the "AMP init."setting appropriately, we should change the probability distribution of the initial candidate solution in our replica analysis.

V. CONCLUSION
We performed a mean-field analysis of the ARA for the CDMA multiuser detection.In the CDMA multiuser detection, the first-order phase transition is encountered in the intermediate pattern ratio.This first-order phase transition degrades the estimation performance.To avoid the first-order phase transition, we applied the ARA to the CDMA multiuser detection.Firstly, we considered the ARA in oracle cases where the initial candidate solution is randomly generated from the original signal with a fixed fraction c.The first-order phase transition can be avoided by the ARA if we prepare for the proper initial condition.In the ARA, the differences in the magnetization between the two local minima at the first-order phase transition were smaller than those in the vanilla QA.The prior information of the original signal avoids or mitigates the firstorder phase transition.To validate our analysis, we performed quantum Monte Carlo simulations.The numerical results were consistent with the RS solutions under the static approximation, except for the intermediate values of the annealing parameter.Although the RS solutions under the static approximation were invalid in these cases, the results obtained from the RS solutions that the ARA can avoid the first-order phase transition were consistent with the numerical results.The RS solutions under the static approximation can be useful for understanding the qualitative behaviors of the ARA.
Next, we considered the ARA in practical cases where we prepare for the initial candidate solution by the practical algorithms.We considered the three algorithms: SA, SQA, and the AMP algorithms.The fraction c obtained by these practical algorithms can exceed the threshold c min to avoid the first-order phase transition.To evaluate the performance of the ARA with the initial candidate solution attained by the practical algorithms, we performed the quantum Monte Carlo simulations.In the "AMP init."setting, we prepared for the initial candidate solution with the AMP algorithm.To compare with the "AMP init."setting, we considered the "random init."setting where the initial candidate solution was randomly generated from Eq.( 19) with a fixed fraction c, which was the same as one obtained by the AMP algorithm for each instance.The ARA in the "random init."setting can utilize the additional information about the original signal not included in the received signal and the spreading codes.Since the free energy barrier is removed by the additional information about the original signal, the ARA in the "random init."setting can avoid the first-order phase transition.Meanwhile,in the "AMP init."setting, the initial candidate solution was correlated with the received signal and the spreading codes.The initial candidate solution only had the same information attained through the received signal and the spreading codes.Because no additional information about the original signal existed, the effective free energy landscape was the same as the original one.Therefore, the ARA in the "AMP init."setting can not avoid the first-order phase transition.The ARA in the "AMP init."setting did not match the RS solutions, whereas the ARA in the "random init."setting matched.The deviations between the numerical result of the ARA in the "AMP init."setting and the RS solutions were due to the assumption of the probability distribution of the initial candidate solution in our replica analysis.
In the "AMP init."setting, the initial candidate solution was correlated with the received signal and the spreading codes.To incorporate the correlation with the received signal and the spreading codes into the initial candidate solution, we need to consider an equilibrium configuration governed by the Gibbs-Boltzmann distribution.Then, the free energy is constrained by the equilibrium configuration.The equilibrium property of the constrained free energy can be analyzed by the Franz-Parisi potential, which is developed to study the metastable state structure for discontinuous mean-field spin glasses [73][74][75][76][77].In a future study, we will analyze the Franz-Parisi potential for the CDMA to investigate the performance of the ARA in practical cases properly.
Although we can not directly apply our theoretical results to the practical cases, we showed that the ARA in the "random init."setting can avoid the first-order phase transition for the CDMA multiuser detection.We exhibited that the probability distribution of the initial candidate solution was crucial in the ARA.Our results indicated that the effective free energy landscape did not change by the ARA with the initial candidate solution obtained by some practical algorithms.In practical cases, the ARA did not enhance the estimation per-formance for the CDMA multiuser detection.The similar behaviors would occur when we apply the ARA to combinatorial optimization problems.
The local stability condition of the RS solutions under the static approximation is expressed as This condition corresponds to the AT condition [78] in the ARA.This result is consistent with the previous result in Ref. [5] for the classical limit s = 1 and λ = 1.We can attain this condition by taking into account the perturbations to the RS solutions [79,80].The detailed calculations for deriving the AT condition in Eq. (A10) are presented in [66].The global instability condition of the RS solutions is related to the negative entropy.The existence of the global instability corresponds to the freezing behavior [81].To detect the freezing behavior, we calculate the RS entropy as follows:

FIG. 1 .FIG. 2 .
FIG. 1. Phase diagram of CDMA model in the ARA with λ = 1.The horizontal axis denotes the pattern ratio.The vertical axis represents the annealing parameter.The experimental settings are (a) T 0 = 0, (b) T 0 = 0.05, and (c) T 0 = 0.1.The "spinodal 1" and "spinodal 2" curves indicate the solutions from the two different branches.The "critical" curve denotes the point at which the RS free energy takes the same value.

FIG. 5 .FIG. 6 .
FIG. 5. Phase diagrams of CDMA model in ARA for four different values of c.The horizontal axis denotes the RA parameter.The vertical axis denotes the annealing parameter.These curves represent the points at which the first-order phase transitions occur.The experimental settings are (a) α = 0.6 and (b) α = 0.5.

FIG. 7 .FIG. 8 .
FIG. 7. Phase diagrams of CDMA model in ARA for four different values of c.Both axes are the same as those in Fig. 5.The experimental settings are (a) α = 0.62 and (b) α = 0.57.

FIG. 9 .
FIG.9.Dependence of order parameters on annealing parameter with λ = 0.8.Both axes are the same as those in Fig.3.

FIG. 10 .
FIG. 10.The dependence of the fraction of the ground state in the estimated signal obtained by SA, SQA and the AMP algorithm on the system size.The experimental settings are as follows: (a) α = 0.5 and T 0 = 0.0 and (b) α = 0.57 and T 0 = 0.05.

FIG. 11 .FIG. 12 .
FIG. 11.Dependence of magnetization on annealing parameter with λ = 0.6 in the "random init." and the "AMP init."settings.The dashed curve for each case is obtained by the saddle-point equations.The experimental settings are as follows: (a) is α = 0.5 and T 0 = 0, and (b) is α = 0.57 and T 0 = 0.05.Both axes are the same as those in Figs.3(a). 1