diff --git a/docs/src/changelog.md b/docs/src/changelog.md index eeae4041f..cc2cf42e5 100644 --- a/docs/src/changelog.md +++ b/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 diff --git a/docs/src/examples/classic2d/1.hard-hexagon/index.md b/docs/src/examples/classic2d/1.hard-hexagon/index.md index 6f103a911..b57ef308a 100644 --- a/docs/src/examples/classic2d/1.hard-hexagon/index.md +++ b/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ξ = $ξ") diff --git a/docs/src/man/algorithms.md b/docs/src/man/algorithms.md index 02b078b61..8d94fa73d 100644 --- a/docs/src/man/algorithms.md +++ b/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 diff --git a/docs/src/man/operators.md b/docs/src/man/operators.md index 05525add2..837bb12f2 100644 --- a/docs/src/man/operators.md +++ b/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] * ψ))) +``` diff --git a/docs/src/man/states.md b/docs/src/man/states.md index 168b2117c..3e4e30706 100644 --- a/docs/src/man/states.md +++ b/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). diff --git a/examples/classic2d/1.hard-hexagon/main.jl b/examples/classic2d/1.hard-hexagon/main.jl index c7b3531b1..436f7603a 100644 --- a/examples/classic2d/1.hard-hexagon/main.jl +++ b/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ξ = $ξ") diff --git a/src/MPSKit.jl b/src/MPSKit.jl index f900884ae..d09a9e473 100644 --- a/src/MPSKit.jl +++ b/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") diff --git a/src/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl index 5cfff2de8..3e17ffbce 100644 --- a/src/algorithms/changebonds/randexpand.jl +++ b/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 diff --git a/src/algorithms/changebonds/svdcut.jl b/src/algorithms/changebonds/svdcut.jl index 295e48557..fc7cbb4ed 100644 --- a/src/algorithms/changebonds/svdcut.jl +++ b/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) diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl index 191d2f93c..c3961d746 100644 --- a/src/algorithms/expval.jl +++ b/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 # -------------- diff --git a/src/algorithms/grassmann.jl b/src/algorithms/grassmann.jl index 7da1a0133..cc73839c0 100644 --- a/src/algorithms/grassmann.jl +++ b/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!, NoTimerOutput, @timeit, default_allocator using TensorOperations: AbstractBackend, DefaultBackend using TensorKit @@ -195,7 +195,7 @@ function fg( scheduler::Scheduler = MPSKit.Defaults.scheduler[], ) where {O <: InfiniteMPO} @timeit timeroutput "envs (parallel)" recalculate!(envs, state, operator, state; timeroutput) - f = @timeit timeroutput "expval" expectation_value(state, operator, envs) + f = @timeit timeroutput "expval" dominant_eigenvalue(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)}) @@ -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) + f = @timeit timeroutput "expval" dominant_eigenvalue(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 diff --git a/src/algorithms/groundstate/vumps.jl b/src/algorithms/groundstate/vumps.jl index 7be8c096a..ed0b030b1 100644 --- a/src/algorithms/groundstate/vumps.jl +++ b/src/algorithms/groundstate/vumps.jl @@ -75,21 +75,24 @@ function dominant_eigsolve( state = VUMPSState(mps, operator, envs, iter, ϵ, which, timeroutput) it = IterativeSolver(alg, state) + # both `find_groundstate` and `leading_boundary` go here + objective = mps isa MultilineMPS ? dominant_eigenvalue : expectation_value + result = LoggingExtras.withlevel(; alg.verbosity) do - @infov 2 loginit!(log, ϵ, sum(expectation_value(mps, operator, envs))) + @infov 2 loginit!(log, ϵ, objective(mps, operator, envs)) for (mps, envs, ϵ) in it if ϵ ≤ alg.tol @infov 4 TimerReport(timeroutput) - @infov 2 logfinish!(log, it.iter, ϵ, expectation_value(mps, operator, envs)) + @infov 2 logfinish!(log, it.iter, ϵ, objective(mps, operator, envs)) return mps, envs, ϵ end if it.iter ≥ alg.maxiter @infov 4 TimerReport(timeroutput) - @warnv 1 logcancel!(log, it.iter, ϵ, expectation_value(mps, operator, envs)) + @warnv 1 logcancel!(log, it.iter, ϵ, objective(mps, operator, envs)) return mps, envs, ϵ end - @infov 3 logiter!(log, it.iter, ϵ, expectation_value(mps, operator, envs)) + @infov 3 logiter!(log, it.iter, ϵ, objective(mps, operator, envs)) end # this should never be reached diff --git a/src/algorithms/statmech/idmrg.jl b/src/algorithms/statmech/idmrg.jl index 84edab763..3f1039505 100644 --- a/src/algorithms/statmech/idmrg.jl +++ b/src/algorithms/statmech/idmrg.jl @@ -7,7 +7,7 @@ function leading_boundary( local iter LoggingExtras.withlevel(; alg.verbosity) do - @infov 2 loginit!(log, ϵ, expectation_value(ψ, operator, envs)) + @infov 2 loginit!(log, ϵ, dominant_eigenvalue(ψ, operator, envs)) for outer iter in 1:(alg.maxiter) alg_eigsolve = adapt_solver(alg.alg_eigsolve; iter, g_global = ϵ) C_current = ψ.C[:, 0] @@ -44,13 +44,13 @@ function leading_boundary( ϵ = norm(C_current - ψ.C[:, 0]) if ϵ < alg.tol - @infov 2 logfinish!(log, iter, ϵ, expectation_value(ψ, operator, envs)) + @infov 2 logfinish!(log, iter, ϵ, dominant_eigenvalue(ψ, operator, envs)) break end if iter == alg.maxiter - @warnv 1 logcancel!(log, iter, ϵ, expectation_value(ψ, operator, envs)) + @warnv 1 logcancel!(log, iter, ϵ, dominant_eigenvalue(ψ, operator, envs)) else - @infov 3 logiter!(log, iter, ϵ, expectation_value(ψ, operator, envs)) + @infov 3 logiter!(log, iter, ϵ, dominant_eigenvalue(ψ, operator, envs)) end end end diff --git a/src/algorithms/statmech/vomps.jl b/src/algorithms/statmech/vomps.jl index 9427f5488..7fd08aef4 100644 --- a/src/algorithms/statmech/vomps.jl +++ b/src/algorithms/statmech/vomps.jl @@ -68,18 +68,18 @@ function dominant_eigsolve( it = IterativeSolver(alg, state) return LoggingExtras.withlevel(; alg.verbosity) do - @infov 2 loginit!(log, ϵ, sum(expectation_value(mps, operator, envs))) + @infov 2 loginit!(log, ϵ, dominant_eigenvalue(mps, operator, envs)) for (mps, envs, ϵ) in it if ϵ ≤ alg.tol - @infov 2 logfinish!(log, it.iter, ϵ, expectation_value(mps, operator, envs)) + @infov 2 logfinish!(log, it.iter, ϵ, dominant_eigenvalue(mps, operator, envs)) return mps, envs, ϵ end if it.iter ≥ alg.maxiter - @warnv 1 logcancel!(log, it.iter, ϵ, expectation_value(mps, operator, envs)) + @warnv 1 logcancel!(log, it.iter, ϵ, dominant_eigenvalue(mps, operator, envs)) return mps, envs, ϵ end - @infov 3 logiter!(log, it.iter, ϵ, expectation_value(mps, operator, envs)) + @infov 3 logiter!(log, it.iter, ϵ, dominant_eigenvalue(mps, operator, envs)) end # this should never be reached diff --git a/src/environments/multiline_envs.jl b/src/environments/multiline_envs.jl index a8e1d73f4..d100009e2 100644 --- a/src/environments/multiline_envs.jl +++ b/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) diff --git a/src/operators/multilinempo.jl b/src/operators/multilinempo.jl index 7432cf91a..81870e192 100644 --- a/src/operators/multilinempo.jl +++ b/src/operators/multilinempo.jl @@ -1,50 +1,57 @@ # MultilineMPO # ------------ +#TODO: add algorithm support for finite MPOs +const _MPOs = Union{InfiniteMPO, FiniteMPO} + """ - const MultilineMPO = Multiline{<:AbstractMPO} + const MultilineMPO = Multiline{<:Union{InfiniteMPO, FiniteMPO}} -Type that represents multiple lines of `MPO` objects. +Type that represents multiple lines of `MPO` objects, i.e. the rows of a two-dimensional +tensor network. Lines are restricted to `InfiniteMPO` or `FiniteMPO` objects as `MultilineMPO` +represents rows of a statistical mechanical transfer operator. +See the manual on [MultilineMPO](@ref) for details. # Constructors - MultilineMPO(mpos::AbstractVector{<:Union{SparseMPO, DenseMPO}}) - MultilineMPO(Os::AbstractMatrix{<:MPOTensor}) + MultilineMPO(mpos::AbstractVector{<:Union{InfiniteMPO, FiniteMPO}}) + MultilineMPO(Os::PeriodicMatrix{<:MPOTensor}) + MultilineMPO(t::MPOTensor) + +!!! note "Finite lines" + Finite lines are accepted by the type and by the constructors, but no algorithm supports + them yet, so they will fail somewhere further down. This is on purpose: there is currently + no support for finite multiline boundaries, but this allows them to be built and inspected. # See also -[`Multiline`](@ref), [`AbstractMPO`](@ref) +[`Multiline`](@ref), [`MultilineMPS`](@ref), [`dominant_eigenvalue`](@ref) """ -const MultilineMPO = Multiline{<:AbstractMPO} +const MultilineMPO = Multiline{<:_MPOs} -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{<:_MPOs}) = 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::_MPOs) = 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) +Base.convert(::Type{FiniteMPO}, t::MultilineMPO{<:FiniteMPO}) = only(t) -function Base.:*(mpo::MultilineMPO, st::MultilineMPS) - size(mpo) == size(st) || throw(ArgumentError("dimension mismatch")) - return Multiline(map(*, zip(mpo, st))) -end - -function Base.:*(mpo1::MultilineMPO, mpo2::MultilineMPO) - size(mpo1) == size(mpo2) || throw(ArgumentError("dimension mismatch")) - return Multiline(map(*, zip(mpo1, mpo2))) +function Base.:*(mpo::MultilineMPO, st::InfiniteMPS) + check_length(mpo[1], st) + for i in 1:size(mpo, 1) + st = mpo[i] * st + end + return st end for f_space in (:physicalspace, :left_virtualspace, :right_virtualspace) diff --git a/src/states/multilinemps.jl b/src/states/multilinemps.jl index 524280ffc..2a0069df6 100644 --- a/src/states/multilinemps.jl +++ b/src/states/multilinemps.jl @@ -1,15 +1,20 @@ # MultilineMPS # ------------ -const MultilineMPS = Multiline{<:InfiniteMPS} +#TODO: add support for finite MPS +const _MPSs = Union{InfiniteMPS, FiniteMPS} +const MultilineMPS = Multiline{<:_MPSs} @doc """ - const MultilineMPS = Multiline{<:InfiniteMPS} + const MultilineMPS = Multiline{<:Union{InfiniteMPS, FiniteMPS}} -Type that represents multiple lines of [`InfiniteMPS`](@ref) objects. +Type that represents multiple lines of MPS objects. When used in the context of +[`leading_boundary`](@ref) with `InfiniteMPS`, this is not to be confused with the fixed point +of a 2D tensor network, which is a single `InfiniteMPS`. +See the manual on [MultilineMPS](@ref) for details. # Constructors - MultilineMPS(mpss::AbstractVector{<:InfiniteMPS}) + MultilineMPS(mpss::AbstractVector{<:Union{InfiniteMPS, FiniteMPS}}) MultilineMPS( [f, eltype], physicalspaces::Matrix{<:Union{S, CompositeSpace{S}}}, virtualspaces::Matrix{<:Union{S, CompositeSpace{S}}} @@ -27,13 +32,27 @@ Type that represents multiple lines of [`InfiniteMPS`](@ref) objects. - `AC`: center-gauged MPS tensors - `C`: gauge (bond) tensors +# Notes + +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. + +Only the first constructor accepts finite lines; the others build `InfiniteMPS` lines from +spaces or tensors. + +!!! note "Finite lines" + Finite lines are accepted by the type and by the first constructor, but no algorithm + supports them yet, so they will fail somewhere further down. This is on purpose: there is currently + no support for finite multiline boundaries, but this allows them to be built and inspected. + # See also -[`Multiline`](@ref) +[`Multiline`](@ref), [`MultilineMPO`](@ref) """ function MultilineMPS end -MultilineMPS(mpss::AbstractVector{<:InfiniteMPS}) = Multiline(mpss) +MultilineMPS(mpss::AbstractVector{<:_MPSs}) = Multiline(mpss) function MultilineMPS( pspaces::AbstractMatrix{S}, Dspaces::AbstractMatrix{S}; kwargs... ) where {S <: VectorSpace} @@ -99,13 +118,11 @@ function TensorKit.dot(a::MultilineMPS, b::MultilineMPS; kwargs...) end 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.convert(::Type{MultilineMPS}, st::_MPSs) = Multiline([st]) +Base.convert(::Type{InfiniteMPS}, st::MultilineMPS{<:InfiniteMPS}) = only(st) +Base.convert(::Type{FiniteMPS}, st::MultilineMPS{<:FiniteMPS}) = only(st) 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)...) diff --git a/src/states/orthoview.jl b/src/states/orthoview.jl index 4bbddf2a1..fb0c82604 100644 --- a/src/states/orthoview.jl +++ b/src/states/orthoview.jl @@ -278,25 +278,28 @@ 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) + return (length(parent(psi.parent)), length(first(parent(psi.parent))) + 1) end function Base.axes(psi::CView{<:Multiline{<:AbstractFiniteMPS}}) - return (Base.OneTo(length(psi.parent.data)), 0:length(first(psi.parent.data))) + return (Base.OneTo(length(parent(psi.parent))), 0:length(first(parent(psi.parent)))) 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! +const MultilineOrthoView{S} = Union{ACView{S}, ALView{S}, ARView{S}, CView{S}} + +# the row index is always periodic +# for infinite lines the column index is periodic as well +Base.checkbounds(::Type{Bool}, ::MultilineOrthoView{<:Multiline{<:InfiniteMPS}}, a, b) = true + +# finite line bounds differ per view, so just delegate to the view's own checkbounds function Base.checkbounds( - ::Type{Bool}, - psi::Union{ACView{<:Multiline}, ALView{<:Multiline}, ARView{<:Multiline}, CView{<:Multiline}}, - a, b + ::Type{Bool}, psi::MultilineOrthoView{<:Multiline{<:AbstractFiniteMPS}}, a, b ) - return if first(psi.parent.data) isa InfiniteMPS - true - else - checkbounds(Bool, CView(first(psi.parent.data)), b) - end + return checkbounds(Bool, _lineview(psi), b) end +_lineview(psi::ACView) = ACView(first(parent(psi.parent))) +_lineview(psi::ALView) = ALView(first(parent(psi.parent))) +_lineview(psi::ARView) = ARView(first(parent(psi.parent))) +_lineview(psi::CView) = CView(first(parent(psi.parent))) # Gauging routines # ---------------- diff --git a/src/utility/multiline.jl b/src/utility/multiline.jl index ed4107cf2..e50e40157 100644 --- a/src/utility/multiline.jl +++ b/src/utility/multiline.jl @@ -2,7 +2,20 @@ $(TYPEDEF) Object that represents multiple lines of objects of type `T`. Typically used to represent -multiple lines of `InfiniteMPS` (`MultilineMPS`) or MPO (`Multiline{<:AbstractMPO}`). +multiple lines of `InfiniteMPS` (`MultilineMPS`) or `InfiniteMPO` (`MultilineMPO`). + +`Multiline` plays two different, orthogonal roles at once, and its Base overloads are split +accordingly: + +- As a sequence of lines, matching what is actually stored: `length`, `eltype`, `iterate` + and `m[i]` (a single integer index) all refer to the `T`-typed lines themselves, i.e. + `length(m) == nrows` and `m[i]::T`. +- As a lattice, describing the 2D shape spanned by the lines together: `size(m)` is + `(nrows, ncols)`, and `axes`/`eachindex` follow `size`. + +These two views disagree on purpose (`length(m) != prod(size(m))`). +Code that wants to work line-by-line should use `m[i]`/`parent(m)`, +while code that wants the lattice shape should use `size`. # Fields @@ -25,14 +38,17 @@ Multiline(data::AbstractVector{T}) where {T} = Multiline{T}(data) # ----------------------- Base.parent(m::Multiline) = m.data Base.size(m::Multiline) = (length(parent(m)), length(parent(m)[1])) -Base.size(m::Multiline, i::Int) = i == 1 ? length(parent(m)) : i == 2 ? length(parent(m)[1]) : error() -Base.length(m::Multiline) = prod(size(m)) -function Base.axes(m::Multiline, i::Int) - return i == 1 ? axes(parent(m), 1) : - i == 2 ? axes(parent(m)[1], 1) : throw(ArgumentError("Invalid index $i")) +function Base.size(m::Multiline, i::Int) # acts like abstract array + return i == 1 ? length(parent(m)) : i == 2 ? length(parent(m)[1]) : 1 +end +Base.length(m::Multiline) = length(parent(m)) +function Base.axes(m::Multiline, d::Int) + return d <= 2 ? axes(m)[d] : Base.OneTo(1) # matches size end Base.eachindex(m::Multiline) = CartesianIndices(size(m)) Base.isfinite(m::Multiline) = isfinite(typeof(m)) +Base.isfinite(::Type{Multiline{T}}) where {T} = isfinite(T) +Base.eltype(::Type{Multiline{T}}) where {T} = T eachsite(m::Multiline) = eachsite(first(parent(m))) @@ -53,7 +69,7 @@ Base.reverse(A::Multiline) = Multiline(reverse(parent(A))) Base.only(A::Multiline) = only(parent(A)) function Base.repeat(A::Multiline, rows::Int, cols::Int) - inner = map(Base.Fix2(repeat, cols), A.data) + inner = map(Base.Fix2(repeat, cols), parent(A)) outer = repeat(inner, rows) return Multiline(outer) end @@ -102,6 +118,7 @@ end VectorInterface.add!!(x::Multiline, y::Multiline, α::Number, β::Number) = add!(x, y, α, β) +# FIXME? is it intentional that a nontrivial multilinemps of normalised rows never has norm 1? function VectorInterface.inner(x::Multiline, y::Multiline) T = VectorInterface.promote_inner(x, y) init = zero(T) @@ -117,6 +134,7 @@ site_type(::Type{Multiline{S}}) where {S} = site_type(S) bond_type(::Type{Multiline{S}}) where {S} = bond_type(S) site_type(st::Multiline) = site_type(typeof(st)) bond_type(st::Multiline) = bond_type(typeof(st)) -TensorKit.sectortype(::Type{Multiline{T}}) where {T} = sectortype(T) -TensorKit.spacetype(::Type{Multiline{T}}) where {T} = spacetype(T) -TensorKit.storagetype(::Type{Multiline{T}}) where {T} = storagetype(T) +for ftype in (:spacetype, :sectortype, :storagetype) + @eval TensorKit.$ftype(::Type{Multiline{T}}) where {T} = $ftype(T) + @eval TensorKit.$ftype(m::Multiline) = $ftype(typeof(m)) +end diff --git a/src/utility/show.jl b/src/utility/show.jl index 2709b7d43..5228c41db 100644 --- a/src/utility/show.jl +++ b/src/utility/show.jl @@ -10,6 +10,23 @@ for T in (:FiniteMPS, :InfiniteMPS, :FiniteMPO, :InfiniteMPO, :FiniteMPOHamilton end end +function Base.summary(io::IO, m::MultilineMPS) + R, C = size(m) + D = maximum(dim, left_virtualspace(m)) + E = scalartype(m) + S = TensorKit.type_repr(spacetype(m)) + print(io, R, "×", C, " MultilineMPS(", E, ", ", S, ") with maximal dimension ", D) + return nothing +end +function Base.summary(io::IO, m::MultilineMPO) + R, C = size(m) + D = maximum(dim, left_virtualspace(m)) + E = scalartype(m) + S = TensorKit.type_repr(spacetype(m)) + print(io, R, "×", C, " MultilineMPO(", E, ", ", S, ") with maximal dimension ", D) + return nothing +end + function Base.show(io::IO, ::MIME"text/plain", ψ::FiniteMPS) summary(io, ψ) get(io, :compact, false)::Bool && return nothing @@ -164,6 +181,40 @@ function Base.show(io::IO, ::MIME"text/plain", mpo::AbstractMPO) return nothing end +function Base.show(io::IO, ::MIME"text/plain", m::Union{MultilineMPS, MultilineMPO}) + summary(io, m) + get(io, :compact, false)::Bool && return nothing + println(io, ":") + + R = size(m, 1) + limit = get(io, :limit, true)::Bool + # rows are typically few so the row budget itself is small and fixed + # rather than computed from `displaysize` + half_rows = (limit && R > 4) ? 2 : R + shown = R <= 2 * half_rows ? (1:R) : [1:half_rows; (R - half_rows + 1):R] + + for (k, i) in enumerate(shown) + println(io, "row ", i, ":") + show(io, MIME"text/plain"(), m[i]) # dispatches to the line's own show for correct connectors + println(io) + if k < length(shown) + next = shown[k + 1] + if next == i + 1 + if m isa MultilineMPO + println(io, " ↓ (row ", i, " maps onto row ", next, ")") + else + println(io, " ⋮") + end + else + println(io, " ⋮\n ⋮") + end + end + end + m isa MultilineMPO && println(io, " ↓ (row ", R, " maps onto row 1)") + + return nothing +end + function Base.summary(io::IO, envs::Union{FiniteEnvironments, InfiniteEnvironments}) print(io, length(envs.GLs), "-site ") Base.showarg(io, envs, true) diff --git a/test/algorithms/statmech.jl b/test/algorithms/statmech.jl index 8a736c02b..2523c95ec 100644 --- a/test/algorithms/statmech.jl +++ b/test/algorithms/statmech.jl @@ -7,6 +7,7 @@ println(" using .TestSetup using Test, TestExtras using MPSKit +using MPSKit: Multiline using TensorKit using TensorKit: ℙ @@ -123,4 +124,17 @@ using TensorKit: ℙ e = expectation_value(ψ, (O_mpo, 1 => E)) @test e ≈ e_th atol = 1.0e-2 end + + @testset "MultilineMPS expectation_value (multi-row)" begin + # reuses the exact permutation-MPO fixed point: O[i]*ψ[i] == ψ[i+1] exactly + ψ = MultilineMPS([product_mps(1), product_mps(2), product_mps(3)]) + O = MultilineMPO([perm_mpo([2, 3, 1]), perm_mpo([1, 3, 2]), perm_mpo([3, 2, 1])]) + @test dominant_eigenvalue(ψ, O) ≈ 1 atol = 1.0e-10 + + # `Multiline([H, H])` of Hamiltonian bypasses the MultilineMPO(...) construction guard + # expectation_value must still reject it + H = transverse_field_ising() + ψ2 = MultilineMPS([InfiniteMPS(ℂ^2, ℂ^4), InfiniteMPS(ℂ^2, ℂ^4)]) + @test_throws MethodError expectation_value(ψ2, Multiline([H, H])) + end end diff --git a/test/operators/multilinempo.jl b/test/operators/multilinempo.jl new file mode 100644 index 000000000..2f561d86a --- /dev/null +++ b/test/operators/multilinempo.jl @@ -0,0 +1,56 @@ +println(" +---------------------------- +| MultilineMPO tests | +---------------------------- +") + +using .TestSetup +using Test, TestExtras +using MPSKit +using TensorKit + +d = 3 +P = ℂ^d + +# ψ[i] = |i⟩ (3 rows, 1 column), O[i]|i⟩ = |i + 1⟩ exactly +# ψ[i] is an InfiniteMPS and is periodic in the line index i +# this properly tests row-shift convention: row i of a MultilineMPO maps row i of the network onto row i + 1 +ψ = MultilineMPS([product_mps(1), product_mps(2), product_mps(3)]) +O = MultilineMPO( + [ + perm_mpo([2, 3, 1]), # row 1: |1⟩ -> |2⟩ + perm_mpo([1, 3, 2]), # row 2: |2⟩ -> |3⟩ + perm_mpo([3, 2, 1]), # row 3: |3⟩ -> |1⟩ + ] +) + +@testset "MultilineMPO * InfiniteMPS" begin + for i in 1:3 + @test abs(dot(ψ[i + 1], O[i] * ψ[i])) ≈ 1 atol = 1.0e-10 + @test abs(dot(ψ[i], O[i] * ψ[i])) ≈ 0 atol = 1.0e-10 + end + + @test abs(dot(ψ[1], O * ψ[1])) ≈ 1 atol = 1.0e-10 + @test abs(dot(ψ[2], O * ψ[2])) ≈ 1 atol = 1.0e-10 # cyclic period + + # the order matters, here testing reverse order not returning to ψ[1] + reversed = foldl((st, i) -> O[i] * st, 3:-1:1; init = ψ[1]) + @test abs(dot(ψ[1], reversed)) ≈ 0 atol = 1.0e-10 +end + +@testset "dominant_eigenvalue" begin + @test dominant_eigenvalue(ψ, O) ≈ 1 atol = 1.0e-10 + + # row-accumulated eigenvalues + scaled(mpo, c) = InfiniteMPO([c * mpo[i] for i in 1:length(mpo)]) + cs = (2.0, 3.0, 5.0) + O_rows = MultilineMPO([scaled(O[i], cs[i]) for i in 1:3]) + @test dominant_eigenvalue(ψ, O_rows) ≈ prod(cs) atol = 1.0e-8 + + # column-accumulated eigenvalues + id_mpo = perm_mpo([1, 2, 3]) # identity operator, so anything's a fixed point + cs = (2.0, 3.0) + O_cols = MultilineMPO([InfiniteMPO([cs[1] * id_mpo[1], cs[2] * id_mpo[1]])]) + ψ_cols = MultilineMPS([InfiniteMPS([P, P], [ℂ^2, ℂ^2])]) + @test dominant_eigenvalue(ψ_cols, O_cols) ≈ prod(cs) atol = 1.0e-8 +end diff --git a/test/setup/testsetup.jl b/test/setup/testsetup.jl index 5aacebd2b..d9f7761a9 100644 --- a/test/setup/testsetup.jl +++ b/test/setup/testsetup.jl @@ -24,6 +24,7 @@ export transverse_field_ising, heisenberg_XXX, bilinear_biquadratic_model, XY_mo export classical_ising_tensors, classical_ising, sixvertex export bad_initial_state export SCHEDULERS, with_scheduler +export perm_mpo, product_mps # using TensorOperations @@ -291,4 +292,23 @@ function bad_initial_state(H, L; T = ComplexF64, n_states = 20, n_fixed = 3) end end +# functions for multiline tests +# bond-dim-1 permutation MPO with O|p⟩ = |π[p]⟩ +function perm_mpo(pi) + d = 3 + P, V = ℂ^d, ℂ^1 + t = zeros(ComplexF64, 1, d, d, 1) + for pin in 1:d + t[1, pi[pin], pin, 1] = 1.0 + end + return InfiniteMPO([TensorMap(t, V ⊗ P ← P ⊗ V)]) +end +function product_mps(k) + d = 3 + P, V = ℂ^d, ℂ^1 + a = zeros(ComplexF64, 1, d, 1) + a[1, k, 1] = 1.0 + return InfiniteMPS([TensorMap(a, V ⊗ P ← V)]) +end + end diff --git a/test/states/multilinemps.jl b/test/states/multilinemps.jl index cf5456f6f..c02e4eb4d 100644 --- a/test/states/multilinemps.jl +++ b/test/states/multilinemps.jl @@ -27,6 +27,7 @@ using TensorKit: ℙ @test TensorKit.sectortype(ψ) == sectortype(D) @test !isfinite(typeof(ψ)) + @test !isfinite(ψ) @test physicalspace(ψ) == fill(d, 2, 2) @test all(x -> x ≾ D, left_virtualspace(ψ))