MPO-dense matrix contraction

[RaimelMedina]

Hi all,

I've been trying find a good way to contract an MPO (constructed by tensorizing and decomposing some given matrix) with a dense tensor. I find, however, that the returned contraction path essentially takes the MPO to its dense form and then performs the contraction with the dense tensor. I wonder if I'm doing something wrong when specifying the contraction path or if in fact that is the best way to carry out the contraction. I would appreciate any input you guys can give me. Thanks!

using OMEinsum, OMEinsumContractionOrders

size_dict = Dict('a'=> 4, 'g'=> 4, 'm'=> 16,
  'b'=> 4, 'h'=> 4, 'n'=> 256, 'c'=> 4, i'=> 4,
  'o'=> 1228, 'd'=> 4, 'j'=> 4, 'p'=> 896, 'e'=> 4,
  'k'=> 7, 'q'=> 32, 'f'=> 4, 'l'=> 8, 'x'=> 65536)
  
code = ein'agm,mbhn,ncio,odjp,pekq,qfl,xabcdef->xghijkl'

optcode = optimize_code(code, size_dict, ExactTreewidth())

The layout of the MPO cores here is (left_bond, in_dim, out_dim, right_bond) and the input tensor is written as (b, in_1, ..., in_N)
Doing this I got:

agbhciekfldj, xabcdef -> xghijkl
├─ agbhcio, ekflodj -> agbhciekfldj
│  ├─ agbhn, ncio -> agbhcio
│  │  ├─ agm, mbhn -> agbhn
│  │  │  ├─ agm
│  │  │  └─ mbhn
│  │  └─ ncio
│  └─ pekfl, odjp -> ekflodj
│     ├─ pekq, qfl -> pekfl
│     │  ├─ pekq
│     │  └─ qfl
│     └─ odjp
└─ xabcdef

So essentially, first reconstruct the dense matrix from the MPO cores and then contract :|

[GiggleLiu]

Try use a different backend, e.g. TreeSA(). If you get consistent result, then it is optimal. BTW, you might want to compare their complexities with contraction_complexity(code).

[RaimelMedina]

Thanks for the reply :) I tried that and indeed, I get a consistent result (plus a warning)

julia> optcode = optimize_code(code, size_dict, TreeSA())
┌ Warning: target space complexity not found, got: 29.807354922057606, with time complexity 41.81454635234436, read-write complexity 30.319728342837152.
└ @ OMEinsumContractionOrders ~/.julia/packages/OMEinsumContractionOrders/2wuNi/src/treesa.jl:229
SlicedEinsum{Char, DynamicNestedEinsum{Char}}(Char[], ekfldjcibhag, xabcdef -> xghijkl
├─ ekflodj, ciobhag -> ekfldjcibhag
│  ├─ pekfl, odjp -> ekflodj
│  │  ├─ pekq, qfl -> pekfl
│  │  │  ├─ pekq
│  │  │  └─ qfl
│  │  └─ odjp
│  └─ ncio, bhnag -> ciobhag
│     ├─ ncio
│     └─ mbhn, agm -> bhnag
│        ├─ mbhn
│        └─ agm
└─ xabcdef
)

And this is the contraction complexity (same for both)

julia> contraction_complexity(optcode, size_dict)
Time complexity: 2^41.81454635234436
Space complexity: 2^29.807354922057606
Read-write complexity: 2^30.319728342837152

I wasn't expecting this, but it is good to know that it is not possible to do better then :) Thanks once again