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

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
``````

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

1 Like

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.

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

1 Like

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.

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.

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).

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?

1 Like

I would suggest leaving pivoted Cholesky for a separate PR.