Factorize() on Symmetric types does not try cholesky()

Internals & Design
1

factorize() on a dense matrix tests several properties to find an efficient factorization. If factorize is called on Symmetric or Hermitian types, it does not try cholesky(), but only considers bunchkaufman():

From the source

function _factorize(A::HermOrSym{T}; check::Bool=true) where T
    TT = typeof(sqrt(oneunit(T)))
    if TT <: BlasFloat
        return bunchkaufman(A; check=check)
    else # fallback
        return lu(A; check=check)
    end
end

Is there any argument against adding (analogous to factorize on dense matrices):

cf = cholesky(A; check = false)
            if cf.info == 0
                return cf
            else
                return bunchkaufman(A; check=check)
           end
2

Or even a SymmetricPositiveDefinite tag type to dispatch to it with factorize?

1 Like
3

I like the idea, though I see the following issue:

Symmetric creates a symmetric view into the original matrix, which is an efficient operation and the symmetry property is enforced.

To enforce positive-definiteness, we would need to project an arbitrary matrix into the space of p.d. matrices, which itself is an O(n^3) operation, similar to the factorization we want to compute. If we did not care about the projection, cholesky() might fail when we call it on the p.d. type. That in turn would require a fall-back, which now makes the solution look very similar to trying cholesky() on a Symmetric type, and falling back to bunchkaufman, and lu.

Another potential type might be diagonally-dominant matrices, which are positive-semi definite, and are easy to check for. Again though, cholesky might fail on them because the current implementation requires strictly p.d. matrices. For that reason, I don’t see a great benefit in adding these types.

4

If you know your matrix is positive definite and want to use Cholesky, why don’t you call cholesky then?

1 Like
5

If the factorize() call occurs in a generic function, which is separate from the piece of code that “knows” about the positive definiteness.

Though my main issue is that factorize() does not try cholesky() on Symmetric or Hermitian types.

6

It is for generic code and consistency. It is nice to be able to exploit the matrix structure without knowing the actual structure. With this one exception, factorize is such a beautiful setup and is the base of \ implementation.

7

I think this is a reasonable suggestion. factorize for dense matrices already tries Cholesky, so the two methods would be consistent. Please consider opening an issue.

That said, in general I agree with @fredrikekre: if you know you can do Cholesky, call cholesky, and similarly for other factorizations. factorize is not type stable (by construction).

8

I experimented locally with a change to make the two factorize methods consistent and could make a pull request soon.

However, I am wondering if we should modify all the cholesky calls (also in the dense factorize) to use the pivoted Cholesky decomposition, which also works on rank deficient matrices and is significantly faster than bunchkaufman in those cases. Opinions?

9

https://github.com/JuliaLang/julia/pull/33652

1 Like
10

I would suggest leaving pivoted Cholesky for a separate PR.