{"id":430,"date":"2016-02-08T16:28:20","date_gmt":"2016-02-08T16:28:20","guid":{"rendered":"http:\/\/drsfenner.org\/blog\/?p=430"},"modified":"2016-02-08T16:29:59","modified_gmt":"2016-02-08T16:29:59","slug":"basic-cholesky-implementation","status":"publish","type":"post","link":"https:\/\/drsfenner.org\/blog\/2016\/02\/basic-cholesky-implementation\/","title":{"rendered":"Basic Cholesky Implementation"},"content":{"rendered":"
\n
\n<\/div>\n
\n
\n

I spent<\/a> a bunch<\/a> of time<\/a> talking about using lower level libraries (LAPACK directly and via LAPACKE or hand wrappers). My next set of posts is going to look at manually implementing these techniques. Mostly, I wanted to spend somem time with the numerical algorithms. I last visited them about 15 years ago in graduate school. There are also several performance lessons along the way (demonstating the limits of a "just compile it" for optimization). For starters, we’ll implement the Cholesky decomposition which is very easy to describe mathematically. <\/p>\n<\/div>\n<\/div>\n<\/div>\n

\n
\n<\/div>\n
\n
\n

Spin Up<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [2]:<\/div>\n
\n
\n
\n
from<\/span> math<\/span> import<\/span> sqrt<\/span>\nimport<\/span> numpy<\/span> as<\/span> np<\/span>\nimport<\/span> numba<\/span>\nimport<\/span> itertools<\/span> as<\/span> it<\/span>\nimport<\/span> functools<\/span> as<\/span> ft<\/span>\nimport<\/span> operator<\/span> as<\/span> op<\/span>\nimport<\/span> timeit<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
Some Testing Helpers<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
I really need to factor these testing helpers out into a testing_utils.py<\/code> module for all these posts. The problem with that is the posts are no longer self-contained for reading and running. *sigh* Nothing’s perfect. Turns out, there was enough "testing" related code that I did<\/em> factor it out into a standalone post<\/a>. In turn, I placed the python code from that post in chol_test_post.py<\/code> and I’m going to import<\/code> it here.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [3]:<\/div>\n
\n
\n
\n
from<\/span> chol_test_post<\/span> import<\/span> (<\/span>time_routines<\/span>,<\/span> same_results<\/span>,<\/span> expand<\/span>,<\/span> compose<\/span>,<\/span> \n                           make_random_vecs_with_cond<\/span>,<\/span> make_spd_opt<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
Basic Cholesky<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
As I said, the mathematical description of Cholesky is pretty simple. For a symmetric, positive-definite matrix (so says Atkinson, Intro. to Numerical Analysis — I forget the corner case of positive-semidefinite), we can factor it: \\(A \\rightarrow LL^T\\)<\/span> where \\(L\\)<\/span> is a lower-triangular matrix. Sometimes, this is called a "matrix square root". Element wise, the diagonal and off diagonal entries of \\(L\\)<\/span> are:<\/p>\n
\\[l_{kk} = \\sqrt{a_{kk} – \\sum^{k-1}_{j=1} l^2_{kj}}\\]<\/span> \\[l_{ik} = \\frac{1}{l_{kk}} \\left( a_{ik} – \\sum^{k-1}_{j=1} l_{ij} l_{kj} \\right)\\]<\/span><\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [4]:<\/div>\n
\n
\n
\n
def<\/span> cholesky0<\/span>(<\/span>A<\/span>):<\/span>\n    """<\/span>\n       Performs a Cholesky decomposition of on symmetric, pos-def A.<\/span>\n       Returns lower-triangular L (full sized, zeroed above diag)<\/span>\n    """<\/span>\n    n<\/span> =<\/span> A<\/span>.<\/span>shape<\/span>[<\/span>0<\/span>]<\/span>\n    L<\/span> =<\/span> np<\/span>.<\/span>zeros_like<\/span>(<\/span>A<\/span>)<\/span>\n\n    # Perform the Cholesky decomposition<\/span>\n    for<\/span> row<\/span> in<\/span> xrange<\/span>(<\/span>n<\/span>):<\/span>\n        for<\/span> col<\/span> in<\/span> xrange<\/span>(<\/span>row<\/span>+<\/span>1<\/span>):<\/span>\n            tmp_sum<\/span> =<\/span> np<\/span>.<\/span>dot<\/span>(<\/span>L<\/span>[<\/span>row<\/span>,:<\/span>col<\/span>],<\/span> L<\/span>[<\/span>col<\/span>,:<\/span>col<\/span>])<\/span>         \n            if<\/span> (<\/span>row<\/span> ==<\/span> col<\/span>):<\/span> # Diagonal elements<\/span>\n                L<\/span>[<\/span>row<\/span>,<\/span> col<\/span>]<\/span> =<\/span> sqrt<\/span>(<\/span>max<\/span>(<\/span>A<\/span>[<\/span>row<\/span>,<\/span>row<\/span>]<\/span> -<\/span> tmp_sum<\/span>,<\/span> 0<\/span>))<\/span>\n            else<\/span>:<\/span>\n                L<\/span>[<\/span>row<\/span>,<\/span>col<\/span>]<\/span> =<\/span> (<\/span>1.0<\/span> \/<\/span> L<\/span>[<\/span>col<\/span>,<\/span>col<\/span>])<\/span> *<\/span> (<\/span>A<\/span>[<\/span>row<\/span>,<\/span>col<\/span>]<\/span> -<\/span> tmp_sum<\/span>)<\/span>\n    return<\/span> L<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
We’re doing as direct an implementation as possible. But, that allows us to call on ContinuumIO’s<\/a> great numba<\/a> project to JIT-compile (just-in-time) our code. It helps (in some ways), as we’ll see below. I’ll throw a note out: I’m not a jit-jedi. More specifically, I’ve only really toyed around with numba<\/code>, so nothing I say here is necessarily current (numba<\/code> moves pretty fast in terms of features) or even best-practice.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [5]:<\/div>\n
\n
\n
\n
@numba.jit<\/span>(<\/span>numba<\/span>.<\/span>double<\/span>[:,:](<\/span>numba<\/span>.<\/span>double<\/span>[:,:]))<\/span>\ndef<\/span> cholesky1<\/span>(<\/span>A<\/span>):<\/span>\n    """<\/span>\n       Performs a Cholesky decomposition of on symmetric, pos-def A.<\/span>\n       Returns lower-triangular L (full sized, zeroed above diag)<\/span>\n    """<\/span>\n    n<\/span> =<\/span> A<\/span>.<\/span>shape<\/span>[<\/span>0<\/span>]<\/span>\n    L<\/span> =<\/span> np<\/span>.<\/span>zeros_like<\/span>(<\/span>A<\/span>)<\/span>\n\n    # Perform the Cholesky decomposition<\/span>\n    for<\/span> row<\/span> in<\/span> xrange<\/span>(<\/span>n<\/span>):<\/span>\n        for<\/span> col<\/span> in<\/span> xrange<\/span>(<\/span>row<\/span>+<\/span>1<\/span>):<\/span>\n            tmp_sum<\/span> =<\/span> np<\/span>.<\/span>dot<\/span>(<\/span>L<\/span>[<\/span>row<\/span>,:<\/span>col<\/span>],<\/span> L<\/span>[<\/span>col<\/span>,:<\/span>col<\/span>])<\/span>         \n            if<\/span> (<\/span>row<\/span> ==<\/span> col<\/span>):<\/span> # Diagonal elements<\/span>\n                L<\/span>[<\/span>row<\/span>,<\/span> col<\/span>]<\/span> =<\/span> sqrt<\/span>(<\/span>max<\/span>(<\/span>A<\/span>[<\/span>row<\/span>,<\/span>row<\/span>]<\/span> -<\/span> tmp_sum<\/span>,<\/span> 0<\/span>))<\/span>\n            else<\/span>:<\/span>\n                L<\/span>[<\/span>row<\/span>,<\/span>col<\/span>]<\/span> =<\/span> (<\/span>1.0<\/span> \/<\/span> L<\/span>[<\/span>col<\/span>,<\/span>col<\/span>])<\/span> *<\/span> (<\/span>A<\/span>[<\/span>row<\/span>,<\/span>col<\/span>]<\/span> -<\/span> tmp_sum<\/span>)<\/span>\n    return<\/span> L<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
I’m not a numba wizard, but it seemed that with my version, I had to eliminate np.dot<\/code> to enable numba<\/code>‘s nopython<\/code>-mode (full-speed-ahead-numba).<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [6]:<\/div>\n
\n
\n
\n
@numba.jit<\/span>(<\/span>numba<\/span>.<\/span>double<\/span>[:,:](<\/span>numba<\/span>.<\/span>double<\/span>[:,:]),<\/span> nopython<\/span>=<\/span>True<\/span>)<\/span>\ndef<\/span> cholesky2<\/span>(<\/span>A<\/span>):<\/span>\n    """<\/span>\n       Performs a Cholesky decomposition of on symmetric, pos-def A.<\/span>\n       Returns lower-triangular L (full sized, zeroed above diag)<\/span>\n    """<\/span>\n    n<\/span> =<\/span> A<\/span>.<\/span>shape<\/span>[<\/span>0<\/span>]<\/span>\n    L<\/span> =<\/span> np<\/span>.<\/span>zeros_like<\/span>(<\/span>A<\/span>)<\/span>\n\n    # Perform the Cholesky decomposition<\/span>\n    for<\/span> row<\/span> in<\/span> xrange<\/span>(<\/span>n<\/span>):<\/span>\n        for<\/span> col<\/span> in<\/span> xrange<\/span>(<\/span>row<\/span>+<\/span>1<\/span>):<\/span>\n\n            # eliminating np.dot allows us to go numba nopython<\/span>\n            # but, the release notes for numba 0.23.0<\/span>\n            # reference:  https:\/\/github.com\/numba\/numba\/pull\/1560<\/span>\n            tmp_sum<\/span> =<\/span> 0.0<\/span>\n            for<\/span> j<\/span> in<\/span> xrange<\/span>(<\/span>col<\/span>):<\/span>\n                tmp_sum<\/span> +=<\/span> L<\/span>[<\/span>row<\/span>,<\/span>j<\/span>]<\/span> *<\/span> L<\/span>[<\/span>col<\/span>,<\/span>j<\/span>]<\/span>\n\n            if<\/span> (<\/span>row<\/span> ==<\/span> col<\/span>):<\/span> \n                # diag elts.<\/span>\n                L<\/span>[<\/span>row<\/span>,<\/span>col<\/span>]<\/span> =<\/span> sqrt<\/span>(<\/span>A<\/span>[<\/span>row<\/span>,<\/span>row<\/span>]<\/span> -<\/span> tmp_sum<\/span>)<\/span>\n            else<\/span>:<\/span>\n                # off-diag elts.<\/span>\n                L<\/span>[<\/span>row<\/span>,<\/span>col<\/span>]<\/span> =<\/span> (<\/span>1.0<\/span> \/<\/span> L<\/span>[<\/span>col<\/span>,<\/span>col<\/span>]<\/span> *<\/span> (<\/span>A<\/span>[<\/span>row<\/span>,<\/span>col<\/span>]<\/span> -<\/span> tmp_sum<\/span>))<\/span>\n    return<\/span> L<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
One of the key rules of numerical linear algebra is that you gain performance when you don’t do work. Ha! Actually, that’s a general rule of life. But, here, it means that (1) when we know the structure of a matrix, then (2) we take advantage of that structure to do less work. Practically speaking, instead of working with a zero-ed starting matrix, we’ll start with raw (arbitrary) memory, and only use values that we know have been set properly. We’ll zero out the super-diagonal pieces separately.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [7]:<\/div>\n
\n
\n
\n
@numba.jit<\/span>(<\/span>numba<\/span>.<\/span>double<\/span>[:,:](<\/span>numba<\/span>.<\/span>double<\/span>[:,:]),<\/span> nopython<\/span>=<\/span>True<\/span>)<\/span>\ndef<\/span> cholesky3<\/span>(<\/span>A<\/span>):<\/span>\n    """<\/span>\n       Performs a Cholesky decomposition of on symmetric, pos-def A.<\/span>\n       Returns lower-triangular L (full sized, zeroed above diag)<\/span>\n    """<\/span>\n    n<\/span> =<\/span> A<\/span>.<\/span>shape<\/span>[<\/span>0<\/span>]<\/span>\n    L<\/span> =<\/span> np<\/span>.<\/span>empty_like<\/span>(<\/span>A<\/span>)<\/span>\n\n    # Perform the Cholesky decomposition<\/span>\n    for<\/span> row<\/span> in<\/span> xrange<\/span>(<\/span>n<\/span>):<\/span>\n        for<\/span> col<\/span> in<\/span> xrange<\/span>(<\/span>row<\/span>+<\/span>1<\/span>):<\/span>\n            tmp_sum<\/span> =<\/span> 0.0<\/span>\n            for<\/span> j<\/span> in<\/span> xrange<\/span>(<\/span>col<\/span>):<\/span>\n                tmp_sum<\/span> +=<\/span> L<\/span>[<\/span>row<\/span>,<\/span>j<\/span>]<\/span> *<\/span> L<\/span>[<\/span>col<\/span>,<\/span>j<\/span>]<\/span>\n            \n            if<\/span> (<\/span>row<\/span> ==<\/span> col<\/span>):<\/span> \n                # diag elts.<\/span>\n                L<\/span>[<\/span>row<\/span>,<\/span>col<\/span>]<\/span> =<\/span> sqrt<\/span>(<\/span>A<\/span>[<\/span>row<\/span>,<\/span>row<\/span>]<\/span> -<\/span> tmp_sum<\/span>)<\/span>\n            else<\/span>:<\/span>\n                # off diag elts.<\/span>\n                L<\/span>[<\/span>row<\/span>,<\/span>col<\/span>]<\/span> =<\/span> (<\/span>1.0<\/span> \/<\/span> L<\/span>[<\/span>col<\/span>,<\/span>col<\/span>]<\/span> *<\/span> (<\/span>A<\/span>[<\/span>row<\/span>,<\/span>col<\/span>]<\/span> -<\/span> tmp_sum<\/span>))<\/span>\n        L<\/span>[<\/span>row<\/span>,<\/span> row<\/span>+<\/span>1<\/span>:]<\/span> =<\/span> 0.0<\/span>\n    return<\/span> L<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
Packed Storage Calculations<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
Turns out, that we really were only scratching the surface of respecting the structure of the matrix. We were still using a \\(n^2\\)<\/span> piece of memory for a triangular matrix when a \\(\\frac{n(n+1)}{2}\\)<\/span> piece of memory will do. Sure, it’s still \\(O(n^2)\\)<\/span>, but savings is savings. And saving 50% is a good thing. We’ll use a long, flat array to store each row of the triangular matrix. To help us out, we’ll definte some helpers to convert from row-and-columns "full" indices to a "straightline" packed index.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [8]:<\/div>\n
\n
\n
\n
#if we were column packing<\/span>\n#@numba.jit(numba.uint(numba.uint, numba.uint, numba.uint), nopython=True)<\/span>\n#def col_packed_lt(r,c,n):<\/span>\n#    return r + (c * (-c + 2*n + 1) * .5)<\/span>\n\n# this would be row-packed upper triangular<\/span>\n#@numba.jit(numba.uint(numba.uint, numba.uint, numba.uint), nopython=True)<\/span>\n#def row_packed(r,c,n):<\/span>\n#    return col_packed(c,r,n)<\/span>\n\n@numba.jit<\/span>(<\/span>numba<\/span>.<\/span>uint<\/span>(<\/span>numba<\/span>.<\/span>uint<\/span>,<\/span> numba<\/span>.<\/span>uint<\/span>,<\/span> numba<\/span>.<\/span>uint<\/span>),<\/span> nopython<\/span>=<\/span>True<\/span>)<\/span>\ndef<\/span> row_packed_lt<\/span>(<\/span>r<\/span>,<\/span>c<\/span>,<\/span>n<\/span>):<\/span>\n    return<\/span> c<\/span> +<\/span> (<\/span>r<\/span> *<\/span> (<\/span>r<\/span>+<\/span>1<\/span>))<\/span> *<\/span> .<\/span>5<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
And, for testing purposes, we may want to expand our packed form back to a full, normal form.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [9]:<\/div>\n
\n
\n
\n
def<\/span> expand_row_packed<\/span>(<\/span>L<\/span>,<\/span> n<\/span>,<\/span> cols<\/span>=<\/span>None<\/span>):<\/span>\n    E<\/span> =<\/span> np<\/span>.<\/span>empty<\/span>((<\/span>n<\/span>,<\/span>cols<\/span>))<\/span> if<\/span> cols<\/span> else<\/span> np<\/span>.<\/span>empty<\/span>((<\/span>n<\/span>,<\/span>n<\/span>))<\/span>\n    E<\/span>[<\/span>np<\/span>.<\/span>tril_indices<\/span>(<\/span>n<\/span>,<\/span>  m<\/span>=<\/span>cols<\/span>)]<\/span> =<\/span> L<\/span>\n    E<\/span>[<\/span>np<\/span>.<\/span>triu_indices<\/span>(<\/span>n<\/span>,<\/span>1<\/span>,<\/span>m<\/span>=<\/span>cols<\/span>)]<\/span> =<\/span> 0.0<\/span>\n    return<\/span> E<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
Here’s a conversion of [0, 1, 2, 3, 4, 5]<\/code> interpreted as a 3×3 lower triangular matrix:<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [10]:<\/div>\n
\n
\n
\n
L<\/span> =<\/span> np<\/span>.<\/span>arange<\/span>(<\/span>6<\/span>)<\/span>\nexpand_row_packed<\/span>(<\/span>L<\/span>,<\/span> 3<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
\n
Out[10]:<\/div>\n
\n
array([[ 0.,  0.,  0.],\n       [ 1.,  2.,  0.],\n       [ 3.,  4.,  5.]])<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
The elements in column 0 of a lower-triangular matrix stored in row-packed form will be at:<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [11]:<\/div>\n
\n
\n
\n
N<\/span> =<\/span> 3<\/span>\nfor<\/span> i<\/span> in<\/span> range<\/span>(<\/span>N<\/span>):<\/span>\n    print<\/span> row_packed_lt<\/span>(<\/span>i<\/span>,<\/span> 0<\/span>,<\/span> N<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
\n
<\/div>\n
\n
0\n1\n3\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
Packed Cholesky Routines<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [12]:<\/div>\n
\n
\n
\n
@numba.jit<\/span>(<\/span>numba<\/span>.<\/span>double<\/span>[:](<\/span>numba<\/span>.<\/span>double<\/span>[:,:]),<\/span> nopython<\/span>=<\/span>True<\/span>)<\/span>\ndef<\/span> cholesky_packed1<\/span>(<\/span>A<\/span>):<\/span>\n    """<\/span>\n       Performs a Cholesky decomposition of on symmetric, pos-def A.<\/span>\n       Returns lower-triangular L packed (row-wise) in a vector<\/span>\n    """<\/span>\n    n<\/span> =<\/span> A<\/span>.<\/span>shape<\/span>[<\/span>0<\/span>]<\/span>\n    numElts<\/span> =<\/span> n<\/span> *<\/span> (<\/span>n<\/span>+<\/span>1<\/span>)<\/span> \/<\/span> 2<\/span>\n    L<\/span> =<\/span> np<\/span>.<\/span>empty<\/span>(<\/span>numElts<\/span>)<\/span>\n    \n    for<\/span> row<\/span> in<\/span> xrange<\/span>(<\/span>n<\/span>):<\/span>\n        for<\/span> col<\/span> in<\/span> xrange<\/span>(<\/span>row<\/span>+<\/span>1<\/span>):<\/span>\n            if<\/span> (<\/span>row<\/span> ==<\/span> col<\/span>):<\/span> \n                # diag elts<\/span>\n                tmp_sum<\/span> =<\/span> 0.0<\/span>\n                base<\/span> =<\/span> row_packed_lt<\/span>(<\/span>col<\/span>,<\/span> 0<\/span>,<\/span> n<\/span>)<\/span>\n                for<\/span> j<\/span> in<\/span> xrange<\/span>(<\/span>col<\/span>):<\/span>\n                    tmp_sum<\/span> +=<\/span> L<\/span>[<\/span>base<\/span> +<\/span> j<\/span>]<\/span> **<\/span> 2<\/span>\n\n                L<\/span>[<\/span>row_packed_lt<\/span>(<\/span>row<\/span>,<\/span>col<\/span>,<\/span>n<\/span>)]<\/span> =<\/span> sqrt<\/span>(<\/span>A<\/span>[<\/span>row<\/span>,<\/span>row<\/span>]<\/span> -<\/span> tmp_sum<\/span>)<\/span>\n            else<\/span>:<\/span>\n                # off diag elts<\/span>\n                tmp_sum<\/span> =<\/span> 0.0<\/span>\n                for<\/span> j<\/span> in<\/span> xrange<\/span>(<\/span>col<\/span>):<\/span>\n                    tmp_sum<\/span> +=<\/span> L<\/span>[<\/span>row_packed_lt<\/span>(<\/span>row<\/span>,<\/span> 0<\/span>,<\/span> n<\/span>)<\/span> +<\/span> j<\/span>]<\/span> *<\/span> L<\/span>[<\/span>row_packed_lt<\/span>(<\/span>col<\/span>,<\/span> 0<\/span>,<\/span> n<\/span>)<\/span> +<\/span> j<\/span>]<\/span>\n                L<\/span>[<\/span>row_packed_lt<\/span>(<\/span>row<\/span>,<\/span>col<\/span>,<\/span>n<\/span>)]<\/span> =<\/span> ((<\/span>1.0<\/span> \/<\/span> L<\/span>[<\/span>row_packed_lt<\/span>(<\/span>col<\/span>,<\/span>col<\/span>,<\/span>n<\/span>)])<\/span> *<\/span> \n                                               (<\/span>A<\/span>[<\/span>row<\/span>][<\/span>col<\/span>]<\/span> -<\/span> tmp_sum<\/span>))<\/span>\n    return<\/span> L<\/span>    \n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
And now, let’s factor out some come computations:<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [13]:<\/div>\n
\n
\n
\n
@numba.jit<\/span>(<\/span>numba<\/span>.<\/span>double<\/span>[:](<\/span>numba<\/span>.<\/span>double<\/span>[:,:]),<\/span> nopython<\/span>=<\/span>True<\/span>)<\/span>\ndef<\/span> cholesky_packed2<\/span>(<\/span>A<\/span>):<\/span>\n    """<\/span>\n       Performs a Cholesky decomposition of on symmetric, pos-def A.<\/span>\n       Returns lower-triangular L packed (row-wise) in a vector<\/span>\n    """<\/span>\n    n<\/span> =<\/span> A<\/span>.<\/span>shape<\/span>[<\/span>0<\/span>]<\/span>\n    numElts<\/span> =<\/span> n<\/span> *<\/span> (<\/span>n<\/span>+<\/span>1<\/span>)<\/span> \/<\/span> 2<\/span>\n    L<\/span> =<\/span> np<\/span>.<\/span>empty<\/span>(<\/span>numElts<\/span>)<\/span>\n    \n    for<\/span> row<\/span> in<\/span> xrange<\/span>(<\/span>n<\/span>):<\/span>\n        for<\/span> col<\/span> in<\/span> xrange<\/span>(<\/span>row<\/span>+<\/span>1<\/span>):<\/span>\n            if<\/span> (<\/span>row<\/span> ==<\/span> col<\/span>):<\/span> # Diagonal elements<\/span>\n                tmp_sum<\/span> =<\/span> 0.0<\/span>\n                base<\/span> =<\/span> row_packed_lt<\/span>(<\/span>col<\/span>,<\/span> 0<\/span>,<\/span> n<\/span>)<\/span>\n                for<\/span> j<\/span> in<\/span> xrange<\/span>(<\/span>col<\/span>):<\/span>\n                    tmp_sum<\/span> +=<\/span> L<\/span>[<\/span>base<\/span> +<\/span> j<\/span>]<\/span> **<\/span> 2<\/span>\n                L<\/span>[<\/span>row_packed_lt<\/span>(<\/span>row<\/span>,<\/span>col<\/span>,<\/span>n<\/span>)]<\/span> =<\/span> sqrt<\/span>(<\/span>A<\/span>[<\/span>row<\/span>,<\/span>row<\/span>]<\/span> -<\/span> tmp_sum<\/span>)<\/span>\n            else<\/span>:<\/span>\n                tmp_sum<\/span> =<\/span> 0.0<\/span>\n                base1<\/span> =<\/span> row_packed_lt<\/span>(<\/span>row<\/span>,<\/span> 0<\/span>,<\/span> n<\/span>)<\/span>\n                base2<\/span> =<\/span> row_packed_lt<\/span>(<\/span>col<\/span>,<\/span> 0<\/span>,<\/span> n<\/span>)<\/span>\n                for<\/span> j<\/span> in<\/span> xrange<\/span>(<\/span>col<\/span>):<\/span>\n                    tmp_sum<\/span> +=<\/span> L<\/span>[<\/span>base1<\/span> +<\/span> j<\/span>]<\/span> *<\/span> L<\/span>[<\/span>base2<\/span> +<\/span> j<\/span>]<\/span>\n                #+= base1 and base2 here *crushed* the timing<\/span>\n                L<\/span>[<\/span>base1<\/span>+<\/span>col<\/span>]<\/span> =<\/span> ((<\/span>1.0<\/span> \/<\/span> L<\/span>[<\/span>base2<\/span>+<\/span>col<\/span>])<\/span> *<\/span> (<\/span>A<\/span>[<\/span>row<\/span>][<\/span>col<\/span>]<\/span> -<\/span> tmp_sum<\/span>))<\/span>\n    return<\/span> L<\/span>    \n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [14]:<\/div>\n
\n
\n
\n
@numba.jit<\/span>(<\/span>numba<\/span>.<\/span>double<\/span>[:](<\/span>numba<\/span>.<\/span>double<\/span>[:,:]),<\/span> nopython<\/span>=<\/span>True<\/span>)<\/span>\ndef<\/span> cholesky_packed3<\/span>(<\/span>A<\/span>):<\/span>\n    """<\/span>\n       Performs a Cholesky decomposition of on symmetric, pos-def A.<\/span>\n       Returns lower-triangular L packed (row-wise) in a vector<\/span>\n    """<\/span>\n    n<\/span> =<\/span> A<\/span>.<\/span>shape<\/span>[<\/span>0<\/span>]<\/span>\n    numElts<\/span> =<\/span> n<\/span> *<\/span> (<\/span>n<\/span>+<\/span>1<\/span>)<\/span> \/<\/span> 2<\/span>\n    L<\/span> =<\/span> np<\/span>.<\/span>empty<\/span>(<\/span>numElts<\/span>)<\/span>\n    \n    for<\/span> row<\/span> in<\/span> xrange<\/span>(<\/span>n<\/span>):<\/span>\n        for<\/span> col<\/span> in<\/span> xrange<\/span>(<\/span>row<\/span>+<\/span>1<\/span>):<\/span>\n            if<\/span> (<\/span>row<\/span> ==<\/span> col<\/span>):<\/span> # Diagonal elements<\/span>\n                base<\/span> =<\/span> row_packed_lt<\/span>(<\/span>col<\/span>,<\/span> 0<\/span>,<\/span> n<\/span>)<\/span>\n                tmp_sum<\/span> =<\/span> np<\/span>.<\/span>sum<\/span>(<\/span>L<\/span>[<\/span>base<\/span>:<\/span>base<\/span>+<\/span>col<\/span>]<\/span>**<\/span>2<\/span>)<\/span>\n                L<\/span>[<\/span>row_packed_lt<\/span>(<\/span>row<\/span>,<\/span>col<\/span>,<\/span>n<\/span>)]<\/span> =<\/span> sqrt<\/span>(<\/span>A<\/span>[<\/span>row<\/span>,<\/span>row<\/span>]<\/span> -<\/span> tmp_sum<\/span>)<\/span>\n            else<\/span>:<\/span>\n                base1<\/span> =<\/span> row_packed_lt<\/span>(<\/span>row<\/span>,<\/span> 0<\/span>,<\/span> n<\/span>)<\/span>\n                base2<\/span> =<\/span> row_packed_lt<\/span>(<\/span>col<\/span>,<\/span> 0<\/span>,<\/span> n<\/span>)<\/span>\n                b1p<\/span> =<\/span> base1<\/span> +<\/span> col<\/span>\n                b2p<\/span> =<\/span> base2<\/span> +<\/span> col<\/span>\n                tmp_sum<\/span> =<\/span> np<\/span>.<\/span>sum<\/span>(<\/span>L<\/span>[<\/span>base1<\/span>:<\/span>b1p<\/span>]<\/span> *<\/span> L<\/span>[<\/span>base2<\/span>:<\/span>b2p<\/span>])<\/span>\n                L<\/span>[<\/span>b1p<\/span>]<\/span> =<\/span> ((<\/span>1.0<\/span> \/<\/span> L<\/span>[<\/span>b2p<\/span>])<\/span> *<\/span> (<\/span>A<\/span>[<\/span>row<\/span>][<\/span>col<\/span>]<\/span> -<\/span> tmp_sum<\/span>))<\/span>\n    return<\/span> L<\/span>    \n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
Testing (Correctness)<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [15]:<\/div>\n
\n
\n
\n
A<\/span> =<\/span> np<\/span>.<\/span>array<\/span>([[<\/span>6<\/span>,<\/span> 3<\/span>,<\/span> 4<\/span>,<\/span> 8<\/span>],<\/span> \n              [<\/span>3<\/span>,<\/span> 6<\/span>,<\/span> 5<\/span>,<\/span> 1<\/span>],<\/span> \n              [<\/span>4<\/span>,<\/span> 5<\/span>,<\/span> 10<\/span>,<\/span> 7<\/span>],<\/span> \n              [<\/span>8<\/span>,<\/span> 1<\/span>,<\/span> 7<\/span>,<\/span> 25<\/span>]],<\/span> dtype<\/span>=<\/span>np<\/span>.<\/span>float64<\/span>)<\/span>\nprint<\/span> "A:"<\/span>\nprint<\/span> A<\/span>\n\nprint<\/span> "L (np.linalg):"<\/span>\nprint<\/span> np<\/span>.<\/span>linalg<\/span>.<\/span>cholesky<\/span>(<\/span>A<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
\n
<\/div>\n
\n
A:\n[[  6.   3.   4.   8.]\n [  3.   6.   5.   1.]\n [  4.   5.  10.   7.]\n [  8.   1.   7.  25.]]\nL (np.linalg):\n[[ 2.44948974  0.          0.          0.        ]\n [ 1.22474487  2.12132034  0.          0.        ]\n [ 1.63299316  1.41421356  2.30940108  0.        ]\n [ 3.26598632 -1.41421356  1.58771324  3.13249102]]\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [16]:<\/div>\n
\n
\n
\n
routines<\/span> =<\/span> [<\/span>np<\/span>.<\/span>linalg<\/span>.<\/span>cholesky<\/span>,<\/span> cholesky0<\/span>,<\/span> cholesky1<\/span>,<\/span> cholesky2<\/span>,<\/span> cholesky3<\/span>]<\/span>\nsame_results<\/span>(<\/span>routines<\/span>,<\/span> (<\/span>A<\/span>,))<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
\n
Out[16]:<\/div>\n
\n
True<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [17]:<\/div>\n
\n
\n
\n
expander<\/span> =<\/span> ft<\/span>.<\/span>partial<\/span>(<\/span>expand_row_packed<\/span>,<\/span> n<\/span>=<\/span>A<\/span>.<\/span>shape<\/span>[<\/span>0<\/span>])<\/span>\nroutines<\/span> =<\/span> [<\/span>np<\/span>.<\/span>linalg<\/span>.<\/span>cholesky<\/span>]<\/span> +<\/span> [<\/span>compose<\/span>(<\/span>expander<\/span>,<\/span> r<\/span>)<\/span>for<\/span> r<\/span> in<\/span> (<\/span>cholesky_packed1<\/span>,<\/span> \n                                                                cholesky_packed2<\/span>,<\/span> \n                                                                cholesky_packed3<\/span>)]<\/span>\nsame_results<\/span>(<\/span>routines<\/span>,<\/span> (<\/span>A<\/span>,))<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
\n
Out[17]:<\/div>\n
\n
True<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
Testing (Performance)<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [18]:<\/div>\n
\n
\n
\n
A<\/span> =<\/span> np<\/span>.<\/span>array<\/span>([[<\/span>6<\/span>,<\/span> 3<\/span>,<\/span> 4<\/span>,<\/span> 8<\/span>],<\/span> \n              [<\/span>3<\/span>,<\/span> 6<\/span>,<\/span> 5<\/span>,<\/span> 1<\/span>],<\/span> \n              [<\/span>4<\/span>,<\/span> 5<\/span>,<\/span> 10<\/span>,<\/span> 7<\/span>],<\/span> \n              [<\/span>8<\/span>,<\/span> 1<\/span>,<\/span> 7<\/span>,<\/span> 25<\/span>]],<\/span> dtype<\/span>=<\/span>np<\/span>.<\/span>float64<\/span>)<\/span>\nprint<\/span> "A:"<\/span>\nprint<\/span> A<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
\n
<\/div>\n
\n
A:\n[[  6.   3.   4.   8.]\n [  3.   6.   5.   1.]\n [  4.   5.  10.   7.]\n [  8.   1.   7.  25.]]\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [19]:<\/div>\n
\n
\n
\n
time_routines<\/span>([<\/span>np<\/span>.<\/span>linalg<\/span>.<\/span>cholesky<\/span>],<\/span> (<\/span>A<\/span>,),<\/span> number<\/span>=<\/span>10000<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
\n
<\/div>\n
\n
cholesky : 0.1380\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [20]:<\/div>\n
\n
\n
\n
time_routines<\/span>([<\/span>cholesky0<\/span>,<\/span> cholesky1<\/span>,<\/span> cholesky2<\/span>,<\/span> cholesky3<\/span>],<\/span> (<\/span>A<\/span>,),<\/span> number<\/span>=<\/span>10000<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
\n
<\/div>\n
\n
cholesky2 : 0.0064\ncholesky3 : 0.0065\ncholesky0 : 0.2947\ncholesky1 : 0.3358\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [21]:<\/div>\n
\n
\n
\n
time_routines<\/span>([<\/span>cholesky_packed1<\/span>,<\/span> cholesky_packed2<\/span>,<\/span> cholesky_packed3<\/span>],<\/span> (<\/span>A<\/span>,),<\/span> number<\/span>=<\/span>10000<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
\n
<\/div>\n
\n
cholesky_packed1 : 0.0073\ncholesky_packed2 : 0.0076\ncholesky_packed3 : 0.0179\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [22]:<\/div>\n
\n
\n
\n
_all<\/span> =<\/span> [<\/span>np<\/span>.<\/span>linalg<\/span>.<\/span>cholesky<\/span>,<\/span> \n        cholesky0<\/span>,<\/span> cholesky1<\/span>,<\/span> cholesky2<\/span>,<\/span> cholesky3<\/span>,<\/span>\n        cholesky_packed1<\/span>,<\/span> cholesky_packed2<\/span>,<\/span> cholesky_packed3<\/span>]<\/span>\ntime_routines<\/span>(<\/span>_all<\/span>,<\/span> (<\/span>A<\/span>,),<\/span> number<\/span>=<\/span>1000<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
\n
<\/div>\n
\n
       cholesky2 : 0.0006\n       cholesky3 : 0.0006\ncholesky_packed1 : 0.0007\ncholesky_packed2 : 0.0007\ncholesky_packed3 : 0.0017\n        cholesky : 0.0138\n       cholesky0 : 0.0296\n       cholesky1 : 0.0336\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
Bigger Fish<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
What happens when we scale up to bigger input matrices? Well, since we need SPD (symmetric, positive-definite) matrices, we have to construct them. Doing that in a naive way has a major issue: we can get a horribly<\/em> conditioned matrix (and Cholesky is sensitive to the square<\/em> of the condition number). We can improve that by adding (relatively) large elements on the diagonal to push the matrix towards being "diagonally dominant".<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [23]:<\/div>\n
\n
\n
\n
n<\/span> =<\/span> np<\/span>.<\/span>random<\/span>.<\/span>random_integers<\/span>(<\/span>1<\/span>,<\/span>10<\/span>,(<\/span>5<\/span>,<\/span>1<\/span>))<\/span>\nouter_n<\/span> =<\/span> np<\/span>.<\/span>outer<\/span>(<\/span>n<\/span>,<\/span>n<\/span>)<\/span>\nprint<\/span> outer_n<\/span>\nprint<\/span> np<\/span>.<\/span>real_if_close<\/span>(<\/span>np<\/span>.<\/span>linalg<\/span>.<\/span>eig<\/span>(<\/span>outer_n<\/span>)[<\/span>0<\/span>])<\/span> # remove small imaginary parts<\/span>\nprint<\/span> np<\/span>.<\/span>linalg<\/span>.<\/span>cond<\/span>(<\/span>outer_n<\/span>)<\/span> # can go to float("inf")<\/span>\n\n# condition normalizer; not-so-regular Tikhonov regularization<\/span>\n# pushes towards diagonally dominant<\/span>\nouter_n<\/span> =<\/span> np<\/span>.<\/span>outer<\/span>(<\/span>n<\/span>,<\/span>n<\/span>)<\/span> +<\/span> 5<\/span> *<\/span> np<\/span>.<\/span>eye<\/span>(<\/span>5<\/span>)<\/span> \nprint<\/span> outer_n<\/span>\nprint<\/span> np<\/span>.<\/span>real_if_close<\/span>(<\/span>np<\/span>.<\/span>linalg<\/span>.<\/span>eig<\/span>(<\/span>outer_n<\/span>)[<\/span>0<\/span>])<\/span>\nprint<\/span> np<\/span>.<\/span>linalg<\/span>.<\/span>cond<\/span>(<\/span>outer_n<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
\n
<\/div>\n
\n
[[ 49  35  70  14  70]\n [ 35  25  50  10  50]\n [ 70  50 100  20 100]\n [ 14  10  20   4  20]\n [ 70  50 100  20 100]]\n[  0.00000000e+00   2.78000000e+02   1.84889275e-32   8.20146785e-17\n   0.00000000e+00]\n3.80552400456e+49\n[[  54.   35.   70.   14.   70.]\n [  35.   30.   50.   10.   50.]\n [  70.   50.  105.   20.  100.]\n [  14.   10.   20.    9.   20.]\n [  70.   50.  100.   20.  105.]]\n[   5.  283.    5.    5.    5.]\n56.6\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
But, we’re going to make use of the make_spd_opt<\/code> routine we defined in our testing helpers. Here’s what it looks like (3 is the length of the diagonal, 10 is the condition number of the result).<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [24]:<\/div>\n
\n
\n
\n
A<\/span> =<\/span> make_spd_opt<\/span>(<\/span>3<\/span>,<\/span> 10<\/span>)<\/span>\nprint<\/span> A<\/span>\nprint<\/span> np<\/span>.<\/span>linalg<\/span>.<\/span>eig<\/span>(<\/span>A<\/span>)[<\/span>0<\/span>]<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
\n
<\/div>\n
\n
[[ 3.02040816  2.75510204  2.20408163]\n [ 2.75510204  8.43877551  0.55102041]\n [ 2.20408163  0.55102041  5.04081633]]\n[  1.   10.    5.5]\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [24]:<\/div>\n
\n
\n
\n
A<\/span> =<\/span> make_spd_opt<\/span>(<\/span>2<\/span>**<\/span>10<\/span>,<\/span> 10<\/span>)<\/span>\n_all<\/span> =<\/span> [<\/span>np<\/span>.<\/span>linalg<\/span>.<\/span>cholesky<\/span>,<\/span> \n        cholesky0<\/span>,<\/span> cholesky1<\/span>,<\/span> cholesky2<\/span>,<\/span> cholesky3<\/span>,<\/span>\n        cholesky_packed1<\/span>,<\/span> cholesky_packed2<\/span>,<\/span> cholesky_packed3<\/span>]<\/span>\ntime_routines<\/span>(<\/span>_all<\/span>,<\/span> (<\/span>A<\/span>,),<\/span> number<\/span>=<\/span>1<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
\n
<\/div>\n
\n
        cholesky : 0.0240\n       cholesky3 : 0.1463\n       cholesky2 : 0.1465\ncholesky_packed2 : 0.2107\ncholesky_packed1 : 0.2144\ncholesky_packed3 : 0.3146\n       cholesky0 : 1.1020\n       cholesky1 : 1.3282\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [25]:<\/div>\n
\n
\n
\n
A<\/span> =<\/span> make_spd_opt<\/span>(<\/span>2<\/span>**<\/span>12<\/span>,<\/span> 10<\/span>)<\/span>\ndrop_pokey<\/span> =<\/span> [<\/span>np<\/span>.<\/span>linalg<\/span>.<\/span>cholesky<\/span>,<\/span>\n              cholesky2<\/span>,<\/span> cholesky3<\/span>]<\/span>\ntime_routines<\/span>(<\/span>drop_pokey<\/span>,<\/span> (<\/span>A<\/span>,),<\/span> number<\/span>=<\/span>1<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
\n
<\/div>\n
\n
 cholesky : 0.8519\ncholesky3 : 10.5039\ncholesky2 : 10.5152\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [26]:<\/div>\n
\n
\n
\n
time_routines<\/span>([<\/span>cholesky_packed2<\/span>],<\/span> [<\/span>A<\/span>],<\/span> number<\/span>=<\/span>1<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
\n
<\/div>\n
\n
cholesky_packed2 : 14.2552\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [27]:<\/div>\n
\n
\n
\n
A<\/span> =<\/span> make_spd_opt<\/span>(<\/span>2<\/span>**<\/span>13<\/span>,<\/span> 10<\/span>)<\/span>\nshowdown<\/span> =<\/span> [<\/span>np<\/span>.<\/span>linalg<\/span>.<\/span>cholesky<\/span>,<\/span> cholesky2<\/span>,<\/span> cholesky_packed2<\/span>]<\/span>\ntime_routines<\/span>(<\/span>showdown<\/span>,<\/span> [<\/span>A<\/span>],<\/span> number<\/span>=<\/span>1<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
\n
<\/div>\n
\n
        cholesky : 5.3662\n       cholesky2 : 82.1272\ncholesky_packed2 : 113.4598\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
When the data gets bigger, we can’t compete. Why not? Because while we are increasing the efficiency of the steps of computation, we’re not increasing the efficiency of data movement. And to solve that, generally we switch to "block-algoirthms" that rely on moving chunks (submatrices) of data around, instead of just single rows and\/or columns. Here’s the fundamental issue laid out nicely<\/a> and here’s how deep that rabbit hole goes<\/a>.<\/p>\n
If you want to implement a block-Cholesky, see:<\/p>\n