From db4c6269e26d9f0bcc150298be5e7a6e029b65d7 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Wed, 4 Jun 2025 18:05:06 -0400 Subject: [PATCH 1/2] [WIP] Fix truncated SVD --- Project.toml | 4 +-- src/factorizations.jl | 23 ------------ test/Project.toml | 2 +- test/test_factorizations.jl | 70 ++++++++++++++++++------------------- 4 files changed, 38 insertions(+), 61 deletions(-) diff --git a/Project.toml b/Project.toml index 669574fa..bdeb193e 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GradedArrays" uuid = "bc96ca6e-b7c8-4bb6-888e-c93f838762c2" authors = ["ITensor developers and contributors"] -version = "0.4.6" +version = "0.4.7" [deps] BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" @@ -25,7 +25,7 @@ GradedArraysTensorAlgebraExt = "TensorAlgebra" [compat] BlockArrays = "1.6.0" -BlockSparseArrays = "0.6.5" +BlockSparseArrays = "0.7.0" Compat = "4.16.0" DerivableInterfaces = "0.4.4" FillArrays = "1.13.0" diff --git a/src/factorizations.jl b/src/factorizations.jl index 70b23714..23fcd572 100644 --- a/src/factorizations.jl +++ b/src/factorizations.jl @@ -42,29 +42,6 @@ end const TGradedUSVᴴ = Tuple{<:GradedMatrix,<:GradedMatrix,<:GradedMatrix} -function BlockSparseArrays.similar_truncate( - ::typeof(svd_trunc!), - (U, S, Vᴴ)::TGradedUSVᴴ, - strategy::BlockPermutedDiagonalTruncationStrategy, - indexmask=MatrixAlgebraKit.findtruncated(diagview(S), strategy), -) - u_axis, v_axis = axes(S) - counter = Base.Fix1(count, Base.Fix1(getindex, indexmask)) - s_lengths = map(counter, blocks(u_axis)) - u_sectors = sectors(u_axis) .=> s_lengths - v_sectors = sectors(v_axis) .=> s_lengths - u_sectors_filtered = filter(>(0) ∘ last, u_sectors) - v_sectors_filtered = filter(>(0) ∘ last, v_sectors) - u_axis′ = gradedrange(u_sectors_filtered) - u_axis = isdual(u_axis) ? dual(u_axis′) : u_axis′ - v_axis′ = gradedrange(v_sectors_filtered) - v_axis = isdual(v_axis) ? dual(v_axis′) : v_axis′ - Ũ = similar(U, axes(U, 1), dual(u_axis)) - S̃ = similar(S, u_axis, v_axis) - Ṽᴴ = similar(Vᴴ, dual(v_axis), axes(Vᴴ, 2)) - return Ũ, S̃, Ṽᴴ -end - function BlockSparseArrays.similar_output( ::typeof(qr_compact!), A::GradedMatrix, R_axis, alg::BlockPermutedDiagonalAlgorithm ) diff --git a/test/Project.toml b/test/Project.toml index 4a3c679a..b3be449f 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -17,7 +17,7 @@ TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a" [compat] Aqua = "0.8.11" BlockArrays = "1.6.0" -BlockSparseArrays = "0.6" +BlockSparseArrays = "0.7.0" GradedArrays = "0.4" LinearAlgebra = "1.10.0" MatrixAlgebraKit = "0.2" diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl index 6ee5c1dd..d401dcef 100644 --- a/test/test_factorizations.jl +++ b/test/test_factorizations.jl @@ -82,41 +82,41 @@ end end end -@testset "svd_trunc (eltype=$elt)" for elt in elts - for i in [2, 3], j in [2, 3], k in [2, 3], l in [2, 3] - r1 = gradedrange([U1(0) => i, U1(1) => j]) - r2 = gradedrange([U1(0) => k, U1(1) => l]) - a = zeros(elt, r1, dual(r2)) - a[Block(2, 2)] = randn(elt, blocksizes(a)[2, 2]) - @test flux(a) == U1(0) - u, s, vᴴ = svd_trunc(a; trunc=(; maxrank=1)) - @test sort(diag(Matrix(s)); rev=true) ≈ svdvals(Matrix(a))[1:size(s, 1)] - @test size(u) == (size(a, 1), 1) - @test size(s) == (1, 1) - @test size(vᴴ) == (1, size(a, 2)) - @test Array(u'u) ≈ I - @test Array(vᴴ * vᴴ') ≈ I - @test flux(u) == trivial(flux(a)) - @test flux(s) == flux(a) - @test flux(vᴴ) == trivial(flux(a)) - - r1 = gradedrange([U1(0) => i, U1(1) => j]) - r2 = gradedrange([U1(0) => k, U1(1) => l]) - a = zeros(elt, r1, dual(r2)) - a[Block(1, 2)] = randn(elt, blocksizes(a)[1, 2]) - @test flux(a) == U1(-1) - u, s, vᴴ = svd_trunc(a; trunc=(; maxrank=1)) - @test sort(diag(Matrix(s)); rev=true) ≈ svdvals(Matrix(a))[1:size(s, 1)] - @test size(u) == (size(a, 1), 1) - @test size(s) == (1, 1) - @test size(vᴴ) == (1, size(a, 2)) - @test Array(u'u) ≈ I - @test Array(vᴴ * vᴴ') ≈ I - @test flux(u) == trivial(flux(a)) - @test flux(s) == flux(a) - @test flux(vᴴ) == trivial(flux(a)) - end -end +## @testset "svd_trunc (eltype=$elt)" for elt in elts +## for i in [2, 3], j in [2, 3], k in [2, 3], l in [2, 3] +## r1 = gradedrange([U1(0) => i, U1(1) => j]) +## r2 = gradedrange([U1(0) => k, U1(1) => l]) +## a = zeros(elt, r1, dual(r2)) +## a[Block(2, 2)] = randn(elt, blocksizes(a)[2, 2]) +## @test flux(a) == U1(0) +## u, s, vᴴ = svd_trunc(a; trunc=(; maxrank=1)) +## @test sort(diag(Matrix(s)); rev=true) ≈ svdvals(Matrix(a))[1:size(s, 1)] +## @test size(u) == (size(a, 1), 1) +## @test size(s) == (1, 1) +## @test size(vᴴ) == (1, size(a, 2)) +## @test Array(u'u) ≈ I +## @test Array(vᴴ * vᴴ') ≈ I +## @test flux(u) == trivial(flux(a)) +## @test flux(s) == flux(a) +## @test flux(vᴴ) == trivial(flux(a)) +## +## r1 = gradedrange([U1(0) => i, U1(1) => j]) +## r2 = gradedrange([U1(0) => k, U1(1) => l]) +## a = zeros(elt, r1, dual(r2)) +## a[Block(1, 2)] = randn(elt, blocksizes(a)[1, 2]) +## @test flux(a) == U1(-1) +## u, s, vᴴ = svd_trunc(a; trunc=(; maxrank=1)) +## @test sort(diag(Matrix(s)); rev=true) ≈ svdvals(Matrix(a))[1:size(s, 1)] +## @test size(u) == (size(a, 1), 1) +## @test size(s) == (1, 1) +## @test size(vᴴ) == (1, size(a, 2)) +## @test Array(u'u) ≈ I +## @test Array(vᴴ * vᴴ') ≈ I +## @test flux(u) == trivial(flux(a)) +## @test flux(s) == flux(a) +## @test flux(vᴴ) == trivial(flux(a)) +## end +## end @testset "qr_compact, left_orth (eltype=$elt)" for elt in elts for i in [2, 3], j in [2, 3], k in [2, 3], l in [2, 3] From 5788a6d5bea101a40f2c24668fcb67d833aca364 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Wed, 4 Jun 2025 18:47:36 -0400 Subject: [PATCH 2/2] Fix logical indexing of graded arrays --- src/gradedunitrange.jl | 18 ++++++++++ test/test_factorizations.jl | 70 ++++++++++++++++++------------------- 2 files changed, 53 insertions(+), 35 deletions(-) diff --git a/src/gradedunitrange.jl b/src/gradedunitrange.jl index 79ab6aeb..046a27eb 100644 --- a/src/gradedunitrange.jl +++ b/src/gradedunitrange.jl @@ -320,3 +320,21 @@ function BlockSparseArrays.blockedunitrange_getindices( new_range = blockedrange(new_first, length.(new_axes)) return GradedUnitRange(new_axes, new_range) end + +using BlockArrays: BlockedVector +function BlockSparseArrays.blockedunitrange_getindices( + a::AbstractGradedUnitRange, indices::AbstractVector{Bool} +) + blocked_indices = BlockedVector(indices, axes(a)) + bs = map(Base.OneTo(blocklength(blocked_indices))) do b + binds = blocked_indices[Block(b)] + bstart = blockfirsts(only(axes(blocked_indices)))[b] + return findall(binds) .+ (bstart - 1) + end + keep = map(!isempty, bs) + secs = sectors(a)[keep] + bs = bs[keep] + r = gradedrange(secs .=> length.(bs); isdual=isdual(a)) + I = mortar(bs, (r,)) + return I +end diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl index d401dcef..6ee5c1dd 100644 --- a/test/test_factorizations.jl +++ b/test/test_factorizations.jl @@ -82,41 +82,41 @@ end end end -## @testset "svd_trunc (eltype=$elt)" for elt in elts -## for i in [2, 3], j in [2, 3], k in [2, 3], l in [2, 3] -## r1 = gradedrange([U1(0) => i, U1(1) => j]) -## r2 = gradedrange([U1(0) => k, U1(1) => l]) -## a = zeros(elt, r1, dual(r2)) -## a[Block(2, 2)] = randn(elt, blocksizes(a)[2, 2]) -## @test flux(a) == U1(0) -## u, s, vᴴ = svd_trunc(a; trunc=(; maxrank=1)) -## @test sort(diag(Matrix(s)); rev=true) ≈ svdvals(Matrix(a))[1:size(s, 1)] -## @test size(u) == (size(a, 1), 1) -## @test size(s) == (1, 1) -## @test size(vᴴ) == (1, size(a, 2)) -## @test Array(u'u) ≈ I -## @test Array(vᴴ * vᴴ') ≈ I -## @test flux(u) == trivial(flux(a)) -## @test flux(s) == flux(a) -## @test flux(vᴴ) == trivial(flux(a)) -## -## r1 = gradedrange([U1(0) => i, U1(1) => j]) -## r2 = gradedrange([U1(0) => k, U1(1) => l]) -## a = zeros(elt, r1, dual(r2)) -## a[Block(1, 2)] = randn(elt, blocksizes(a)[1, 2]) -## @test flux(a) == U1(-1) -## u, s, vᴴ = svd_trunc(a; trunc=(; maxrank=1)) -## @test sort(diag(Matrix(s)); rev=true) ≈ svdvals(Matrix(a))[1:size(s, 1)] -## @test size(u) == (size(a, 1), 1) -## @test size(s) == (1, 1) -## @test size(vᴴ) == (1, size(a, 2)) -## @test Array(u'u) ≈ I -## @test Array(vᴴ * vᴴ') ≈ I -## @test flux(u) == trivial(flux(a)) -## @test flux(s) == flux(a) -## @test flux(vᴴ) == trivial(flux(a)) -## end -## end +@testset "svd_trunc (eltype=$elt)" for elt in elts + for i in [2, 3], j in [2, 3], k in [2, 3], l in [2, 3] + r1 = gradedrange([U1(0) => i, U1(1) => j]) + r2 = gradedrange([U1(0) => k, U1(1) => l]) + a = zeros(elt, r1, dual(r2)) + a[Block(2, 2)] = randn(elt, blocksizes(a)[2, 2]) + @test flux(a) == U1(0) + u, s, vᴴ = svd_trunc(a; trunc=(; maxrank=1)) + @test sort(diag(Matrix(s)); rev=true) ≈ svdvals(Matrix(a))[1:size(s, 1)] + @test size(u) == (size(a, 1), 1) + @test size(s) == (1, 1) + @test size(vᴴ) == (1, size(a, 2)) + @test Array(u'u) ≈ I + @test Array(vᴴ * vᴴ') ≈ I + @test flux(u) == trivial(flux(a)) + @test flux(s) == flux(a) + @test flux(vᴴ) == trivial(flux(a)) + + r1 = gradedrange([U1(0) => i, U1(1) => j]) + r2 = gradedrange([U1(0) => k, U1(1) => l]) + a = zeros(elt, r1, dual(r2)) + a[Block(1, 2)] = randn(elt, blocksizes(a)[1, 2]) + @test flux(a) == U1(-1) + u, s, vᴴ = svd_trunc(a; trunc=(; maxrank=1)) + @test sort(diag(Matrix(s)); rev=true) ≈ svdvals(Matrix(a))[1:size(s, 1)] + @test size(u) == (size(a, 1), 1) + @test size(s) == (1, 1) + @test size(vᴴ) == (1, size(a, 2)) + @test Array(u'u) ≈ I + @test Array(vᴴ * vᴴ') ≈ I + @test flux(u) == trivial(flux(a)) + @test flux(s) == flux(a) + @test flux(vᴴ) == trivial(flux(a)) + end +end @testset "qr_compact, left_orth (eltype=$elt)" for elt in elts for i in [2, 3], j in [2, 3], k in [2, 3], l in [2, 3]