<<MaTeX` (* For LaTeX typesetting in Mathematica *)

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

Working with qubits

A general state X is specified by the amplitudes for being in each of the 16 basis states. In the standard basis, these amplitudes are ordered as follows:

X@@@Tuples[{{1,2},{1,2},{1,2},{1,2}}]

The following functions convert states to vectors and vice versa:

stateToVec[state_]:=UnitVector[16,FromDigits[Last[state]-1,2]+1]

stateToVec[state1_+state2_]:=stateToVec[state1]+stateToVec[state2]

stateToVec[c_*state_]:=c*stateToVec[state]

stateToVec[state[{1,1,1,2}]]

stateToVec[(1/2)*state[{1,1,1,2}]+state[{1,2,1,2}]]

vecToState[vec_]:=vec.state/@Tuples[{{1,2},{1,2},{1,2},{1,2}}]

vecToState[RandomReal[1,16]]

There are 12 single-qubit spin operators: (Subscript[σ, j]^i) with i=1,2,3,4 and j=x,y,z. Each of these operators can be represented by a matrix in the 16-dimensional Hilbert space. For example, the spin operator (Subscript[σ, x]^2)=1⊗(Subscript[σ, x])⊗1⊗1 is represented, in the standard basis, by the Kronecker product of the respective single-particle matrices:

MatrixForm[σ[2,1]=KroneckerProduct[IdentityMatrix[2],PauliMatrix[1],IdentityMatrix[2],IdentityMatrix[2]]]

We can verify that (Subscript[σ, x]^2) flips the spin of the 2nd qubit, without affecting the other qubits:

vecToState[σ[2,1].stateToVec[state[{1,1,2,1}]-state[{2,2,1,1}]]]

Similarly, we can construct the matrices for all other single-qubit spin operators:

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

Multi-particle operators can be simply constructed by taking products of the single-qubit operators. Thus, the operator (Subscript[σ, x]^1)(Subscript[σ, y]^3) is given by the following matrix:

σ[1,1].σ[3,2]//MatrixForm

As can be seen, the matrix elements are zero except for a single ±i is each row and column. The positions of the ±i entries are given below:

Position[σ[1,1].σ[3,2],I]

Position[σ[1,1].σ[3,2],-I]

Pizza for (up to) five

The equipartition projection is implemented by appropriately shifting the tensor elements. For a tensor X of dimension (Subscript[n, 1])×(Subscript[n, 2])×...×(Subscript[n, N]), whose elements should sum to S, the equipartition projection (Subscript[P, E])(X) is given as follows:

where  and . Below we define this projection:

equipartitionProj[x_,S_]:=
Block[{dimensions,avgoverslices,scalarterm},
dimensions=Dimensions[x];
avgoverslices=Sum[TensorProduct@@ReplacePart[Table[1,{i,#}]&/@dimensions,j->Mean/@Flatten/@Flatten[x,{j}]],{j,Length[dimensions]}];
scalarterm=(Length[dimensions]-1)*Mean[Flatten[x]]+S/Times@@dimensions;
x-avgoverslices+scalarterm
]

Check this projection on a 3×4 tensor:

With[{S=12},MatrixForm[equipartitionCheck=equipartitionProj[RandomReal[{-1,1},{3,4}],S]]]

All rows add up to 4 and all columns add up to 3:

{Total/@equipartitionCheck,Total@equipartitionCheck}

Next we define the sparsity projection (Subscript[P, S]) which yields a tensor containing a maximum of p nonzero elements which are all positive:

sparsityProj[x_,p_]:=
Block[{y},
y=Flatten[Clip[x,{0,Infinity}]];
y[[Ordering[y,Times@@Dimensions[x]-p]]]=0.;
ArrayReshape[y,Dimensions[x]]]

Check this projection on a 3×4 tensor:

With[{p=6},Block[{y=RandomReal[{-1,0.5},{3,4}]},{MatrixForm[y],MatrixForm[sparsityProj[y,p]]}]]

We will first apply the RRR algorithm to the simpler “pizza for up to four” problem:

β=0.5;

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

nmax=4;
size=Range[Floor[nmax/2]+1,nmax];
S=LCM@@size;
p=6;

x=RandomReal[1,size];
error={};

Do[
proj1=equipartitionProj[x,S];
proj2=sparsityProj[2*proj1-x,p];
diff=proj2-proj1;
x+=β*diff;
AppendTo[error,norm[diff]];
If[Last[error]<(10^(-12)),Break[]];
,{1000}]//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]

We find the expected solution:

Row[MatrixForm/@Chop[{proj1,proj2,x}],"    "]

Now it is time to tackle the original problem:

1st run:

nmax=5;
size=Range[Floor[nmax/2]+1,nmax];
S=LCM@@size;
p=9;

x=RandomReal[1,size];
error={};

Do[
proj1=equipartitionProj[x,S];
proj2=sparsityProj[2*proj1-x,p];
diff=proj2-proj1;
x+=β*diff;
AppendTo[error,norm[diff]];
If[Last[error]<(10^(-12)),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]

We find the following solution:

pizzaproj1=Chop@proj1

Piece sizes:

pieces1=SparseArray[pizzaproj1]["NonzeroValues"]

How the pieces are shared by 3/4/5 people:

divisions1[3]=pieces1[[#]]&/@Flatten[Position[SparseArray[pizzaproj1]["NonzeroPositions"],{#,_,_}]]&/@Range[3]

divisions1[4]=pieces1[[#]]&/@Flatten[Position[SparseArray[pizzaproj1]["NonzeroPositions"],{_,#,_}]]&/@Range[4]

divisions1[5]=pieces1[[#]]&/@Flatten[Position[SparseArray[pizzaproj1]["NonzeroPositions"],{_,_,#}]]&/@Range[5]

Everyone gets the same amount of pizza:

{Plus@@@divisions1[3],Plus@@@divisions1[4],Plus@@@divisions1[5]}

2nd run:

x=RandomReal[1,size];
error={};

Do[
proj1=equipartitionProj[x,S];
proj2=sparsityProj[2*proj1-x,p];
diff=proj2-proj1;
x+=β*diff;
AppendTo[error,norm[diff]];
If[Last[error]<(10^(-12)),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]

We find the dynamical system is trapped in a cycle.

3rd run:

x=RandomReal[1,size];
error={};

Do[
proj1=equipartitionProj[x,S];
proj2=sparsityProj[2*proj1-x,p];
diff=proj2-proj1;
x+=β*diff;
AppendTo[error,norm[diff]];
If[Last[error]<(10^(-12)),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]

Now we get a different solution:

pizzaproj2=Chop@proj1

Piece sizes:

pieces2=SparseArray[pizzaproj2]["NonzeroValues"]

How the pieces are shared by 3/4/5 people:

divisions2[3]=pieces2[[#]]&/@Flatten[Position[SparseArray[pizzaproj2]["NonzeroPositions"],{#,_,_}]]&/@Range[3]

divisions2[4]=pieces2[[#]]&/@Flatten[Position[SparseArray[pizzaproj2]["NonzeroPositions"],{_,#,_}]]&/@Range[4]

divisions2[5]=pieces2[[#]]&/@Flatten[Position[SparseArray[pizzaproj2]["NonzeroPositions"],{_,_,#}]]&/@Range[5]

Everyone gets the same amount of pizza:

{Plus@@@divisions2[3],Plus@@@divisions2[4],Plus@@@divisions2[5]}

Latin squares by divide and concur

We’ll solve the problem in several stages. First, we consider the a simple 4×4 latin square problem, to which we know the solution.

incompleteSquare={{1,0,0,2},{0,2,0,0},{0,4,3,0},{3,1,0,4}};
num=Length[incompleteSquare];
givencoords=Flatten[Table[Append[#,j]&/@Position[incompleteSquare,j],{j,num}],1];
sqgrid=Tuples[{Range[num],Range[num]}];
keys=
Table[Block[{positionrules},
positionrules=ArrayRules[Flatten[SparseArray[Drop[#,-1]->Last[#]&/@RotateLeft[givencoords,{0,j-1}],num]]];
Through[{Keys,Values}[Drop[positionrules,-1]]]
],{j,3}];

incompleteSquare//Grid

To implement the “divide” projection (Subscript[P, D]), we set the largest element in each row/column/stack equal to 1 and set the other elements to 0:

divideProj[y_]:=
SparseArray[
#->1&/@
Flatten[Table[
Block[{maxpositions,coordinates},
maxpositions=Flatten[Ordering[#,-1]&/@Flatten[TensorTranspose[y[[j]],RotateRight[Range[3],j-1]],1]];
coordinates=RotateRight[MapThread[Append[#1,#2]&,{sqgrid,maxpositions}],{0,j-1}];
Prepend[#,j]&/@coordinates],
{j,3}],1]
]

To implement the “concur” projection (Subscript[P, C]), we compute the average of the 3 copies and set each element to 0 if it is less than 1/2 and to 1 otherwise:

concurProj[y_]:=ConstantArray[UnitStep[Total[y]-1.5],3]

Next we modify each projection to incorporate the positions of the given elements. We simply set the given coordinates to 1:

divideProjConstrained[y_]:=
SparseArray[
#->1&/@
Flatten[Table[
Block[{maxpositions,unitypositions,coordinates},
maxpositions=Flatten[Ordering[#,-1]&/@Delete[Flatten[TensorTranspose[y[[j]],RotateRight[Range[3],j-1]],1],keys[[j,1]]]];
unitypositions=multiInsert[maxpositions,keys[[j,2]],Flatten[keys[[j,1]]]];
coordinates=RotateRight[MapThread[Append[#1,#2]&,{sqgrid,unitypositions}],{0,j-1}];
Prepend[#,j]&/@coordinates],
{j,3}],1]
]

concurProjConstrained[y_]:=ConstantArray[MapAt[1&,UnitStep[Total[y]-1.5],givencoords],3]

We apply the RRR algorithm to complete the 4×4 latin square:

y=RandomReal[1,Join[{3},Table[num,{j,3}]]];
error={};

Do[
proj1=divideProjConstrained[y];
proj2=concurProj[2*proj1-y];
diff=proj2-proj1;
y+=β*diff;
AppendTo[error,norm[diff]];
If[Last[error]<(10^(-12)),Break[]];
,{1000}]//AbsoluteTiming

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

Here is the completed square:

Grid[Normal[SparseArray[Drop[#,-1]->Last[#]&/@SparseArray[First[proj2]]["NonzeroPositions"]]]]

Next we move up to a 7×7 square. We’ll use this example to compare the different possible RRR schemes and study the variation with the time step β.

incompleteSquare=Normal[SparseArray[{{2,2}->2,{2,5}->7,{3,3}->5,{3,5}->4,{5,4}->5,{7,2}->4},7]];
num=Length[incompleteSquare];
givencoords=Flatten[Table[Append[#,j]&/@Position[incompleteSquare,j],{j,num}],1];
sqgrid=Tuples[{Range[num],Range[num]}];
keys=
Table[Block[{positionrules},
positionrules=ArrayRules[Flatten[SparseArray[Drop[#,-1]->Last[#]&/@RotateLeft[givencoords,{0,j-1}],num]]];
Through[{Keys,Values}[Drop[positionrules,-1]]]
],{j,3}];

incompleteSquare//Grid

We’ll apply the RRR scheme to complete the above square. We have 4 different projections at our disposal: (i) divide (Subscript[P, D]), (ii) concur (Subscript[P, C]), (iii) constrained divide (Subscript[(OverTilde[P]), D]), and (iv) constrained concur (Subscript[(OverTilde[P]), C]). We can combine them in 4 different ways:

1. (OverDot[y])=(Subscript[P, C])(2(Subscript[(OverTilde[P]), D])(y)-y)-(Subscript[(OverTilde[P]), D])(y)

β=0.1;

y=RandomReal[1,Join[{3},Table[num,{j,3}]]];
error={};

Do[
proj1=divideProjConstrained[y];
proj2=concurProj[2*proj1-y];
diff=proj2-proj1;
y+=β*diff;
AppendTo[error,norm[diff]];
If[Last[error]<(10^(-5)),Break[]];
,{1000}]//AbsoluteTiming

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

Row[{
Column[{Style["Incomplete",12],Grid[incompleteSquare]},Alignment->Center],
Column[{Style["Complete",12],Grid[Normal[SparseArray[Drop[#,-1]->Last[#]&/@SparseArray[First[proj1]]["NonzeroPositions"]]]]},Alignment->Center]},
"         "]

2. (OverDot[y])=(Subscript[(OverTilde[P]), C])(2(Subscript[P, D])(y)-y)-(Subscript[P, D])(y)

y=RandomReal[1,Join[{3},Table[num,{j,3}]]];
error={};

Do[
proj1=divideProj[y];
proj2=concurProjConstrained[2*proj1-y];
diff=proj2-proj1;
y+=β*diff;
AppendTo[error,norm[diff]];
If[Last[error]<(10^(-5)),Break[]];
,{1000}]//AbsoluteTiming

The error decreases to zero in the same timescale:

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

We find a different solution -- this is simply because there are multiple solutions and we start from a random point:

Row[{
Column[{Style["Incomplete",12],Grid[incompleteSquare]},Alignment->Center],
Column[{Style["Complete",12],Grid[Normal[SparseArray[Drop[#,-1]->Last[#]&/@SparseArray[First[proj1]]["NonzeroPositions"]]]]},Alignment->Center]},
"         "]

3. (OverDot[y])=(Subscript[(OverTilde[P]), D])(2(Subscript[P, C])(y)-y)-(Subscript[P, C])(y)

y=RandomReal[1,Join[{3},Table[num,{j,3}]]];
error={};

Do[
proj1=concurProj[y];
proj2=divideProjConstrained[2*proj1-y];
diff=proj2-proj1;
y+=β*diff;
AppendTo[error,norm[diff]];
If[Last[error]<(10^(-5)),Break[]];
,{1000}]//AbsoluteTiming

We see a different behavior -- the error keeps fluctuating and does not approach zero (in the same timescale as before):

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

4. (OverDot[y])=(Subscript[P, D])(2(Subscript[(OverTilde[P]), C])(y)-y)-(Subscript[(OverTilde[P]), C])(y)

y=RandomReal[1,Join[{3},Table[num,{j,3}]]];
error={};

Do[
proj1=concurProjConstrained[y];
proj2=divideProj[2*proj1-y];
diff=proj2-proj1;
y+=β*diff;
AppendTo[error,norm[diff]];
If[Last[error]<(10^(-5)),Break[]];
,{1000}]//AbsoluteTiming

Again, the error keeps fluctuating:

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

Thus, we find that the order of the projections make a big difference here. The first two combinations are much better at finding solutions.

Effect of changing the time step parameter β

A larger time step causes more fluctuations, requiring more iterations to find the solution:

β=0.8;

y=RandomReal[1,Join[{3},Table[num,{j,3}]]];
error={};

Do[
proj1=divideProjConstrained[y];
proj2=concurProj[2*proj1-y];
diff=proj2-proj1;
y+=β*diff;
AppendTo[error,norm[diff]];
If[Last[error]<(10^(-5)),Break[]];
,{10000}]//AbsoluteTiming

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

Choosing a very small β reduces fluctuations, but also needs more time steps to reach the solution:

β=0.02;

y=RandomReal[1,Join[{3},Table[num,{j,3}]]];
error={};

Do[
proj1=divideProjConstrained[y];
proj2=concurProj[2*proj1-y];
diff=proj2-proj1;
y+=β*diff;
AppendTo[error,norm[diff]];
If[Last[error]<(10^(-5)),Break[]];
,{10000}]//AbsoluteTiming

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

The sweet spot seems to be around β=0.1:

β=0.1;

y=RandomReal[1,Join[{3},Table[num,{j,3}]]];
error={};

Do[
proj1=divideProjConstrained[y];
proj2=concurProj[2*proj1-y];
diff=proj2-proj1;
y+=β*diff;
AppendTo[error,norm[diff]];
If[Last[error]<(10^(-5)),Break[]];
,{1000}]//AbsoluteTiming

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

The 25×25 latin square completion

incompleteSquare=Import["https://raw.githubusercontent.com/veitelser/RRR_mini_course/master/data/25x25_latin_square","Table"];
num=Length[incompleteSquare];
givencoords=Flatten[Table[Append[#,j]&/@Position[incompleteSquare,j],{j,num}],1];
sqgrid=Tuples[{Range[num],Range[num]}];
keys=
Table[Block[{positionrules},
positionrules=ArrayRules[Flatten[SparseArray[Drop[#,-1]->Last[#]&/@RotateLeft[givencoords,{0,j-1}],num]]];
Through[{Keys,Values}[Drop[positionrules,-1]]]
],{j,3}];

incompleteSquare//Grid

We try the RRR scheme  y'=y+β*((Subscript[P, C])(2(Subscript[(OverTilde[P]), D])(y)-y)-(Subscript[(OverTilde[P]), D])(y))  with β=0.1:

y=RandomReal[1,Join[{3},Table[num,{j,3}]]];
error={};

Do[
proj1=divideProjConstrained[y];
proj2=concurProj[2*proj1-y];
diff=proj2-proj1;
y+=β*diff;
AppendTo[error,norm[diff]];
If[Last[error]<(10^(-5)),Break[]];
,{25000}]//AbsoluteTiming

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

Here is the solution:

Column[
{Style["Completed Latin Square",15],Grid[latinsquare1=Normal[SparseArray[Drop[#,-1]->Last[#]&/@SparseArray[First[proj1]]["NonzeroPositions"]]]]
},Alignment->Center,Spacings->1]

We also try  y'=y+β*((Subscript[(OverTilde[P]), C])(2(Subscript[P, D])(y)-y)-(Subscript[P, D])(y))  with β=0.1:

y=RandomReal[1,Join[{3},Table[num,{j,3}]]];
error={};

Do[
proj1=divideProj[y];
proj2=concurProjConstrained[2*proj1-y];
diff=proj2-proj1;
y+=β*diff;
AppendTo[error,norm[diff]];
If[Last[error]<(10^(-5)),Break[]];
,{25000}]//AbsoluteTiming

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

And we find another solution:

Column[
{Style["Completed Latin Square",15],Grid[latinsquare2=Normal[SparseArray[Drop[#,-1]->Last[#]&/@SparseArray[First[proj1]]["NonzeroPositions"]]]]
},Alignment->Center,Spacings->1]

The two solutions are different!

norm[latinsquare1-latinsquare2]

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