Fix Multiline (MultilineMPS/MultilineMPO) correctness and API inconsistencies

docs/src/changelog.md

@@ -38,6 +38,7 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
   by `MPSKit.default_allocator`, instead of leaving them to the garbage collector
   (two-site DMRG: -64% allocations, -57% GC time, -23% wall time).
   Disable with `MPSKit.Defaults.set_buffering!(false)`. ([#467](https://github.com/QuantumKitHub/MPSKit.jl/pull/467))
+- Custom `show`/`summary` for `MultilineMPS`/`MultilineMPO`. Each row is now rendered via each row's own display, and row shifting is shown explicitly for `MultilineMPO`.
 
 ### Changed
 
@@ -70,11 +71,23 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 - `correlator` now throws an `ArgumentError` when the sites are not ordered as `i < j`.
   Previously such a call only logged an `@error` and then continued into a contraction that is
   not the requested correlator. ([#489](https://github.com/QuantumKitHub/MPSKit.jl/pull/489))
+- `Multiline` (and therefore `MultilineMPS`/`MultilineMPO`) now consistently treats
+  `length`/`eltype`/`iterate`/`m[i]` as referring to the individual lines it stores
+  (`length(m) == nrows`), while `size`/`axes`/`eachindex` refer to the `(nrows, ncols)` lattice
+  shape.
+- `MultilineMPO` construction is now restricted to `InfiniteMPO` lines. Constructing one from
+  Hamiltonians or finite MPOs now throws a `MethodError` at the construction site instead of producing an object that fails
+  later. The `AbstractMatrix` constructor that silently built finite-line `MultilineMPO`s was
+  removed.
 
 ### Deprecated
 
 ### Removed
 
+- `expectation_value(::MultilineMPS, ::MultilineMPO, envs...)` fallback method, which silently
+  computed a meaningless value (`prod` instead of `sum`, no row shift, `envs` ignored) for any
+  `MultilineMPO` line type not covered by the guarded method. Most notably this prevents a fallback for `InfiniteMPOHamiltonian`, a legal but never-meaningful `Multiline` line type.
+
 ### Fixed
 
 - `isfinite(::WindowMPOHamiltonian)` was undefined. ([#489](https://github.com/QuantumKitHub/MPSKit.jl/pull/489))
@@ -96,6 +109,19 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
   and the right virtual leg of `mpo[end]` and contracted the two — which at length 1 is the *same*
   tensor, so it returned `O * O` on twice the physical space instead of `O`.
   ([#484](https://github.com/QuantumKitHub/MPSKit.jl/pull/484))
+- `MultilineMPO * MultilineMPS` and `MultilineMPO * MultilineMPO` (`*`) never worked at all
+  (`MethodError`, from `map(*, zip(...))` calling `*` on a 2-tuple). Fixing the immediate error
+  still left a shape bug (`O * ψ` produced an `nrows * ncols`-line object with `#undef` entries
+  past `nrows`) and a row-shift bug (`O * ψ` had fidelity 0 with the expected state on every
+  row). `*` now correctly applies the convention that row `i` of a `MultilineMPO` maps row `i`
+  of the state onto row `i + 1`.
+- `isfinite(::MultilineMPO)` threw (`isfinite(typeof(m))` had no matching type-level method for
+  `Multiline`).
+- `changebonds(::MultilineMPO, ::SvdCut)` threw (`convert(MultilineMPS, ::MultilineMPO)` has no
+  method).
+- `axes(m::Multiline, i)` threw for `i > 2`, but now returns `Base.OneTo(1)`, matching Base's own
+  out-of-range convention (already the case for `size(m, i)`).
+- `spacetype`/`sectortype`/`storagetype` on a `Multiline` instance were undefined. Only the type-level methods existed.
 
 ### Performance
 

docs/src/man/operators.md

@@ -274,3 +274,22 @@ a collection (direct sum) of spaces, one for each row/column.
 ```@example operators
 left_virtualspace(H_ising, 1), right_virtualspace(H_ising, 1), physicalspace(H_ising, 1)
 ```
+
+## MultilineMPO
+
+Much like a [`MultilineMPS`](@ref) is a repeating set of `InfiniteMPS` lines, a [`MultilineMPO`](@ref) is a repeating set of `InfiniteMPO` lines.
+This is typically the row-to-row or column-to-column transfer matrix of a 2D classical partition function, or a boundary MPO in the context of PEPS.
+See the `MultilineMPS` section of the [states](@ref um_states) page for the row-shift convention relating a `MultilineMPO` to the `MultilineMPS` it acts on: row `i` of the operator maps row `i` of the state onto row `i + 1`.
+
+```@example operators
+mpo_multi = MultilineMPO([mpo, mpo])
+```
+
+Only `InfiniteMPO` lines are supported; neither Hamiltonians nor finite MPOs are allowed currently.
+This is enforced by the constructors rather than by the type itself, so `MultilineMPO([H, H])` for a Hamiltonian `H` throws a `MethodError` right at the construction site, instead of producing an object that only fails once used.
+
+!!! note "`*` between `MultilineMPO`s does not compose the way you might expect"
+    Because of the `+1` row shift, multiplying two `MultilineMPO`s produces an operator that spans two rows of the lattice: `(O1 * O2)` maps row `i` onto row `i + 2`.
+    As a consequence, `(O1 * O2) * ψ != O1 * (O2 * ψ)` in general, for any definition of the product.
+    The two sides are genuinely different operations (one fused two-row step, versus two sequential one-row steps).
+    Treat `O1 * O2` as a "fused double row", not as an operator that can be composed further.

docs/src/man/states.md

@@ -170,5 +170,35 @@ col = 2
 al = state.AL[row, col];
 ```
 
+### The row-shift convention
+
+The row direction and the column direction of a `MultilineMPS` play different roles.
+Within a row, a `MultilineMPS` behaves exactly like the `InfiniteMPS` it repeats: columns are periodic, and `state.AL[row, col]`, `state.AR[row, col]`, `state.C[row, col]` and `state.AC[row, col]` behave exactly as they would for `state[row]::InfiniteMPS`.
+
+The row direction is the direction along which a transfer matrix, represented as a [`MultilineMPO`](@ref), is applied.
+By convention, row `i` of a `MultilineMPO` maps row `i` of the state onto row `i + 1`, so that
+
+```julia
+(O * ψ)[i] == O[i - 1] * ψ[i - 1]
+```
+
+i.e. applying `O` shifts every row up by one.
+This convention is not just a property of `*`: it is baked into `environments`, `expectation_value`, and the derivative machinery as well, all of which pair `ψ[row + 1]` against `O[row]`.
+It is also why `O * ψ` requires `size(O, 1) == size(ψ, 1)`, rather than the row counts of operator and state being independent.
+
+### What is currently supported
+
+Only lines that are themselves infinite are supported: a `MultilineMPS` is always built out of
+[`InfiniteMPS`](@ref) lines, never [`FiniteMPS`](@ref) or [`WindowMPS`](@ref) ones.
+This matches the intended use case of rows/columns in an infinite 2D network.
+The constructors reject finite lines.
+
+### Subtleties
+
+- **`size` vs. iteration:** `size(state)` returns `(nrows, ncols)`, describing the lattice shape.
+Iterating over a `MultilineMPS` (or indexing it with a single integer, `state[i]`) instead yields the individual `InfiniteMPS` *lines*, so `length(state) == nrows`, not `nrows * ncols`.
+Code that wants to operate line-by-line should use `state[i]`/`parent(state)`, while code that wants the lattice shape should use `size`.
+- **Norms and inner products:** `dot`/`norm` sum the contribution of every row, so a `MultilineMPS` whose individual rows are each normalized to 1 does *not* itself have norm 1: `norm(state) == sqrt(nrows)` for `nrows` identical normalized rows.
+
 These objects are also used extensively in the context of [PEPSKit.jl](https://github.com/QuantumKitHub/PEPSKit.jl).
 

src/algorithms/changebonds/randexpand.jl

@@ -56,7 +56,7 @@ function changebonds!(ψ::InfiniteMPS, alg::RandExpand)
 end
 
 function changebonds!(ψ::MultilineMPS, alg::RandExpand)
-    foreach(Base.Fix2(changebonds!, alg), ψ.data)
+    foreach(Base.Fix2(changebonds!, alg), parent(ψ))
     return ψ
 end
 

src/algorithms/changebonds/svdcut.jl

@@ -86,14 +86,14 @@ function changebonds!(mpo::FiniteMPO, alg::SvdCut)
 end
 
 # TODO: this assumes the MPO is infinite, and does weird things for finite MPOs.
-function changebonds(ψ::InfiniteMPO, alg::SvdCut)
-    return convert(InfiniteMPO, changebonds(convert(InfiniteMPS, ψ), alg))
+function changebonds(mpo::InfiniteMPO, alg::SvdCut)
+    return convert(InfiniteMPO, changebonds(convert(InfiniteMPS, mpo), alg))
 end
-function changebonds(ψ::MultilineMPO, alg::SvdCut)
-    return convert(MultilineMPO, changebonds(convert(MultilineMPS, ψ), alg))
+function changebonds(mpo::MultilineMPO, alg::SvdCut)
+    return Multiline(map(Base.Fix2(changebonds, alg), parent(mpo)))
 end
 function changebonds(ψ::MultilineMPS, alg::SvdCut)
-    return Multiline(map(x -> changebonds(x, alg), ψ.data))
+    return Multiline(map(Base.Fix2(changebonds, alg), parent(ψ)))
 end
 function changebonds(ψ::InfiniteMPS, alg::SvdCut)
     copied = copy.(ψ.AL)

src/algorithms/expval.jl

@@ -32,6 +32,13 @@ the operator is a `AbstractTensorMap` that acts on the physical space of a singl
     return value is the total over one unit cell; divide by `length(ψ)` to obtain a
     per-site value.
 
+!!! note "MultilineMPS/MultilineMPO"
+    For a `MultilineMPS`/`MultilineMPO{<:InfiniteMPO}` pair, `O` is applied with the usual
+    row-shift convention (row `i` of `O` pairs `ψ[i]` against `ψ[i + 1]`).
+    The returned value is the fixed-point contraction, not a normalized `⟨ψ|O|ψ⟩ / ⟨ψ|ψ⟩` in the sense of
+    the single-line methods.
+    `MultilineMPO` lines of `InfiniteMPOHamiltonian` (or any non-`InfiniteMPO` line type) are not supported.
+
 # Examples
 
 ```jldoctest
@@ -180,8 +187,13 @@ end
 function expectation_value(ψ::FiniteQP, mpo::FiniteMPO)
     return expectation_value(convert(FiniteMPS, ψ), mpo)
 end
-function expectation_value(ψ::InfiniteMPS, mpo::InfiniteMPO, envs...)
-    return expectation_value(convert(MultilineMPS, ψ), convert(MultilineMPO, mpo), envs...)
+function expectation_value(ψ::InfiniteMPS, mpo::InfiniteMPO)
+    return expectation_value(convert(MultilineMPS, ψ), convert(MultilineMPO, mpo))
+end
+function expectation_value(ψ::InfiniteMPS, mpo::InfiniteMPO, envs::AbstractMPSEnvironments)
+    return expectation_value(
+        convert(MultilineMPS, ψ), convert(MultilineMPO, mpo), convert(MultilineEnvironments, envs)
+    )
 end
 function expectation_value(
         ψ::MultilineMPS, O::MultilineMPO{<:InfiniteMPO},
@@ -193,10 +205,6 @@ function expectation_value(
         return contract_mpo_expval(ψ.AC[i, j], GL, O[i, j], GR, ψ.AC[i + 1, j])
     end
 end
-function expectation_value(ψ::MultilineMPS, mpo::MultilineMPO, envs...)
-    # TODO: fix environments
-    return prod(x -> expectation_value(x...), zip(parent(ψ), parent(mpo)))
-end
 # fallback
 function expectation_value(ψ::AbstractMPS, mpo::AbstractMPO, envs...)
     return dot(ψ, mpo, ψ) / dot(ψ, ψ)

src/algorithms/grassmann.jl

@@ -11,7 +11,7 @@ The module exports nothing, and all references to it should be qualified, e.g.
 module GrassmannMPS
 
 using ..MPSKit
-using ..MPSKit: AbstractMPSEnvironments, InfiniteEnvironments, MultilineEnvironments,
+using ..MPSKit: AbstractMPSEnvironments, InfiniteEnvironments, MultilineEnvironments, site_type,
     AC_projection, recalculate!, TimerOutput, DISABLED_TIMER, @timeit, default_allocator
 using TensorOperations: AbstractBackend, DefaultBackend
 using TensorKit
@@ -216,12 +216,12 @@ function fg(
         backend::AbstractBackend = DefaultBackend(),
         scheduler::Scheduler = MPSKit.Defaults.scheduler[],
     )
-    @assert length(state) == 1 "not implemented"
+    @assert size(state, 1) == 1 "not implemented"
     @timeit timeroutput "envs (parallel)" recalculate!(envs, state, operator, state; timeroutput)
     f = @timeit timeroutput "expval" expectation_value(state, operator, envs)
     isapprox(imag(f), 0; atol = eps(abs(f))^(3 / 4)) || @warn "MPO might not be Hermitian: $f"
 
-    A = Core.Compiler.return_type(Grassmann.project, Tuple{eltype(state), eltype(state)})
+    A = Core.Compiler.return_type(Grassmann.project, Tuple{site_type(state), site_type(state)})
     gs = Matrix{A}(undef, size(state))
     allocator = default_allocator(state, scheduler)
     @timeit timeroutput "gradient" tforeach(eachindex(state); scheduler) do i

src/environments/multiline_envs.jl

@@ -98,3 +98,7 @@ function transfer_rightenv!(
     )
     return transfer_rightenv!(envs, below, O, above, site)
 end
+
+# converters
+Base.convert(::Type{MultilineEnvironments}, envs::InfiniteEnvironments) = Multiline([envs])
+Base.convert(::Type{InfiniteEnvironments}, envs::MultilineEnvironments) = only(envs)

src/operators/multilinempo.jl

@@ -5,46 +5,49 @@
 
 Type that represents multiple lines of `MPO` objects.
 
+Only `InfiniteMPO` lines are supported currently.
+Hamiltonians and finite MPOs are rejected by the constructors.
+
 # Constructors
 
-    MultilineMPO(mpos::AbstractVector{<:Union{SparseMPO, DenseMPO}})
-    MultilineMPO(Os::AbstractMatrix{<:MPOTensor})
+    MultilineMPO(mpos::AbstractVector{<:InfiniteMPO})
+    MultilineMPO(Os::PeriodicMatrix{<:MPOTensor})
+    MultilineMPO(t::MPOTensor)
 
 # See also
 
 [`Multiline`](@ref), [`AbstractMPO`](@ref)
 """
 const MultilineMPO = Multiline{<:AbstractMPO}
 
-function MultilineMPO(Os::AbstractMatrix)
-    return MultilineMPO(map(FiniteMPO, eachrow(Os)))
-end
 function MultilineMPO(Os::PeriodicMatrix)
     return MultilineMPO(map(InfiniteMPO, eachrow(Os)))
 end
-MultilineMPO(mpos::AbstractVector{<:AbstractMPO}) = Multiline(mpos)
+MultilineMPO(mpos::AbstractVector{<:InfiniteMPO}) = Multiline(mpos)
 MultilineMPO(t::MPOTensor) = MultilineMPO(PeriodicMatrix(fill(t, 1, 1)))
 
 # allow indexing with two indices
-Base.getindex(t::MultilineMPO, ::Colon, j::Int) = Base.getindex.(t.data, j)
+Base.getindex(t::MultilineMPO, ::Colon, j::Int) = Base.getindex.(parent(t), j)
 Base.getindex(t::MultilineMPO, i::Int, j) = Base.getindex(t[i], j)
 Base.getindex(t::MultilineMPO, I::CartesianIndex{2}) = t[I.I...]
 
 # converters
-Base.convert(::Type{MultilineMPO}, t::AbstractMPO) = Multiline([t])
+Base.convert(::Type{MultilineMPO}, t::InfiniteMPO) = Multiline([t])
 Base.convert(::Type{DenseMPO}, t::MultilineMPO{<:DenseMPO}) = only(t)
 Base.convert(::Type{SparseMPO}, t::MultilineMPO{<:SparseMPO}) = only(t)
-Base.convert(::Type{FiniteMPO}, t::MultilineMPO{<:FiniteMPO}) = only(t)
 Base.convert(::Type{InfiniteMPO}, t::MultilineMPO{<:InfiniteMPO}) = only(t)
 
 function Base.:*(mpo::MultilineMPO, st::MultilineMPS)
     size(mpo) == size(st) || throw(ArgumentError("dimension mismatch"))
-    return Multiline(map(*, zip(mpo, st)))
+    # row i of the operator maps row i of the state onto row i + 1
+    return Multiline(circshift(map(*, parent(mpo), parent(st)), 1))
 end
 
 function Base.:*(mpo1::MultilineMPO, mpo2::MultilineMPO)
     size(mpo1) == size(mpo2) || throw(ArgumentError("dimension mismatch"))
-    return Multiline(map(*, zip(mpo1, mpo2)))
+    # `mpo2[i]` maps row i onto row i + 1, where `mpo1[i + 1]` picks up: the resulting
+    # operator spans two rows of the lattice and maps row i onto row i + 2
+    return Multiline(map(*, circshift(parent(mpo1), -1), parent(mpo2)))
 end
 
 for f_space in (:physicalspace, :left_virtualspace, :right_virtualspace)

src/states/multilinemps.jl

@@ -27,6 +27,11 @@ Type that represents multiple lines of [`InfiniteMPS`](@ref) objects.
 - `AC`: center-gauged MPS tensors
 - `C`: gauge (bond) tensors
 
+Only `InfiniteMPS` lines are supported.
+Note that `length`, `eltype` and iteration refer to the lines (so e.g. `length(ψ) == nrows`),
+while `size` refers to the `(nrows, ncols)` lattice shape.
+See [`Multiline`](@ref) for details.
+
 # See also
 
 [`Multiline`](@ref)
@@ -101,11 +106,8 @@ TensorKit.normalize!(a::MultilineMPS) = (normalize!.(parent(a)); return a)
 
 Base.convert(::Type{MultilineMPS}, st::InfiniteMPS) = Multiline([st])
 Base.convert(::Type{InfiniteMPS}, st::MultilineMPS) = only(st)
-Base.eltype(t::MultilineMPS) = eltype(t[1])
 Base.copy!(ψ::MultilineMPS, ϕ::MultilineMPS) = (copy!.(parent(ψ), parent(ϕ)); ψ)
 
-Base.isfinite(::Type{<:MultilineMPS}) = false
-
 for f_space in (:physicalspace, :left_virtualspace, :right_virtualspace)
     @eval $f_space(t::MultilineMPS, i::Int, j::Int) = $f_space(t[i], j)
     @eval $f_space(t::MultilineMPS, I::CartesianIndex{2}) = $f_space(t, Tuple(I)...)

src/states/orthoview.jl

@@ -277,12 +277,6 @@ Base.size(psi::CView{<:AbstractFiniteMPS}) = (length(psi.parent) + 1,)
 Base.axes(psi::CView{<:AbstractFiniteMPS}) = map(n -> 0:(n - 1), size(psi))
 
 Base.size(psi::CView{<:Multiline{<:InfiniteMPS}}) = size(psi.parent)
-function Base.size(psi::CView{<:Multiline{<:AbstractFiniteMPS}})
-    return (length(psi.parent.data), length(first(psi.parent.data)) + 1)
-end
-function Base.axes(psi::CView{<:Multiline{<:AbstractFiniteMPS}})
-    return (Base.OneTo(length(psi.parent.data)), 0:length(first(psi.parent.data)))
-end
 
 #the checkbounds for multiline objects needs to be changed, as the first index is periodic
 #however if it is a Multiline(Infinitemps), then the second index is also periodic!
@@ -291,10 +285,10 @@ function Base.checkbounds(
         psi::Union{ACView{<:Multiline}, ALView{<:Multiline}, ARView{<:Multiline}, CView{<:Multiline}},
         a, b
     )
-    return if first(psi.parent.data) isa InfiniteMPS
+    return if first(parent(psi.parent)) isa InfiniteMPS
         true
     else
-        checkbounds(Bool, CView(first(psi.parent.data)), b)
+        checkbounds(Bool, CView(first(parent(psi.parent))), b)
     end
 end