Support TruncationMultiplets

Project.toml

@@ -57,7 +57,7 @@ FiniteDifferences = "0.12"
 GPUArrays = "11.4.1"
 LRUCache = "1.6"
 LinearAlgebra = "1"
-MatrixAlgebraKit = "0.6.9"
+MatrixAlgebraKit = "0.6.13"
 Mooncake = "0.5.27"
 OhMyThreads = "0.8.0"
 PrecompileTools = "1"

ext/TensorKitGPUArraysExt.jl

@@ -85,7 +85,7 @@ function MatrixAlgebraKit.findtruncated_svd(values::GPUSectorVector, strategy::S
     return Adapt.adapt(Vector, MatrixAlgebraKit.findtruncated(values, strategy))
 end
 
-for strat in (:(MatrixAlgebraKit.TruncationByOrder), :(MatrixAlgebraKit.TruncationByError), :(MatrixAlgebraKit.TruncationIntersection), :(TensorKit.Factorizations.TruncationSpace))
+for strat in (:(MatrixAlgebraKit.TruncationByOrder), :(MatrixAlgebraKit.TruncationByError), :(MatrixAlgebraKit.TruncationIntersection), :(MatrixAlgebraKit.TruncationMultiplets), :(TensorKit.Factorizations.TruncationSpace))
     @eval function MatrixAlgebraKit.findtruncated_svd(values::GPUSectorVector, strategy::$strat)
         # returning a GPUSectorVector wrecks things in truncate_{co}domain
         # because of scalar indexing

src/TensorKit.jl

@@ -99,7 +99,7 @@ export @tensor, @tensoropt, @ncon, ncon, @planar, @plansor
 export scalar, add!, contract!
 
 # truncation schemes
-export notrunc, truncrank, trunctol, truncfilter, truncspace, truncerror
+export notrunc, truncrank, trunctol, truncfilter, truncspace, truncerror, truncmultiplets
 
 # cache management
 export empty_globalcaches!

src/factorizations/factorizations.jl

@@ -23,7 +23,8 @@ using MatrixAlgebraKit
 import MatrixAlgebraKit as MAK
 using MatrixAlgebraKit: AbstractAlgorithm, TruncatedAlgorithm, DiagonalAlgorithm
 using MatrixAlgebraKit: TruncationStrategy, NoTruncation, TruncationByValue, TruncationByError,
-    TruncationIntersection, TruncationUnion, TruncationByFilter, TruncationByOrder
+    TruncationIntersection, TruncationUnion, TruncationByFilter, TruncationByOrder,
+    TruncationMultiplets
 using MatrixAlgebraKit: diagview
 
 include("utility.jl")

src/factorizations/truncation.jl

@@ -337,6 +337,37 @@ function MAK.findtruncated_svd(values::SectorVector, strategy::TruncationUnion)
     return SectorDict{eltype(sectors), eltype(vals)}(sectors, vals)
 end
 
+# Degenerate values can occur in different blocks, for example when the sectors are related by
+# a symmetry, so groups of degenerate values are determined across all blocks simultaneously.
+function MAK.findtruncated(values::SectorVector, strategy::TruncationMultiplets)
+    return _findtruncated_multiplets(values, MAK.findtruncated(values, strategy.trunc), strategy)
+end
+function MAK.findtruncated_svd(values::SectorVector, strategy::TruncationMultiplets)
+    return _findtruncated_multiplets(values, MAK.findtruncated_svd(values, strategy.trunc), strategy)
+end
+
+function _findtruncated_multiplets(values::SectorVector, inds, strategy::TruncationMultiplets)
+    keep = similar(values, Bool)
+    to_keep = collect(parent(_ind_to_mask!(keep, inds)))
+    kept_values = collect(parent(values))
+    perm = sortperm(kept_values; by = abs, rev = true)
+    MAK._keep_multiplets!(to_keep, kept_values, perm, strategy)
+    copyto!(parent(keep), to_keep)
+    return keep
+end
+
+function _ind_to_mask!(keep::SectorVector{Bool}, inds::SectorVector{Bool})
+    copyto!(parent(keep), parent(inds))
+    return keep
+end
+function _ind_to_mask!(keep::SectorVector{Bool}, inds::SectorDict)
+    fill!(parent(keep), false)
+    for (c, ind) in pairs(inds)
+        keep[c][ind] .= true
+    end
+    return keep
+end
+
 # Truncation error
 # ----------------
 MAK.truncation_error(values::SectorVector, ind) = MAK.truncation_error!(copy(values), ind)

test/amd/factorizations.jl

@@ -434,6 +434,31 @@ for V in spacelist
                 @test norm(t - U5 * S5 * Vᴴ5) ≈ ϵ5 atol = eps(real(T))^(4 / 5)
                 @test minimum(diagview(S5)) >= λ
                 @test abs(dim(domain(S5)) - nvals) ≤ maximum(c -> blockdim(domain(t), c), blocksectors(t); init = 1)
+
+                U0, S0, Vᴴ0 = svd_compact(t, alg)
+                Sd = similar(S0)
+                copyto!(Sd.data, real(T)[3 - mod(i, 3) for i in 0:(length(Sd.data) - 1)])
+                td = U0 * Sd * Vᴴ0
+                rtol = sqrt(eps(real(T)))
+                _, Sr, _, _ = svd_trunc(td; trunc = truncrank(nvals), alg)
+                for mode in (:down, :up)
+                    trunc = truncmultiplets(truncrank(nvals); rtol, mode)
+                    U6, S6, Vᴴ6, ϵ6 = @constinferred svd_trunc(td; trunc, alg)
+                    @test td * Vᴴ6' ≈ U6 * S6
+                    @test isisometric(U6)
+                    @test isisometric(Vᴴ6; side = :right)
+                    @test norm(td - U6 * S6 * Vᴴ6) ≈ ϵ6 atol = eps(real(T))^(4 / 5)
+                    if mode == :down
+                        @test space(S6, 1) ≾ space(Sr, 1)
+                    else
+                        @test space(Sr, 1) ≾ space(S6, 1)
+                    end
+                    if dim(domain(S6)) > 0
+                        λ6 = minimum(diagview(S6))
+                        _, S7, _, _ = svd_trunc(td; trunc = trunctol(; atol = λ6 - 1 // 2), alg)
+                        @test space(S6, 1) == space(S7, 1)
+                    end
+                end
             end
         end
 

test/cuda/factorizations.jl

@@ -453,6 +453,31 @@ for V in spacelist
                 @test norm(t - U5 * S5 * Vᴴ5) ≈ ϵ5 atol = eps(real(T))^(4 / 5)
                 @test minimum(diagview(S5)) >= λ
                 @test abs(dim(domain(S5)) - nvals) ≤ maximum(c -> blockdim(domain(t), c), blocksectors(t); init = 1)
+
+                U0, S0, Vᴴ0 = svd_compact(t)
+                Sd = similar(S0)
+                copyto!(Sd.data, real(T)[3 - mod(i, 3) for i in 0:(length(Sd.data) - 1)])
+                td = U0 * Sd * Vᴴ0
+                rtol = sqrt(eps(real(T)))
+                _, Sr, _, _ = svd_trunc(td; trunc = truncrank(nvals))
+                for mode in (:down, :up)
+                    trunc = truncmultiplets(truncrank(nvals); rtol, mode)
+                    U6, S6, Vᴴ6, ϵ6 = @constinferred svd_trunc(td; trunc)
+                    @test td * Vᴴ6' ≈ U6 * S6
+                    @test isisometric(U6)
+                    @test isisometric(Vᴴ6; side = :right)
+                    @test norm(td - U6 * S6 * Vᴴ6) ≈ ϵ6 atol = eps(real(T))^(4 / 5)
+                    if mode == :down
+                        @test space(S6, 1) ≾ space(Sr, 1)
+                    else
+                        @test space(Sr, 1) ≾ space(S6, 1)
+                    end
+                    if dim(domain(S6)) > 0
+                        λ6 = minimum(diagview(S6))
+                        _, S7, _, _ = svd_trunc(td; trunc = trunctol(; atol = λ6 - 1 // 2))
+                        @test space(S6, 1) == space(S7, 1)
+                    end
+                end
             end
         end
 

test/factorizations/svd.jl

@@ -290,6 +290,31 @@ for V in spacelist
                 @test norm(t - U5 * S5 * Vᴴ5) ≈ ϵ5 atol = eps(real(T))^(4 / 5)
                 @test minimum(diagview(S5)) >= λ
                 @test abs(dim(domain(S5)) - nvals) ≤ maximum(c -> blockdim(domain(t), c), blocksectors(t); init = 1)
+
+                U0, S0, Vᴴ0 = svd_compact(t)
+                Sd = similar(S0)
+                copyto!(Sd.data, real(T)[3 - mod(i, 3) for i in 0:(length(Sd.data) - 1)])
+                td = U0 * Sd * Vᴴ0
+                rtol = sqrt(eps(real(T)))
+                _, Sr, _, _ = svd_trunc(td; trunc = truncrank(nvals))
+                for mode in (:down, :up)
+                    trunc = truncmultiplets(truncrank(nvals); rtol, mode)
+                    U6, S6, Vᴴ6, ϵ6 = @constinferred svd_trunc(td; trunc)
+                    @test td * Vᴴ6' ≈ U6 * S6
+                    @test isisometric(U6)
+                    @test isisometric(Vᴴ6; side = :right)
+                    @test norm(td - U6 * S6 * Vᴴ6) ≈ ϵ6 atol = eps(real(T))^(4 / 5)
+                    if mode == :down
+                        @test space(S6, 1) ≾ space(Sr, 1)
+                    else
+                        @test space(Sr, 1) ≾ space(S6, 1)
+                    end
+                    if dim(domain(S6)) > 0
+                        λ6 = minimum(diagview(S6))
+                        _, S7, _, _ = svd_trunc(td; trunc = trunctol(; atol = λ6 - 1 // 2))
+                        @test space(S6, 1) == space(S7, 1)
+                    end
+                end
             end
         end
     end

test/setup.jl

@@ -362,8 +362,8 @@ VIBMRepA4 = (
 for V in (Vtr, VRepℤ₂, VRepℤ₃, VRepU₁, VfRepℤ₂, VTwistedVecℤ₄, VTwistedVecℤ₃, VRepCU₁, VRepSU₂, VfRepCU₁, VfRepSU₂, VfHubbard, Vfib, VRepA4, VfRepA4, VRepA4Twistedℤ₄, VIBM, VIBMRepA4)
     V1, V2, V3, V4, V5 = V
     for Vi in (V1, V2, V3, V4, V5)
-        @test allequal(leftunit, sectors(Vi))
-        @test allequal(rightunit, sectors(Vi))
+        @test TensorKit._allequal(leftunit, sectors(Vi))
+        @test TensorKit._allequal(rightunit, sectors(Vi))
     end
     @test dim(V1 ⊗ V2 ⊗ V3 ← (V4 ⊗ V5)') > 0
     @test V1 ⊗ V2 ≾ (V3 ⊗ V4 ⊗ V5)' # wide blocks for QR