Fix Multiline (MultilineMPS/MultilineMPO) correctness and API inconsistencies

docs/src/changelog.md

@@ -38,7 +38,11 @@ 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`.
+- `*(::MultilineMPO, ::InfiniteMPS)`, which pushes the boundary MPS through every row of the
+  network in turn, advancing it by one full period.
+- `dominant_eigenvalue(ψ, O, [environments])`, the eigenvalue of the transfer operator `O` for the
+  boundary MPS `ψ`. `expectation_value(::InfiniteMPS, ::InfiniteMPO)` forwards here.
 ### Changed
 
 - Renormalization during time evolution is now controlled by an explicit `normalize` keyword on
@@ -72,16 +76,38 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
   not the requested correlator. ([#489](https://github.com/QuantumKitHub/MPSKit.jl/pull/489))
 - TimerOutputs 1.x is now required. The timing tables printed at `verbosity > 3` use the new
   layout (tree guides, heat bars) and additionally report per-section GC time.
+- `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` and `MultilineMPS` lines are now restricted by the type to
+  `Union{InfiniteMPO, FiniteMPO}` and `Union{InfiniteMPS, FiniteMPS}` respectively, rather than to
+  any `AbstractMPO`/`InfiniteMPS`. Hamiltonian lines are excluded outright. Finite lines are
+  accepted by both the type and the constructors so that finite multiline networks can be built
+  and inspected. No algorithm supports them yet, so they fail further down. The `AbstractMatrix`
+  constructor that silently built finite-line `MultilineMPO`s was removed.
 
 ### Deprecated
 
 ### Removed
 
 - Support for TimerOutputs 0.5.
+- `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.
+- `expectation_value` for a `MultilineMPS`/`MultilineMPO` pair entirely, replaced by
+  `dominant_eigenvalue`.
+- `*(::MultilineMPO, ::MultilineMPS)` and `*(::MultilineMPO, ::MultilineMPO)`, as these
+  were not meaningful operations. Neither method had ever been callable previously.
 
 ### Fixed
 
 - `isfinite(::WindowMPOHamiltonian)` was undefined. ([#489](https://github.com/QuantumKitHub/MPSKit.jl/pull/489))
+- `checkbounds` on the `AL`/`AR`/`AC`/`C` views of a `Multiline` now delegates to the
+  matching view, and dispatches on `Multiline{<:InfiniteMPS}` versus `Multiline{<:AbstractFiniteMPS}`.
+  The row index remains unchecked in both cases due to periodicity.
+- `size`/`axes` for a `CView` over a `Multiline` with finite lines were missing.
 - `excitations(::InfiniteMPO, ::QuasiparticleAnsatz, ::InfiniteQP, lenvs, renvs)` referenced `H_eff`  before assigning. ([#489](https://github.com/QuantumKitHub/MPSKit.jl/pull/489))
 - `Base.:+`/`-` on `FiniteMPS` returned a wrong state for near-parallel operands carried by
   different tensor networks, e.g. `norm(E₀ * gs - H * gs)` coming out as `2 * norm(gs) * E₀`
@@ -104,6 +130,13 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 - Fix hardcoding of number of physical spaces in the `changebonds` implementations for
   `FiniteMPS`, enabling its use for systems with composite physical spaces
   ([#514](https://github.com/QuantumKitHub/MPSKit.jl/pull/514))
+- `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/examples/classic2d/1.hard-hexagon/index.md

@@ -41,7 +41,7 @@ end
 
 One way to study statistical mechanics in infinite systems with tensor networks is by approximating the dominant eigenvector of the transfer matrix by an MPS.
 This dominant eigenvector contains a lot of hidden information.
-For example, the free energy can be extracted by computing the expectation value of the mpo.
+For example, the partition function per site is the dominant eigenvalue of the transfer matrix, from which the free energy follows.
 Additionally, we can compute the entanglement entropy as well as the correlation length of the state:
 
 ````julia
@@ -52,7 +52,7 @@ V = virtual_space(D)
     ψ₀, mpo,
     VUMPS(; verbosity = 0, alg_eigsolve = MPSKit.Defaults.alg_eigsolve(; ishermitian = false))
 ) # use non-hermitian eigensolver
-F = real(expectation_value(ψ, mpo))
+F = real(dominant_eigenvalue(ψ, mpo))
 S = real(first(entropy(ψ)))
 ξ = correlation_length(ψ; sector = leftunit(ψ))
 println("F = $F\tS = $S\tξ = $ξ")

docs/src/man/algorithms.md

@@ -22,7 +22,6 @@ Many of these algorithms have different advantages and disadvantages, and figuri
 Here, we enumerate some of their properties in hopes of pointing you in the right direction. For convenience, the full list of algorithms is:
 
 - [DMRG](@ref)
-- [DMRG2](@ref)
 - [VUMPS](@ref)
 - [Gradient descent](@ref)
 - [TDVP](@ref)
@@ -329,6 +328,14 @@ use GradientGrassmann.
 leading_boundary
 ```
 
+The associated eigenvalue, i.e. the partition function per unit cell, is obtained with [`dominant_eigenvalue`](@ref).
+This is deliberately not `expectation_value`: the boundary contracts against a different line than the one the operator acts on, so it is not an overlap.
+For a single-line `InfiniteMPS`/`InfiniteMPO` pair the two coincide and `expectation_value` forwards accordingly.
+
+```@docs; canonical=false
+dominant_eigenvalue
+```
+
 ## `approximate`
 
 Often, it is useful to approximate a given MPS by another, typically by one of a different

docs/src/man/operators.md

@@ -53,8 +53,7 @@ MPOs also support a range of linear algebra operations, such as addition, subtra
 multiplication, either among themselves or with a finite MPS. Here, it is important to note
 that these operations will increase the virtual dimension of the resulting MPO or MPS, and
 this naive application is thus typically not optimal. For approximate operations that do not
-increase the virtual dimension, the more advanced algorithms in the [um_algorithms](@ref)
-sections should be used.
+increase the virtual dimension, the more advanced algorithms in the [algorithms](@ref um_algorithms) sections should be used.
 
 ```@example operators
 O_xzx² = O_xzx * O_xzx
@@ -274,3 +273,23 @@ 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
+
+A [`MultilineMPO`](@ref) is a stack of MPO lines making up the rows of a two-dimensional tensor network.
+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.
+
+```@example operators
+mpo_multi = MultilineMPO([mpo, mpo])
+```
+
+The type restricts lines to `InfiniteMPO` and `FiniteMPO`, since a `MultilineMPO` is a statistical mechanical transfer operator.
+The quantity extracted from it is a [`dominant_eigenvalue`](@ref) rather than an energy.
+Finite lines are accepted so that finite networks can be built and inspected, but no algorithm supports them yet.
+
+Applying a `MultilineMPO` pushes an [`InfiniteMPS`](@ref) through every row in turn, advancing the boundary by one full period of the network:
+
+```julia
+O * ψ == O[end] * (… * (O[2] * (O[1] * ψ)))
+```

docs/src/man/states.md

@@ -155,20 +155,48 @@ A two-dimensional classical partition function can often be represented by an in
 There are many ways to evaluate such a network, but here we focus on the so-called boundary MPS methods.
 These first reduce the problem from contracting a two-dimensional network to the contraction of a one-dimensional MPS, by finding the fixed point of the row-to-row (or column-to-column) transfer matrix.
 In these cases however, there might be a non-trivial periodicity in both the horizontal as well as vertical direction.
-Therefore, in MPSKit they are represented by [`MultilineMPS`](@ref), which are simply a repeating set of [`InfiniteMPS`](@ref).
+Therefore, in MPSKit they are represented by [`MultilineMPS`](@ref), which are simply a repeating set of MPS lines, one per row of the network.
 
 ```@example states
 state = MultilineMPS(fill(infinite_state, 2))
 ```
 
-They offer some convenience functionality for using cartesian indexing (row - column):
+Properties are accessed with cartesian (row, column) indexing:
 
-You can access properties by calling
 ```@example states
-row = 2
-col = 2
-al = state.AL[row, col];
+al = state.AL[2, 1];
 ```
 
+### The row-shift convention
+
+Row `i` of a [`MultilineMPO`](@ref) maps row `i` of the network onto row `i + 1`.
+`environments`, the derivative operators and [`dominant_eigenvalue`](@ref) all follow it, pairing `state[i + 1]` as the bra against `O[i]` acting on the ket `state[i]`.
+Within a row, `state.AL[i, j]` and friends behave exactly as they would for `state[i]::InfiniteMPS`.
+
+Bra and ket are therefore different lines, so a `MultilineMPS` is not a state whose expectation value makes sense to take.
+
+### One fixed point, many lines
+
+It is worth mentioning explicitly what the point of the extra rows is, since it is easy to read a `MultilineMPS` as a stack of independent states.
+It is not.
+However many rows the operator has, the network has exactly one boundary fixed point, namely the MPS that returns to itself after being pushed through all of the rows.
+The lines of a `MultilineMPS` are bookkeeping for that single problem.
+Line `i + 1` is the boundary after row `i` has been applied, so this lets [`leading_boundary`](@ref) cut it into `nrows` coupled subproblems instead of contracting every row into one operator.
+Each subproblem contributes a partial factor, and only their product is the eigenvalue of the fixed point.
+Independently of that, the `ncols` sites of the unit cell contribute factors of their own on top.
+[`dominant_eigenvalue`](@ref) accumulates both.
+
+Applying an operator therefore acts on an ordinary [`InfiniteMPS`](@ref), advancing it by one full period:
+
+```julia
+O * ψ == O[end] * (… * (O[2] * (O[1] * ψ)))
+```
+
+### Subtleties
+
+- **`size` vs. iteration:** `size(state)` is the `(nrows, ncols)` lattice shape, while `length`, iteration and `state[i]` refer to the lines, so `length(state) == nrows`. See [`Multiline`](@ref MPSKit.Multiline).
+- **Norms:** `dot`/`norm` sum over rows, so `norm(state) == sqrt(nrows)` for `nrows` normalized rows.
+- **Finite lines** are accepted by the type and the vector constructor so that finite networks can be built and inspected, but no algorithm supports them yet.
+
 These objects are also used extensively in the context of [PEPSKit.jl](https://github.com/QuantumKitHub/PEPSKit.jl).
 

examples/classic2d/1.hard-hexagon/main.jl

@@ -37,7 +37,7 @@ md"""
 
 One way to study statistical mechanics in infinite systems with tensor networks is by approximating the dominant eigenvector of the transfer matrix by an MPS. 
 This dominant eigenvector contains a lot of hidden information.
-For example, the free energy can be extracted by computing the expectation value of the mpo.
+For example, the partition function per site is the dominant eigenvalue of the transfer matrix, from which the free energy follows.
 Additionally, we can compute the entanglement entropy as well as the correlation length of the state:
 """
 
@@ -48,7 +48,7 @@ V = virtual_space(D)
     ψ₀, mpo,
     VUMPS(; verbosity = 0, alg_eigsolve = MPSKit.Defaults.alg_eigsolve(; ishermitian = false))
 ) # use non-hermitian eigensolver
-F = real(expectation_value(ψ, mpo))
+F = real(dominant_eigenvalue(ψ, mpo))
 S = real(first(entropy(ψ)))
 ξ = correlation_length(ψ; sector = leftunit(ψ))
 println("F = $F\tS = $S\tξ = $ξ")

src/MPSKit.jl

@@ -46,6 +46,7 @@ export exact_diagonalization, fidelity_susceptibility
 
 # toolbox:
 export expectation_value, correlator, variance
+export dominant_eigenvalue
 export correlation_length, marek_gap, transfer_spectrum
 export entropy, entanglement_spectrum
 export open_boundary_conditions, periodic_boundary_conditions
@@ -116,8 +117,8 @@ include("utility/linearcombination.jl")
 # maybe we should introduce an abstract state type
 include("states/abstractmps.jl")
 include("states/infinitemps.jl")
-include("states/multilinemps.jl")
 include("states/finitemps.jl")
+include("states/multilinemps.jl")
 include("states/windowmps.jl")
 include("states/orthoview.jl")
 include("states/quasiparticle_state.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,16 @@ 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 "Multiline operators"
+    There is no method for a `MultilineMPS`/`MultilineMPO` pair. Such an operator pairs row
+    `i` of `O` with line `i` of `ψ` as the ket and line `i + 1` as the bra. A well-defined
+    overlap requires sandwiching a top fixed point with a bottom fixed point of a
+    two-dimensional tensor network. The quantity of interest for just one fixed point
+    is a dominant eigenvalue, see [`dominant_eigenvalue`](@ref).
+
+    For an `InfiniteMPS` and an `InfiniteMPO` there is only a single line. The bra and ket
+    coincide, and the two notions agree.
+
 # Examples
 
 ```jldoctest
@@ -181,26 +191,68 @@ 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...)
+    return dominant_eigenvalue(ψ, mpo, envs...)
 end
-function expectation_value(
-        ψ::MultilineMPS, O::MultilineMPO{<:InfiniteMPO},
+# fallback
+function expectation_value(ψ::AbstractMPS, mpo::AbstractMPO, envs...)
+    return dot(ψ, mpo, ψ) / dot(ψ, ψ)
+end
+
+# Dominant eigenvalues
+# --------------------
+"""
+    dominant_eigenvalue(ψ, O, [environments]) -> λ
+
+Eigenvalue of the transfer operator `O` for the boundary MPS `ψ`, accumulated over one unit
+cell of the network. See the manual at [One fixed point, many lines](@ref) on how to
+interpret the role of the lines in a `MultilineMPS` and how the eigenvalue is
+accumulated over them.
+
+# Arguments
+
+- `ψ::Union{InfiniteMPS, MultilineMPS}`: the boundary MPS
+- `O::Union{InfiniteMPO, MultilineMPO}`: the transfer operator
+- `environments`: the environments to use, calculated if not given
+
+# Returns
+
+- `λ::Number`: the eigenvalue for one unit cell. Complex in general, since a transfer
+    operator need not be Hermitian.
+
+!!! note "This is not an expectation value"
+    For a [`MultilineMPO`](@ref) with more than one row this is not an overlap. The
+    contraction pairs row `i` of the operator with line `i` of `ψ` as the ket and line
+    `i + 1` as the bra, so no `⟨ϕ|O|ϕ⟩` is computed anywhere.
+    A genuine expectation value of a two-dimensional operator additionally requires the
+    fixed point on the other side of the network, which is a separate
+    [`leading_boundary`](@ref) run, contracted through mixed environments.
+
+    For a single line the bra and the ket coincide and the two notions agree, which is why
+    [`expectation_value`](@ref) forwards here for `InfiniteMPS`/`InfiniteMPO`.
+
+# See also
+
+[`leading_boundary`](@ref), [`expectation_value`](@ref)
+"""
+function dominant_eigenvalue(ψ::InfiniteMPS, O::InfiniteMPO)
+    return dominant_eigenvalue(convert(MultilineMPS, ψ), convert(MultilineMPO, O))
+end
+function dominant_eigenvalue(ψ::InfiniteMPS, O::InfiniteMPO, envs::AbstractMPSEnvironments)
+    return dominant_eigenvalue(
+        convert(MultilineMPS, ψ), convert(MultilineMPO, O), convert(MultilineEnvironments, envs)
+    )
+end
+function dominant_eigenvalue(
+        ψ::MultilineMPS{<:InfiniteMPS}, O::MultilineMPO{<:InfiniteMPO},
         envs::MultilineEnvironments = environments(ψ, O, ψ)
     )
+    #TODO: a true overlap needs the top and bottom fixed points with mixed environments
     return prod(product(1:size(ψ, 1), 1:size(ψ, 2))) do (i, j)
         GL = envs[i].GLs[j]
         GR = envs[i].GRs[j]
         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(ψ, ψ)
-end
 
 # Lazy operators
 # --------------