diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index b57b6361f..a9bab8511 100644 --- a/src/PEPSKit.jl +++ b/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") diff --git a/src/algorithms/contractions/ctmrg/projector.jl b/src/algorithms/contractions/ctmrg/projector.jl index f8791d77f..0af18b268 100644 --- a/src/algorithms/contractions/ctmrg/projector.jl +++ b/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 diff --git a/src/algorithms/ctmrg/projectors/projectors.jl b/src/algorithms/ctmrg/projectors/projectors.jl index 6815669f8..76c4ddf6a 100644 --- a/src/algorithms/ctmrg/projectors/projectors.jl +++ b/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) diff --git a/src/algorithms/ctmrg/projectors/randomized.jl b/src/algorithms/ctmrg/projectors/randomized.jl new file mode 100644 index 000000000..afaa27d54 --- /dev/null +++ b/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 diff --git a/src/algorithms/ctmrg/simultaneous.jl b/src/algorithms/ctmrg/simultaneous.jl index ecaf52c7b..fcfee42c5 100644 --- a/src/algorithms/ctmrg/simultaneous.jl +++ b/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) diff --git a/src/environments/ctmrg_environments.jl b/src/environments/ctmrg_environments.jl index 1b9308244..f5f57ca78 100644 --- a/src/environments/ctmrg_environments.jl +++ b/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')