<<MaTeX`

SetOptions[MaTeX,"Preamble"->{"\\usepackage{color,txfonts}"}];

A generalization of RRR

As before, we decompose the dynamical map in the neighborhood of a solution into mutually orthogonal spaces (Subscript[X, ‖])≡span((OverTilde[A]),(OverTilde[B])) and (Subscript[X, ⊥])≡span((Subscript[c, A]),(Subscript[c, B])), where x∈(Subscript[X, ‖])⊕(Subscript[X, ⊥])=(ℝ^N). We thus write x=((Subscript[x, ‖]),(Subscript[x, ⊥])) where (Subscript[x, ‖,⊥])∈(Subscript[X, ‖,⊥]) and (Subscript[P, A,B])(x)=((Subscript[P, (OverTilde[A]),(OverTilde[B])])((Subscript[x, ‖])),(Subscript[c, A,B])). Then the equation of motion takes the form

We see that (Subscript[(OverDot[x]), ⊥])=(Subscript[c, A])-(Subscript[c, B]). Thus, when (Subscript[c, A])≠(Subscript[c, B]), (Subscript[(OverDot[x]), ⊥])≠0 and the system flows away monotonically from the non-solution. In the feasible case ((Subscript[c, A])=(Subscript[c, B])), (Subscript[(OverDot[x]), ⊥])=0. Then the problem reduces to the dynamics in (Subscript[X, ‖]), given by

To see whether the system approaches the desired solution, we look at the how the norm of (Subscript[x, ‖])evolves:

This equation can be simplified by decomposing (Subscript[x, ‖])into parts belonging to (OverTilde[B]) and orthogonal to (OverTilde[B]). Thus, we write (Subscript[x, ‖])=(Subscript[x, (OverTilde[B])])+(Subscript[x, (OverBar[B])]) where (Subscript[x, (OverBar[B])]^T)(Subscript[x, (OverTilde[B])])=0 and (Subscript[X, ‖])=(OverTilde[B])⊕(OverBar[B]). Substituting this expression into the above equation and using the relations (Subscript[P, (OverTilde[B])])(Subscript[x, (OverTilde[B])])=(Subscript[x, (OverTilde[B])]), (Subscript[P, (OverTilde[B])])(Subscript[x, (OverBar[B])])=0, and (P^T)=P, we find

where (OverBar[A]) denotes the subspace of (Subscript[X, ‖]) that is orthogonal to (OverTilde[A]). This is same rate we had obtained in class apart from the factor of γ. Using the property (P^2)=P, the rate can be expressed as a sum of squares,

Thus, the norm decreases monotonically until both (Subscript[P, (OverBar[A])])(Subscript[x, (OverTilde[B])])=0 and (Subscript[P, (OverTilde[A])])(Subscript[x, (OverBar[B])])=0, which is true if and only if (Subscript[P, (OverTilde[A])])(Subscript[x, ‖])=(Subscript[P, (OverTilde[B])])(Subscript[x, ‖]). Therefore the system approaches a hyperplane of fixed points corresponding to the desired solution (QED).

MaTeX codes for the displayed equations:

MaTeX["\\color{blue} \\partial_t \\lVert{x_{\\parallel}}\\rVert^2 = -2\\gamma \\left(\\lVert P_{\\bar{A}} x_{\\tilde{B}} \\rVert^2 + \\lVert P_{\\tilde{A}} x_{\\bar{B}} \\rVert^2 \\right) \\leq 0\\;.",Magnification->2]

MaTeX["
 \\color{blue}\\begin{aligned}
 \\frac{1}{2}\\partial_t \\lVert{x_{\\parallel}}\\rVert^2 &= \\left(x_{\\tilde{B}}^{T} + x_{\\bar{B}}^{T}\\right) \\left[ (1+\\gamma) P_{\\tilde{A}} x_{\\tilde{B}} - (1-\\gamma) P_{\\tilde{B}} P_{\\tilde{A}} (x_{\\tilde{B}} + x_{\\bar{B}}) - \\gamma P_{\\tilde{A}} (x_{\\tilde{B}} + x_{\\bar{B}}) - \\gamma x_{\\tilde{B}} \\right]\\\\
 &= \\left(x_{\\tilde{B}}^{T} + x_{\\bar{B}}^{T}\\right) \\left[- \\gamma x_{\\tilde{B}} + P_{\\tilde{A}} x_{\\tilde{B}} - \\gamma P_{\\tilde{A}} x_{\\bar{B}} - (1-\\gamma) P_{\\tilde{B}} P_{\\tilde{A}} (x_{\\tilde{B}} + x_{\\bar{B}}) \\right] \\\\
 &= -\\gamma x_{\\tilde{B}}^{T} x_{\\tilde{B}} + \\left(x_{\\tilde{B}}^{T} + x_{\\bar{B}}^{T}\\right)  \\left(P_{\\tilde{A}} x_{\\tilde{B}} - \\gamma P_{\\tilde{A}} x_{\\bar{B}}\\right) - (1-\\gamma) x_{\\tilde{B}}^{T} P_{\\tilde{A}} (x_{\\tilde{B}} + x_{\\bar{B}}) \\\\
 &= -\\gamma x_{\\tilde{B}}^{T} \\left(1 - P_{\\tilde{A}}\\right) x_{\\tilde{B}} - \\gamma x_{\\bar{B}}^{T} P_{\\tilde{A}} x_{\\bar{B}} + x_{\\bar{B}}^{T} P_{\\tilde{A}} x_{\\tilde{B}} - x_{\\tilde{B}}^{T} P_{\\tilde{A}} x_{\\bar{B}} \\\\
 &= -\\gamma \\left( x_{\\tilde{B}}^{T} P_{\\bar{A}} x_{\\tilde{B}} + x_{\\bar{B}}^{T} P_{\\tilde{A}} x_{\\bar{B}} \\right)
 \\end{aligned}
",Magnification->2]

MaTeX["\\color{blue}\\frac{1}{2}\\partial_t \\lVert{x_{\\parallel}}\\rVert^2 =x_{\\parallel}^{T} \\dot{x}_{\\parallel} = x_{\\parallel}^{T} \\left[(1+\\gamma) P_{\\tilde{A}} P_{\\tilde{B}} - (1-\\gamma) P_{\\tilde{B}} P_{\\tilde{A}} -\\gamma (P_{\\tilde{A}} + P_{\\tilde{B}})\\right]x_{\\parallel}\\;.",Magnification->2]

MaTeX["\\color{blue}\\dot{x}_{\\parallel}=\\left[(1+\\gamma) P_{\\tilde{A}} P_{\\tilde{B}} - (1-\\gamma) P_{\\tilde{B}} P_{\\tilde{A}} -\\gamma (P_{\\tilde{A}} + P_{\\tilde{B}})\\right](x_{\\parallel})\\;.",Magnification->2]

MaTeX["\\color{blue}(\\dot{x}_{\\parallel},\\dot{x}_{\\perp})=\\left[P_{\\tilde{A}}\\left((1+\\gamma) P_{\\tilde{B}}(x_{\\parallel}) - \\gamma x_{\\parallel}\\right) - P_{\\tilde{B}}\\left((1-\\gamma) P_{\\tilde{A}}(x_{\\parallel}) + \\gamma x_{\\parallel}\\right), c_A - c_B \\right].",Magnification->2]

He who must not be retrieved

But retrieve we will! Let’s first import the diffraction data:

fourierMagData1=Import["https://raw.githubusercontent.com/veitelser/RRR_mini_course/master/data/DarkLordData","Table"];

gridsize1=Length[fourierMagData1]

Here’s what it looks like:

ListDensityPlot[RotateLeft[fourierMagData1,{1/2,1/2}*gridsize1],FrameTicksStyle->Directive[Black,12],PlotLegends->Automatic]

Next we will determine the support for the Dark Lord from the autocorrelation:

MatrixPlot[autocorrelation1=RotateLeft[Chop@InverseFourier[(fourierMagData1^2)],{1/2,1/2}*gridsize1],FrameTicksStyle->Directive[Black,12],Epilog->{EdgeForm[Black],Opacity[0],Rectangle[{39,15},{218,240}]}]

The autocorrelation is contained well within a rectangle of size (240-15+1) × (218-39+1)=226×180. Thus, we choose a rectangular support of size 113 × 90:

L1={113,90};

Now let’s define the two projections. First, the Fourier magnitude projection (note that if x is real, the projection will also be real):

fourierProj1[x_]:=Chop[InverseFourier[fourierMagData1*Exp[I*Arg[Fourier[x]]]]]

Second, the fixed support constraint. This is implemented by first dropping all entries outside the support, then clipping all negative elements, then padding:

supportProj1[x_]:=PadRight[Clip[Drop[x,First[L1]-gridsize1,Last[L1]-gridsize1],{0,Infinity}],{1,1}*gridsize1]

We are now ready to use the RRR algorithm. We’ll iterate with β=0.5, starting from a random array:

β=0.5;

norm[x_]:=Sqrt[Mean[((Flatten[x])^2)]]

x1=PadRight[RandomReal[{0,1},L1],{1,1}*gridsize1];
error1={};

Do[
fourierProjection1=fourierProj1[x1];
supportProjection1=supportProj1[2*fourierProjection1-x1];
diff1=supportProjection1-fourierProjection1;
x1+=β*diff1;
AppendTo[error1,norm[diff1]];
If[Last[error1]<0.0001,Break[]];
,{2000}]//AbsoluteTiming

The error bounces around before finally becoming zero:

ListPlot[error1,Joined->True,PlotRange->{0,0.1},Frame->True,FrameTicksStyle->Directive[Black,12],FrameLabel->{Style["Iteration",Black,16],MaTeX["\\text{Error $\\equiv$ }\\lVert \\dot{x} \\rVert / \\sqrt{N}",Magnification->2]},ImageSize->450]

Finally, here is our Dark Lord:

Row[{ArrayPlot[Take[fourierProjection1,First[L1],Last[L1]],ImageSize->300,PlotLabel->Style["Fourier projection",Black,16]],ArrayPlot[Take[supportProjection1,First[L1],Last[L1]],ImageSize->300,PlotLabel->Style["Support projection",Black,16]],ArrayPlot[Take[x1,First[L1],Last[L1]],ImageSize->300,PlotLabel->Style["x",Black,16]]},"         "]

Crystallographic phase retrieval

Sparsity projections

Suppose we have a real-valued signal ρ on N pixels. The projection  should map ρ to the nearest signal ρ' that is nonzero on at most |S| of the pixels. This is most easily implemented by setting the smallest N - |S| entries (by magnitude) in ρ to zero. Let’s do an example for N=64 and |S|=10:

nonzeros=10;

ProjS[ρ_]:=Block[{y=Flatten[ρ]},y[[Ordering[Abs[y],Length[y]-nonzeros]]]=0.;Partition[y,Length[ρ]]]

MatrixForm[ρ=RandomReal[{-1,1},{8,8}]]

MatrixForm[ProjS[ρ]]

If the projected signal ρ' is also required to be non-negative, one has to first clip all negative elements in ρ (to zero), then set the smallest N - |S| (non-negative) entries to zero. One could also check whether they are already zero after clipping, but it is not necessary. Here is an example for N=64 and |S|=10:

ProjSplus[ρ_]:=Block[{y=Flatten[Clip[ρ,{0,Infinity}]]},y[[Ordering[y,Length[y]-nonzeros]]]=0.;Partition[y,Length[ρ]]]

MatrixForm[ρ=RandomReal[{-1,1},{8,8}]]

MatrixForm[ProjSplus[ρ]]

Note that if the original signal ρ does not have enough positive entries, the projected signal will have less than |S| (=10) nonzero elements:

MatrixForm[ρ=RandomReal[{-15,1},{8,8}]]

MatrixForm[ProjSplus[ρ]]

Reconstructing atom distribution from diffraction data:

First, import the data for photon counts:

photoncounts=Import["https://raw.githubusercontent.com/veitelser/phase-retrieval-benchmarks/master/data/data100E"];

gridsize2=Length[photoncounts]

Next we construct the Fourier magnitudes throughout the N×N grid (here N=128). To do this, we note that the Fourier magnitudes are proportional to the square root of the photon counts and use the inversion symmetry and periodicity of the Fourier magnitudes. Thus, e.g., |(OverHat[ρ])(π/N+(Subscript[q, x]),π/N+(Subscript[q, y]))| = |(OverHat[ρ])(π/N-(Subscript[q, x]),π/N-(Subscript[q, y]))| .

fourierMagData2=
Block[{lefthalf,righthalf},
lefthalf=N@Sqrt[photoncounts];
righthalf=ArrayPad[RotateRight[Partition[Reverse[Flatten[lefthalf[[All,2;;]]]],gridsize2/2-1]],{{0},{1,0}}];
ArrayFlatten[{{lefthalf,righthalf}}]
];

Here’s how it looks:

ListDensityPlot[RotateLeft[fourierMagData2,{1/2,1/2}*gridsize2],FrameTicksStyle->Directive[Black,12],PlotLegends->Automatic]

Although we don’t need it here, here’s what the autocorrelation of the electron density looks like (within the unit cell):

MatrixPlot[autocorrelation2=RotateLeft[Chop@InverseFourier[(fourierMagData2^2)],{1/2,1/2}*gridsize2],FrameTicksStyle->Directive[Black,12]]

We’ll now define the projections. The Fourier magnitude projection needs a slight modification. Since forward-scattered photons are not detected, there is no measurement of the zeroth Fourier component (OverHat[ρ])(0,0), which is simply the average electron density. This is a free parameter, and thus, we simply leave this Fourier coefficient unchanged during projection.

fourierProj2[x_]:=Block[{fourier=Fourier[x]},Chop[InverseFourier[ReplacePart[fourierMagData2*Exp[I*Arg[fourier]],{1,1}->fourier[[1,1]]]]]]

Next, we use the projection  from above to implement the sparsity constraint. We know that there are about 100 atoms which should take up 8 pixels each on average. Thus, we impose the constraint that the signal can have at most 800 nonzero entries which are all non-negative.

zeros=(gridsize2^2)-800

supportProj2[x_]:=Block[{y=Flatten[Clip[x,{0,Infinity}]]},y[[Ordering[y,zeros]]]=0.;Partition[y,gridsize2]]

Now it’s time to have the RRR machinery run its course:

x2=RandomReal[{0,1},{1,1}*gridsize2];
error2={};

Do[
fourierProjection2=fourierProj2[x2];
supportProjection2=supportProj2[2*fourierProjection2-x2];
diff2=supportProjection2-fourierProjection2;
x2+=β*diff2;
AppendTo[error2,norm[diff2]];
If[Last[error2]<0.0001,Break[]];
,{1000}]//AbsoluteTiming

The error decreases at first, but never vanishes due to the noise introduced by the finite photon counts. It keeps fluctuating around a stable value:

ListPlot[error2,Joined->True,PlotRange->{0,6.5},Frame->True,FrameTicksStyle->Directive[Black,13],FrameLabel->{Style["Iteration",Black,16],MaTeX["\\text{Error $\\equiv$ }\\lVert \\dot{x} \\rVert / \\sqrt{N}",Magnification->2]},ImageSize->450]

Finally, here’s what we find for the reconstructed atom distribution within the unit cell:

Row[{ArrayPlot[fourierProjection2,ImageSize->300,PlotLabel->Style["Fourier projection",Black,16]],ArrayPlot[supportProjection2,ImageSize->300,PlotLabel->Style["Support projection",Black,16]],ArrayPlot[x2,ImageSize->300,PlotLabel->Style["x",Black,16]]},"         "]

---------------The End---------------