A convex optimization approach for image recovery from nonlinear measurements in optical interferometry Anna Aur´ ıa * , Rafael E. Carrillo * , and Yves Wiaux *†‡ * Institute of Electrical Engineering, Ecole Polytechnique F´ ed´ erale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland. † Department of Radiology and Medical Informatics, University of Geneva (UniGE), CH-1211 Geneva, Switzerland ‡ Department of Radiology, Lausanne University Hospital (CHUV), CH-1011 Lausanne, Switzerland. Abstract—Image recovery in optical interferometry is an ill-posed nonlinear inverse problem arising from incomplete power spectrum and bi-spectrum measurements. We formulate a linear version of the problem for the order-3 tensor formed by the tensor product of the signal with itself. This linear problem is regularized by standard convex ℓ 1 -relaxations of sparsity and low rank constraints and solved using the most advanced algorithms in convex optimization. We show preliminary results on small size synthetic images as a proof of concept. Interferometry is a unique tool to image the sky at otherwise inaccessible resolutions. The set of visibilities measured provides an incomplete Fourier coverage of the brightness intensity as a function of angular direction, in vector form x ∈ R N with components xi . However, at optical wavelengths, the phase of the complex visibilities is affected by atmospheric turbulence. The measurable quantities are power spectrum data | ˆ xi | 2 , and phases associated with the bi- spectrum ˆ xi · ˆ xj · ˆ x k [1]. This poses a nonlinear inverse problem for image reconstruction, which is sensitive to the optimization strategy. Generalizing the Phase Lift approach [2], we formulate a linear ver- sion of the problem for the order-3 tensor X = x⊗x⊗x ∈ R N×N×N with components X ijk , which arises from the tensor product of the signal with itself. The total flux is measured independently and we consider a normalized signal such that ∑ i xi =ˆ x0 =1, ˆ x0 representing the zero-frequency. Thus, the linear measurement model y = A(X ) ∈ C M encompasses both power spectrum and bi- spectrum measurements, for a measurement operator A defined as a selection operator after Fourier transform along all tensor dimensions. In this setting, prior information is essential to regularize the ill- posed inverse problem in the perspective of image reconstruction. Firstly, we adopt a sparsity prior on the tensor to acknowledge some signal sparsity K ≪ N . While ℓ0−minimization would promote sparsity explicitly we adopt the common ℓ1−convex relaxation. Secondly, convex reality and positivity constraints are also enforced to acknowledge the fact that we deal with intensity images. Thirdly, in order to counter-balance the large increase of dimensionality when solving for X instead of x, we rely explicitly on the fact that X is formed by the tensor product of x, whose components sum to unity, so that summations over one or two indices respectively lead to the order-2 tensor C(X )= x ⊗ x ∈ R N×N , and to the signal x itself. We enforce the semi-definite positivity of C(X ) and its rank- 1 structure by resorting to nuclear norm minimization. The nuclear norm of an order-2 tensor is defined as the ℓ1−norm of its singular values vector. Its minimization should be understood as the convex relaxation of the minimization of the rank function counting the number of singular values. Finally, we also resort to a re-weighting scheme consisting in approaching both ℓ0−minimization on X and rank minimization on C(X ) by solving a sequence of weighted ℓ1 and nuclear norm minimizations, each initialized to the solution of the previous problem. The fundamental symmetry of the tensor X over index permutation is also enforced by ensuring that any operation performed preserves the symmetry. The weighted ℓ1 and nuclear norm minimization problem thus reads as: min X ||C(X )|| W * + λ||X || W 1 such that ||y −A(X )||2 ≤ ǫ, and C(X )  0, X≥ 0, (1) where the symbols || · || W * and || · || W 1 respectively denote weighted nuclear and ℓ1 norms. In the weighted nuclear norm, the singular values of C(X ) are essentially divided by their value at the previous iteration, in order to approximate the rank function. In the weighted ℓ1−norm, each tensor component X ijk is essentially divided by some robust estimation of its value from the previous iteration, in order to promote ℓ0−minimization. This estimation is obtained by symmetrized sums over two dimensions in order to promote structure in the tensor sparsity. In the first iteration, no weighting is applied. A non-weighted ℓ1−norm is not a meaningful prior function as the tensor values sum to unity. We thus set λ =0 at the first iteration. We solve this complex problem taking advantage of the versatility of convex optimization, using a combination of the Douglas-Rachford and dual forward-backward algorithms. The solution is low rank and we extract x as the principal eigenvector of C(X ). We want to highlight the fact that in this framework the results do not depend on the initialization of the algorithm, in stark contrast with common non-convex approaches. Figure 1 shows an example of reconstruction of a 16 ×16 synthetic image from M =0.75N measurements affected by 30dB of input noise, along with a phase transition diagram for random 8 ×8 images, representing the probability of good reconstruction in the sparsity- undersampling plane in a noiseless setting. In both cases equal numbers of random power and bi-spectrum data are considered. Figure 1. Left and centre: 16 × 16 image and reconstruction with M = 0.75N : SNR= 37.2dB. Right: phase transition diagram for 8 × 8 images. The principal drawback of this approach is the dimension of the problem, leading to computation time and memory requirements issues. Advanced algorithmic, coding and hardware solutions need to be investigated in this perspective. REFERENCES [1] E. Thi´ ebaut and J. Giovanelli, “Image reconstruction in optical interfer- ometry,” IEEE Signal Processing Magazine, vol. 21, pp. 97–109, 2010. [2] E. Cand` es, T. Strohmer, and V. Voroninski, “Phaselift: Exact and stable signal recovery from magnitude measurements via convex programming,” Communications on Pure and Applied Mathematics, 2011, preprint in arXiv:1109.4499v1.