Quantum state tomography

We first reconstruct the single-qubit spin operators in the standard basis from the last assignment:

Do[
σ[i,j]=KroneckerProduct@@RotateRight[{PauliMatrix[j],IdentityMatrix[2],IdentityMatrix[2],IdentityMatrix[2]},i-1],
{i,4},{j,3}]

Next, we form a set of linear equations which describe the polarization measurements. For a density-matrix X, the measurement μ of an operator H is given by the equation Tr(X H)=μ,  or  (Subscript[Σ, i,j])(Subscript[H, j,i])(Subscript[X, i,j])=μ,  where the sum is over all basis states. This equation describes a linear operation on a vector  (Overscript[X, →])  composed of the flattened matrix elements (Subscript[X, i,j]). Thus, we  can represent all such measurements by a single matrix equation  A.(Overscript[X, →])=(Overscript[μ, →]),  where the rows of  A  are composed of the flattened matrix elements (Subscript[H, j,i]), as below:

A=
SparseArray[
Join[
{Flatten[IdentityMatrix[Length[σ[1,1]]]]},
Flatten[Table[Flatten[Transpose[σ[α,j]]],{α,4},{j,3}],1],
Flatten[Table[Flatten[Transpose[σ[1,i].σ[α,j]]],{i,3},{α,2,4},{j,3}],2]
]]

We form the vector  (Overscript[μ, →])  from the given measurements:

μ=Flatten[{1,{-1/4,0,0},{-1/4,0,0},{0,0,0},{1/4,1/4,0},{0,1/4,-1/4},{-1/4,1/4,0},{0,0,0},{0,1/4,0},{-1/4,1/4,3/4},{0,0,1/4},{0,-1/4,0},{-1/4,-1/4,0},{0,0,0}}];

In the RRR algorithm, we will be searching for the actual density matrix in the space of all 16×16 Hermitian matrices (note, there are 16 basis states). Given such a matrix X, the projection to the nearest (flattened) Hermitian matrix which satisfies the measurement constraint  A.(Overscript[X, →])=(Overscript[μ, →])  is given by  (Subscript[P, M])(X)=( I-(((SuperDagger[A])(A (SuperDagger[A])))^(-1))A).(Overscript[X, →])+(((SuperDagger[A])(A (SuperDagger[A])))^(-1)).(Overscript[μ, →]) :

scale=SparseArray[IdentityMatrix[Last[Dimensions[A]]]-ConjugateTranspose[A].Inverse[A.ConjugateTranspose[A]].A]

shift=SparseArray[ConjugateTranspose[A].Inverse[A.ConjugateTranspose[A]].μ]

measProj[x_]:=Partition[scale.Flatten[x]+shift,Length[x]]

For a (non-null) pure state ψ, the density matrix is given by (Subscript[X, i,j])=(Subscript[ψ, i])(Subscript[ψ, j]^( * )). Thus, it is a positive-definite Hermitian matrix of rank 1. The projection to this constraint is implemented by forming the eigen-decomposition of  X  and retaining only the term with largest positive eigenvalue:

pureStateProj[x_]:=
Block[{val,vec},
{val,vec}=Eigensystem[x,1,Method->{"Arnoldi","Criteria"->"RealPart"}];
If[Chop@First[val]>0,Transpose[vec].DiagonalMatrix[val].Conjugate[vec],ConstantArray[0,Dimensions[x]]]
]

We are now ready to run the RRR algorithm:

norm[list_]:=Sqrt[Mean[((Abs[Flatten[list]])^2)]];

β=0.5;

x=Block[{y=RandomComplex[{-1-I,1+I},{16,16}]},y+ConjugateTranspose[y]];
error={};

Do[
proj1=pureStateProj[x];
proj2=measProj[2*proj1-x];
diff=proj2-proj1;
x+=β*diff;
AppendTo[error,norm[diff]];
If[Last[error]<(10^(-12)),Break[]];
,{10000}]//AbsoluteTiming

The error decreases exponentially with iterations:

ListLogPlot[error,Joined->True,PlotRange->All,Frame->True,FrameTicksStyle->Directive[Black,12],FrameLabel->{Style["Iteration",Black,16],Style["Error",Black,16]},ImageSize->400]

Chop@Eigenvalues[proj1]

The system is in an equal superposition of all the basis states!

Normalize@Rationalize[Chop@First[proj1],(10^(-6))]

Complex Hadamard matrices

We will be searching for 7×7 complex Hadamard matrices, X, and see if we find any that are different from the unitary Fourier matrix F. Let us first define F:

n=7;

fourierMat=Table[Exp[2*π*I*p*q/n],{p,0,n-1},{q,0,n-1}];

There are two constraints to satisfy: (i) X must be unitary and (ii) all the elements of X must have the same magnitude. The projection to the nearest unitary matrix, (Subscript[P, U])(X), is obtained by forming the singular value decomposition, X=V D W, and setting the diagonal matrix D to identity, i.e., (Subscript[P, U])(X)=V W.

unitaryProj[x_]:=Block[{svd=SingularValueDecomposition[x]},First[svd].ConjugateTranspose[Last[svd]]]

The equal-magnitude constraint is implemented simply by normalizing each element. In addition, both constraints are unaffected when we multiply rows or columns of X by arbitrary phases. Thus, without loss of generality, we adopt the “Fourier convention,” which amounts choosing the phases such that all the elements in the first row and first column are equal to 1/(Sqrt[n]).

magnitudeProj[x_]:=Block[{norm=1/Sqrt[Length[x]]},ArrayPad[norm*Exp[I*Arg@x[[2;;,2;;]]],{1,0},N@norm]]

Below we present four different runs with the RRR algorithm. We find Hadamard matrices which are not found by permuting rows/columns of the Fourier matrix.

norm[list_]:=Sqrt[Mean[((Abs[Flatten[list]])^2)]];

β=0.5;

RRR - 1st run:

x=RandomComplex[{-1-I,1+I},{n,n}];
error={};

Do[
proj1=unitaryProj[x];
proj2=magnitudeProj[2*proj1-x];
diff=proj2-proj1;
x+=β*diff;
AppendTo[error,norm[diff]];
If[Last[error]<(10^(-10)),Break[]];
,{50000}]//AbsoluteTiming

The error decreases exponentially but has sustained oscillations:

ListLogPlot[error,Joined->True,PlotRange->All,Frame->True,FrameTicksStyle->Directive[Black,12],FrameLabel->{Style["Iteration",Black,16],Style["Error",Black,16]},ImageSize->400]

Test the constraints:

MatrixForm/@{Chop[proj1.ConjugateTranspose[proj1]],Chop[ConjugateTranspose[proj1].proj1],Sqrt[n]*Abs[proj1]}

The Hadamard matrix we found:

Sqrt[n]*(hadamard1=Chop@proj1)//MatrixForm

One way to check if the above matrix is identical to some row/column permutation of the unitary Fourier matrix is to plot the distribution of the phases:

Row[{Histogram[Flatten[Mod[Arg[hadamard1]/2*π,1]],{0.01},Frame->True,FrameTicksStyle->Directive[Black,14],FrameLabel->{Style["Normalized phases",Black,16],Style["Multiplicity",Black,16]},ImageSize->400,PlotLabel->Style["What we found",18]],
Histogram[Flatten[Mod[Arg[fourierMat]/2*π,1]],{0.01},Frame->True,FrameTicksStyle->Directive[Black,14],FrameLabel->{Style["Normalized phases",Black,16],Style["Multiplicity",Black,16]},ImageSize->400,PlotLabel->Style["Fourier",18]]},"        "]

Since the two distributions are different, the Hadamard matrix we found is not equivalent to the Fourier matrix.

RRR - 2nd run:

x=RandomComplex[{-1-I,1+I},{n,n}];
error={};

Do[
proj1=unitaryProj[x];
proj2=magnitudeProj[2*proj1-x];
diff=proj2-proj1;
x+=β*diff;
AppendTo[error,norm[diff]];
If[Last[error]<(10^(-10)),Break[]];
,{50000}]//AbsoluteTiming

The error decreases exponentially:

ListLogPlot[error,Joined->True,PlotRange->All,Frame->True,FrameTicksStyle->Directive[Black,12],FrameLabel->{Style["Iteration",Black,16],Style["Error",Black,16]},ImageSize->400]

The Hadamard matrix we found:

Sqrt[n]*(hadamard2=Chop@proj1)//MatrixForm

Phase distribution:

Row[{Histogram[Flatten[Mod[Arg[hadamard2]/2*π,1]],{0.01},Frame->True,FrameTicksStyle->Directive[Black,14],FrameLabel->{Style["Normalized phases",Black,16],Style["Multiplicity",Black,16]},ImageSize->400,PlotLabel->Style["What we found",18]],
Histogram[Flatten[Mod[Arg[fourierMat]/2*π,1]],{0.01},Frame->True,FrameTicksStyle->Directive[Black,14],FrameLabel->{Style["Normalized phases",Black,16],Style["Multiplicity",Black,16]},ImageSize->400,PlotLabel->Style["Fourier",18]]},"        "]

RRR - 3rd run:

x=RandomComplex[{-1-I,1+I},{n,n}];
error={};

Do[
proj1=unitaryProj[x];
proj2=magnitudeProj[2*proj1-x];
diff=proj2-proj1;
x+=β*diff;
AppendTo[error,norm[diff]];
If[Last[error]<(10^(-10)),Break[]];
,{50000}]//AbsoluteTiming

The error decreases exponentially:

ListLogPlot[error,Joined->True,PlotRange->All,Frame->True,FrameTicksStyle->Directive[Black,12],FrameLabel->{Style["Iteration",Black,16],Style["Error",Black,16]},ImageSize->400]

The Hadamard matrix we found:

Sqrt[n]*(hadamard3=Chop@proj1)//MatrixForm

Phase distribution:

Row[{Histogram[Flatten[Mod[Arg[hadamard3]/2*π,1]],{0.01},Frame->True,FrameTicksStyle->Directive[Black,14],FrameLabel->{Style["Normalized phases",Black,16],Style["Multiplicity",Black,16]},ImageSize->400,PlotLabel->Style["What we found",18]],
Histogram[Flatten[Mod[Arg[fourierMat]/2*π,1]],{0.01},Frame->True,FrameTicksStyle->Directive[Black,14],FrameLabel->{Style["Normalized phases",Black,16],Style["Multiplicity",Black,16]},ImageSize->400,PlotLabel->Style["Fourier",18]]},"        "]

RRR - 4th run:

x=RandomComplex[{-1-I,1+I},{n,n}];
error={};

Do[
proj1=unitaryProj[x];
proj2=magnitudeProj[2*proj1-x];
diff=proj2-proj1;
x+=β*diff;
AppendTo[error,norm[diff]];
If[Last[error]<(10^(-10)),Break[]];
,{50000}]//AbsoluteTiming

The error decreases exponentially:

ListLogPlot[error,Joined->True,PlotRange->All,Frame->True,FrameTicksStyle->Directive[Black,12],FrameLabel->{Style["Iteration",Black,16],Style["Error",Black,16]},ImageSize->400]

The Hadamard matrix we found:

Sqrt[n]*(hadamard4=Chop@proj1)//MatrixForm

Phase distribution:

Row[{Histogram[Flatten[Mod[Arg[hadamard4]/2*π,1]],{0.01},Frame->True,FrameTicksStyle->Directive[Black,14],FrameLabel->{Style["Normalized phases",Black,16],Style["Multiplicity",Black,16]},ImageSize->400,PlotLabel->Style["What we found",18]],
Histogram[Flatten[Mod[Arg[fourierMat]/2*π,1]],{0.01},Frame->True,FrameTicksStyle->Directive[Black,14],FrameLabel->{Style["Normalized phases",Black,16],Style["Multiplicity",Black,16]},ImageSize->400,PlotLabel->Style["Fourier",18]]},"        "]

Not only are the above phase distributions different from that of the Fourier matrix, they are also different from one another:

Row[Histogram[Flatten[Mod[Arg[ToExpression["hadamard"<>ToString[#]]]/2*π,1]],{0.01},Frame->True,FrameTicksStyle->Directive[Black,14],FrameLabel->{Style["Normalized phases",Black,16],Style["Multiplicity",Black,16]},ImageSize->400,PlotLabel->Style["Hadamard-"<>ToString[#],18]]&/@Range[4],"        "]

Coloring the queens graph

We have (n^2) hostile queens placed on an n×n checkerboard. They can attack one another by moving any number of places in the vertical, horizontal, or diagonal direction. But, for some unknown reason, queen A would only attack queen B if A and B have the same color. Thus, to maintain peace, we need to color the checkerboard in such a way that no two queens in a mutually attacking position have the same color.

This is equivalent to coloring a graph with (n^2) vertices, where each vertex represents the position of a queen. Two vertices are connected by an edge if they are in a mutually attacking position. Let us denote the set of all such pairs by  E. We need to assign colors such that for all  (i,j)ϵ E,  vertices  i and  j  have different colors.

We have k different colors available, each of which can be represented by a vector. In order to preserve full permutation symmetry between the colors, we encode them by (k-1)-dimensional vectors which give the vertices of a simplex. We denote the color assigned to vertex j by  (Subscript[(Overscript[v, →]), j]) and consider the matrix (Subscript[X, i,j])=(Subscript[(Overscript[v, →]), i]).(Subscript[(Overscript[v, →]), j]).

The color vectors are normalized, so the diagonal elements must be unity, (Subscript[X, i,i])=1. We also must have  (Subscript[X, i,j])=-1/(k-1)  ∀(i,j)ϵ E  to prevent the queens from starting a fight. Below we define a projection to satisfy these two constraints.

n=6;
grid=Flatten[Table[{i,j},{i,n},{j,n}],1];
diag={#,#}&/@Range[(n^2)];
edges=Complement[Position[Outer[Block[{s=#1-#2},Count[s,0]>0||SameQ@@Abs[s]]&,grid,grid,1],True],diag];
others=Complement[Flatten[Table[{i,j},{i,(n^2)},{j,(n^2)}],1],diag,edges];

k=7;c=-1/(k-1);avg=(1+c)/2;

queenProj[x_]:=Block[{y=MapAt[If[#>avg,1,c]&,x,others]},ReplacePart[y,{diag->1,edges->c}]]

In addition, since (Subscript[X, i,j])=(Subscript[(Overscript[v, →]), i]).(Subscript[(Overscript[v, →]), j])∀i,j, where (Subscript[(Overscript[v, →]), j])’s are (k-1)-dimensional real vectors, X must be a real, symmetric, positive-definite matrix with rank ≤k-1. We implement this rank projection by retaining the (k-1) largest positive eigenvalues in the eigen-decomposition of X:

rankProj[x_]:=
Block[{val,vec},
{val,vec}=Transpose[Select[Take[Sort[Transpose[Eigensystem[x]]],-k+1],First[#]>0&]];
Transpose[vec].DiagonalMatrix[val].vec
]

We are ready to run RRR:

norm[list_]:=Sqrt[Mean[((Flatten[list])^2)]];

β=0.5;

x=Block[{y=RandomReal[{-1,1},{(n^2),(n^2)}]},y+Transpose[y]];
error={};

Do[
proj1=queenProj[x];
proj2=rankProj[2*proj1-x];
diff=proj2-proj1;
x+=β*diff;
AppendTo[error,norm[diff]];
If[Last[error]<(10^(-10)),Break[]];
,{5000}]//AbsoluteTiming

ListLogPlot[error,Joined->True,PlotRange->All,Frame->True,FrameTicksStyle->Directive[Black,12],FrameLabel->{Style["Iteration",Black,16],Style["Error",Black,16]},ImageSize->400]

Check that the two constraints are indeed satisfied:

{Norm[proj1[[First[#],Last[#]]]&/@diag-1],Norm[proj1[[First[#],Last[#]]]&/@edges-c]}

Eigenvalues[proj1]

The color vectors at each vertex:

colorvecs=Transpose[Eigenvectors[proj1,k-1]]

The k distinct color vectors:

DeleteDuplicates[colorvecs]

We simply label the above colors 1,2,...,k and find the labels at each vertex:

colors=Position[DeleteDuplicates[colorvecs],#][[1,1]]&/@colorvecs

Here is how the queens are colored!

Graphics[Join[{Raster[Partition[colors/k ,n],ColorFunction->"Rainbow"]},
Text[Style[ToString@Partition[colors ,n][[Last[#],First[#]]],18],#-{1,1}/2]&/@grid]]

Next we try a 7×7 grid:

n=7;
grid=Flatten[Table[{i,j},{i,n},{j,n}],1];
diag={#,#}&/@Range[(n^2)];
edges=Complement[Position[Outer[Block[{s=#1-#2},Count[s,0]>0||SameQ@@Abs[s]]&,grid,grid,1],True],diag];
others=Complement[Flatten[Table[{i,j},{i,(n^2)},{j,(n^2)}],1],diag,edges];

β=0.5;

x=Block[{y=RandomReal[{-1,1},{(n^2),(n^2)}]},y+Transpose[y]];
error={};

Do[
proj1=queenProj[x];
proj2=rankProj[2*proj1-x];
diff=proj2-proj1;
x+=β*diff;
AppendTo[error,norm[diff]];
If[Last[error]<(10^(-10)),Break[]];
,{5000}]//AbsoluteTiming

ListLogPlot[error,Joined->True,PlotRange->All,Frame->True,FrameTicksStyle->Directive[Black,12],FrameLabel->{Style["Iteration",Black,16],Style["Error",Black,16]},ImageSize->400]

Check that the two constraints are indeed satisfied:

{Norm[proj1[[First[#],Last[#]]]&/@diag-1],Norm[proj1[[First[#],Last[#]]]&/@edges-c]}

Eigenvalues[proj1]

Colors at each vertex:

colorvecs=Transpose[Eigenvectors[proj1,k-1]];
colors=Position[DeleteDuplicates[colorvecs],#][[1,1]]&/@colorvecs

How the queens are colored:

Graphics[Join[{Raster[Partition[colors/k ,n],ColorFunction->"Rainbow"]},
Text[Style[ToString@Partition[colors ,n][[Last[#],First[#]]],18],#-{1,1}/2]&/@grid]]

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

Survey questions:

1. I prefer Mathematica demonstrations, since that’s the only language I’m comfortable with!
2. I’m not particularly motivated to learn any programming language unless I need it to solve problems.
3. Origami, perhaps! I thought we had a good exposure to a diverse set of applications.
4. Maybe “projpy,” just because “proxpy” has already been taken:  (HyperlinkURL["https://github.com/yeukhon/proxpy", https://github.com/yeukhon/proxpy]).