From ef39e510fdad6690497df5b322cd2f5d0618ac9e Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Wed, 21 Oct 2026 15:27:09 +0200 Subject: [PATCH 1/6] New truncation strategy for degenerate clusters --- Project.toml | 5 ++++- src/TensorKit.jl | 2 +- src/factorizations/factorizations.jl | 3 ++- src/factorizations/truncation.jl | 28 ++++++++++++++++++++++++++++ 4 files changed, 35 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index a759bdaef..f0545b027 100644 --- a/Project.toml +++ b/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.11" Mooncake = "0.5.27" OhMyThreads = "0.8.0" PrecompileTools = "1" @@ -72,3 +72,6 @@ TimerOutputs = "1" TupleTools = "1.5" VectorInterface = "0.6, 0.7" julia = "1.10" + +[sources] +MatrixAlgebraKit = {url = "https://github.com/quantumkithub/matrixalgebrakit.jl/", rev = "main"} diff --git a/src/TensorKit.jl b/src/TensorKit.jl index 6bdbebfb5..2d1d806db 100644 --- a/src/TensorKit.jl +++ b/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! diff --git a/src/factorizations/factorizations.jl b/src/factorizations/factorizations.jl index 7bada491c..a1753fa6a 100644 --- a/src/factorizations/factorizations.jl +++ b/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") diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index bd6956440..bf4a7ddfc 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -337,6 +337,34 @@ 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 + +_ind_to_mask!(keep::SectorVector{Bool}, inds::SectorVector{Bool}) = copyto!(keep, inds) +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) From 093101c44aa22dc3b842c018fc3a2ee299933f16 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Tue, 27 Oct 2026 06:57:57 -0400 Subject: [PATCH 2/6] Fixes and tests --- ext/TensorKitGPUArraysExt.jl | 2 +- test/amd/factorizations.jl | 25 +++++++++++++++++++++++++ test/cuda/factorizations.jl | 25 +++++++++++++++++++++++++ test/factorizations/svd.jl | 25 +++++++++++++++++++++++++ 4 files changed, 76 insertions(+), 1 deletion(-) diff --git a/ext/TensorKitGPUArraysExt.jl b/ext/TensorKitGPUArraysExt.jl index 5e013865d..45f6df640 100644 --- a/ext/TensorKitGPUArraysExt.jl +++ b/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 diff --git a/test/amd/factorizations.jl b/test/amd/factorizations.jl index b65823964..d45d165d1 100644 --- a/test/amd/factorizations.jl +++ b/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 diff --git a/test/cuda/factorizations.jl b/test/cuda/factorizations.jl index 3c482a04c..d560de296 100644 --- a/test/cuda/factorizations.jl +++ b/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 diff --git a/test/factorizations/svd.jl b/test/factorizations/svd.jl index e2eb19076..e6e8e9cf2 100644 --- a/test/factorizations/svd.jl +++ b/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 From 14d5db0c9ec20c0b5ff5d357a34cdb553c30e480 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Tue, 27 Oct 2026 08:41:03 -0400 Subject: [PATCH 3/6] Fix --- Project.toml | 5 +---- src/factorizations/truncation.jl | 5 ++++- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index f0545b027..e8d17e10f 100644 --- a/Project.toml +++ b/Project.toml @@ -57,7 +57,7 @@ FiniteDifferences = "0.12" GPUArrays = "11.4.1" LRUCache = "1.6" LinearAlgebra = "1" -MatrixAlgebraKit = "0.6.11" +MatrixAlgebraKit = "0.6.13" Mooncake = "0.5.27" OhMyThreads = "0.8.0" PrecompileTools = "1" @@ -72,6 +72,3 @@ TimerOutputs = "1" TupleTools = "1.5" VectorInterface = "0.6, 0.7" julia = "1.10" - -[sources] -MatrixAlgebraKit = {url = "https://github.com/quantumkithub/matrixalgebrakit.jl/", rev = "main"} diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index bf4a7ddfc..426e14927 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -356,7 +356,10 @@ function _findtruncated_multiplets(values::SectorVector, inds, strategy::Truncat return keep end -_ind_to_mask!(keep::SectorVector{Bool}, inds::SectorVector{Bool}) = copyto!(keep, inds) +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) From d7ec9962781a9d37cdf98f5b540eedecece316cc Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Fri, 30 Oct 2026 09:05:01 +0100 Subject: [PATCH 4/6] Drop some collects --- src/factorizations/truncation.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index 426e14927..4f6ab618e 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -348,8 +348,8 @@ 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)) + to_keep = parent(_ind_to_mask!(keep, inds)) + kept_values = parent(values) perm = sortperm(kept_values; by = abs, rev = true) MAK._keep_multiplets!(to_keep, kept_values, perm, strategy) copyto!(parent(keep), to_keep) From 53aecb00a822a1d7f933e005167b85b4b9710db5 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Fri, 30 Oct 2026 05:21:48 -0400 Subject: [PATCH 5/6] Restore collects --- src/factorizations/truncation.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index 4f6ab618e..426e14927 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -348,8 +348,8 @@ end function _findtruncated_multiplets(values::SectorVector, inds, strategy::TruncationMultiplets) keep = similar(values, Bool) - to_keep = parent(_ind_to_mask!(keep, inds)) - kept_values = parent(values) + 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) From 8c34020c5db789ffcca6f7ee87b032043dacd6e6 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Fri, 30 Oct 2026 07:39:57 -0400 Subject: [PATCH 6/6] Fix for 1.10 --- test/setup.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/setup.jl b/test/setup.jl index ef8f1896f..dca3a6c4e 100644 --- a/test/setup.jl +++ b/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