From 021e860796f3638d7168fbc7cee2e3eb9740d6f3 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Fri, 19 Dec 2025 12:29:49 -0500 Subject: [PATCH 1/8] use `SectorVector` as output for `svd_vals(::DiagonalTensorMap)` --- src/factorizations/diagonal.jl | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/factorizations/diagonal.jl b/src/factorizations/diagonal.jl index a101b1f21..daf8b9fff 100644 --- a/src/factorizations/diagonal.jl +++ b/src/factorizations/diagonal.jl @@ -1,6 +1,7 @@ # DiagonalTensorMap # ----------------- _repack_diagonal(d::DiagonalTensorMap) = Diagonal(d.data) +_repack_diagonal(d::SectorVector) = Diagonal(parent(d)) MAK.diagview(t::DiagonalTensorMap) = SectorVector(t.data, TensorKit.diagonalblockstructure(space(t))) @@ -100,10 +101,12 @@ for f! in (:eig_vals!, :eigh_vals!, :svd_vals!) $f!(_repack_diagonal(d), diagview(_repack_diagonal(V)), alg) return V end - @eval function MAK.initialize_output( - ::typeof($f!), d::DiagonalTensorMap, alg::DiagonalAlgorithm - ) - data = MAK.initialize_output($f!, _repack_diagonal(d), alg) - return DiagonalTensorMap(data, d.domain) - end +end + +# For diagonal inputs we don't have to promote the scalartype since we know they are symmetric +function MAK.initialize_output(::typeof(eig_vals!), t::AbstractTensorMap, alg::DiagonalAlgorithm) + V_D = fuse(domain(t)) + Tc = scalartype(t) + A = similarstoragetype(t, Tc) + return SectorVector{Tc, sectortype(t), A}(undef, V_D) end From efd44420f6786148c581ba4bc71eb940c63e6f3d Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Fri, 19 Dec 2025 12:30:17 -0500 Subject: [PATCH 2/8] Don't sort diagonal svd_vals between blocks --- src/factorizations/diagonal.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/factorizations/diagonal.jl b/src/factorizations/diagonal.jl index daf8b9fff..5e4bd9e6c 100644 --- a/src/factorizations/diagonal.jl +++ b/src/factorizations/diagonal.jl @@ -96,7 +96,8 @@ end # f_vals # ------ -for f! in (:eig_vals!, :eigh_vals!, :svd_vals!) +# this is invalid for svd_vals since we cannot sort between blocks +for f! in (:eig_vals!, :eigh_vals!) @eval function MAK.$f!(d::AbstractTensorMap, V, alg::DiagonalAlgorithm) $f!(_repack_diagonal(d), diagview(_repack_diagonal(V)), alg) return V From 813e825137bfce58c9bae3157af49f25be287d79 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Fri, 19 Dec 2025 12:38:50 -0500 Subject: [PATCH 3/8] update tests --- test/tensors/factorizations.jl | 55 ++++++++++++++++------------------ 1 file changed, 26 insertions(+), 29 deletions(-) diff --git a/test/tensors/factorizations.jl b/test/tensors/factorizations.jl index b017b746b..176d62657 100644 --- a/test/tensors/factorizations.jl +++ b/test/tensors/factorizations.jl @@ -1,6 +1,7 @@ using Test, TestExtras using TensorKit using LinearAlgebra: LinearAlgebra +using MatrixAlgebraKit: diagview @isdefined(TestSetup) || include("../setup.jl") using .TestSetup @@ -191,10 +192,9 @@ for V in spacelist @test isposdef(s) @test isisometric(vᴴ; side = :right) - s′ = LinearAlgebra.diag(s) - for (c, b) in pairs(LinearAlgebra.svdvals(t)) - @test b ≈ s′[c] - end + s′ = @constinferred svd_vals(t) + @test s′ ≈ diagview(s) + @test s′ isa TensorKit.SectorVector v, c = @constinferred left_orth(t; alg = :svd) @test v * c ≈ t @@ -261,14 +261,14 @@ for V in spacelist @test norm(t - U1 * S1 * Vᴴ1) ≈ ϵ1 atol = eps(real(T))^(4 / 5) @test dim(domain(S1)) <= nvals - λ = minimum(minimum, values(LinearAlgebra.diag(S1))) + λ = minimum(diagview(S1)) trunc = trunctol(; atol = λ - 10eps(λ)) U2, S2, Vᴴ2, ϵ2 = @constinferred svd_trunc(t; trunc) @test t * Vᴴ2' ≈ U2 * S2 @test isisometric(U2) @test isisometric(Vᴴ2; side = :right) @test norm(t - U2 * S2 * Vᴴ2) ≈ ϵ2 atol = eps(real(T))^(4 / 5) - @test minimum(minimum, values(LinearAlgebra.diag(S1))) >= λ + @test minimum(diagview(S1)) >= λ @test U2 ≈ U1 @test S2 ≈ S1 @test Vᴴ2 ≈ Vᴴ1 @@ -297,7 +297,7 @@ for V in spacelist @test isisometric(U5) @test isisometric(Vᴴ5; side = :right) @test norm(t - U5 * S5 * Vᴴ5) ≈ ϵ5 atol = eps(real(T))^(4 / 5) - @test minimum(minimum, values(LinearAlgebra.diag(S5))) >= λ + @test minimum(diagview(S5)) >= λ @test dim(domain(S5)) ≤ nvals end end @@ -312,13 +312,11 @@ for V in spacelist d, v = @constinferred eig_full(t) @test t * v ≈ v * d - d′ = LinearAlgebra.diag(d) - for (c, b) in pairs(LinearAlgebra.eigvals(t)) - @test sort(b; by = abs) ≈ sort(d′[c]; by = abs) - end + d′ = @constinferred eig_vals(t) + @test d′ ≈ diagview(d) + @test d′ isa TensorKit.SectorVector - vdv = v' * v - vdv = (vdv + vdv') / 2 + vdv = project_hermitian!(v' * v) @test @constinferred isposdef(vdv) t isa DiagonalTensorMap || @test !isposdef(t) # unlikely for non-hermitian map @@ -327,35 +325,34 @@ for V in spacelist @test t * v ≈ v * d @test dim(domain(d)) ≤ nvals - t2 = (t + t') + t2 = @constinferred project_hermitian(t) D, V = eigen(t2) @test isisometric(V) D̃, Ṽ = @constinferred eigh_full(t2) @test D ≈ D̃ @test V ≈ Ṽ - λ = minimum( - minimum(real(LinearAlgebra.diag(b))) - for (c, b) in blocks(D) - ) + λ = minimum(real, diagview(D)) @test cond(Ṽ) ≈ one(real(T)) @test isposdef(t2) == isposdef(λ) @test isposdef(t2 - λ * one(t2) + 0.1 * one(t2)) @test !isposdef(t2 - λ * one(t2) - 0.1 * one(t2)) - add!(t, t') - - d, v = @constinferred eigh_full(t) - @test t * v ≈ v * d + d, v = @constinferred eigh_full(t2) + @test t2 * v ≈ v * d @test isunitary(v) - λ = minimum(minimum(real(LinearAlgebra.diag(b))) for (c, b) in blocks(d)) + d′ = @constinferred eigh_vals(t2) + @test d′ ≈ diagview(d) + @test d′ isa TensorKit.SectorVector + + λ = minimum(real, diagview(d)) @test cond(v) ≈ one(real(T)) - @test isposdef(t) == isposdef(λ) - @test isposdef(t - λ * one(t) + 0.1 * one(t)) - @test !isposdef(t - λ * one(t) - 0.1 * one(t)) + @test isposdef(t2) == isposdef(λ) + @test isposdef(t2 - λ * one(t) + 0.1 * one(t2)) + @test !isposdef(t2 - λ * one(t) - 0.1 * one(t2)) - d, v = @constinferred eigh_trunc(t; trunc = truncrank(nvals)) - @test t * v ≈ v * d + d, v = @constinferred eigh_trunc(t2; trunc = truncrank(nvals)) + @test t2 * v ≈ v * d @test dim(domain(d)) ≤ nvals end end @@ -390,7 +387,7 @@ for V in spacelist @test cond(t2) == 0.0 end for T in eltypes, t in (rand(T, W, W), rand(T, W, W)') - add!(t, t') + project_hermitian!(t) vals = @constinferred LinearAlgebra.eigvals(t) λmax = maximum(s -> maximum(abs, s), values(vals)) λmin = minimum(s -> minimum(abs, s), values(vals)) From 32d2320aa0bcb68be13e4059101ab4de07a14224 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Fri, 19 Dec 2025 12:42:33 -0500 Subject: [PATCH 4/8] update changelog --- docs/src/Changelog.md | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/src/Changelog.md b/docs/src/Changelog.md index 6309d0c03..d2ff74b51 100644 --- a/docs/src/Changelog.md +++ b/docs/src/Changelog.md @@ -31,6 +31,7 @@ When releasing a new version, move the "Unreleased" changes to a new version sec - Issue with using relative tolerances in truncation schemes ([#314](https://github.com/QuantumKitHub/TensorKit.jl/issues/314)) - Using `scalartype` instead of `eltype` in BLAS contraction ([#326](https://github.com/QuantumKitHub/TensorKit.jl/pull/326)) - Divide by zero error in `show` for empty tensors ([#329](https://github.com/QuantumKitHub/TensorKit.jl/pull/329)) +- `svd_vals(::DiagonalTensorMap)` correctly outputs `SectorVector` and implementation fix. ([#333](https://github.com/QuantumKitHub/TensorKit.jl/pull/329)) ## [0.16.0](https://github.com/QuantumKitHub/TensorKit.jl/releases/tag/v0.16.0) - 2025-12-08 From ed459475213321cfcd80e1a91bef0561529abc4b Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sun, 21 Dec 2025 15:47:10 -0500 Subject: [PATCH 5/8] be careful about not mixing eigenvalues through blocks --- src/factorizations/diagonal.jl | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/src/factorizations/diagonal.jl b/src/factorizations/diagonal.jl index 5e4bd9e6c..cec4a9061 100644 --- a/src/factorizations/diagonal.jl +++ b/src/factorizations/diagonal.jl @@ -94,16 +94,6 @@ function MAK.svd_compact!(t::AbstractTensorMap, USVᴴ, alg::DiagonalAlgorithm) return svd_full!(t, USVᴴ, alg) end -# f_vals -# ------ -# this is invalid for svd_vals since we cannot sort between blocks -for f! in (:eig_vals!, :eigh_vals!) - @eval function MAK.$f!(d::AbstractTensorMap, V, alg::DiagonalAlgorithm) - $f!(_repack_diagonal(d), diagview(_repack_diagonal(V)), alg) - return V - end -end - # For diagonal inputs we don't have to promote the scalartype since we know they are symmetric function MAK.initialize_output(::typeof(eig_vals!), t::AbstractTensorMap, alg::DiagonalAlgorithm) V_D = fuse(domain(t)) From b640d024555b96cf9effb61ac6a3f0f168af997f Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Mon, 22 Dec 2025 13:58:32 +0100 Subject: [PATCH 6/8] LA svd comment --- src/factorizations/factorizations.jl | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/factorizations/factorizations.jl b/src/factorizations/factorizations.jl index 6a23ef045..897af7396 100644 --- a/src/factorizations/factorizations.jl +++ b/src/factorizations/factorizations.jl @@ -49,10 +49,14 @@ function LinearAlgebra.eigvals(t::AbstractTensorMap; kwargs...) end LinearAlgebra.eigvals!(t::AbstractTensorMap; kwargs...) = eig_vals!(t) -LinearAlgebra.svd(t::AbstractTensorMap; full::Bool = false) = - full ? svd_full(t) : svd_compact(t) -LinearAlgebra.svd!(t::AbstractTensorMap; full::Bool = false) = - full ? svd_full!(t) : svd_compact!(t) +function LinearAlgebra.svd(t::AbstractTensorMap; full::Bool = false) + U, S, Vᴴ = full ? svd_full(t) : svd_compact(t) + return U, diagview(S), adjoint(Vᴴ) +end +function LinearAlgebra.svd!(t::AbstractTensorMap; full::Bool = false) + U, S, Vᴴ = full ? svd_full!(t) : svd_compact!(t) + return U, diagview(S), adjoint(Vᴴ) +end function LinearAlgebra.svdvals(t::AbstractTensorMap) tcopy = copy_oftype(t, factorisation_scalartype(svd_vals!, t)) From 867485b800991ed17c4e47e9eadf89f3e987d500 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Mon, 22 Dec 2025 20:42:00 +0100 Subject: [PATCH 7/8] Don't overload LinearAlgebra's svd for now --- src/factorizations/factorizations.jl | 9 --------- 1 file changed, 9 deletions(-) diff --git a/src/factorizations/factorizations.jl b/src/factorizations/factorizations.jl index 897af7396..0c2272ac8 100644 --- a/src/factorizations/factorizations.jl +++ b/src/factorizations/factorizations.jl @@ -49,15 +49,6 @@ function LinearAlgebra.eigvals(t::AbstractTensorMap; kwargs...) end LinearAlgebra.eigvals!(t::AbstractTensorMap; kwargs...) = eig_vals!(t) -function LinearAlgebra.svd(t::AbstractTensorMap; full::Bool = false) - U, S, Vᴴ = full ? svd_full(t) : svd_compact(t) - return U, diagview(S), adjoint(Vᴴ) -end -function LinearAlgebra.svd!(t::AbstractTensorMap; full::Bool = false) - U, S, Vᴴ = full ? svd_full!(t) : svd_compact!(t) - return U, diagview(S), adjoint(Vᴴ) -end - function LinearAlgebra.svdvals(t::AbstractTensorMap) tcopy = copy_oftype(t, factorisation_scalartype(svd_vals!, t)) return LinearAlgebra.svdvals!(tcopy) From b6f58f3def76a329703fe2488b5bc0a2c4100071 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Tue, 23 Dec 2025 07:20:13 +0100 Subject: [PATCH 8/8] Fix Changelog --- docs/src/Changelog.md | 1 - 1 file changed, 1 deletion(-) diff --git a/docs/src/Changelog.md b/docs/src/Changelog.md index d2ff74b51..2b5d6ca87 100644 --- a/docs/src/Changelog.md +++ b/docs/src/Changelog.md @@ -24,7 +24,6 @@ When releasing a new version, move the "Unreleased" changes to a new version sec - Extended support for selecting storage types in the `TensorMap` constructors ([#327](https://github.com/QuantumKitHub/TensorKit.jl/pull/327)) - `similar_diagonal` to handle storage types when constructing diagonals ([#330](https://github.com/QuantumKitHub/TensorKit.jl/pull/330)) -- `LinearAlgebra.svd` overloads ### Fixed