Implement `evolve` of an `MPS` with an `MPO`

[jofrevalles]

Summary

This PR introduces the evolve function, which evolves a Matrix Product State (MPS) using a Matrix Product Operator (MPO). The idea from @starsfordummies is that we contract the tensors of each corresponding site, and we then create the new λ with a kron product with the old λ (which come from the MPS) and the identity (comes from the MPO). Right now, this function only works for mps in the canonical form.

Example

Here we show how can we use this function:

julia> using Qrochet

julia> ψ = rand(Chain, Open, State; n=8, χ=10); canonize!(ψ)
MPS (inputs=0, outputs=8)

julia> mpo = rand(Chain, Open, Operator; n=8, χ=10)

MPO (inputs=8, outputs=8)

julia> Qrochet.@reindex! inputs(mpo) => outputs(ψ)

MPS (inputs=0, outputs=8)

julia> ψ = Qrochet.evolve(ψ, mpo)
MPS (inputs=0, outputs=8)

julia> size.(tensors(ψ))
15-element Vector{Tuple{Int64, Vararg{Int64}}}:
 (2, 8)
 (2, 8, 40)
 (2, 40, 80)
 (2, 80, 100)
 (2, 100, 80)
 (2, 80, 40)
 (2, 40, 8)
 (2, 8)
 (8,)
 (40,)
 (80,)
 (100,)
 (80,)
 (40,)
 (8,)

[jofrevalles]

This is not the optimal way to evolve a MPS through a MPO, but it's good enough. Thanks for your PR!

What is more optimal? Here we just contract and update the bond dimensions.

[mofeing]

What is more optimal? Here we just contract and update the bond dimensions.

Well, you can do MPS compression algorithm with the partial traces but inserting the MPO in between the original MPS and the target MPS.

But just like we did with evolve!, it's good to have this method; we just need to keep in mind but in the future we might need to extend the function to support this other method.

[jofrevalles]

@mofeing Now I managed to make this function in place, but if we followed my suggestion in https://github.com/bsc-quantic/Tenet.jl/issues/147, these could be much more clean.

[mofeing]

Great! Then let's fix bsc-quantic/Tenet.jl#147 first

[jofrevalles]

@mofeing Now this works in place and the code is much better with Tenet.@unsafe_region.

julia> mps = rand(Chain, Open, State; n=8, χ=10)
MPS (inputs=0, outputs=8)

julia> mpo = rand(Chain, Open, Operator; n=8, χ=10)
MPO (inputs=8, outputs=8)

julia> canonize!(mps)
MPS (inputs=0, outputs=8)

julia> Qrochet.@reindex! inputs(mpo) => outputs(mps)
MPS (inputs=0, outputs=8)

julia> evolve!(mps, mpo)
MPS (inputs=0, outputs=8)

julia> size.(tensors(mps))
15-element Vector{Tuple{Int64, Vararg{Int64}}}:
 (8,)
 (40,)
 (80,)
 (100,)
 (80,)
 (40,)
 (8,)
 (2, 8)
 (2, 8, 40)
 (2, 40, 80)
 (2, 80, 100)
 (2, 100, 80)
 (2, 80, 40)
 (2, 40, 8)
 (2, 8)

[mofeing]

Wait, why are there vectors of size 80? Also the sizes are too big for a 8-site MPS. The maximum theoretical bond dimension for a physical dim of 2 should be $2^4 = 16$.

[jofrevalles]

Wait, why are there vectors of size 80? Also the sizes are too big for a 8-site MPS. The maximum theoretical bond dimension for a physical dim of 2 should be 24=16.

The bond dimensions go up to 100 since for both mps and mpo, the maximum bond dimension is 10, and since we have to contract both bond dimensions we get 10^2.