Clojure Numerics, Part 6 - More Linear Algebra Fun with Least Squares

Need help with your custom Clojure software? I'm open to (selected) contract work.

December 27, 2017

Please share: .

These books fund my work! Please check them out.

This time, we look at a few important variants of the least squares problem: least squares with equality constraints, and generalized least squares. They show how we can solve problems where there is some norm to be minimized with additional constraints to be satisfied. A very cool and useful thing to know!

Before I continue, a few reminders:

The namespaces we'll use: 

(require '[uncomplicate.neanderthal
           [native :refer [dge dv]]
           [core :refer [mm mv]]
           [linalg :refer [lse gls]]])

Least Squares With Equality Constraints

We have a system \(Ax=c\) that we need to find the least squares solution, as in the last article. But, additionally, we have the constraint \(Bx=d\). There is much theory connected to how to do this in a stable way, but the programmer in me just want to skip it and show you the code. So, here it is. In Clojure, just use the lse function! 

(let [a (dge 6 4 [-0.57 -1.28 -0.39 0.25
                  -1.93 1.08 -0.31 -2.14
                  2.30 0.24 0.40 -0.35
                  -1.93 0.64 -0.66 0.08
                  0.15 0.30 0.15 -2.13
                  -0.02 1.03 -1.43 0.50]
             {:layout :row})
      b (dge 2 4 [1 0 -1 0
                  0 1 0 -1]
             {:layout :row})
      c (dv [-1.50 -2.14 1.23 -0.54 -1.68 0.82])
      d (dv [7 2])]
  (lse a b c d))

#RealBlockVector[double, n:4, offset: 0, stride:1]
[   2.57    1.12   -4.43   -0.88 ]

Our solution is \(x=(2.57, 1.12, -4.43, -0.88)\).

Is it correct? Let's see if it satisfies the constraint: 

(mv (dge 2 4 [1 0 -1 0
              0 1 0 -1]
         {:layout :row})
    (dv [2.57 1.12 -4.43 -0.88 ]))

#RealBlockVector[double, n:2, offset: 0, stride:1]
[   7.00    2.00 ]

Yes!

I won't bother you longer with this, other than repeat that the search key phrase for this kind of problems is Least Sqares with Equality Constraints (LSE) and that the lse and its impure cousin lse! is what you need in those cases.

Generalized Least Squares

Here's another stock problem from the Least Squares family. This one is connected to problems commonly found in statistics, data analysis, machine learning, etc. such as regularization, weighting and similar stuff.

Let's see this case similar to the stock problem solvable by least squares: \(Ax+w=c\). We added \(w\) to the basic equation. It can describe, for example, normally distributed noise in the system, with zero mean, and covariance matirx \(\sigma^2 W\), where \(W=BB^T\). Now, we could have proceeded by doing some matrix algebra and concluding that the solution can be found by minimizing \({\lVert B^{-1}(Ax-c) \rVert}_2\). We know (since we read that in the previous articles in this series! :) that computing and working with inverse matrices is bad, bad, bad: a slow procedure that often gives wrong answers.

It was shown by smart mathematicians that this problem is equivalent to the problem of minimizing the norm \({\lVert y^Ty \rVert}\) with the constraint \(Ax+By=c\). Fair enough, but how does it help us?

That problem is called the Generalized Least Sqares problem (GLS) and it is numerically stable! It is done by a flavor of QR factorization that is too exotic for us to look in detail. The neat thing is that we don't have to, because Neanderthal comes with the function that does everything for us: gls. 

(let [a (dge 4 3 [-0.57 -1.28 -0.39
                  -1.93 1.08 -0.31
                  2.30 0.24 -0.40
                  -0.02 1.03 -1.43]
             {:layout :row})
      b (dge 4 4 [0.5 0 0 0
                  0 1 0 0
                  0 0 2 0
                  0 0 0 5]
             {:layout :row})
      c (dv [1.32 -4 5.52 3.24])]
  (gls a b c))

'(#RealBlockVector(double  n:3  offset: 0  stride:1)
(   1.99   -1.01   -2.99 )
 #RealBlockVector(double  n:4  offset: 0  stride:1)
(  -0.00   -0.00   -0.00    0.01 )
)

The solution is given as a Clojure vector; a pair containing \(x\) and \(y\). In this example, GLS solution is \(x=(1.99, -1.01, -2.99), y=(0, 0, 0, 0.01)\) (roughly; here, only the first two decimals are printed).

Enough of the least squares

Huh. This time I managed to be brief. For better or for worse :)

There are quite a few problems that can be described and solved as one of the variants of least squares. I hope that these typical flavors that I showed can give you enough foundations to find your own way when there's something more exotic to be solved. Don't think that what we'd covered is a small thing, though. Quite a lot of problems could be solved with : ls, lse, and gls. A lot!

Maybe it's time to give some attention to eigen problems now. Remember the Clojure Linear Algebra Refresher (2) - Eigenvalues and Eigenvectors article from the Linear Algebra Refresher series? In the upcoming articles, we'll look into it a bit more seriously.

Clojure Numerics, Part 6 - More Linear Algebra Fun with Least Squares - December 27, 2017 - Dragan Djuric