2

Often in modeling, data processing, and optimization; the intensive portions of the code can often reduce to solving a lot of linear equations.

e.g.

solve for x1 in A x1 = b1

solve for x2 in A x2 = b2

... and so on

Sometimes the matrix 'A' doesn't change between computations, and my question is how (and when) can I take advantage of this with some pre-computation? This assumes that you're already using linear algebra libraries (e.g. BLAS, LAPACK) optimized for your system, and can't find a way to solve it in a single batch operation for which (LAPACK, Matlab, etc) already have specialized functions.

For example, one strategy is to compute and store an appropriate matrix decomposition (LUP or QR) once that would otherwise be (internally) called by the linear algebra library each time. However, I can find little guidance on which ones to use that work well with the intermediate solvers (in particular when working with LAPACK), or the relative merits of the decompositions in terms of speed or running into bad edge cases.

Note: A "bad" strategy is to compute the inverse matrix inv(A)due to speed and accuracy issues. (MatLab's documentation for inv() discusses this in greater detail )

J. Paulsen
  • 23
  • 4
  • [Memoization](https://en.wikipedia.org/wiki/Memoization) comes to mind. – Robert Harvey Jun 01 '16 at 23:44
  • @Robert Harvey. Perhaps instead of 'Pre-computation' I should say 'Memoization strategies...' in the title. There are a lot of approaches, for instance I mentioned determining either QR and LUP matrix decompositions ahead of time. The fact that there are still a fair number of options in these libraries after decade(s?) hints that the answer is 'it depends' but I'd like to have a better idea on what. – J. Paulsen Jun 02 '16 at 00:56

1 Answers1

2

Matrix decomposition is definitely the way to go here. Which decomposition you use will be determined by the structure of the systems you are are trying to solve: solving using Cholesky decomposition requires a square, symmetric, positive definite, matrix. You can solve a general square matrix using LU or QR. LU is typically used because it usually takes fewer steps than QR. QR is used for under and over determined systems like least squares, and finding eigenvalues.

If you are using LAPACK, you can call dgetrf to LU factorize a general, square matrix, and dgetrs to find the solution, using the input vector, and the factor matrices you created with dgetrf.

Charles E. Grant
  • 16,612
  • 1
  • 46
  • 73
  • Ah, this is very helpful. So would I be right in interpreting this as: "In the 'general case' where you can't rely on the matrix structure use LU and add checks in the off chance the matrix is poorly conditioned and then switch to QR. For cases when you expect to see poorly conditioned, or over and underdetermined matrices frequently or need more info (e.g. eigenvalues), use QR straight off. The rest are for when you can ensure special properties of the matrix." – J. Paulsen Jun 02 '16 at 05:29
  • Yes, but in addition I'd note that the numerical linear algebra textbook I have on hand, [Numerical Linear Algebra by Tbefethen and Bau](http://bookstore.siam.org/ot50/), assures us that while theoretically less stable, LUP is almost always fine in practice. Look at [this question](http://scicomp.stackexchange.com/questions/3262/how-to-choose-a-method-for-solving-linear-equations/3302#3302) on the Scientific Computation Stack Exchange for links to more resources. – Charles E. Grant Jun 05 '16 at 01:53