Representing a symmetric matrix factorization

Specific Domains Numerics
,
1

This question is about implementing a matrix factorization with some operations in Julia, ideally relying on existing optimized implementations, in a generic way do that I don’t have to commit to a particular implementation (dense/sparse, etc) and have it work with AD.

What I need

Suppose J is an n \times m matrix (n \ge m possible). We form B = J'J, and want to support the following operations:

  1. products x'Bx for any conformable x vectors,

  2. testing if B is positive definite (it is always positive semidefinite by construction),

  3. B^{-1} y for conformable y vectors.

  4. Forming the actual matrix B (for input into hvcat).

Implementation questions

Ideally these should be fast, work for dense and sparse B (J is a Jacobian obtained from AD, may be sparse). The code should be generic and work with AD (ForwardDiff.Dual, Enzyme.jl, etc).

I thought of making a J=QR decomposition and storing R. This is not unlike a Cholesky decomposition, except that LinearAlgebra.Cholesky requires positive definiteness if I understand correctly, and dot(x, B, x) does not specialize for Cholesky (but the generic fallback may not be so bad). And I did not check AD’ability of the sparse representation.

I am wondering if there is a better way to do it though — can one arrive at an LDL' representation from J without forming B? Or is there an better factorization?

2

Sounds like you might want ldlt? With a wrapper that’s a lot like PDMat except with mathematics appropriate for Bunch-Kauffman?

(it is always positive semidefinite by construction)

Though if this is true, why not just use PDMat?

1 Like
3

Possibly, but recall that the input is J, and I am not sure if I want to form B = J'J. Is that the only way to get B = LDL'?

Because it forms both the decomposition and the full matrix unconditionally. I rarely need the full matrix.

4

Ah, I didn’t understand that J was supplied. I agree with computing the QR decomposition and storing R. Then

struct SymProduct{T,RT<:AbstractMatrix{T}} <: Factorization{T}
    R::RT
end

Base.size(B::SymProduct) = (m = size(B.R, 1); return (m, m))
Base.:*(B::SymProduct, x) = (y = B.R*x; return y'*y)
function LinearAlgebra.isposdef(B::SymProduct)
    (; R = B)
    m, n = size(R)
    m > n && return false  # fast path
    # check for zeros in R
end
...

If you want SymProduct{T,RT<:AbstractMatrix{T}} <: AbstractMatrix{T} instead, then you need to store B because B[i, j] needs to be fast for it to truly be considered an AbstractMatrix (that’s basically the defining difference between Factorization and AbstractMatrix).

It would be useful to have this as a standalone package, as this does come up in multiple applications. JuliaLinearAlgebra could host, ping me here or on slack if you need permissions. Or perhaps it could be added to PDMats.jl.

1 Like