Fix polar

Author: mtfishman

ext/MatrixAlgebraKitFillArraysExt.jl

@@ -14,7 +14,7 @@ end
 
 struct ZerosAlgorithm <: AbstractAlgorithm end
 
-for f in [:eig, :eigh, :lq, :qr, :svd]
+for f in [:eig, :eigh, :lq, :polar, :qr, :svd]
     ff = Symbol("default_", f, "_algorithm")
     @eval begin
         function MatrixAlgebraKit.$ff(::Type{<:AbstractZerosMatrix}; kwargs...)
@@ -170,6 +170,30 @@ function MatrixAlgebraKit.svd_vals!(A::AbstractZerosMatrix, F, alg::ZerosAlgorit
     return diagview(A)
 end
 
+function MatrixAlgebraKit.check_input(::typeof(left_polar!), A::AbstractZerosMatrix, F)
+    m, n = size(A)
+    m >= n ||
+        throw(ArgumentError("input matrix needs at least as many rows as columns"))
+    return nothing
+end
+function MatrixAlgebraKit.left_polar!(A::AbstractZerosMatrix, F, alg::ZerosAlgorithm)
+    check_input(left_polar!, A, F)
+    U, S, Vᴴ = svd_compact(A)
+    return (Eye((axes(U, 1), axes(Vᴴ, 2))), Vᴴ' * S * Vᴴ)
+end
+
+function MatrixAlgebraKit.check_input(::typeof(right_polar!), A::AbstractZerosMatrix, F)
+    m, n = size(A)
+    n >= m ||
+        throw(ArgumentError("input matrix needs at least as many columns as rows"))
+    return nothing
+end
+function MatrixAlgebraKit.right_polar!(A::AbstractZerosMatrix, F, alg::ZerosAlgorithm)
+    check_input(right_polar!, A, F)
+    U, S, Vᴴ = svd_compact(A)
+    return (U * S * U', Eye((axes(U, 1), axes(Vᴴ, 2))))
+end
+
 struct EyeAlgorithm <: AbstractAlgorithm end
 
 for f in [:eig, :eigh, :lq, :qr, :polar, :svd]

test/fillarrays.jl

@@ -64,15 +64,15 @@ using FillArrays: SquareEye
     @test iszero(R)
     @test R isa Zeros
 
-    # A = Zeros(4, 3)
-    # Q, R = @constinferred left_polar(A)
-    # @test Q * R == A
-    # @test size(Q) == (4, 3)
-    # @test size(R) == (3, 3)
-    # @test Q == Matrix(I, (4, 3))
-    # @test Q isa Eye
-    # @test iszero(R)
-    # @test R isa Zeros
+    A = Zeros(4, 3)
+    Q, R = @constinferred left_polar(A)
+    @test Q * R == A
+    @test size(Q) == (4, 3)
+    @test size(R) == (3, 3)
+    @test Q == Matrix(I, (4, 3))
+    @test Q isa Eye
+    @test iszero(R)
+    @test R isa Zeros
 
     A = Zeros(3, 4)
     L, Q = @constinferred lq_compact(A)
@@ -109,15 +109,15 @@ using FillArrays: SquareEye
     @test Q == Eye(4)
     @test Q isa Eye
 
-    # A = Zeros(3, 4)
-    # L, Q = @constinferred right_polar(A)
-    # @test L * Q == A
-    # @test size(L) == (3, 3)
-    # @test size(Q) == (3, 4)
-    # @test iszero(L)
-    # @test L isa Zeros
-    # @test Q == Matrix(I, (3, 4))
-    # @test Q isa Eye
+    A = Zeros(3, 4)
+    L, Q = @constinferred right_polar(A)
+    @test L * Q == A
+    @test size(L) == (3, 3)
+    @test size(Q) == (3, 4)
+    @test iszero(L)
+    @test L isa Zeros
+    @test Q == Matrix(I, (3, 4))
+    @test Q isa Eye
 
     A = Zeros(3, 4)
     U, S, V = @constinferred svd_compact(A)