numerical solution of an under-determined linear equation in high dimensions

Question

I need to solve a linear regression problem $$Ax=y$$ which is hugely underdetermined. I have around $10^6$ features but only $10^3$ equations. So $A$ is a $1,000\times 1,000,000$ matrix and $y$ a vector of length $1,000$, both are given and I need to find $x$ (of length $1,000,000$).

The solution is of course $x^*=A^\dagger y$ where $A^\dagger$ is the pseudo inverse. $x^*$ minimizes the least-square expression $(Ax-y)^T(Ax-y)$.

All this is well known and I can do very efficiently using standard libraries (lstsq of numpy for example). The problem is that I want to add a regularization term of the form

$$\min_x\Big\{ (Ax-y)^T(Ax-y)+x^T R x\Big\}$$ where $R$ is some $10^6\times 10^6$ matrix that I have (it is of course very sparse). The analytic solution this equation is

$$x=\left(A^TA+R\right)^\dagger y\ ,$$ but it is completely impractical to even instantiate the matrix $A^TA+R$ (a $10^6\times10^6$ matrix). What's the best way to obtain the solution in a numerically stable manner? If that helps, I know how to express $R$ as $R=B^TB$ with b is a matrix of the same shape as $A$.

--- Preemtive apology:: I realize that this might be a very basic question, but I can't find a standard way to do it in any of the packages I use (scikit and the likes). I also know that when $R$ is the identity this is the usual Tikhonov regularization, and there is a trick there of evaluating $A A^T$ instead of $A^TA$, but I'm not sure if that trick applies for arbitrary $R$.

Answer 1

You want to minimize

$\min \| Ax -y \|_{2}^{2} + x^{T}B^{T}Bx=\| Ax -y \|_{2}^{2} + \| Bx \|_{2}^{2}$

Recall that

$\| u \|_{2}^{2} + \| v \|_{2}^{2}= \left\| \left[ \begin{array}{c} u \\ v \end{array} \right] \right\|_{2}^{2}$.

Thus your problem can be written as

$\min \| Hx - g \|_{2}^{2}$

where

$H=\left[ \begin{array}{c} A \\ B \end{array} \right]$

and

$g=\left[ \begin{array}{c} y \\ 0 \end{array} \right]$.

In using an iterative algorithm such as LSQR to solve the least squares problem, you'll need to be able to compute products of the form

$w=Hv$.

This can be done as

$w=\left[ \begin{array}{c} Av \\ Bv \end{array} \right]$.

Similarly, you'll need to be able to compute products of the form

$v=H^{T}v$.

This can be done as

$v=\left[ A^{T} B^{T} \right]w$.