# Sparse Cholesky of Gram Matrix (SuiteSparse?)

**URL:** <https://discourse.julialang.org/t/sparse-cholesky-of-gram-matrix-suitesparse/37303>\
**Category:** Numerics\
**Tags:** question, linearalgebra\
**Created:** [April 9, 2020, 10:48pm UTC](https://discourse.julialang.org/t/sparse-cholesky-of-gram-matrix-suitesparse/37303 "2020-04-09T22:48:27Z")\
**Posts on this page:** 6\
**Page:** 2

<div class="post-metadata">

**Author:** ![guille](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/guille/32/7505_2.png) [@guille](https://discourse.julialang.org/u/guille)\
**Post date:** [April 11, 2020, 12:14am UTC](https://discourse.julialang.org/t/sparse-cholesky-of-gram-matrix-suitesparse/37303/21 "2020-04-11T00:14:59Z")

</div>

I’m curious as well. I discarded this “expanded system” solution some time ago since the matrices are quite large, but I can definitely give it a try to see how quickly the problem is solved!

> [@dpo](#):
>
> Note that the HSL codes are typically only free for academia. MUMPS is open source and can use MPI if your systems are that large. If you have access to it, I find that HSL\_MA57 is often the fastest. You can store the factorization and reuse it to solve for multiple right-hand sides.

I actually haven’t tried HSL at all for this (I actually forgot it existed…), so I will do that at some point in the near future. Thanks for the heads up!

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 11, 2020, 2:09am UTC](https://discourse.julialang.org/t/sparse-cholesky-of-gram-matrix-suitesparse/37303/22 "2020-04-11T02:09:54Z")

</div>

> [@guille](#):
>
> actually computing AA^T and taking the Cholesky factorization of that seems to be roughly the same performance

You’re right. And if the matrix is very non-square, then forming AA^T first is much faster:

```julia
A = sprand(1000, 10^6, 1e-6);
A += SparseMatrixCSC(I, size(A)...)
cholgram2(A) = cholesky(A*A')
@btime cholgram($A);
@btime cholgram2($A);

```

gives

```julia
21.030 ms (38 allocations: 23.15 MiB)
1.880 ms (47 allocations: 339.10 KiB)

```

whereas a square matrix

```julia
A = sprand(10^6, 10^6, 1e-9);
A += SparseMatrixCSC(I, size(A)...)
@btime cholgram($A);
@btime cholgram2($A);

```

gives

```julia
207.157 ms (38 allocations: 213.74 MiB)
200.582 ms (52 allocations: 240.73 MiB)

```

I wonder what the point of this feature of CHOLMOD is, if it isn’t faster? Maybe it is more accurate if A is badly conditioned, since it avoids squaring the condition number (I guess)?

---

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [April 11, 2020, 12:19pm UTC](https://discourse.julialang.org/t/sparse-cholesky-of-gram-matrix-suitesparse/37303/23 "2020-04-11T12:19:04Z")

</div>

no need for an interpolation mark (`$`) with `@btime` ?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 11, 2020, 12:37pm UTC](https://discourse.julialang.org/t/sparse-cholesky-of-gram-matrix-suitesparse/37303/24 "2020-04-11T12:37:09Z")

</div>

> [@LaurentPlagne](#):
>
> no need for an interpolation mark ( `$` ) with `@btime` ?

Interpolating global arguments to functions is [always a good idea](https://github.com/JuliaCI/BenchmarkTools.jl/blob/master/doc/manual.md#interpolating-values-into-benchmark-expressions) with `@btime`, so that there is no dynamic-dispatch overhead.

---

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [April 11, 2020, 12:49pm UTC](https://discourse.julialang.org/t/sparse-cholesky-of-gram-matrix-suitesparse/37303/25 "2020-04-11T12:49:53Z")

</div>

That is why I wonder about the impact of missing $ marks in @guille measurements.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 11, 2020, 3:07pm UTC](https://discourse.julialang.org/t/sparse-cholesky-of-gram-matrix-suitesparse/37303/26 "2020-04-11T15:07:47Z")

</div>

These matrices are big enough that the overhead of dynamic dispatch is probably negligible, so `$` interpolation is just a good habit.

[Previous page](https://discourse.julialang.org/t/sparse-cholesky-of-gram-matrix-suitesparse/37303.md?page=1)