From d0df70b2366dd31f55977b71c9ab6f8ab8f65a20 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Wed, 28 May 2025 09:31:18 -0400 Subject: [PATCH] Implement graded LQ --- Project.toml | 4 +- src/factorizations.jl | 25 ++++++- test/test_factorizations.jl | 145 ++++++++++++++++++++++++++++++++---- 3 files changed, 158 insertions(+), 16 deletions(-) diff --git a/Project.toml b/Project.toml index 2846bba..669574f 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.5" +version = "0.4.6" [deps] BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" @@ -25,7 +25,7 @@ GradedArraysTensorAlgebraExt = "TensorAlgebra" [compat] BlockArrays = "1.6.0" -BlockSparseArrays = "0.6.2" +BlockSparseArrays = "0.6.5" Compat = "4.16.0" DerivableInterfaces = "0.4.4" FillArrays = "1.13.0" diff --git a/src/factorizations.jl b/src/factorizations.jl index c16f44f..70b2371 100644 --- a/src/factorizations.jl +++ b/src/factorizations.jl @@ -9,7 +9,14 @@ using BlockSparseArrays: mortar_axis using LinearAlgebra: Diagonal using MatrixAlgebraKit: - MatrixAlgebraKit, qr_compact!, qr_full!, svd_compact!, svd_full!, svd_trunc! + MatrixAlgebraKit, + lq_compact!, + lq_full!, + qr_compact!, + qr_full!, + svd_compact!, + svd_full!, + svd_trunc! function BlockSparseArrays.similar_output( ::typeof(svd_compact!), A::GradedMatrix, S_axes, alg::BlockPermutedDiagonalAlgorithm @@ -73,3 +80,19 @@ function BlockSparseArrays.similar_output( R = similar(A, R_axis, axes(A, 2)) return Q, R end + +function BlockSparseArrays.similar_output( + ::typeof(lq_compact!), A::GradedMatrix, L_axis, alg::BlockPermutedDiagonalAlgorithm +) + L = similar(A, axes(A, 1), L_axis) + Q = similar(A, dual(L_axis), axes(A, 2)) + return L, Q +end + +function BlockSparseArrays.similar_output( + ::typeof(lq_full!), A::GradedMatrix, L_axis, alg::BlockPermutedDiagonalAlgorithm +) + L = similar(A, axes(A, 1), L_axis) + Q = similar(A, dual(L_axis), axes(A, 2)) + return L, Q +end diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl index 9ef8d97..6ee5c1d 100644 --- a/test/test_factorizations.jl +++ b/test/test_factorizations.jl @@ -1,8 +1,19 @@ using BlockArrays: Block, blocksizes using GradedArrays: U1, dual, flux, gradedrange, trivial using LinearAlgebra: I, diag, svdvals -using MatrixAlgebraKit: qr_compact, qr_full, svd_compact, svd_full, svd_trunc -using Test: @test, @testset +using MatrixAlgebraKit: + left_orth, + left_polar, + lq_compact, + lq_full, + qr_compact, + qr_full, + right_orth, + right_polar, + svd_compact, + svd_full, + svd_trunc +using Test: @test, @test_broken, @testset const elts = (Float32, Float64, ComplexF32, ComplexF64) @testset "svd_compact (eltype=$elt)" for elt in elts @@ -107,29 +118,33 @@ end end end -@testset "qr_compact (eltype=$elt)" for elt in elts +@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] 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) - q, r = qr_compact(a) - @test q * r ≈ a - @test Array(q'q) ≈ I - @test flux(q) == trivial(flux(a)) - @test flux(r) == flux(a) + for f in (qr_compact, left_orth) + q, r = f(a) + @test q * r ≈ a + @test Array(q'q) ≈ I + @test flux(q) == trivial(flux(a)) + @test flux(r) == flux(a) + end 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) - q, r = qr_compact(a) - @test q * r ≈ a - @test Array(q'q) ≈ I - @test flux(q) == trivial(flux(a)) - @test flux(r) == flux(a) + for f in (qr_compact, left_orth) + q, r = f(a) + @test q * r ≈ a + @test Array(q'q) ≈ I + @test flux(q) == trivial(flux(a)) + @test flux(r) == flux(a) + end end end @@ -160,3 +175,107 @@ end @test flux(r) == flux(a) end end + +@testset "left_polar (eltype=$elt)" for elt in elts + r1 = gradedrange([U1(0) => 3, U1(1) => 4]) + r2 = gradedrange([U1(0) => 2, U1(1) => 3]) + a = zeros(elt, r1, dual(r2)) + a[Block(2, 2)] = randn(elt, blocksizes(a)[2, 2]) + @test flux(a) == U1(0) + q, r = left_polar(a) + @test q * r ≈ a + @test Array(q'q) ≈ I + @test flux(q) == trivial(flux(a)) + @test flux(r) == flux(a) + + r1 = gradedrange([U1(0) => 3, U1(1) => 4]) + r2 = gradedrange([U1(0) => 2, U1(1) => 3]) + a = zeros(elt, r1, dual(r2)) + a[Block(1, 2)] = randn(elt, blocksizes(a)[1, 2]) + @test flux(a) == U1(-1) + q, r = left_polar(a) + @test q * r ≈ a + @test Array(q'q) ≈ I + @test_broken flux(q) == trivial(flux(a)) + @test_broken flux(r) == flux(a) +end + +@testset "lq_compact, right_orth (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) + l, q = lq_compact(a) + @test l * q ≈ a + @test Array(q * q') ≈ I + @test flux(l) == flux(a) + @test flux(q) == trivial(flux(a)) + end + 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(1, 2)] = randn(elt, blocksizes(a)[1, 2]) + @test flux(a) == U1(-1) + l, q = lq_compact(a) + @test l * q ≈ a + @test Array(q * q') ≈ I + @test flux(l) == flux(a) + @test flux(q) == trivial(flux(a)) + end +end + +@testset "lq_full (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) + l, q = lq_full(a) + @test l * q ≈ a + @test Array(q * q') ≈ I + @test Array(q'q) ≈ I + @test flux(l) == flux(a) + @test flux(q) == trivial(flux(a)) + end + 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(1, 2)] = randn(elt, blocksizes(a)[1, 2]) + @test flux(a) == U1(-1) + l, q = lq_full(a) + @test l * q ≈ a + @test Array(q * q') ≈ I + @test Array(q'q) ≈ I + @test flux(l) == flux(a) + @test flux(q) == trivial(flux(a)) + end +end + +@testset "right_polar (eltype=$elt)" for elt in elts + r1 = gradedrange([U1(0) => 2, U1(1) => 3]) + r2 = gradedrange([U1(0) => 3, U1(1) => 4]) + a = zeros(elt, r1, dual(r2)) + a[Block(2, 2)] = randn(elt, blocksizes(a)[2, 2]) + @test flux(a) == U1(0) + l, q = right_polar(a) + @test l * q ≈ a + @test Array(q * q') ≈ I + @test flux(l) == flux(a) + @test flux(q) == trivial(flux(a)) + + r1 = gradedrange([U1(0) => 2, U1(1) => 3]) + r2 = gradedrange([U1(0) => 3, U1(1) => 4]) + a = zeros(elt, r1, dual(r2)) + a[Block(1, 2)] = randn(elt, blocksizes(a)[1, 2]) + @test flux(a) == U1(-1) + l, q = right_polar(a) + @test l * q ≈ a + @test Array(q * q') ≈ I + @test_broken flux(l) == flux(a) + @test_broken flux(q) == trivial(flux(a)) +end