diff --git a/.gitignore b/.gitignore index b8339a652..6d6e49a18 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,4 @@ examples/.ipynb_checkpoints .vscode/ .DS_Store **/dev/ +tmp/ diff --git a/docs/make.jl b/docs/make.jl index 3e7b253eb..51832519a 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -64,7 +64,12 @@ makedocs(; ), pages = [ "Home" => "index.md", - "Manual" => ["man/models.md", "man/multithreading.md", "man/precompilation.md"], + "Manual" => [ + "man/models.md", + "Subspace iteration" => "man/subspace_iteration.md", + "man/multithreading.md", + "man/precompilation.md", + ], "Examples" => [ "examples/index.md", "Optimization" => joinpath.(Ref("examples"), examples_optimization), diff --git a/docs/src/man/subspace_iteration.md b/docs/src/man/subspace_iteration.md new file mode 100644 index 000000000..e3e2ffbf1 --- /dev/null +++ b/docs/src/man/subspace_iteration.md @@ -0,0 +1,308 @@ +# [Subspace-iteration CTMRG projectors](@id subspace_iteration) + +## Goal and compression + +Given a large CTMRG full-infinite environment ``M \in \mathbb{C}^{m \times m}`` and ``\chi \leq k \leq m``, we seek **range finders** ``X,Y`` satisfying + +```math +\begin{aligned} + X &\in \mathbb{C}^{m \times k}, & + X^{\dagger} X &= I_k, + \\ + Y &\in \mathbb{C}^{m \times k}, & + Y^{\dagger} Y &= I_k, +\end{aligned} +``` + +such that the leading ``\chi`` singular directions of ``M`` can be recovered from the compressed matrix + +```math +\rho = Y^{\dagger} M X \in \mathbb{C}^{k \times k}. +``` + +``\mathrm{range}(X)`` approximates the dominant right singular subspace; ``\mathrm{range}(Y)`` approximates the dominant left singular subspace. + +!!! note + + In CTMRG, ``M`` is made from two half-infinite environments, + + ```math + M = L R, \qquad L, R \in \mathbb{C}^{m \times m}. + ``` + + Since ``m`` is usually much larger than ``k``, it is cheaper to compute ``M X = L(R X)``, ``M^{\dagger} Y = R^{\dagger}(L^{\dagger} Y)``, and ``\rho = (Y^{\dagger} L)(R X)``. + +## Trial-space selection + +| `subspace` | Charge allocation | +|:--|:--| +| Elementary space | Explicit multiplicities, clipped to the available spaces. | +| Positive integer | A total dimension budget allocated using the sampled singular spectrum. | + +`trunc` independently selects the retained space: `subspace=24, trunc=truncrank(16)` allows both charge distributions to change. +Complete multiplets are preserved, so integer budgets may not be filled exactly. +Fixed-space truncation reserves its required multiplicities in the trial allocation. + +For integer budgets: + +1. Initially and every `charge_refresh=10` sweeps, probe every shared charge, reusing compatible cached directions and adding orthogonal Gaussian directions. + Exploration temporarily uses up to the trial budget per sector. +2. Apply the selected range update and reduced SVD, then keep the strongest trial directions within the total budget. +3. Between refreshes, recycle the selected space normally. + +Refreshes select directions in reduced SVD coordinates; convergence is confirmed on a refresh sweep. +`charge_probe=nothing` preserves broad exploration. Set `charge_probe=p` to add at +most `p` multiplicities per shared sector on refreshes of a compatible cache, +clipped to the same per-sector budget cap. Previously absent charges receive up to +`p` trial directions. Required fixed-space multiplicities remain reserved. +Initialization and incompatible-cache recovery always use broad exploration. +A sufficiently large `p` recovers the broad probe; small values can require several +refreshes to move many directions into a newly important charge. + +Set `charge_refresh=0` for initial and incompatible-cache exploration only. + +## Updating range finders + +Introduce notations to describe how ``X, Y`` are updated during CTMRG: + +- **Step:** ``n`` labels the CTMRG step, with incoming cached bases ``X_{n-1},Y_{n-1}`` from the previous step. + The environment ``M_n`` approaches a fixed point ``M`` as CTMRG converges. +- **Update schemes**: *Independent* schemes update ``X_n, Y_n`` independently, while in *alternating* schemes ``X_n, Y_n`` are related. +- **Initialization:** orthogonalize Gaussian samples to obtain ``X_0`` and, for the independent updates, an independent ``Y_0``. +- **Orthogonalization:** for a full-rank map ``A``, ``\mathrm{leftorth}(A)`` and ``\mathrm{rightorth}(A)`` are the orthonormal factors obtained from thin QR and LQ. + +Within each CTMRG step, ``q \geq 0`` counts extra inner iterations with ``M_n`` held fixed. +Superscript ``(j)`` labels the inner iteration. +All updates use the factorized contractions above without forming ``M``. + +### Independent scheme + +- **Initialize:** + + ```math + X_n^{(0)} = X_{n-1}, \qquad Y_n^{(0)} = Y_{n-1}. + ``` + +- **Iterate**, ``j = 1,\ldots,q + 1``: + + ```math + \begin{aligned} + X_n^{(j)} &= \mathrm{leftorth}(M_n^{\dagger} M_n X_n^{(j - 1)}), \\ + Y_n^{(j)} &= \mathrm{leftorth}(M_n M_n^{\dagger} Y_n^{(j - 1)}). + \end{aligned} + ``` + +- **Compress:** + + ```math + X_n = X_n^{(q + 1)}, \qquad Y_n = Y_n^{(q + 1)}, \qquad + \rho_n = Y_n^{\dagger} M_n X_n. + ``` + +**Complete trial spaces.** If ``k = m`` in every symmetry sector, we skip the power steps and use fixed canonical unitary bases ``X_n,Y_n``: + +```math +X_n X_n^{\dagger} = I, \qquad Y_n Y_n^{\dagger} = I, \qquad M_n = Y_n \rho_n X_n^{\dagger}. +``` + +Thus ``\rho_n`` is an exact change of coordinates: its singular values are those of ``M_n``, and lifting recovers its singular vectors. +All singular directions are already represented, so power filtering adds no information; skipping it avoids roundoff from normal products and redundant basis changes. + +### Alternating scheme + +- **Initialize:** + + ```math + Q_n^{(0)} = \mathrm{leftorth}(M_n X_{n-1}). + ``` + +- **Iterate**, ``j = 1,\ldots,q`` (skip if ``q = 0``): + + ```math + \begin{aligned} + X_n^{(j)} &= \mathrm{leftorth}(M_n^{\dagger} Q_n^{(j - 1)}), \\ + Q_n^{(j)} &= \mathrm{leftorth}(M_n X_n^{(j)}). + \end{aligned} + ``` + +- **``X`` update**: With ``Q_n = Q^{(q)}_n``, perform QR decomposition + + ```math + M^\dagger_n Q_n = X_n T_n + \qquad (X^\dagger_n X_n = I_k). + ``` + +There are two ways to update ``Y_n``. + +- **Simple ``Y`` update**: Directly set ``Y_n = Q_n``. + Then ``\rho_n`` is already given by ``T^\dagger_n``, because + + ```math + \rho_n = Y_n^{\dagger} M_n X_n + = Q^\dagger_n M_n X_n + = T^\dagger_n X^\dagger_n X_n = T^\dagger_n. + ``` + +- **Refined ``Y`` update:** Perform QR decomposition + + ```math + M_n X_n = Y_n \widetilde{T}_n \qquad + (Y_n^{\dagger} Y_n = I_k). + ``` + + Then ``\rho_n`` is already given by ``\widetilde{T}_n``, because + + ```math + \rho_n = Y_n^{\dagger} M_n X_n + = Y_n^{\dagger} Y_n \widetilde{T}_n = \widetilde{T}_n. + ``` + +### Range comparison + +All three give the same right-range recurrence in sampled sectors: + +```math +\mathrm{range}(X_n) += \mathrm{range}\!\left((M_n^{\dagger} M_n)^{q + 1} X_{n-1}\right). +``` + +Their left-range recurrences differ: + +| Variant | ``\mathrm{range}(Y_n)`` | +| --- | --- | +| Independent | ``\mathrm{range}\!\left((M_n M_n^{\dagger})^{q + 1} Y_{n-1}\right)`` | +| Alternating | ``\mathrm{range}\!\left(M_n (M_n^{\dagger} M_n)^q X_{n-1}\right)`` | +| Refined alternating | ``\mathrm{range}\!\left(M_n (M_n^{\dagger} M_n)^{q + 1} X_{n-1}\right)`` | + +## Why the dominant subspaces are recovered + +Write the SVD ``M = U \Sigma V^{\dagger}``, with singular vectors ``u_i,v_i`` and singular values ``\sigma_1 \geq \sigma_2 \geq \cdots``. +Let ``U_\chi, V_\chi \in \mathbb{C}^{m \times \chi}`` denote the leading ``\chi`` singular directions, with singular values collected in ``\Sigma_\chi``. + +### Range recurrences filter the singular directions + +For a fixed map ``M_n = M``, repeated active updates give + +```math +\mathrm{range}(X_n) += \mathrm{range}\!\left((M^{\dagger} M)^{n(q + 1)} X_0\right). +``` + +The reason this favors the dominant directions is + +```math +(M^{\dagger} M) v_i = \sigma_i^2 v_i, \qquad +(M M^{\dagger}) u_i = \sigma_i^2 u_i, \qquad +M v_i = \sigma_i u_i. +``` + +- **Right space:** each application of ``M^{\dagger} M`` suppresses a component along ``v_j`` relative to one along ``v_i`` by ``(\sigma_j/\sigma_i)^2`` when ``\sigma_j < \sigma_i``. + With full-rank initial overlap and ``\sigma_k > \sigma_{k+1}``, the range converges to the leading ``k`` right directions. +- **Left space:** the independent recurrence performs the same filtering with ``M M^{\dagger}``. + The alternating recurrences use ``M v_i = \sigma_i u_i`` to turn the filtered right directions into left directions; refinement adds one further right-to-left application. +- **During CTMRG:** ``M_n`` changes, so the fixed-map convergence statement does not apply directly. + As ``M_n`` settles, these recurrences let the cached spaces approach its dominant subspaces, provided the gap persists and updates remain active. + +Let ``U_k, V_k \in \mathbb{C}^{m \times k}`` collect the leading ``k`` singular directions of the limiting map ``M``. +Under the overlap, gap, and active-update assumptions above, as ``M_n`` converges towards ``M``, we expect + +```math +\mathrm{range}(X_n) \to \mathrm{range}(V_k), +\quad +\mathrm{range}(Y_n) \to \mathrm{range}(U_k). +``` + +Since ``k \geq \chi``, these limiting ranges contain the leading ``\chi`` directions we want to retain. + +### Compression preserves the dominant directions + +Ideally, after ``X, Y`` successfully find the dominant ``k`` singular directions, + +```math +\mathrm{range}(V_\chi) \subseteq \mathrm{range}(X), \quad +\mathrm{range}(U_\chi) \subseteq \mathrm{range}(Y). +``` + +Since ``X^{\dagger} X = Y^{\dagger} Y = I_k``, when restricting to the dominant singular subspace, + +```math +X X^{\dagger} V_\chi = V_\chi, \qquad Y Y^{\dagger} U_\chi = U_\chi. +``` + +Using ``M V_\chi = U_\chi \Sigma_\chi`` and ``M^{\dagger} U_\chi = V_\chi \Sigma_\chi`` in ``\rho = Y^{\dagger} M X`` gives + +```math +\begin{aligned} + \rho(X^{\dagger} V_\chi) &= (Y^{\dagger} U_\chi)\Sigma_\chi, \\ + \rho^{\dagger}(Y^{\dagger} U_\chi) &= (X^{\dagger} V_\chi)\Sigma_\chi. +\end{aligned} +``` + +Thus ``X^{\dagger} V_\chi`` and ``Y^{\dagger} U_\chi`` collect orthonormal right and left singular directions of ``\rho``, with the original singular values ``\Sigma_\chi``. + +Compression cannot increase any singular value, so these captured leading directions remain leading in ``\rho``. + +### Lifting returns to the full space + +Compute ``\rho = u s v^{\dagger}`` and embed all ``k`` directions back through the rangefinders: + +```math +\widehat{U} = Y u, \qquad \widehat{V} = X v. +``` + +- **Exact containment:** the leading ``\chi`` lifted directions span the target subspaces when ``\sigma_\chi > \sigma_{\chi+1}``. +- **Approximate containment:** as the range recurrences improve the captured spaces, the lifted subspaces approach the targets, provided this gap stays nonzero. + +## Rotating range finders + +After extracting the singular directions, optionally rotate all ``k`` trial directions, including oversampling directions, before caching them for the next step: + +```math +X_n \leftarrow X_n v_n = \widehat{V}_n, \qquad Y_n \leftarrow Y_n u_n = \widehat{U}_n. +``` + +The rotated range finders directly give ``Y_n^{\dagger} M_n X_n = s_n``. +The range recurrence relations are unchanged. + +- **Cache:** independent updates retain ``X_n,Y_n``; alternating variants retain only ``X_n``. + +## Convergence and differentiation + +Range finders are updated every sweep. Corner and edge spectra provide the initial +convergence check, as in ordinary CTMRG. Once these converge, scrambled bond-gauge +alignment is applied on the remaining sweeps. Convergence additionally requires +elementwise changes of the aligned environment and changes of the complete trial +subspaces at every position and direction to satisfy the CTMRG tolerance: + +```math +\epsilon(B_{n-1}, B_n) += \frac{\|B_n-B_{n-1}(B_{n-1}^{\dagger}B_n)\|_F}{\sqrt{k}}. +``` + +For orthonormal bases this is the RMS sine of the principal angles. The residual +form avoids subtracting nearly equal norms. Its maximum over the cached bases +must satisfy the CTMRG tolerance, including oversampling directions. + +`FixedPointGradient` differentiates the joint fixed point of the environment and +cached range finders. Charge allocations are held fixed during differentiation. +The new range finders are expressed in the input environment's enlarged spaces; +aligning the output environment back to that frame keeps the pair consistent. +The adjoint uses the same bond-gauge alignment and fixes trial-basis frames before +solving the coupled equations. `gradient_alg=nothing` instead differentiates through +all sweeps. + +## Computational and memory costs + +For square half environments of dimension ``m``, the table counts active range updates and construction of ``\rho`` per unit-cell position and CTMRG move direction, at fixed ``k,q`` in sampled sectors. +One half-environment application to ``k`` directions costs ``O(m^2 k)``; thin QR/LQ costs ``O(m k^2)``; the reduced SVD costs ``O(k^3)``. + +| Variant | Half-environment applications | QR/LQ count | Basis cache | +| --- | ---: | ---: | ---: | +| Independent | ``8(q + 1) + 2`` | ``2(q + 1)`` | ``2mk`` | +| Alternating | ``4(q + 1)`` | ``2(q + 1)`` | ``mk`` | +| Refined alternating | ``4(q + 1) + 3`` | ``2(q + 1) + 1`` | ``mk`` | + +- **Projector construction:** excluded above; standard alternating additionally applies ``R`` to ``\chi`` directions, costing ``O(m^2\chi)``. +- **Memory:** all require ``O(mk + k^2)`` temporary workspace, excluding shared corner tensors. +- **Comparison:** alternating reduces operation counts and the basis cache; all retain the same asymptotic scaling. + Iteration counts at matched accuracy can differ. diff --git a/src/Defaults.jl b/src/Defaults.jl index 78fc88ab2..1fcf1b60e 100644 --- a/src/Defaults.jl +++ b/src/Defaults.jl @@ -71,7 +71,11 @@ Module containing default algorithm parameter values and arguments. * `projector_alg=:$(Defaults.projector_alg)` : Default variant of the CTMRG projector algorithm. - `:HalfInfiniteProjector` : Projection via SVDs of half-infinite (two enlarged corners) CTMRG environments. - `:FullInfiniteProjector` : Projection via SVDs of full-infinite (all four enlarged corners) CTMRG environments. + - `:SubspaceIterationProjector` : Projection via recycled subspace iteration of the full-infinite environment. * `projector_verbosity=$(Defaults.projector_verbosity)` : Projector output information verbosity. +* `orth_alg=$(Defaults.orth_alg)` : QR/LQ algorithm configuration used by the subspace-iteration projector. +* `charge_probe=nothing` : Positive integers limit added trial multiplicities per sector on compatible-cache refreshes; `nothing` probes broadly. +* `charge_refresh=$(Defaults.charge_refresh)` : Sweep interval for exploring trial charges when `subspace` is an integer; zero disables periodic exploration. * `projector_alg_c4v=:$(Defaults.projector_alg_c4v)` : Default variant of the C4v CTMRG projector algorithm. - `:C4vEighProjector` : Projection via truncated Eigh of an enlarged corner. - `:C4vQRProjector` : Projection via QR decomposition of a column-enlarged corner. @@ -169,8 +173,10 @@ const orth_rrule_alg = :FullPullback const orth_rrule_verbosity = 0 # Projectors -const projector_alg = :HalfInfiniteProjector # ∈ {:HalfInfiniteProjector, :FullInfiniteProjector} +const projector_alg = :HalfInfiniteProjector # ∈ {:HalfInfiniteProjector, :FullInfiniteProjector, :SubspaceIterationProjector} const projector_verbosity = 0 +const orth_alg = (;) +const charge_refresh = 10 const projector_alg_c4v = :C4vEighProjector # ∈ {:C4vEighProjector, :C4vQRProjector} # Fixed-point gradient diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index b715dc2ec..6c6de439c 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -133,8 +133,11 @@ include("algorithms/contractions/correlator/pepo_1layer.jl") include("algorithms/ctmrg/sparse_environments.jl") include("algorithms/ctmrg/ctmrg.jl") include("algorithms/ctmrg/projectors/projectors.jl") +include("algorithms/ctmrg/projectors/cache.jl") include("algorithms/ctmrg/projectors/halfinfinite.jl") include("algorithms/ctmrg/projectors/fullinfinite.jl") +include("algorithms/ctmrg/projectors/subspaceiter.jl") +include("algorithms/ctmrg/projectors/subspaceiter_spaces.jl") include("algorithms/ctmrg/projectors/c4v_eigh.jl") include("algorithms/ctmrg/projectors/c4v_qr.jl") include("algorithms/ctmrg/simultaneous.jl") @@ -185,7 +188,7 @@ export initialize_ctmrg_environment, RandomInitialization, ProductStateInitialization, ApplicationInitialization, IdentityInitialization export corner, edge, setcorner!, setedge! export FixedSpaceTruncation, SiteDependentTruncation -export HalfInfiniteProjector, FullInfiniteProjector +export HalfInfiniteProjector, FullInfiniteProjector, SubspaceIterationProjector export C4vCTMRG, C4vEighProjector, C4vQRProjector export initialize_random_c4v_env, initialize_singlet_c4v_env export LocalOperator, physicalspace diff --git a/src/algorithms/ctmrg/ctmrg.jl b/src/algorithms/ctmrg/ctmrg.jl index 9263eb89f..2d76cb1ef 100644 --- a/src/algorithms/ctmrg/ctmrg.jl +++ b/src/algorithms/ctmrg/ctmrg.jl @@ -35,19 +35,29 @@ function CTMRGAlgorithm(; if decomposition_alg isa NamedTuple decomposition_alg = (; fwd_alg = decomposition_alg) end - projector_algorithm = ProjectorAlgorithm(; - alg = projector_alg, decomposition_alg, trunc, verbosity - ) + projector_alg = projector_alg isa Symbol ? (; alg = projector_alg) : projector_alg + projector_alg isa NamedTuple || + throw(ArgumentError("projector_alg must be a Symbol or NamedTuple.")) + projector_kwargs = merge((; decomposition_alg, trunc, verbosity), projector_alg) + projector_algorithm = ProjectorAlgorithm(; projector_kwargs...) return alg_type(tol, maxiter, miniter, verbosity, projector_algorithm) end """ ctmrg_iteration(network, env, alg::CTMRGAlgorithm) -> env′, info + ctmrg_iteration(network, env, alg::CTMRGAlgorithm, cache) -> env′, info, cache′ Perform a single CTMRG iteration in which all directions are being grown and renormalized. + +The four-argument form threads projector-specific runtime data between iterations. """ -function ctmrg_iteration(network, env, alg::CTMRGAlgorithm) end +function ctmrg_iteration(network, env, alg::CTMRGAlgorithm) + cache = initialize_projector_cache(network, env, alg) + env′, info, = ctmrg_iteration(network, env, alg, cache) + return env′, info +end +function ctmrg_iteration(network, env, alg::CTMRGAlgorithm, cache) end """ leading_boundary(env₀, network; kwargs...) -> env, info @@ -84,11 +94,16 @@ supplied via the keyword arguments or directly as an [`CTMRGAlgorithm`](@ref) st - `:truncrank` : Additionally supply truncation dimension `η`; truncate to the maximal virtual dimension of `η` - `:truncspace` : Additionally supply truncation space `η`; truncate according to the supplied vector space - `:trunctol` : Additionally supply singular value cutoff `η`; truncate such that every retained singular value is larger than `η` -* `projector_alg::Symbol=:$(Defaults.projector_alg)` : Variant of the projector algorithm. See also [`ProjectorAlgorithm`](@ref). +* `projector_alg::Union{Symbol,NamedTuple}=:$(Defaults.projector_alg)` : Variant of the projector algorithm, optionally supplied as a `NamedTuple` containing projector-specific keyword arguments. See also [`ProjectorAlgorithm`](@ref). - `:HalfInfiniteProjector` : Projection via SVDs of half-infinite (two enlarged corners) CTMRG environments. - `:FullInfiniteProjector` : Projection via SVDs of full-infinite (all four enlarged corners) CTMRG environments. + - `:SubspaceIterationProjector` : Stateful projection via recycled subspace iteration of the full-infinite environment, with `subspace` specified as a symmetry-resolved space or an integer trial-dimension budget. Supports fixed-point differentiation or backpropagation through the CTMRG sweeps. - `:C4vEighProjector` : Projection via `eigh` of the Hermitian enlarged corner, works only for [`C4vCTMRG`](@ref). - `:C4vQRProjector` : Projection via QR decomposition of the lower-rank column-enlarged corner, works only for [`C4vCTMRG`](@ref). + Its `update` selects `:independent` (default, independent updates), `:alternating` (alternating QR/LQ), or `:alternating_refined` (with a final left-basis refinement). + For example, use `projector_alg=(; alg=:SubspaceIterationProjector, subspace=ℂ^10, update=:alternating, n_subspace_iter=1)` for active alternating refinement with one extra power step. + Use `subspace=24, trunc=truncrank(16)` to adapt both trial and retained charge distributions; `charge_refresh` controls the exploration interval. + `rotate_basis=true` expresses all trial directions in reduced singular-vector coordinates. * `decomposition_alg::Union{NamedTuple,<:SVDAdjoint,<:EighAdjoint,<:OrthAdjoint}` : Tensor decomposition algorithm used for computing projectors. When specified as a `NamedTuple`, the settings are passed as the forward algorithm to the appropriate decomposition for the given projector algorithm. For information on which forward algorithms are available, and how to specify them, see [`SVDAdjoint`](@ref), [`EighAdjoint`](@ref) and [`OrthAdjoint`](@ref). @@ -108,6 +123,14 @@ end function leading_boundary( env₀::AbstractCTMRGEnv, network::InfiniteSquareNetwork, alg::CTMRGAlgorithm ) + env, info, = _leading_boundary_with_cache(env₀, network, alg) + return env, info +end + +"""Contract the boundary and retain its final projector cache for fixed-point differentiation.""" +function _leading_boundary_with_cache( + env₀::AbstractCTMRGEnv, network::InfiniteSquareNetwork, alg::CTMRGAlgorithm + ) check_input(leading_boundary, network, env₀, alg) log = ignore_derivatives(() -> MPSKit.IterLog("CTMRG")) return with_verbosity(; alg.verbosity) do @@ -118,12 +141,25 @@ function leading_boundary( η = one(real(scalartype(network))) ctmrg_loginit!(log, η, network, env₀) local info_iter + projector_cache = initialize_projector_cache(network, env, alg) converged = false + align_gauge = false for iter in 1:(alg.maxiter) - env, info_iter = ctmrg_iteration(network, env, alg) + previous = env + env, info_iter, projector_cache = ctmrg_iteration(network, env, alg, projector_cache) η, CS, TS = calc_convergence(env, CS, TS, alg) + ready = iter ≥ alg.miniter && _projector_convergence_ready(projector_cache, alg.projector_alg) + if projector_cache isa SubspaceIterationCache && + (align_gauge || (η ≤ alg.tol && ready)) + align_gauge = true + signs = ignore_derivatives(() -> first(compute_gauge_fix_gauge(env, previous, _scrambling_env_gauge(alg)))) + # New rangefinders live in the input frame; align the output to that frame. + env = fix_global_phases(fix_relative_phases(env, signs), previous) + η = max(η, ignore_derivatives(() -> norm(env - previous) / norm(previous))) + end + η = max(η, get(info_iter.contraction_metrics, :subspace_error, zero(η))) - if η ≤ alg.tol && iter ≥ alg.miniter + if η ≤ alg.tol && ready ctmrg_logfinish!(log, iter, η, network, env) converged = true break @@ -139,7 +175,7 @@ function leading_boundary( convergence_error = η, contraction_metrics = info_iter.contraction_metrics, ) - return env, info + return env, info, projector_cache end end function leading_boundary(env₀::AbstractCTMRGEnv, state, args...; kwargs...) diff --git a/src/algorithms/ctmrg/projectors/cache.jl b/src/algorithms/ctmrg/projectors/cache.jl new file mode 100644 index 000000000..9eb10f2e1 --- /dev/null +++ b/src/algorithms/ctmrg/projectors/cache.jl @@ -0,0 +1,89 @@ +"""Abstract runtime cache threaded between projector iterations.""" +abstract type ProjectorCache end + +"""Empty runtime cache used by projector algorithms without recycled data.""" +struct EmptyProjectorCache <: ProjectorCache end + +"""Allow convergence unless an adaptive projector still needs a charge-allocation check.""" +_projector_convergence_ready(cache, alg) = true + +""" +Runtime rangefinder cache for SI-CTMRG iterations. +Both `X` and `Y` embed the trial space into the full right and left spaces, respectively. +""" +struct SubspaceIterationCache{X, Y} <: ProjectorCache + iteration::Int + X::X + Y::Y +end + +"""Copy all SI tensor grids while preserving cache metadata.""" +function Base.copy(cache::SubspaceIterationCache) + return SubspaceIterationCache( + cache.iteration, + map(_copy_subspace_grid, (cache.X, cache.Y))..., + ) +end + +"""Initialize an empty runtime projector cache for a fresh boundary contraction.""" +initialize_projector_cache(network, env, alg::CTMRGAlgorithm) = EmptyProjectorCache() + +"""Thread an unchanged empty projector cache through a three-argument CTMRG iteration.""" +function ctmrg_iteration(network, env, alg::CTMRGAlgorithm, cache::EmptyProjectorCache) + env′, info = ctmrg_iteration(network, env, alg) + return env′, info, cache +end + +"""Report whether cached SI rangefinders match the current map spaces and scalar type.""" +function _compatible_rangefinders( + X, Y, halfinf_left, halfinf_right, + subspace::ElementarySpace, require_left::Bool = true, + ) + isnothing(X) && return false + right = scalartype(X) == scalartype(halfinf_left) && + space(X) == (domain(halfinf_right) ← subspace) + left = !require_left || ( + !isnothing(Y) && + scalartype(Y) == scalartype(halfinf_left) && + space(Y) == (codomain(halfinf_left) ← subspace) + ) + return right && left +end + +"""Copy an optional cache grid while leaving absent state unallocated.""" +_copy_subspace_grid(grid) = isnothing(grid) ? nothing : copy(grid) + +"""Index a trial grid, returning no left basis when that grid is absent.""" +_subspace_grid_entry(grid, indices...) = isnothing(grid) ? nothing : grid[indices...] + +"""Check grid layout, scalar type, and state required by the configured update.""" +function _compatible_subspace_cache(cache, env, alg) + return all(g -> isnothing(g) || size(g) == size(env), (cache.X, cache.Y)) && + scalartype(first(cache.X)) == scalartype(env) && + (!_retain_subspace_left(alg) || !isnothing(cache.Y)) +end + +"""Rotate a directional tensor grid without mutating arrays needed by backpropagation.""" +_rotate_subspace_grid(grid) = isnothing(grid) ? nothing : reverse(permutedims(circshift(grid, (-1, 0, 0)), (1, 3, 2)); dims = 2) + +"""Rotate all directional SI cache arrays counter-clockwise.""" +function Base.rotl90(cache::SubspaceIterationCache) + X′, Y′ = map(_rotate_subspace_grid, (cache.X, cache.Y)) + return SubspaceIterationCache(cache.iteration, X′, Y′) +end + +"""RMS sine of trial-space changes, computed without subtracting nearly equal norms.""" +_subspace_error(::Nothing, ::Nothing) = 0.0 +function _subspace_error(old, new) + size(old) == size(new) || return Inf + return maximum(zip(old, new)) do (a, b) + space(a) == space(b) || return Inf + norm(b - a * (a' * b)) / sqrt(dim(domain(b))) + end +end +@non_differentiable _subspace_error(args...) + +function _update_subspace_cache(cache, X, Y) + error = max(_subspace_error(cache.X, X), _subspace_error(cache.Y, Y)) + return SubspaceIterationCache(cache.iteration + 1, X, Y), error +end diff --git a/src/algorithms/ctmrg/projectors/projectors.jl b/src/algorithms/ctmrg/projectors/projectors.jl index 8d1f6c637..e58defd92 100644 --- a/src/algorithms/ctmrg/projectors/projectors.jl +++ b/src/algorithms/ctmrg/projectors/projectors.jl @@ -15,6 +15,14 @@ Keyword argument parser returning the appropriate `ProjectorAlgorithm` algorithm function ProjectorAlgorithm(; alg = Defaults.projector_alg, decomposition_alg = (;), + orth_alg = Defaults.orth_alg, + subspace = nothing, + charge_refresh::Int = Defaults.charge_refresh, + charge_probe = nothing, + update::Symbol = :independent, + n_subspace_iter::Int = 0, + rotate_basis::Bool = true, + rng = Random.default_rng(), trunc = (;), verbosity = Defaults.projector_verbosity, ) @@ -25,7 +33,7 @@ function ProjectorAlgorithm(; # parse SVD forward & rrule algorithm - decomposition_algorithm = if alg in [:HalfInfiniteProjector, :FullInfiniteProjector] + decomposition_algorithm = if alg in [:HalfInfiniteProjector, :FullInfiniteProjector, :SubspaceIterationProjector] _alg_or_nt(SVDAdjoint, decomposition_alg) elseif alg in [:C4vEighProjector] _alg_or_nt(EighAdjoint, decomposition_alg) @@ -47,6 +55,42 @@ function ProjectorAlgorithm(; throw(ArgumentError("unknown trunc $trunc")) end + if alg == :SubspaceIterationProjector + truncation_strategy isa Union{FixedSpaceTruncation, TruncationSpace, MatrixAlgebraKit.TruncationByOrder} || + throw(ArgumentError("SubspaceIterationProjector supports fixed-space, truncspace, or truncrank truncation.")) + subspace isa ElementarySpace || (subspace isa Integer && !(subspace isa Bool) && 0 < subspace <= typemax(Int)) || + throw(ArgumentError("subspace must be an elementary space or a positive integer dimension.")) + charge_refresh >= 0 || throw(ArgumentError("charge_refresh must be nonnegative.")) + isnothing(charge_probe) || (charge_probe isa Int && charge_probe > 0) || + throw(ArgumentError("charge_probe must be nothing or a positive Int.")) + if truncation_strategy isa Union{TruncationSpace, MatrixAlgebraKit.TruncationByOrder} + truncation_strategy.by === abs && truncation_strategy.rev || + throw(ArgumentError("SubspaceIterationProjector retains the largest singular values.")) + end + if truncation_strategy isa MatrixAlgebraKit.TruncationByOrder + truncation_strategy.howmany > 0 || throw(ArgumentError("The retained rank must be positive.")) + subspace isa Integer && subspace < truncation_strategy.howmany && + throw(ArgumentError("The trial dimension must be at least the retained rank.")) + end + update in (:independent, :alternating, :alternating_refined) || + throw(ArgumentError("unknown subspace update: $update")) + n_subspace_iter >= 0 || throw(ArgumentError("n_subspace_iter must be nonnegative.")) + orth_alg = _alg_or_nt(OrthAdjoint, orth_alg) + return alg_type( + decomposition_algorithm, + truncation_strategy, + orth_alg, + subspace, + charge_refresh, + charge_probe, + verbosity, + update, + n_subspace_iter, + rotate_basis, + rng, + ) + end + return alg_type(decomposition_algorithm, truncation_strategy, verbosity) end @@ -66,6 +110,54 @@ function truncation_strategy(alg::ProjectorAlgorithm, edge) end end +""" +Collect the cyclic coordinate layout and target truncation for a projector. + +``N`` is the number of enlarged corners involved in the construction of each projector. +- ``N = 2`` for half-infinite environment. +- ``N = 4`` for full-infinite environment. +""" +function _enlarged_corner_layout( + coordinate::NTuple{3, Int}, env::CTMRGEnv, + alg::ProjectorAlgorithm, ::Val{N}, + ) where {N} + rowsize, colsize = size(env)[2:3] + coordinate2 = _next_coordinate(coordinate, rowsize, colsize) + trunc = truncation_strategy(alg, env.edges[coordinate[1], coordinate2[2:3]...]) + coordinates = if N == 2 + (coordinate, coordinate2) + elseif N == 4 + coordinate3 = _next_coordinate(coordinate2, rowsize, colsize) + coordinate4 = _next_coordinate(coordinate3, rowsize, colsize) + (coordinate4, coordinate, coordinate2, coordinate3) + else + throw(ArgumentError("projectors require either two or four enlarged corners")) + end + return coordinates, trunc +end + +"""Collect precomputed enlarged corners and the target truncation for simultaneous CTMRG.""" +function enlarged_corner_inputs( + corner_grid::AbstractArray{E, 3}, coordinate::NTuple{3, Int}, + env::CTMRGEnv, alg::ProjectorAlgorithm, count::Val{N}, + ) where {E, N} + coordinates, trunc = _enlarged_corner_layout(coordinate, env, alg, count) + enlarged_corners = map(coordinate -> corner_grid[coordinate...], coordinates) + return enlarged_corners, trunc +end + +"""Construct enlarged corners and collect the target truncation for sequential CTMRG.""" +function enlarged_corner_inputs( + network::InfiniteSquareNetwork, coordinate::NTuple{3, Int}, + env::CTMRGEnv, alg::ProjectorAlgorithm, count::Val{N}, + ) where {N} + coordinates, trunc = _enlarged_corner_layout(coordinate, env, alg, count) + enlarged_corners = map(coordinates) do coordinate + return TensorMap(EnlargedCorner(network, env, coordinate)) + end + return enlarged_corners, trunc +end + """ _set_truncation(alg::ProjectorAlgorithm, trunc::TruncationStrategy) diff --git a/src/algorithms/ctmrg/projectors/subspaceiter.jl b/src/algorithms/ctmrg/projectors/subspaceiter.jl new file mode 100644 index 000000000..bca8ca89d --- /dev/null +++ b/src/algorithms/ctmrg/projectors/subspaceiter.jl @@ -0,0 +1,306 @@ +""" +$(TYPEDEF) + +Recycled subspace approximation of the full-infinite CTMRG environment. +Supply `subspace` as an elementary space with explicit charge multiplicities, or as a positive integer specifying the total trial dimension. +An integer budget selects its charge distribution from the sampled singular spectrum, preserving complete multiplets. +`trunc` independently selects the retained environment: [`FixedSpaceTruncation`](@ref), `truncspace`, or `truncrank`. +With fixed target sectors, the trial allocation reserves their required multiplicities. +`charge_probe=nothing` uses broad exploration; a positive integer limits added multiplicities per sector on compatible-cache refreshes. +Initialization and incompatible-cache recovery remain broad. +Integer allocations are explored initially and every `charge_refresh=$(Defaults.charge_refresh)` sweeps. +Set `charge_refresh=0` to disable periodic exploration; incompatible cached spaces still trigger exploration. +Refreshes select trial directions in reduced SVD coordinates even when ordinary basis rotations are disabled. +See [Subspace-iteration CTMRG projectors](@ref subspace_iteration) for update equations and convergence. +Both trial bases embed the reduced space into the full environment spaces: ``X^{\\dagger} X = Y^{\\dagger} Y = I`` and ``ρ = Y^{\\dagger} M X``. + +The `update` selector chooses `:independent` (default, the paper's two-factorization independent update), `:alternating` (alternating QR/LQ with a reused reduced factor), or `:alternating_refined` (three factorizations with a final left-basis refresh). +`n_subspace_iter=q` adds `q` independent updates, or `q` stabilized right-basis power steps before the mandatory alternating update. +`rotate_basis=true` expresses the complete trial bases in the reduced SVD coordinates, retaining all oversampling directions. + +[`FixedPointGradient`](@ref) differentiates the joint environment/rangefinder fixed point, holding only the charge allocation fixed. +Use `gradient_alg=nothing` in [`PEPSOptimize`](@ref) to differentiate through all sweeps and recycled bases. + +## Fields + +$(TYPEDFIELDS) + +## Constructors + + SubspaceIterationProjector(; subspace, kwargs...) + +Defaults are `update=:independent`, `n_subspace_iter=0`, `rotate_basis=true`, `rng=Random.default_rng()`, `decomposition_alg=SVDAdjoint()`, `orth_alg=OrthAdjoint()`, and `trunc=FixedSpaceTruncation()`. +""" +struct SubspaceIterationProjector{M, S <: SVDAdjoint, T, Q <: OrthAdjoint, E <: Union{ElementarySpace, Integer}, R} <: ProjectorAlgorithm + decomposition_alg::S + trunc::T + orth_alg::Q + subspace::E + charge_refresh::Int + charge_probe::Union{Nothing, Int} + verbosity::Int + update::Symbol + n_subspace_iter::Int + rotate_basis::Bool + rng::R +end + +function SubspaceIterationProjector(; kwargs...) + return ProjectorAlgorithm(; alg = :SubspaceIterationProjector, trunc = FixedSpaceTruncation(), kwargs...) +end + +"""Specialize the range kernel on the selected update without additional public algorithm types.""" +function SubspaceIterationProjector(svd, trunc, orth, space, refresh, probe, verbosity, update::Symbol, q::Int, rotate::Bool, rng) + return SubspaceIterationProjector{update, typeof(svd), typeof(trunc), typeof(orth), typeof(space), typeof(rng)}( + svd, trunc, orth, space, refresh, probe, verbosity, update, q, rotate, rng + ) +end + +PROJECTOR_SYMBOLS[:SubspaceIterationProjector] = SubspaceIterationProjector + +_retain_subspace_left(alg::SubspaceIterationProjector{M}) where {M} = M == :independent + +"""Apply a half environment through its enlarged corners without forming their product.""" +struct SIHalfEnvironment{L, R, S} + left::L + right::R + scale::S +end + +"""Normalize an implicit half without copying its enlarged corners.""" +SIHalfEnvironment(left, right) = SIHalfEnvironment(left, right, inv(norm(left) * norm(right))) + +TensorKit.domain(h::SIHalfEnvironment) = domain(h.right) +TensorKit.codomain(h::SIHalfEnvironment) = codomain(h.left) +TensorKit.scalartype(h::SIHalfEnvironment) = scalartype(h.left) +Base.adjoint(h::SIHalfEnvironment) = SIHalfEnvironment(h.right', h.left', h.scale) +Base.:*(h::SIHalfEnvironment, x::AbstractTensorMap) = scale!!(h.left * (h.right * x), h.scale) +Base.:*(x::AbstractTensorMap, h::SIHalfEnvironment) = scale!!((x * h.left) * h.right, h.scale) +Base.:*(h::SIHalfEnvironment{<:Any, <:Any, Nothing}, x::AbstractTensorMap) = h.left * (h.right * x) +Base.:*(x::AbstractTensorMap, h::SIHalfEnvironment{<:Any, <:Any, Nothing}) = (x * h.left) * h.right + +"""Normalize unscaled corner halves while preserving explicitly supplied maps.""" +_normalize_subspace_half(h) = h +_normalize_subspace_half(h::SIHalfEnvironment{<:Any, <:Any, Nothing}) = SIHalfEnvironment(h.left, h.right) + +"""Check compatible network, environment, and explicitly specified trial sectors.""" +function _check_subspace_iteration_input(network, env, alg::SubspaceIterationProjector) + sector = sectortype(spacetype(network)) + sectortype(spacetype(env)) === sector && + (!(alg.subspace isa ElementarySpace) || sectortype(alg.subspace) === sector) || + throw(ArgumentError("Network, environment, and subspace must use the same symmetry sectors.")) + return nothing +end + +"""Select trial multiplicities from retained charges or an explicit space, clipped to both map dimensions.""" +function _subspace_iteration_space(left, right, alg, trunc) + available = infimum(fuse(left), fuse(right)) + _adaptive_subspace(alg) && return _subspace_probe_space(available, alg, trunc) + K = infimum(alg.subspace, available) + !(trunc isa TruncationSpace) || infimum(trunc.space, K) == trunc.space || + throw(ArgumentError("Subspace must contain the target truncation space.")) + return K +end + +"""Initialize right and required left isometries from the trial space into the full spaces.""" +function rangefinder_initialize(::Type{T}, left::ProductSpace, right::ProductSpace, K::ElementarySpace, alg::SubspaceIterationProjector{M}) where {T, M} + X, _ = left_orth!(randn(alg.rng, T, right ← K), alg.orth_alg) + Y = if M == :independent + last(right_orth!(randn(alg.rng, T, K ← left), alg.orth_alg))' + else + nothing + end + return X, Y +end + +"""Initialize trial bases from factorized half-environment spaces.""" +function rangefinder_initialize(L, R, K::ElementarySpace, alg::SubspaceIterationProjector) + return rangefinder_initialize(scalartype(L), codomain(L), domain(R), K, alg) +end + +"""Initialize one directional cache entry without forming enlarged corners.""" +function _initialize_subspace_iteration_entry(network::InfiniteSquareNetwork, env::CTMRGEnv, coordinate::NTuple{3, Int}, alg::SubspaceIterationProjector) + coordinates, trunc = _enlarged_corner_layout(coordinate, env, alg, Val(4)) + left = codomain(EnlargedCorner(network, env, first(coordinates))) + right = domain(EnlargedCorner(network, env, last(coordinates))) + K = _subspace_iteration_space(left, right, alg, trunc) + T = scalartype(network) + X, Y = rangefinder_initialize(T, left, right, K, alg) + return (; X, Y) +end + +"""Initialize only the trial grids required by the selected update.""" +function initialize_subspace_iteration_cache(network::InfiniteSquareNetwork, env::CTMRGEnv, alg::SubspaceIterationProjector) + entries = map(c -> _initialize_subspace_iteration_entry(network, env, c, alg), eachcoordinate(env, 1:4)) + Y = _retain_subspace_left(alg) ? map(x -> x.Y, entries) : nothing + return SubspaceIterationCache(0, map(x -> x.X, entries), Y) +end + +@non_differentiable initialize_subspace_iteration_cache(args...) + +"""Construct a sampled left range and exact bases for complete left sectors.""" +function _subspace_alternating_range(L, R, X, K, alg) + W = ignore_derivatives(() -> fuse(codomain(L))) + sampled = ignore_derivatives(() -> Vect[sectortype(K)](s => dim(K, s) for s in sectors(K) if dim(K, s) < dim(W, s))) + dim(sampled) == 0 && return ignore_derivatives(() -> isometry(scalartype(L), codomain(L) ← K)), false + E = sampled == K ? nothing : ignore_derivatives(() -> isometry(scalartype(X), K ← sampled)) + Xsampled = isnothing(E) ? X : X * E + Q, _ = left_orth!(L * (R * Xsampled), alg.orth_alg) + isnothing(E) && return Q, true + Q0, E0 = ignore_derivatives() do + exact = Vect[sectortype(K)](s => dim(K, s) for s in sectors(K) if dim(sampled, s) == 0) + isometry(scalartype(L), codomain(L) ← exact), isometry(scalartype(L), K ← exact) + end + return Q0 * E0' + Q * E', true +end + +""" +Update the trial bases and construct `core`, the reduced matrix ``ρ = Y^{\\dagger} L R X``. +Use independent normal-map updates or alternating QR steps with a final LQ or refined QR factor reuse. +""" +function rangefinder_update(X::AbstractTensorMap, Y, L, R, K, alg::SubspaceIterationProjector{M}) where {M} + if M == :independent + if ignore_derivatives(() -> fuse(domain(X)) == fuse(domain(R)) == fuse(codomain(L))) + X, Y = ignore_derivatives() do + isometry(storagetype(X), domain(R) ← K), isometry(storagetype(X), codomain(L) ← K) + end + else + for _ in 0:alg.n_subspace_iter + MX = L * (R * X) + Xnew = R' * (L' * MX) + # Adjoint LQ evaluates leftorth(M M'Y) in the cached row layout. + YadjM = (Y' * L) * R + _, Yrow = right_orth!((YadjM * R') * L', alg.orth_alg) + Y = Yrow' + X, _ = left_orth!(Xnew, alg.orth_alg) + end + end + right = R * X + left = Y' * L + core = left * right + else + # Initialize Q, then perform q right/left update pairs. + Q, sampled = if M == :alternating + _subspace_alternating_range(L, R, X, K, alg) + else + Q, _ = left_orth!(L * (R * X), alg.orth_alg) + Q, true + end + for _ in 1:alg.n_subspace_iter + sampled || break + X, _ = left_orth!(R' * (L' * Q), alg.orth_alg) + Q, _ = if M == :alternating + _subspace_alternating_range(L, R, X, K, alg) + else + left_orth!(L * (R * X), alg.orth_alg) + end + end + if M == :alternating + # The LQ of Q'M is the adjoint of the final QR; its triangular factor is ρ. + Y = Q + left = Y' * L + B = left * R + core, Z = sampled ? right_orth!(B, alg.orth_alg) : (B, nothing) + X = isnothing(Z) ? nothing : Z' + right = nothing + else + X, _ = left_orth!(R' * (L' * Q), alg.orth_alg) + right = R * X + Y, core = left_orth!(L * right, alg.orth_alg) + left = Y' * L + end + end + return (; X, Y, left, right, core) +end + +""" +Apply MatrixAlgebraKit's SVD phase convention with a differentiable correction. + +The forward SVD fixes this gauge by default, but its pullback assumes gauge-invariant cotangents. +Keep the phase correction in the AD graph to support phase-sensitive singular-vector outputs. +""" +function _subspace_svd_phase_correction(U::AbstractTensorMap, V::AbstractTensorMap) + layout = ignore_derivatives(() -> collect(values(TensorKit.blockstructure(space(U))))) + phase_blocks = map(layout) do (dims, inds) + u = reshape(view(U.data, inds), dims) + pivots = ignore_derivatives(() -> [last(findmax(abs, col)) for col in eachcol(u)]) + return vec(diagm([sign(u[pivots[j], j]) for j in axes(u, 2)])) + end + phases = TensorMap{scalartype(U)}(reduce(vcat, phase_blocks), domain(U) ← domain(U)) + return U * phases', phases * V +end + +""" +Factor `core`, the reduced matrix ``ρ = Y^{\\dagger} M X``, retaining complete singular directions when recycling needs them. +The `:alternating` update normalizes this matrix before its SVD. +""" +function _subspace_reduced_svd(core, trunc, alg, complete::Bool) + svd_alg = decomposition_algorithm(alg) + fwd = svd_alg.fwd_alg + fwd = fwd isa TruncatedAlgorithm ? fwd.alg : fwd + # TODO: Unify alternating core normalization with the other updates' half-environment scaling; + # verify projector scaling, convergence, and gradients before changing the conventions. + rho = alg.update == :alternating ? core / norm(core) : core + if complete + svd_alg = @set svd_alg.fwd_alg = TruncatedAlgorithm(fwd, notrunc()) + u, s, v = svd_trunc_no_error!(rho, svd_alg) + u, v = _subspace_svd_phase_correction(u, v) + E = ignore_derivatives(() -> _subspace_truncation_isometry(s, trunc, scalartype(core))) + uχ, sχ, vχ = u * E, DiagonalTensorMap(E' * s * E), E' * v + error = ignore_derivatives(() -> norm(s - E * sχ * E') / norm(sχ)) + return (; u, s, v, uχ, sχ, vχ, error) + end + svd_alg = @set svd_alg.fwd_alg = TruncatedAlgorithm(fwd, trunc) + u, s, v, error = svd_trunc!(rho, svd_alg) + u, v = _subspace_svd_phase_correction(u, v) + error = ignore_derivatives(() -> error / norm(s)) + return (; u, s, v, uχ = u, sχ = s, vχ = v, error) +end + +""" +Compute shared projectors, singular frames, and the minimal next-sweep trial state. +The internal `core` is the reduced matrix ``ρ = Y^{\\dagger} M X`` with ``M = L R``, obtained by contraction or a reused QR/LQ factor. +""" +function _subspace_iteration_decomposition(L, R, alg::SubspaceIterationProjector{M}, ranges) where {M} + trunc = ranges.trunc + K, refresh = ignore_derivatives() do + return _adaptive_subspace(alg) ? _adaptive_subspace_space(L, R, alg, ranges) : + (_subspace_iteration_space(codomain(L), domain(R), alg, trunc), false) + end + compatible = _compatible_rangefinders(ranges.X, ranges.Y, L, R, K, _retain_subspace_left(alg)) + X, Y = if refresh + X = _resize_subspace_basis(ranges.X, domain(R), K, alg, scalartype(L)) + Y = _retain_subspace_left(alg) ? copy(_resize_subspace_basis(ranges.Y, codomain(L), K, alg, scalartype(L))')' : nothing + X, Y + else + compatible ? (ranges.X, ranges.Y) : ignore_derivatives(() -> rangefinder_initialize(L, R, K, alg)) + end + if M != :alternating + L, R = _normalize_subspace_half(L), _normalize_subspace_half(R) + end + state = rangefinder_update(X, Y, L, R, K, alg) + full = _subspace_reduced_svd(state.core, trunc, alg, refresh || alg.rotate_basis || isnothing(state.X)) + U = state.Y * full.uχ + V = isnothing(state.X) ? full.vχ : full.vχ * state.X' + isqS = sdiag_pow(full.sχ, -0.5) + P_left = isnothing(state.right) ? R * (V' * isqS) : state.right * (full.vχ' * isqS) + P_right = (isqS * full.uχ') * state.left + rotate = refresh || alg.rotate_basis + X = isnothing(state.X) ? copy(full.v') : (rotate ? state.X * full.v' : copy(state.X)) + # Keep the left cache as a lazy adjoint of the existing row-oriented contraction. + Y = _retain_subspace_left(alg) ? (rotate ? full.u' * state.Y' : copy(state.Y'))' : nothing + if refresh + E = ignore_derivatives(() -> _adaptive_subspace_isometry(full.s, alg, trunc, scalartype(state.core))) + X = X * E + Y = isnothing(Y) ? nothing : (E' * Y')' + end + return (; U, S = full.sχ, V, X, Y, truncation_error = full.error, projectors = (P_left, P_right)) +end + +"""Build factorized half environments using the established scaling for each update.""" +function compute_projector(corners, alg::SubspaceIterationProjector, ranges) + L = SIHalfEnvironment(corners[1], corners[2], nothing) + R = SIHalfEnvironment(corners[3], corners[4], nothing) + info = _subspace_iteration_decomposition(L, R, alg, ranges) + return info.projectors, info +end diff --git a/src/algorithms/ctmrg/projectors/subspaceiter_spaces.jl b/src/algorithms/ctmrg/projectors/subspaceiter_spaces.jl new file mode 100644 index 000000000..df15c07c2 --- /dev/null +++ b/src/algorithms/ctmrg/projectors/subspaceiter_spaces.jl @@ -0,0 +1,87 @@ +"""Identify trial spaces whose charge allocation is selected from a sampled spectrum.""" +_adaptive_subspace(alg::SubspaceIterationProjector) = alg.subspace isa Integer + +"""Probe shared sectors broadly or add bounded multiplicities to a compatible allocation.""" +function _subspace_probe_space(available, alg, trunc, current = nothing) + budget = Int(alg.subspace) + required = trunc isa TruncationSpace ? trunc.space : zero(available) + infimum(required, available) == required || throw(ArgumentError("The environment cannot contain the target truncation space.")) + dim(required) <= budget || throw(ArgumentError("The trial dimension must contain the target truncation space.")) + K = typeof(available)( + s => begin + cap = min(dim(available, s), floor(Int, budget / dim(s))) + if isnothing(current) || isnothing(alg.charge_probe) + cap + else + n = dim(current, s) + max(dim(required, s), n + min(alg.charge_probe, cap - n)) + end + end for s in sectors(available) + ) + dim(K) > 0 || throw(ArgumentError("The trial budget cannot fit an available symmetry multiplet.")) + return K +end + +"""Refresh shared charges periodically.""" +_subspace_charge_refresh(cache, alg) = _adaptive_subspace(alg) && (cache.iteration == 0 || (alg.charge_refresh > 0 && cache.iteration % alg.charge_refresh == 0)) + +"""Confirm convergence on a charge-refresh sweep when periodic exploration is enabled.""" +function _projector_convergence_ready(cache::SubspaceIterationCache, alg::SubspaceIterationProjector) + return !_adaptive_subspace(alg) || alg.charge_refresh == 0 || (cache.iteration - 1) % alg.charge_refresh == 0 +end + +"""Reuse a valid allocation or probe every shared charge up to the total trial budget per sector.""" +function _adaptive_subspace_space(L, R, alg, ranges) + available = infimum(fuse(codomain(L)), fuse(domain(R))) + budget = Int(alg.subspace) + K = isnothing(ranges.X) ? zero(available) : fuse(domain(ranges.X)) + valid = _compatible_rangefinders(ranges.X, ranges.Y, L, R, K, _retain_subspace_left(alg)) && + infimum(K, available) == K && dim(K) <= budget + if ranges.trunc isa TruncationSpace + valid &= infimum(K, ranges.trunc.space) == ranges.trunc.space + end + refresh = !valid || get(ranges, :refresh, true) + if refresh + K = _subspace_probe_space(available, alg, ranges.trunc, valid ? K : nothing) + end + return K, refresh +end + +"""Preserve compatible cached columns and add orthogonal Gaussian directions in enlarged sectors.""" +function _resize_subspace_basis(X, output, K, alg, ::Type{T}) where {T} + if isnothing(X) || scalartype(X) != T || codomain(X) != output + return ignore_derivatives(() -> first(left_orth!(randn(alg.rng, T, output ← K), alg.orth_alg))) + end + fuse(domain(X)) == K && return X + old, keep, added, extra = ignore_derivatives() do + common = infimum(fuse(domain(X)), K) + extra = K ⊖ common + isometry(T, domain(X) ← common), isometry(T, K ← common), left_null(isometry(T, K ← common)), extra + end + U = X * old + dim(extra) == 0 && return U * keep' + Z = ignore_derivatives(() -> randn(alg.rng, T, output ← extra)) + Q, _ = left_orth!(Z - U * (U' * Z), alg.orth_alg) + return U * keep' + Q * added' +end + +"""Select reduced singular directions with fixed sector counts or a global truncation strategy.""" +function _subspace_truncation_isometry(s, trunc, ::Type{T}) where {T} + W = space(s, 1) + trunc isa TruncationSpace && return isometry(T, W ← trunc.space) + indices = MatrixAlgebraKit.findtruncated_svd(diagview(s), trunc) + target = TensorKit.Factorizations.truncate_space(W, indices) + dim(target) > 0 || throw(ArgumentError("Truncation retained no symmetry multiplet.")) + return TensorKit.Factorizations.truncate_domain!(zeros(T, W ← target), id(T, W), indices) +end + +"""Keep the strongest sampled directions within the budget, reserving fixed target multiplicities.""" +function _adaptive_subspace_isometry(s, alg, trunc, ::Type{T}) where {T} + scores = trunc isa TruncationSpace ? copy(s) : s + if trunc isa TruncationSpace + for c in sectors(trunc.space) + diagview(block(scores, c))[1:dim(trunc.space, c)] .= Inf + end + end + return _subspace_truncation_isometry(scores, truncrank(Int(alg.subspace)), T) +end diff --git a/src/algorithms/ctmrg/sequential.jl b/src/algorithms/ctmrg/sequential.jl index 3ce75c4c1..bb74060fd 100644 --- a/src/algorithms/ctmrg/sequential.jl +++ b/src/algorithms/ctmrg/sequential.jl @@ -22,7 +22,7 @@ For a full description, see [`leading_boundary`](@ref). The supported keywords a * `verbosity::Int=$(Defaults.ctmrg_verbosity)` * `trunc::Union{TruncationStrategy,NamedTuple}=(; alg::Symbol=:$(Defaults.trunc))` * `decomposition_alg::Union{<:SVDAdjoint,NamedTuple}` -* `projector_alg::Symbol=:$(Defaults.projector_alg)` +* `projector_alg::Union{Symbol,NamedTuple}=:$(Defaults.projector_alg)` """ struct SequentialCTMRG{P <: ProjectorAlgorithm} <: CTMRGAlgorithm tol::Float64 @@ -37,12 +37,25 @@ end CTMRG_SYMBOLS[:SequentialCTMRG] = SequentialCTMRG +function initialize_projector_cache( + network, env, alg::Union{SimultaneousCTMRG{<:SubspaceIterationProjector}, SequentialCTMRG{<:SubspaceIterationProjector}} + ) + return initialize_subspace_iteration_cache(network, env, alg.projector_alg) +end + +function check_input( + ::typeof(leading_boundary), network, env, + alg::Union{SimultaneousCTMRG{<:SubspaceIterationProjector}, SequentialCTMRG{<:SubspaceIterationProjector}}, + ) + return _check_subspace_iteration_input(network, env, alg.projector_alg) +end + """ ctmrg_leftmove(col::Int, network, env::CTMRGEnv, alg::SequentialCTMRG) Perform sequential CTMRG left move on the `col`-th column. """ -function ctmrg_leftmove(col::Int, network, env::CTMRGEnv, alg::SequentialCTMRG) +function ctmrg_leftmove(col::Int, network, env::CTMRGEnv, alg::SequentialCTMRG, args...) #= ----> left move C1 ← T1 ← r-1 @@ -52,12 +65,15 @@ function ctmrg_leftmove(col::Int, network, env::CTMRGEnv, alg::SequentialCTMRG) C4 → T3 → r+1 c-1 c =# - projectors, info = sequential_projectors(col, network, env, alg.projector_alg) + projectors, info = sequential_projectors(col, network, env, alg.projector_alg, args...) env = renormalize_sequentially(col, projectors, network, env) return env, info end -function ctmrg_iteration(network, env::CTMRGEnv, alg::SequentialCTMRG) +function ctmrg_iteration( + network, env::CTMRGEnv, alg::SequentialCTMRG, + cache::EmptyProjectorCache, + ) truncation_error = zero(real(scalartype(network))) for _ in 1:4 # rotate for col in 1:size(network, 2) # left move column-wise @@ -69,7 +85,62 @@ function ctmrg_iteration(network, env::CTMRGEnv, alg::SequentialCTMRG) network = rotate_north(network, EAST) env = rotate_north(env, EAST) end - return env, (; contraction_metrics = (; truncation_error)) + return env, (; contraction_metrics = (; truncation_error)), cache +end + +"""Replace the updated direction using projector columns without mutating the old cache.""" +function _replace_subspace_phase(grid, columns, key::Symbol) + isnothing(grid) && return nothing + return map(CartesianIndices(grid)) do I + return I[1] == WEST ? getproperty(columns[I[3]], key)[I[2]] : grid[I] + end +end + +function ctmrg_iteration( + network, env::CTMRGEnv, + alg::SequentialCTMRG{<:SubspaceIterationProjector}, + cache::SubspaceIterationCache, + ) + truncation_error = zero(real(scalartype(network))) + + if !_compatible_subspace_cache(cache, env, alg.projector_alg) + cache = initialize_projector_cache(network, env, alg) + end + grids = (cache.X, cache.Y) + + # Perform four directional moves, expressing each one as a column left-move. + for _ in 1:4 + # Select the rangefinders aligned with the current left-move direction. + Xdirection = grids[1][WEST, :, :] + Ydirection = _subspace_grid_entry(grids[2], WEST, :, :) + + columns = () + for col in 1:size(network, 2) + env, info = ctmrg_leftmove( + col, network, env, alg, Xdirection, Ydirection, + _subspace_charge_refresh(cache, alg.projector_alg) + ) + columns = (columns..., info) + ignore_derivatives() do + truncation_error = max(truncation_error, info.truncation_error) + end + end + + grids = map( + (grid, key) -> _rotate_subspace_grid(_replace_subspace_phase(grid, columns, key)), + grids, (:X, :Y), + ) + + # Rotate the network, environment, and cache into the frame of the next left move. + network = rotate_north(network, EAST) + env = rotate_north(env, EAST) + end + + cache′, subspace_error = _update_subspace_cache( + cache, grids... + ) + info = (; contraction_metrics = (; truncation_error, subspace_error)) + return env, info, cache′ end """ @@ -79,45 +150,32 @@ end Compute CTMRG projectors in the `:SequentialCTMRG` scheme either for an entire column `col` or for a specific `coordinate` (where `dir=WEST` is already implied in the `:SequentialCTMRG` scheme). """ -function sequential_projectors(col::Int, network, env::CTMRGEnv, alg::ProjectorAlgorithm) +function sequential_projectors(col::Int, network, env::CTMRGEnv, alg::ProjectorAlgorithm, args...) coordinates = eachcoordinate(env)[:, col] T_dst = Base.promote_op( - sequential_projectors, NTuple{3, Int}, typeof(network), typeof(env), typeof(alg) + sequential_projectors, NTuple{3, Int}, typeof(network), typeof(env), typeof(alg), + map(typeof, _projector_arguments(first(coordinates), args...))..., ) proj_and_info = similar(coordinates, T_dst) proj_and_info′::typeof(proj_and_info) = dtmap!!(proj_and_info, coordinates) do (r, c) - proj, info = sequential_projectors((WEST, r, c), network, env, alg) - return proj, info + return sequential_projectors((WEST, r, c), network, env, alg, _projector_arguments((r, c), args...)...) end return _split_proj_and_info(proj_and_info′) end function sequential_projectors( - coordinate::NTuple{3, Int}, network, env::CTMRGEnv, alg::HalfInfiniteProjector + coordinate::NTuple{3, Int}, network, env::CTMRGEnv, alg::Union{HalfInfiniteProjector, FullInfiniteProjector} ) - _, r, c = coordinate - r′ = r - 1 - trunc = truncation_strategy(alg, edge(env, WEST, r′, c)) + ec, trunc = enlarged_corner_inputs(network, coordinate, env, alg, Val(alg isa HalfInfiniteProjector ? 2 : 4)) alg´ = _set_decomposition_truncation(alg, trunc) - Q1 = TensorMap(EnlargedCorner(network, env, (SOUTHWEST, r, c))) - Q2 = TensorMap(EnlargedCorner(network, env, (NORTHWEST, r′, c))) - return compute_projector((Q1, Q2), alg´) + return compute_projector(ec, alg´) end + function sequential_projectors( - coordinate::NTuple{3, Int}, network, env::CTMRGEnv, alg::FullInfiniteProjector - ) - rowsize, colsize = size(env)[2:3] - coordinate_nw = _next_coordinate(coordinate, rowsize, colsize) - coordinate_ne = _next_coordinate(coordinate_nw, rowsize, colsize) - coordinate_se = _next_coordinate(coordinate_ne, rowsize, colsize) - trunc = truncation_strategy(alg, env.edges[WEST, coordinate_nw[2:3]...]) - alg´ = _set_decomposition_truncation(alg, trunc) - ec = ( - TensorMap(EnlargedCorner(network, env, coordinate_se)), - TensorMap(EnlargedCorner(network, env, coordinate)), - TensorMap(EnlargedCorner(network, env, coordinate_nw)), - TensorMap(EnlargedCorner(network, env, coordinate_ne)), + coordinate::NTuple{3, Int}, network, env::CTMRGEnv, + alg::SubspaceIterationProjector, X, Y, refresh::Bool = false, ) - return compute_projector(ec, alg´) + ec, trunc = enlarged_corner_inputs(network, coordinate, env, alg, Val(4)) + return compute_projector(ec, alg, (; X, Y, trunc, refresh)) end """ diff --git a/src/algorithms/ctmrg/simultaneous.jl b/src/algorithms/ctmrg/simultaneous.jl index ecaf52c7b..8037c0ecd 100644 --- a/src/algorithms/ctmrg/simultaneous.jl +++ b/src/algorithms/ctmrg/simultaneous.jl @@ -21,7 +21,7 @@ For a full description, see [`leading_boundary`](@ref). The supported keywords a * `verbosity::Int=$(Defaults.ctmrg_verbosity)` * `trunc::Union{TruncationStrategy,NamedTuple}=(; alg::Symbol=:$(Defaults.trunc))` * `decomposition_alg::Union{<:SVDAdjoint,NamedTuple}` -* `projector_alg::Symbol=:$(Defaults.projector_alg)` +* `projector_alg::Union{Symbol,NamedTuple}=:$(Defaults.projector_alg)` """ struct SimultaneousCTMRG{P <: ProjectorAlgorithm} <: CTMRGAlgorithm tol::Float64 @@ -36,7 +36,8 @@ end CTMRG_SYMBOLS[:SimultaneousCTMRG] = SimultaneousCTMRG -function ctmrg_iteration(network, env::CTMRGEnv, alg::SimultaneousCTMRG) +"""Expand every corner tensor for a simultaneous CTMRG iteration.""" +function _expand_corners_simultaneously(network, env::CTMRGEnv) coordinates = eachcoordinate(network, 1:4) T_corners = Base.promote_op( TensorMap ∘ EnlargedCorner, typeof(network), typeof(env), eltype(coordinates) @@ -45,14 +46,49 @@ function ctmrg_iteration(network, env::CTMRGEnv, alg::SimultaneousCTMRG) enlarged_corners::typeof(enlarged_corners′) = dtmap!!(enlarged_corners′, eachcoordinate(network, 1:4)) do idx return TensorMap(EnlargedCorner(network, env, idx)) - end # expand environment + end + return enlarged_corners +end + +function ctmrg_iteration( + network, env::CTMRGEnv, alg::SimultaneousCTMRG, + cache::EmptyProjectorCache, + ) + enlarged_corners = _expand_corners_simultaneously(network, env) projectors, info = simultaneous_projectors(enlarged_corners, env, alg.projector_alg) # compute projectors on all coordinates env′ = renormalize_simultaneously(enlarged_corners, projectors, network, env) # renormalize enlarged corners info = (; contraction_metrics = (; info.truncation_error), info.U, info.S, info.V, ) - return env′, info + return env′, info, cache +end + +function ctmrg_iteration( + network, env::CTMRGEnv, + alg::SimultaneousCTMRG{<:SubspaceIterationProjector}, + cache::SubspaceIterationCache, + ) + if !_compatible_subspace_cache(cache, env, alg.projector_alg) + cache = initialize_projector_cache(network, env, alg) + end + enlarged_corners = _expand_corners_simultaneously(network, env) + projectors, projector_info = simultaneous_projectors( + enlarged_corners, env, alg.projector_alg, cache.X, cache.Y, + _subspace_charge_refresh(cache, alg.projector_alg) + ) + env′ = renormalize_simultaneously(enlarged_corners, projectors, network, env) + cache′, subspace_error = _update_subspace_cache( + cache, + projector_info.X, projector_info.Y + ) + info = (; + contraction_metrics = (; + projector_info.truncation_error, subspace_error, + ), + projector_info.U, projector_info.S, projector_info.V, + ) + return env′, info, cache′ end # Work-around to stop Zygote from choking on first execution (sometimes) @@ -65,9 +101,18 @@ function _split_proj_and_info(proj_and_info) S = map(x -> x[2].S, proj_and_info) V = map(x -> x[2].V, proj_and_info) info = (; truncation_error, U, S, V) + if hasproperty(first(proj_and_info)[2], :X) + X = map(x -> x[2].X, proj_and_info) + Y = isnothing(first(proj_and_info)[2].Y) ? nothing : map(x -> x[2].Y, proj_and_info) + info = (; info..., X, Y) + end return (P_left, P_right), info end +"""Select optional recycled bases for one projector coordinate.""" +_projector_arguments(coordinate, X, Y, refresh::Bool) = (X[coordinate...], _subspace_grid_entry(Y, coordinate...), refresh) +_projector_arguments(coordinate) = () + """ simultaneous_projectors(enlarged_corners::Array{E,3}, env::CTMRGEnv, alg::ProjectorAlgorithm) simultaneous_projectors(coordinate, enlarged_corners::Array{E,3}, alg::ProjectorAlgorithm) @@ -76,45 +121,36 @@ Compute CTMRG projectors in the `:SimultaneousCTMRG` scheme either for all provi enlarged corners or on a specific `coordinate`. """ function simultaneous_projectors( - enlarged_corners::Array{E, 3}, env::CTMRGEnv, alg::ProjectorAlgorithm + enlarged_corners::Array{E, 3}, env::CTMRGEnv, alg::ProjectorAlgorithm, args... ) where {E} coordinates = eachcoordinate(env, 1:4) T_dst = Base.promote_op( simultaneous_projectors, NTuple{3, Int}, typeof(enlarged_corners), typeof(env), typeof(alg), + map(typeof, _projector_arguments(first(coordinates), args...))..., ) proj_and_info′ = similar(coordinates, T_dst) proj_and_info::typeof(proj_and_info′) = dtmap!!(proj_and_info′, coordinates) do coordinate - return simultaneous_projectors(coordinate, enlarged_corners, env, alg) + return simultaneous_projectors(coordinate, enlarged_corners, env, alg, _projector_arguments(coordinate, args...)...) end return _split_proj_and_info(proj_and_info) end function simultaneous_projectors( - coordinate, enlarged_corners::Array{E, 3}, env, alg::HalfInfiniteProjector + coordinate::NTuple{3, Int}, enlarged_corners::Array{E, 3}, env, + alg::Union{HalfInfiniteProjector, FullInfiniteProjector}, ) where {E} - coordinate′ = _next_coordinate(coordinate, size(env)[2:3]...) - trunc = truncation_strategy(alg, env.edges[coordinate[1], coordinate′[2:3]...]) + ec, trunc = enlarged_corner_inputs(enlarged_corners, coordinate, env, alg, Val(alg isa HalfInfiniteProjector ? 2 : 4)) alg′ = _set_decomposition_truncation(alg, trunc) - ec = (enlarged_corners[coordinate...], enlarged_corners[coordinate′...]) return compute_projector(ec, alg′) end + function simultaneous_projectors( - coordinate, enlarged_corners::Array{E, 3}, env, alg::FullInfiniteProjector + coordinate::NTuple{3, Int}, enlarged_corners::Array{E, 3}, env, + alg::SubspaceIterationProjector, X, Y, refresh::Bool = false, ) where {E} - rowsize, colsize = size(env)[2:3] - coordinate2 = _next_coordinate(coordinate, rowsize, colsize) - coordinate3 = _next_coordinate(coordinate2, rowsize, colsize) - coordinate4 = _next_coordinate(coordinate3, rowsize, colsize) - trunc = truncation_strategy(alg, env.edges[coordinate[1], coordinate2[2:3]...]) - alg′ = _set_decomposition_truncation(alg, trunc) - ec = ( - enlarged_corners[coordinate4...], - enlarged_corners[coordinate...], - enlarged_corners[coordinate2...], - enlarged_corners[coordinate3...], - ) - return compute_projector(ec, alg′) + ec, trunc = enlarged_corner_inputs(enlarged_corners, coordinate, env, alg, Val(4)) + return compute_projector(ec, alg, (; X, Y, trunc, refresh)) end """ diff --git a/src/algorithms/ctmrg/sparse_environments.jl b/src/algorithms/ctmrg/sparse_environments.jl index 59f7efea1..27489ecf9 100644 --- a/src/algorithms/ctmrg/sparse_environments.jl +++ b/src/algorithms/ctmrg/sparse_environments.jl @@ -61,6 +61,18 @@ function EnlargedCorner(network::InfiniteSquareNetwork, env, coordinates) end end +"""Return the domain of an enlarged corner without materializing its contraction.""" +function TensorKit.domain(Q::EnlargedCorner) + A = _rotate_north_localsandwich(Q.A, Q.dir) + return domain(Q.E_2) * _elementwise_dual(east_virtualspace(A)) +end + +"""Return the codomain of an enlarged corner without materializing its contraction.""" +function TensorKit.codomain(Q::EnlargedCorner) + A = _rotate_north_localsandwich(Q.A, Q.dir) + return first(codomain(Q.E_1)) * south_virtualspace(A) +end + """ TensorMap(Q::EnlargedCorner) diff --git a/src/algorithms/optimization/implicit_differentiation.jl b/src/algorithms/optimization/implicit_differentiation.jl index 37ae322b4..6aa90fad9 100644 --- a/src/algorithms/optimization/implicit_differentiation.jl +++ b/src/algorithms/optimization/implicit_differentiation.jl @@ -205,10 +205,51 @@ Evaluating the gradient of the cost function for CTMRG: _scrambling_env_gauge(::CTMRGAlgorithm) = ScramblingEnvGauge() function _set_fixed_truncation(alg::CTMRGAlgorithm) - alg_fixed = @set alg.projector_alg = _set_truncation(alg.projector_alg, FixedSpaceTruncation()) + projector = _set_truncation(alg.projector_alg, FixedSpaceTruncation()) + if projector isa SubspaceIterationProjector && _adaptive_subspace(alg.projector_alg) + # Hold the discrete charge allocation fixed during the joint response solve. + projector = @set projector.charge_refresh = 0 + end + alg_fixed = @set alg.projector_alg = projector return alg_fixed end +"""Convert structured AD cotangents to the vectors used by the fixed-point solver.""" +function _fixedpoint_tangent(x, dx) + dx = unthunk(dx) + dx isa AbstractZero && return VI.zerovector(x) + return x isa Tuple ? ntuple(i -> _fixedpoint_tangent(x[i], dx[i]), length(x)) : dx +end + +"""Prepare the differentiable fixed-point state and gauge-fixed iteration.""" +function _fixedpoint_iteration(env, alg, cache::EmptyProjectorCache, signs) + iteration(A, x) = fix_global_phases( + fix_relative_phases(ctmrg_iteration(InfiniteSquareNetwork(A), x, alg, cache)[1], signs), env + ) + return env, iteration, identity +end + +function _fixedpoint_iteration(env, alg, cache::SubspaceIterationCache, signs) + # Krylov vectors and AD cotangents must use the same concrete tensor storage. + bases = isnothing(cache.Y) ? (cache.X,) : (cache.X, map(copy, cache.Y)) + function iteration(A, x) + incoming = SubspaceIterationCache(cache.iteration, x[2][1], isnothing(cache.Y) ? nothing : x[2][2]) + en, _, next = ctmrg_iteration(InfiniteSquareNetwork(A), x[1], alg, incoming) + en = fix_global_phases(fix_relative_phases(en, signs), env) + nextbases = isnothing(next.Y) ? (next.X,) : (next.X, next.Y) + # Fix column-frame freedom and remove normal-to-isometry perturbations. + aligned = map(nextbases, bases) do grid, refs + map(grid, refs) do b, ref + out = b * (b' * ref) + ref + out - ref * (ref' * out) + end + end + return en, aligned + end + rhs(Δenv) = (Δenv, VI.zerovector(bases)) + return (env, bases), iteration, rhs +end + # compute the CTMRG gradient through fixed-point differentiation function _rrule( gradmode::FixedPointGradient, @@ -220,31 +261,26 @@ function _rrule( ) _check_algorithm_combination(alg, gradmode) - env, = leading_boundary(envinit, state, alg) + env, info, projector_cache = _leading_boundary_with_cache(envinit, InfiniteSquareNetwork(state), alg) # prepare iterating function corresponding to a single gauge-fixed CTMRG iteration alg_fixed = _set_fixed_truncation(alg) # fix spaces during differentiation alg_gauge = _scrambling_env_gauge(alg) # select appropriate gauge-fixing algorithm - env_conv, info = ctmrg_iteration(InfiniteSquareNetwork(state), env, alg_fixed) - signs, = compute_gauge_fix_gauge(env_conv, env, alg_gauge) - # fix global phases relative to the fixed point inside the iteration - function gauge_fixed_iteration(A, x) - return fix_global_phases( - fix_relative_phases(ctmrg_iteration(InfiniteSquareNetwork(A), x, alg_fixed)[1], signs), - env, - ) - end + env_conv, _, = ctmrg_iteration(InfiniteSquareNetwork(state), env, alg_fixed, projector_cache) + signs = first(compute_gauge_fix_gauge(env_conv, env, alg_gauge)) + x, gauge_fixed_iteration, rhs = _fixedpoint_iteration(env, alg_fixed, projector_cache, signs) # prepare its pullback - _, env_vjp = rrule_via_ad(config, gauge_fixed_iteration, state, env) + _, env_vjp = rrule_via_ad(config, gauge_fixed_iteration, state, x) # split off state and environment parts ∂f∂A(x)::typeof(state) = env_vjp(x)[2] - ∂f∂x(x)::typeof(env) = env_vjp(x)[3] + ∂f∂x(y) = _fixedpoint_tangent(x, env_vjp(y)[3]) function leading_boundary_fixed_pullback((Δenv′, Δinfo)) Δenv = unthunk(Δenv′) # evaluate the geometric sum - ∂F∂env = fixedpoint_gradient(Δenv, ∂f∂x, ∂f∂A, Δenv, gradmode.solver_alg) + Δx = rhs(Δenv) + ∂F∂env = fixedpoint_gradient(Δx, ∂f∂x, ∂f∂A, Δx, gradmode.solver_alg) return NoTangent(), ZeroTangent(), ∂F∂env, NoTangent() end @@ -310,7 +346,7 @@ function fixedpoint_gradient(∂E∂x, ∂f∂x, ∂f∂A, y₀, alg::ManualIter dx = ∂f∂A(y) ϵ = 1.0 for i in 1:(alg.maxiter) - y′ = ∂E∂x + ∂f∂x(y) + y′ = VI.add(∂E∂x, ∂f∂x(y)) dxnew = ∂f∂A(y′) ϵnew = norm(dxnew - dx) @@ -345,7 +381,8 @@ function fixedpoint_gradient(∂E∂x, ∂f∂x, ∂f∂A, x₀, alg::KrylovKit. return (VI.add!!(y, ∂E∂x, X[2]), X[2]) end X₀ = (x₀, one(scalartype(x₀))) - _, vecs, info = realeigsolve(f, X₀, 1, :LM, alg) + which = EigSorter(λ -> abs(λ - one(λ))) + _, vecs, info = realeigsolve(f, X₀, 1, which, alg) if alg.verbosity > 0 && info.converged < 1 @warn("gradient fixed-point iteration reached maximal number of iterations:", info) end @@ -426,6 +463,21 @@ function _check_algorithm_combination(::SimultaneousCTMRG{<:FullInfiniteProjecto throw(ArgumentError(msg)) end +"""Reject characteristic-equation rules that are not defined for SI projectors.""" +function _reject_subspace_iteration_gradient() + throw(ArgumentError("SubspaceIterationProjector supports FixedPointGradient or `gradient_alg=nothing`; its ImplicitGradient characteristic equations are not implemented.")) +end +function _check_algorithm_combination( + ::SequentialCTMRG{<:SubspaceIterationProjector}, ::ImplicitGradient + ) + return _reject_subspace_iteration_gradient() +end +function _check_algorithm_combination( + ::SimultaneousCTMRG{<:SubspaceIterationProjector}, ::ImplicitGradient + ) + return _reject_subspace_iteration_gradient() +end + @doc raw""" implicit_gradient(x̆, ∂ₓF, ∂ₚF, y₀, alg) diff --git a/test/ctmrg/subspace_dimensions.jl b/test/ctmrg/subspace_dimensions.jl new file mode 100644 index 000000000..c387c367c --- /dev/null +++ b/test/ctmrg/subspace_dimensions.jl @@ -0,0 +1,10 @@ +using PEPSKit + +@isdefined(TestSuite) || include("../testsuite/TestSuite.jl") +using .TestSuite + +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +if !is_buildkite + TestSuite.ctmrg_subspace_dimensions(Vector) +end diff --git a/test/ctmrg/subspace_iteration.jl b/test/ctmrg/subspace_iteration.jl new file mode 100644 index 000000000..85221a29c --- /dev/null +++ b/test/ctmrg/subspace_iteration.jl @@ -0,0 +1,10 @@ +using PEPSKit + +@isdefined(TestSuite) || include("../testsuite/TestSuite.jl") +using .TestSuite + +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +if !is_buildkite + TestSuite.ctmrg_subspace_iteration(Vector) +end diff --git a/test/ctmrg/subspace_updates.jl b/test/ctmrg/subspace_updates.jl new file mode 100644 index 000000000..40b9974f8 --- /dev/null +++ b/test/ctmrg/subspace_updates.jl @@ -0,0 +1,10 @@ +using PEPSKit + +@isdefined(TestSuite) || include("../testsuite/TestSuite.jl") +using .TestSuite + +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +if !is_buildkite + TestSuite.ctmrg_subspace_updates(Vector) +end diff --git a/test/gradients/subspace_dimensions.jl b/test/gradients/subspace_dimensions.jl new file mode 100644 index 000000000..71f748ef7 --- /dev/null +++ b/test/gradients/subspace_dimensions.jl @@ -0,0 +1,10 @@ +using PEPSKit + +@isdefined(TestSuite) || include("../testsuite/TestSuite.jl") +using .TestSuite + +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +if !is_buildkite + TestSuite.gradients_subspace_dimensions(Vector) +end diff --git a/test/gradients/subspace_iteration.jl b/test/gradients/subspace_iteration.jl new file mode 100644 index 000000000..837311dc2 --- /dev/null +++ b/test/gradients/subspace_iteration.jl @@ -0,0 +1,10 @@ +using PEPSKit + +@isdefined(TestSuite) || include("../testsuite/TestSuite.jl") +using .TestSuite + +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +if !is_buildkite + TestSuite.gradients_subspace_iteration(Vector) +end diff --git a/test/gradients/subspace_updates.jl b/test/gradients/subspace_updates.jl new file mode 100644 index 000000000..9e4e8366c --- /dev/null +++ b/test/gradients/subspace_updates.jl @@ -0,0 +1,10 @@ +using PEPSKit + +@isdefined(TestSuite) || include("../testsuite/TestSuite.jl") +using .TestSuite + +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +if !is_buildkite + TestSuite.gradients_subspace_updates(Vector) +end diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl index 3e576386d..157723e68 100644 --- a/test/testsuite/TestSuite.jl +++ b/test/testsuite/TestSuite.jl @@ -164,6 +164,24 @@ module CTMRGUnitCell end using .CTMRGUnitCell +module CTMRGSubspaceDimensions + include("ctmrg/subspace_dimensions.jl") + export ctmrg_subspace_dimensions +end +using .CTMRGSubspaceDimensions + +module CTMRGSubspaceIteration + include("ctmrg/subspace_iteration.jl") + export ctmrg_subspace_iteration +end +using .CTMRGSubspaceIteration + +module CTMRGSubspaceUpdates + include("ctmrg/subspace_updates.jl") + export ctmrg_subspace_updates +end +using .CTMRGSubspaceUpdates + # Examples # -------- module ExamplesBoseHubbard @@ -210,6 +228,24 @@ module GradientsCTMRG end using .GradientsCTMRG +module GradientsSubspaceDimensions + include("gradients/subspace_dimensions.jl") + export gradients_subspace_dimensions +end +using .GradientsSubspaceDimensions + +module GradientsSubspaceIteration + include("gradients/subspace_iteration.jl") + export gradients_subspace_iteration +end +using .GradientsSubspaceIteration + +module GradientsSubspaceUpdates + include("gradients/subspace_updates.jl") + export gradients_subspace_updates +end +using .GradientsSubspaceUpdates + # Optimization # ------------ module OptimizationNormPreservation diff --git a/test/testsuite/ctmrg/subspace_dimensions.jl b/test/testsuite/ctmrg/subspace_dimensions.jl new file mode 100644 index 000000000..d460e6eac --- /dev/null +++ b/test/testsuite/ctmrg/subspace_dimensions.jl @@ -0,0 +1,141 @@ +using Test +using Random +using LinearAlgebra +using MatrixAlgebraKit +using MatrixAlgebraKit: TruncatedAlgorithm +using TensorKit +using PEPSKit +using Adapt + +"""Check integer budgets, new-charge discovery, and fixed target multiplicities.""" +function ctmrg_subspace_dimensions_selection(AT) + return @testset "Integer subspace allocation" begin + W = U1Space(0 => 3, 1 => 3, 2 => 3) + A, changed = (zeros(ComplexF64, W ← W) for _ in 1:2) + for (hot, M) in ((2, A), (1, changed)), (c, b) in blocks(M) + b[diagind(b)] .= c == U1Irrep(hot) ? [9, 8, 7] : c == U1Irrep(0) ? [6, 5, 4] : [3, 2, 1] + end + A, changed, trunc = adapt(AT, A), adapt(AT, changed), truncrank(4) + @testset "$update" for update in (:independent, :alternating, :alternating_refined) + Random.seed!(273) + alg = SubspaceIterationProjector(; subspace = 6, trunc, update) + first_step = PEPSKit._subspace_iteration_decomposition(A, id(W), alg, (; X = nothing, Y = nothing, trunc)) + @test fuse(domain(first_step.X)) == U1Space(0 => 3, 2 => 3) + incoming = (; first_step.X, first_step.Y, trunc) + result = PEPSKit._subspace_iteration_decomposition(changed, id(W), alg, (; incoming..., refresh = true)) + @test fuse(domain(result.X)) == U1Space(0 => 3, 1 => 3) + @test result.X' * result.X ≈ id(ComplexF64, domain(result.X)) + u, s, v, = svd_trunc(changed, TruncatedAlgorithm(DefaultAlgorithm(), trunc)) + actual, expected = result.U * result.S * result.V, u * s * v + @test actual / norm(actual) ≈ expected / norm(expected) atol = 1.0e-10 + fixed = truncspace(U1Space(0 => 2, 1 => 2)) + result = PEPSKit._subspace_iteration_decomposition(A, id(W), alg, (; X = nothing, Y = nothing, trunc = fixed)) + @test fuse(domain(result.X)) == U1Space(0 => 2, 1 => 2, 2 => 2) + end + @testset "$S" for (S, charges) in ((SU2Irrep, (0, 1 // 2)), (A4Irrep, (0, 3)), (FibonacciAnyon, (:I, :τ))) + W = Vect[S](charges[1] => 3, charges[2] => 2) + M = zeros(ComplexF64, W ← W) + for (c, b) in blocks(M) + b[diagind(b)] .= (c == S(charges[1]) ? 3 : 10) .- (0:(size(b, 1) - 1)) + end + M = adapt(AT, M) + alg = SubspaceIterationProjector(; subspace = 6, trunc = truncrank(3)) + result = PEPSKit._subspace_iteration_decomposition(M, id(W), alg, (; X = nothing, Y = nothing, trunc = alg.trunc)) + _, s, _, = svd_trunc(M, TruncatedAlgorithm(DefaultAlgorithm(), alg.trunc)) + @test space(result.S, 1) == space(s, 1) + @test dim(domain(result.X)) <= 6 + end + end +end + +"""Compare rank-changing sweeps with full SVD and check refresh and fixed-point policies.""" +function ctmrg_subspace_dimensions_sweeps(AT) + return @testset "Adaptive CTMRG sweeps" begin + @testset "$flavor, $update, probe=$probe" for flavor in (SimultaneousCTMRG, SequentialCTMRG), update in (:independent, :alternating, :alternating_refined), probe in (nothing, 1) + Random.seed!(274) + P, V = U1Space(-1 => 1, 0 => 1, 1 => 1), U1Space(0 => 1, 1 => 1) + ψ = adapt(AT, InfinitePEPS(P, V; unitcell = (2, 1))) + network = InfiniteSquareNetwork(ψ) + env = CTMRGEnv(ψ, U1Space(0 => 2)) + projector = SubspaceIterationProjector(; subspace = 4, trunc = truncrank(2), update, charge_refresh = 1, charge_probe = probe) + alg = flavor(0.0, 2, 2, 0, projector) + reference = flavor(; projector_alg = :FullInfiniteProjector, trunc = truncrank(2), verbosity = 0) + cache = PEPSKit.initialize_projector_cache(network, env, alg) + env_ref = env + for _ in 1:2 + env, _, cache = PEPSKit.ctmrg_iteration(network, env, alg, cache) + env_ref, = PEPSKit.ctmrg_iteration(network, env_ref, reference) + end + @test norm(ψ, env) ≈ norm(ψ, env_ref) rtol = 1.0e-7 + @test all(X -> dim(domain(X)) <= 4, cache.X) + fixed = PEPSKit._set_fixed_truncation(alg).projector_alg + @test fixed.trunc isa FixedSpaceTruncation && !PEPSKit._subspace_charge_refresh(cache, fixed) + end + periodic = SubspaceIterationProjector(; subspace = 24, charge_refresh = 10) + cache = PEPSKit.SubspaceIterationCache(3, nothing, nothing) + @test !PEPSKit._projector_convergence_ready(cache, periodic) + @test PEPSKit._projector_convergence_ready(PEPSKit.@set(cache.iteration = 11), periodic) + @test PEPSKit._subspace_charge_refresh(PEPSKit.@set(cache.iteration = 10), periodic) + end +end + +"""Discover absent charges incrementally and recover broad probing when requested.""" +function ctmrg_subspace_dimensions_hybrid(AT) + return @testset "Hybrid charge probes" begin + @test isnothing(SubspaceIterationProjector(; subspace = 6).charge_probe) + for probe in (0, -1, true, 1.5) + @test_throws ArgumentError SubspaceIterationProjector(; subspace = 6, charge_probe = probe) + end + W = U1Space(0 => 3, 1 => 3, 2 => 3) + A = zeros(ComplexF64, W ← W) + for (c, b) in blocks(A) + b[diagind(b)] .= c == U1Irrep(1) ? [100, 90, 80] : c == U1Irrep(0) ? [6, 5, 4] : [3, 2, 1] + end + A = adapt(AT, A) + initial = U1Space(0 => 3, 2 => 3) + @testset "$update" for update in (:independent, :alternating, :alternating_refined) + alg = SubspaceIterationProjector(; subspace = 6, trunc = truncrank(4), charge_probe = 1, update) + Random.seed!(277) + X, Y = PEPSKit.rangefinder_initialize(A, id(W), initial, alg) + incoming = (; X, Y, trunc = alg.trunc, refresh = true) + K, refresh = PEPSKit._adaptive_subspace_space(A, id(W), alg, incoming) + @test refresh && K == U1Space(0 => 3, 1 => 1, 2 => 3) + broad = PEPSKit.@set alg.charge_probe = nothing + large = PEPSKit.@set alg.charge_probe = 6 + @test first(PEPSKit._adaptive_subspace_space(A, id(W), broad, incoming)) == W + @test first(PEPSKit._adaptive_subspace_space(A, id(W), large, incoming)) == W + @test first(PEPSKit._adaptive_subspace_space(A, id(W), alg, (; incoming..., X = nothing))) == W + @test first(PEPSKit._adaptive_subspace_space(A, id(W), alg, (; incoming..., X = real(X)))) == W + local result + for n in 1:3 + probe, = PEPSKit._adaptive_subspace_space(A, id(W), alg, incoming) + expanded = PEPSKit._resize_subspace_basis(incoming.X, codomain(A), probe, alg, ComplexF64) + @test expanded' * expanded ≈ id(ComplexF64, domain(expanded)) atol = 1.0e-12 + result = PEPSKit._subspace_iteration_decomposition(A, id(W), alg, incoming) + selected = fuse(domain(result.X)) + @test dim(selected, U1Irrep(1)) == n + @test dim(selected) <= 6 + incoming = (; result.X, result.Y, trunc = alg.trunc, refresh = true) + end + @test fuse(domain(result.X)) == U1Space(0 => 3, 1 => 3) + u, s, v, = svd_trunc(A, TruncatedAlgorithm(DefaultAlgorithm(), alg.trunc)) + actual, expected = result.U * result.S * result.V, u * s * v + @test actual / norm(actual) ≈ expected / norm(expected) atol = 1.0e-10 + fixed = truncspace(U1Space(0 => 2, 1 => 2)) + @test first(PEPSKit._adaptive_subspace_space(A, id(W), alg, (; incoming..., X, Y, trunc = fixed))) == W + end + W = SU2Space(0 => 6, 1 // 2 => 4) + alg = SubspaceIterationProjector(; subspace = 6, charge_probe = 1) + current = SU2Space(0 => 2, 1 // 2 => 1) + @test PEPSKit._subspace_probe_space(W, alg, truncrank(3), current) == SU2Space(0 => 3, 1 // 2 => 2) + end +end + +"""Run focused checks for automatically allocated integer trial dimensions.""" +function ctmrg_subspace_dimensions(AT) + return @testset "Subspace dimensions ($AT)" begin + ctmrg_subspace_dimensions_hybrid(AT) + ctmrg_subspace_dimensions_selection(AT) + ctmrg_subspace_dimensions_sweeps(AT) + end +end diff --git a/test/testsuite/ctmrg/subspace_iteration.jl b/test/testsuite/ctmrg/subspace_iteration.jl new file mode 100644 index 000000000..86c4d8c4e --- /dev/null +++ b/test/testsuite/ctmrg/subspace_iteration.jl @@ -0,0 +1,345 @@ +using Test +using Random +using LinearAlgebra +using MatrixAlgebraKit +using TensorKit +using PEPSKit +using Adapt +using PEPSKit: ctmrg_iteration, initialize_projector_cache + +"""Construct a CTMRG flavor configured with the subspace iteration projector for testing.""" +function si_ctmrg(flavor; subspace, projector_kwargs = (;), kwargs...) + projector_alg = merge( + (; alg = :SubspaceIterationProjector, subspace), projector_kwargs + ) + return flavor(; projector_alg, verbosity = 0, kwargs...) +end + +"""Check configuration for the subspace iteration projector.""" +function ctmrg_subspace_iteration_configuration(AT) + return @testset "SubspaceIterationProjector configuration" begin + alg = SubspaceIterationProjector(; subspace = ℂ^10) + @test alg.update == :independent + @test alg.n_subspace_iter == 0 && alg.rotate_basis + + configured = SimultaneousCTMRG(; + projector_alg = (; + alg = :SubspaceIterationProjector, + subspace = ℂ^9, + orth_alg = (; fwd_alg = (; alg = :Householder)), + ), + ) + @test configured.projector_alg.orth_alg.fwd_alg isa Householder + + @test_throws ArgumentError SimultaneousCTMRG(; projector_alg = :SubspaceIterationProjector) + @test SubspaceIterationProjector(; subspace = 10).subspace == 10 + for subspace in (0, -1, true, truncrank(24)) + @test_throws ArgumentError SubspaceIterationProjector(; subspace) + end + @test_throws ArgumentError SubspaceIterationProjector(; subspace = 12, trunc = truncrank(16)) + @test_throws ArgumentError SubspaceIterationProjector(; subspace = 24, trunc = truncrank(16; rev = false)) + @test_throws ArgumentError SubspaceIterationProjector(; subspace = 24, charge_refresh = -1) + @test_throws ArgumentError SubspaceIterationProjector(; subspace = ℂ^10, update = :unknown) + @test_throws ArgumentError SubspaceIterationProjector(; subspace = ℂ^10, n_subspace_iter = -1) + @test_throws ArgumentError SubspaceIterationProjector(; + subspace = ℂ^10, trunc = truncerror(; rtol = 1.0e-8) + ) + end +end + +"""Check convergence for the subspace iteration projector.""" +function ctmrg_subspace_iteration_convergence(AT) + return @testset "Dominant subspace convergence" begin + Random.seed!(0x51ab) + V = ComplexSpace(4) + L = adapt(AT, TensorMap(Matrix(Diagonal(ComplexF64[5, 3, 0.05, 0.01])), V ← V)) + R = id(V) + L /= norm(L) + R /= norm(R) + trunc = truncspace(ComplexSpace(2)) + full_alg = MatrixAlgebraKit.TruncatedAlgorithm(DefaultAlgorithm(), trunc) + Ufull, Sfull, Vfull, = svd_trunc(L * R, full_alg) + alg = SubspaceIterationProjector(; subspace = ℂ^3) + rangefinders = (; X = nothing, Y = nothing, trunc) + local si + for _ in 1:5 + si = PEPSKit._subspace_iteration_decomposition(L, R, alg, rangefinders) + rangefinders = (; X = si.X, Y = si.Y, trunc) + end + + @test norm(si.U * si.U' - Ufull * Ufull') < 1.0e-4 + @test norm(si.V' * si.V - Vfull' * Vfull) < 1.0e-4 + @test si.S / norm(si.S) ≈ Sfull / norm(Sfull) rtol = 1.0e-6 + end +end + +"""Check half environments for the subspace iteration projector.""" +function ctmrg_subspace_iteration_half_environments(AT) + return @testset "Factorized half environments" begin + Random.seed!(0x51fa) + Q = [adapt(AT, randn(ComplexF64, ℂ^5 ← ℂ^5)) for _ in 1:4] + left, right = PEPSKit.SIHalfEnvironment(Q[1], Q[2]), PEPSKit.SIHalfEnvironment(Q[3], Q[4]) + alg = SubspaceIterationProjector(; subspace = ℂ^3) + X, Y = PEPSKit.rangefinder_initialize(left, right, alg.subspace, alg) + inputs = (; X, Y, trunc = truncspace(ℂ^2)) + implicit = PEPSKit._subspace_iteration_decomposition(left, right, alg, inputs) + L = Q[1] * Q[2] / (norm(Q[1]) * norm(Q[2])) + R = Q[3] * Q[4] / (norm(Q[3]) * norm(Q[4])) + explicit = PEPSKit._subspace_iteration_decomposition(L, R, alg, inputs) + @test implicit.S ≈ explicit.S + @test implicit.U * implicit.S * implicit.V ≈ explicit.U * explicit.S * explicit.V + P = PEPSKit.contract_projectors(implicit.U, implicit.S, implicit.V, left, right) + P′ = PEPSKit.contract_projectors(explicit.U, explicit.S, explicit.V, L, R) + @test P[1] * P[2] ≈ P′[1] * P′[2] + end +end + +"""Check rangefinders for the subspace iteration projector.""" +function ctmrg_subspace_iteration_rangefinders(AT) + return @testset "TensorMap rangefinders" begin + Random.seed!(0x51c7) + psi = adapt(AT, InfinitePEPS(ComplexSpace(2), ComplexSpace(2); unitcell = (1, 1))) + network = InfiniteSquareNetwork(psi) + env = CTMRGEnv(psi, ComplexSpace(4)) + alg = si_ctmrg(SimultaneousCTMRG; subspace = ℂ^104) + + for coordinate in PEPSKit.eachcoordinate(env, 1:4) + corner = PEPSKit.EnlargedCorner(network, env, coordinate) + materialized_corner = adapt(AT, TensorMap(corner)) + @test domain(corner) == domain(materialized_corner) + @test codomain(corner) == codomain(materialized_corner) + end + + cache = initialize_projector_cache(network, env, alg) + _, _, cache = ctmrg_iteration(network, env, alg, cache) + X = first(cache.X) + @test dim(domain(X)) < dim(alg.projector_alg.subspace) + + larger_env = CTMRGEnv(psi, ComplexSpace(5)) + _, _, reset_cache = ctmrg_iteration(network, larger_env, alg, cache) + @test space(first(reset_cache.X)) != space(first(cache.X)) + end +end + +"""Check phase caches for the subspace iteration projector.""" +function ctmrg_subspace_iteration_phase_caches(AT) + return @testset "Sequential subspace phase caches" begin + Random.seed!(0x5e90) + psi = adapt(AT, InfinitePEPS(ComplexSpace(2), ComplexSpace(1); unitcell = (2, 3))) + network = InfiniteSquareNetwork(psi) + env = CTMRGEnv(psi, ComplexSpace(2)) + alg = si_ctmrg(SequentialCTMRG; subspace = ℂ^7) + cache = initialize_projector_cache(network, env, alg) + env, _, cache = ctmrg_iteration(network, env, alg, cache) + _, _, cache = ctmrg_iteration(network, env, alg, cache) + @test all(grid -> size(grid) == (4, 2, 3), (cache.X, cache.Y)) + + basis = adapt(AT, isometry(ComplexF64, ℂ^3 ← ℂ^2)) + rotated = copy(basis) + block(rotated, Trivial())[1, 1] = cos(1.0e-8) + block(rotated, Trivial())[3, 1] = sin(1.0e-8) + @test PEPSKit._subspace_error([basis], [rotated]) ≈ sin(1.0e-8) / sqrt(2) rtol = 1.0e-8 + rotation = first(left_orth(adapt(AT, randn(ComplexF64, ℂ^2 ← ℂ^2)))) + @test PEPSKit._subspace_error([basis], [basis * rotation]) < 1.0e-14 + + values = reshape(collect(1:24), 4, 2, 3) + original = PEPSKit.SubspaceIterationCache(2, values, -values) + copied = copy(original) + @test all( + pair -> pair[1] == pair[2] && pair[1] !== pair[2], zip( + (copied.X, copied.Y), + (original.X, original.Y), + ) + ) + @test copied.iteration == 2 + rotated = original + for _ in 1:4 + rotated = rotl90(rotated) + end + @test (rotated.X, rotated.Y) == + (values, -values) + @test rotated.iteration == 2 + end +end + +"""Check unequal spaces for the subspace iteration projector.""" +function ctmrg_subspace_iteration_unequal_spaces(AT) + return @testset "Unequal dense spaces" begin + Pspaces = fill(ComplexSpace(2), 2, 2) + Nspaces = ComplexSpace.([2 3; 3 2]) + Espaces = ComplexSpace.([3 2; 2 3]) + chis_north = ComplexSpace.([3 4; 4 3]) + chis_east = ComplexSpace.([4 3; 3 4]) + chis_south = ComplexSpace.([3 4; 4 3]) + chis_west = ComplexSpace.([4 3; 3 4]) + psi = adapt(AT, InfinitePEPS(randn, ComplexF64, Pspaces, Nspaces, Espaces)) + env = CTMRGEnv( + randn, + ComplexF64, + psi, + chis_north, + chis_east, + chis_south, + chis_west, + ) + @testset "$flavor" for flavor in (SimultaneousCTMRG, SequentialCTMRG) + alg = si_ctmrg(flavor; subspace = ℂ^9) + env′, = ctmrg_iteration(InfiniteSquareNetwork(psi), env, alg) + @test all(space.(env′.corners) .== space.(env.corners)) + @test all(space.(env′.edges) .== space.(env.edges)) + end + end +end + +"""Check unitcells for the subspace iteration projector.""" +function ctmrg_subspace_iteration_unitcells(AT) + return @testset "Dense unit cells" begin + flavors = (SimultaneousCTMRG, SequentialCTMRG) + unitcells = ((1, 1), (2, 3)) + @testset "$flavor, $unitcell" for (flavor_index, flavor) in enumerate(flavors), + (unitcell_index, unitcell) in enumerate(unitcells) + Random.seed!(0x51ce + 10 * flavor_index + unitcell_index) + psi = adapt(AT, InfinitePEPS(ComplexSpace(2), ComplexSpace(1); unitcell)) + env = CTMRGEnv(psi, ComplexSpace(2)) + common = (; maxiter = 2, miniter = 2, tol = 0.0) + env_full, = leading_boundary( + env, + psi, + flavor(; projector_alg = :FullInfiniteProjector, verbosity = 0, common...), + ) + env_si, = leading_boundary( + env, psi, si_ctmrg(flavor; subspace = ℂ^7, common...) + ) + + @test abs(norm(psi, env_si)) ≈ abs(norm(psi, env_full)) rtol = 1.0e-7 + corner_distance = maximum( + splat(PEPSKit._singular_value_distance), + zip(map(svd_vals, env_si.corners), map(svd_vals, env_full.corners)), + ) + edge_distance = maximum( + splat(PEPSKit._singular_value_distance), + zip(map(svd_vals, env_si.edges), map(svd_vals, env_full.edges)), + ) + @test corner_distance < 1.0e-7 + @test edge_distance < 1.0e-7 + + H = adapt(AT, heisenberg_XYZ(InfiniteSquare(unitcell...))) + @test cost_function(psi, env_si, H) ≈ cost_function(psi, env_full, H) rtol = 1.0e-7 + end + end +end + +"""Check symmetry sectors for the subspace iteration projector.""" +function ctmrg_subspace_iteration_symmetry_sectors(AT) + return @testset "Symmetry sectors" for (S, charges) in ( + (U1Irrep, (-1 // 2, 1 // 2, 3 // 2)), + (Z2Irrep, (0, 1)), + (Z3Irrep, (0, 1, 2)), + (SU2Irrep, (0, 1 // 2, 1)), + (A4Irrep, (0, 1, 3)), + (FermionParity, (0, 1)), + (FermionParity ⊠ U1Irrep, ((0, 0), (1, 1), (1, -1))), + (FibonacciAnyon, (:I, :τ)), + ) + Random.seed!(0x51c4) + W = Vect[S](c => d for (c, d) in zip(charges, (5, 6, 4))) + K = Vect[S](c => d for (c, d) in zip(charges, (3, 4, 3))) + χ = Vect[S](c => d for (c, d) in zip(charges, (2, 2, 1))) + L = adapt(AT, zeros(ComplexF64, W ← W)) + for (c, b) in blocks(L) + r = dim(χ, c) + b[diagind(b)] .= [i <= r ? 10.0 - i : 0.1 / i for i in axes(b, 1)] + end + L = first(left_orth(adapt(AT, randn(ComplexF64, space(L))))) * L + R = first(left_orth(adapt(AT, randn(ComplexF64, space(L))))) + alg = SubspaceIterationProjector(; subspace = K) + X, Y = PEPSKit.rangefinder_initialize(L, R, K, alg) + inputs = (; X, Y, trunc = truncspace(χ)) + local result + for _ in 1:3 + result = PEPSKit._subspace_iteration_decomposition(L, R, alg, inputs) + inputs = (; result.X, result.Y, inputs.trunc) + end + U, S, V, = svd_trunc(L * R, MatrixAlgebraKit.TruncatedAlgorithm(DefaultAlgorithm(), truncspace(χ))) + @test space(result.S, 1) == χ + @test domain(result.X) == domain(result.Y) == ProductSpace(K) + @test result.U * result.S * result.V ≈ U * S * V rtol = 1.0e-7 + @test result.U' * result.U ≈ id(ComplexF64, χ) + @test result.V * result.V' ≈ id(ComplexF64, χ) + + end +end + +"""Check symmetric unitcells for the subspace iteration projector.""" +function ctmrg_subspace_iteration_symmetric_unitcells(AT) + return @testset "Symmetric unit cells" begin + cases = ( + (ComplexF64, U1Space(-1 => 1, 0 => 1, 1 => 1), U1Space(0 => 1, 1 => 1), U1Space(-1 => 1, 0 => 2, 1 => 2)), + (Float64, Z2Space(0 => 1, 1 => 1), Z2Space(0 => 1, 1 => 1), Z2Space(0 => 2, 1 => 1)), + (Float64, Z3Space(0 => 1, 1 => 1, 2 => 1), Z3Space(0 => 1, 1 => 1), Z3Space(0 => 2, 1 => 1, 2 => 1)), + (ComplexF64, SU2Space(0 => 1, 1 => 1), SU2Space(1 // 2 => 1), SU2Space(0 => 1, 1 => 1)), + (ComplexF64, Vect[A4Irrep](3 => 1), Vect[A4Irrep](3 => 1), Vect[A4Irrep](0 => 1, 3 => 1)), + (Float64, Vect[FermionParity](0 => 1, 1 => 1), Vect[FermionParity](0 => 1, 1 => 1), Vect[FermionParity](0 => 2, 1 => 1)), + ( + ComplexF64, Vect[FermionParity ⊠ U1Irrep]((0, 0) => 1, (1, 1) => 1, (1, -1) => 1), + Vect[FermionParity ⊠ U1Irrep]((0, 0) => 1, (1, 1) => 1, (1, -1) => 1), + Vect[FermionParity ⊠ U1Irrep]((0, 0) => 1, (1, 1) => 1, (1, -1) => 1), + ), + ) + @testset "$T, $P, $V, $χ, $flavor, $unitcell" for (T, P, V, χ) in cases, flavor in (SimultaneousCTMRG, SequentialCTMRG), unitcell in ((1, 1), (2, 1)) + Random.seed!(0x51c5) + K = fuse(χ ⊗ V ⊗ V') + psi = adapt(AT, InfinitePEPS(randn, T, fill(P, unitcell), fill(V, unitcell), fill(V, unitcell))) + env = CTMRGEnv(randn, T, psi, χ) + alg = si_ctmrg(flavor; subspace = K, maxiter = 2, miniter = 2, tol = 0.0) + env′, info = leading_boundary(env, psi, alg) + @test all(space.(env′.corners) .== space.(env.corners)) + @test all(space.(env′.edges) .== space.(env.edges)) + @test isfinite(norm(env′)) && isfinite(info.convergence_error) + end + end +end + +"""Check compatibility for the subspace iteration projector.""" +function ctmrg_subspace_iteration_compatibility(AT) + return @testset "Space compatibility and gradient configuration" begin + P = Z2Space(0 => 1) + V = Z2Space(0 => 1, 1 => 1) + psi = adapt(AT, InfinitePEPS(randn, ComplexF64, fill(P, 1, 1), fill(V, 1, 1), fill(V, 1, 1))) + env = CTMRGEnv(psi, V) + alg = SimultaneousCTMRG(; + projector_alg = (; alg = :SubspaceIterationProjector, subspace = V) + ) + @test_throws ArgumentError leading_boundary( + env, psi, si_ctmrg(SimultaneousCTMRG; subspace = Z2Space(0 => 10)) + ) + @test_throws ArgumentError PEPSKit.check_input( + leading_boundary, InfiniteSquareNetwork(psi), env, + si_ctmrg(SimultaneousCTMRG; subspace = ℂ^10), + ) + @test_throws ArgumentError PEPSKit.check_input( + leading_boundary, InfiniteSquareNetwork(psi), CTMRGEnv(adapt(AT, InfinitePEPS(ℂ^2, ℂ^2)), ℂ^2), alg, + ) + + @test_throws ArgumentError PEPSKit._check_algorithm_combination( + SimultaneousCTMRG(; projector_alg = (; alg = :SubspaceIterationProjector, subspace = ℂ^7)), + ImplicitGradient(), + ) + end +end + +"""Check subspace iteration for the requested storage type.""" +function ctmrg_subspace_iteration(AT) + return @testset "Subspace iteration ($AT)" begin + ctmrg_subspace_iteration_configuration(AT) + ctmrg_subspace_iteration_convergence(AT) + ctmrg_subspace_iteration_half_environments(AT) + ctmrg_subspace_iteration_rangefinders(AT) + ctmrg_subspace_iteration_phase_caches(AT) + ctmrg_subspace_iteration_unequal_spaces(AT) + ctmrg_subspace_iteration_unitcells(AT) + ctmrg_subspace_iteration_symmetry_sectors(AT) + ctmrg_subspace_iteration_symmetric_unitcells(AT) + ctmrg_subspace_iteration_compatibility(AT) + end +end diff --git a/test/testsuite/ctmrg/subspace_updates.jl b/test/testsuite/ctmrg/subspace_updates.jl new file mode 100644 index 000000000..c6fecaec4 --- /dev/null +++ b/test/testsuite/ctmrg/subspace_updates.jl @@ -0,0 +1,179 @@ +using Test +using Random +using LinearAlgebra +using MatrixAlgebraKit +using TensorKit +using PEPSKit +using Adapt +using MatrixAlgebraKit: TruncatedAlgorithm + +"""Count the orthogonalizations performed by a production range update.""" +struct SubspaceCountingOrth{A} <: MatrixAlgebraKit.AbstractAlgorithm + alg::A + qr::Base.RefValue{Int} + lq::Base.RefValue{Int} +end + +"""Count QR calls while using the ordinary forward factorization.""" +function MatrixAlgebraKit.left_orth!(t::AbstractTensorMap, a::SubspaceCountingOrth) + a.qr[] += 1 + return left_orth!(t, a.alg) +end + +"""Count LQ calls while using the ordinary forward factorization.""" +function MatrixAlgebraKit.right_orth!(t::AbstractTensorMap, a::SubspaceCountingOrth) + a.lq[] += 1 + return right_orth!(t, a.alg) +end + +"""Evaluate the defining recurrences using an explicit full environment.""" +function expected_subspace_ranges(M, X, Y, update, q) + if update == :independent + for _ in 0:q + X = first(left_orth(M' * M * X)) + Y = first(left_orth(M * M' * Y)) + end + else + for _ in 1:q + X = first(left_orth(M' * M * X)) + end + Q = first(left_orth(M * X)) + X = first(left_orth(M' * Q)) + Y = update == :alternating_refined ? first(left_orth(M * X)) : Q + end + return X, Y +end + +"""Check recurrences for the subspace iteration projector.""" +function ctmrg_subspace_updates_recurrences(AT) + return @testset "Range recurrences and factorization counts" begin + @testset "$T, $update, q=$q" for T in (Float64, ComplexF64), update in (:independent, :alternating, :alternating_refined), q in (0, 1) + Random.seed!(264) + L, R = adapt(AT, randn(T, ℂ^9 ← ℂ^9)), adapt(AT, randn(T, ℂ^9 ← ℂ^7)) + X, Y = first(left_orth(adapt(AT, randn(T, ℂ^7 ← ℂ^4)))), last(right_orth(adapt(AT, randn(T, ℂ^4 ← ℂ^9))))' + counter = SubspaceCountingOrth(OrthAdjoint().fwd_alg, Ref(0), Ref(0)) + alg = SubspaceIterationProjector(; + subspace = ℂ^4, update, n_subspace_iter = q, + rotate_basis = false, orth_alg = OrthAdjoint(; fwd_alg = counter) + ) + before = deepcopy((X, Y)) + trunc = truncspace(ℂ^2) + result = PEPSKit._subspace_iteration_decomposition(L, R, alg, (; X, Y, trunc)) + expected_count = update == :alternating_refined ? 2(q + 1) + 1 : 2(q + 1) + @test counter.qr[] + counter.lq[] == expected_count + if update == :alternating + @test counter.lq[] == 1 + end + @test (X, Y) == before + Xref, Yref = expected_subspace_ranges(L * R, X, Y, update, q) + @test result.X * result.X' ≈ Xref * Xref' atol = 1.0e-10 + Ynew = PEPSKit.rangefinder_update(X, Y, L, R, alg.subspace, alg).Y + @test Ynew * Ynew' ≈ Yref * Yref' atol = 1.0e-10 + @test Ynew' * Ynew ≈ id(T, domain(Ynew)) atol = 1.0e-10 + @test isnothing(result.Y) == (update != :independent) + core = Yref' * L * R * Xref + u, s, v, = svd_trunc(core, TruncatedAlgorithm(DefaultAlgorithm(), trunc)) + expected = Yref * u * s * v * Xref' + actual = result.U * result.S * result.V + @test actual / norm(actual) ≈ expected / norm(expected) atol = 1.0e-10 + + end + end +end + +"""Check complete sectors for the subspace iteration projector.""" +function ctmrg_subspace_updates_complete_sectors(AT) + return @testset "Complete and rectangular sectors" begin + @testset "$update, $m ← $n, q=$q" for update in (:independent, :alternating, :alternating_refined), (m, n) in ((7, 5), (5, 7), (5, 5)), q in (0, 1) + Random.seed!(265) + M = adapt(AT, randn(ComplexF64, ℂ^m ← ℂ^n)) + alg = SubspaceIterationProjector(; subspace = ℂ^20, update, n_subspace_iter = q, rotate_basis = false) + trunc = truncspace(ℂ^2) + result = PEPSKit._subspace_iteration_decomposition( + M, id(ℂ^n), alg, + (; X = nothing, Y = nothing, trunc) + ) + u, s, v, = svd_trunc(M, TruncatedAlgorithm(DefaultAlgorithm(), trunc)) + actual, expected = result.U * result.S * result.V, u * s * v + @test actual / norm(actual) ≈ expected / norm(expected) atol = 1.0e-10 + @test result.X' * result.X ≈ id(ComplexF64, domain(result.X)) atol = 1.0e-10 + + end + @testset "$S, $update" for (S, charges) in ( + (U1Irrep, (0, 1)), + (SU2Irrep, (0, 1 // 2)), + (A4Irrep, (0, 3)), + (FermionParity ⊠ U1Irrep, ((0, 0), (1, 1))), + (FibonacciAnyon, (:I, :τ)), + ), update in (:independent, :alternating, :alternating_refined) + W = Vect[S](charges[1] => 3, charges[2] => 7) + K = Vect[S](charges[1] => 3, charges[2] => 4) + χ = Vect[S](charges[1] => 2, charges[2] => 2) + Random.seed!(266) + A = adapt(AT, zeros(ComplexF64, W ← W)) + for (_, b) in blocks(A) + b[diagind(b)] .= 0.1 .^ (0:(size(b, 1) - 1)) + end + A = first(left_orth(adapt(AT, randn(ComplexF64, space(A))))) * A * last(right_orth(adapt(AT, randn(ComplexF64, space(A))))) + alg = SubspaceIterationProjector(; subspace = K, update) + state = (; X = nothing, Y = nothing, trunc = truncspace(χ)) + local result + for _ in 1:5 + result = PEPSKit._subspace_iteration_decomposition(A, id(ComplexF64, W), alg, state) + state = (; result.X, result.Y, state.trunc) + end + u, s, v, = svd_trunc(A, TruncatedAlgorithm(DefaultAlgorithm(), state.trunc)) + actual, expected = result.U * result.S * result.V, u * s * v + @test actual / norm(actual) ≈ expected / norm(expected) atol = 1.0e-8 + @test domain(result.X) == ProductSpace(K) + @test result.U' * result.U ≈ id(ComplexF64, χ) atol = 1.0e-10 + end + end +end + +"""Check recovery for the subspace iteration projector.""" +function ctmrg_subspace_updates_recovery(AT) + return @testset "Cache configuration and scalar recovery" begin + @testset "$flavor" for flavor in (SimultaneousCTMRG, SequentialCTMRG) + Random.seed!(268) + ψ = adapt(AT, InfinitePEPS(ℂ^2, ℂ^2)) + network, env = InfiniteSquareNetwork(ψ), CTMRGEnv(ψ, ℂ^3) + alternating_alg = flavor(0.0, 1, 1, 0, SubspaceIterationProjector(; subspace = ℂ^5, update = :alternating)) + independent_alg = flavor(0.0, 1, 1, 0, SubspaceIterationProjector(; subspace = ℂ^5)) + cache = PEPSKit.initialize_projector_cache(network, env, alternating_alg) + _, _, reset = PEPSKit.ctmrg_iteration(network, env, independent_alg, cache) + @test !isnothing(reset.Y) + real_ψ = adapt(AT, InfinitePEPS(randn, Float64, ℂ^2, ℂ^2)) + real_network, real_env = InfiniteSquareNetwork(real_ψ), CTMRGEnv(real_ψ, ℂ^3) + real_cache = PEPSKit.initialize_projector_cache(real_network, real_env, alternating_alg) + _, _, reset = PEPSKit.ctmrg_iteration(network, env, alternating_alg, real_cache) + @test scalartype(first(reset.X)) == ComplexF64 + end + end +end + +"""Check initialization for the subspace iteration projector.""" +function ctmrg_subspace_updates_initialization(AT) + return @testset "Seeded initialization" begin + @testset "$update" for update in (:independent, :alternating, :alternating_refined) + a, b = [SubspaceIterationProjector(; subspace = ℂ^4, update, rng = Random.Xoshiro(269)) for _ in 1:2] + first_bases = PEPSKit.rangefinder_initialize(ComplexF64, ProductSpace(ℂ^7), ProductSpace(ℂ^6), a.subspace, a) + second_bases = PEPSKit.rangefinder_initialize(ComplexF64, ProductSpace(ℂ^7), ProductSpace(ℂ^6), b.subspace, b) + @test first_bases == second_bases + @test first_bases[1]' * first_bases[1] ≈ id(ComplexF64, ℂ^4) atol = 1.0e-12 + if !isnothing(first_bases[2]) + @test first_bases[2]' * first_bases[2] ≈ id(ComplexF64, ℂ^4) atol = 1.0e-12 + end + end + end +end + +"""Check subspace range updates for the requested storage type.""" +function ctmrg_subspace_updates(AT) + return @testset "Subspace range updates ($AT)" begin + ctmrg_subspace_updates_recurrences(AT) + ctmrg_subspace_updates_complete_sectors(AT) + ctmrg_subspace_updates_recovery(AT) + ctmrg_subspace_updates_initialization(AT) + end +end diff --git a/test/testsuite/gradients/subspace_dimensions.jl b/test/testsuite/gradients/subspace_dimensions.jl new file mode 100644 index 000000000..eaeab1cda --- /dev/null +++ b/test/testsuite/gradients/subspace_dimensions.jl @@ -0,0 +1,33 @@ +using Test +using Random +using LinearAlgebra +using TensorKit +using PEPSKit +using Zygote +using Adapt +using ChainRulesCore: ignore_derivatives + +"""Check cache feedback through active and charge-refresh sweeps.""" +function gradients_subspace_dimensions(AT) + return @testset "Adaptive subspace gradients ($AT)" begin + @testset "$update, q=$q, $mode, probe=$probe" for update in (:independent, :alternating, :alternating_refined), q in (0, 1), mode in (:active, :refresh), probe in (nothing, 1) + Random.seed!(275) + W = U1Space(0 => 6, 1 => 6, 2 => 6) ⊗ U1Space(0 => 1) + A, direction, target = (adapt(AT, randn(ComplexF64, W ← W)) for _ in 1:3) + direction /= norm(direction) + alg = SubspaceIterationProjector(; subspace = 8, trunc = truncrank(4), charge_probe = probe, update, n_subspace_iter = q, rotate_basis = q == 0) + X, Y = PEPSKit.rangefinder_initialize(A, id(W), fuse(W), alg) + function loss(t) + ignore_derivatives(() -> Random.seed!(276)) + first_step = PEPSKit._subspace_iteration_decomposition(t, id(W), alg, (; X, Y, trunc = alg.trunc)) + result = PEPSKit._subspace_iteration_decomposition(t, id(W), alg, (; first_step.X, first_step.Y, trunc = alg.trunc, refresh = mode == :refresh)) + return real(dot(target, result.U * result.S * result.V)) + end + gradient = only(Zygote.gradient(loss, A)) + h = 1.0e-6 + central(δ) = (loss(A + δ * direction) - loss(A - δ * direction)) / (2δ) + finite = (4central(h) - central(2h)) / 3 + @test real(dot(gradient, direction)) ≈ finite rtol = 3.0e-5 atol = 3.0e-7 + end + end +end diff --git a/test/testsuite/gradients/subspace_iteration.jl b/test/testsuite/gradients/subspace_iteration.jl new file mode 100644 index 000000000..6502175cb --- /dev/null +++ b/test/testsuite/gradients/subspace_iteration.jl @@ -0,0 +1,178 @@ +using Test +using Random +using LinearAlgebra +using MatrixAlgebraKit +using TensorKit +using Zygote +using PEPSKit +using Adapt +using ChainRulesCore: ignore_derivatives + +"""Check a directional derivative with Richardson-extrapolated central differences.""" +function si_directional_difference(f, x, direction; h = 1.0e-4) + central(δ) = (f(x + δ * direction) - f(x - δ * direction)) / (2δ) + return (4central(h) - central(2h)) / 3 +end + +"""Build a Hermitian nearest-neighbor operator compatible with the physical symmetry.""" +function si_test_operator(AT, T, P) + h = adapt(AT, randn(T, P ⊗ P ← P ⊗ P)) + h = h + h' + h /= norm(h) + return LocalOperator(fill(P, 1, 1), ((1, 1), (1, 2)) => h, ((1, 1), (2, 1)) => h) +end + +"""Supply small PEPS spaces spanning trivial, Abelian, non-Abelian, and fermionic sectors.""" +function si_energy_test_spaces() + F = Vect[FermionParity ⊠ U1Irrep]((0, 0) => 1, (1, 1) => 1, (1, -1) => 1) + A = Vect[A4Irrep](3 => 1) + return ( + (ℂ^2, ℂ^2, ℂ^2), + (U1Space(-1 => 1, 0 => 1, 1 => 1), U1Space(0 => 1, 1 => 1), U1Space(-1 => 1, 0 => 1, 1 => 1)), + (SU2Space(0 => 1, 1 => 1), SU2Space(1 // 2 => 1), SU2Space(0 => 1, 1 => 1)), + (F, F, F), + (A, A, Vect[A4Irrep](0 => 1, 3 => 1)), + ) +end + +"""Check orthogonalization for the subspace iteration projector.""" +function gradients_subspace_iteration_orthogonalization(AT) + return @testset "Right orthogonalization pullback" begin + spaces = ( + (ℂ^3, ℂ^6), + (U1Space(-1 // 2 => 2, 1 // 2 => 3), U1Space(-1 // 2 => 4, 1 // 2 => 5)), + (SU2Space(0 => 2, 1 // 2 => 3), SU2Space(0 => 4, 1 // 2 => 5)), + ( + Vect[FermionParity ⊠ U1Irrep]((0, 0) => 2, (1, 1) => 3), + Vect[FermionParity ⊠ U1Irrep]((0, 0) => 4, (1, 1) => 5), + ), + ) + @testset "$T, $K ← $W" for T in (Float64, ComplexF64), (K, W) in spaces + Random.seed!(731) + A, direction = adapt(AT, randn(T, K ← W)), adapt(AT, randn(T, K ← W)) + target = adapt(AT, randn(T, W ← W)) + loss(t) = real(dot(target, last(right_orth(t, OrthAdjoint()))' * last(right_orth(t, OrthAdjoint())))) + g = only(Zygote.gradient(loss, A)) + @test real(dot(g, direction)) ≈ si_directional_difference(loss, A, direction) rtol = 1.0e-6 atol = 1.0e-8 + end + end +end + +"""Check recycled basis for the subspace iteration projector.""" +function gradients_subspace_iteration_recycled_basis(AT) + return @testset "Recycled-basis pullback" begin + spaces = ( + (ℂ^6, ℂ^4, ℂ^2), + ( + U1Space(-1 // 2 => 4, 1 // 2 => 5, 3 // 2 => 3), + U1Space(-1 // 2 => 3, 1 // 2 => 3, 3 // 2 => 2), + U1Space(-1 // 2 => 2, 1 // 2 => 2, 3 // 2 => 1), + ), + ( + SU2Space(0 => 4, 1 // 2 => 5, 1 => 3), + SU2Space(0 => 3, 1 // 2 => 3, 1 => 2), SU2Space(0 => 2, 1 // 2 => 2, 1 => 1), + ), + ( + Vect[A4Irrep](0 => 4, 1 => 5, 3 => 3), + Vect[A4Irrep](0 => 3, 1 => 3, 3 => 2), Vect[A4Irrep](0 => 2, 1 => 2, 3 => 1), + ), + ( + Vect[FermionParity ⊠ U1Irrep]((0, 0) => 4, (1, 1) => 5, (1, -1) => 3), + Vect[FermionParity ⊠ U1Irrep]((0, 0) => 3, (1, 1) => 3, (1, -1) => 2), + Vect[FermionParity ⊠ U1Irrep]((0, 0) => 2, (1, 1) => 2, (1, -1) => 1), + ), + ( + Vect[FibonacciAnyon](:I => 4, :τ => 5), + Vect[FibonacciAnyon](:I => 3, :τ => 3), Vect[FibonacciAnyon](:I => 2, :τ => 2), + ), + ) + @testset "$W, $K, $χ" for (W, K, χ) in spaces + Random.seed!(732) + L, R = adapt(AT, randn(ComplexF64, W ← W)), adapt(AT, randn(ComplexF64, W ← W)) + direction, target = adapt(AT, randn(ComplexF64, space(L))), adapt(AT, randn(ComplexF64, space(L))) + alg = SubspaceIterationProjector(; subspace = K) + X, Y = PEPSKit.rangefinder_initialize(L, R, alg.subspace, alg) + inputs = (; X, Y, trunc = truncspace(χ)) + function loss(t) + first_step = PEPSKit._subspace_iteration_decomposition(t, R, alg, inputs) + next_inputs = (; first_step.X, first_step.Y, inputs.trunc) + result = PEPSKit._subspace_iteration_decomposition(t, R, alg, next_inputs) + return real(dot(target, result.U * result.S * result.V)) + end + g = only(Zygote.gradient(loss, L)) + @test real(dot(g, direction)) ≈ si_directional_difference(loss, L, direction) rtol = 1.0e-5 atol = 1.0e-7 + # A phase-sensitive loss checks the gauge convention separately from reconstruction. + target_U = adapt(AT, randn(ComplexF64, codomain(L) ← χ)) + phase_loss(t) = real(dot(target_U, PEPSKit._subspace_iteration_decomposition(t, R, alg, inputs).U)) + gphase = only(Zygote.gradient(phase_loss, L)) + @test real(dot(gphase, direction)) ≈ si_directional_difference(phase_loss, L, direction) rtol = 1.0e-5 atol = 1.0e-7 + end + end +end + +"""Check energy for the subspace iteration projector.""" +function gradients_subspace_iteration_energy(AT) + return @testset "Energy backpropagation" begin + @testset "$T, $flavor, $P, $V, $χ" for T in (Float64, ComplexF64), flavor in (SimultaneousCTMRG, SequentialCTMRG), (P, V, χ) in si_energy_test_spaces() + Random.seed!(733) + K = infimum(fuse(χ ⊗ V ⊗ V'), χ ⊕ χ) + psi = adapt(AT, InfinitePEPS(randn, T, P, V)) + psi = psi / norm(psi) + direction = adapt(AT, InfinitePEPS(randn, T, P, V)) + direction = direction / norm(direction) + env = CTMRGEnv(randn, T, psi, χ) + H = si_test_operator(AT, T, P) + alg = flavor(; + miniter = 2, maxiter = 2, tol = 0.0, verbosity = 0, + projector_alg = (; alg = :SubspaceIterationProjector, subspace = K), + ) + function loss(ψ) + ignore_derivatives(() -> Random.seed!(734)) + return cost_function(ψ, first(leading_boundary(env, ψ, alg)), H) + end + g = only(Zygote.gradient(loss, psi)) + @test real(dot(g, direction)) ≈ si_directional_difference(loss, psi, direction; h = 1.0e-6) rtol = 1.0e-5 atol = 1.0e-7 + end + end +end + +"""Check fixed point for the subspace iteration projector.""" +function gradients_subspace_iteration_fixed_point(AT) + return @testset "Fixed-point gradient" begin + @testset "$flavor, $P, $V, $χ" for flavor in (SimultaneousCTMRG, SequentialCTMRG), (P, V, χ) in si_energy_test_spaces() + Random.seed!(738) + K = fuse(χ ⊗ V ⊗ V') + psi = adapt(AT, InfinitePEPS(P, V)) + psi = psi / norm(psi) + direction = adapt(AT, InfinitePEPS(P, V)) + direction = direction / norm(direction) + env = CTMRGEnv(psi, χ) + H = si_test_operator(AT, ComplexF64, P) + alg = flavor(; + miniter = 20, maxiter = 200, tol = 1.0e-12, verbosity = 0, + projector_alg = (; alg = :SubspaceIterationProjector, subspace = K), + ) + gradmode = FixedPointGradient(; solver_alg = (; alg = :GMRES, tol = 1.0e-10, maxiter = 100)) + function loss(ψ) + ignore_derivatives(() -> Random.seed!(736)) + env′, = PEPSKit.hook_pullback(leading_boundary, env, ψ, alg; alg_rrule = gradmode) + return cost_function(ψ, env′, H) + end + Random.seed!(736) + _, info = leading_boundary(env, psi, alg) + @test info.converged + g = only(Zygote.gradient(loss, psi)) + @test real(dot(g, direction)) ≈ si_directional_difference(loss, psi, direction) rtol = 1.0e-4 atol = 1.0e-7 + end + end +end + +"""Check subspace iteration gradients for the requested storage type.""" +function gradients_subspace_iteration(AT) + return @testset "Subspace iteration gradients ($AT)" begin + gradients_subspace_iteration_orthogonalization(AT) + gradients_subspace_iteration_recycled_basis(AT) + gradients_subspace_iteration_energy(AT) + gradients_subspace_iteration_fixed_point(AT) + end +end diff --git a/test/testsuite/gradients/subspace_updates.jl b/test/testsuite/gradients/subspace_updates.jl new file mode 100644 index 000000000..b2595c820 --- /dev/null +++ b/test/testsuite/gradients/subspace_updates.jl @@ -0,0 +1,175 @@ +using Test +using Random +using LinearAlgebra +using MatrixAlgebraKit +using TensorKit +using Zygote +using PEPSKit +using Adapt +using ChainRulesCore: ignore_derivatives + +"""Estimate a directional derivative using extrapolated central differences.""" +function subspace_update_difference(f, x, direction; h = 1.0e-6) + central(δ) = (f(x + δ * direction) - f(x - δ * direction)) / (2δ) + return (4central(h) - central(2h)) / 3 +end + +"""Contract two production sweeps with the first sweep's cache retained in the graph.""" +function subspace_two_sweep_energy(ψ, env, alg, cache, H) + network = InfiniteSquareNetwork(ψ) + env1, _, cache1 = PEPSKit.ctmrg_iteration(network, env, alg, cache) + env2, info, _ = PEPSKit.ctmrg_iteration(network, env1, alg, cache1) + return cost_function(ψ, env2, H) +end + +"""Check matrix for the subspace iteration projector.""" +function gradients_subspace_updates_matrix(AT) + return @testset "Two-step matrix pullbacks" begin + @testset "$T, $update, q=$q, rotate=$rotate" for T in (Float64, ComplexF64), update in (:independent, :alternating, :alternating_refined), q in (0, 1), rotate in (false, true) + Random.seed!(84264) + W, K, χ = ℂ^7, ℂ^4, ℂ^2 + A, direction, target = adapt(AT, randn(T, W ← W)), adapt(AT, randn(T, W ← W)), adapt(AT, randn(T, W ← W)) + direction /= norm(direction) + alg = SubspaceIterationProjector(; + subspace = K, update, n_subspace_iter = q, + rotate_basis = rotate, + ) + X, Y = PEPSKit.rangefinder_initialize(A, id(T, W), K, alg) + initial = (; X, Y, trunc = truncspace(χ)) + function loss(t) + step = PEPSKit._subspace_iteration_decomposition(t, id(T, W), alg, initial) + result = PEPSKit._subspace_iteration_decomposition( + t, id(T, W), alg, + (; step.X, step.Y, initial.trunc) + ) + return real(dot(target, result.U * result.S * result.V)) + end + gradient = only(Zygote.gradient(loss, A)) + @test real(dot(gradient, direction)) ≈ subspace_update_difference(loss, A, direction) rtol = 3.0e-5 atol = 3.0e-7 + target_U = adapt(AT, randn(T, W ← χ)) + phase_loss(t) = real(dot(target_U, PEPSKit._subspace_iteration_decomposition(t, id(T, W), alg, initial).U)) + phase_gradient = only(Zygote.gradient(phase_loss, A)) + @test real(dot(phase_gradient, direction)) ≈ subspace_update_difference(phase_loss, A, direction) rtol = 3.0e-5 atol = 3.0e-7 + end + end +end + +"""Check symmetric for the subspace iteration projector.""" +function gradients_subspace_updates_symmetric(AT) + return @testset "Symmetric recycled pullbacks" begin + F = FermionParity ⊠ U1Irrep + spaces = ( + (SU2Space(0 => 3, 1 // 2 => 7), SU2Space(0 => 3, 1 // 2 => 4), SU2Space(0 => 2, 1 // 2 => 2)), + (Vect[F]((0, 0) => 3, (1, 1) => 7), Vect[F]((0, 0) => 3, (1, 1) => 4), Vect[F]((0, 0) => 2, (1, 1) => 2)), + ) + @testset "$W, $update" for (W, K, χ) in spaces, update in (:independent, :alternating, :alternating_refined) + Random.seed!(84265) + A, direction, target = adapt(AT, randn(ComplexF64, W ← W)), adapt(AT, randn(ComplexF64, W ← W)), adapt(AT, randn(ComplexF64, W ← W)) + direction /= norm(direction) + alg = SubspaceIterationProjector(; subspace = K, update, n_subspace_iter = 1) + X, Y = PEPSKit.rangefinder_initialize(A, id(ComplexF64, W), K, alg) + initial = (; X, Y, trunc = truncspace(χ)) + function loss(t) + step = PEPSKit._subspace_iteration_decomposition(t, id(ComplexF64, W), alg, initial) + result = PEPSKit._subspace_iteration_decomposition( + t, id(ComplexF64, W), alg, + (; step.X, step.Y, initial.trunc) + ) + return real(dot(target, result.U * result.S * result.V)) + end + gradient = only(Zygote.gradient(loss, A)) + @test real(dot(gradient, direction)) ≈ subspace_update_difference(loss, A, direction) rtol = 3.0e-5 atol = 3.0e-7 + end + end +end + +"""Check differentiation through consecutive production sweeps.""" +function gradients_subspace_updates_sweeps(AT) + return @testset "Two-sweep production pullbacks" begin + settings = ((:independent, 0), (:alternating_refined, 0), (:alternating, 0), (:alternating, 1)) + @testset "$update, q=$q, $flavor, rotate=$rotate" for (update, q) in settings, flavor in (SimultaneousCTMRG, SequentialCTMRG), rotate in (false, true) + Random.seed!(84266) + ψ = adapt(AT, InfinitePEPS(ℂ^2, ℂ^2)) + ψ /= norm(ψ) + direction = adapt(AT, InfinitePEPS(ℂ^2, ℂ^2)) + direction /= norm(direction) + env = CTMRGEnv(ψ, ℂ^3) + H = adapt(AT, heisenberg_XYZ(InfiniteSquare())) + projector = SubspaceIterationProjector(; + subspace = ℂ^4, update, n_subspace_iter = q, + rotate_basis = rotate + ) + alg = flavor(0.0, 2, 2, 0, projector) + cache = PEPSKit.initialize_projector_cache(InfiniteSquareNetwork(ψ), env, alg) + loss(x) = subspace_two_sweep_energy(x, env, alg, cache, H) + gradient = only(Zygote.gradient(loss, ψ)) + @test real(dot(gradient, direction)) ≈ subspace_update_difference(loss, ψ, direction; h = 1.0e-7) rtol = 3.0e-5 atol = 3.0e-7 + end + end +end + +"""Check fixed point for the subspace iteration projector.""" +function gradients_subspace_updates_fixed_point(AT) + references = Dict() + cases = [(update, flavor, unitcell, false) for update in (:independent, :alternating, :alternating_refined) for flavor in (SimultaneousCTMRG, SequentialCTMRG) for unitcell in ((1, 1), (2, 1))] + # Complete-range regressions for the two updates with drifting bond gauges. + append!(cases, ((update, SimultaneousCTMRG, (1, 1), true) for update in (:independent, :alternating))) + return @testset "Subspace fixed-point gradients" begin + @testset "$update, $flavor, $unitcell, complete=$complete" for (update, flavor, unitcell, complete) in cases + Random.seed!(complete ? 5 : 84267) + ψ = adapt(AT, InfinitePEPS(ℂ^2, ℂ^2; unitcell)) + ψ /= norm(ψ) + direction = adapt(AT, InfinitePEPS(ℂ^2, ℂ^2; unitcell)) + direction /= norm(direction) + env = CTMRGEnv(ψ, ℂ^4) + H = adapt(AT, heisenberg_XYZ(InfiniteSquare(unitcell...))) + alg = flavor(; + tol = 1.0e-11, miniter = 20, maxiter = 500, verbosity = 0, + projector_alg = (; + alg = :SubspaceIterationProjector, subspace = ℂ^(complete ? 16 : 6), update, + n_subspace_iter = update == :alternating ? 1 : 0, + ) + ) + mode = FixedPointGradient(; solver_alg = (; alg = :GMRES, tol = 1.0e-10, maxiter = 100)) + function loss(x) + ignore_derivatives(() -> Random.seed!(84268)) + boundary, = PEPSKit.hook_pullback(leading_boundary, env, x, alg; alg_rrule = mode) + return cost_function(x, boundary, H) + end + _, info = leading_boundary(env, ψ, alg) + @test info.converged + gradient = only(Zygote.gradient(loss, ψ)) + reference = get!(references, (flavor, unitcell, complete)) do + full = flavor(; tol = 1.0e-12, maxiter = 500, verbosity = 0, projector_alg = :FullInfiniteProjector) + f(x) = cost_function(x, first(PEPSKit.hook_pullback(leading_boundary, env, x, full; alg_rrule = mode)), H) + only(Zygote.gradient(f, ψ)) + end + @test norm(gradient - reference) / norm(reference) < 1.0e-7 + if !complete && update == :independent && flavor == SimultaneousCTMRG && unitcell == (1, 1) + for solver in (:Arnoldi, :BiCGStab, :GeomSum, :ManualIter) + other = FixedPointGradient(; solver_alg = (; alg = solver, tol = 1.0e-10, maxiter = 200)) + f(x) = cost_function(x, first(PEPSKit.hook_pullback(leading_boundary, env, x, alg; alg_rrule = other)), H) + @test norm(only(Zygote.gradient(f, ψ)) - reference) / norm(reference) < 1.0e-7 + end + end + if complete && update == :independent + function unrolled(x) + ignore_derivatives(() -> Random.seed!(84268)) + return cost_function(x, first(leading_boundary(env, x, alg)), H) + end + @test norm(only(Zygote.gradient(unrolled, ψ)) - reference) / norm(reference) < 1.0e-7 + end + @test real(dot(gradient, direction)) ≈ subspace_update_difference(loss, ψ, direction; h = 1.0e-4) rtol = 1.0e-7 atol = 1.0e-9 + end + end +end + +"""Check subspace update gradients for the requested storage type.""" +function gradients_subspace_updates(AT) + return @testset "Subspace update gradients ($AT)" begin + gradients_subspace_updates_matrix(AT) + gradients_subspace_updates_symmetric(AT) + gradients_subspace_updates_sweeps(AT) + gradients_subspace_updates_fixed_point(AT) + end +end