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

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.

• Do successive subproblems differ just by the insertion/deletion of one row? If so, your best option may be QR factorization updates (you can find them described e.g. in the Golub-Van Loan book). May 3 at 18:10
• Thank you for the suggestion. Unfortunately not, they may differ by deletion of several rows. May 3 at 19:11
• If you know all the subproblems at the start, you could use a batched routine to parallelize this well on a GPU. May 3 at 23:35
• Do the matrices of the sub problems overlap? May 4 at 6:53
• Yes, some will, but only a few, say ~ 10. May 5 at 15:08

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

• Thank you, it is indeed a projection. I will try it out! May 9 at 13:05
• @Raibyo: yes, but I applied it from the wrong side. It got to act from the left, in order to reduce the rows. May 9 at 13:40