Continuous-variable quantum approximate optimization on a programmable photonic quantum processor

Variational quantum algorithms (VQAs) provide a promising approach to achieving quantum advantage for practical problems on near-term noisy intermediate-scale quantum (NISQ) devices. Thus far, most studies on VQAs have focused on qubit-based systems, but the power of VQAs can be potentially boosted by exploiting infinite-dimensional continuous-variable (CV) systems. Here, we implement the CV version of one VQA, a quantum approximate optimization algorithm by developing an automated collaborative computing system between a programmable photonic quantum computer and a classical computer. We experimentally demonstrate that this algorithm solves the minimization problem of simple continuous functions by implementing the quantum version of gradient descent to localize an initially broadly-distributed wavefunction to the minimum. This method allows the execution of a practical CV quantum algorithm on a physical platform. Our work can be extended to the minimization of more general functions, providing an alternative to achieve the quantum advantage in practical problems.


I. INTRODUCTION
Variational quantum algorithms (VQAs) have recently emerged as the leading approach to achieving quantum advantage for practical problems under the constraints of near-term noisy intermediate-scale quantum (NISQ) devices [1,2].In VQAs, such constraints are avoided by the common strategy to repeatedly run shallow-depth quantum circuits with the circuit parameters updated by classical optimizers.This strategy enables us to mitigate the accumulation of errors and fully exploit the computational space offered by the limited-scale devices.Thus far, a wide variety of VQAs have been proposed theoretically for qubit-based systems, such as ones for combinatorial optimization [3], chemistry simulation [4], and machine learning [5].They have already been demonstrated experimentally on several physical platforms [4,[6][7][8][9][10].
In contrast, there have been much fewer proposals [11][12][13][14][15] and no experimental implementations on continuousvariable (CV) VQAs, although CV quantum computing can potentially offer superior computational power in the NISQ era.The potential of CV systems lies in the ability to process infinite-dimensional quantum information even on single-mode devices, while in qubit-based systems each qubit provides only two-dimensional computational space.Furthermore, CV systems natively and efficiently handle continuous real parameters that often appear in real-world problems.In general, fully exploiting such infinite dimensionality and continuous degree of freedom in CV systems for quantum computation has been regarded as impractical due to their noise sensitivity and difficulty in error correction [16].However, this can be in turn a promising approach to extracting high computational power in the NISQ era when the error correction is not assumed.
In this article, we implement the CV version [11] of one of the most typical VQAs, a quantum approximate optimization algorithm (QAOA) [3] by developing an automated collaborative computing system between a programmable photonic quantum computer and a classical computer.We choose the simplest problem for the CV-QAOA and experimentally demonstrate it, where this algorithm minimizes one-variable continuous quadratic functions by working as the quantum version of gradient descent and localizing an initially broadly-distributed wavefunction to the minimum.The algorithm is shown to robustly find approximate answers to the problem with a noisy shallow-depth quantum circuit, thus confirming that a VQA works also in CV systems.This method showcases a CV quantum algorithm that can be used to solve practical problems, except for Gaussian boson sampling [17][18][19], which is designed for achieving quantum supremacy and has been partially linked to some practical problems [20,21].We also show that our implementation can be extended to the minimization problems of arbitrary-order functions by adding optical resources [22].It can also be extended to multivariable functions by incorporating multi-mode interactions.Thus our work highlights a new approach using CV quantum computing in the NISQ era, offering an alternative to achieve the quantum advantage in practical problems.

II. THEORY OF CV-QAOA
Many practical problems in various fields come down to optimization or minimization, which often require high computational costs for classical computers.The QAOA is a heuristic algorithm that could potentially offer a quantum speed-up to solve such problems on NISQ devices [3].Theoretical aspects and experimental imple-mentations of the QAOA have been recently studied intensively on qubit-based systems to solve discrete combinatorial optimization problems [6-9, 23, 24].
Later, a CV version of the QAOA was proposed to solve continuous optimization problems on CV systems [11].This algorithm is designed for minimizations of continuous real-valued functions, which have many practical applications in finance [25], machine learning [26], and engineering [27].This proposal [11] indicates that the CV-QAOA has potential of a quantum speed-up as its circuits can encode CV Grover's search algorithm [28], which achieves a quadratic speed-up over the classical algorithms.
The CV-QAOA is formulated as follows.The goal of the algorithm is to find an approximate minimum of a real-valued function f (x) are position and momentum of the particle, respectively.The initial state is |p = 0⟩, which is an eigenstate of p and thus equally-weighted superposition of |x⟩ for all x.The unitary operator given by where ĤC = f ( x) and ĤM = p2 /2, transforms the initial state to the final state.P determines the circuit depth.ĤC and ĤM are called cost and mixer Hamiltonians, respectively, and η = (η 1 , η 2 , • • • , η P ) and γ = (γ 1 , γ 2 , • • • , γ P ) are the tunable real positive parameters.The final state is then measured in the x-basis, the measurement outcome x being a candidate of the minimum of the function.This algorithm can be considered as a quantum version of the gradient descent method; a pair of the cost and mixer Hamiltonians transforms x as and indeed Û (η, γ) represents the Trotterized approximation of the time-evolution of a particle trapped in a potential of f ( x) in an N -dimensional space [3,11].Thus, while the distribution of x is initially uniform, it moves under the influence of the potential f ( x) and then localizes around the minimum if the parameters (η, γ) are properly chosen.
For the demonstration of the algorithm, we adopt the simplest problem setup, minimizing a quadratic function of one variable with the shallowest depth.Specifically, we choose f (x) = (x − a) 2 with a ∈ R (N = 1) and set P = 1.

III. PHOTONIC-CIRCUIT IMPLEMENTATION
We implement the CV-QAOA on a programmable photonic quantum computer, where optical amplitude and phase are identified with position and momentum in the algorithm, respectively.Our implementation is enabled by recent technological advances in photonic CV quantum computing, including state preparations, gate operations, and measurements [29].Especially, programmable and multi-step CV gate operations have been demonstrated very recently, which are indispensable for implementing the CV-QAOA [30][31][32].Such previous studies have been limited to the proof-of-principle demonstrations of predetermined quantum gates.Here, we implement the CV-QAOA by developing an automated collaborative computing system between such a programmable CV photonic quantum computer and a classical computer, where the latter assesses the output of the former and automatically updates the program of the former in real time.

A. Concept of our implementation
Our implementation can be conceptually shown in Fig. 1(a).The input state is a p-squeezed state, which approximates |p = 0⟩.After the parametrized operation of the cost and mixer on the input state, the output state is measured in the x-basis.According to Eq. ( 2), Û (η, γ) in our specific problem setting transforms x as (x in , pin ) and (x out , pout ) being the quadrature amplitude of the input and output states, respectively.The parameters η and γ are iteratively updated according to the sampling results of the output state.
To realize the concept of Fig. 1(a), we design and implement an optical and electrical setup shown in Fig. 1(b).To perform a gate operation of Eq. (3) in a programmable way, we use a measurement-induced gate composed of an ancillary state, measurement, and feedforward.In fact, this setup can be regarded as the simplest version of a more general setup to implement the CV-QAOA for arbitrary-order functions [22], as we show in Appendix B. The following is a specific description of the current setup.The input state is produced by an optical parametric oscillator (OPO) named OPO-1.The wavefunction of the state is squeezed in the p-direction in the sense that ⟨x 2 in ⟩ = e 2r+ /2 and ⟨p 2 in ⟩ = e −2r− /2 with r + , r − > 0. The ancillary state from OPO-2 is an orthogonally squeezed state with ⟨x 2 a ⟩ = e −2r− /2 and ⟨p 2 a ⟩ = e 2r+ /2.These two fields interfere at the beam splitter having variable transmissivity T and then they are sent to two homodyne detectors.The homodyne detectors measure the quadrature amplitudes x1,θ = x1 cos θ + p1 sin θ and x2,ϕ = x2 cos ϕ + p2 sin ϕ.Here (x 1 , p1 ) and (x 2 , p2 ) are the quadrature amplitudes of the upper and lower beams coming out of the variable beam splitter, respectively.The parameters of the processor T , θ, and ϕ can be varied via the applied voltage to the corresponding electro-optic modulators (EOMs).Then x1,θ is fed forward to x2,ϕ with a certain gain.To let The circuit parameters (T, θ, ϕ), the feedforward gain g, and the constant displacement x d are determined according to a set of (η, γ) so that the measurement-induced gate is performed in a consistent way.EOM, electro-optic modulator; PBS, polarizing beam splitter; QWP, quarter-wave plate.
the above-described circuit act as the gate designated by (η, γ) as Eq. ( 3), the parameters of the optical circuit are determined by and tan ϕ = γ, and correspondingly the feedforward gain is set by g = γ 2 − 4ηγ + 4η 2 (1 + γ 2 ), in which settings the contribution of the anti-squeezed quadrature of the ancilla pa to the gate output is canceled so that the measurement-induced gate operates properly.The constant displacement of x d = 2aηγ is also performed for the constant term in Eq. ( 3).The feedforward and constant displacing are numerically done by post-processing in the classical computer, which is justified because it gives the same results as those given by optical displacing.In fact, feedforward operations were performed by postprocessing in recent demonstrations of one-way quantum computation [30,31].After the above classical postprocessing, the gate output xout becomes which asymptotically coincides with Eq. (3) in the high squeezing limit r − → ∞ in the sense that the variance of the noise term −2ηx a approaches zero (see Appendix A for the derivations).Therefore, our photonic processor depicted in Fig. 1(b) is capable of obtaining the xmeasurement result of the output state, xout , with the gate parameters η and γ varied.

B. Experimental setup
Here we describe the details of our experimental setup in Fig. 1(b).We use a continuous-wave laser of wavelength 1545.3 nm.Two OPOs are pumped by the second harmonic fields with wavelength of 772.7 nm.The pump power is set to 200 mW.The full width at the halfmaximum of the OPOs is 60 MHz.The variable beam splitter is composed of a bulk electro-optic modulator named EOM-1, a quarter-wave plate, and a pair of polarizing beam splitters.EOM-1 serves as a variable polarization rotator and thus works as a variable beam splitter with polarization optics.We inserted the quarter-wave plate so that the transmissivity is 50 % when no voltage is applied to EOM-1, which makes it easy to lock the relative phase between the input and the ancillary beams.
Each homodyne detection is performed by interfering the local oscillator field with the signal field at a 50:50 beam splitter.Two beams from the beam splitter are received by two photodiodes, the photocurrents of which are subtracted with each other and amplified in the electric circuit.The bandwidth of the circuit is about 200 MHz.The optical power of the local oscillator field is set to 10 mW.A fiber-coupled electro-optic modulator EOM-2 or 3 shifts the optical phase of each local oscillator for the control of the homodyne angle θ or ϕ.
The outcome of the homodyne detection is acquired by an oscilloscope and then sent to the classical computer.The time series from the oscilloscope is converted to a quadrature amplitude by convoluting it with a mode function h(t) defined by where Γ = 3 × 10 7 /s and t 1 = 50 ns.The purpose of using this mode function is to eliminate undesirable effect from low-frequency electrical noise from homodyne detectors [33].Based on the measured quadratures and the subsequent analysis, the classical computer can automatically reprogram the photonic quantum circuit by changing the voltages applied to EOM-1, 2, and 3.This enables the classical computer to collaborate with the photonic quantum computer updated in real time and perform the CV-QAOA.As a preliminary measurement, the outputs of OPO-1 and OPO-2 are measured by the homodyne detectors with the transmissivity of the variable beam splitter set to zero.The squeezing level and the anti-squeezing level of these modes are measured to be −5.3 dB and +9.0 dB on average.This measurement result indicates that the overall optical loss of the experimental setup is estimated to be 22 %.

IV. OPTIMIZATION LANDSCAPE AND ALGORITHM PERFORMANCE
The CV-QAOA is performed on the processor as follows; first, we repeatedly run the circuit and sample xout with the parameters η and γ fixed to calculate the mean value ⟨f (x out )⟩; then, an outer-loop optimizer in the classical computer suggests new parameters to decrease ⟨f (x out )⟩.These two stages are repeated alternately.To experimentally demonstrate the capability of the above, we operate the whole system in three different conditions for the parameter update: (i) the parameters η and γ are scanned like a grid search; (ii) the parameters are fixed at their optimum; (iii) the parameters are updated according to the protocol of the Bayesian optimization.Each result is shown in the following.
First, to qualitatively diagnose that the processor operates properly for the full search range of the outerloop optimizer, ⟨f (x out )⟩ is evaluated with the parameters η and γ scanned like a grid search.The landscapes of ⟨f (x out )⟩ as a function of η and γ are obtained from the experiment or the numerical simulation as shown in Fig. 2. The landscape structure of the experiment reasonably well agrees with that of the simulation.Specifically, ⟨f (x out )⟩ is small around ηγ = 1/2 and the smallest around (η, γ) = (1/2, 1), which is the theoretical optimum in the high squeezing limit indicated by white stars in the figure.Note that the optimal point in the finite-squeezing case is , where δ = e −2r− /(e 2r+ + 2a 2 ), and thus the deviation from (1/2, 1) is less than 2 % with the values of r + and r − for our experiment no matter the value of a.The overall landscape structure is also independent of the value of a because it only affects the steepness of the valley of ηγ ∼ 1/2.
Next, to visualize how the algorithm works like the gradient descent method to reshape the wavefunction of x, we run the processor with the parameters fixed at the optimal point (η, γ) = (1/2, 1). Figure 3 shows the histogram of the sampled xout for a specific a.The sampling results of the input state xin are also plotted.It can be found that the input state has a broad distribution, and the optimal gate operation localizes the distribution around a, the exact solution.In fact, the standard deviation of xout is 0.5767 (4), smaller than that of the vacuum state 1/ √ 2 thanks to the quantumness of the processor.Also, the gap between the mean output and the exact x out = a 2XWSXWRSWSDUDPV ,QSXWx out = xin solution is ⟨x out ⟩ − a = −7.1(5)× 10 −3 , the absolute value of which is much smaller than the standard deviation.This indicates that the systematic shift from the exact solution is negligible with respect to the statistical broadening of the distribution.These results visually prove that the experimentally implemented gate for Û (η, γ) in our problem setting properly localizes the distribution around the minimum like the gradient descent method when the parameters are optimal.
Finally, to evaluate the performance of the CV-QAOA in a realistic condition where the optimum of the parameters is unknown, we perform the algorithm with the parameter updated by the Bayesian optimization.The optimizer suggests new parameters every 1000 samples of xout so that ln⟨f (x out )⟩ is minimized (see Appendix C).We repeat such a process many times with the value of a changed to see the statistical behavior of the algorithm (see Appendix D).In Fig. 2, the overlaid red dots show the distribution of the classically optimized pair of (η, γ).Each dot corresponds to the pair of the parameters that gives the smallest ⟨f (x out )⟩ among 100 suggestions by the optimizer for each execution of the CV-QAOA.These figures show that the classical optimizer reaches the point around the optimum within the 100 steps as the simulation predicts.Figure 4(a) shows the typical trace of the parameter update by the Bayesian optimization, where a wide area is initially explored but gradually the search becomes concentrated around the optimum.The behavior of ln⟨f (x out )⟩, which is the target function of the Bayesian optimization, is shown in Fig. 4(b).The traces denote the average of ln⟨f (x out )⟩ for all the CV-QAOA trials.The observed decreasing trend of the average of the target function is comparable to the simulation expectation.The small differences between the experimental and simulated traces can be attributed to imperfections of the optical setup such as intensity fluctuations and/or alignment drifts of the local oscillator beams for the homodyne detectors.
To quantify the performance in minimizing the function, the success probability of finding the minimum is evaluated experimentally and numerically.Figure 4(c) shows the cumulative success probability of finding the minimum of f (x) up to a certain step of the parameter updates.Here we set the criterion of the success by f (x out ) = (x out − a) 2 < 1 × 10 −9 .The success probability is calculated by the success frequency for 30 different values of a in f (x) = (x−a) 2 .The solid line is derived by averaging the success probability for eleven sets of the trials while the shaded area denotes the ±1σ-region around the average derived from the eleven trials.These results show that the experimental results of the increase of the cumulative success probability coincide with the numerically simulated ones.They also show that the success probability obtained by the CV-QAOA is significantly better than that by random sampling, where the input state is directly measured.

V. DISCUSSION
In conclusion, we demonstrate the successful implementation of the CV-QAOA with the parameters of the programmable photonic quantum processor updated.The demonstration experimentally shows that the performance of the algorithm for the minimization of quadratic functions is significantly better than that of the random sampling and comparable to what the numerical simulation predicts.Even though the experimental system is influenced by various imperfections including optical loss, the CV-QAOA still successfully finds the minimum of the function, indicating that the algorithm works robustly.
Such imperfections limit the effective squeezing levels.The optimum effective squeezing level for the overall performance is nontrivial.In fact, with higher squeezing levels, the width of the output distribution in Fig. 3 gets narrower and this increases the probability of finding the optimal point of x.On the other hand, it is also observed in our numerical simulation that the parameters η and γ are optimized more slowly with higher squeezing levels due to the landscape change in Fig. 2. The former (the latter) is advantageous (disadvantageous) to the algorithm.Optimum squeezing level should be investigated, but its detailed analysis is left for future work.
Our experiment in this paper can be simulated efficiently with classical computers because our system is entirely built with Gaussian building blocks [34].However, we show that our implementation can be extended to arbitrary-order functions by adding non-Gaussian ancillary states other than squeezed states [22] (Appendix B).This may push the CV-QAOA to the non-Gaussian regime beyond efficient classical simulation. .It can also be extended to multivariable functions by using beam splitters for multi-mode interactions.In this way, more complex functions can be minimized to address the practical problems.For such more complex functions, we may set P > 1 with a deeper circuit at the cost of the increased number of the classical parameters to be optimized (η, γ).The current system took ∼100 seconds to perform one CV-QAOA trial with 100 optimization steps.This runtime was mainly dominated by the measurement time, which can be shortened by increasing the bandwidths of the OPOs and the homodyne detectors.
This work is the first experimental realization of a quantum algorithm using CV information, which proves the usefulness of CV quantum systems for natively solving CV problems.In general, CV problems can also be solved with qubit-based quantum computers [35], but it requires more resources to represent discretized CV parameters with many qubits and also introduces quantization errors.CV quantum computers can avoid such difficulties and efficiently encode CV problems.This demonstration sheds light on the advantages of CV quantum computing that its infinite dimensional space can be exploited in NISQ applications.It stimulates the realizations of other VQAs in CV systems such as quantum machine learning [12,13] and thus opens a new promising way toward quantum advantage.
If we set the feedforward gain by g ′ = (1 − T )/T / sin θ, the output described by Eq. (A3) becomes Note that g ′ is chosen so that the anti-squeezed quadrature pa disappears in this expression.This expression indicates that, apart from a noise term proportional to the squeezed quadrature xa , the circuit in Fig. 5 is capable of performing various linear transformations by changing T , θ, ϕ, and x ′ d .In fact, by setting ) which asymptotically coincides with Eq. (A1) in the high squeezing limit of xa → 0. Note that the upper row of this equation is identical to Eq. ( 4).
Given that we eventually measure the quadrature amplitude xout , the phase rotation by ϕ can be achieved instead by changing the homodyne angle for the measurement of the output state.In addition, the displacing operation can be replaced by numerical post-processing after the homodyne measurement as the displacing operation only shifts the mean value of the measurement outcome.For this reason, the circuit in Fig. 1(b) is equivalent to the one in Fig. 5 as long as the quadrature amplitude of the output state is measured.In the circuit in Fig. 1(b), xout is in fact provided using the two homodyne-measurement outcomes by x2,ϕ + gx 1,θ + x d , which corresponds to the upper row of Eq. (A3).Here we redefined g := g ′ sin ϕ = γ 2 − 4ηγ + 4η 2 (1 + γ 2 ) and x d := x ′ d sin ϕ = 2aηγ.

APPENDIX B: IMPLEMENTATION OF HIGHER-ORDER FUNCTIONS
Let us show that our setup can be extended to a quantum circuit for minimizing a higher-order function.As can be seen from Eq. ( 2), if the minimized function f (x) is an n-th order polynomial, the transformation by a pair of the cost and mixer Hamiltonians in the QAOA involves an (n − 1)-th order polynomial of x: where we express f (x) = n k=0 a k x k with real coefficients a k .The above transformation can be implemented by the circuit shown in Fig. 6.In contrast to Fig. 5, where one beam from the beam splitter is measured by simple homodyne detection, the beam is nonlinearly measured with the help of the ancillary states.The working principle can be made clear by dividing the circuit into two part: the nonlinear measurement followed by the feedforward, and the phase rotation by ϕ.
The combination of the nonlinear measurement and feedforward provides an (n − 1)-th order polynomial of x in the following way.The ancillary states , where (x k , pk ) are canonical pairs of the quadrature operators of the corresponding modes.The beam is sequentially coupled with the ancillary states A k by the beam splitters having transmissivity of T k .Then, one outgoing mode from each beam splitter is measured by homodyne detection in x quadrature, the outcome of which is labeled as q k .The transmissivity T k is adaptively controlled depending on the past outcomes q j , where k < j ≤ n.The rightmost homodyne measurement is done with the homodyne angle of θ, which is also controlled depending on the past outcomes.The quadratures after appropriate feedforwards from the nonlinear measurement block are expressed, in the limit that the ancillae are ideal, as where C k depends on q j (k < j ≤ n), T l (k ≤ l ≤ n) and θ (see Eqs.( 5) and ( 9) in Ref. [22]).Here, by recurrently determining T n , T n−1 , • • • , T 3 in the order from n to 3 and then θ, the coefficient C k (2 ≤ k ≤ n) can be arbitrarily chosen.As the coefficient C 1 corresponds to constant displacing for p, it is also arbitrarily determined at the displacing operation.Therefore, the above process realizes n-th order nonlinearity with arbitrary coefficients C k .
In this way, the transformation by a pair of the cost and mixer Hamiltonians for a given polynomial f (x) can be implemented.Thus, in the viewpoint of implementing generic higher-order functions, our experimental setup can be regarded as the simplest case of Fig. 6.

APPENDIX C: PARAMETER SEARCH CONDITION WITH THE BAYESIAN OPTIMIZATION
For the update of the circuit parameters (η, γ) in Fig. 1, we adopt the Bayesian optimization, which has been commonly used among the derivative-free optimization methods in the parameter optimization of the QAOA [6,7].The following conditions of the Bayesian optimization are set after checking the convergence of the circuit parameters in numerical simulations.log 10 η and log 10 γ are handed to the Bayesian optimizer as free parameters.The target function to be optimized is ln⟨f (x out )⟩.We use a python package for the Bayesian optimization [37].The Matérn kernel with ν = 2.5 is chosen as the covariance kernel function.The acquisition function is a type called the upper confidence bound given by µ + κ(t)σ, where µ and σ are the estimated mean and standard deviation of the target function, respectively.We set κ(t) = κ 0 × 0.97 t , where κ 0 = 2.576 and t is the number of optimization steps.
As for the QAOA of the qubit system, the parameter range of the unitary operations is often limited by [0, 2π).However, the CV case does not have such a limit.Since it is difficult to optimize parameters with unlimited search range, we limit the range by estimating the order of the optimum of the parameters using a generic prescription described below.The search range of the parameters is set by 0.1 < η < 10 and 0.1 < γ < 10 because the optimum is estimated to be (η opt , γ opt ) ∼ (1, 1).
Nonlinear measurement FIG. 6. Implementation of a higher-order function.The outgoing beam from the first beam splitter is nonlinearly measured with the help of the ancillary nonlinear phase states, and the measurement outcomes are fed forward before the phase rotation by ϕ.Instead of also preparing a quadratic phase state and an additional beam splitter at the right end of nonlinear measurement block as proposed in Fig. 1d in Ref [22], the homodyne angle θ of the rightmost homodyne measurement is made variable.In fact, the implementation of the fourth-order nonlinear phase gate is proposed in this way (Fig. 2 in Ref. [22]).
Here we explain the prescription for the estimation.Let us consider the CV-QAOA with P = 1 for a onevariable function f (x).If implemented by a combination of measurement-induced gates [22], a pair of the cost and mixer operations transforms the quadratures as where N is a noise term arising from non-ideal ancillary states.xa and pa formally denote the quadrature amplitudes of the ancillary states.Let us assume that the ancillary states are linearly or non-linearly squeezed sufficiently and thus |⟨N ⟩| ≪ |x min |, where x min = argminf (x).We can also assume that the optimal parameters, namely (η opt , γ opt ), localize the distribution of xout around x min , and thus provide two conditions: ⟨x out ⟩ ∼ x min and ⟨∆x 2 out ⟩ is minimized.By assuming that the optimal parameters (η opt , γ opt ) are adopted and taking the expectation value of Eq. (C1), the assumption ⟨x out ⟩ ∼ x min reduces to Here we assume ⟨x in ⟩ = 0 and ⟨p in ⟩ = 0 since the input state is the squeezed vacuum.In this way, if the prior information on the order of magnitude of x min and ⟨∇f (x in )⟩ is available, the product η opt γ opt can be inferred.Let us next consider the variance of xout .Since Eq. (C2) estimates the product of the optimal parameters, we evaluate the variance under the constraint of ηγ = c, where c denotes the estimated value of the product η opt γ opt .By substituting ηγ = c into Eq.(C1) to eliminate η, we have As the first term does not depend on γ, γ opt can be estimated by minimizing the sum of the second and third terms.The third term can be calculated if the specific form of N (η, γ; xa , pa ) is known.Once the third term is calculated, the sum of the second and third terms becomes a one-variable function of γ, which should be able to be minimized.Let us consider the specific case of f (x) = (x−a) 2 , and estimate the order of magnitude of the optimal parameters by using the typical magnitude of xin .We denote ⟨x 2 in ⟩ = σ, which we regard as the typical magnitude of xin .As the function f (x) is quadratic in x and its leading term is x 2 , the order of magnitude of the variation of f (x) in the range of ±σ can be estimated by σ 2 .We then estimate the gradient of the function by dividing the typical range of the function output σ 2 by the typical range of the function input σ: |⟨∇f (x)⟩| ∼ σ 2 /σ = σ.We can also assume that |x min | ∼ σ.Using these estimations, we can estimate η opt γ opt as η opt γ opt ∼ 1 (C4) from Eq. (C2).Let us then calculate the variance under ηγ = 1.Since N (η, γ; xa , pa ) = −2ηx a from Eq. ( 4), the sum of the second and third terms of Eq. (C3) becomes γ 2 ⟨p 2 in ⟩ + (4/γ 2 )⟨x 2 a ⟩.By minimizing this in terms of γ, γ opt can be estimated as In this way, the optimal parameters are estimated as (η opt , γ opt ) ∼ (1, 1), which is why we set the search range by [0.1, 10] ⊗2 with a margin of about a factor of ten.The above discussion is only a rough estimation of the order of magnitude, and although it works in our case, there is no guarantee that the optimal parameters exist in the search range.For this reason, a general strategy may be to perform the CV-QAOA within the initial search range and, if there seems to be the optimal point outside the search range, expand or change the search range by observing the behavior of the parameter optimization.

APPENDIX D: REPEATED EXECUTION OF THE CIRCUITS
In the demonstration of the CV-QAOA with parameters updated (Fig. 4), the iterations of the circuit execu-tion are hierarchical.For clarity, the conditions for that hierarchical execution are summarized here.For each pair of the parameters, we repeatedly run the circuit and sample xout 1000 times with the parameters fixed to obtain ⟨f (x out )⟩.The Bayesian optimizer suggests 100 pairs of the parameters by using the results of ⟨f (x out )⟩.Such 100-step classical optimization is repeated eleven times for the same value of a (in f (x) = (x − a) 2 ).Finally, this repeat is done for 30 different values of a.
The values of a are randomly sampled from the Gaussian distribution of zero mean and the standard deviation equal to 1.99, which corresponds to ⟨x 2 in ⟩.This is because we intend to mimic the following situation.Suppose that the range of x that gives the minimum of f (x) (x min ) can be roughly estimated by some conditions in the problem settings such as, for example, physical conditions or features of f (x).In this case, the distribution of the initial state can be set so that it covers the estimated range for x min .In this demonstration, mimicking the situation where the range estimation is correct and the solution a is in fact in the estimated range, we repeatedly execute the CV-QAOA with many different a sampled from that range to statistically evaluate the performance of the algorithm.Note that the range for x min can always be matched to the range of the initial state distribution by rescaling and translating x, and resultantly redefining f (x).Generally, if the estimation is incorrect, one can perform the algorithm again with an effectively broader initial state.

APPENDIX E: SIMULATION CONDITION
In the simulation, the optical circuit is numerically simulated by expressing the quadratures of the initial and ancillary states using Gaussian random numbers.Based on the measured squeezing level, we set ⟨x 2 in ⟩ = ⟨p 2 a ⟩ = 10 9.0/10 /2 and ⟨p 2 in ⟩ = ⟨x 2 a ⟩ = 10 −5.3/10 /2.The asymmetry between the squeezing and anti-squeezing implies that optical loss is taken into account.No other imperfections are included in the simulation.

FIG. 2 .
FIG. 2. Experimental and simulated landscapes of the CV-QAOA.The landscape structure of the experiment (left) agrees with that of the simulation (right).Each grid is evaluated by 1000 samples of xout.f (x) = (x − 1) 2 is used for the landscape evaluations.Each red dot on the landscapes shows the parameter set that records the smallest ⟨f (xout)⟩ during the 100-step Bayesian optimization.The white stars denote the theoretical optimum of (η, γ) = (1/2, 1).See Appendix E for the simulation condition.

FIG. 3 .
FIG. 3. Histogram of the output distribution.The normalized frequencies of the sampled xout in the experiment are shown.The green bars are the output distribution with the optimal parameters of (η, γ) = (1/2, 1) while the red ones are the distribution of the input state, which corresponds to (η, γ) = (0, 0).The number of the sampled xout is 1.1 × 10 6 for each.The curve of f (xout) = (xout − a) 2 is overlaid.We use a = 2.745 for these measurements.The comparison of two distributions demonstrates how the CV-QAOA works.

FIG. 4 .
FIG. 4. Algorithm performance with parameter updates.(a) Typical trace of the parameters updated by the Bayesian optimization.The initial search uniformly spans the full range but gradually the search becomes concentrated around the optimum, which is depicted by the dashed lines.(b) Convergence of the classical optimization.The solid traces show the average of the target function ln⟨f (xout)⟩ from the experiment and the simulation.The black dashed line is the theoretical values for the case the parameters are fixed (η, γ) = (1/2, 1).(c) Success probability to sample xout such that f (xout) < 1 × 10 −9 .The solid line is derived by averaging the probability for the eleven sets of such trials.The shaded area shows the ±1σ region around the average.The green and blue plots are the experimental and simulated results of the CV-QAOA, respectively.As a reference, the red and orange plots show the experimental and simulated results of random sampling.They indicate that the CV-QAOA finds the minimum of f (x) significantly more efficiently than the random sampling.

FIG. 5 .
FIG.5.Optical circuit for a measurement-induced linear transformation.The input state interferes with the ancilla at the beamsplitter having the transmissivity of T .One outgoing beam from the beamsplitter is measured by homodyne detection, and the measurement outcome is fed forward to the other beam.The quadrature p2 is displaced by g ′ x1,θ + x ′ d , and then a phase rotation by ϕ is applied.