3
$\begingroup$

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.

$\endgroup$
5
  • $\begingroup$ 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). $\endgroup$ May 3 at 18:10
  • $\begingroup$ Thank you for the suggestion. Unfortunately not, they may differ by deletion of several rows. $\endgroup$
    – Raibyo
    May 3 at 19:11
  • $\begingroup$ If you know all the subproblems at the start, you could use a batched routine to parallelize this well on a GPU. $\endgroup$ May 3 at 23:35
  • $\begingroup$ Do the matrices of the sub problems overlap? $\endgroup$
    – Richard
    May 4 at 6:53
  • $\begingroup$ Yes, some will, but only a few, say ~ 10. $\endgroup$
    – Raibyo
    May 5 at 15:08

1 Answer 1

2
$\begingroup$

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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