{"id":390,"date":"2015-12-17T19:33:59","date_gmt":"2015-12-17T19:33:59","guid":{"rendered":"http:\/\/drsfenner.org\/blog\/?p=390"},"modified":"2016-02-05T16:28:51","modified_gmt":"2016-02-05T16:28:51","slug":"into-the-weeds-i-lapack-lapacke-and-three-paths","status":"publish","type":"post","link":"http:\/\/drsfenner.org\/blog\/2015\/12\/into-the-weeds-i-lapack-lapacke-and-three-paths\/","title":{"rendered":"Into the Weeds I – LAPACK, LAPACKE, and Three+ Paths"},"content":{"rendered":"
\n
\n<\/div>\n
\n
\n

In the last post<\/a>, we looked at three ways of solving the least-squares regression problem (Cholesky, QR, and SVD decompositions) using LAPACK routines in a round-about fashion. The different decompositions (i.e., factorizations) were performed by (1) calling a NumPy np.linalg<\/code> routine that called (2) an LAPACK driver routine<\/a> for a "bigger" problem like linear equation solving<\/a> or least-squares<\/a>. The driver routine called one of the decompositions internally. <\/p>\n<\/div>\n<\/div>\n<\/div>\n

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

LAPACK Routines<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n

Since the QR and SVD solutions both used a copy<\/em> of our \\(X\\)<\/span> matrix (or at least the memory used hinted at that), I want to dive into the lower level of LAPACK (avoiding NumPy) and call the driver routines directly. We’ll try to avoid doing the decompositions by hand, if we can!. The most immediate issue, in this series of experiments is: how do we call the LAPACK routines?<\/em> Here is a table of the LAPACK routines used by the solver\/least squares routines in NumPy, SciPy, and cvxopt.<\/p>\n\n\n\n\n\n\n\n\n\n
Python
Library
Routine<\/th>\n		Underlying
LAPACK
Routine<\/th>\n		Matrix Type<\/th>\nAlgorithm<\/th>\n<\/tr>\n<\/thead>\n
sla.solve<\/code>
w\/ sym_pos=True<\/code><\/td>\n		posv<\/code><\/td>\nsymmetric positive
definite matrix<\/td>\n		solve eqn using Cholesky<\/td>\n<\/tr>\n
nla.solve<\/code>
sla.solve<\/code><\/td>\n		gesv<\/code><\/td>\nsquare matrix<\/td>\nsolve eqn using LU<\/td>\n<\/tr>\n
cla.solve<\/code><\/td>\ngels<\/code><\/td>\ngeneral matrix (full rank)<\/td>\nleast-squares via QR<\/td>\n<\/tr>\n
sla.lstsq<\/code><\/td>\ngelss<\/code><\/td>\ngeneral matrix<\/td>\nleast-squares via SVD<\/td>\n<\/tr>\n
nla.lstsq<\/code><\/td>\ngelsd<\/code><\/td>\ngeneral matrix<\/td>\nleast-squares via Divide-and-Conquer SVD<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n

The different libraries make use of different LAPACK routines. And, in general, the libraries only expose the LAPACK routines that they need. For example, scipy.linalg<\/code> does not provide dgels<\/code> (as of scipy-0.16.1) it is not listed in the scipy.linalg.lapack function list<\/a>. It is almost possible to pick and choose from each library to get the desired, underlying LAPACK routines, but I want to do something slightly different.<\/p>\n<\/div>\n<\/div>\n<\/div>\n

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

LAPACK: Fortran to C (Sharks\/Jets? Montague\/Capulet?)<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n

One of the biggest obstacles to using LAPACK in Python is the fact that Python (at least CPython) and NumPy (and friends) are written in C. LAPACK is written in Fortran. This isn’t a show stopper. There are a number of interoperabiliity options out there: (1) converting Fortran to C (f2c<\/code> program) and (2) linking Fortran binaries into a C program. The first is just a fancy pure-C program. The second generally requires declaring a C function prototype for the Fortran function you want to call (with something like: extern "C" void dgels(arg1, arg2, etc.)<\/code>) and then simply calling the function. Here’s a quick overview<\/a>.<\/p>\n

Declaring your own prototypes is a real drag and it can also leave some holes in the conversion between C-land and Fortran-land. The biggest hole is the fact that Fortran stores data in column-major<\/em> order, while C stores data in row-major<\/em> order (wikipedia:Row-major<\/a>). There are other issues as well and having a well-defined interface between the two layers is a really good idea. BLAS (lower level routines than LAPACK) has had one for years called CBLAS<\/a>. However, it took significantly longer for LAPACK to get a standard C-interface: LAPACKE<\/a>. A quick side-note: CLAPACK<\/a> is a f2c<\/code>‘d version of LAPACK, so it provides similar functionality in a different way.<\/p>\n<\/div>\n<\/div>\n<\/div>\n

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

LAPACKE and Python<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n

After all that lead-in, my goal is to use LAPACKE to call LAPACK routines — from Python. It turns out that this is really easy. py-liblapack<\/a> is a tiny shim built using a Python CFFI library called cffi<\/a>. Here, CFFI stands for C foreign function interface<\/em> and that acronym shows up all over the place. One example is LISP’s CFFI<\/a> where the C stands for Common<\/em>. Anyway, cffi<\/code> let’s us do some really cool things:<\/p>\n

    \n
  1. Define and allocate C-level objects (variables). For example, workspace = ffi.new("int[]", n)<\/code> creates a Python reference to a length-n<\/code> int<\/code> array<\/li>\n
  2. Read in declarations from a .h<\/code> file: header_file = open("liblapack.h"); ffi.cdef(header_file.read())<\/code><\/li>\n
  3. Access functions defined in a C library. For example, liblapack = ffi.dlopen("libreflapacke.so")<\/code> gives us dot-access (via Python attributes) to functions in libreflapacke<\/code>. We can use dgels<\/code> by coding: liblapack.LAPACKE_dgels<\/code>.<\/li>\n<\/ol>\n

    These three capabilities let mikefc<\/em> write an ~10 line script that exposes LAPACKE for us in Python. Granted, more work went into making a cffi<\/code>-compatible set of .h<\/code> files for the LAPACKE routines, but that was done by another set of Python scripts. py-liblapack<\/code> exposes all of the single<\/code> and double<\/code> (s<\/code> and d<\/code>) routines from LAPACKE; it would be straight-forward to also expose the complex<\/code> and double-precision complex<\/code> (c<\/code> and z<\/code>) routines.<\/p>\n

    To use the py-liblapack<\/code> code, all you need are two files from the source repository: liblapack.py<\/code> and liblapack.h<\/code>. Then, put the proper file names in liblapack.py<\/code>:<\/p>\n

    \n
    import<\/span> os\r\nfrom<\/span> cffi import<\/span> FFI\r\nffi =<\/span> FFI()\r\n\r\nthis_dir =<\/span> os.path.dirname(os.path.abspath(__file__<\/span>)) +<\/span> "\/"<\/span>\r\nwith<\/span> open<\/span>(this_dir +<\/span> "liblapack.h"<\/span>) as<\/span> header_file:\r\n    ffi.cdef(header_file.read())\r\nliblapack =<\/span> ffi.dlopen("libreflapacke.so"<\/span>)<\/code><\/pre>\n<\/div>\n
    In the examples below, I’ll assume that liblapack.py<\/code> lives in the same directory as my notebook and I’ll be import<\/code>ing it, shortly.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Generic Setup and Imports<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [1]:<\/div>\n
    \n
    \n
    \n
    import<\/span> operator<\/span> as<\/span> op<\/span>\r\n\r\nimport<\/span> numpy<\/span> as<\/span> np<\/span>\r\nimport<\/span> numpy.linalg<\/span> as<\/span> nla<\/span>\r\nimport<\/span> scipy.linalg<\/span> as<\/span> sla<\/span>\r\n\r\nfrom<\/span> cvxopt<\/span> import<\/span> matrix<\/span> as<\/span> cvx_matrix<\/span>\r\nimport<\/span> cvxopt.lapack<\/span> as<\/span> cla<\/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
    In [2]:<\/div>\n
    \n
    \n
    \n
    def<\/span> create_data<\/span>(<\/span>n<\/span>=<\/span>5000<\/span>,<\/span> p<\/span>=<\/span>100<\/span>):<\/span>\r\n    ''' <\/span>\r\n        n is number cases\/observations\/examples<\/span>\r\n        p is number of features\/attributes\/variables<\/span>\r\n    '''<\/span>\r\n    X<\/span>      =<\/span> np<\/span>.<\/span>random<\/span>.<\/span>uniform<\/span>(<\/span>0<\/span>,<\/span>10<\/span>,<\/span> (<\/span>n<\/span>,<\/span>p<\/span>))<\/span>\r\n    coeffs<\/span> =<\/span> np<\/span>.<\/span>random<\/span>.<\/span>uniform<\/span>(<\/span>0<\/span>,<\/span>10<\/span>,<\/span> p<\/span>)<\/span>\r\n    def<\/span> f<\/span>(<\/span>X<\/span>):<\/span> return<\/span> X<\/span>.<\/span>dot<\/span>(<\/span>coeffs<\/span>)<\/span>\r\n\r\n    noise<\/span> =<\/span> np<\/span>.<\/span>random<\/span>.<\/span>normal<\/span>(<\/span>0.0<\/span>,<\/span> 2.0<\/span>,<\/span> n<\/span>)<\/span>\r\n    Y<\/span> =<\/span> f<\/span>(<\/span>X<\/span>)<\/span> +<\/span> noise<\/span>\r\n    Y<\/span> =<\/span> Y<\/span>.<\/span>reshape<\/span>((<\/span>n<\/span>,<\/span>1<\/span>))<\/span>\r\n    \r\n    return<\/span> X<\/span>,<\/span>Y<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [3]:<\/div>\n
    \n
    \n
    \n
    n_obs<\/span> =<\/span> 2<\/span>**<\/span>16<\/span>\r\nn_ftrs<\/span> =<\/span> 2<\/span>**<\/span>10<\/span>\r\nX<\/span>,<\/span>Y<\/span> =<\/span> create_data<\/span>(<\/span>n_obs<\/span>,<\/span>n_ftrs<\/span>)<\/span>\r\nprint<\/span> "about"<\/span>,<\/span> 10<\/span>**<\/span>int<\/span>(<\/span>np<\/span>.<\/span>log10<\/span>(<\/span>n_obs<\/span> *<\/span> n_ftrs<\/span>)),<\/span> "entries"<\/span>\r\nprint<\/span> "about"<\/span>,<\/span> (<\/span>X<\/span>.<\/span>itemsize<\/span> *<\/span> X<\/span>.<\/span>size<\/span>)<\/span> \/<\/span> float<\/span>(<\/span>1024<\/span>**<\/span>2<\/span>),<\/span> "MB"<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    \n
    <\/div>\n
    \n
    about 10000000 entries\r\nabout 512.0 MB\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    liblapack<\/code> and Utilities<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [4]:<\/div>\n
    \n
    \n
    \n
    from<\/span> liblapack<\/span> import<\/span> liblapack<\/span> as<\/span> lap<\/span>\r\nfrom<\/span> cffi<\/span> import<\/span> FFI<\/span>\r\nffi<\/span> =<\/span> FFI<\/span>()<\/span>\r\n\r\n#<\/span>\r\n# numpy array to underlying C array<\/span>\r\n#<\/span>\r\ndef<\/span> get_dptr<\/span>(<\/span>arr<\/span>):<\/span>\r\n    return<\/span> ffi<\/span>.<\/span>cast<\/span>(<\/span>"double *"<\/span>,<\/span> arr<\/span>.<\/span>__array_interface__<\/span>[<\/span>'data'<\/span>][<\/span>0<\/span>])<\/span>\r\n\r\ndef<\/span> get_iptr<\/span>(<\/span>arr<\/span>):<\/span>\r\n    return<\/span> ffi<\/span>.<\/span>cast<\/span>(<\/span>"int *"<\/span>,<\/span> arr<\/span>.<\/span>__array_interface__<\/span>[<\/span>'data'<\/span>][<\/span>0<\/span>])<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Cholesky via dposv<\/code><\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [5]:<\/div>\n
    \n
    \n
    \n
    # posv<\/span>\r\n# http:\/\/www.netlib.org\/lapack\/explore-html\/dc\/d7f\/lapacke__dposv_8c.html<\/span>\r\n# https:\/\/software.intel.com\/en-us\/node\/520982<\/span>\r\n# <\/span>\r\n# lapack_int LAPACKE_dposv(int matrix_layout,<\/span>\r\n#                          char uplo,<\/span>\r\n#                          lapack_int n, lapack_int nrhs,<\/span>\r\n#                          double *a, lapack_int lda,<\/span>\r\n#                          double *b, lapack_int ldb)<\/span>\r\ndef<\/span> sls_cholesky_orig<\/span>(<\/span>X<\/span>,<\/span> Y<\/span>):<\/span>\r\n    XtX<\/span> =<\/span> X<\/span>.<\/span>T<\/span>.<\/span>dot<\/span>(<\/span>X<\/span>)<\/span>\r\n    XtY<\/span> =<\/span> X<\/span>.<\/span>T<\/span>.<\/span>dot<\/span>(<\/span>Y<\/span>)<\/span>\r\n\r\n    m<\/span>,<\/span>n<\/span> =<\/span> XtX<\/span>.<\/span>shape<\/span>\r\n    nrhs<\/span> =<\/span> XtY<\/span>.<\/span>shape<\/span>[<\/span>1<\/span>]<\/span> if<\/span> XtY<\/span>.<\/span>ndim<\/span> ==<\/span> 2<\/span> else<\/span> 1<\/span>\r\n\r\n    info<\/span> =<\/span> lap<\/span>.<\/span>LAPACKE_dposv<\/span>(<\/span>lap<\/span>.<\/span>LAPACK_ROW_MAJOR<\/span>,<\/span>\r\n                             "L"<\/span>,<\/span> # prefer L.dot(L.T) factorization<\/span>\r\n                             n<\/span>,<\/span> nrhs<\/span>,<\/span>\r\n                             get_dptr<\/span>(<\/span>XtX<\/span>),<\/span> n<\/span>,<\/span>    # ldXtX = n<\/span>\r\n                             get_dptr<\/span>(<\/span>XtY<\/span>),<\/span> nrhs<\/span>)<\/span> # ldXtY = nrhs<\/span>\r\n    assert<\/span> not<\/span> info<\/span>\r\n    return<\/span> XtY<\/span>[:<\/span>X<\/span>.<\/span>shape<\/span>[<\/span>1<\/span>]]<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [6]:<\/div>\n
    \n
    \n
    \n
    # %memit sls_cholesky_orig(X,Y)<\/span>\r\n# peak memory: 581.86 MiB, increment: 11.87 MiB<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    LU via dgesv<\/code><\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [7]:<\/div>\n
    \n
    \n
    \n
    # gesv<\/span>\r\n# http:\/\/www.netlib.org\/lapack\/explore-html\/d6\/de8\/lapacke__dgesv_8c.html<\/span>\r\n# https:\/\/software.intel.com\/en-us\/node\/520973<\/span>\r\n#<\/span>\r\n#lapack_int LAPACKE_dgesv(int matrix_layout,<\/span>\r\n#                         lapack_int n, lapack_int nrhs,<\/span>\r\n#                         double *a, lapack_int lda,<\/span>\r\n#                         lapack_int *ipiv,<\/span>\r\n#                         double *b, lapack_int ldb)<\/span>\r\ndef<\/span> sls_lu_orig<\/span>(<\/span>X<\/span>,<\/span>Y<\/span>):<\/span>\r\n    m<\/span>,<\/span>n<\/span> =<\/span> X<\/span>.<\/span>shape<\/span>\r\n    nrhs<\/span> =<\/span> Y<\/span>.<\/span>shape<\/span>[<\/span>1<\/span>]<\/span> if<\/span> Y<\/span>.<\/span>ndim<\/span> ==<\/span> 2<\/span> else<\/span> 1<\/span>\r\n\r\n    workspace<\/span> =<\/span> ffi<\/span>.<\/span>new<\/span>(<\/span>"int[]"<\/span>,<\/span> n<\/span>)<\/span> # workspace array<\/span>\r\n\r\n    info<\/span> =<\/span> lap<\/span>.<\/span>LAPACKE_dgesv<\/span>(<\/span>lap<\/span>.<\/span>LAPACK_ROW_MAJOR<\/span>,<\/span>\r\n                             n<\/span>,<\/span> nrhs<\/span>,<\/span>\r\n                             get_dptr<\/span>(<\/span>X<\/span>),<\/span> n<\/span>,<\/span>    # ldX = n<\/span>\r\n                             workspace<\/span>,<\/span>\r\n                             get_dptr<\/span>(<\/span>Y<\/span>),<\/span> nrhs<\/span>)<\/span> # ldY = nrhs<\/span>\r\n\r\n    assert<\/span> not<\/span> info<\/span>\r\n    return<\/span> Y<\/span>[:<\/span>n<\/span>]<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [8]:<\/div>\n
    \n
    \n
    \n
    # %memit sls_orig(X,Y)<\/span>\r\n# peak memory: 576.93 MiB, increment: 6.71 MiB<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    QR via dgels<\/code><\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [9]:<\/div>\n
    \n
    \n
    \n
    # gels<\/span>\r\n# http:\/\/www.netlib.org\/lapack\/explore-html\/d9\/da2\/lapacke__dgels_8c.html<\/span>\r\n# https:\/\/software.intel.com\/en-us\/node\/521112<\/span>\r\n#<\/span>\r\n#<\/span>\r\n#lapack_int LAPACKE_dgels(int matrix_layout,<\/span>\r\n#                         char trans,<\/span>\r\n#                         lapack_int m, lapack_int n, lapack_int nrhs,<\/span>\r\n#                         double* a, lapack_int lda,<\/span>\r\n#                         double* b, lapack_int ldb);<\/span>\r\ndef<\/span> sls_qr_orig<\/span>(<\/span>X<\/span>,<\/span>Y<\/span>):<\/span>\r\n    m<\/span>,<\/span>n<\/span>  =<\/span> X<\/span>.<\/span>shape<\/span>\r\n    nrhs<\/span> =<\/span> Y<\/span>.<\/span>shape<\/span>[<\/span>1<\/span>]<\/span> if<\/span> Y<\/span>.<\/span>ndim<\/span> ==<\/span> 2<\/span> else<\/span> 1<\/span>\r\n    info<\/span> =<\/span> lap<\/span>.<\/span>LAPACKE_dgels<\/span>(<\/span>lap<\/span>.<\/span>LAPACK_ROW_MAJOR<\/span>,<\/span>\r\n                             "N"<\/span>,<\/span> # do not transpose<\/span>\r\n                             m<\/span>,<\/span> n<\/span>,<\/span> nrhs<\/span>,<\/span>\r\n                             get_dptr<\/span>(<\/span>X<\/span>),<\/span> n<\/span>,<\/span>\r\n                             get_dptr<\/span>(<\/span>Y<\/span>),<\/span> nrhs<\/span>)<\/span>\r\n\r\n    assert<\/span> not<\/span> info<\/span>\r\n    return<\/span> Y<\/span>[:<\/span>n<\/span>]<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [10]:<\/div>\n
    \n
    \n
    \n
    # %memit sls_qr_orig(X,Y)<\/span>\r\n# peak memory: 1083.94 MiB, increment: 514.06 MiB<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    QR via dgels<\/code>: Take Two<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [11]:<\/div>\n
    \n
    \n
    \n
    def<\/span> sls_qr_minwork<\/span>(<\/span>X<\/span>,<\/span>Y<\/span>):<\/span>\r\n    m<\/span>,<\/span>n<\/span>  =<\/span> X<\/span>.<\/span>shape<\/span>\r\n    nrhs<\/span> =<\/span> Y<\/span>.<\/span>shape<\/span>[<\/span>1<\/span>]<\/span> if<\/span> Y<\/span>.<\/span>ndim<\/span> ==<\/span> 2<\/span> else<\/span> 1<\/span>\r\n    \r\n\r\n    # http:\/\/www.netlib.org\/lapack\/explore-html\/d8\/dde\/dgels_8f_source.html<\/span>\r\n    # min LWORK >= max( 1, MN + max( MN, NRHS ) )<\/span>\r\n    work_size<\/span> =<\/span> X<\/span>.<\/span>size<\/span> +<\/span> max<\/span>(<\/span>X<\/span>.<\/span>size<\/span>,<\/span>nrhs<\/span>)<\/span>\r\n    work<\/span> =<\/span> np<\/span>.<\/span>zeros<\/span>(<\/span>work_size<\/span>)<\/span>\r\n    \r\n    info<\/span> =<\/span> lap<\/span>.<\/span>LAPACKE_dgels_work<\/span>(<\/span>lap<\/span>.<\/span>LAPACK_ROW_MAJOR<\/span>,<\/span>\r\n                                 "N"<\/span>,<\/span> # do not transpose<\/span>\r\n                                 m<\/span>,<\/span> n<\/span>,<\/span> nrhs<\/span>,<\/span>\r\n                                 get_dptr<\/span>(<\/span>X<\/span>),<\/span> n<\/span>,<\/span>\r\n                                 get_dptr<\/span>(<\/span>Y<\/span>),<\/span> nrhs<\/span>,<\/span>\r\n                                 get_dptr<\/span>(<\/span>work<\/span>),<\/span> work_size<\/span>)<\/span>\r\n\r\n    assert<\/span> not<\/span> info<\/span>\r\n    return<\/span> Y<\/span>[:<\/span>n<\/span>]<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [12]:<\/div>\n
    \n
    \n
    \n
    # %memit sls_qr_minwork(X,Y)<\/span>\r\n#peak memory: 1084.91 MiB, increment: 517.84 MiB<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    SVD via dgelss<\/code><\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [13]:<\/div>\n
    \n
    \n
    \n
    # gelss<\/span>\r\n#<\/span>\r\n#<\/span>\r\n# http:\/\/www.netlib.org\/lapack\/explore-html\/db\/d77\/lapacke__dgelss_8c.html<\/span>\r\n# https:\/\/software.intel.com\/en-us\/node\/521114<\/span>\r\n#<\/span>\r\n#<\/span>\r\n#lapack_int LAPACKE_dgelss(int matrix_layout, <\/span>\r\n#                          lapack_int m, lapack_int n, lapack_int nrhs, <\/span>\r\n#                          double* a, lapack_int lda, <\/span>\r\n#                          double* b, lapack_int ldb, <\/span>\r\n#                          double* s, <\/span>\r\n#                          double rcond, <\/span>\r\n#                          lapack_int* rank )<\/span>\r\n#<\/span>\r\ndef<\/span> sls_svd_dgelss<\/span>(<\/span>X<\/span>,<\/span>Y<\/span>):<\/span>\r\n    m<\/span>,<\/span>n<\/span> =<\/span> X<\/span>.<\/span>shape<\/span>\r\n    nrhs<\/span> =<\/span> Y<\/span>.<\/span>shape<\/span>[<\/span>1<\/span>]<\/span> if<\/span> Y<\/span>.<\/span>ndim<\/span> ==<\/span> 2<\/span> else<\/span> 1<\/span>\r\n\r\n    singular_values<\/span> =<\/span> np<\/span>.<\/span>empty<\/span>(<\/span>min<\/span>(<\/span>m<\/span>,<\/span>n<\/span>))<\/span>\r\n    rank<\/span> =<\/span> np<\/span>.<\/span>zeros<\/span>(<\/span>1<\/span>,<\/span> dtype<\/span>=<\/span>np<\/span>.<\/span>int64<\/span>)<\/span>\r\n\r\n    info<\/span> =<\/span> lap<\/span>.<\/span>LAPACKE_dgelss<\/span>(<\/span>lap<\/span>.<\/span>LAPACK_ROW_MAJOR<\/span>,<\/span>\r\n                              m<\/span>,<\/span>n<\/span>,<\/span>nrhs<\/span>,<\/span>\r\n                              get_dptr<\/span>(<\/span>X<\/span>),<\/span> n<\/span>,<\/span>\r\n                              get_dptr<\/span>(<\/span>Y<\/span>),<\/span> nrhs<\/span>,<\/span>\r\n                              get_dptr<\/span>(<\/span>singular_values<\/span>),<\/span>\r\n                              -<\/span>1<\/span>,<\/span> \r\n                              get_iptr<\/span>(<\/span>rank<\/span>))<\/span>\r\n    return<\/span> Y<\/span>[:<\/span>n<\/span>]<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [14]:<\/div>\n
    \n
    \n
    \n
    # %memit sls_svd_dgelss(X,Y)<\/span>\r\n# peak memory: 1083.97 MiB, increment: 513.99 MiB<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    SVD via dgelsd<\/code><\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [15]:<\/div>\n
    \n
    \n
    \n
    # gelsd<\/span>\r\n#<\/span>\r\n# http:\/\/www.netlib.org\/lapack\/explore-html\/df\/de3\/lapacke__dgelsd_8c.html<\/span>\r\n# https:\/\/software.intel.com\/en-us\/node\/521115<\/span>\r\n#<\/span>\r\n#<\/span>\r\n#lapack_int LAPACKE_dgelsd(int matrix_layout,<\/span>\r\n#                          lapack_int m, lapack_int n, lapack_int nrhs,<\/span>\r\n#                          double* a, lapack_int lda,<\/span>\r\n#                          double* b, lapack_int ldb,<\/span>\r\n#                          double* s,<\/span>\r\n#                          double rcond,<\/span>\r\n#                          lapack_int* rank );<\/span>\r\n#<\/span>\r\ndef<\/span> sls_svd_orig<\/span>(<\/span>X<\/span>,<\/span>Y<\/span>):<\/span>\r\n    m<\/span>,<\/span>n<\/span> =<\/span> X<\/span>.<\/span>shape<\/span>\r\n    nrhs<\/span> =<\/span> Y<\/span>.<\/span>shape<\/span>[<\/span>1<\/span>]<\/span> if<\/span> Y<\/span>.<\/span>ndim<\/span> ==<\/span> 2<\/span> else<\/span> 1<\/span>\r\n\r\n    singular_values<\/span> =<\/span> np<\/span>.<\/span>empty<\/span>(<\/span>min<\/span>(<\/span>m<\/span>,<\/span>n<\/span>))<\/span>\r\n    rank<\/span> =<\/span> np<\/span>.<\/span>zeros<\/span>(<\/span>1<\/span>,<\/span> dtype<\/span>=<\/span>np<\/span>.<\/span>int64<\/span>)<\/span>\r\n\r\n    info<\/span> =<\/span> lap<\/span>.<\/span>LAPACKE_dgelsd<\/span>(<\/span>lap<\/span>.<\/span>LAPACK_ROW_MAJOR<\/span>,<\/span>\r\n                              m<\/span>,<\/span>n<\/span>,<\/span>nrhs<\/span>,<\/span>\r\n                              get_dptr<\/span>(<\/span>X<\/span>),<\/span> n<\/span>,<\/span>\r\n                              get_dptr<\/span>(<\/span>Y<\/span>),<\/span> nrhs<\/span>,<\/span>\r\n                              get_dptr<\/span>(<\/span>singular_values<\/span>),<\/span>\r\n                              -<\/span>1<\/span>,<\/span> \r\n                              get_iptr<\/span>(<\/span>rank<\/span>))<\/span>\r\n    return<\/span> Y<\/span>[:<\/span>n<\/span>]<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [16]:<\/div>\n
    \n
    \n
    \n
    # %memit sls_svd_orig(X,Y)<\/span>\r\n#peak memory: 1084.62 MiB, increment: 514.18 MiB<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Check for Correctness<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Here, we test the code to make sure they all give the same \\(\\beta\\)<\/span>s. Note that we make a copy because most of these routines modify X<\/code> adn Y<\/code> in-place. This differs a bit from our semantics of early posts, but it isn’t terribly important because within this<\/em> post, we are still comparing apples to apples. When we go back to the outside world, we may still need to make copies, depending on our use-cases for X<\/code> and Y<\/code>.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [17]:<\/div>\n
    \n
    \n
    \n
    _X<\/span>,<\/span>_Y<\/span> =<\/span> create_data<\/span>(<\/span>n_obs<\/span>,<\/span>n_ftrs<\/span>)<\/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\nbetas_chol<\/span> =<\/span> sls_cholesky_orig<\/span>(<\/span>X<\/span>,<\/span>Y<\/span>)<\/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\nbetas_qr<\/span> =<\/span> sls_qr_orig<\/span>(<\/span>X<\/span>,<\/span>Y<\/span>)<\/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\nbetas_svd<\/span> =<\/span> sls_svd_orig<\/span>(<\/span>X<\/span>,<\/span>Y<\/span>)<\/span>\r\n\r\nbetas<\/span> =<\/span> [<\/span>betas_chol<\/span>,<\/span> betas_qr<\/span>,<\/span> betas_svd<\/span>]<\/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
    And Do a Quick Timing Run<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [18]:<\/div>\n
    \n
    \n
    \n
    %<\/span>timeit<\/span> sls_cholesky_orig(X,Y)\r\n%<\/span>timeit<\/span> sls_lu_orig(X,Y)\r\n%<\/span>timeit<\/span> sls_qr_orig(X,Y)\r\n%<\/span>timeit<\/span> sls_qr_minwork(X,Y)\r\n%<\/span>timeit<\/span> sls_svd_dgelss(X,Y)\r\n%<\/span>timeit<\/span> sls_svd_orig(X,Y)\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    \n
    <\/div>\n
    \n
    1 loops, best of 3: 1.54 s per loop\r\n10 loops, best of 3: 32.6 ms per loop\r\n1 loops, best of 3: 8.99 s per loop\r\n1 loops, best of 3: 8.77 s per loop\r\n1 loops, best of 3: 9.17 s per loop\r\n1 loops, best of 3: 9.19 s per loop\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Summary Results<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n\n\n\n\n\n\n\n\n\n\n
    Algorithm<\/th>\nCalled Via<\/th>\nMemory
    Increment(MB)<\/th>\n    		~Time(s)<\/th>\n<\/tr>\n<\/thead>\n
    Cholesky<\/td>\ndposv<\/code><\/td>\n12<\/td>\n1.5<\/td>\n<\/tr>\n
    LU<\/td>\ndgesv<\/code><\/td>\n7<\/td>\n>0.0<\/td>\n<\/tr>\n
    QR<\/td>\ndgels<\/code><\/td>\n514<\/td>\n9<\/td>\n<\/tr>\n
    QR (min LWORK)<\/td>\ndgels<\/code><\/td>\n517<\/td>\n9<\/td>\n<\/tr>\n
    SVD<\/td>\ndgelss<\/code><\/td>\n514<\/td>\n9<\/td>\n<\/tr>\n
    SVD<\/td>\ndgelsd<\/code><\/td>\n514<\/td>\n9<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    As you can see, things are still not "right". And there are some hints running around in the code and discussion above as to what the remaining issue is. I won’t spill the beans, but stay tuned for next time.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    \nAdditional Resources
    \n<\/h3>\n
    \nYou can grab a copy of this notebook<\/a>.\n<\/p>\n
    \nEven better, you can view it using nbviewer<\/a>.\n<\/p>\n
    \nLicense
    \n<\/h3>\n
    \nUnless otherwise noted, the contents of this notebook are under the following license. The code in the notebook should be considered part of the text (i.e., licensed and treated as as follows).\n<\/p>\n
    \"Creative<\/a>
    DrsFenner.org Blog And Notebooks<\/span> by     Mark and Barbara Fenner<\/a> is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License<\/a>.
    Permissions beyond the scope of this license may be available at     drsfenner.org\/blog\/about-and-contacts<\/a>.<\/p>\n<\/div>\n<\/div>\n<\/div>\n","protected":false},"excerpt":{"rendered":"
    In the last post, we looked at three ways of solving the least-squares regression problem (Cholesky, QR, and SVD decompositions) using LAPACK routines in a round-about fashion. The different decompositions (i.e., factorizations) were performed by (1) calling a NumPy np.linalg routine that called (2) an LAPACK driver routine for a "bigger" problem like linear equation […]<\/p>\n","protected":false},"author":3,"featured_media":0,"comment_status":"closed","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[6,7],"tags":[],"class_list":["post-390","post","type-post","status-publish","format-standard","hentry","category-mrdr","category-sci-math-stat-python"],"_links":{"self":[{"href":"http:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/posts\/390","targetHints":{"allow":["GET"]}}],"collection":[{"href":"http:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"http:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"http:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/users\/3"}],"replies":[{"embeddable":true,"href":"http:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/comments?post=390"}],"version-history":[{"count":2,"href":"http:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/posts\/390\/revisions"}],"predecessor-version":[{"id":427,"href":"http:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/posts\/390\/revisions\/427"}],"wp:attachment":[{"href":"http:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/media?parent=390"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"http:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/categories?post=390"},{"taxonomy":"post_tag","embeddable":true,"href":"http:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/tags?post=390"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}