{"id":407,"date":"2016-01-27T16:14:59","date_gmt":"2016-01-27T16:14:59","guid":{"rendered":"http:\/\/drsfenner.org\/blog\/?p=407"},"modified":"2016-02-05T16:30:49","modified_gmt":"2016-02-05T16:30:49","slug":"into-the-weeds-iii-hacking-lapack-calls-to-save-memory","status":"publish","type":"post","link":"http:\/\/drsfenner.org\/blog\/2016\/01\/into-the-weeds-iii-hacking-lapack-calls-to-save-memory\/","title":{"rendered":"Into the Weeds III – Hacking LAPACK Calls To Save Memory"},"content":{"rendered":"
\n
\n<\/div>\n
\n
\n

In the last post<\/a>, we tried to use lower-level LAPACK\/LAPACKE calls to remove some of the copying overhead of our operations. Well, it turns out that it didn’t work out so well (see the table at the end of the post). The fundamental reason why is that LAPACKE took the easy way out, when they decided to wrap LAPACK functions. In particular, whenever there is a C-\/Fortran-ordering issue (for a reminder<\/a>), LAPACKE simply makes a copy.<\/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 [1]:<\/div>\n
\n
\n
\n
import<\/span> numpy<\/span> as<\/span> np<\/span>\r\n\r\nfrom<\/span> shared<\/span> import<\/span> (<\/span>lap<\/span>,<\/span> get_dptr<\/span>,<\/span> create_data1<\/span>,<\/span> create_data2<\/span>,<\/span> simple_args<\/span>,<\/span> my_dlalsd<\/span>)<\/span>\r\nfrom<\/span> shared_dgelsd<\/span> import<\/span> make_structures_dgelsd<\/span>\r\n\r\n%<\/span>load_ext<\/span> memory_profiler\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
And now for something completely different …<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
Here we have two re-worked methods for solving least squares using either Fortran (col-major) data or, here’s the fun part folks, solving the least squares problem with C (row-major) data by calling the LAPACK\/E routines and lying<\/em>. We tell the Fortran routines that we have Fortran data (remember, in memory it is in row-major<\/em> order) and we solve a "transposed" problem<\/em>. That turns out to be just about trivial for using the \\(QR\\)<\/span> based method because the function call accepts a "T" (transposed) argument.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [2]:<\/div>\n
\n
\n
\n
def<\/span> sls_qr_via_dgels<\/span>(<\/span>X<\/span>,<\/span>Y<\/span>,<\/span>trans<\/span>):<\/span>\r\n    assert<\/span> X<\/span>.<\/span>flags<\/span>.<\/span>f_contiguous<\/span> or<\/span> trans<\/span> # don't allow the copy case although it should work<\/span>\r\n    # assert np.linalg.matrix_rank(X) == np.shape[1] # "full column rank"<\/span>\r\n\r\n    m<\/span>,<\/span>n<\/span>,<\/span>nrhs<\/span>,<\/span>ldX<\/span>,<\/span>ldY<\/span>,<\/span>la_order<\/span> =<\/span> simple_args<\/span>(<\/span>X<\/span>,<\/span>Y<\/span>,<\/span>trans<\/span>)<\/span>\r\n    X_p<\/span>,<\/span> Y_p<\/span> =<\/span> get_dptr<\/span>(<\/span>X<\/span>),<\/span> get_dptr<\/span>(<\/span>Y<\/span>)<\/span>\r\n\r\n    TN<\/span> =<\/span> "N"<\/span> if<\/span> not<\/span> trans<\/span> else<\/span> "T"<\/span>\r\n    info<\/span> =<\/span> lap<\/span>.<\/span>LAPACKE_dgels<\/span>(<\/span>la_order<\/span>,<\/span> TN<\/span>,<\/span> m<\/span>,<\/span>n<\/span>,<\/span>nrhs<\/span>,<\/span>\r\n                             X_p<\/span>,<\/span> ldX<\/span>,<\/span> Y_p<\/span>,<\/span> ldY<\/span>)<\/span>    \r\n\r\n    beta_shape<\/span> =<\/span> n<\/span> if<\/span> not<\/span> trans<\/span> else<\/span> m<\/span>\r\n    return<\/span> Y<\/span>[:<\/span>beta_shape<\/span>]<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
Manual Implementation of dgelsd<\/code> for C-order without Copies<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
However, for \\(SVD\\)<\/span>, things are a bit more tricky. I based this off of dgelsd<\/code> (the more "modern", slightly faster, slightly more memory hungry) \\(SVD\\)<\/span> least squares solver. Since the outer dgelsd<\/code> call doesn’t know about solving a transposed problem, we have to do it by hand<\/em>.<\/p>\n
Here’s a small aside. The entire<\/em> LAPACKE wrapper for LAPACK could adopt the techniques I’m using in this post and in many cases<\/em> (if not all?), could remove unnecessary copies. The problem is that copying and transposing gives a direct<\/em> way to use the preexisting implementation. The techniques I use here require rewriting<\/em> a (mathematically) equivalent implementation. LAPACK is unweidly (from a software engineering perspective), as is. It’s an exact case study for the labor savings of OOP (the tradeoff for performance is the main competingn concern). Adding another set of implementations is sure to causes headaches.<\/p>\n
The basic algorithm of dgelsd<\/code> given the regression problem \\(X \\beta = Y\\)<\/span> is:<\/p>\n
    \n
  1. Decompose X to two orthogonal matrices and an (upper) bidiagonal: \\(X \\rightarrow QBP\\)<\/span><\/li>\n
  2. Rotate \\(Y\\)<\/span> with \\(Q^T\\)<\/span> (this is done inplace, so it updates \\(Y \\leftarrow Q^T Y\\)<\/span>)<\/li>\n
  3. Solve \\(B\\beta_{tmp} = Q^T Y\\)<\/span> (written in original terms, don’t rotate \\(Y\\)<\/span> twice). Note, this process gives us both the \\(\\beta_{tmp}\\)<\/span> and<\/em> the singular values of \\(X\\)<\/span>.<\/li>\n
  4. Rotate \\(\\beta_{tmp}\\)<\/span> with \\(P\\)<\/span>: \\(\\beta \\leftarrow P\\beta_{tmp}\\)<\/span><\/li>\n<\/ol>\n
    These steps are performed with these Fortran LAPACK calls:<\/p>\n
      \n
    1. dgebrd<\/code><\/li>\n
    2. dormbr<\/code><\/li>\n
    3. dlalsd<\/code><\/li>\n
    4. dormbr<\/code><\/li>\n<\/ol>\n
      The combination of dgebrd<\/code> and dlalsd<\/code> performs an "applied" \\(SVD\\)<\/span>. dgebrd<\/code> actually computes \\(Q\\)<\/span> and \\(P^T\\)<\/span> (which is status quo in many numerical systems — probably because they either call LAPACK or want to be compatible).<\/p>\n
      One other step in dgelsd<\/code> is dealing with the case when there are many more rows than columns. Actually, dgelsd<\/code> handles under- and over-determined systems differently, but we’re assuming (and simplifying) to consider that we are either (1) square or (2) tall-and-skinny (over-determined). With many more rows than columsn, dgelsd<\/code> first performs a \\(QR\\)<\/span> decomposition (using dgeqrf<\/code>) on \\(X\\)<\/span> and then proceeds as above using \\(R\\)<\/span> in place of \\(X\\)<\/span> and \\(Q_{QR}^T Y\\)<\/span> for \\(Y\\)<\/span> (the rotation is applied using dormqr<\/code>).<\/p>\n<\/div>\n<\/div>\n<\/div>\n
      \n
      \n<\/div>\n
      \n
      \n
      That is all the "easy" case. We still have talked about what it means to do dgelsd<\/code> in a transposed sense. Remember, dgels<\/code> had an explicit T<\/code>\/N<\/code> flag that could control whether to transpose the system or not. Without futher ado, here’s the mathematical outline of steps:<\/p>\n
        \n
      1. \\(X^T \\rightarrow (QBP)^T = P^TB^TQ^T\\)<\/span> (as above with \\(Q \\leftrightarrow P\\)<\/span>; \\(B\\)<\/span> was upper bidiagonal, \\(B^T\\)<\/span> should be lower bidiagonal.<\/li>\n
      2. Rotate \\(Y\\)<\/span> with \\(P^T\\)<\/span>: \\(Y \\leftarrow P^TY\\)<\/span><\/li>\n
      3. Solve \\(B^T \\beta_{tmp} = P^TY\\)<\/span> giving \\(\\beta_{tmp}\\)<\/span> and the singular values<\/li>\n
      4. Rotate \\(\\beta_{tmp}\\)<\/span> with \\(Q\\)<\/span>: \\(\\beta \\leftarrow Q\\beta_{tmp}\\)<\/span><\/li>\n<\/ol>\n
        The analogue (tranpose) of the \\(QR\\)<\/span> step for the many rows case, is to perform an \\(LQ\\)<\/span> decomposition (the decomposition is done with dgelqf<\/code> and the rotation is applied with dormlq<\/code>). That is: \\(QR(X)\\)<\/span> and \\(LQ(X.T)\\)<\/span> give equivalent \\(Q\\)<\/span> values.<\/p>\n
        Getting all of this encoded properly was a royal pain in the tookis. The row-\/col-major ordering, the transposed problem, the transposes (or not) on the orthogonal matrices, the upper-\/lower- bidiagonalization (which is only computed as an upper because we have more rows than columns), and the differences in how LAPACK’s \\(QR\\)<\/span> and \\(LQ\\)<\/span> store their results made me feel like I was caught in an Escher painting.<\/p>\n
        In additional to encoding the mathematics (properly), it turned out that dlalsd<\/code> is not wrapped by the current version of LAPACKE. So, I had to write a small wrapper for it.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n<\/div>\n
        \n
        \n
        dlalsd Aside<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n<\/div>\n
        \n
        \n
        To make use of dlalsd<\/code>, I wrote a tiny C library. Then, I exposed it in Python using the cffi<\/code> technique<\/a>. I then gave it a Python function name of my_dlalsd<\/code>. My wrapper is to the LAPACK (FORTRAN) function, so it doesn’t have the C-ORDER "capabilities" (automatic copy and transpose) of other LAPACKE routines.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n<\/div>\n
        \n
        \n
        Here’s a quick look at the .h, .c, and Makefile and the Python helper. The dlalsd_<\/code> declaration is necesssary because it is defined externally (in libreflapack.so<\/code>).<\/p>\n
        \n
        \/* mylapack.h *\/<\/span>\r\nint<\/span> MY_dlalsd_work(char<\/span> UPLO, int<\/span> SMLSIZ, int<\/span> N, int<\/span> NRHS, \r\n                   double<\/span> *D, double<\/span> *E, double<\/span> *B,int<\/span> LDB,   \r\n                   double<\/span> RCOND, int<\/span> *rank, double<\/span> *work, int<\/span> *iwork);   \r\n\r\n\/* mylapack.c *\/<\/span>\r\n#include <stdio.h><\/span>\r\nvoid<\/span> dlalsd_(char<\/span> *UPLO, int<\/span> *SMLSIZ, int<\/span> *N, int<\/span> *NRHS,\r\n             double<\/span> *D, double<\/span> *E, double<\/span> *B,int<\/span> *LDB,\r\n             double<\/span> *RCOND, int<\/span> *rank, double<\/span> *work, int<\/span> *iwork, int<\/span> *info);\r\n\r\nint<\/span> MY_dlalsd_work(char<\/span> UPLO, int<\/span> SMLSIZ, int<\/span> N, int<\/span> NRHS,\r\n                   double<\/span> *D, double<\/span> *E, double<\/span> *B,int<\/span> LDB,\r\n                   double<\/span> RCOND, int<\/span> *rank, double<\/span> *work, int<\/span> *iwork){\r\n  int<\/span> info = 0<\/span>;\r\n  dlalsd_(&UPLO, &SMLSIZ, &N, &NRHS, D, E, B, &LDB, &RCOND, rank, work, iwork, &info);\r\n  if<\/span> (info < 0<\/span>)\r\n    info = info - 1<\/span>;\r\n\r\n  return<\/span> info;\r\n}<\/code><\/pre>\n<\/div>\n
        This Makefile worked on my Linux box with gcc.<\/p>\n
        \n
        # Makefile<\/span>\r\nall<\/span>:\r\n    gcc<\/span> -shared -fPIC -o libmylapack.so mylapack.c -lreflapack<\/code><\/pre>\n<\/div>\n
        And a little Python shim.<\/p>\n
        \n
        from<\/span> cffi import<\/span> FFI\r\nffi =<\/span> FFI()\r\nwith<\/span> open<\/span>(".\/mylapack.h"<\/span>) as<\/span> header_file:\r\n    ffi.cdef(header_file.read())\r\nmylapack =<\/span> ffi.dlopen(".\/libmylapack.so"<\/span>)\r\nmy_dlalsd =<\/span> mylapack.MY_dlalsd_work<\/code><\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n<\/div>\n
        \n
        \n
        The Code<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n
        In [3]:<\/div>\n
        \n
        \n
        \n
        def<\/span> sls_svd_via_manual_dgelsd<\/span>(<\/span>X<\/span>,<\/span>Y<\/span>,<\/span>trans<\/span>):<\/span>\r\n    ''' we essentially reimplement dgelsd and add the capability of solving a transposed problem<\/span>\r\n        i.e., $X \\Beta = Y$ --or-- $X.T \\Beta = Y$ '''<\/span>\r\n    assert<\/span> X<\/span>.<\/span>flags<\/span>.<\/span>f_contiguous<\/span> or<\/span> trans<\/span>    \r\n    common_args<\/span>,<\/span> arrays<\/span>,<\/span> pointers<\/span> =<\/span> make_structures_dgelsd<\/span>(<\/span>X<\/span>,<\/span>Y<\/span>,<\/span>trans<\/span>=<\/span>trans<\/span>)<\/span>\r\n\r\n    m<\/span>,<\/span>n<\/span>,<\/span>nrhs<\/span>,<\/span>ldX<\/span>,<\/span>ldY<\/span>,<\/span>la_order<\/span>                       =<\/span> common_args<\/span>\r\n    singular_values<\/span>,<\/span> E<\/span>,<\/span> Q<\/span>,<\/span> P<\/span>,<\/span> rank<\/span>,<\/span> M<\/span>,<\/span> I<\/span>            =<\/span> arrays<\/span>\r\n    X_p<\/span>,<\/span> Y_p<\/span>,<\/span> SV_p<\/span>,<\/span> E_p<\/span>,<\/span> Q_p<\/span>,<\/span> P_p<\/span>,<\/span> rank_p<\/span>,<\/span> M_p<\/span>,<\/span> I_p<\/span> =<\/span> pointers<\/span>\r\n\r\n    # helpful names<\/span>\r\n    minmn<\/span> =<\/span> min<\/span>(<\/span>m<\/span>,<\/span>n<\/span>)<\/span>\r\n    real_rows<\/span>,<\/span> real_cols<\/span> =<\/span> (<\/span>m<\/span>,<\/span>n<\/span>)<\/span> if<\/span> not<\/span> trans<\/span> else<\/span> (<\/span>n<\/span>,<\/span>m<\/span>)<\/span>\r\n    y_shape<\/span>,<\/span> beta_shape<\/span> =<\/span> real_rows<\/span>,<\/span> real_cols<\/span>\r\n\r\n    UL_arg<\/span> =<\/span> 'U'<\/span> # without preprocessing, both trans={T,F} need 'U'<\/span>\r\n\r\n    # cross over for pre-QRing is defined at:<\/span>\r\n    # tsiee.f line 576 ... thresh = min(m,n) * 1.6<\/span>\r\n    lotsaRows<\/span> =<\/span> real_rows<\/span> >=<\/span> real_cols<\/span> *<\/span> 1.6<\/span>\r\n    if<\/span> lotsaRows<\/span>:<\/span>\r\n        if<\/span> not<\/span> trans<\/span>:<\/span>\r\n            decomp<\/span>,<\/span> mult<\/span>,<\/span> TN<\/span> =<\/span> lap<\/span>.<\/span>LAPACKE_dgeqrf_work<\/span>,<\/span> lap<\/span>.<\/span>LAPACKE_dormqr_work<\/span>,<\/span> 'T'<\/span>\r\n        else<\/span>:<\/span>\r\n            decomp<\/span>,<\/span> mult<\/span>,<\/span> TN<\/span> =<\/span> lap<\/span>.<\/span>LAPACKE_dgelqf_work<\/span>,<\/span> lap<\/span>.<\/span>LAPACKE_dormlq_work<\/span>,<\/span> 'N'<\/span>\r\n            UL_arg<\/span> =<\/span> 'L'<\/span> # B\/E need to be used to define a lower bi-di<\/span>\r\n\r\n        info<\/span> =<\/span> decomp<\/span>(<\/span>la_order<\/span>,<\/span> m<\/span>,<\/span> n<\/span>,<\/span> X_p<\/span>,<\/span> ldX<\/span>,<\/span> E_p<\/span>,<\/span> M_p<\/span>,<\/span> M<\/span>.<\/span>size<\/span>)<\/span>\r\n        assert<\/span> not<\/span> info<\/span>,<\/span> "<\/span>%s<\/span>:<\/span>%s<\/span>"<\/span> %<\/span> (<\/span>decomp<\/span>.<\/span>__name__<\/span>,<\/span> info<\/span>)<\/span>\r\n\r\n        # this can be done at the end (so says dgelsd.f) <\/span>\r\n        info<\/span> =<\/span> mult<\/span>(<\/span>la_order<\/span>,<\/span> 'L'<\/span>,<\/span> TN<\/span>,<\/span> y_shape<\/span>,<\/span> nrhs<\/span>,<\/span> real_cols<\/span>,<\/span> X_p<\/span>,<\/span> ldX<\/span>,<\/span> E_p<\/span>,<\/span> Y_p<\/span>,<\/span> ldY<\/span>,<\/span> M_p<\/span>,<\/span> M<\/span>.<\/span>size<\/span>)<\/span>\r\n        assert<\/span> not<\/span> info<\/span>,<\/span> "<\/span>%s<\/span>:<\/span>%s<\/span>"<\/span> %<\/span> (<\/span>decomp<\/span>.<\/span>__name__<\/span>,<\/span> info<\/span>)<\/span>\r\n            \r\n        y_shape<\/span> =<\/span> beta_shape<\/span> =<\/span> n<\/span> =<\/span> m<\/span> =<\/span> minmn<\/span>\r\n        X<\/span>[<\/span>np<\/span>.<\/span>tril_indices<\/span>(<\/span>minmn<\/span>,<\/span> -<\/span>1<\/span>)]<\/span> =<\/span> 0.0<\/span>\r\n\r\n    # print "calling dgebrd"<\/span>\r\n    assert<\/span> M<\/span>.<\/span>size<\/span> >=<\/span> max<\/span>(<\/span>m<\/span>,<\/span>n<\/span>)<\/span> # minimal workspace<\/span>\r\n    info<\/span> =<\/span> lap<\/span>.<\/span>LAPACKE_dgebrd_work<\/span>(<\/span>la_order<\/span>,<\/span> m<\/span>,<\/span> n<\/span>,<\/span> X_p<\/span>,<\/span> ldX<\/span>,<\/span> SV_p<\/span>,<\/span> E_p<\/span>,<\/span> Q_p<\/span>,<\/span> P_p<\/span>,<\/span> M_p<\/span>,<\/span> M<\/span>.<\/span>size<\/span>)<\/span>\r\n    assert<\/span> not<\/span> info<\/span>,<\/span> "dgebrd bad arg <\/span>%d<\/span>"<\/span> %<\/span> info<\/span>\r\n\r\n    pq<\/span>,<\/span> PQ_p<\/span>,<\/span> real_k<\/span> =<\/span> (<\/span>'Q'<\/span>,<\/span> Q_p<\/span>,<\/span> n<\/span>)<\/span>  if<\/span> not<\/span> trans<\/span> else<\/span> (<\/span>'P'<\/span>,<\/span> P_p<\/span>,<\/span> m<\/span>)<\/span>\r\n    # print "calling dormbr (I)"  <\/span>\r\n    info<\/span> =<\/span> lap<\/span>.<\/span>LAPACKE_dormbr_work<\/span>(<\/span>la_order<\/span>,<\/span> pq<\/span>,<\/span> 'L'<\/span>,<\/span> 'T'<\/span>,<\/span> y_shape<\/span>,<\/span> nrhs<\/span>,<\/span> real_cols<\/span>,<\/span>\r\n                                   X_p<\/span>,<\/span> ldX<\/span>,<\/span> PQ_p<\/span>,<\/span> Y_p<\/span>,<\/span> ldY<\/span>,<\/span> M_p<\/span>,<\/span> M<\/span>.<\/span>size<\/span>)<\/span>\r\n    assert<\/span> not<\/span> info<\/span>,<\/span> "dormbr (I) bad arg <\/span>%d<\/span>"<\/span> %<\/span> info<\/span>\r\n    \r\n    # print "calling dlalsd"<\/span>\r\n    # FIXME: this wraps a direct call to LAPACK (broken for C-ORDER)<\/span>\r\n    # 25 is "small size" for this dlalsd<\/span>\r\n    # -1 says to use machine eps for rcond parameter (det. small singular values)<\/span>\r\n    info<\/span> =<\/span> my_dlalsd<\/span>(<\/span>UL_arg<\/span>,<\/span> 25<\/span>,<\/span> minmn<\/span>,<\/span> nrhs<\/span>,<\/span> SV_p<\/span>,<\/span> E_p<\/span>,<\/span> Y_p<\/span>,<\/span> ldY<\/span>,<\/span> -<\/span>1<\/span>,<\/span> rank_p<\/span>,<\/span> M_p<\/span>,<\/span> I_p<\/span>)<\/span>\r\n    assert<\/span> not<\/span> info<\/span>,<\/span> "dlalsd bad arg <\/span>%d<\/span>"<\/span> %<\/span> info<\/span>\r\n    \r\n    pq<\/span>,<\/span> PQ_p<\/span> =<\/span> (<\/span>'P'<\/span>,<\/span> P_p<\/span>)<\/span> if<\/span> not<\/span> trans<\/span> else<\/span> (<\/span>'Q'<\/span>,<\/span> Q_p<\/span>)<\/span>\r\n    # print "calling dormbr (II)"<\/span>\r\n    info<\/span> =<\/span> lap<\/span>.<\/span>LAPACKE_dormbr_work<\/span>(<\/span>la_order<\/span>,<\/span> pq<\/span>,<\/span> 'L'<\/span>,<\/span> 'N'<\/span>,<\/span> beta_shape<\/span>,<\/span> nrhs<\/span>,<\/span> n<\/span>,<\/span>\r\n                                   X_p<\/span>,<\/span> ldX<\/span>,<\/span> PQ_p<\/span>,<\/span> Y_p<\/span>,<\/span> ldY<\/span>,<\/span> M_p<\/span>,<\/span> M<\/span>.<\/span>size<\/span>)<\/span>\r\n    assert<\/span> not<\/span> info<\/span>,<\/span> "dormbr (II) bad arg <\/span>%d<\/span>"<\/span> %<\/span> info<\/span>\r\n\r\n    return<\/span> Y<\/span>[:<\/span>beta_shape<\/span>],<\/span> singular_values<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n<\/div>\n
        \n
        \n
        Testing<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n<\/div>\n
        \n
        \n
        There were more than enough corner cases (and missteps) while I was coding to actually want test cases so I didn’t "fix" problems with regressions (the software engineering kind) in other parts of the code. So, some unittests<\/code> to the rescue:<\/p>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n
        In [4]:<\/div>\n
        \n
        \n
        \n
        !<\/span>python test_dgelsd.py\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n
        \n
        <\/div>\n
        \n
        test_heavy_pure_fortran_1 (__main__.TestManual_dgelsd) ... ok\r\ntest_heavy_pure_fortran_2 (__main__.TestManual_dgelsd) ... ok\r\ntest_heavy_tranposed_corder_1 (__main__.TestManual_dgelsd) ... ok\r\ntest_heavy_tranposed_corder_2 (__main__.TestManual_dgelsd) ... ok\r\ntest_medium_pure_fortran_1 (__main__.TestManual_dgelsd) ... ok\r\ntest_medium_pure_fortran_2 (__main__.TestManual_dgelsd) ... ok\r\ntest_medium_tranposed_corder_1 (__main__.TestManual_dgelsd) ... ok\r\ntest_medium_tranposed_corder_2 (__main__.TestManual_dgelsd) ... ok\r\ntest_simple_pure_fortran (__main__.TestManual_dgelsd) ... ok\r\ntest_simple_transposed_corder (__main__.TestManual_dgelsd) ... ok\r\n\r\n----------------------------------------------------------------------\r\nRan 10 tests in 40.326s\r\n\r\nOK\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n<\/div>\n
        \n
        \n
        Same Results?<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n<\/div>\n
        \n
        \n
        So, let’s run each of these routines with row- and column-major data and verify we get the same output.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n
        In [5]:<\/div>\n
        \n
        \n
        \n
        #<\/span>\r\n# we have to copy the data so we can use the same data in the next call<\/span>\r\n# lapack\/e overwrite its inputs with the results (decomps and betas)<\/span>\r\n#<\/span>\r\ndef<\/span> copy_and_do<\/span>(<\/span>func<\/span>,<\/span> X<\/span>,<\/span> Y<\/span>,<\/span> trans<\/span>,<\/span> order<\/span>):<\/span>\r\n    _X<\/span>,<\/span>_Y<\/span> =<\/span> np<\/span>.<\/span>copy<\/span>(<\/span>X<\/span>,<\/span>order<\/span>=<\/span>order<\/span>),<\/span> np<\/span>.<\/span>copy<\/span>(<\/span>Y<\/span>,<\/span>order<\/span>=<\/span>order<\/span>)<\/span>\r\n    return<\/span> func<\/span>(<\/span>_X<\/span>,<\/span>_Y<\/span>,<\/span>trans<\/span>)<\/span>\r\n\r\nn_obs<\/span>,<\/span> n_ftrs<\/span> =<\/span> 2<\/span>**<\/span>16<\/span>,<\/span> 2<\/span>**<\/span>10<\/span>\r\n#n_obs, n_ftrs = 10, 2<\/span>\r\n#n_obs, n_ftrs = 1024, 50<\/span>\r\n\r\n# X,Y = create_data2()<\/span>\r\nX<\/span>,<\/span>Y<\/span> =<\/span> create_data1<\/span>(<\/span>n_obs<\/span>,<\/span> n_ftrs<\/span>)<\/span>\r\nbetas_qr_f<\/span>  =<\/span> copy_and_do<\/span>(<\/span>sls_qr_via_dgels<\/span>,<\/span>          X<\/span>,<\/span> Y<\/span>,<\/span> False<\/span>,<\/span> "F"<\/span>)<\/span>\r\nbetas_svd_f<\/span> =<\/span> copy_and_do<\/span>(<\/span>sls_svd_via_manual_dgelsd<\/span>,<\/span> X<\/span>,<\/span> Y<\/span>,<\/span> False<\/span>,<\/span> "F"<\/span>)[<\/span>0<\/span>]<\/span>\r\n\r\nbetas_qr_c<\/span>  =<\/span> copy_and_do<\/span>(<\/span>sls_qr_via_dgels<\/span>,<\/span>          X<\/span>,<\/span> Y<\/span>,<\/span> True<\/span>,<\/span> "C"<\/span>)<\/span>\r\nbetas_svd_c<\/span> =<\/span> copy_and_do<\/span>(<\/span>sls_svd_via_manual_dgelsd<\/span>,<\/span> X<\/span>,<\/span> Y<\/span>,<\/span> True<\/span>,<\/span> "C"<\/span>)[<\/span>0<\/span>]<\/span>\r\n\r\nbetas_expected<\/span> =<\/span> np<\/span>.<\/span>linalg<\/span>.<\/span>lstsq<\/span>(<\/span>X<\/span>,<\/span>Y<\/span>)[<\/span>0<\/span>]<\/span>\r\n\r\n# betas = [betas_qr_f, betas_svd_f[0], betas_qr_c, betas_svd_c[0]]<\/span>\r\nbetas<\/span> =<\/span> [<\/span>betas_expected<\/span>,<\/span> betas_qr_f<\/span>,<\/span> betas_qr_c<\/span>,<\/span> betas_svd_f<\/span>,<\/span> betas_svd_c<\/span>]<\/span>\r\n#for b in betas:<\/span>\r\n#    print b<\/span>\r\nprint<\/span> "All same betas?:"<\/span>,<\/span> all<\/span>(<\/span>np<\/span>.<\/span>allclose<\/span>(<\/span>b1<\/span>,<\/span> b2<\/span>)<\/span> for<\/span> b1<\/span> in<\/span> betas<\/span>[:<\/span>-<\/span>1<\/span>]<\/span> for<\/span> b2<\/span> in<\/span> betas<\/span>[<\/span>1<\/span>:])<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n
        \n
        <\/div>\n
        \n
        All same betas?: True\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n<\/div>\n
        \n
        \n
        Memory Usage<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n
        In [6]:<\/div>\n
        \n
        \n
        \n
        n_obs<\/span>,<\/span> n_ftrs<\/span> =<\/span> 2<\/span>**<\/span>16<\/span>,<\/span> 2<\/span>**<\/span>10<\/span>\r\nX<\/span>,<\/span>Y<\/span> =<\/span> create_data1<\/span>(<\/span>n_obs<\/span>,<\/span> n_ftrs<\/span>,<\/span> "C"<\/span>)<\/span>\r\n# %memit sls_qr_via_dgels(X,Y,True)           # peak memory: 600.68 MiB, increment: 3.19 MiB<\/span>\r\n# %memit sls_svd_via_manual_dgelsd(X,Y,True)  # peak memory: 584.77 MiB, increment: 20.05 MiB<\/span>\r\n# for comparison:<\/span>\r\n# %memit np.linalg.lstsq(X,Y)                 # peak memory: 1082.57 MiB, increment: 517.70 MiB<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n<\/div>\n
        \n
        \n
        And now, it’s time for a big happy dance with music: we avoided copies AND<\/em> we stayed in row-major (C-land) order!<\/p>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n<\/div>\n
        \n
        \n
        Timing<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n
        In [7]:<\/div>\n
        \n
        \n
        \n
        n_obs<\/span>,<\/span> n_ftrs<\/span> =<\/span> 2<\/span>**<\/span>16<\/span>,<\/span> 2<\/span>**<\/span>10<\/span>\r\n_X<\/span>,<\/span>_Y<\/span> =<\/span> create_data1<\/span>(<\/span>n_obs<\/span>,<\/span> n_ftrs<\/span>,<\/span> "C"<\/span>)<\/span>\r\n\r\n%<\/span>timeit<\/span> np.linalg.lstsq(_X,_Y)\r\n\r\n# forcing copies since all of these calls place the results in the input (overwrite our data)<\/span>\r\nX<\/span>,<\/span>Y<\/span> =<\/span> np<\/span>.<\/span>copy<\/span>(<\/span>_X<\/span>),<\/span> np<\/span>.<\/span>copy<\/span>(<\/span>_Y<\/span>)<\/span>\r\n%<\/span>timeit<\/span> sls_qr_via_dgels(X,Y,True)\r\n\r\nX<\/span>,<\/span>Y<\/span> =<\/span> np<\/span>.<\/span>copy<\/span>(<\/span>_X<\/span>),<\/span> np<\/span>.<\/span>copy<\/span>(<\/span>_Y<\/span>)<\/span>\r\n%<\/span>timeit<\/span> sls_svd_via_manual_dgelsd(X,Y,True)\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n
        \n
        <\/div>\n
        \n
        1 loops, best of 3: 8.72 s per loop\r\n1 loops, best of 3: 11.9 s per loop\r\n1 loops, best of 3: 11.7 s per loop\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n
        In [8]:<\/div>\n
        \n
        \n
        \n
        import<\/span> scipy.linalg<\/span> as<\/span> sla<\/span>\r\n\r\nX<\/span>,<\/span>Y<\/span> =<\/span> np<\/span>.<\/span>copy<\/span>(<\/span>_X<\/span>),<\/span> np<\/span>.<\/span>copy<\/span>(<\/span>_Y<\/span>)<\/span>\r\n# NOTE: sla.lstsq uses dgelss insead of dgelsd<\/span>\r\n%<\/span>timeit<\/span> sla.lstsq(X,Y,overwrite_a=True,overwrite_b=True)\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n
        \n
        <\/div>\n
        \n
        1 loops, best of 3: 28.7 s per loop\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n<\/div>\n
        \n
        \n
        A few things stand out here.<\/p>\n
        Even though np.linalg.lstsq<\/code> makes a copy (and it also has additional memory overhead to return its other results – it doesn’t overwrite), it is still faster at this problem size. I’ll tinker with finding a tradeoff point, but remember the point here was (largely) about memory and less about processing time (although the two are intimately intertwined).<\/p>\n
        My Python sls<\/code> functions are not particularly optimized. The sls_qr_via_dgels<\/code> function is very small. The sls_svd_via_manual_dgelsd<\/code>, however, has a lot of machinery on the inside. I’ve left some of that machinery out of this post (which is already quite long and detailed). But, there is memory allocation, workspace (temporary) size calculation, and a variety of utility routines that I’ve only lightly tuned. Most of those are only executed once, so we’re not talking about inner loop performance hits.<\/p>\n
        Next, while scipy.linalg.lstsq<\/code> can operate with overwriting, it uses a different underlying algorithm dgelss<\/code> instead of dgelsd<\/code>. That algorithm is a bit slower (than dgelsd<\/code>), but it uses less memory.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n<\/div>\n
        \n
        \n
        LAPACK Glossary<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n<\/div>\n
        \n
        \n
        There are many LAPACK routine names in the post. Here’s a quick guide to interpreting them.<\/p>\n
        Important abbreviations:<\/p>\n