{"id":388,"date":"2015-12-11T15:54:10","date_gmt":"2015-12-11T15:54:10","guid":{"rendered":"http:\/\/drsfenner.org\/blog\/?p=388"},"modified":"2016-02-05T16:28:27","modified_gmt":"2016-02-05T16:28:27","slug":"three-paths-to-least-squares-linear-regression","status":"publish","type":"post","link":"https:\/\/drsfenner.org\/blog\/2015\/12\/three-paths-to-least-squares-linear-regression\/","title":{"rendered":"Three Paths to Least Squares Linear Regression"},"content":{"rendered":"
\n
\n<\/div>\n
\n
\n

In the prior installment<\/a>, we talked about solving the basic problem of linear regression using a standard least-squares approach. Along the way, I showed a few methods that gave the same results to the least squares problem. Today, I want to look at three methods for solving linear regression:<\/p>\n

    \n
  1. solving the full normal equations via Cholesky decomposition<\/li>\n
  2. solving a least-squares problem via QR decomposition<\/li>\n
  3. solving a least-squares problem via SVD (singular value decomposition)<\/li>\n<\/ol>\n

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

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

    The following table may not make complete sense to start with, but consider it an outline of what we are going to discuss. nla<\/code>\/sla<\/code>\/cla<\/code> refer to numpy.linalg<\/code>, scipy.linalg<\/code>, and cvxopt.lapack<\/code>, respectively. Each of the top-level routines makes use of high quality LAPACK<\/a> routines under the hood.<\/p>\n

    While we’ll use abbreviations for the LAPACK routine names, these abbreviated names don’t really exist in LAPACK. For example, there is no posv<\/code> in LAPACK. posv<\/code> is our shorthand for dposv<\/code> which is the double<\/code> (double precision floating-point) specific version of _posv<\/code>. There is also sposv<\/code>, cposv<\/code>, and zposv<\/code>. For example, see the pair of rows labelled "symmetric\/Hermiatian positive definite" (the two row headers belong together, it’s just bad formatting) in Table 2.2 of the LUG<\/a>.<\/p>\n

    The LAPACK routines for these solutions are described generally under (Solving) Linear Equations<\/a> and Linear Least Squares<\/a>. The first set of routines are for square matrices. In some cases, the Intel documentation for the MKL (which includes an LAPACK implementation) has better descriptions of some of these routines. See: (Solving) Linear Equations<\/a> and Linear Least Squares<\/a>.<\/p>\n<\/div>\n<\/div>\n<\/div>\n

    \n
    \n<\/div>\n
    \n
    \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<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n

    Imports<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n

    Not too surprising, but we’ll do some imports to get us started. If you haven’t seen it before, cvxopt<\/code> is a nice package for coding optimization problems. Because these problems often involve linear algebra, cvxopt<\/code> has an (partial) interface to LAPACK and it exposes some routines that numpy<\/code> and scipy<\/code>‘s LAPACK interfaces miss. Annoyingly, none of these three LAPACK interfaces are complete. Even if you put them all together, there is still missing functionality. But, we’ll deal with a workaround next time.<\/p>\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<\/div>\n
    \n
    \n
    Data Creation and a Memory Profiling Note<\/h3>\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<\/div>\n
    \n
    \n
    You’ll notice that I used an iPython magic (%load_ext memory_profiler<\/code>) above in the import<\/code> section. I’m not going to jump into the details of installing the memory profiler (but pip install memory_profiler<\/code> might be all you need). However, we do have to discuss a few issues that come up with its use. If you want to recreate the memory use values I gather below, you’ll need to be a bit careful.<\/p>\n
    If I say:<\/p>\n
        X,Y = create_data()\r\n    %memit X.T.dot(X)<\/code><\/pre>\n
    and look at the result, I might see a memory increment of ~8MB. That is: the process of creating X.T.dot(X)<\/code> required an additional 8MB of memory. But, if I then go and run the dot-product again in the next cell:<\/p>\n
        %memit X.T.dot(X)<\/code><\/pre>\n
    it may report an increment of ~0MB. What happened?<\/p>\n
    The first X.T.dot(X)<\/code> keeps no references to the result and the result can be free()<\/code>‘d – garbage collected – back to the OS. However, the OS might not automatically reclaim that memory and\/or Python’s garbage collector might not automatically release it (that’s a completely post worthy topic of its own). Thus, the amount of memory that our process is using is still at the +8MB level. So, when we run another dot-product, our process uses that unused-by-a-variable, but not reclaimed-by-the-OS, memory and we see no increment.<\/p>\n
    The moral of that story is that if you want to see consistent increments, you need to either (1) keep a reference to the resulting object (XtX = X.T.dot(X)<\/code>) so the memory is both not freed and not reclaimed (i.e., it is in use) or (2) after running the %memit<\/code> line, restart your kernel before you do the next %memit<\/code> (foring a freed and reclaimed state for -all- of the process memory).<\/p>\n
    The problem is even a little more severe: any<\/em> large computation (not just memit<\/code> expressions) that is evaluated and not retained via a reference will cause this memory use masking — and it makes perfect sense when you think about it. So, all of my memit<\/code>s will be done in a fresh interpreter and I’ll simply leave comments about the results.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Create Some Data and Investigate Memory Use<\/h3>\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\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    \n
    <\/div>\n
    \n
    about 10000000 entries\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [4]:<\/div>\n
    \n
    \n
    \n
    print<\/span> X<\/span>.<\/span>shape<\/span>\r\nprint<\/span> X<\/span>.<\/span>itemsize<\/span>,<\/span> X<\/span>.<\/span>size<\/span>,<\/span> reduce<\/span>(<\/span>op<\/span>.<\/span>mul<\/span>,<\/span> X<\/span>.<\/span>shape<\/span>)<\/span>\r\n# memory(X) : 2**(16 + 10 + 3) [2**3 (8 byte)]<\/span>\r\nprint<\/span> (<\/span>X<\/span>.<\/span>itemsize<\/span> *<\/span> X<\/span>.<\/span>size<\/span>)<\/span> \/<\/span> float<\/span>(<\/span>1024<\/span>**<\/span>2<\/span>)<\/span>  # 1024**2 == 2**20 == 1 MiB<\/span>\r\nprint<\/span> 2<\/span>**<\/span>(<\/span>16<\/span> +<\/span> 10<\/span> +<\/span> 3<\/span>)<\/span> \/<\/span> 2<\/span>**<\/span>20<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    \n
    <\/div>\n
    \n
    (65536, 1024)\r\n8 67108864 67108864\r\n512.0\r\n512\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [5]:<\/div>\n
    \n
    \n
    \n
    print<\/span> X<\/span>.<\/span>shape<\/span>\r\nXtX<\/span> =<\/span> X<\/span>.<\/span>T<\/span>.<\/span>dot<\/span>(<\/span>X<\/span>)<\/span>\r\nprint<\/span> XtX<\/span>.<\/span>shape<\/span>\r\nprint<\/span> (<\/span>XtX<\/span>.<\/span>itemsize<\/span> *<\/span> XtX<\/span>.<\/span>size<\/span>)<\/span> \/<\/span> float<\/span>(<\/span>1024<\/span>**<\/span>2<\/span>)<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    \n
    <\/div>\n
    \n
    (65536, 1024)\r\n(1024, 1024)\r\n8.0\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    np.dot()<\/code> Tests<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Let’s do two simple tests to get an idea of (1) the cost of some basic operations and (2) understand the portion of resources used in the Cholesky that are not<\/em> due to the factorization and backwards solution.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [6]:<\/div>\n
    \n
    \n
    \n
    #%memit X.T.dot(X)<\/span>\r\n#increment ~27.5MB<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    It’s possibly surprising that \\(X^TX\\)<\/span> is so much smaller than \\(X\\)<\/span>. However, we’re dealing with a "tall-skinny" matrix — there are many more rows than colums (\\(n_{obs} >> n_{ftrs}\\)<\/span>). So, when we take the transpose and multiply we have a \\((n_{ftr},n_{obs}) \\times (n_{obs}, n_{ftr}) \\rightarrow (n_{ftr},n_{ftr})\\)<\/span> which has far fewer elements than \\((n_{obs}, n_{ftr})\\)<\/span>.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [7]:<\/div>\n
    \n
    \n
    \n
    %<\/span>timeit<\/span> X.T.dot(X)\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    \n
    <\/div>\n
    \n
    1 loops, best of 3: 1.39 s per loop\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Compressing data requires a trade-off in pre-processing. Here, we see that the dot-product takes a fair bit of time. Remember that it is doing its work (a vector multiplication) in the long (\\(n_{obs}\\)<\/span>\/column) dimension.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Via Cholesky<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Recall that that solving linear regression using least squares means that we want to perform the following optimization: \\(\\min RSS(\\beta)=\\min(\\bf{y}-\\bf{X}\\beta)^T(\\bf{y}-\\bf{X}\\beta)\\)<\/span>. That means that we want: \\(\\frac{\\partial RSS(\\beta)}{\\partial \\beta} = -2 X^T(Y-X\\beta) = 0\\)<\/span>. We saw that this gives: \\(X^TX\\beta = X^TY\\)<\/span>, which we want to solve for \\(\\beta\\)<\/span>.<\/p>\n
    \\(Cholesky(X^TX)\\)<\/span> gives \\(X^TX=LL^T\\)<\/span> with \\(L\\)<\/span> lower triangular.<\/p>\n
    Then, we have \\(L^TL\\beta = X^TY\\)<\/span>. We let \\(w=L\\beta\\)<\/span> and solve \\(Lw=X^TY\\)<\/span> for \\(w\\)<\/span> using simple backward substitution. Now, we return to \\(L\\beta=w\\)<\/span> and we solve it for \\(\\beta\\)<\/span> using simple forward substitution.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [8]:<\/div>\n
    \n
    \n
    \n
    #%%memit<\/span>\r\n#XtX, XtY = X.T.dot(X), X.T.dot(Y)<\/span>\r\n#betas = sla.solve(XtX,XtY, sym_pos=True, check_finite=False)<\/span>\r\n# increment ~28MB<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [9]:<\/div>\n
    \n
    \n
    \n
    %%<\/span>timeit<\/span> \r\nXtX, XtY = X.T.dot(X), X.T.dot(Y)\r\nbetas = sla.solve(XtX,XtY, sym_pos=True, check_finite=False)\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    \n
    <\/div>\n
    \n
    1 loops, best of 3: 1.45 s per loop\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    For comparison sake, here’s the same effective call but using cvxopt<\/code>‘s LAPACK interface to access posv<\/code> directly. Because cla.posv<\/code> is so low-level it has the same "modify the input" semantics of LAPACK’s posv<\/code>. Thus, if we want to keep the comparison between the two methods (sla.solve(sym_pos=True)<\/code> and cla.posv<\/code>) fair, we have to explicitly make our own copy. There’s a second reason we have to make that copy: the cvxopt.lapack<\/code> interface expects cvxopt.matrix<\/code> objects.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [10]:<\/div>\n
    \n
    \n
    \n
    %%<\/span>timeit<\/span>\r\nXtX, XtY = X.T.dot(X), X.T.dot(Y)\r\ncXtX, cXtY = cvx_matrix(XtX), cvx_matrix(XtY)\r\ncla.posv(cXtX, cXtY)\r\nbetas = cXtY[:n_ftrs]\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    \n
    <\/div>\n
    \n
    1 loops, best of 3: 1.46 s per loop\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Via QR<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    The QR decomposition let’s us write: \\(X=QR\\)<\/span> so that \\(X\\beta=Y\\)<\/span> becomes \\(QR\\beta=Y\\)<\/span>. The QR decomposition gives a \\(Q\\)<\/span> that is orthogonal which means \\(Q^T=Q^{-1}\\)<\/span>. So, we can multiply on the left by \\(Q^{-1}\\)<\/span>: \\(R \\beta = Q^TY\\)<\/span>. Lastly, \\(R\\)<\/span> is upper triangular so it is efficient and stable to use backwards substitution to solve that system for \\(\\beta\\)<\/span> giving (the numerically stable equivalent to): \\(\\beta = R^{-1}Q^TY\\)<\/span>.<\/p>\n
    Here, we will use cla.gels<\/code>. As with cla.posv<\/code>, we need to make our own copy to keep both the result and the input matrix available.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [11]:<\/div>\n
    \n
    \n
    \n
    #%%memit<\/span>\r\n#cX, cY = cvx_matrix(X),cvx_matrix(Y)  # makes a copy<\/span>\r\n# increment:  ~512MB<\/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<\/span>\r\n#cla.gels(cX,cY) # betas are in cY[:p]<\/span>\r\n# increment: ~28 MiB<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [13]:<\/div>\n
    \n
    \n
    \n
    %%<\/span>timeit<\/span>\r\ncX, cY = cvx_matrix(X),cvx_matrix(Y)  # makes a copy\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    \n
    <\/div>\n
    \n
    1 loops, best of 3: 954 ms per loop\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [14]:<\/div>\n
    \n
    \n
    \n
    %%<\/span>timeit<\/span>\r\ncX, cY = cvx_matrix(X),cvx_matrix(Y)  # makes a copy\r\ncla.gels(cX,cY) # betas are in cY[:p]\r\nbetas = cY[:n_ftrs]\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    \n
    <\/div>\n
    \n
    1 loops, best of 3: 8.23 s per loop\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Via SVD<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Performing linear regression with SVD is quite easy – even if the theory behind it is slightly more complex. If \\(X=U \\Sigma V\\)<\/span>, then the pseudo-inverse of \\(X\\)<\/span> is \\(X^{\\dagger} = V {\\Sigma}^{\\dagger} U^T\\)<\/span> and we can solve \\(X\\beta=Y\\)<\/span> with \\(\\beta = X^{\\dagger} X \\beta = X^\\dagger Y = V{\\Sigma}^{\\dagger} U^T Y\\)<\/span><\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [15]:<\/div>\n
    \n
    \n
    \n
    #%%memit <\/span>\r\n#betas = nla.lstsq(X,Y)[0]<\/span>\r\n# increment ~540MB<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [16]:<\/div>\n
    \n
    \n
    \n
    %%<\/span>timeit<\/span>\r\nbetas = nla.lstsq(X,Y)[0]\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    \n
    <\/div>\n
    \n
    1 loops, best of 3: 8.16 s per loop\r\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Comparing Betas<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Perhaps we should have done this first, but we should make certain that all of our methods for producing the \\(\\beta\\)<\/span> coefficients give us the same results. Fortunately, it’s easy enough. You might wonder why we don’t just use the values we computed above. Well, like %memit<\/code>, %timeit<\/code> has quirks. And one of those quirks is that variables defined in a %%timeit<\/code> cell don’t leak out into the global environment (they are local to the %timeit<\/code> execution). Darn.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [17]:<\/div>\n
    \n
    \n
    \n
    #<\/span>\r\n# Cholesky <\/span>\r\n#<\/span>\r\nXtX<\/span>,<\/span> XtY<\/span> =<\/span> X<\/span>.<\/span>T<\/span>.<\/span>dot<\/span>(<\/span>X<\/span>),<\/span> X<\/span>.<\/span>T<\/span>.<\/span>dot<\/span>(<\/span>Y<\/span>)<\/span>\r\nbetas_chol_1<\/span> =<\/span> sla<\/span>.<\/span>solve<\/span>(<\/span>XtX<\/span>,<\/span>XtY<\/span>,<\/span> sym_pos<\/span>=<\/span>True<\/span>,<\/span> check_finite<\/span>=<\/span>False<\/span>)<\/span>\r\n\r\ncXtX<\/span>,<\/span> cXtY<\/span> =<\/span> cvx_matrix<\/span>(<\/span>XtX<\/span>),<\/span> cvx_matrix<\/span>(<\/span>XtY<\/span>)<\/span>\r\ncla<\/span>.<\/span>posv<\/span>(<\/span>cXtX<\/span>,<\/span> cXtY<\/span>)<\/span>\r\nbetas_chol_2<\/span> =<\/span> cXtY<\/span>[:<\/span>n_ftrs<\/span>]<\/span>\r\n\r\n\r\n#<\/span>\r\n# QR<\/span>\r\n#<\/span>\r\ncX<\/span>,<\/span> cY<\/span> =<\/span> cvx_matrix<\/span>(<\/span>X<\/span>),<\/span>cvx_matrix<\/span>(<\/span>Y<\/span>)<\/span>  # makes a copy<\/span>\r\ncla<\/span>.<\/span>gels<\/span>(<\/span>cX<\/span>,<\/span>cY<\/span>)<\/span> # betas are in cY[:p]<\/span>\r\nbetas_qr<\/span> =<\/span> cY<\/span>[:<\/span>n_ftrs<\/span>]<\/span>\r\n\r\n#<\/span>\r\n# SVD<\/span>\r\n#<\/span>\r\nbetas_svd<\/span> =<\/span> nla<\/span>.<\/span>lstsq<\/span>(<\/span>X<\/span>,<\/span>Y<\/span>)[<\/span>0<\/span>]<\/span>\r\n\r\nbetas<\/span> =<\/span> [<\/span>betas_chol_1<\/span>,<\/span> betas_chol_2<\/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
    Summary Results<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    For an np.array<\/code> of ~10M entries (and 512MB) and with non-overwriting semantics, we have:<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\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
    dot-product<\/td>\nnp.dot()<\/code><\/td>\n28<\/td>\n1.5<\/td>\n<\/tr>\n
    Cholesky<\/td>\nsla.solve(sym_pos=True)<\/code>
    cla.posv()<\/code><\/td>\n    		28<\/td>\n1.5<\/td>\n<\/tr>\n
    QR<\/td>\ncla.gels()<\/code><\/td>\n540<\/td>\n8<\/td>\n<\/tr>\n
    SVD<\/td>\nnla.lstsq()<\/code><\/td>\n540<\/td>\n8<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Now, if you’ve read about the Cholesky, QR, and SVD methods as applied to linear regression, you might be a bit concerned. You probably read that in-order they go from fastest to slowest and from least to most numerically stable. We haven’t talked at all about the stability issue, but we don’t see any difference between QR and SVD here. We also see that both QR and SVD are making full copies of the X<\/code> matrix. In the case of QR, we saw that explicitly (we made the copy) because cla.gels<\/code> is low enough level (close enough to the LAPACK gels<\/code> call) that it overwrites its inputs.<\/p>\n
    We’ll dig into these methods by directly using the LAPACK routines, including a (hopefully) lower overhead gels<\/code> 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 prior installment, we talked about solving the basic problem of linear regression using a standard least-squares approach. Along the way, I showed a few methods that gave the same results to the least squares problem. Today, I want to look at three methods for solving linear regression: solving the full normal equations via […]<\/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-388","post","type-post","status-publish","format-standard","hentry","category-mrdr","category-sci-math-stat-python"],"_links":{"self":[{"href":"https:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/posts\/388","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/users\/3"}],"replies":[{"embeddable":true,"href":"https:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/comments?post=388"}],"version-history":[{"count":2,"href":"https:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/posts\/388\/revisions"}],"predecessor-version":[{"id":426,"href":"https:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/posts\/388\/revisions\/426"}],"wp:attachment":[{"href":"https:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/media?parent=388"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/categories?post=388"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/tags?post=388"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}