Working with Tensors III

This notebook continues on developing the tools needed to work with Tensors in Mathematica

Routines from first two notebooks

Reordering Indices

NamedTranspose[mat_,orig_,final_]:=Module[{ord=Flatten[Position[final,#]&/@orig]},Transpose[mat,ord]]

NamedTranspose::usage="NamedTranspose[mat_,orig_,final_] takes a tensor `mat' whose indices are named `orig' and shuffles the order, so they match `final'.  For example NamedTranspose[SparseArray[{1,2,3}->5,{3,3,3}],{'a','b','c'},{'b','c','a'}]";

ReOrder[tensordescription_List,finalnames_]:=Module[{
tensor="tensor"/.tensordescription, (* Extract Tensor *)
names="names"/.tensordescription, (* Extract names of indices *)
newtensor (* Declare local variables *)
},
(* Calculate new Tensor *)
newtensor=NamedTranspose[tensor,names,finalnames];
(* Assemble into rules *)
{"tensor"->newtensor,"names"->finalnames}
]

ReOrder::usage="ReOrder[tensordescription_List,finalnames_] takes a NamedTensor object of the form {'tensor'->tensor,'names'->names} and returns a NamedTensor, with the indices permuted to match finalnames"

Combining Indices

CombineIndices[namedtensor_List,indicestocombine_List,newindex_]:=
Module[{oldnames="names"/.namedtensor,names,targpos,newnamedtensor,dims,tensor,newdims,newtensor,finalnamedtensor,ind,indexkeys,refine,newnames,targnames},
(* Shift the indices we want to the right *)
targpos=First@Flatten@Position[oldnames,Alternatives@@indicestocombine];
names=Join[Complement[oldnames,indicestocombine],indicestocombine];
newnamedtensor=ReOrder[namedtensor,names];
(* Next extract the tensor *)
tensor="tensor"/.newnamedtensor;
(* Calculate the lengths of each of the indices *)
dims=Dimensions[tensor];
newdims=Append[dims[[;;-Length[indicestocombine]-1]],Times@@dims[[-Length[indicestocombine];;]]];
(* Do the reshaping *)
newtensor=ArrayReshape[tensor,newdims];
newnames=Append[names[[;;-Length[indicestocombine]-1]],newindex];
(* Do a final reordering *)
targnames=Join[newnames[[;;targpos-1]],{newnames[[-1]]},newnames[[targpos;;-2]]];
ReOrder[{"tensor"->newtensor,"names"->newnames},targnames]]

rightCombineIndices[namedtensor_List,indicestocombine_List,newindex_]:=
Module[{oldnames="names"/.namedtensor,names,newnamedtensor,dims,tensor,newdims,newtensor,finalnamedtensor,ind,indexkeys,refine,newnames,targnames},
(* Shift the indices we want to the right *)
names=Join[SortBy[Complement[oldnames,indicestocombine],Position[oldnames,#]&],indicestocombine];
newnamedtensor=ReOrder[namedtensor,names];
(* Next extract the tensor *)
tensor="tensor"/.newnamedtensor;
(* Calculate the lengths of each of the indices *)
dims=Dimensions[tensor];
newdims=Append[dims[[;;-Length[indicestocombine]-1]],Times@@dims[[-Length[indicestocombine];;]]];
(* Do the reshaping *)
newtensor=ArrayReshape[tensor,newdims];
newnames=Append[names[[;;-Length[indicestocombine]-1]],newindex];
{"tensor"->newtensor,"names"->newnames}]

leftCombineIndices[namedtensor_List,indicestocombine_List,newindex_]:=
Module[{oldnames="names"/.namedtensor,names,newnamedtensor,dims,tensor,newdims,newtensor,finalnamedtensor,ind,indexkeys,refine,newnames,targnames},
(* Shift the indices we want to the right *)
names=Join[indicestocombine,SortBy[Complement[oldnames,indicestocombine],Position[oldnames,#]&]];
newnamedtensor=ReOrder[namedtensor,names];
(* Next extract the tensor *)
tensor="tensor"/.newnamedtensor;
(* Calculate the lengths of each of the indices *)
dims=Dimensions[tensor];
newdims=Prepend[dims[[Length[indicestocombine]+1;;]],Times@@dims[[;;Length[indicestocombine]]]];
(* Do the reshaping *)
newtensor=ArrayReshape[tensor,newdims];
newnames=Prepend[names[[Length[indicestocombine]+1;;]],newindex];
{"tensor"->newtensor,"names"->newnames}]

makename[indicestocombine_List]:=StringJoin@@Riffle[ToString/@indicestocombine,Table["_",{Length[indicestocombine]-1}]]

breakname[compoundname_String]:=StringSplit[compoundname,"_"]

CombineIndices[namedtensor_List,indicestocombine_List]:=
Module[{newindex=makename[indicestocombine]},
CombineIndices[namedtensor,indicestocombine,newindex]]

leftCombineIndices[namedtensor_List,indicestocombine_List]:=
Module[{newindex=makename[indicestocombine]},
leftCombineIndices[namedtensor,indicestocombine,newindex]]

rightCombineIndices[namedtensor_List,indicestocombine_List]:=
Module[{newindex=makename[indicestocombine]},
rightCombineIndices[namedtensor,indicestocombine,newindex]]

BreakIndices[namedtensor_List,indextobreak_,newnames_List,newdimensions_List]:=
Module[{tensor="tensor"/.namedtensor,
names="names"/.namedtensor,
olddims,
pos,newtensor,newdims,finalnames,dimtobereplaced
},
olddims=Dimensions[tensor];
pos=First@Flatten@Position[names,indextobreak];
(* Do a quick consistency check *)
dimtobereplaced=olddims[[pos]];
If[dimtobereplaced≠Times@@newdimensions,
Return[Failure["Inconsistent Dimensions",
Association["MessageTemplate"->"Trying to break an index of length `l1` into indices of length `dims`",
"MessageParameters"->Association[
"l1"->dimtobereplaced,
"dims"->newdimensions
]]]];
];
newdims=Join[olddims[[;;pos-1]],newdimensions,olddims[[pos+1;;]]];
newtensor=ArrayReshape[tensor,newdims];
finalnames=Join[names[[;;pos-1]],newnames,names[[pos+1;;]]];
{"tensor"->newtensor,"names"->finalnames}
]

Outer Products

disambiguate[names_List]:=Module[{seive,count},
count[_]=0;
seive[a_]:=Module[{astr=ToString[a],root,numticks=0},
root=astr;
While[StringTake[root,-1]=="'",
numticks+=1;
root=StringTake[root,{1,-2}];
];
count[root]+=1;
root<>StringRepeat["'",count[root]-1]];
seive/@names]

Combine[namedtensors__List]:={"tensor"->TensorProduct@@("tensor"/.{namedtensors}),"names"->disambiguate[Join@@("names"/.{namedtensors})]}

Pauli Matrices

pauliZ=SparseArray@{{1,0},{0,-1}};
pauliX=SparseArray@{{0,1},{1,0}};
pauliY=SparseArray@{{0,-I},{I,0}};
pauliI = SparseArray@{{1,0},{0,1}};

id[left_,right_,len_]:={"tensor"->IdentityMatrix[len],"names"->{left,right}}

New Routines

Contractions (Inner Products)

To Contract two tensors, the most natural way is to use “Dot” -- one puts the indices that you want contracted at the beginning and end of the two tensor.

Mb={"tensor"->(M=SparseArray[{{1,1,1}->a,{2,2,1}->b,{2,3,2}->c,{3,4,2}->d},{3,4,2}]),"names"->{"α","β","σ"}}

Mc=ReOrder[Mb,{"σ","α","β"}]

Md=ReOrder[Mb,{"α","β","σ"}]

Dot["tensor"/.Md,"tensor"/.Mc]

Again we can encapsulate this

Contract[namedtensor1_List,namedtensor2_List,indextocontract_]:=
Module[{t1, 
t2,
names1="names"/.namedtensor1,
names2="names"/.namedtensor2,
newnames1,newnames2,newnamedtensor1,newnamedtensor2,newt,newnames},
(*Get new ordering of names*)
newnames1=RotateLeft[names1,Flatten@Position[names1,indextocontract]];
newnames2=RotateLeft[names2,Flatten@Position[names2,indextocontract]-1];
(*Reorder tensor indices*)
newnamedtensor1=ReOrder[namedtensor1,newnames1];
newnamedtensor2=ReOrder[namedtensor2,newnames2];
(*Extract Tensors*)
t1="tensor"/.newnamedtensor1;
t2="tensor"/.newnamedtensor2;
(*Do Dot Product *)
newt=Dot[t1,t2];
newnames=Join[newnames1[[;;-2]],newnames2[[2;;]]];
{"tensor"->newt,"names"->newnames}
]

R1=Contract[Mb,Mb/.{"α"->"α`","β"->"β`"},"σ"]

Another thing we might want to do is to trace over two indices --  which can be done with “TensorContract”

TensorContract["tensor"/.R1,{1,3}]

We can overload the same command for doing this -- if you give one tensor to “Contract” it contracts indices, if you give two, it does a dot product:

Contract[namedtensor_List,indicestocontract_List]:=
Module[{tensor="tensor"/.namedtensor,
names="names"/.namedtensor,
positions,newtensor,keep,newnames},
(* Find out which indices we need to contract *)
positions=Flatten[Position[names,#]&/@indicestocontract];
(* Do contraction *)
newtensor=TensorContract[tensor,positions];
(* Find out names of indices we are keeping *)
keep=Complement[Range[Length[names]],positions];
newnames=names[[keep]];
(* Reassemble *)
{"tensor"->newtensor,"names"->newnames}
]

Usage:

Contract[R1,{"α","α`"}]

Fancier contractions

Another way to do things (which is the notation of ITensor for example) is to use “Einstein” notation -- and simply contract repeated indices.  We first make the single-tensor version

autoContract[namedtensor_List]:=
Module[{tensor="tensor"/.namedtensor,newtensor,
names="names"/.namedtensor,signature,repeats,newindices},
signature=Last@Reap[MapIndexed[Sow[#2,#1]&,names],_,Rule];
repeats=Cases[signature[[All,2]],(locs:{_,_}):>Flatten[locs]];
newtensor=TensorContract[tensor,repeats];
newindices=Delete[names,Transpose[{Flatten[repeats]}]];
{"tensor"-> newtensor,"names"->newindices}
]

R1

autoContract[R1/.{"α`"->"α"}]

Now we do two tensors

autoContract[namedtensor1_List,namedtensor2_List]:=
Module[{names1="names"/.namedtensor1,
names2="names"/.namedtensor2,commonnames,left,right,leftnames,rightnames,
ltens,rtens,newtensor,newnames,pos1,pos2},
commonnames=Intersection[names1,names2];
left=rightCombineIndices[namedtensor1,commonnames];
right=leftCombineIndices[namedtensor2,commonnames];
ltens="tensor"/.left;
rtens="tensor"/.right;
newtensor=Dot[ltens,rtens];
leftnames="names"/.left;
rightnames="names"/.right;
newnames=Join[leftnames[[;;-2]],rightnames[[2;;]]];
{"tensor"->newtensor,"names"->newnames}]

|

autoContract[R1,R1/."α"->"γ"]

autoContract[R1, R1]

Now more

autoContract[namedtensor1_List,namedtensor2_List,moretensors__List]:=
autoContract[autoContract[namedtensor1,namedtensor2],moretensors]

To take the trace of the product of three matrices, we could do

sig1={"tensor"->pauliX,"names"->{"a","b"}};
sig2={"tensor"->pauliX,"names"->{"b","c"}};
sig3={"tensor"->IdentityMatrix[2],"names"->{"c","a"}};

autoContract[sig1,sig2,sig3]

Another example

autoContract[Mb ,Mb]

Note:  If the indices are all distinct, autoContract does the same thing as Combine

autoContract[id["a","a'",2],id["b","b'",2],id["c","c'",2]]

Adding Tensors

addTensor[namedtensor1_List,namedtensor2_List]:=
Module[{tens1="tensor"/.namedtensor1,
tens2="tensor"/.namedtensor2,
names1="names"/.namedtensor1,
names2="names"/.namedtensor2,ordnamedtens2},
If[names1==names2,
Return[{"tensor"->tens1+tens2,"names"->names1}]];
If[Sort[names1]==Sort[names2],
ordnamedtens2=ReOrder[namedtensor2,names1];
Return[{"tensor"->tens1+("tensor"/.ordnamedtens2),"names"->names1}]];
Failure["Shape Mismatch",Association["MessageTemplate"->
"Trying to add a tensor with indices `ind1` to one with indices `ind2`",
"MessageParameters"->Association["ind1"->names1,"ind2"->names2]]]
]

addTensor[R1,R1]

Applying operations to tensors

It is useful to be able to apply an operation (such as complex conjugation) to a tensor in one of these constructions

ApplytoTensor[namedtensor_List,op_]:=namedtensor /.("tensor"->t_):>("tensor"->op[t])

ctens={"tensor"->SparseArray[{{1,1,1}->I,{2,2,1}->2-I,{2,3,2}->0.3,{3,4,2}->0.4},{3,4,2}],"names"->{"a","b","c"}}

cctens=ApplytoTensor[ctens,Conjugate]

ArrayRules["tensor"/.cctens ]

ArrayRules["tensor"/.multiply[ctens,2]]

Multiplying a tensor by a number

multiply[namedtensor_List,number_]:=namedtensor /.("tensor"->t_):>("tensor"->number t)

Exercises

In the lecture 1 and Homework 1, you work with a Matrix Product state which is the product of tensors of the form

Md2={"tensor"->SparseArray[{{1,1,1}->1,{1,2,2}->ϵ,{2,1,1}->ϵ,{2,2,2}->1}],"names"->{"s'","s","σ"}}

MatrixForm["tensor"/.Md2]

Generate the transfer matrix EI, as a rank 4 tensor.  Hint -- you first need to make M* with the correct indices, and then use autoContract

Combine the indices of EI to express it as a 4x4 matrix

Use MatrixForm to Display it

Find the right eigenvectors and eigenvalues

Take the right eigenvector corresponding to the largest eigenvalue, and re-express it as a 2x2 matrix in the s-t space.