Improving infinite mpo matrix

[LHerviou]

As described, some redefinitions and improvements on the form of InfiniteMPOMatrix, and modifications that follow.
The core idea:

• Arbitrary filled InfiniteMPOMatrix are allowed
• Instead of waiting for nrange(mpo) to evaluate the terms in the local MPO, we do not delay evaluation. T
• Lines and columns in an InfiniteMPOMatrix have well defined indices, shared for all elements.

All of that normally significantly simplifies the manipulation and the optimization of the InfiniteMPOMatrix and the corresponding InfiniteMPO
I did not have time to fully debug, so I will keep updating in the next few days, and confirm when one can merge.

[LHerviou]

All tests work on my branch without any measurable slowdown.

Probably a lot of simplification/optimization can be done. Up to you if you integrate it right now, or if we wait a bit til I can have an overview of possible changes.

[mtfishman]

@LHerviou I'm fine waiting to merge this if you think there are simplifications or optimizations to make once you have the time, unless @JanReimers needs this for something and therefore it's useful to have merged.

[JanReimers]

I have a branch with

• InfiniteMPO compression outputting InfiniteCanonicalMPO
• VUMPS suppert for InfiniteMPO
• Simplification of VUMPS generic code
• VUMPS suppert for InfiniteCanonicalMPO (which may require some discussion)
  I was volunteering to let Loic merge first and then I can deal with any resulting conflicts.

[LHerviou]

On my computer, this runs nicely.

To tackle some missing methods in NDTensors, I explicitely defined an addition formula for EmptyITensors. I would be interested in someone having a look at it for type stability. This touches part of the Pkg I am not familiar with.

[LHerviou]

Sounds perfect. Does it allow also to remove the extension of datatype I had to do? Thanks a lot, Mat.

…

On Tue, Jun 20, 2023 at 8:23 PM Matt Fishman ***@***.***> wrote: ***@***.**** commented on this pull request. ------------------------------ In src/infinitempomatrix.jl <#77 (comment)> : > + new_tags = get(kwargs, :left_tags, tags(old_ind)) + top, middle, bottom = build_three_projectors_from_index(old_ind; tags=new_tags) + vector = fill(op("Zero", local_sit), 3, 1) + else + new_tags = get(kwargs, :right_tags, tags(old_ind)) + top, middle, bottom = build_three_projectors_from_index(old_ind; tags=new_tags) + vector = fill(op("Zero", local_sit), 1, 3) + end + for (idx, proj) in enumerate([top, middle, bottom]) + vector[idx] = proj * tensor + end + return vector +end + +#Matrix multiplication for Matrix of ITensors +function Base.:*(A::Matrix{ITensor}, B::Matrix{ITensor}) @LHerviou <https://github.com/LHerviou> I have a simple fix for this issue which makes sure the contraction code outputs the correct types, so then there is no issue in subsequent operations like addition. I'll make a PR to NDTensors/ITensors directly. I can also add the definition Base.zero(T::ITensor) = false * T so then with those two fixes Matrix{ITensor} * Matrix{ITensor} will "just work". — Reply to this email directly, view it on GitHub <#77 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AXMRNR4CA3F7F3G2L2OG343XMHTBHANCNFSM6AAAAAAZIEX4ZI> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[mtfishman]

I don't seem to need to overload NDTensors.datatype for Matrix{ITensor} * Matrix{ITensor} to work. Is it possible another operation needed that?

[LHerviou]

I will check again tomorrow. It is possible it came from something else.

…

On Tue, Jun 20, 2023 at 11:23 PM Matt Fishman ***@***.***> wrote: I don't seem to need to overload NDTensors.datatype for Matrix{ITensor} * Matrix{ITensor} to work. Is it possible another operation needed that? — Reply to this email directly, view it on GitHub <#77 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AXMRNR653Z3QAD5V6T675KTXMIIEDANCNFSM6AAAAAAZIEX4ZI> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[mtfishman]

@LHerviou the issues we discussed have been implemented in the latest version of ITensors.jl (ITensor/ITensors.jl#1135).

[LHerviou]

Apparently, there is still one definition missing with how you implemented zero.

julia> T = ITensor()
ITensor ord=0
NDTensors.EmptyStorage{NDTensors.EmptyNumber, NDTensors.Dense{NDTensors.EmptyNumber, Vector{NDTensors.EmptyNumber}}}

julia> zero(T)
ERROR: MethodError: no method matching fill!(::NDTensors.NoData, ::NDTensors.EmptyNumber)

Adding "Base.fill!(::NDTensors.NoData, ::Any) = 0 " is enough to solve the issue.

[mtfishman]

@LHerviou what do you think of changing the name of the type from InfiniteMPOMatrix to InfiniteBlockMPO? I think that is a bit more descriptive, and captures that the MPO is composed of block matrices.

I do not mind that change, as long as I don't have to change it again after that.

Ideally we would come up with a name that works well that we won't have to change going forward. Also we could do that name change in a future PR.

[LHerviou]

I think InfiniteBlockMPO is a good name. It is a MPO by block in the same sense block matrices exist.
So I have no objections.

I have no idea how many people (beyond me and some of my collaborators) used the InfiniteMPOMatrix anyway.

[mtfishman]

@LHerviou this looks good to me now. Do you think this is ready to merge?

[LHerviou]

Hi I believe so.

…

On Tue, 5 Sept 2023, 19:28 Matt Fishman, ***@***.***> wrote: @LHerviou <https://github.com/LHerviou> this looks good to me now. Do you think this is ready to merge? — Reply to this email directly, view it on GitHub <#77 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AXMRNRYSBNUSG3QI2ZCVKDDXY5OL7ANCNFSM6AAAAAAZIEX4ZI> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[LHerviou]

For the future: I guess I am letting @JanReimers finally do his own updates, then I will add all the iDMRG and so on using his version of compression.

[mtfishman]

Thanks again!

[mtfishman]

@LHerviou not sure what's going on but it looks like this PR broke some tests on main (https://github.com/ITensor/ITensorInfiniteMPS.jl/commits/main), could you take a look?

[LHerviou]

I will have a look either tonight or tomorrow.
The errors seem to imply something is missing in ITensors.NDTensors, but it might be something else

[mtfishman]

Ah we have been changing some internals of NDTensors, and I guess this PR relied on some of those internals, so I could imagine there might be an issue there.

We are in the process of refactoring and simplifying NDTensors so hopefully it will mean some of the NDTensors overloads made in this PR can just be removed in the near future.

[LHerviou]

@mtfishman Ok, I found the source of the bug.

l269 and l272 of [tensor_operations/tensor_algebra.jl](https://github.com/ITensor/ITensors.jl/blob/main/src/tensor_operations/tensor_algebra.jl),
you changed the way you define the directsum_projectors! from

 D1[...] = true

to

D1[...] = one(eltype(D1))

In the end, the problem is that one(EmptyNumber) is not defined. I am testing some fixes right now.

[mtfishman]

Ok interesting, thanks for looking into it. Does that mean an EmptyStorage is getting direct-summed somewhere?

[LHerviou]

Yes, when I am fusing some indices.
To be more precise: when I am moving from a big MPOMatrix to the InfiniteMPO

[LHerviou]

For info:
NDTensors.convert(::Type{NDTensors.EmptyNumber}, x::Float64) = NDTensors.EmptyNumber()
fixes the problem. Should I update ITensorInfiniteMPS or do you want to wait to fix NDTensors.

I am not sure it is a good workaround though, as the projectors are not projectors anymore, and could be used at some point later.

[mtfishman]

I'm not sure that is something we want to define, since conversion shouldn't lose information (i.e. Float64 -> ComplexF64 conversion is fine since it widens the type, but ComplexF64 -> Float64 is only ok if the imaginary part is zero).

I wonder if one(::EmptyNumber) = true would fix it? I think that makes some sense as a definition anyway, since EmptyNumber is kind of a universal zero value and true is the smallest one value in Julia.

[mtfishman]

We are replacing the EmptyStorage type with a better and simpler design in ITensor/ITensors.jl#1161 so hopefully things like this will be easier to deal with after that is merged.