6
$\begingroup$

I have a large number of systems of the form:

$Ax=b_i$

To solve for a large numbers of such $b_i\;1\leq i \leq k$ but where $A$ is fixed (A is a rank $p$ general --i.e. non sparse, non PSD-- matrix).

I can solve them individually using an LU decomposition (costs $O(p^3)$) but was wondering whether there is a more efficient way to get all of the $k$ vector of solution $x^*_i$ than solving these $k$ systems independently?

The typical range of values of $p$, $k$ i'm considering are in the 10-100 (typically, i'm expecting $k\approx p$).

A pointer to a c++ implementation of whatever method is proposed would also be greatly appreciated.

Best regards,

Improve this question
$\endgroup$
2
  • $\begingroup$ Could you clarify what you mean by "jump" from one solution to the next? $\endgroup$ – Paul♦ Jun 17 '12 at 15:43
  • $\begingroup$ @Paul: i shouldn't have used that expression. It's now corrected. $\endgroup$ – user189035 Jun 17 '12 at 15:45
5
$\begingroup$

If you know all right hand sides in advance, you can combine them into a matrix and solve $AX=B$. Routines for this (in C and Fortran) can be found in the LAPACK library, which represents the state of the art.

If $k\gg p$ it is probably more efficient (than using Paul's standard way of proceeding) to compute the inverse matrix and then multiply each right hand side with it. This vectorizes (and parallelizes) much better, and the small initial overhead in computing the inverse may be ignored. Unlike in the case $k\le p$. For example, if $k=p$ then the overhead in arithmetic operations is 50%, and the gain from vectorization may not be sufficient to make up for this.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Neumaiser: yes all the right hand sides are known in advance. Can you explicit the dimensions of X and B in your answer and how you construct them from the original matrix A and the vectors $b_i$? $\endgroup$ – user189035 Jun 17 '12 at 16:53
  • 1
    $\begingroup$ With your notation, $B$ is $p \times k$, and has the $b_i$ as columns. The same for $X$ and the $x_i^*$. Then form $X=A^{-1}B$ (simple, and better for $k\gg p$) or factor $PA=LR$ and then solve $LY=PB$ by forward substitution and $RX=Y$ by back substitution. $\endgroup$ – Arnold Neumaier Jun 17 '12 at 16:55
  • $\begingroup$ I personally like this solution better for the subjective reason that it makes the program simpler to write (i.e. the loop is managed by LAPACK). $\endgroup$ – user189035 Jun 17 '12 at 16:57
  • 1
    $\begingroup$ @user189035: Well, LAPACK will do the three parts of the second version with a single call, while it needs two calls for the first version. $\endgroup$ – Arnold Neumaier Jun 17 '12 at 17:02
7
$\begingroup$

First, you need to use LU factorization for the first system of equations. This is an $\Theta(p^3)$ operation. But since you are not changing the matrix coefficients and only changing the RHS vector, you can reuse the factors L and U to "resolve" the new systems in $\Theta(p^2)$ operations.

Algorithm:

  1. Find the LU factorization of A.
  2. For i = 1 to N
  3. Solve the system $Ly_i=b_i$ (forward substitution)
  4. Use the solution $y_i$ to solve the system $Ux_i=y_i$
  5. The solution to the system with RHS vector $b_i$ is then $x_i$.
  6. End For
Improve this answer
$\endgroup$
0
$\begingroup$

If $k<p$ then it might be worthwhile to use one of the "subspace recycling" methods you get by recording the iterates of a Krylov subspace method such as GMRES. If $k\approx p$ or $k>p$ I think it should be simple to show that these subspace recycling methods won't be effective.

Improve this answer
$\endgroup$

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.