Merge branch 'main' into copilot/update-docstrings-eigenvalue-decompositions

Author: lkdvos

ext/MatrixAlgebraKitChainRulesCoreExt.jl

@@ -14,7 +14,7 @@ MatrixAlgebraKit.iszerotangent(::AbstractZero) = true
 @non_differentiable MatrixAlgebraKit.select_algorithm(args...)
 @non_differentiable MatrixAlgebraKit.initialize_output(args...)
 @non_differentiable MatrixAlgebraKit.check_input(args...)
-@non_differentiable MatrixAlgebraKit.isisometry(args...)
+@non_differentiable MatrixAlgebraKit.isisometric(args...)
 @non_differentiable MatrixAlgebraKit.isunitary(args...)
 
 function ChainRulesCore.rrule(::typeof(copy_input), f, A)

src/MatrixAlgebraKit.jl

@@ -3,16 +3,16 @@ module MatrixAlgebraKit
 using LinearAlgebra: LinearAlgebra
 using LinearAlgebra: norm # TODO: eleminate if we use VectorInterface.jl?
 using LinearAlgebra: mul!, rmul!, lmul!, adjoint!, rdiv!, ldiv!
-using LinearAlgebra: sylvester
+using LinearAlgebra: sylvester, lu!
 using LinearAlgebra: isposdef, issymmetric
 using LinearAlgebra: Diagonal, diag, diagind, isdiag
 using LinearAlgebra: UpperTriangular, LowerTriangular
 using LinearAlgebra: BlasFloat, BlasReal, BlasComplex, BlasInt
 
-export isisometry, isunitary, ishermitian, isantihermitian
+export isisometric, isunitary, ishermitian, isantihermitian
 
-export project_hermitian, project_antihermitian
-export project_hermitian!, project_antihermitian!
+export project_hermitian, project_antihermitian, project_isometric
+export project_hermitian!, project_antihermitian!, project_isometric!
 export qr_compact, qr_full, qr_null, lq_compact, lq_full, lq_null
 export qr_compact!, qr_full!, qr_null!, lq_compact!, lq_full!, lq_null!
 export svd_compact, svd_full, svd_vals, svd_trunc
@@ -34,6 +34,7 @@ export LAPACK_HouseholderQR, LAPACK_HouseholderLQ, LAPACK_Simple, LAPACK_Expert,
     LAPACK_QRIteration, LAPACK_Bisection, LAPACK_MultipleRelativelyRobustRepresentations,
     LAPACK_DivideAndConquer, LAPACK_Jacobi
 export LQViaTransposedQR
+export PolarViaSVD, PolarNewton
 export DiagonalAlgorithm
 export NativeBlocked
 export CUSOLVER_Simple, CUSOLVER_HouseholderQR, CUSOLVER_QRIteration, CUSOLVER_SVDPolar,
@@ -64,6 +65,7 @@ export notrunc, truncrank, trunctol, truncerror, truncfilter
             :svd_pullback!, :svd_trunc_pullback!
         )
     )
+    eval(Expr(:public, :is_left_isometric, :is_right_isometric))
 end
 
 include("common/defaults.jl")

src/common/matrixproperties.jl

@@ -1,5 +1,5 @@
 """
-    isisometry(A; side=:left, isapprox_kwargs...) -> Bool
+    isisometric(A; side=:left, isapprox_kwargs...) -> Bool
 
 Test whether a linear map is an isometry, where the type of isometry is controlled by `kind`:
 
@@ -8,13 +8,14 @@ Test whether a linear map is an isometry, where the type of isometry is controll
 
 The `isapprox_kwargs` are passed on to `isapprox` to control the tolerances.
 
-New specializations should overload [`is_left_isometry`](@ref) and [`is_right_isometry`](@ref).
+New specializations should overload [`MatrixAlgebraKit.is_left_isometric`](@ref) and
+[`MatrixAlgebraKit.is_right_isometric`](@ref).
 
 See also [`isunitary`](@ref).
 """
-function isisometry(A; side::Symbol = :left, isapprox_kwargs...)
-    side === :left && return is_left_isometry(A; isapprox_kwargs...)
-    side === :right && return is_right_isometry(A; isapprox_kwargs...)
+function isisometric(A; side::Symbol = :left, isapprox_kwargs...)
+    side === :left && return is_left_isometric(A; isapprox_kwargs...)
+    side === :right && return is_right_isometric(A; isapprox_kwargs...)
 
     throw(ArgumentError(lazy"Invalid isometry side: $side"))
 end
@@ -25,48 +26,42 @@ end
 Test whether a linear map is unitary, i.e. `A * A' ≈ I ≈ A' * A`.
 The `isapprox_kwargs` are passed on to `isapprox` to control the tolerances.
 
-See also [`isisometry`](@ref).
+See also [`isisometric`](@ref).
 """
 function isunitary(A; isapprox_kwargs...)
-    return is_left_isometry(A; isapprox_kwargs...) &&
-        is_right_isometry(A; isapprox_kwargs...)
+    return is_left_isometric(A; isapprox_kwargs...) &&
+        is_right_isometric(A; isapprox_kwargs...)
 end
 function isunitary(A::AbstractMatrix; isapprox_kwargs...)
     size(A, 1) == size(A, 2) || return false
-    return is_left_isometry(A; isapprox_kwargs...)
+    return is_left_isometric(A; isapprox_kwargs...)
 end
 
 @doc """
-    is_left_isometry(A; isapprox_kwargs...) -> Bool
+    is_left_isometric(A; isapprox_kwargs...) -> Bool
 
-Test whether a linear map is a left isometry, i.e. `A' * A ≈ I`.
+Test whether a linear map is a (left) isometry, i.e. `A' * A ≈ I`.
 The `isapprox_kwargs` can be used to control the tolerances of the equality.
 
-See also [`isisometry`](@ref) and [`is_right_isometry`](@ref).
-""" is_left_isometry
+See also [`isisometric`](@ref) and [`MatrixAlgebraKit.is_right_isometric`](@ref).
+""" is_left_isometric
 
-function is_left_isometry(A::AbstractMatrix; atol::Real = 0, rtol::Real = defaulttol(A), norm = LinearAlgebra.norm)
+function is_left_isometric(A::AbstractMatrix; atol::Real = 0, rtol::Real = defaulttol(A), norm = LinearAlgebra.norm)
     P = A' * A
     nP = norm(P) # isapprox would use `rtol * max(norm(P), norm(I))`
     diagview(P) .-= 1
     return norm(P) <= max(atol, rtol * nP) # assume that the norm of I is `sqrt(n)`
 end
 
 @doc """
-    is_right_isometry(A; isapprox_kwargs...) -> Bool
+    is_right_isometric(A; isapprox_kwargs...) -> Bool
 
-Test whether a linear map is a right isometry, i.e. `A * A' ≈ I`.
+Test whether a linear map is a (right) isometry, i.e. `A * A' ≈ I`.
 The `isapprox_kwargs` can be used to control the tolerances of the equality.
 
-See also [`isisometry`](@ref) and [`is_left_isometry`](@ref).
-""" is_right_isometry
-
-function is_right_isometry(A::AbstractMatrix; atol::Real = 0, rtol::Real = defaulttol(A), norm = LinearAlgebra.norm)
-    P = A * A'
-    nP = norm(P) # isapprox would use `rtol * max(norm(P), norm(I))`
-    diagview(P) .-= 1
-    return norm(P) <= max(atol, rtol * nP) # assume that the norm of I is `sqrt(n)`
-end
+See also [`isisometric`](@ref) and [`MatrixAlgebraKit.is_left_isometric`](@ref).
+""" is_right_isometric
+is_right_isometric(A; kwargs...) = is_left_isometric(A'; kwargs...)
 
 """
     ishermitian(A; isapprox_kwargs...)

src/implementations/polar.jl

@@ -11,7 +11,7 @@ function check_input(::typeof(left_polar!), A::AbstractMatrix, WP, ::AbstractAlg
     @assert W isa AbstractMatrix && P isa AbstractMatrix
     @check_size(W, (m, n))
     @check_scalar(W, A)
-    @check_size(P, (n, n))
+    isempty(P) || @check_size(P, (n, n))
     @check_scalar(P, A)
     return nothing
 end
@@ -21,7 +21,7 @@ function check_input(::typeof(right_polar!), A::AbstractMatrix, PWᴴ, ::Abstrac
     n >= m ||
         throw(ArgumentError("input matrix needs at least as many columns as rows"))
     @assert P isa AbstractMatrix && Wᴴ isa AbstractMatrix
-    @check_size(P, (m, m))
+    isempty(P) || @check_size(P, (m, m))
     @check_scalar(P, A)
     @check_size(Wᴴ, (m, n))
     @check_scalar(Wᴴ, A)
@@ -43,25 +43,154 @@ function initialize_output(::typeof(right_polar!), A::AbstractMatrix, ::Abstract
     return (P, Wᴴ)
 end
 
-# Implementation
-# --------------
+# Implementation via SVD
+# -----------------------
 function left_polar!(A::AbstractMatrix, WP, alg::PolarViaSVD)
     check_input(left_polar!, A, WP, alg)
-    U, S, Vᴴ = svd_compact!(A, alg.svdalg)
+    U, S, Vᴴ = svd_compact!(A, alg.svd_alg)
     W, P = WP
     W = mul!(W, U, Vᴴ)
-    S .= sqrt.(S)
-    SsqrtVᴴ = lmul!(S, Vᴴ)
-    P = mul!(P, SsqrtVᴴ', SsqrtVᴴ)
+    if !isempty(P)
+        S .= sqrt.(S)
+        SsqrtVᴴ = lmul!(S, Vᴴ)
+        P = mul!(P, SsqrtVᴴ', SsqrtVᴴ)
+    end
     return (W, P)
 end
 function right_polar!(A::AbstractMatrix, PWᴴ, alg::PolarViaSVD)
     check_input(right_polar!, A, PWᴴ, alg)
-    U, S, Vᴴ = svd_compact!(A, alg.svdalg)
+    U, S, Vᴴ = svd_compact!(A, alg.svd_alg)
     P, Wᴴ = PWᴴ
     Wᴴ = mul!(Wᴴ, U, Vᴴ)
-    S .= sqrt.(S)
-    USsqrt = rmul!(U, S)
-    P = mul!(P, USsqrt, USsqrt')
+    if !isempty(P)
+        S .= sqrt.(S)
+        USsqrt = rmul!(U, S)
+        P = mul!(P, USsqrt, USsqrt')
+    end
     return (P, Wᴴ)
 end
+
+# Implementation via Newton
+# --------------------------
+function left_polar!(A::AbstractMatrix, WP, alg::PolarNewton)
+    check_input(left_polar!, A, WP, alg)
+    W, P = WP
+    if isempty(P)
+        W = _left_polarnewton!(A, W, P; alg.kwargs...)
+        return W, P
+    else
+        W = _left_polarnewton!(copy(A), W, P; alg.kwargs...)
+        # we still need `A` to compute `P`
+        P = project_hermitian!(mul!(P, W', A))
+        return W, P
+    end
+end
+
+function right_polar!(A::AbstractMatrix, PWᴴ, alg::PolarNewton)
+    check_input(right_polar!, A, PWᴴ, alg)
+    P, Wᴴ = PWᴴ
+    if isempty(P)
+        Wᴴ = _right_polarnewton!(A, Wᴴ, P; alg.kwargs...)
+        return P, Wᴴ
+    else
+        Wᴴ = _right_polarnewton!(copy(A), Wᴴ, P; alg.kwargs...)
+        # we still need `A` to compute `P`
+        P = project_hermitian!(mul!(P, A, Wᴴ'))
+        return P, Wᴴ
+    end
+end
+
+# these methods only compute W and destroy A in the process
+function _left_polarnewton!(A::AbstractMatrix, W, P = similar(A, (0, 0)); tol = defaulttol(A), maxiter = 10)
+    m, n = size(A) # we must have m >= n
+    Rᴴinv = isempty(P) ? similar(P, (n, n)) : P # use P as workspace when available
+    if m > n # initial QR
+        Q, R = qr_compact!(A)
+        Rc = view(A, 1:n, 1:n)
+        copy!(Rc, R)
+        Rᴴinv = ldiv!(UpperTriangular(Rc)', one!(Rᴴinv))
+    else # m == n
+        R = A
+        Rc = view(W, 1:n, 1:n)
+        copy!(Rc, R)
+        Rᴴinv = ldiv!(lu!(Rc)', one!(Rᴴinv))
+    end
+    γ = sqrt(norm(Rᴴinv) / norm(R)) # scaling factor
+    rmul!(R, γ)
+    rmul!(Rᴴinv, 1 / γ)
+    R, Rᴴinv = _avgdiff!(R, Rᴴinv)
+    copy!(Rc, R)
+    i = 1
+    conv = norm(Rᴴinv, Inf)
+    while i < maxiter && conv > tol
+        Rᴴinv = ldiv!(lu!(Rc)', one!(Rᴴinv))
+        γ = sqrt(norm(Rᴴinv) / norm(R)) # scaling factor
+        rmul!(R, γ)
+        rmul!(Rᴴinv, 1 / γ)
+        R, Rᴴinv = _avgdiff!(R, Rᴴinv)
+        copy!(Rc, R)
+        conv = norm(Rᴴinv, Inf)
+        i += 1
+    end
+    if conv > tol
+        @warn "`left_polar!` via Newton iteration did not converge within $maxiter iterations (final residual: $conv)"
+    end
+    if m > n
+        return mul!(W, Q, Rc)
+    end
+    return W
+end
+
+function _right_polarnewton!(A::AbstractMatrix, Wᴴ, P = similar(A, (0, 0)); tol = defaulttol(A), maxiter = 10)
+    m, n = size(A) # we must have m <= n
+    Lᴴinv = isempty(P) ? similar(P, (m, m)) : P # use P as workspace when available
+    if m < n # initial QR
+        L, Q = lq_compact!(A)
+        Lc = view(A, 1:m, 1:m)
+        copy!(Lc, L)
+        Lᴴinv = ldiv!(LowerTriangular(Lc)', one!(Lᴴinv))
+    else # m == n
+        L = A
+        Lc = view(Wᴴ, 1:m, 1:m)
+        copy!(Lc, L)
+        Lᴴinv = ldiv!(lu!(Lc)', one!(Lᴴinv))
+    end
+    γ = sqrt(norm(Lᴴinv) / norm(L)) # scaling factor
+    rmul!(L, γ)
+    rmul!(Lᴴinv, 1 / γ)
+    L, Lᴴinv = _avgdiff!(L, Lᴴinv)
+    copy!(Lc, L)
+    i = 1
+    conv = norm(Lᴴinv, Inf)
+    while i < maxiter && conv > tol
+        Lᴴinv = ldiv!(lu!(Lc)', one!(Lᴴinv))
+        γ = sqrt(norm(Lᴴinv) / norm(L)) # scaling factor
+        rmul!(L, γ)
+        rmul!(Lᴴinv, 1 / γ)
+        L, Lᴴinv = _avgdiff!(L, Lᴴinv)
+        copy!(Lc, L)
+        conv = norm(Lᴴinv, Inf)
+        i += 1
+    end
+    if conv > tol
+        @warn "`right_polar!` via Newton iteration did not converge within $maxiter iterations (final residual: $conv)"
+    end
+    if m < n
+        return mul!(Wᴴ, Lc, Q)
+    end
+    return Wᴴ
+end
+
+# in place computation of the average and difference of two arrays
+function _avgdiff!(A::AbstractArray, B::AbstractArray)
+    axes(A) == axes(B) || throw(DimensionMismatch())
+    @simd for I in eachindex(A, B)
+        @inbounds begin
+            a = A[I]
+            b = B[I]
+            A[I] = (a + b) / 2
+            B[I] = b - a
+        end
+    end
+    return A, B
+end

src/implementations/projections.jl

@@ -5,6 +5,8 @@ function copy_input(::typeof(project_hermitian), A::AbstractMatrix)
 end
 copy_input(::typeof(project_antihermitian), A) = copy_input(project_hermitian, A)
 
+copy_input(::typeof(project_isometric), A) = copy_input(left_polar, A)
+
 function check_input(::typeof(project_hermitian!), A::AbstractMatrix, B::AbstractMatrix, ::AbstractAlgorithm)
     LinearAlgebra.checksquare(A)
     n = Base.require_one_based_indexing(A)
@@ -18,6 +20,16 @@ function check_input(::typeof(project_antihermitian!), A::AbstractMatrix, B::Abs
     return nothing
 end
 
+function check_input(::typeof(project_isometric!), A::AbstractMatrix, W::AbstractMatrix, ::AbstractAlgorithm)
+    m, n = size(A)
+    m >= n ||
+        throw(ArgumentError("input matrix needs at least as many rows as columns"))
+    @assert W isa AbstractMatrix
+    @check_size(W, (m, n))
+    @check_scalar(W, A)
+    return nothing
+end
+
 # Outputs
 # -------
 function initialize_output(::typeof(project_hermitian!), A::AbstractMatrix, ::NativeBlocked)
@@ -27,15 +39,26 @@ function initialize_output(::typeof(project_antihermitian!), A::AbstractMatrix,
     return A
 end
 
+function initialize_output(::typeof(project_isometric!), A::AbstractMatrix, ::AbstractAlgorithm)
+    return similar(A)
+end
+
 # Implementation
 # --------------
-function project_hermitian!(A::AbstractMatrix, B, alg::NativeBlocked)
-    check_input(project_hermitian!, A, B, alg)
-    return project_hermitian_native!(A, B, Val(false); alg.kwargs...)
+function project_hermitian!(A::AbstractMatrix, Aₕ, alg::NativeBlocked)
+    check_input(project_hermitian!, A, Aₕ, alg)
+    return project_hermitian_native!(A, Aₕ, Val(false); alg.kwargs...)
 end
-function project_antihermitian!(A::AbstractMatrix, B, alg::NativeBlocked)
-    check_input(project_antihermitian!, A, B, alg)
-    return project_hermitian_native!(A, B, Val(true); alg.kwargs...)
+function project_antihermitian!(A::AbstractMatrix, Aₐ, alg::NativeBlocked)
+    check_input(project_antihermitian!, A, Aₐ, alg)
+    return project_hermitian_native!(A, Aₐ, Val(true); alg.kwargs...)
+end
+
+function project_isometric!(A::AbstractMatrix, W, alg::AbstractAlgorithm)
+    check_input(project_isometric!, A, W, alg)
+    noP = similar(W, (0, 0))
+    W, _ = left_polar!(A, (W, noP), alg)
+    return W
 end
 
 function project_hermitian_native!(A::AbstractMatrix, B::AbstractMatrix, anti::Val; blocksize = 32)