Simple, reliable and noise-resilient continuous-variable quantum state tomography with convex optimization

Precise reconstruction of unknown quantum states from measurement data, a process commonly called quantum state tomography, is a crucial component in the development of quantum information processing technologies. Many different tomography methods have been proposed over the years. Maximum likelihood estimation is a prominent example, being the most popular method for a long period of time. Recently, more advanced neural network methods have started to emerge. Here, we go back to basics and present a method for continuous variable state reconstruction that is both conceptually and practically simple, based on convex optimization. Convex optimization has been used for process tomography and qubit state tomography, but seems to have been overlooked for continuous variable quantum state tomography. We demonstrate high-fidelity reconstruction of an underlying state from data corrupted by thermal noise and imperfect detection, for both homodyne and heterodyne measurements. A major advantage over other methods is that convex optimization algorithms are guaranteed to converge to the optimal solution.


I. INTRODUCTION
Quantum state tomography has a rich history, and interest in it is only increasing due to the current development of many quantum technologies such as quantum sensing, communication, cryptography and computing. The goal of quantum state tomography is to infer information about a quantum state from measurement data. Calculating probability distributions for observables given a specific quantum state is straightforward, but the inverse problem of calculating a quantum state given measured probability distributions is fraught with issues due to the problem being ill-posed, meaning slight variations or noise in the measurement data can heavily affect the result [1].
While nowadays multi-qubit state tomography is an area of great interest, the development of quantum tomography started with continuous variable (CV) states in 1989 when Vogel and Risken proposed that the Wigner function can be obtained from the inverse Radon transform of homodyne measurement data [2]. Ever since the first experimental quantum tomography was performed using such measurements [3], homodyne tomography has been a mainstay in the quantum optics community. At the beginning, reconstruction techniques had the issue that the resulting states could be unphysical [4]. New numerical methods were developed to handle quantum state reconstruction from homodyne measurements, ranging from pattern function methods in the mid 90s [5,6] to maximum likelihood estimation (MLE) emerging in the late 90s [7,8]. An iterative MLE algorithm was introduced in the 2000s [9], and it has remained very popular over time, with different variations offering improvements being developed [10][11][12][13]. Currently, even more advanced methods are emerging. Bayesian methods have been proposed [14][15][16], and neural networks for CV state tomography have started to appear [17][18][19][20]. There are downsides to all above mentioned methods: Bayesian methods require the use of Monte Carlo sampling algorithms, iterations of an MLE algorithm have to be terminated at some arbitrary point, and neural networks are practically black boxes where the result can depend strongly on the chosen cost function.
Quantum state tomography is an essential tool in the field of quantum information. Not only is it needed to characterize states and verify state preparation, having a reliable method of state tomography also plays an important role in several methods for quantum process tomography, which characterizes a quantum gate or channel [21][22][23][24][25][26].
Here, we demonstrate a method for continuous variable state reconstruction that is both conceptually and practically simple, based on convex optimization. Convex optimization comprise a special class of mathematical optimization problems that can be solved numerically very reliably and efficiently [27]. It has been applied to process tomography [28,29], detector tomography [30,31] and discrete variable state-tomography (though often assuming low-rank density matrices) [32][33][34][35]. Still, it seems to have been overlooked for CV state tomography. We use it to demonstrate high-fidelity reconstruction of a CV quantum state with no assumptions of the rank, from data contaminated by thermal noise or imperfect detection. The primary advantage of this method compared to iterative maximum likelihood and neural network methods is that convex optimization algorithms are guaranteed to converge to the optimal solution.
The paper is structured as follows. We begin by introducing the basic concepts of quantum states and measurements in Section II, and how their properties are suitable for forming the tomography problem as a convex optimization problem. In Section III our state tomography method is demonstrated first on simulated homodyne data and then heterodyne data. For heterodyne measurement, its performance is compared to the ubiquitous maximum likelihood estimation and a new neural network tomography method. We also test our method on real experimental data with good results. Finally, we conclude and summarize in Section IV.

II. QUANTUM STATE TOMOGRAPHY
A. The quantum state The state of a quantum system can be described by a N ×N density matrix ρ in a Hilbert space, where N is the dimension of the system in question. For a CV state, ρ is often presented in the Fock (number state) basis. In theory the corresponding state space is infinite dimensional, but for practical reasons one may truncate the Hilbert space to a finite dimension that is sufficient to contain all non-zero elements of the density matrix [36]. There are a few conditions a density matrix must fulfill in order to describe a physical system: it must be Hermitian, have unit trace, and have no negative eigenvalues [37].

B. Measurements
A general quantum measurement can be described by a set of operators {Π k } called positive operator valued measure (POVM), where each operator is associated with a possible measurement outcome k [38] such that the probability of obtaining k is [39] p k = Tr[Π k ρ]. (1) When the density matrix ρ is reshaped into a vector ρ (see Appendix A), the tomography problem can be expressed as a linear equation where the state ρ is the unknown to be solved for. The elements b k of b are the measurement records p k [40]. here placed in a onedimensional array. We use the Fock basis with operators Ω i×N +j = |i j| , i, j = 0, . . . , N − 1.
The density matrix is then represented as ρ = The matrix A is given by the M × N 2 matrix C. The tomography problem Direct linear inversion of Eq. (2) to obtain the density matrix elements ρ k is not feasible since it does not guarantee positive semidefiniteness of the state, so any noise in the measurement data b can result in an unphysical density matrix [41]. Also, in this case it is likely that no solution exists at all if the system is overdetermined (M > N ). For this reason it makes sense to determine an approximate solution that renders the residual A ρ − b as small as possible. This means the problem of finding a solution to the set of linear equations A ρ = b can be formulated as the optimization problem When the vector norm · is chosen to be the Euclidean or 2 -norm [42], this is well-known as the least-squares problem [43]. Unlike the linear problem, an analytical solution can be obtained for this problem. However, the ill-posedness of an inverse problem in the form of a linear system of equations (2) is reflected in an ill-conditioned matrix A. This means an analytical solution by the normal equations is numerically unstable due to the high condition number of A, and the result is highly sensitive to noise in the measurement results [44]. A stable solution to an ill-posed problem can be obtained via numerical regularization methods [45]. Typically this entails adding a term to the expression (5) that is optimized. In our case this would still not be adequate since an unconstrained optimization does not guarantee a physical state. Fortunately, the quantum tomography problem has a number of nice features that we can utilize as described below.

D. Convex optimization
The set of all quantum states is the set of Hermitian matrices with non-negative eigenvalues and unit trace. This is a convex set [46]. The combination of being Hermitian and having non-negative eigenvalues means that a density matrix describing a physical quantum state is positive semidefinite (ρ 0). As such, we can write our optimization problem as a semidefinite program subject to ρ 0, (7) Tr ρ = 1.
The constraints that ensure a physical density matrix also act as regularization, meaning that the problem is now well-posed [47][48][49].
Since the feasible solution set {ρ} is convex, the cost function (6) is a convex function over this set, and the trace equality constraint is linear, it is a convex optimization problem. This means we can utilize efficient convex optimization methods. Convex optimization problems are easier to solve than general nonlinear optimization problems, but most importantly, it can be shown that every local minimum to a convex optimization problem is also a global minimum, guaranteeing an optimal solution [27].
The least-squares problem (6) will be the cost function for our convex optimization. The optimization variable corresponding to the unknown state is defined as an N × N Hermitian matrix X; it is only vectorized when calculating the cost, as indicated by the vector notation ρ in Eq. (6) as opposed to ρ in Eqs. (7) and (8). Having the variable in matrix form makes stating the constraints X 0, Tr X = 1 (9) very straightforward. There are efficient numerical methods that solve these types of problems, and software packages that provide them. We will use the open source Python package cvxpy [50][51][52]. Being able to utilize already existing software makes this state reconstruction procedure very easy to implement.

III. APPLICATION WITH COMMON MEASUREMENT SCHEMES
Homodyne and heterodyne detection are practical experimental techniques for investigating the quantum properties of light. Below we test the convex optimization state reconstruction method on data from these two measurement schemes, starting with simulated homodyne measurements having different levels of efficiency in section III A. In section III B we move on to heterodyne measurement. First we test the method on simulated noisy data in subsection III B 1. Then we compare the performance of our method to a neural network and a maximum likelihood method in subsection III B 2, using simulated, ideal data as the available code for those two methods cannot handle noisy data. Finally, we use our method on experimental heterodyne data and compare it to a maximum likelihood method based on moments in subsection III B 3. Another example of a reconstruction using experimental data in the form of Wigner tomography can be found in Appendix B.

A. Homodyne measurement
In an optical homodyne experiment, a current proportional to the generalized field quadrature x θ = (a † e iθ + ae −iθ )/ √ 2 is measured, where a and a † are the bosonic field annihilation and creation operators, respectively. The rotation angle θ is set by the phase of a local oscillator. The full manifold of quadratures can be obtained by varying θ between 0 and π [53]. The underlying quantum state of light can be reconstructed from a number of such measurements [54]. Here we test the convex reconstruction method with simulated homodyne measurement data generated by solving a stochastic master equation as implemented in Ref. [55]. Simulated homodyne currents are integrated and the data is discretized by sorting it into bins, here denoted by an index j. Fock basis matrix elements for the measurement operator corresponding to rotation angle θ and bin j are where the integral from x j to x j+1 is over photocurrents in bin j, and the overlap between the number and quadrature eigenstates is given by the harmonic oscillator eigenfunctions in the position basis multiplied by a phase [56]: with H n being the nth Hermite polynomial. We use a test state because it is a state that is asymmetric in phase space and has Wigner negativity, and it could also be prepared with only small modifications of the code in Ref. [55].

State reconstruction with noisy homodyne data
Noise and losses are an inevitable part of any experiment. Common sources are amplification noise and nonideal detection efficiency, and these types of imperfections can be related to each other as shown in Appendix C. In traditional quantum optics literature [57] the dominant imperfection is the measurement efficiency η of homodyne detection, which is in realistic cases typically less than 1, where η = 1 would correspond to 100 % efficiency. Inefficiency in homodyne tomography can be accounted for in the reconstruction by modifying the measurement operator as [9] Π = m,n,k B m+k,m B n+k,n Π mn (θ, j) |n + k m + k| .
Using this, our reconstruction of the test state (12)  well even with very inefficient measurements, with fidelities over 0.9 for efficiencies as low as η = 0.3 (see Table I), where the fidelity between the reconstructed state ρ and the ideal state ρ ideal is defined as The simulated measurement used 20 rotation angles, with data binned into 20 bins per angle. Due to fluctuations in the simulated data which consists of integrated currents for 2000 stochastic trajectories [58], the values are averaged over six measurement rounds. A reconstruction of the test state from data with efficiency η = 0.5 is displayed in Fig. 1.
To the best of our knowledge, there is unfortunately no publicly available code for other reconstruction methods using noisy homodyne data to compare with. The code and simulated data used to produce these results are provided in Ref. [59].

B. Heterodyne measurement
Noiseless heterodyne measurement corresponds to measuring the Husimi Q-function. The Q-function for a state ρ is defined as [60] The corresponding measurement operators are projections onto coherent states Π k = |α k α k | where α k ∈ C denotes a point in phase space. We test the method on a coherent superposition of coherent states |β , also called Schrödinger's cat state: It is an interesting looking state [cf. Fig. 2(c)] of interest for error-corrected quantum computing [61]. The code used to produce the results below are provided in Ref. [59].

State reconstruction with noisy heterodyne data
As an example of a noisy system, consider superconducting circuits, which provide a promising platform for quantum computing experiments. Since those circuits must be cooled to extremely low temperatures, thermal noise is a disruptive factor. It can be generated at room temperature and transmitted to the chip, but additionally, amplification of a weak measurement signal inevitably adds noise. This noise can be modeled as a thermal state [62]. When a measurement contains amplifier noise described by a thermal state ρ th n with mean number of photons n, the measured histogram corresponds to the convolution [63] where Q(β) corresponds to the Husimi Q-function [64] of the underlying wanted state, and is the P -function of the thermal state [65]. The noisecompensated measurement operators are given by [66] Π For testing, simulated noiseless data was obtained by calculating the Q-function of the test state at discrete points on a grid in phase space. This was then convolved with the thermal P -function to generate simulated noisy histograms. Using the noise-compensated operators (20) to construct the matrix A defined in Eq. (4), we show in Fig. 2 that the underlying quantum state can be reconstructed with very high fidelity from noisy data. The clean histogram is displayed in Fig. 2(a) for comparison to the noisy histogram Fig. 2(b) which is smeared out due to thermal noise with n th = 5 photons on average. The state was reconstructed using the noisy data and the knowledge that n th = 5. The high-fidelity result can be seen in Fig. 2(c). This is a reconstruction the neural network in Ref. [20] could not do.

Comparison to maximum likelihood estimation and neural network
We compare our convex optimization method (CVX) to a standard maximum likelihood estimation method (MLE) and a recent method based on a type of neural network called conditional generative adversarial network (CGAN) [19,20] with code provided in Ref. [67]. We observe that MLE and the CGAN can behave erratically depending on reconstruction parameters for the same underlying state, even without noise. For example, we show a test with different phase space limits α max = x max = p max and two different grid sizes in Fig. 3. The MLE is set to terminate when the Hilbert-Schmidt distance between two consecutive estimations is smaller than 10 −6 . The fidelity is fairly stable, but as seen in Fig. 3(b), it is generally the slowest method with the larger grid. The CGAN is initialized with random parameters for each training, leading to different results for each run. In Fig. 3(a) we show the result for two different random seeds. Not only does the fidelity vary for different seeds, it fluctuates significantly for different α max . Notably, our CVX method consistently reconstructs the state with fidelity 1, almost always with the fastest runtime. For more information on the time complexity of the CVX method, see Appendix D. We also test reconstruction using different Fock space

State reconstruction with experimental heterodyne data
Here we demonstrate the CVX state reconstruction method on experimental data produced for the publication Ref. [68]. The experiment consisted of selecting a particular mode from a propagating field which was emitted by a coherently driven qubit in front of a mirror. Particular modes were chosen by the shape of a temporal filter. The procedure for state reconstruction with noisy data used in this paper, and commonly used in circuit QED, was to extract moments of the bosonic operatorsâ,â † from the heterodyne histogram and then doing a maximum likelihood estimation of the state based on the extracted moments [69][70][71]. For states with large photon numbers, high-order moments are required, since nth order moments only contain information about number states up to n − 1. This means full information of a density matrix in a Fock space of dimension N = n − 1 can be obtained. This method is however limited to low photon-number states, since the standard deviations of the moments increase rapidly with higher and higher moments. For this reason, the Fock space was truncated to N = 4 in Ref. [68]. Notably, this is not a restriction for the CVX method. In Fig. 5(a) we show the fidelities between the moments-based MLE with N = 4 and CVX with N = 4 and N = 8, where the fidelity for the latter was calculated by truncating the reconstructed density matrix. The x-axis in Fig. 5 corresponds to the filter width, where larger widths are expected to correspond to states with higher photon numbers. Each filter width corresponds to separate states with separate measurement histograms. In Fig. 5(a) it can be seen that the fidelity between the MLE and CVX reconstructions with both N = 4 is close to 1 for all measurements, indicating that CVX performs at least as well as MLE. However, if the CVX reconstruction Fock space is increased to N = 8, fidelity starts to decrease for longer filter widths, i.e. states with higher photon numbers. In Fig. 5(b) we see that for the more severely truncated N = 4 density matrix, the three-photon population ρ 3 increases in lieu of higher photon-number populations. But with N = 8, we see that the four-photon probability ρ 4 (and to an even lesser degree higher photon-number states not shown in the plot) is slightly increasing. Because the CVX reconstruction with N = 8 matches the reconstruction with N = 4 with high fidelity for the shortest temporal filters, it is reasonable to suspect that the smaller density matrix from the moments-MLE method fails to completely correctly capture the states with a slightly larger photon number while the CVX method has no such issue.

IV. SUMMARY AND CONCLUSIONS
Convex optimization for quantum tomography is not a new idea in itself, but it does not seem to have been implemented for single-mode continuous variable state tomography previously.Unlike emerging machine learning methods, the convex optimization method is easy to understand and to use. Here we have demonstrated such a method on different types of realistic, noisy data. We showed that the convex optimization method gives stable results even when measurement parameters are varied, unlike an implementation of maximum likelihood estimation as well as a neural network, which both gave aberrant reconstructions for certain parameter combinations seemingly without reason. Besides the reliability of the convex optimization reconstruction, a major benefit is that unlike iterative maximum likelihood and neural network methods, there is no arbitrary stopping criteria for when to terminate the iterations, and the optimal solution is guaranteed to be reached. The method is easy to implement with openly available Python packages, and example code is provided in Ref. [59].
While continuous-variable states was the focus of this study, the presented framework and method is general and can also be used for discrete-variable systems. The state reconstruction takes less than 30 seconds for states occupying a Hilbert space of dimension N up to 32 [cf. Fig. 4(b)], and under ten minutes with N = 64 on a desk- for experimental data with different filter widths, corresponding to states with different photon numbers. (a) Fidelity between reconstructions. When CVX is truncated to the same Fock space dimension (N = 4) as the MLE, fidelity is close to 1, meaning the two reconstruction methods agree. When CVX is allowed a larger Fock space, fidelity decreases with longer temporal filters, corresponding to states with higher photon numbers. This indicates that N = 4 is not a large enough Fock space to reconstruct the correct state. Note that the y-axis begins at 0.9. (b) Three-photon populations ρ3 for the N = 4 and N = 8 reconstructions, and four-photon population ρ4 for the latter. It can be seen that for N = 4, the largest possible photon number state ρ3 increases excessively when higher photon number populations instead start appearing with the larger Fock space.
top computer with 16 cores. Based on this, the current method is expected to be practical for states with up to six qubits. It is notable that the method is not limited by the computational cost of the convex optimization itself, but rather by the construction of the necessary operators (see Appendix D). This was not the main concern of this work, so there is possibly room for improvement.

ACKNOWLEDGMENTS
Huge thanks to Shahnawaz Ahmed for valuable discussions about all things tomography as well as programming. Also big thanks to Yong Lu and Marina Kudra for providing experimental data. Thanks to Christopher Eichler for an informative discussion. Finally I thank Fernando Quijandría for discussions and proofreading the manuscript, and Göran Johansson for proofreading as well. This work was supported by Chalmers Excellence Initiative Nano.

Appendix A: Vectorization of ρ
The density matrix is cast to a vector in row-major order, meaning that the second column is placed under the first, and so on. As an example, we show the vectorization of a density matrix with N = 2: The vectorized version is Appendix B: State reconstruction with Wigner tomography The convex optimization state reconstruction can also be performed for Wigner tomography. The Wigner function can be measured as the expectation value of the displaced parity operator P [72,73]: with the displacement operator D(α) = e αâ † −α * â , wherê a andâ † are the annihilation and creation operators of the bosonic mode. From Eq. (B1) the corresponding measurement operators are seen to be Π k = D(α k )P D † (α k ). We test the method on measurement data of the target state (|0 + |4 )/ √ 2 as shown in Fig. 2(b) in Ref. [74], recreated in Fig. 6(a). The Wigner function and photon populations of the reconstructed state are shown in Figs. 6(b) and 6(c), respectively. The phase space grid is 61 × 61 with limit α max = 2.32, and the reconstruction was performed with Fock space dimension N = 10. Unit fidelity is not expected due to losses, and our reconstructed state has 0.95 fidelity to the target state, which is close to 0.94 as reported in the paper. The reconstruction in [74], was produced by the CGAN [19,20]- The expression for a heterodyne histogram corrupted by n thermal photons is given by inserting Eq. (19) into Eq. (18), which results in the Gaussian convolution We now introduce the s-parametrized quasiprobability distributionW (α, a) of which the Q-function, P -function and the Wigner function are special cases of for s = −1, 1, 0, respectively. For different values of s the distributions are related to each other via a Gaussian convolution [64] W (α, s) = 1 π To show the equivalency between thermal noise and detection inefficiency we utilize that there is a relation between the parameter s and the quantum efficiency η [75]: Inserting this and also setting t = −1 in Eq. (C2) gives Comparing Eqs. (C1) and (C4) we can after some simple algebra see that there is a relation or equivalently The equivalence between the models of loss and noise can intuitively be understood by imagining a fictitious beam splitter in front of a perfect detector. The transmissivity η of the beam splitter corresponds to the detection efficiency. Any and all types of losses can be modelled in this fashion; the overall detection efficiency η is simply the product of efficiencies of different parts of the system. Also, dissipation is intimately related to fluctuations [76], and this statistical noise is formally described by vacuum entering via the second port of the beam splitter. In the case of additional thermal noise, a thermal field can be imagined to enter the second port of the beam splitter [57].
Appendix D: Notes on the computational cost CVXPY is a modeling language for convex optimization problems [50,52], not a solver in itself. CVXPY supports several solvers, and the computational complexity may be different for different solvers. The results in this paper were obtained with the Splitting Conic Solver (SCS) [77,78]. The SCS algorithm comprises several steps, but the most time-consuming operation is projection on the semidefinite cone which has cubic complexity with respect to the matrix size [79]. The SCS method is constructed to handle very large problems and is indeed very efficient for the tomography problem. In fact, constructing the matrix A (4) is by far the most costly part of the CVX tomography program rather than the convex optimization itself. This is exemplified in Table II plication between the N × N matrices Π, Ω, which scales as O(N 3 ) [80], and taking the trace which scales linearly. As such, the overall time complexity of this step is O(N 3 ). This is to be done for each of the M N 2 elements of A, where M is the number of measurements. As such, the cost of constructing the entire matrix is O(N 4 M ).
The process can be parallelized to improve the speed, but even using 16 cores the matrix construction is the limiting factor of the CVX tomography program as can be seen in Table. II. It is possible this can be sped-up further, for example by using GPUs. It is difficult to ascertain the exact number of measurement samples that is needed to obtain sufficient statistics, since this is state dependent as can be seen in Fig. 7. Nevertheless, a comparison can be made between reconstruction methods. The figure shows reconstruction fidelities as a function of the number of samples for three different states: squeezed vacuum, a so-called binomial code state (|0 + |4 )/