{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "In this example, we relax the phase retrieval problem similar to the classical [MaxCut](http://www-math.mit.edu/~goemans/PAPERS/maxcut-jacm.pdf) semidefinite program and recover the phase of the signal given the magnitude of the linear measurements.\n", "\n", "Phase recovery has wide applications such as in X-ray and crystallography imaging, diffraction imaging or microscopy and audio signal processing. In all these applications, the detectors cannot measure the phase of the incoming wave and only record its amplitude i.e complex measurements of a signal $x \\in \\mathbb{C}^p$ are obtained from a linear injective operator A, **but we can only measure the magnitude vector Ax, not the phase fo Ax**.\n", "\n", "Recovering the phase of Ax from |Ax| is a **nonconvex optimization problem**. Using results from [this paper](https://arxiv.org/abs/1206.0102), the problem can be relaxed to a (complex) semidefinite program (complex SDP).\n", "\n", "The original reprsentation of the problem is as follows:\n", "\n", ">>>> find x\n", "\n", ">>>> such that |Ax| = b\n", "\n", ">>>> where $x \\in \\mathbb{C}^p$, $A \\in \\mathbb{C}^{n \\times p}$ and $b \\in \\mathbb{R}^n$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this example, **the problem is to find the phase of Ax given the value |Ax|**. Given a linear operator $A$ and a vector $b= |Ax|$ of measured amplitudes, in the noiseless case, we can write Ax = diag(b)u where $u \\in \\mathbb{C}^n$ is a phase vector, satisfying |$\\mathbb{u}_i$| = 1 for i = 1,. . . , n. \n", "\n", "We relax this problem as Complex Semidefinite Programming.\n", "\n", "### Relaxed Problem similar to [MaxCut](http://www-math.mit.edu/~goemans/PAPERS/maxcut-jacm.pdf)\n", "\n", "Define the positive semidefinite hermitian matrix $M = \\text{diag}(b) (I - A A^*) \\text{diag}(b)$. The problem is:\n", "\n", " minimize < U,M >\n", " subject to \n", " diag(U) = 1\n", " U in :HermitianSemiDefinite\n", " \n", "Here the variable $U$ must be hermitian ($U \\in \\mathbb{H}_n $), and we have a solution to the phase recovery problem if $U = u u^*$ has rank one. Otherwise, the leading singular vector of $U$ can be used to approximate the solution." ] }, { "cell_type": "code", "execution_count": 101, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(size(coeff),size(var)) = ((40,40),(40,40))\n", "----------------------------------------------------------------------------\n", "\tSCS v1.1.8 - Splitting Conic Solver\n", "\t(c) Brendan O'Donoghue, Stanford University, 2012-2015\n", "----------------------------------------------------------------------------\n", "Lin-sys: sparse-direct, nnz in A = 3181\n", "eps = 1.00e-04, alpha = 1.80, max_iters = 20000, normalize = 1, scale = 5.00\n", "Variables n = 801, constraints m = 1641\n", "Cones:\tprimal zero / dual free vars: 821\n", "\tsd vars: 820, sd blks: 1\n", "Setup time: 1.03e-03s\n", "----------------------------------------------------------------------------\n", " Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time (s)\n", "----------------------------------------------------------------------------\n", " 0| inf inf -nan -inf -inf inf 1.17e-03 \n", " 100| inf inf -nan -inf -inf inf 1.48e-01 \n", " 200| inf inf -nan -inf -inf inf 2.76e-01 \n", " 300| inf inf -nan -inf -inf inf 4.01e-01 \n", " 400| 5.94e-02 1.47e-03 1.06e-05 -1.33e+03 -1.33e+03 2.40e-17 5.62e-01 \n", " 500| 1.22e-02 2.87e-04 2.14e-06 -1.33e+03 -1.33e+03 8.67e-18 6.90e-01 \n", " 600| 2.64e-03 6.29e-05 4.60e-07 -1.33e+03 -1.33e+03 7.28e-17 8.22e-01 \n", " 700| 5.73e-04 1.37e-05 9.98e-08 -1.33e+03 -1.33e+03 7.30e-17 9.63e-01 \n", " 800| 1.25e-04 2.98e-06 2.17e-08 -1.33e+03 -1.33e+03 7.73e-18 1.09e+00 \n", " 820| 9.18e-05 2.20e-06 1.60e-08 -1.33e+03 -1.33e+03 2.45e-17 1.11e+00 \n", "----------------------------------------------------------------------------\n", "Status: Solved\n", "Timing: Solve time: 1.11e+00s\n", "\tLin-sys: nnz in L factor: 6763, avg solve time: 5.92e-05s\n", "\tCones: avg projection time: 1.26e-03s\n", "----------------------------------------------------------------------------\n", "Error metrics:\n", "dist(s, K) = 6.5244e-09, dist(y, K*) = 1.8637e-09, s'y/m = -5.3273e-10\n", "|Ax + s - b|_2 / (1 + |b|_2) = 9.1773e-05\n", "|A'y + c|_2 / (1 + |c|_2) = 2.1952e-06\n", "|c'x + b'y| / (1 + |c'x| + |b'y|) = 1.6034e-08\n", "----------------------------------------------------------------------------\n", "c'x = -1326.1472, -b'y = -1326.1471\n", "============================================================================\n" ] }, { "data": { "text/plain": [ "20×20 Array{Complex{Float64},2}:\n", " 0.999998+1.49475e-15im … 0.641021-0.766111im \n", " 0.446918-0.760165im -0.279396-0.844012im \n", " 0.540963+0.840695im 0.991945+0.124332im \n", " 0.310127-0.919005im -0.498066-0.83538im \n", " 0.601155+0.79323im 0.997432+0.0489258im\n", " 0.318761+0.886608im … 0.899099+0.324567im \n", " -0.273286-0.917372im -0.866077-0.384836im \n", " 0.940279-0.213429im 0.435288-0.868812im \n", " -0.72136-0.548241im -0.863905+0.207653im \n", " 0.535909+0.814248im 0.975065+0.118254im \n", " 0.766546+0.641349im … 0.981572-0.175106im \n", " 0.858843-0.134508im 0.44869-0.767066im \n", " 0.491533-0.848226im -0.335281-0.929426im \n", " 0.939298+0.337107im 0.862322-0.505741im \n", " -0.289044-0.902153im -0.862436-0.361781im \n", " 0.872244+0.474994im … 0.923208-0.369149im \n", " 0.532716+0.829509im 0.984631+0.122242im \n", " -0.446504-0.715904im -0.81144-0.10799im \n", " 0.876944+0.475247im 0.927712-0.364238im \n", " 0.641021+0.766111im 0.999993+3.07926e-15im " ] }, "execution_count": 101, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using Convex\n", "n = 20\n", "p = 2\n", "A = rand(n,p) + im*randn(n,p)\n", "x = rand(p) + im*randn(p)\n", "b = abs(A*x) + rand(n)\n", "\n", "M = diagm(b)*(eye(n)-A*ctranspose(A))*diagm(b)\n", "U = ComplexVariable(n,n)\n", "objective = inner_product(U,M)\n", "c1 = diag(U) == 1 \n", "c2 = U in :SDP\n", "p = minimize(objective,c1,c2)\n", "solve!(p)\n", "U.value" ] }, { "cell_type": "code", "execution_count": 107, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2\n" ] }, { "data": { "text/plain": [ "20-element Array{Complex{Float64},1}:\n", " 0.642427-0.766347im \n", " -0.320032-0.947407im \n", " 0.992174+0.124859im \n", " -0.513966-0.857811im \n", " 0.998852+0.0479087im\n", " 0.942305+0.334756im \n", " -0.912508-0.40906im \n", " 0.446383-0.894842im \n", " -0.974714+0.223458im \n", " 0.993151+0.116834im \n", " 0.9847-0.174256im \n", " 0.498451-0.866918im \n", " -0.342385-0.93956im \n", " 0.86251-0.506039im \n", " -0.920369-0.391052im \n", " 0.92876-0.370682im \n", " 0.992626+0.121218im \n", " -0.98923-0.146366im \n", " 0.93096-0.365121im \n", " 1.0+0.0im " ] }, "execution_count": 107, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Verify if the rank of U is 1\n", "B, C = eig(U.value);\n", "println(length([e for e in B if(abs(real(e))>1e-4)]))\n", "#Decompose U = uu*\n", "# u is the phase of Ax\n", "u = C[:,1];\n", "for i in 1:n\n", " u[i] = u[i]/abs(u[i])\n", "end\n", "u" ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Julia 0.5.0-rc3", "language": "julia", "name": "julia-0.5" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.5.1" } }, "nbformat": 4, "nbformat_minor": 0 }