From d245da12451d9366d43a20d9936acbc57bbe90d9 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Wed, 21 Oct 2026 15:14:02 +0200 Subject: [PATCH 1/6] Add a new truncation strategy to protect degenerate clusters --- docs/src/user_interface/truncations.md | 13 +++++ src/MatrixAlgebraKit.jl | 5 +- src/implementations/truncation.jl | 34 +++++++++++ src/interface/truncation.jl | 39 +++++++++++++ test/common/truncate.jl | 79 +++++++++++++++++++++++++- 5 files changed, 167 insertions(+), 3 deletions(-) diff --git a/docs/src/user_interface/truncations.md b/docs/src/user_interface/truncations.md index 9c4ed0c8a..7a1cc8a29 100644 --- a/docs/src/user_interface/truncations.md +++ b/docs/src/user_interface/truncations.md @@ -145,6 +145,19 @@ combined_trunc = truncrank(10) & trunctol(; atol = 1e-6); combined_trunc = trunctol(; atol = 1e-6) | truncrank(3); ``` +Truncating with a fixed number of values can split a group of (near-)degenerate values, such as +when these degeneracies are due to a symmetry. Wrapping a strategy with [`truncdegenerate`](@ref) +ensures that such groups are either kept or discarded as a whole cluster: + +```@docs; canonical=false +truncdegenerate +``` + +```julia +# Keep at most 10 values, without splitting degenerate groups +degenerate_trunc = truncdegenerate(truncrank(10); rtol = 1e-10); +``` + ## Truncation Error When using truncated decompositions such as [`svd_trunc`](@ref), [`eig_trunc`](@ref), or [`eigh_trunc`](@ref), an additional truncation error value is returned. diff --git a/src/MatrixAlgebraKit.jl b/src/MatrixAlgebraKit.jl index 78d36bfe5..af11d6aff 100644 --- a/src/MatrixAlgebraKit.jl +++ b/src/MatrixAlgebraKit.jl @@ -52,7 +52,7 @@ export CUSOLVER_Simple, CUSOLVER_HouseholderQR, CUSOLVER_QRIteration, CUSOLVER_S export ROCSOLVER_HouseholderQR, ROCSOLVER_QRIteration, ROCSOLVER_Jacobi, ROCSOLVER_DivideAndConquer, ROCSOLVER_Bisection -export notrunc, truncrank, trunctol, truncerror, truncfilter +export notrunc, truncrank, trunctol, truncerror, truncfilter, truncdegenerate @static if VERSION >= v"1.11.0-DEV.469" eval( @@ -64,7 +64,8 @@ export notrunc, truncrank, trunctol, truncerror, truncfilter eval( Expr( :public, :TruncationByOrder, :TruncationByFilter, :TruncationByValue, - :TruncationByError, :TruncationIntersection, :TruncationUnion, :truncate + :TruncationByError, :TruncationIntersection, :TruncationUnion, :TruncationDegenerate, + :truncate ) ) eval( diff --git a/src/implementations/truncation.jl b/src/implementations/truncation.jl index 18f564130..15e54229a 100644 --- a/src/implementations/truncation.jl +++ b/src/implementations/truncation.jl @@ -163,6 +163,40 @@ _ind_union(A::Base.OneTo, B::Base.OneTo) = Base.OneTo(max(length(A), length(B))) _ind_union(A::AbstractUnitRange, B::AbstractUnitRange) = union(A, B) _ind_union(A, B) = union(A, B) +function findtruncated(values::AbstractVector, strategy::TruncationDegenerate) + keep = _ind_to_mask(values, findtruncated(values, strategy.strategy)) + perm = sortperm(values; by = abs, rev = true) + return findall(_keep_degenerate!(keep, values, perm, strategy)) +end +function findtruncated_svd(values::AbstractVector, strategy::TruncationDegenerate) + # values are already sorted by decreasing magnitude + keep = _ind_to_mask(values, findtruncated_svd(values, strategy.strategy)) + return findall(_keep_degenerate!(keep, values, eachindex(values), strategy)) +end + +_ind_to_mask(values::AbstractVector, ::Colon) = trues(length(values)) +_ind_to_mask(values::AbstractVector, ind::AbstractVector{Bool}) = BitVector(ind) +function _ind_to_mask(values::AbstractVector, ind) + keep = falses(length(values)) + keep[ind] .= true + return keep +end + +# walk the values in sorted order `perm`, and make sure each group of consecutive degenerate +# values is either kept or discarded as a whole cluster +function _keep_degenerate!(keep::AbstractVector{Bool}, values::AbstractVector, perm, strategy::TruncationDegenerate) + σs = abs.(collect(view(values, perm))) + start = firstindex(σs) + for k in eachindex(σs) + if k == lastindex(σs) || σs[k] - σs[k + 1] > max(strategy.atol, strategy.rtol * σs[k]) + group = view(perm, start:k) + allequal(view(keep, group)) || (keep[group] .= strategy.mode === :up) + start = k + 1 + end + end + return keep +end + # Truncation error # ---------------- truncation_error(values::AbstractVector, ind) = truncation_error!(copy(values), ind) diff --git a/src/interface/truncation.jl b/src/interface/truncation.jl index 4d4e9d10c..68c359b38 100644 --- a/src/interface/truncation.jl +++ b/src/interface/truncation.jl @@ -270,6 +270,45 @@ Base.:|(::NoTruncation, ::NoTruncation) = notrunc() Base.:|(::NoTruncation, ::TruncationUnion) = notrunc() Base.:|(::TruncationUnion, ::NoTruncation) = notrunc() +""" + TruncationDegenerate(strategy::TruncationStrategy, atol::Real, rtol::Real, mode::Symbol) + +Truncation strategy that applies `strategy` and then adjusts the selection such that groups of +(near-)degenerate values are never split. + +See also [`truncdegenerate`](@ref). +""" +struct TruncationDegenerate{S <: TruncationStrategy, T <: Real} <: TruncationStrategy + strategy::S + atol::T + rtol::T + mode::Symbol +end +function TruncationDegenerate(strategy::TruncationStrategy, atol::Real, rtol::Real, mode::Symbol) + return TruncationDegenerate(strategy, promote(atol, rtol)..., mode) +end + +""" + truncdegenerate(strategy::TruncationStrategy; atol::Real=0, rtol::Real=1e-10, mode::Symbol=:down) + +Truncation strategy that applies `strategy`, and then makes sure that groups of (near-)degenerate +values are either kept or discarded as a whole. + +Values are sorted by decreasing magnitude, and two consecutive values `abs(vᵢ) ≥ abs(vᵢ₊₁)` are +considered degenerate whenever `abs(vᵢ) - abs(vᵢ₊₁) ≤ max(atol, rtol * abs(vᵢ))`. Groups of +degenerate values that are only partially kept by `strategy` are discarded entirely if +`mode = :down`, or kept entirely if `mode = :up`. + +This is useful when the degeneracies are due to a symmetry, where splitting a +degenerate group would break the symmetry due to rounding errors. +In particular, `truncdegenerate(truncrank(howmany))` keeps at most `howmany` values for +`mode = :down`, and at least `min(howmany, length(values))` values for `mode = :up`. +""" +function truncdegenerate(strategy::TruncationStrategy; atol::Real = 0, rtol::Real = 1.0e-10, mode::Symbol = :down) + mode in (:down, :up) || throw(ArgumentError(lazy"mode must be `:down` or `:up`, got `$mode`")) + return TruncationDegenerate(strategy, atol, rtol, mode) +end + @doc """ truncation_error(values, ind) Compute the truncation error as the 2-norm of the values that are not kept by `ind`. diff --git a/test/common/truncate.jl b/test/common/truncate.jl index 8be5e7000..bc28669fe 100644 --- a/test/common/truncate.jl +++ b/test/common/truncate.jl @@ -2,7 +2,9 @@ using MatrixAlgebraKit using Test using TestExtras using MatrixAlgebraKit: NoTruncation, TruncationIntersection, TruncationUnion, - TruncationByOrder, TruncationByValue, TruncationStrategy, findtruncated, findtruncated_svd + TruncationByOrder, TruncationByValue, TruncationDegenerate, TruncationStrategy, + findtruncated, findtruncated_svd +using LinearAlgebra: Diagonal, norm @testset "truncate" begin trunc = @constinferred TruncationStrategy() @@ -94,3 +96,78 @@ using MatrixAlgebraKit: NoTruncation, TruncationIntersection, TruncationUnion, @test union2 isa TruncationUnion @test length(union2.components) == 3 end + +@testset "truncdegenerate" begin + trunc = @constinferred truncdegenerate(truncrank(2)) + @test trunc isa TruncationDegenerate + @test trunc.strategy == truncrank(2) + @test trunc.atol == 0 && trunc.rtol == 1.0e-10 + @test trunc.mode === :down + @test trunc == truncdegenerate(truncrank(2)) + @test trunc != truncdegenerate(truncrank(2); mode = :up) + @test_throws ArgumentError truncdegenerate(truncrank(2); mode = :sideways) + + # a degenerate group of three values around 0.5, split by truncrank(2) and truncrank(3) + ϵ = 1.0e-14 + vals_sorted = [1.0, 0.5 + ϵ, 0.5, -0.5 + ϵ, 0.1] + vals_unsorted = vals_sorted[[3, 5, 1, 4, 2]] + for howmany in (2, 3) + down = truncdegenerate(truncrank(howmany)) + up = truncdegenerate(truncrank(howmany); mode = :up) + @test @constinferred(findtruncated_svd(vals_sorted, down)) == [1] + @test @constinferred(findtruncated_svd(vals_sorted, up)) == 1:4 + @test sort(vals_unsorted[@constinferred(findtruncated(vals_unsorted, down))]) == [1.0] + @test sort(vals_unsorted[@constinferred(findtruncated(vals_unsorted, up))]) == + sort(vals_sorted[1:4]) + end + + # nothing changes if the cut does not split a degenerate cluster + for howmany in (0, 1, 4, 5, 6) + for mode in (:down, :up) + strategy = truncdegenerate(truncrank(howmany); mode) + @test findtruncated_svd(vals_sorted, strategy) == + findtruncated_svd(vals_sorted, truncrank(howmany)) + @test sort(findtruncated(vals_unsorted, strategy)) == + sort(findtruncated(vals_unsorted, truncrank(howmany))) + end + end + @test findtruncated_svd(vals_sorted, truncdegenerate(notrunc())) == 1:5 + @test isempty(findtruncated_svd(Float64[], truncdegenerate(truncrank(2)))) + + # without tolerance, only exactly degenerate values are grouped + @test findtruncated_svd(vals_sorted, truncdegenerate(truncrank(2); rtol = 0)) == 1:2 + @test findtruncated_svd([1.0, 0.5, 0.5, 0.1], truncdegenerate(truncrank(2); rtol = 0)) == [1] + + vals_chain = [1.0, 0.9, 0.8, 0.7, 0.1] + @test findtruncated_svd(vals_chain, truncdegenerate(truncrank(2); atol = 0.15)) == Int[] + @test findtruncated_svd(vals_chain, truncdegenerate(truncrank(2); atol = 0.15, mode = :up)) == 1:4 + @test findtruncated_svd(vals_chain, truncdegenerate(truncrank(2); rtol = 0.12)) == Int[] + @test findtruncated_svd(vals_chain, truncdegenerate(truncrank(2); rtol = 0.12, mode = :up)) == 1:3 + + strategy = truncdegenerate(truncrank(3) & trunctol(; atol = 0.2)) + @test findtruncated_svd(vals_sorted, strategy) == [1] + strategy = truncdegenerate(trunctol(; atol = 0.5) | truncrank(1); mode = :up) + @test findtruncated_svd(vals_sorted, strategy) == 1:4 + strategy = truncdegenerate(truncfilter(x -> x < 0.5 + ϵ / 2); mode = :up) + @test sort(findtruncated(vals_unsorted, strategy)) == sort(findtruncated(vals_unsorted, truncfilter(x -> x < 0.9))) + @test truncdegenerate(truncrank(2)) & trunctol(; atol = 0.2) isa TruncationIntersection + + # exactly degenerate singular values + for T in (Float64, ComplexF64) + n = 6 + Q₁ = qr_compact(randn(T, n, n))[1] + Q₂ = qr_compact(randn(T, n, n))[1] + Sref = [3.0, 2.0, 2.0, 2.0, 1.0, 0.5] + A = Q₁ * Diagonal(Sref) * Q₂' + for howmany in 2:3 + U, S, Vᴴ, ϵ = svd_trunc(A; trunc = truncdegenerate(truncrank(howmany); rtol = 1.0e-8)) + @test diagview(S) ≈ Sref[1:1] + @test ϵ ≈ norm(Sref[2:end]) + @test U * S * Vᴴ ≈ Sref[1] * Q₁[:, 1] * Q₂[:, 1]' + U, S, Vᴴ, ϵ = svd_trunc(A; trunc = truncdegenerate(truncrank(howmany); rtol = 1.0e-8, mode = :up)) + @test diagview(S) ≈ Sref[1:4] + @test ϵ ≈ norm(Sref[5:end]) + @test U * S * Vᴴ ≈ Q₁[:, 1:4] * Diagonal(Sref[1:4]) * Q₂[:, 1:4]' + end + end +end From 59ec1619d88ac3acf5b58b06bbad964b0fff6628 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Thu, 22 Oct 2026 02:37:44 -0400 Subject: [PATCH 2/6] Formatter --- test/common/truncate.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/common/truncate.jl b/test/common/truncate.jl index bc28669fe..11f04cf46 100644 --- a/test/common/truncate.jl +++ b/test/common/truncate.jl @@ -121,7 +121,7 @@ end sort(vals_sorted[1:4]) end - # nothing changes if the cut does not split a degenerate cluster + # nothing changes if the cut does not split a degenerate cluster for howmany in (0, 1, 4, 5, 6) for mode in (:down, :up) strategy = truncdegenerate(truncrank(howmany); mode) From fd198afdb507fad496ed1a2598a144f5b90b04ac Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Thu, 22 Oct 2026 09:05:11 -0400 Subject: [PATCH 3/6] Test truncdegenerate for SVD also --- test/testsuite/decompositions/svd.jl | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/test/testsuite/decompositions/svd.jl b/test/testsuite/decompositions/svd.jl index a5fd80c9b..6723338b6 100644 --- a/test/testsuite/decompositions/svd.jl +++ b/test/testsuite/decompositions/svd.jl @@ -686,6 +686,21 @@ function test_svd_trunc_algs( @test ϵ2 ≈ norm(diagview(S)[3:4]) atol = atol @test_throws ArgumentError svd_trunc(A; alg = truncalg, trunc = (; maxrank = 2)) @test_throws ArgumentError svd_trunc_no_error(A; alg = truncalg, trunc = (; maxrank = 2)) + + # truncdegenerate + Sdiag = similar(A, real(eltype(T)), m4) + copyto!(Sdiag, real(eltype(T))[0.9, 0.3, 0.3, 0.01]) + S = Diagonal(Sdiag) + Vᴴ = instantiate_unitary(T, A, m4) + A = U * S * Vᴴ + truncalg = TruncatedAlgorithm(alg, truncdegenerate(truncrank(2); atol = 0.2)) + U2, S2, V2ᴴ, ϵ2 = @testinferred svd_trunc(A; alg = truncalg) + @test collect(diagview(S2)) ≈ collect(diagview(S)[1:1]) + @test ϵ2 ≈ norm(diagview(S)[2:4]) atol = atol + truncalg = TruncatedAlgorithm(alg, truncdegenerate(truncrank(3); atol = 0.2)) + U2, S2, V2ᴴ, ϵ2 = @testinferred svd_trunc(A; alg = truncalg) + @test collect(diagview(S2)) ≈ collect(diagview(S)[1:3]) + @test ϵ2 ≈ norm(diagview(S)[4:4]) atol = atol end end end From 33ec0df39815ee1444f6edb6ca3cc93a95f0852e Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Thu, 22 Oct 2026 15:36:17 +0200 Subject: [PATCH 4/6] Formatter --- test/testsuite/decompositions/svd.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/testsuite/decompositions/svd.jl b/test/testsuite/decompositions/svd.jl index 6723338b6..d5a3da1b0 100644 --- a/test/testsuite/decompositions/svd.jl +++ b/test/testsuite/decompositions/svd.jl @@ -686,7 +686,7 @@ function test_svd_trunc_algs( @test ϵ2 ≈ norm(diagview(S)[3:4]) atol = atol @test_throws ArgumentError svd_trunc(A; alg = truncalg, trunc = (; maxrank = 2)) @test_throws ArgumentError svd_trunc_no_error(A; alg = truncalg, trunc = (; maxrank = 2)) - + # truncdegenerate Sdiag = similar(A, real(eltype(T)), m4) copyto!(Sdiag, real(eltype(T))[0.9, 0.3, 0.3, 0.01]) From 873be726603ce18bfedafbad90bee68ff44f9034 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Fri, 23 Oct 2026 04:21:21 -0400 Subject: [PATCH 5/6] Rename the strategy to multiplets --- docs/src/user_interface/truncations.md | 6 ++-- src/MatrixAlgebraKit.jl | 4 +-- src/implementations/truncation.jl | 10 +++--- src/interface/truncation.jl | 18 +++++----- test/common/truncate.jl | 48 +++++++++++++------------- test/testsuite/decompositions/svd.jl | 6 ++-- 6 files changed, 46 insertions(+), 46 deletions(-) diff --git a/docs/src/user_interface/truncations.md b/docs/src/user_interface/truncations.md index 7a1cc8a29..b829b5936 100644 --- a/docs/src/user_interface/truncations.md +++ b/docs/src/user_interface/truncations.md @@ -146,16 +146,16 @@ combined_trunc = trunctol(; atol = 1e-6) | truncrank(3); ``` Truncating with a fixed number of values can split a group of (near-)degenerate values, such as -when these degeneracies are due to a symmetry. Wrapping a strategy with [`truncdegenerate`](@ref) +when these degeneracies are due to a symmetry. Wrapping a strategy with [`truncmultiplets`](@ref) ensures that such groups are either kept or discarded as a whole cluster: ```@docs; canonical=false -truncdegenerate +truncmultiplets ``` ```julia # Keep at most 10 values, without splitting degenerate groups -degenerate_trunc = truncdegenerate(truncrank(10); rtol = 1e-10); +degenerate_trunc = truncmultiplets(truncrank(10); rtol = 1e-10); ``` ## Truncation Error diff --git a/src/MatrixAlgebraKit.jl b/src/MatrixAlgebraKit.jl index af11d6aff..3a9f4b413 100644 --- a/src/MatrixAlgebraKit.jl +++ b/src/MatrixAlgebraKit.jl @@ -52,7 +52,7 @@ export CUSOLVER_Simple, CUSOLVER_HouseholderQR, CUSOLVER_QRIteration, CUSOLVER_S export ROCSOLVER_HouseholderQR, ROCSOLVER_QRIteration, ROCSOLVER_Jacobi, ROCSOLVER_DivideAndConquer, ROCSOLVER_Bisection -export notrunc, truncrank, trunctol, truncerror, truncfilter, truncdegenerate +export notrunc, truncrank, trunctol, truncerror, truncfilter, truncmultiplets @static if VERSION >= v"1.11.0-DEV.469" eval( @@ -64,7 +64,7 @@ export notrunc, truncrank, trunctol, truncerror, truncfilter, truncdegenerate eval( Expr( :public, :TruncationByOrder, :TruncationByFilter, :TruncationByValue, - :TruncationByError, :TruncationIntersection, :TruncationUnion, :TruncationDegenerate, + :TruncationByError, :TruncationIntersection, :TruncationUnion, :TruncationMultiplets, :truncate ) ) diff --git a/src/implementations/truncation.jl b/src/implementations/truncation.jl index 15e54229a..55d87be06 100644 --- a/src/implementations/truncation.jl +++ b/src/implementations/truncation.jl @@ -163,15 +163,15 @@ _ind_union(A::Base.OneTo, B::Base.OneTo) = Base.OneTo(max(length(A), length(B))) _ind_union(A::AbstractUnitRange, B::AbstractUnitRange) = union(A, B) _ind_union(A, B) = union(A, B) -function findtruncated(values::AbstractVector, strategy::TruncationDegenerate) +function findtruncated(values::AbstractVector, strategy::TruncationMultiplets) keep = _ind_to_mask(values, findtruncated(values, strategy.strategy)) perm = sortperm(values; by = abs, rev = true) - return findall(_keep_degenerate!(keep, values, perm, strategy)) + return findall(_keep_multiplets!(keep, values, perm, strategy)) end -function findtruncated_svd(values::AbstractVector, strategy::TruncationDegenerate) +function findtruncated_svd(values::AbstractVector, strategy::TruncationMultiplets) # values are already sorted by decreasing magnitude keep = _ind_to_mask(values, findtruncated_svd(values, strategy.strategy)) - return findall(_keep_degenerate!(keep, values, eachindex(values), strategy)) + return findall(_keep_multiplets!(keep, values, eachindex(values), strategy)) end _ind_to_mask(values::AbstractVector, ::Colon) = trues(length(values)) @@ -184,7 +184,7 @@ end # walk the values in sorted order `perm`, and make sure each group of consecutive degenerate # values is either kept or discarded as a whole cluster -function _keep_degenerate!(keep::AbstractVector{Bool}, values::AbstractVector, perm, strategy::TruncationDegenerate) +function _keep_multiplets!(keep::AbstractVector{Bool}, values::AbstractVector, perm, strategy::TruncationMultiplets) σs = abs.(collect(view(values, perm))) start = firstindex(σs) for k in eachindex(σs) diff --git a/src/interface/truncation.jl b/src/interface/truncation.jl index 68c359b38..508344472 100644 --- a/src/interface/truncation.jl +++ b/src/interface/truncation.jl @@ -271,25 +271,25 @@ Base.:|(::NoTruncation, ::TruncationUnion) = notrunc() Base.:|(::TruncationUnion, ::NoTruncation) = notrunc() """ - TruncationDegenerate(strategy::TruncationStrategy, atol::Real, rtol::Real, mode::Symbol) + TruncationMultiplets(strategy::TruncationStrategy, atol::Real, rtol::Real, mode::Symbol) Truncation strategy that applies `strategy` and then adjusts the selection such that groups of (near-)degenerate values are never split. -See also [`truncdegenerate`](@ref). +See also [`truncmultiplets`](@ref). """ -struct TruncationDegenerate{S <: TruncationStrategy, T <: Real} <: TruncationStrategy +struct TruncationMultiplets{S <: TruncationStrategy, T <: Real} <: TruncationStrategy strategy::S atol::T rtol::T mode::Symbol end -function TruncationDegenerate(strategy::TruncationStrategy, atol::Real, rtol::Real, mode::Symbol) - return TruncationDegenerate(strategy, promote(atol, rtol)..., mode) +function TruncationMultiplets(strategy::TruncationStrategy, atol::Real, rtol::Real, mode::Symbol) + return TruncationMultiplets(strategy, promote(atol, rtol)..., mode) end """ - truncdegenerate(strategy::TruncationStrategy; atol::Real=0, rtol::Real=1e-10, mode::Symbol=:down) + truncmultiplets(strategy::TruncationStrategy; atol::Real=0, rtol::Real=1e-10, mode::Symbol=:down) Truncation strategy that applies `strategy`, and then makes sure that groups of (near-)degenerate values are either kept or discarded as a whole. @@ -301,12 +301,12 @@ degenerate values that are only partially kept by `strategy` are discarded entir This is useful when the degeneracies are due to a symmetry, where splitting a degenerate group would break the symmetry due to rounding errors. -In particular, `truncdegenerate(truncrank(howmany))` keeps at most `howmany` values for +In particular, `truncmultiplets(truncrank(howmany))` keeps at most `howmany` values for `mode = :down`, and at least `min(howmany, length(values))` values for `mode = :up`. """ -function truncdegenerate(strategy::TruncationStrategy; atol::Real = 0, rtol::Real = 1.0e-10, mode::Symbol = :down) +function truncmultiplets(strategy::TruncationStrategy; atol::Real = 0, rtol::Real = 1.0e-10, mode::Symbol = :down) mode in (:down, :up) || throw(ArgumentError(lazy"mode must be `:down` or `:up`, got `$mode`")) - return TruncationDegenerate(strategy, atol, rtol, mode) + return TruncationMultiplets(strategy, atol, rtol, mode) end @doc """ diff --git a/test/common/truncate.jl b/test/common/truncate.jl index 11f04cf46..2c68d2bb0 100644 --- a/test/common/truncate.jl +++ b/test/common/truncate.jl @@ -2,7 +2,7 @@ using MatrixAlgebraKit using Test using TestExtras using MatrixAlgebraKit: NoTruncation, TruncationIntersection, TruncationUnion, - TruncationByOrder, TruncationByValue, TruncationDegenerate, TruncationStrategy, + TruncationByOrder, TruncationByValue, TruncationMultiplets, TruncationStrategy, findtruncated, findtruncated_svd using LinearAlgebra: Diagonal, norm @@ -97,23 +97,23 @@ using LinearAlgebra: Diagonal, norm @test length(union2.components) == 3 end -@testset "truncdegenerate" begin - trunc = @constinferred truncdegenerate(truncrank(2)) - @test trunc isa TruncationDegenerate +@testset "truncmultiplets" begin + trunc = @constinferred truncmultiplets(truncrank(2)) + @test trunc isa TruncationMultiplets @test trunc.strategy == truncrank(2) @test trunc.atol == 0 && trunc.rtol == 1.0e-10 @test trunc.mode === :down - @test trunc == truncdegenerate(truncrank(2)) - @test trunc != truncdegenerate(truncrank(2); mode = :up) - @test_throws ArgumentError truncdegenerate(truncrank(2); mode = :sideways) + @test trunc == truncmultiplets(truncrank(2)) + @test trunc != truncmultiplets(truncrank(2); mode = :up) + @test_throws ArgumentError truncmultiplets(truncrank(2); mode = :sideways) # a degenerate group of three values around 0.5, split by truncrank(2) and truncrank(3) ϵ = 1.0e-14 vals_sorted = [1.0, 0.5 + ϵ, 0.5, -0.5 + ϵ, 0.1] vals_unsorted = vals_sorted[[3, 5, 1, 4, 2]] for howmany in (2, 3) - down = truncdegenerate(truncrank(howmany)) - up = truncdegenerate(truncrank(howmany); mode = :up) + down = truncmultiplets(truncrank(howmany)) + up = truncmultiplets(truncrank(howmany); mode = :up) @test @constinferred(findtruncated_svd(vals_sorted, down)) == [1] @test @constinferred(findtruncated_svd(vals_sorted, up)) == 1:4 @test sort(vals_unsorted[@constinferred(findtruncated(vals_unsorted, down))]) == [1.0] @@ -124,33 +124,33 @@ end # nothing changes if the cut does not split a degenerate cluster for howmany in (0, 1, 4, 5, 6) for mode in (:down, :up) - strategy = truncdegenerate(truncrank(howmany); mode) + strategy = truncmultiplets(truncrank(howmany); mode) @test findtruncated_svd(vals_sorted, strategy) == findtruncated_svd(vals_sorted, truncrank(howmany)) @test sort(findtruncated(vals_unsorted, strategy)) == sort(findtruncated(vals_unsorted, truncrank(howmany))) end end - @test findtruncated_svd(vals_sorted, truncdegenerate(notrunc())) == 1:5 - @test isempty(findtruncated_svd(Float64[], truncdegenerate(truncrank(2)))) + @test findtruncated_svd(vals_sorted, truncmultiplets(notrunc())) == 1:5 + @test isempty(findtruncated_svd(Float64[], truncmultiplets(truncrank(2)))) # without tolerance, only exactly degenerate values are grouped - @test findtruncated_svd(vals_sorted, truncdegenerate(truncrank(2); rtol = 0)) == 1:2 - @test findtruncated_svd([1.0, 0.5, 0.5, 0.1], truncdegenerate(truncrank(2); rtol = 0)) == [1] + @test findtruncated_svd(vals_sorted, truncmultiplets(truncrank(2); rtol = 0)) == 1:2 + @test findtruncated_svd([1.0, 0.5, 0.5, 0.1], truncmultiplets(truncrank(2); rtol = 0)) == [1] vals_chain = [1.0, 0.9, 0.8, 0.7, 0.1] - @test findtruncated_svd(vals_chain, truncdegenerate(truncrank(2); atol = 0.15)) == Int[] - @test findtruncated_svd(vals_chain, truncdegenerate(truncrank(2); atol = 0.15, mode = :up)) == 1:4 - @test findtruncated_svd(vals_chain, truncdegenerate(truncrank(2); rtol = 0.12)) == Int[] - @test findtruncated_svd(vals_chain, truncdegenerate(truncrank(2); rtol = 0.12, mode = :up)) == 1:3 + @test findtruncated_svd(vals_chain, truncmultiplets(truncrank(2); atol = 0.15)) == Int[] + @test findtruncated_svd(vals_chain, truncmultiplets(truncrank(2); atol = 0.15, mode = :up)) == 1:4 + @test findtruncated_svd(vals_chain, truncmultiplets(truncrank(2); rtol = 0.12)) == Int[] + @test findtruncated_svd(vals_chain, truncmultiplets(truncrank(2); rtol = 0.12, mode = :up)) == 1:3 - strategy = truncdegenerate(truncrank(3) & trunctol(; atol = 0.2)) + strategy = truncmultiplets(truncrank(3) & trunctol(; atol = 0.2)) @test findtruncated_svd(vals_sorted, strategy) == [1] - strategy = truncdegenerate(trunctol(; atol = 0.5) | truncrank(1); mode = :up) + strategy = truncmultiplets(trunctol(; atol = 0.5) | truncrank(1); mode = :up) @test findtruncated_svd(vals_sorted, strategy) == 1:4 - strategy = truncdegenerate(truncfilter(x -> x < 0.5 + ϵ / 2); mode = :up) + strategy = truncmultiplets(truncfilter(x -> x < 0.5 + ϵ / 2); mode = :up) @test sort(findtruncated(vals_unsorted, strategy)) == sort(findtruncated(vals_unsorted, truncfilter(x -> x < 0.9))) - @test truncdegenerate(truncrank(2)) & trunctol(; atol = 0.2) isa TruncationIntersection + @test truncmultiplets(truncrank(2)) & trunctol(; atol = 0.2) isa TruncationIntersection # exactly degenerate singular values for T in (Float64, ComplexF64) @@ -160,11 +160,11 @@ end Sref = [3.0, 2.0, 2.0, 2.0, 1.0, 0.5] A = Q₁ * Diagonal(Sref) * Q₂' for howmany in 2:3 - U, S, Vᴴ, ϵ = svd_trunc(A; trunc = truncdegenerate(truncrank(howmany); rtol = 1.0e-8)) + U, S, Vᴴ, ϵ = svd_trunc(A; trunc = truncmultiplets(truncrank(howmany); rtol = 1.0e-8)) @test diagview(S) ≈ Sref[1:1] @test ϵ ≈ norm(Sref[2:end]) @test U * S * Vᴴ ≈ Sref[1] * Q₁[:, 1] * Q₂[:, 1]' - U, S, Vᴴ, ϵ = svd_trunc(A; trunc = truncdegenerate(truncrank(howmany); rtol = 1.0e-8, mode = :up)) + U, S, Vᴴ, ϵ = svd_trunc(A; trunc = truncmultiplets(truncrank(howmany); rtol = 1.0e-8, mode = :up)) @test diagview(S) ≈ Sref[1:4] @test ϵ ≈ norm(Sref[5:end]) @test U * S * Vᴴ ≈ Q₁[:, 1:4] * Diagonal(Sref[1:4]) * Q₂[:, 1:4]' diff --git a/test/testsuite/decompositions/svd.jl b/test/testsuite/decompositions/svd.jl index d5a3da1b0..1baeebbec 100644 --- a/test/testsuite/decompositions/svd.jl +++ b/test/testsuite/decompositions/svd.jl @@ -687,17 +687,17 @@ function test_svd_trunc_algs( @test_throws ArgumentError svd_trunc(A; alg = truncalg, trunc = (; maxrank = 2)) @test_throws ArgumentError svd_trunc_no_error(A; alg = truncalg, trunc = (; maxrank = 2)) - # truncdegenerate + # truncmultiplets Sdiag = similar(A, real(eltype(T)), m4) copyto!(Sdiag, real(eltype(T))[0.9, 0.3, 0.3, 0.01]) S = Diagonal(Sdiag) Vᴴ = instantiate_unitary(T, A, m4) A = U * S * Vᴴ - truncalg = TruncatedAlgorithm(alg, truncdegenerate(truncrank(2); atol = 0.2)) + truncalg = TruncatedAlgorithm(alg, truncmultiplets(truncrank(2); atol = 0.2)) U2, S2, V2ᴴ, ϵ2 = @testinferred svd_trunc(A; alg = truncalg) @test collect(diagview(S2)) ≈ collect(diagview(S)[1:1]) @test ϵ2 ≈ norm(diagview(S)[2:4]) atol = atol - truncalg = TruncatedAlgorithm(alg, truncdegenerate(truncrank(3); atol = 0.2)) + truncalg = TruncatedAlgorithm(alg, truncmultiplets(truncrank(3); atol = 0.2)) U2, S2, V2ᴴ, ϵ2 = @testinferred svd_trunc(A; alg = truncalg) @test collect(diagview(S2)) ≈ collect(diagview(S)[1:3]) @test ϵ2 ≈ norm(diagview(S)[4:4]) atol = atol From fe7bb6f45a5e2ec29f4e52488899f6763037eef0 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Fri, 23 Oct 2026 11:16:52 -0400 Subject: [PATCH 6/6] Comments --- src/implementations/truncation.jl | 6 +++--- src/interface/truncation.jl | 10 +++++----- test/common/truncate.jl | 2 +- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/implementations/truncation.jl b/src/implementations/truncation.jl index 55d87be06..44a4c110d 100644 --- a/src/implementations/truncation.jl +++ b/src/implementations/truncation.jl @@ -164,13 +164,13 @@ _ind_union(A::AbstractUnitRange, B::AbstractUnitRange) = union(A, B) _ind_union(A, B) = union(A, B) function findtruncated(values::AbstractVector, strategy::TruncationMultiplets) - keep = _ind_to_mask(values, findtruncated(values, strategy.strategy)) + keep = _ind_to_mask(values, findtruncated(values, strategy.trunc)) perm = sortperm(values; by = abs, rev = true) return findall(_keep_multiplets!(keep, values, perm, strategy)) end function findtruncated_svd(values::AbstractVector, strategy::TruncationMultiplets) # values are already sorted by decreasing magnitude - keep = _ind_to_mask(values, findtruncated_svd(values, strategy.strategy)) + keep = _ind_to_mask(values, findtruncated_svd(values, strategy.trunc)) return findall(_keep_multiplets!(keep, values, eachindex(values), strategy)) end @@ -190,7 +190,7 @@ function _keep_multiplets!(keep::AbstractVector{Bool}, values::AbstractVector, p for k in eachindex(σs) if k == lastindex(σs) || σs[k] - σs[k + 1] > max(strategy.atol, strategy.rtol * σs[k]) group = view(perm, start:k) - allequal(view(keep, group)) || (keep[group] .= strategy.mode === :up) + allequal(view(keep, group)) || (keep[group] .= strategy.mode == :up) start = k + 1 end end diff --git a/src/interface/truncation.jl b/src/interface/truncation.jl index 508344472..0eaa79dff 100644 --- a/src/interface/truncation.jl +++ b/src/interface/truncation.jl @@ -271,21 +271,21 @@ Base.:|(::NoTruncation, ::TruncationUnion) = notrunc() Base.:|(::TruncationUnion, ::NoTruncation) = notrunc() """ - TruncationMultiplets(strategy::TruncationStrategy, atol::Real, rtol::Real, mode::Symbol) + TruncationMultiplets(trunc::TruncationStrategy, atol::Real, rtol::Real, mode::Symbol) -Truncation strategy that applies `strategy` and then adjusts the selection such that groups of +Truncation strategy that applies `trunc` and then adjusts the selection such that groups of (near-)degenerate values are never split. See also [`truncmultiplets`](@ref). """ struct TruncationMultiplets{S <: TruncationStrategy, T <: Real} <: TruncationStrategy - strategy::S + trunc::S atol::T rtol::T mode::Symbol end -function TruncationMultiplets(strategy::TruncationStrategy, atol::Real, rtol::Real, mode::Symbol) - return TruncationMultiplets(strategy, promote(atol, rtol)..., mode) +function TruncationMultiplets(trunc::TruncationStrategy, atol::Real, rtol::Real, mode::Symbol) + return TruncationMultiplets(trunc, promote(atol, rtol)..., mode) end """ diff --git a/test/common/truncate.jl b/test/common/truncate.jl index 2c68d2bb0..0c71b18aa 100644 --- a/test/common/truncate.jl +++ b/test/common/truncate.jl @@ -100,7 +100,7 @@ end @testset "truncmultiplets" begin trunc = @constinferred truncmultiplets(truncrank(2)) @test trunc isa TruncationMultiplets - @test trunc.strategy == truncrank(2) + @test trunc.trunc == truncrank(2) @test trunc.atol == 0 && trunc.rtol == 1.0e-10 @test trunc.mode === :down @test trunc == truncmultiplets(truncrank(2))