Implement graded LQ (#38)

Author: mtfishman

Project.toml

@@ -1,7 +1,7 @@
 name = "GradedArrays"
 uuid = "bc96ca6e-b7c8-4bb6-888e-c93f838762c2"
 authors = ["ITensor developers <[email protected]> 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"

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

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