Lazy ITensor contraction

[mtfishman]

This is a reminder that it would be great to have lazy ITensor contraction support and broadcasting working with contractions.

This could be done through the broadcast syntax, for example:

C .= a .* A .* B   # mul!(C, A, B, a)
C .= a .* A .* B .+ b .* C  # mul!(C, A, B, a, b)
D .= A .* B .* C  # lazily contract A, B, and C, optimizing contraction order
D .= A .+ B .+ C  # Fuse into a single call

[mtfishman]

Some of these were addressed by #251. It would still be nice to have lazy contraction of arbitrary numbers of ITensors, with contraction order optimization.

I think that we could even include contraction order optimization with regular tensor contraction, like A*B*C, since:

julia> i = Index(2)
(dim=2|id=53)

julia> A = randomITensor(i,i')
ITensor ord=2 (dim=2|id=53) (dim=2|id=53)' 
Dense{Float64,Array{Float64,1}}

julia> B = randomITensor(i',i'')
ITensor ord=2 (dim=2|id=53)' (dim=2|id=53)'' 
Dense{Float64,Array{Float64,1}}

julia> C = randomITensor(i'', i)
ITensor ord=2 (dim=2|id=53)'' (dim=2|id=53) 
Dense{Float64,Array{Float64,1}}

julia> @which A*B*C
*(a, b, c, xs...) in Base at operators.jl:529

so it looks like we can overload *(::ITensor, ::ITensor, ::ITensor, ::ITensor...) and do contraction order optimization (and other optimizations, for example ones involving delta tensors).

[emstoudenmire]

Yes, let's definitely look into this. It could be a really killer feature of ITensor. I'm really glad we might be able to do this just with * in addition to broadcasted operations.

One issue to think through is how we can deal with views versus copies. Copying the tensors is safe but costly; not copying is more efficient but trickier for ensuring correctness.

So I imagine what you're thinking of here is a sort of semi-eager mode, where we are lazy about evaluating the separate contractions within R = ABC*D but as soon as we have all of those tensors, we eagerly compute the result R before going on to the next line of code. Correct?

[mtfishman]

In both cases (the broadcasting and non-broadcasting versions), I was thinking it would work semi-eagerly as you say.

In the case of broadcasting, Julia rewrites a call like R = A .* B .* C into a Broadcast.Broadcasted object that stores views of A, B and C (and the operations between them). Then you can control the final operation that is done with that object which gets stored in R, so the idea would be that it would do the contraction of A, B, and C into normal ITensor that gets stored in R. In theory we could output a fancier lazy ITensor storage that doesn't actually do the contraction yet, but I don't see the need right now.

In the case of R = A*B*C, it appears to be a built-in Julia feature that this gets rewritten into a call R = *(A,B,C), which means that we can overload *(A::ITensor, B::ITensor, C::ITensor). Again, in theory we could overload that function to output a lazy ITensor that stores those tensors without doing the actually operation, but it seems simpler for now to have that actually perform the contraction with possible order optimization.

I think both of those should cover many cases of lazy operations that we are interested in (we can just tell people that if they want operations treated lazily, put them on the same line and possibly use the dot syntax for more complicated expressions).