Examples file for QI Mathematica package (https://github.com/rogercolbeck/QI)

Generate a maximally entangled state of two qubits

ent=(1/2)^(1/2)*(KetV[0,2]⊗KetV[0,2]+KetV[1,2]⊗KetV[1,2])

Make a density operator from this state

rho12=DM[ent];rho12//MatrixForm

Compute the reduced density operator on the first subsystem

rho1=PT[rho12,{1,0},{2,2}];rho1//MatrixForm

And on the second

rho2=PT[rho12,{0,1},{2,2}];rho2//MatrixForm

This can also be done using QubitPartialTrace

QubitPartialTrace[rho12,{2}]//MatrixForm

Purify rho1

Purify[rho1]

Apply the tensor product of two Pauli operators to the state

(σ[3]⊗σ[2]).ent//MatrixForm

Take another pure state of dimension 8

psi={{(Sqrt[(2/33)])},{(Sqrt[(2/33)])},{2 (Sqrt[(2/33)])},{(Sqrt[(3/22)])},{(Sqrt[(3/22)])},{(Sqrt[(2/33)])},{(Sqrt[(2/33)])},{2 (Sqrt[(2/33)])}};

Find its Schmidt decomposition assuming it is a 2D system tensored with a 4D system

{c,vecs1,vecs2}=FullSimplify[SchmidtDecomposition[psi,{2,4}]]

Confirm that this is correct

FullSimplify[Sum[c[[i]]*vecs1[[i]]⊗vecs2[[i]],{i,1,Length[c]}]-psi]

Take the tensor product of psi with ent, where psi is on systems 1 and 3 and ent is 2 and 4

tensored=FullSimplify[Tensor[psi,ent,{1,2,1,2},{2,2,4,2}]]

Tracing out systems 1 and 3 should recover the density matrix of ent

FullSimplify[PT[DM[tensored],{0,1,0,1},{2,2,4,2}]]

We can also swap the systems in tensored so as to get ent⊗psi

ExchangeSystems[tensored,{3,1,4,2},{2,2,4,2}]-ent⊗psi

Simultaneously diagonalize two commuting matrices

a=σ[1]⊗σ[1];b=σ[2]⊗σ[2];a.b-b.a

u=SimultaneouslyDiagonalize[a,b]

CT[u].a.u//MatrixForm

CT[u].b.u//MatrixForm

Treat ent as an isometry and complete to a unitary

unitary=AppendCols[ent];unitary//MatrixForm

unitary.CT[unitary]

Find the RQ decomposition of a random density operator:

{r,q}=RQDecompNeg[rho=PickRandomRho[4]];

Check the result, and that the matrices have the correct form:

Chop[r.CT[q]-rho]//MatrixForm

Chop[r]//MatrixForm

Chop[q.CT[q]]//MatrixForm

Pick a random channel from dimension 2 to 4 with 10 Kraus operators

channel=PickRandomChannel[2,4,10];

Length[channel]

Check that it is a valid channel

Chop[Sum[CT[channel[[i]]].channel[[i]],{i,1,10}]]//MatrixForm

Act on density operator rho

rho={{0.25,0},{0,0.75}};

rhop=Chop[Sum[channel[[i]].rho.CT[channel[[i]]],{i,1,Length[channel]}]]

Find the Choi state of the channel

choist=ChoiState[channel];

Returning to the channel will give a smaller Kraus rank (8 is the maximum needed)

channel2=ChoiChannel[choist,2,4][[1]];

Length[channel2]

Check that it has the same action on rho

rhopp=Chop[Sum[channel2[[i]].rho.CT[channel2[[i]]],{i,1,Length[channel2]}]]

Chop[rhopp-rhop]

This is the same as ChannelCompress

Chop[ChannelCompress[channel]-channel2]

Find the conditional von Neumann entropy of the first qubit with respect to the second in rhop

vNEntropy[rhop,{1,2},{2,2}]

This is the same as the joint entropy minus the reduced entropy

vNEntropy[rhop,{1},{4}]-vNEntropy[PT[rhop,{0,1},{2,2}],{1},{2}]

Pick a random state of 4 qubits with real coefficients

psi=RPickRandomPsi[16]

Compute the conditional mutual information between the 3rd and 4th systems given the first

Chop[vNInfo[DM[psi],{3,0,2,1},{2,2,2,2}]]

There is symmetry in the labels 1 and 2

Chop[vNInfo[DM[psi],{3,0,1,2},{2,2,2,2}]]

It is also invariant under unitaries on the system being conditioned on

Chop[vNInfo[DM[(PickRandomUnitary[2]⊗IdentityMatrix[8]).psi],{3,0,1,2},{2,2,2,2}]]

Since psi is pure, it should have no entropy

Chop[vNEntropy[DM[psi],{1},{16}]]

Pick two random density matrices and compute the relative entropy between them

rho1=PickRandomRho[4];

rho2=PickRandomRho[4];

Chop[RelEnt[rho1,rho2]]

Now choose rho2 without full support

rho2=PickRandomRho[2]⊗{{1,0},{0,0}};

RelEnt[rho1,rho2]

Chop[RelEnt[rho2,rho1]]

Now compute the alpha=1.5 Renyi entropy of the same states.

Chop[RelEnt[1.5,rho2,rho1]]

Create an 8x8 matrix with symbolic values, treating it as three qubits and compute the partial trace where the first and third systems are traced out

rho=CreateMatrix[m,8,8];rho//MatrixForm

PT[rho,{0,1,0},{2,2,2}]//MatrixForm

Pick two random 16 dimensional state (4 qubits) with exact values

psi1=FPickRandomPsi[16,2];psi2=FPickRandomPsi[16,2]

Compute their trace distance and fidelity

Simplify[Dist[DM[psi1],DM[psi2]]]

N[%]

Simplify[Fidelity[DM[psi1],DM[psi2]]]

N[%]

On tracing out the first, third and fourth systems, the distance should decrease and the fidelity increase

FullSimplify[Dist[Simplify[PT[DM[psi1],{0,1,0,0},{2,2,2,2}]],Simplify[PT[DM[psi2],{0,1,0,0},{2,2,2,2}]]]]

N[%]

FullSimplify[Fidelity[Simplify[PT[DM[psi1],{0,1,0,0},{2,2,2,2}]],Simplify[PT[DM[psi2],{0,1,0,0},{2,2,2,2}]]]]

N[%]

Note that QuickDist is faster than Dist in general

rho1={{(162/713),(105/(713 (Sqrt[521]))),(2565/(31 (Sqrt[3966373]))),-(180/(31 (Sqrt[7613])))},{(105/(713 (Sqrt[521]))),(107602/371473),-(5985/(16151 (Sqrt[7613]))),(420/(31 (Sqrt[3966373])))},{(2565/(31 (Sqrt[3966373]))),-(5985/(16151 (Sqrt[7613]))),(1405854/5345981),(10260/(10261 (Sqrt[521])))},{-(180/(31 (Sqrt[7613]))),(420/(31 (Sqrt[3966373]))),(10260/(10261 (Sqrt[521]))),(2259/10261)}};rho2={{(5/26),(363/(650 (Sqrt[41]))),-((266 (Sqrt[(2/41)]))/1625),((22 (Sqrt[2]))/1625)},{(363/(650 (Sqrt[41]))),(9837/26650),-((2142 (Sqrt[2]))/66625),-((486 (Sqrt[(2/41)]))/1625)},{-((266 (Sqrt[(2/41)]))/1625),-((2142 (Sqrt[2]))/66625),(1188/5125),(672/(1625 (Sqrt[41])))},{((22 (Sqrt[2]))/1625),-((486 (Sqrt[(2/41)]))/1625),(672/(1625 (Sqrt[41]))),(336/1625)}};

Timing[Chop[N[Dist[rho1,rho2]]]]

Timing[N[QuickDist[rho1,rho2]]]

Choose a random density operator in dimension 4 and a random measurement with 2 outcomes in dimension 4

rho=PickRandomRho[4];

meas=PickRandomMeasurement[4,2];

Check that it is a valid measurement:

PositiveSemidefiniteMatrixQ[CT[meas[[1]]].meas[[1]]]

PositiveSemidefiniteMatrixQ[CT[meas[[2]]].meas[[2]]]

Chop[Sum[CT[meas[[i]]].meas[[i]],{i,1,2}]]

The post-measurement state after applying the measurement to rho and getting the first outcome is

Chop[meas[[1]].rho.CT[meas[[1]]]]

Its trace is the probability of this outcome

Tr[%]

Repeat for the second outcome:

Chop[meas[[2]].rho.CT[meas[[2]]]]

Tr[%]

We can also define a POVM from these

POVM1=Table[CT[meas[[i]]].meas[[i]],{i,1,2}];

Chop[MeasurePOVM[rho,POVM1,{1},{4}]]

We might also want to apply this POVM to half of a maximally entangled state

ent=(1/4)^(1/2)*Sum[KetV[i,4]⊗KetV[i,4],{i,0,3}]

stateOut=Chop[MeasurePOVM[DM[ent],POVM1,{1,2},{4,4}]];stateOut//MatrixForm

The probability of the first outcome is the trace of the upper left 4x4 submatrix:

Tr[Take[stateOut,4,4]]

Chop[Tr[(POVM1[[1]]⊗IdentityMatrix[4]).DM[ent]]]

See the entangled state in basis form, considering it to be a state of 4 qubits

BasisFormS[ent,{2,2,2,2}]

Or the same but treating it as two 4-dimensional systems

BasisFormS[ent,{4,4}]

Or as a 2 and an 8

BasisFormS[ent,{2,8}]

Find the optimum POVM to distinguish two states

rho1=PickRandomRho[3];rho2=PickRandomRho[3];

POVM=Chop[OptimumPOVM[rho1,rho2]]

Calculate the probability of correctly guessing the state using this POVM

Chop[(Tr[POVM[[1]].rho1]+Tr[POVM[[2]].rho2])/2]

This should be the same as 1/2 ( 1 + Dist[rho1, rho2]) :

Chop[1/2*(1+Dist[rho1,rho2])]