Using submatrices of matrix decomposition for solving a large number least-squares problems

Question

I want to decrease the computational time for solving a large number (>1000) of least-squares problems. Given a matrix, the system matrix for each least-squares problem is a submatrix of the given matrix. My idea is to precompute a matrix decomposition of the given matrix, then reuse it for each least-squares problem. I need help with choosing a decomposition. Other suggestions of how to speed up the computations are welcome.

To be more precise. Given a real, non-symmetric full matrix $A\in \mathbb{R}^{M\times N}$, where $M > N$ (usually $M\sim 100$ and $N\sim 30$), I have to solve a large number P of least-squares problems of the type $A_{S_{i}}x_{i} = y_{i}$ for $x_{i}$ given $y_{i}$, where $i = 1,\ldots,P$. Here, $A_{S_{i}}$ is a submatrix of $A$, in the meaning that certain rows have been extracted from $A$. Thus $A_{S_{i}}\in \mathbb{R}^{M_{i}\times N}$, where $M_{i} \leq M$.

To speed up the computations my idea is to compute a matrix decomposition of $A$, then reuse it for each least-squares problem. This is done by extracting the corresponding submatrices in the decomposition such that they form $A_{S_{i}}$.

I have tried this with QR-factorization. With $A = QR$, then $A_{S_{i}} = Q_{S_{i}}R$, where $Q_{S_{i}}$ is a submatrix of $Q$. However, $Q_{S_{i}}$ is not necessarily unitary. Thus, it is no longer true that $x_{i} = R \backslash Q_{S_{i}}^{T}y_{i}$, but rather $x_{i} = R \backslash Q_{S_{i}}\backslash y_{i}$. This includes two solves and is, according to my timings, slower than $x_{i} = A_{S_{i}} \backslash y_{i}$.

Is there any matrix decomposition as described above that would speed up my computations? Thank you.

Answer 1

The row extraction (as you call it) is basically a projection $P$, which acts from the right on the coefficient matrix

$$ A_{S_i} = A \cdot P_{S_i} $$

where $P_{S_i}$ is a matrix of dimension $N \times |S_i|$, which contains unit vectors if the corresponding row is included.

$$ P_{S_i} = \begin{pmatrix} e_{s_1}, e_{s_2}, \ldots , e_{s_{|S_i|}} \end{pmatrix} $$

With this, the task is to solve the systems

$$ (A P_{S_i}) x = y $$

for which you can you can successively solve the two problems $$ A \underbrace{(P_{S_i} x)}_{=z} = y $$ and $$ P_{S_i} x = z $$

For the first equation, you can decompose the system matrix $A$ once by a QR or SVD decomposition.

By decomposing the projection as 1 = P + (1 - P) and using the fact that a projection is idempotent $PP=P$, the second equation can be decomposed into two equations (i'll drop the subscript)

\begin{align} Px &= Pz\,,\\ 0 &= (1-P)z \end{align}

That is, you solve for the unknown x within the given subspace, and you can't solve for $z$ outside the projected space. An exact solution therefore only exists if $0 = (1-P)z$ exactly holds. In the least squares case, you should ignore this equation, and use $x = Pz$.