{"id":414,"date":"2016-02-05T16:06:07","date_gmt":"2016-02-05T16:06:07","guid":{"rendered":"http:\/\/drsfenner.org\/blog\/?p=414"},"modified":"2016-02-05T16:33:45","modified_gmt":"2016-02-05T16:33:45","slug":"testing-a-cholesky-implementation","status":"publish","type":"post","link":"https:\/\/drsfenner.org\/blog\/2016\/02\/testing-a-cholesky-implementation\/","title":{"rendered":"Testing a Cholesky Implementation"},"content":{"rendered":"
\n
\n<\/div>\n
\n
\n

In working through the Cholesky, QR, and SVD methods for solving the normal regression equations, I got curious about the underlying implementations of each. I did QR and Cholesky in grad school (a great Numerical Analysis sequence). I never got around to SVD (I think we did some variation of tridiagonalization). Anyway, Cholesky is so simple that I think testing Cholesky<\/em> generated more interesting points for me. There were two big issues and the "testing" post is (I think) longer than the "implementation" post.<\/p>\n

Executing the tests brought up an interesting software engineering (code reuse) issue and creating good test inputs brought up an intersting mathematical issue. Read on! <\/p>\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
In [1]:<\/div>\n
\n
\n
\n
import<\/span> numpy<\/span> as<\/span> np<\/span>\nimport<\/span> itertools<\/span> as<\/span> it<\/span>\nimport<\/span> functools<\/span> as<\/span> ft<\/span>\nimport<\/span> timeit<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [2]:<\/div>\n
\n
\n
\n
# from python itertools docs:<\/span>\ndef<\/span> pairwise<\/span>(<\/span>iterable<\/span>):<\/span>\n    "s -> (s0,s1), (s1,s2), (s2, s3), ..."<\/span>\n    a<\/span>,<\/span> b<\/span> =<\/span> it<\/span>.<\/span>tee<\/span>(<\/span>iterable<\/span>)<\/span>\n    next<\/span>(<\/span>b<\/span>,<\/span> None<\/span>)<\/span>\n    return<\/span> it<\/span>.<\/span>izip<\/span>(<\/span>a<\/span>,<\/span> b<\/span>)<\/span>\n\ndef<\/span> time_routines<\/span>(<\/span>routines<\/span>,<\/span> inputs<\/span>,<\/span> number<\/span>=<\/span>10<\/span>):<\/span>\n    name_len<\/span> =<\/span> max<\/span>(<\/span>len<\/span>(<\/span>r<\/span>.<\/span>__name__<\/span>)<\/span> for<\/span> r<\/span> in<\/span> routines<\/span>)<\/span>\n    line_template<\/span> =<\/span> "%"<\/span> +<\/span> str<\/span>(<\/span>name_len<\/span>)<\/span> +<\/span> "s : <\/span>%6.4f<\/span>"<\/span>\n\n    times<\/span> =<\/span> {}<\/span>\n    for<\/span> r<\/span> in<\/span> routines<\/span>:<\/span>\n        call<\/span> =<\/span> ft<\/span>.<\/span>partial<\/span>(<\/span>r<\/span>,<\/span> *<\/span>inputs<\/span>)<\/span>\n        times<\/span>[<\/span>r<\/span>.<\/span>__name__<\/span>]<\/span> =<\/span> min<\/span>(<\/span>timeit<\/span>.<\/span>Timer<\/span>(<\/span>call<\/span>)<\/span>.<\/span>repeat<\/span>(<\/span>repeat<\/span>=<\/span>3<\/span>,<\/span> number<\/span>=<\/span>number<\/span>))<\/span>\n    \n    print<\/span> "<\/span>\\n<\/span>"<\/span>.<\/span>join<\/span>(<\/span>line_template<\/span> %<\/span> (<\/span>k<\/span>,<\/span> times<\/span>[<\/span>k<\/span>])<\/span> for<\/span> k<\/span> in<\/span> sorted<\/span>(<\/span>times<\/span>,<\/span> key<\/span>=<\/span>times<\/span>.<\/span>get<\/span>))<\/span>\n\ndef<\/span> same_results<\/span>(<\/span>routines<\/span>,<\/span> inputs<\/span>):<\/span>\n    try<\/span>:<\/span>\n        results<\/span> =<\/span> {<\/span>r<\/span>.<\/span>__name__<\/span>:<\/span>r<\/span>(<\/span>*<\/span>inputs<\/span>)<\/span> for<\/span> r<\/span> in<\/span> routines<\/span>}<\/span>\n        return<\/span> all<\/span>(<\/span>np<\/span>.<\/span>allclose<\/span>(<\/span>r1<\/span>,<\/span>r2<\/span>)<\/span> for<\/span> r1<\/span>,<\/span>r2<\/span> in<\/span> pairwise<\/span>(<\/span>results<\/span>.<\/span>values<\/span>()))<\/span>\n    except<\/span>:<\/span> # quick-and-dirty<\/span>\n        return<\/span> False<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
So, with these helper routines, we can take several different implementations and compare them two ways: (1) do they all give the same answer (same_results<\/code>) and (2) what are their run times (time_routines<\/code>). Note that if the inputs are mutated in-place, we will be solving different problems in both of these. We’ll address that in a second. same_results<\/code> plays a little fast-and-loose because, in general, floating-point closeness is not transitive<\/em> — but closeness treats us far<\/em> better than floating-point equality. However, we’ll be content enough check that if close(a,b)<\/code> and close(b,c)<\/code>, then we’re "all close enough".<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
The Reuse Issue<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
Now, there are two problems that come up when we use this testing code. Since our implementations are not mutating, we don’t have to worry about that. However, we do have a similar issue. We have two different data structures<\/em> being returned. Our "naive" method returns a full NumPy<\/code> array; we also have "packed" versions that return a flat array that is a packed version of the triangular matrix being returned. So, to check for equality of results we would either (1) define a custom comparison routine or (2) provide a converter that will massage the packed versions into a full version. Here’s code that shows the issue:<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [3]:<\/div>\n
\n
\n
\n
def<\/span> makeFullLowerTri<\/span>(<\/span>n<\/span>):<\/span>\n    X<\/span> =<\/span> np<\/span>.<\/span>empty<\/span>((<\/span>n<\/span>,<\/span>n<\/span>))<\/span>\n    X<\/span>[<\/span>np<\/span>.<\/span>triu_indices<\/span>(<\/span>n<\/span>)]<\/span> =<\/span> 0.0<\/span>\n    X<\/span>[<\/span>np<\/span>.<\/span>tril_indices<\/span>(<\/span>n<\/span>)]<\/span> =<\/span> np<\/span>.<\/span>arange<\/span>((<\/span>n<\/span>*<\/span>(<\/span>n<\/span>+<\/span>1<\/span>))<\/span>\/<\/span>2.0<\/span>)<\/span>\n    return<\/span> X<\/span>\nprint<\/span> makeFullLowerTri<\/span>(<\/span>3<\/span>)<\/span>\n\ndef<\/span> makePackedLowerTri<\/span>(<\/span>n<\/span>):<\/span>\n    return<\/span> np<\/span>.<\/span>arange<\/span>((<\/span>n<\/span>*<\/span>(<\/span>n<\/span>+<\/span>1<\/span>))<\/span>\/<\/span>2.0<\/span>)<\/span>\nprint<\/span> makePackedLowerTri<\/span>(<\/span>3<\/span>)<\/span>\n\nprint<\/span> same_results<\/span>([<\/span>makeFullLowerTri<\/span>,<\/span> makePackedLowerTri<\/span>],<\/span> (<\/span>3<\/span>,))<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
\n
<\/div>\n
\n
[[ 0.  0.  0.]\n [ 1.  2.  0.]\n [ 3.  4.  5.]]\n[ 0.  1.  2.  3.  4.  5.]\nFalse\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
BTW, you might see the 2.0<\/code> in the np.arange<\/code> expression. If arange<\/code> takes a float<\/code> as its input, it produces an array of np.float64<\/code> (more technically, np.float_<\/code> … "machine floating point size"). If the input is an int<\/code>, the dtype<\/code> is also np.int_<\/code>. Anyway, it’s just a trick to get the dtype<\/code> squared away nicely.<\/p>\n
So, here’s a quick "expander" which turns packed triangular into full.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [4]:<\/div>\n
\n
\n
\n
def<\/span> expand<\/span>(<\/span>packedTri<\/span>,<\/span> n<\/span>):<\/span>\n    X<\/span> =<\/span> np<\/span>.<\/span>empty<\/span>((<\/span>n<\/span>,<\/span>n<\/span>))<\/span>\n    X<\/span>[<\/span>np<\/span>.<\/span>triu_indices<\/span>(<\/span>n<\/span>)]<\/span> =<\/span> 0.0<\/span>\n    X<\/span>[<\/span>np<\/span>.<\/span>tril_indices<\/span>(<\/span>n<\/span>)]<\/span> =<\/span> packedTri<\/span>\n    return<\/span> X<\/span>\nprint<\/span> np<\/span>.<\/span>allclose<\/span>(<\/span>expand<\/span>(<\/span>makePackedLowerTri<\/span>(<\/span>3<\/span>),<\/span> 3<\/span>),<\/span> makeFullLowerTri<\/span>(<\/span>3<\/span>))<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
\n
<\/div>\n
\n
True\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
But, our same_results<\/code> routine take one function to apply (not nested applications). That process of nesting applications of a function has a name: function composition<\/em>. And, you might remember it from studying such wonders as the chain rule of derivatives. But, unless you are into functional programming and languages, you might not have played around with it much from a software engineering perspective. What we need is (1) a single function that is the result of (2) applying expand<\/code> on makePackedLowerTri<\/code>. Fortunately, Python makes it relatively easy to create such a critter.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [5]:<\/div>\n
\n
\n
\n
def<\/span> helper<\/span>(<\/span>main<\/span>,<\/span> post<\/span>):<\/span>\n    ''' help out by creating\/returning a process that calls *main* to do the "work"<\/span>\n        and then *post* to do the "post-processing"'''<\/span>\n    def<\/span> inner<\/span>(<\/span>A<\/span>):<\/span>\n        ''' compute the composition of two functions on A'''<\/span>\n        return<\/span> post<\/span>(<\/span>main<\/span>(<\/span>A<\/span>))<\/span>\n    return<\/span> inner<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
Now, the next issue is that we don’t want to modify main<\/code> – which will be our makePackedLowerTri<\/code>. In particular, it returns one array and we want to keep it that way. That means that post<\/code> (expand<\/code>) is going to be missing an argument – the n<\/code> value. If only there were a way to "partially" apply a function to some arguments and leave the rest "to be filled in later". But wait! There is:<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [6]:<\/div>\n
\n
\n
\n
# above, we did:  import functools as ft <\/span>\nexpand3<\/span> =<\/span> ft<\/span>.<\/span>partial<\/span>(<\/span>expand<\/span>,<\/span> n<\/span>=<\/span>3<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
expand3<\/code> is a function of one<\/em> argument that has its n<\/code> argument fixed at 3<\/code>. It is as if we had typed (without the typing and with many more resuse possibilities):<\/p>\n
\n
def<\/span> expand3(packedTri):\n    X =<\/span> np.empty((3<\/span>,3<\/span>))\n    X[np.triu_indices(3<\/span>)] =<\/span> 0.0<\/span>\n    X[np.tril_indices(3<\/span>)] =<\/span> packedTri\n    return<\/span> X<\/code><\/pre>\n<\/div>\n
A couple of quick notes. (1) We’ll see that we don’t have<\/em> to lock down the value of n<\/code>. We could compute it inside of expand<\/code> or we could hack on helper<\/code> or makePackedLowerTri<\/code>. (2) Adding compositions is a problem, if we want to directly compare the work<\/em> (memory, time, etc.) of the base functions. We’ll ignore that for now because we are comparing correctness.<\/p>\n
Here’s how the pieces work together. makeFullLowerTri<\/code> is<\/em> a "direct" function that we give an arg to and get our usable result. For the packed version, we do a bit more work. We call helper<\/code> and it makes a function for us<\/em>. That function that it creates is the one we want to call. So, we have two calls: (1) we call helper<\/code> and (2) we call the returned value (from helper<\/code>). The second call gives us a comparable result. Here they are:<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [7]:<\/div>\n
\n
\n
\n
# simple case<\/span>\nprint<\/span> makeFullLowerTri<\/span>(<\/span>3<\/span>)<\/span>\n\n# complex case<\/span>\nr<\/span> =<\/span> helper<\/span>(<\/span>makePackedLowerTri<\/span>,<\/span> expand3<\/span>)<\/span>\nprint<\/span> r<\/span>(<\/span>3<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
\n
<\/div>\n
\n
[[ 0.  0.  0.]\n [ 1.  2.  0.]\n [ 3.  4.  5.]]\n[[ 0.  0.  0.]\n [ 1.  2.  0.]\n [ 3.  4.  5.]]\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
And here’s how we use it in same_results<\/code>:<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [8]:<\/div>\n
\n
\n
\n
routines<\/span> =<\/span> [<\/span>makeFullLowerTri<\/span>,<\/span> helper<\/span>(<\/span>makePackedLowerTri<\/span>,<\/span> expand3<\/span>)]<\/span>\nprint<\/span> same_results<\/span>(<\/span>routines<\/span>,<\/span> (<\/span>3<\/span>,))<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
\n
<\/div>\n
\n
True\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
Generalizing<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
The process of composing functions – and composing more than two functions – seems like a nicely repeatable process. Can we make the above more general than just helper<\/code>? The answer is "yes!". There used to be a python module (third-party) called functional<\/code> that had a compose<\/code> routine. However, I can’t find it anymore (it is reference in the python 3.1 functional HOWTO document). But, a little hacking and Googling found some resources and resulted in the following code:<\/p>\n