Reconstructing high-dimensional two-photon entangled states via compressive sensing

Accurately establishing the state of large-scale quantum systems is an important tool in quantum information science; however, the large number of unknown parameters hinders the rapid characterisation of such states, and reconstruction procedures can become prohibitively time-consuming. Compressive sensing, a procedure for solving inverse problems by incorporating prior knowledge about the form of the solution, provides an attractive alternative to the problem of high-dimensional quantum state characterisation. Using a modified version of compressive sensing that incorporates the principles of singular value thresholding, we reconstruct the density matrix of a high-dimensional two-photon entangled system. The dimension of each photon is equal to d = 17, corresponding to a system of 83521 unknown real parameters. Accurate reconstruction is achieved with approximately 2500 measurements, only 3% of the total number of unknown parameters in the state. The algorithm we develop is fast, computationally inexpensive, and applicable to a wide range of quantum states, thus demonstrating compressive sensing as an effective technique for measuring the state of large-scale quantum systems.


Introduction
Many areas of quantum mechanics require the efficient and accurate measurement of entangled states.Perhaps the most traditional and widely adopted way of doing so is through full tomographic reconstruction [1], a technique that performs a series of independent measurements on the system in order to uniquely identify its nature.However, the complexity of such a method dramatically increases with increasing dimension of the system, and fully measuring the state of two entangled objects, each of d dimensions, requires at least d 4 measurements [2].As a result, full tomographic reconstruction is effective only at low dimensions and is otherwise prohibitively time consuming and computationally expensive.
Large-dimensional states are necessary for quantum computation and for certain quantum information protocols.Monz et al. reported the generation of a 14-qubit entangled state using trapped ions [3], and Yao et al. reported the generation of an 8-photon entangled state [4], although neither reported the density matrix for their respective states.Zhang et al. performed quantum tomography of a hybrid optical detector with over a million free parameters [5].However, to date, the largest density matrix reported for an entangled state is that of Häffner et al., who recorded the density matrix of 8 trapped ions [6].
Compressive sensing, which originates from the field of signal processing, provides a very efficient mechanism to establish properties of an unknown system with limited observations (see, e.g., Candès [7] and references therein).Compressive sensing uses prior assumptions in order to reduce the number of possible solutions, which can drastically reduce both measurement and processing time.Consequently, it is possible to establish descriptions of very large systems that could previously not be explored.This principle is used extensively in the fields of image reconstruction [8] and medical tomography [9], and it has recently been adopted in various areas of quantum science [10][11][12][13][14][15][16].
In this paper we propose and outline a compressive sensing technique that is able to successfully reconstruct the density matrix of near-pure entangled states of high dimensions.We implement this method to reconstruct a pure state of two 17-dimensional photons entangled in their orbital angular momentum.The recovery of the state is achieved by employing only 3% of the measurements that full tomographic reconstruction would require.The full procedure, including measuring and postprocessing, takes approximately three hours on a standard desktop computer.Our data processing algorithm is similar to the singular value thresholding algorithm detailed in [17]; however, its design is specifically adapted for near-pure entangled state reconstruction.The procedure is fast, computationally inexpensive, and robust to noise.

Theoretical description of compressive sensing and quantum tomography.
Compressive sensing is a data-processing technique widely used in different signal reconstruction applications.Its aim is to find the solution to underdetermined linear systems, under the assumption that such a solution is sparse in some basis.Such problems can be posed in the following way: where x ∈ C N ×1 represents a vector describing the measured object; A ∈ C M ×N is the matrix of measurements, with M N ; b ∈ C M ×1 is the vector of measurement results; • 1 denotes the 1 norm of the vector; and f is arXiv:1407.7426v2[quant-ph] 29 Jul 2014 a transformation to a space in which f (x) has a sparse representation.
In the specific case of quantum state tomography, the aim is to reconstruct an unknown near-pure density matrix, using an under-sampled set of measurements, under the assumption that such a matrix is low rank.The problem to be solved is then [10,18] Here, ρ is the density matrix to be reconstructed, while ρ ∈ C N ×1 is the density matrix in vector form; A ∈ C M ×N is the matrix of measurements; p ∈ C M ×1 is the vector of resulting probabilities; and • Tr stands for the trace norm of the matrix.The rows of the measurement matrix A are individual measurement vectors A i , and the elements of the vector p are the corresponding probabilities p i .

The algorithm.
We develop an operation-projection method similar to the singular value thresholding algorithm shown in [17] and implemented in [10]; however, we significantly modify its design in order to make use of the known features of near-pure entangled states.The algorithm requires an initial guess matrix to begin the procedure.The protocol then has two main stages: (i) the operations on the current matrix ρ to impose the desired characteristics and (ii) the projection of the resulting answer in vector form ρ onto the solution space.Applying these steps repeatedly constitutes an iterative procedure to approach the target solution.We interchange between the matrix form and vector form when implementing the operation and projection stages respectively.
In the operations stage, two steps are performed: First, the rank of the matrix is reduced by thresholding the eigenvalues below a certain level.This is achieved by decomposing the matrix into its eigenvalues and eigenvectors, setting the eigenvalues below the chosen threshold to zero, and then recomposing the matrix using where λ i is the i th eigenvalue and ϕ i the corresponding eigenvector.Second, to make use of the known sparsity characteristic and trace properties of the solution, we apply a thresholding operation on the individual matrix elements and normalise the result to have trace equal to unity.We achieve this by setting the elements that have modulus smaller than a chosen value to zero and dividing the matrix by its trace.
After the operation, the resultant matrix ρ0 has the desired characteristics of the solution; however, it no longer belongs to the linear space defined by the measurements A and probabilities p.To return the matrix ρ0 to the space defined by A ρ = p, we then implement the projection stage of the procedure.In order to describe the projection stage, we first introduce a geometrical formalism of the problem.
Each measurement vector A i and corresponding probability p i represents a hyperplane in a space of N dimensions, where N is the number of elements in ρ.The intersection of these hyperplanes represents the set of all solutions to the system A ρ = p.A simplified version of this concept is shown in Fig. 1, where two intersecting hyperplanes are represented as two-dimensional planes, and their intersection as a line.
After the operations stage, the matrix ρ0 is reshaped into vector form ρ 0 so that it can be projected onto the intersection of the hyperplanes defined by the linear system.The projection procedure is simple and computationally inexpensive if the hyperplanes are all perpendicular to each other.
However, although the matrix of random measurements A is nearly orthonormal, there is small non-zero overlap between any two measurements A i and A j (i = j).This is due to the physical limitations of the measurement procedure.For this reason, we transform the system A ρ = p into a new system A ρ = p , where A is an orthogonal matrix.This is achieved by multiplying both A and p by a matrix B such that BA = A .It is important to note that the system A ρ = p is a mathematical construct and no longer directly relates to the measurements and their corresponding probabilities; however, the set of solutions it defines is exactly the same as that of the original system.
In order to obtain a solution ρ s from the initial point ρ 0 , we progressively project ρ 0 on each hyperplane in turn.This procedure is initiated by projecting the initial point ρ 0 onto the first hyperplane, given by A 1 ρ = p 1 , to find a new point ρ 1 .This new point is then projected onto the second hyperplane, and we continue in this fashion until the desired solution ρ s is found.This occurs after M projections, where M is the number of measurements.Details of this projection procedure can be found in the Supplementary Materials.
Applying this operation-projection procedure repeatedly constitutes an iterative method that provides a solution exhibiting the desired characteristics and belongs to the linear system A ρ = p.The schematic outline of the algorithm is shown in Fig. 1, where the orange and red arrows represent the operation and projection steps, respectively.The method is considered complete when the distance between consecutive iterates is below a predetermined tolerance.

Noise correction.
In our system, noise manifests itself as errors in the measured probabilities.Such noise is unavoidable, and consequently, the density matrix that we recover ρ r will not correspond to the desired solution to the problem A ρ = p; instead, it will be a solution to the system A ρ = p + ∆ p, Space of solutions defined by the measurement A i and the resulting probability p i FIG.1: A schematic representation of the compressive sensing problem.The two green planes A i ρ = p i and A j ρ = p j each correspond to individual measurements and represent two different solution spaces.The intersection of the two planes corresponds to the set of all solutions belonging to the combined space A ρ = p, indicated by the red line.The curved line represents a set of potential solutions in the space that retain the desired characteristics of our answer.Our algorithm works by iterating between the set of solutions with the desired characteristics (orange line) and the set of solutions belonging to A ρ = p (red line).After a number of iterations, the algorithm converges to the solution of A ρ = p that possesses the desired characteristics.
where ∆ p is a vector of errors on the true probabilities.This error in probabilities results in a solution ρ r that is in fact some distance ∆ ρ from the desired solution ρ d in the space in which the algorithm operates.There are many methods for finding the solution in the presence of error [17,19].In our case, we determine ∆ ρ and subtract it from ρ r .This corresponds to the operation We use a priori knowledge of the desired state's characteristics to find ∆ ρ and systematically correct for noise in the system.Further details of our method can be found in the Supplementary Information.

Experimental implementation.
We have performed an experimental recovery of the density matrix of a 17-dimensional two-photon state in the orbital angular momentum (OAM) degree of freedom, produced by parametric downconversion (see Methods for details).The dimension of each photon is equal to d = 17, so the dimension of the entire state is 83521.The reconstruction is performed after 2506 projective measurements, which corresponds to only 3% of the total number of unknown parameters in the state.The reconstruction of the state is shown in Fig. 2.
The state that we measure exhibits strong anticorrelations in the OAM degree of freedom; that is to say that the OAM state | S in the signal photon is correlated with the state | − I in the idler photon.Additionally, the existence of the non-zero off-diagonal elements in the density matrix indicates a high degree of purity in the obtained state.These two features combine to suggest a high degree of entanglement of the OAM modes.
To characterise the entanglement, we use the fidelity of the reconstructed state ρ with the ideal, maximally entangled pure state The fidelity is then given by For the density matrix shown in Fig. 2, this fidelity was found to be 83.1%.In order to characterise the effectiveness of the reconstruction method, we reconstructed the density matrix of a 7-dimensional two-photon state with varying number of measurements.The resultant fidelities are shown in Fig. 3.We show the results both with and without our error correction procedure.For both cases, the fidelity increases as the number of measurements increases, indicating that more information produces a more accurate reconstruction.
However, for the case without error correction, the fidelity gradually decreases beyond 20% of the measurements.Because the measurements performed are nearly orthogonal to each other and are of insufficient number   The orange points correspond to the fidelities without error compensation; the green points correspond to the fidelities with our error compensation.The maximum value for the green points is 0.97 ± 0.01.
to yield a fully determined system, the errors contained within each measurement result do not average out to reduce the uncertainty, but instead sum to increase the uncertainty.Equivalently, every measurement taken into account restricts by one dimension the space of possible solutions to the underdetermined system: fidelity increases with increasing measurements at a low number of samples because the space is large enough to be very close to the desired solution, but the space gets smaller with increasing measurements, progressively excluding other low-rank sparse objects.At the high fidelity peak, the space is small enough so that the lowest rank and sparsest solution it contains is approximately the desired one and the algorithm will converge towards it.As the number of samples increases, the accumulation of errors results in a solution space that is far from the desired one; however, with additional samples, the dimension of the space is reduced.As a result, its distance from the sampled object increases and the algorithm yields an answer that diverges from the desired one.

Discussion
We have developed and tested an efficient method for determining the state of a quantum system based on a few simple assumptions.In this case, we use the prior knowledge of the sparsity of the density matrix associated with the system to achieve high-fidelity recovery from a small number of independent measurements of that system.Thus the state that we report corresponds to that which satisfies the set of measurements and the initial assumption of purity.One way to look at this is to say that we have answered the following question: "What is the purest state that is compatible with the set of measure-ments?"However, a feature of our method is that, using the same measurements and different prior knowledge, it can be readily refashioned to recover many states with a variety of desired characteristics.
In the case of a two-photon entangled state, where each photon exists in a 17-dimensional space with 83521 corresponding unknowns, we are able to recover the system with 3% of the measurements required for informational completeness of an unknown general quantum state.Our result corresponds to one of the largest discrete quantum states yet to be reported.We anticipate that the techniques implemented in this work will have impact in a wide range of areas in quantum science, including implementation and verification of quantum information protocols using high-dimensional states.

Methods
We use a 100-mW diode laser with wavelength 405 nm, along with a 3-mm-thick BBO crystal, to generate entangled photons through the process of parametric downconversion; see Fig. 4. The two-photon state that is generated in this process is given by where |c | 2 indicates the probability of finding the signal photon with OAM − and the idler photon with OAM .In our experiment, we limit the range of OAM states to values between = −8 and = 8.
The signal and idler photons are each incident on a separate half of a spatial light modulator (SLM), displaying computer-generated holograms, and then collected by a single-mode fibre connected to a single-photon detector.This results in a projective measurement on the twophoton mode.The result of the projection is measured by the coincidence detection with a coincidence window of 25 ns.
One of the keys to successful compressive sensing is to ensure that the measurement states are unstructured with respect to the basis in which the sampled state is sparse.For this experiment, that corresponds to measurement settings that are random superpositions of OAM modes.Therefore the measurement states |φ i S and |φ i I are generated from superpositions of OAM states where the coefficients a are generated at random We generate the matrix A, of Eq. ( 2), by performing a number of random separable projective measurements.Each row of A corresponds to the vector form of the individual projectors Âi .The projection operator Âi is given by Âi where the states |φ i S and |φ i I are the modes for the signal and idler arms respectively.The coincidence rate c i for each measurement Âi can be normalised to obtain the equivalent probability p i .Each probability p i constitutes the result of the corresponding measurement Âi .The probability of recording a coincidence count is given by Consequently, the linear system that is defined by the set of measurement operators { Âi } and the corresponding probabilities {p i } is where ρ is the vector form of the density matrix ρ.After performing an appropriate number of measurements, the task is then to solve the inverse problem under the constraints given by Eq. (2).

End
Here orth(•) is the orthogonalizing operator, Γ τ,τ (•) is the operator that enforces the desired characteristics described in the results section, vec(•) is the operator that rearranges the elements of a matrix into a vector, P (v, V ) denotes the projection of a vector v onto a hyperplane having normal V , and mat(•) is the operator that rearranges the elements of a vector into a square matrix.

Projection onto a hyperplane
To project a point ρ i−1 onto a hyperplane A i ρ = p i , it is necessary to find the vector v i that has direction n i normal to the hyperplane A i ρ = p i and size k i , where k i is The desired projection ρ i is then Finding and correcting ∆ ρ The error vector ∆ ρ depends on the experimental error ∆p i associated with each of the probabilities and the measurements made to perform the reconstruction.Instead of extending the search to a non-linear space, we make use of the low rank and sparsity information to estimate the error direction, that is, to find a close approximation for the vector ∆ ρ = ρ r − ρ d , where ρ r is the projection of ρ d onto the linear space intersection of all the hyperplanes.In order to do this we divide the measurements and the corresponding probabilities in subsets A si and p si of sufficient size for our compressive sensing technique to yield convergence to a solution (the size of the subsets depends on the purity of the state and the estimate of the error ∆p i on each of the probabilities).We then perform the operation-projection algorithm on each subset separately to find the vectors ρ si , low rank and sparse solutions to the corresponding subsets systems A si ρ = p si .We hence define a new set of vectors ∆ ρ i = ρ si − ρ di where ρ di is the solution resulting from applying one last set of thresholding operations to ρ si .The sum of the vectors ∆ ρ i is taken as a close approximation for the error vector ∆ ρ.We finally compute a correction vector for the probabilities ∆p that can be subtracted from the measured probabilities p.The correction vector is defined as ∆p = A∆ ρ.
Once the probabilities have been corrected as described above, the compressive sensing algorithm can be performed, making use of the set of performed measurements A and the corrected probabilities p − ∆p.

FIG. 2 :
FIG. 2: (a) The real part of the recovered density matrix.The dimension of each photon is equal to d = 17, so that the total dimension of the combined space is equal to 83521.The index i runs from i = −8 to 8. (b) The imaginary part of the recovered density matrix.

FIG. 3 :
FIG.3: Fidelity with the maximally entangled state vs. the percentage of measurements used for reconstruction for dimension d = 7 (100% corresponds to d 4 random projective measurements).The orange points correspond to the fidelities without error compensation; the green points correspond to the fidelities with our error compensation.The maximum value for the green points is 0.97 ± 0.01.