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).