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

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

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

• Where does $\| u \|_{2}^{2} + \| v \|_{2}^{2}= \left[ \begin{array}{c} u \\ v \end{array} \right]_{2}^{2}$ come from? – Richard Feb 13 at 18:20
• $\| u \|_{2}^{2}+ \| v \|_{2}^{2}=u_{1}^{2} + u_{2}^{2} + \ldots u_{n}^{2} + v_{1}^{2}+\ldots + v_{m}^{2}=\| \left[ \begin{array}{c} u \\ v \end{array} \right] \|_{2}^{2}$. – Brian Borchers Feb 13 at 18:22