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


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


$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


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


This can be done as

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

  • $\begingroup$ Where does $\| u \|_{2}^{2} + \| v \|_{2}^{2}= \left[ \begin{array}{c} u \\ v \end{array} \right]_{2}^{2}$ come from? $\endgroup$ – Richard Feb 13 at 18:20
  • $\begingroup$ $\| 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}$. $\endgroup$ – Brian Borchers Feb 13 at 18:22
  • $\begingroup$ Shoot, saw it a moment later. Thanks! $\endgroup$ – Richard Feb 13 at 18:24

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.