Add polar, add tests

Author: mtfishman

src/BlockSparseArrays.jl

@@ -47,6 +47,7 @@ include("BlockArraysSparseArraysBaseExt/BlockArraysSparseArraysBaseExt.jl")
 include("factorizations/svd.jl")
 include("factorizations/truncation.jl")
 include("factorizations/qr.jl")
+include("factorizations/polar.jl")
 include("factorizations/orthnull.jl")
 
 end

src/factorizations/orthnull.jl

@@ -1,5 +1,11 @@
 using MatrixAlgebraKit:
-  MatrixAlgebraKit, left_polar, qr_compact, select_algorithm, svd_compact
+  MatrixAlgebraKit,
+  left_polar,
+  lq_compact,
+  qr_compact,
+  right_polar,
+  select_algorithm,
+  svd_compact
 
 function MatrixAlgebraKit.left_orth(
   A::AbstractBlockSparseMatrix;
@@ -15,9 +21,7 @@ function MatrixAlgebraKit.left_orth(
   if kind == :qr
     return left_orth_qr(A, alg_qr)
   elseif kind == :polar
-    # TODO: Implement this.
-    # return left_orth_polar(A, alg_polar)
-    return left_orth_svd(A, alg_svd)
+    return left_orth_polar(A, alg_polar)
   elseif kind == :svd
     return left_orth_svd(A, alg_svd, trunc)
   else
@@ -39,7 +43,7 @@ function left_orth_svd(A, alg, trunc::Nothing=nothing)
 end
 
 function MatrixAlgebraKit.right_orth(
-  A;
+  A::AbstractBlockSparseMatrix;
   trunc=nothing,
   kind=isnothing(trunc) ? :lq : :svd,
   alg_lq=(; positive=true),
@@ -54,9 +58,7 @@ function MatrixAlgebraKit.right_orth(
     # return right_orth_lq(A, alg_lq)
     return right_orth_svd(A, alg_svd)
   elseif kind == :polar
-    # TODO: Implement this.
-    # return right_orth_polar(A, alg_polar)
-    return right_orth_svd(A, alg_svd)
+    return right_orth_polar(A, alg_polar)
   elseif kind == :svd
     return right_orth_svd(A, alg_svd, trunc)
   else

src/factorizations/polar.jl

@@ -0,0 +1,49 @@
+using MatrixAlgebraKit:
+  MatrixAlgebraKit,
+  PolarViaSVD,
+  check_input,
+  default_algorithm,
+  left_polar!,
+  right_polar!,
+  svd_compact!
+
+function MatrixAlgebraKit.check_input(
+  ::typeof(left_polar!), A::AbstractBlockSparseMatrix, WP
+)
+  W, P = WP
+  @views for I in eachblockstoredindex(A)
+    m, n = size(A[I])
+    m >= n ||
+      throw(ArgumentError("each input matrix block needs at least as many rows as columns"))
+    # check_input(left_polar!, A[I], (W[I1], P[I2]))
+  end
+  return nothing
+end
+
+function MatrixAlgebraKit.left_polar!(A::AbstractBlockSparseMatrix, WP, alg::PolarViaSVD)
+  check_input(left_polar!, A, WP)
+  U, S, Vᴴ = svd_compact!(A, alg.svdalg)
+  # TODO: Use more in-place operations here, avoid `copy`.
+  W = U * Vᴴ
+  P = copy(Vᴴ') * S * Vᴴ
+  return (W, P)
+end
+function MatrixAlgebraKit.right_polar!(A::AbstractBlockSparseMatrix, PWᴴ, alg::PolarViaSVD)
+  check_input(right_polar!, A, PWᴴ)
+  U, S, Vᴴ = svd_compact!(A, alg.svdalg)
+  # TODO: Use more in-place operations here, avoid `copy`.
+  Wᴴ = U * Vᴴ
+  P = U * S * copy(U')
+  return (P, Wᴴ)
+end
+
+function MatrixAlgebraKit.default_algorithm(
+  ::typeof(left_polar!), a::AbstractBlockSparseMatrix; kwargs...
+)
+  return PolarViaSVD(default_algorithm(svd_compact!, a; kwargs...))
+end
+function MatrixAlgebraKit.default_algorithm(
+  ::typeof(right_polar!), a::AbstractBlockSparseMatrix; kwargs...
+)
+  return PolarViaSVD(default_algorithm(svd_compact!, a; kwargs...))
+end

test/test_factorizations.jl

@@ -1,7 +1,17 @@
 using BlockArrays: Block, BlockedMatrix, BlockedVector, blocks, mortar
 using BlockSparseArrays: BlockSparseArray, BlockDiagonal, eachblockstoredindex
 using MatrixAlgebraKit:
-  qr_compact, qr_full, svd_compact, svd_full, svd_trunc, truncrank, trunctol
+  left_orth,
+  left_polar,
+  qr_compact,
+  qr_full,
+  right_orth,
+  right_polar,
+  svd_compact,
+  svd_full,
+  svd_trunc,
+  truncrank,
+  trunctol
 using LinearAlgebra: LinearAlgebra
 using Random: Random
 using Test: @inferred, @testset, @test
@@ -156,7 +166,7 @@ end
   end
 end
 
-@testset "qr_compact" for T in (Float32, Float64, ComplexF32, ComplexF64)
+@testset "qr_compact (T=$T)" for T in (Float32, Float64, ComplexF32, ComplexF64)
   for i in [1, 2], j in [1, 2], k in [1, 2], l in [1, 2]
     A = BlockSparseArray{T}(undef, ([i, j], [k, l]))
     A[Block(1, 1)] = randn(T, i, k)
@@ -167,7 +177,7 @@ end
   end
 end
 
-@testset "qr_full" for T in (Float32, Float64, ComplexF32, ComplexF64)
+@testset "qr_full (T=$T)" for T in (Float32, Float64, ComplexF32, ComplexF64)
   for i in [2, 3], j in [2, 3], k in [2, 3], l in [2, 3]
     A = BlockSparseArray{T}(undef, ([i, j], [k, l]))
     A[Block(1, 1)] = randn(T, i, k)
@@ -181,3 +191,47 @@ end
     @test A ≈ Q * R
   end
 end
+
+@testset "left_polar (T=$T)" for T in (Float32, Float64, ComplexF32, ComplexF64)
+  A = BlockSparseArray{T}(undef, ([3, 4], [2, 3]))
+  A[Block(1, 1)] = randn(T, 3, 2)
+  A[Block(2, 2)] = randn(T, 4, 3)
+
+  U, C = left_polar(A)
+  @test U * C ≈ A
+  @test Matrix(U'U) ≈ LinearAlgebra.I
+end
+
+@testset "right_polar (T=$T)" for T in (Float32, Float64, ComplexF32, ComplexF64)
+  A = BlockSparseArray{T}(undef, ([2, 3], [3, 4]))
+  A[Block(1, 1)] = randn(T, 2, 3)
+  A[Block(2, 2)] = randn(T, 3, 4)
+
+  C, U = right_polar(A)
+  @test C * U ≈ A
+  @test Matrix(U * U') ≈ LinearAlgebra.I
+end
+
+@testset "left_orth (T=$T)" for T in (Float32, Float64, ComplexF32, ComplexF64)
+  A = BlockSparseArray{T}(undef, ([3, 4], [2, 3]))
+  A[Block(1, 1)] = randn(T, 3, 2)
+  A[Block(2, 2)] = randn(T, 4, 3)
+
+  for kind in (:qr, :polar, :svd)
+    U, C = left_orth(A; kind)
+    @test U * C ≈ A
+    @test Matrix(U'U) ≈ LinearAlgebra.I
+  end
+end
+
+@testset "right_orth (T=$T)" for T in (Float32, Float64, ComplexF32, ComplexF64)
+  A = BlockSparseArray{T}(undef, ([2, 3], [3, 4]))
+  A[Block(1, 1)] = randn(T, 2, 3)
+  A[Block(2, 2)] = randn(T, 3, 4)
+
+  for kind in (:qr, :polar, :svd)
+    C, U = right_orth(A; kind)
+    @test C * U ≈ A
+    @test Matrix(U * U') ≈ LinearAlgebra.I
+  end
+end