Add randomized full-environment CTMRG projector

src/PEPSKit.jl

@@ -6,6 +6,7 @@ using Compat
 using Accessors: @set, @reset
 using VectorInterface
 import VectorInterface as VI
+import Random
 
 using MatrixAlgebraKit
 using MatrixAlgebraKit:
@@ -128,6 +129,7 @@ include("algorithms/ctmrg/ctmrg.jl")
 include("algorithms/ctmrg/projectors/projectors.jl")
 include("algorithms/ctmrg/projectors/halfinfinite.jl")
 include("algorithms/ctmrg/projectors/fullinfinite.jl")
+include("algorithms/ctmrg/projectors/randomized.jl")
 include("algorithms/ctmrg/projectors/c4v_eigh.jl")
 include("algorithms/ctmrg/projectors/c4v_qr.jl")
 include("algorithms/ctmrg/simultaneous.jl")

src/algorithms/contractions/ctmrg/projector.jl

@@ -244,3 +244,14 @@ function contract_projectors(U, S, V, Q, Q_next)
     P_right = isqS * U' * Q
     return P_left, P_right
 end
+
+function contract_projectors(U, S, V, fq::FourQuadrants)
+    isqS = sdiag_pow(S, -0.5)
+    # use * to respect fermionic case
+    p1 = (codomainind(fq), domainind(fq))
+    p2 = (codomainind(fq), (numout(fq) + 1,))
+    P_left = tensorcontract(fq.Q3, p1, false, fq.Q4 * V' * isqS, p2, false, p2)
+    p3 = ((1,), codomainind(fq) .+ 1)
+    P_right = tensorcontract(isqS * U' * fq.Q1, p3, false, fq.Q2, p1, false, p3)
+    return P_left, P_right
+end

src/algorithms/ctmrg/projectors/projectors.jl

@@ -25,7 +25,7 @@ function ProjectorAlgorithm(;
 
     # parse SVD forward & rrule algorithm
 
-    decomposition_algorithm = if alg in [:HalfInfiniteProjector, :FullInfiniteProjector]
+    decomposition_algorithm = if alg in [:HalfInfiniteProjector, :FullInfiniteProjector, :randomized]
         _alg_or_nt(SVDAdjoint, decomposition_alg)
     elseif alg in [:C4vEighProjector]
         _alg_or_nt(EighAdjoint, decomposition_alg)

src/algorithms/ctmrg/projectors/randomized.jl

@@ -0,0 +1,95 @@
+"""Full-environment randomized SVD projector from PR 263."""
+struct RandomizedProjector{S, T, R} <: ProjectorAlgorithm
+    decomposition_alg::S
+    trunc::T
+    rng::R
+    oversampling::Int
+    max_full::Int
+    n_subspace_iter::Int
+    verbosity::Int
+end
+
+PROJECTOR_SYMBOLS[:randomized] = RandomizedProjector
+
+decomposition_algorithm(alg::RandomizedProjector) = alg.decomposition_alg
+
+_default_randomized_oversampling = 10
+_default_randomized_max_full = 100
+_default_n_subspace_iter = 2
+
+# needed as default interface in PEPSKit.ProjectorAlgorithm
+"""Construct the randomized projector with the original PR defaults."""
+function RandomizedProjector(decomposition_alg, trunc, verbosity)
+    return RandomizedProjector(
+        decomposition_alg, trunc, Random.default_rng(), _default_randomized_oversampling,
+        _default_randomized_max_full, _default_n_subspace_iter, verbosity
+    )
+end
+
+"""Choose exact or oversampled sector dimensions using the previous environment space."""
+function random_domain(alg::RandomizedProjector, full_space, last_space)
+    sector_dims = map(sectors(full_space)) do s
+        if dim(full_space, s) <= alg.max_full
+            n = dim(full_space, s)
+        else
+            n = dim(last_space, s) + alg.oversampling
+        end
+        return s => n
+    end
+    return Vect[sectortype(last_space)](sector_dims)
+end
+
+
+"""Find each block range using exact bases for small blocks and randomized iteration otherwise."""
+function randomized_range_finder(A::AbstractTensorMap, alg::RandomizedProjector, randomized_space)
+    Q = TensorMap{eltype(A)}(undef, domain(A) ← randomized_space)
+    foreach(blocks(Q)) do (s, b)
+        Aad = A'
+        m, n = size(b)
+        if m <= alg.max_full
+            b .= LinearAlgebra.I(m)
+        else
+            Ω = randn(alg.rng, eltype(A), domain(A) ← Vect[sectortype(A)](s => n))
+            for _ in 1:alg.n_subspace_iter
+                Y = A * Ω
+                Qs, _ = left_orth!(Y)
+                Ω, _ = left_orth!(Aad * Qs)
+            end
+            Y = A * Ω
+            Qs, _ = left_orth!(Y)
+            b .= block(Qs, s)
+        end
+    end
+    return Q
+end
+
+# impose full env, could also be defined for half_infinite_environment, little gain
+"""Build the CTMRG projectors from the lazy randomized full-environment SVD."""
+function compute_projector(fq::FourQuadrants, last_space, alg::RandomizedProjector)
+    full_space = fuse(domain(fq))
+    randomized_space = random_domain(alg, full_space, last_space)
+    Q = randomized_range_finder(fq, alg, randomized_space)
+    B = Q' * fq
+    normalize!(B)  # TODO better way?
+
+    decomposition_alg = decomposition_algorithm(alg)
+    U′, S, V, truncation_error = svd_trunc!(B, decomposition_alg)
+    U = Q * U′
+    foreach(blocks(S)) do (s, b)
+        if size(b, 1) == dim(randomized_space, s) && size(b, 1) < dim(full_space, s)
+            @warn("Sector is too small, kept all computed values: ", s)
+        end
+    end
+
+    # Check for degenerate singular values; still needed for exact blocks
+    Zygote.isderiving() && ignore_derivatives() do
+        if alg.verbosity > 0 && is_degenerate_spectrum(S)
+            svals = TensorKit.SectorDict(c => diag(b) for (c, b) in blocks(S))
+            @warn("degenerate singular values detected: ", svals)
+        end
+    end
+
+    truncation_error = truncation_error / norm(S) # normalize truncation error
+    P_left, P_right = contract_projectors(U, S, V, fq)
+    return (P_left, P_right), (; U, S, V, truncation_error)
+end

src/algorithms/ctmrg/simultaneous.jl

@@ -117,6 +117,27 @@ function simultaneous_projectors(
     return compute_projector(ec, alg′)
 end
 
+"""Compute the randomized projector with the current truncation and previous sector dimensions."""
+function simultaneous_projectors(
+        coordinate, enlarged_corners::Array{E, 3}, env, alg::RandomizedProjector
+    ) where {E}
+    coordinate′ = _next_coordinate(coordinate, size(env)[2:3]...)
+    trunc = truncation_strategy(alg, env.edges[coordinate[1], coordinate′[2:3]...])
+    alg′ = _set_decomposition_truncation(alg, trunc)
+    rowsize, colsize = size(enlarged_corners)[2:3]
+    coordinate2 = _next_coordinate(coordinate, rowsize, colsize)
+    coordinate3 = _next_coordinate(coordinate2, rowsize, colsize)
+    coordinate4 = _next_coordinate(coordinate3, rowsize, colsize)
+    fq = FourQuadrants(
+        enlarged_corners[coordinate4...],
+        enlarged_corners[coordinate...],
+        enlarged_corners[coordinate2...],
+        enlarged_corners[coordinate3...],
+    )
+    last_space = space(env.edges[coordinate[1], coordinate′[2:3]...], 1)
+    return compute_projector(fq, last_space, alg′)
+end
+
 """
 $(SIGNATURES)
 

src/environments/ctmrg_environments.jl

@@ -480,3 +480,34 @@ function VI.inner(env₁::CTMRGEnv, env₂::CTMRGEnv)
     return inner((env₁.corners, env₁.edges), (env₂.corners, env₂.edges))
 end
 VI.norm(env::CTMRGEnv) = norm((env.corners, env.edges))
+
+# ==========================================================================================
+
+struct FourQuadrants{S, T, N, TM} <: AbstractTensorMap{S, T, N, N}
+    Q1::TM
+    Q2::TM
+    Q3::TM
+    Q4::TM
+
+    function FourQuadrants{T, S, N, TM}(
+            Q1, Q2, Q3, Q4
+        ) where {S, T, N, TM <: AbstractTensorMap{T, S, N, N}}
+        return new{T, S, N, TM}(Q1, Q2, Q3, Q4)
+    end
+end
+
+function FourQuadrants(Q1, Q2, Q3, Q4)
+    return FourQuadrants{eltype(Q1), spacetype(Q1), numin(Q1), typeof(Q1)}(Q1, Q2, Q3, Q4)
+end
+
+TensorKit.TensorMap(fq::FourQuadrants) = fq.Q1 * fq.Q2 * fq.Q3 * fq.Q4
+
+TensorKit.space(fq::FourQuadrants) = codomain(fq.Q1) ← domain(fq.Q4)
+
+# TBD use tensorcontract to handle fermions?
+(fq::FourQuadrants)(m) = fq.Q1 * (fq.Q2 * (fq.Q3 * (fq.Q4 * m)))
+
+Base.:*(m::AbstractTensorMap, fq::FourQuadrants) = m * fq.Q1 * fq.Q2 * fq.Q3 * fq.Q4
+Base.:*(fq::FourQuadrants, m::AbstractTensorMap) = fq.Q1 * (fq.Q2 * (fq.Q3 * (fq.Q4 * m)))
+
+Base.adjoint(fq::FourQuadrants) = FourQuadrants(fq.Q4', fq.Q3', fq.Q2', fq.Q1')