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

**URL:** https://discourse.julialang.org/t/factorize-on-symmetric-types-does-not-try-cholesky/30177
**Category:** Internals & Design
**Tags:** linearalgebra
**Created:** [October 22, 2019, 3:07pm UTC](https://discourse.julialang.org/t/factorize-on-symmetric-types-does-not-try-cholesky/30177 "2019-10-22T15:07:51Z")
**Posts on this page:** 10
**Page:** 1

<div class="post-metadata">

### Author: ![SEA](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sea/32/8956_2.png) [@SEA](https://discourse.julialang.org/u/SEA)
#### Post date: [October 22, 2019, 3:07pm UTC](https://discourse.julialang.org/t/factorize-on-symmetric-types-does-not-try-cholesky/30177/1 "2019-10-22T15:07:51Z")

</div>

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](https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/symmetric.jl)

```julia
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):

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

```

---

<div class="post-metadata">

### Author: ![jlperla](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlperla/32/34332_2.png) [@jlperla](https://discourse.julialang.org/u/jlperla)
#### Post date: [October 22, 2019, 4:44pm UTC](https://discourse.julialang.org/t/factorize-on-symmetric-types-does-not-try-cholesky/30177/2 "2019-10-22T16:44:37Z")

</div>

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

---

<div class="post-metadata">

### Author: ![SEA](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sea/32/8956_2.png) [@SEA](https://discourse.julialang.org/u/SEA)
#### Post date: [October 22, 2019, 7:03pm UTC](https://discourse.julialang.org/t/factorize-on-symmetric-types-does-not-try-cholesky/30177/3 "2019-10-22T19:03:05Z")

</div>

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.

---

<div class="post-metadata">

### Author: ![fredrikekre](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fredrikekre/32/1688_2.png) [@fredrikekre](https://discourse.julialang.org/u/fredrikekre)
#### Post date: [October 22, 2019, 7:16pm UTC](https://discourse.julialang.org/t/factorize-on-symmetric-types-does-not-try-cholesky/30177/4 "2019-10-22T19:16:08Z")

</div>

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

---

<div class="post-metadata">

### Author: ![SEA](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sea/32/8956_2.png) [@SEA](https://discourse.julialang.org/u/SEA)
#### Post date: [October 22, 2019, 7:24pm UTC](https://discourse.julialang.org/t/factorize-on-symmetric-types-does-not-try-cholesky/30177/5 "2019-10-22T19:24:18Z")

</div>

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.

---

<div class="post-metadata">

### Author: ![jlperla](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlperla/32/34332_2.png) [@jlperla](https://discourse.julialang.org/u/jlperla)
#### Post date: [October 23, 2019, 2:17am UTC](https://discourse.julialang.org/t/factorize-on-symmetric-types-does-not-try-cholesky/30177/6 "2019-10-23T02:17:36Z")

</div>

> [@fredrikekre](#):
>
> If you know your matrix is positive definite and want to use Cholesky, why don’t you call `cholesky` then?

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.

---

<div class="post-metadata">

### Author: ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [October 23, 2019, 6:04am UTC](https://discourse.julialang.org/t/factorize-on-symmetric-types-does-not-try-cholesky/30177/7 "2019-10-23T06:04:42Z")

</div>

I think this is a reasonable suggestion. `factorize` for dense matrices [already tries Cholesky](https://github.com/JuliaLang/julia/blob/e5c9a43ba1b3bc67af40513b13481eab27471944/stdlib/LinearAlgebra/src/dense.jl#L1258), 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).

---

<div class="post-metadata">

### Author: ![SEA](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sea/32/8956_2.png) [@SEA](https://discourse.julialang.org/u/SEA)
#### Post date: [October 23, 2019, 9:15pm UTC](https://discourse.julialang.org/t/factorize-on-symmetric-types-does-not-try-cholesky/30177/8 "2019-10-23T21:15:16Z")

</div>

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?

---

<div class="post-metadata">

### Author: ![SEA](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sea/32/8956_2.png) [@SEA](https://discourse.julialang.org/u/SEA)
#### Post date: [October 24, 2019, 12:17am UTC](https://discourse.julialang.org/t/factorize-on-symmetric-types-does-not-try-cholesky/30177/9 "2019-10-24T00:17:20Z")

</div>

[https://github.com/JuliaLang/julia/pull/33652](https://github.com/JuliaLang/julia/pull/33652)

---

<div class="post-metadata">

### Author: ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [October 24, 2019, 5:02am UTC](https://discourse.julialang.org/t/factorize-on-symmetric-types-does-not-try-cholesky/30177/10 "2019-10-24T05:02:09Z")

</div>

I would suggest leaving pivoted Cholesky for a separate PR.
