GSL/BLAS: Multiply a matrix with an inverse matrix

math matrix linear-algebra gsl
Tom picture Tom · Aug 23, 2010 · Viewed 7.4k times · Source

I'm using the GNU GSL to do some matrix calculations. I'm trying to multiply a matrix B with the inverse of a matrix A.

Now I've noticed that the BLAS-part of GSL has a function to do this, but only if A is triangular. Is there a specific reason to this? Also, what would be the fastest way to do this computation? Should I invert A using LU decomposition, or is there a better way?

FWIW, A has the form P'GP, where P is a normal matrix, P' is its inverse, and G is a diagonal matrix.

Thanks a bunch :)

Answer

Alejandro Cámara picture Alejandro Cámara · Sep 5, 2010

I believe Adrien is right with the reason for BLAS not having an inverse function for square matrices. It depends on the matrix you are using to optimize the calculus of its inverse.

In general you can use the LU decomposition, which is valid for any square matrix. I. e., something like:

gsl_linalg_LU_decomp(A, p, signum);
gsl_linalg_LU_invert(A, p, invA);

where A is the square matrix you want its inverse, p is a gsl_permutation object (the permutation object where the permutation matrix is encoded), signum is the sign of the permutation and invA is the inverse of A.

Since you state that A = P' G P being P normal and G diagonal, probably A is a normal matrix. I haven't used them in a while, but there must be a factorization theorem for them, which you can find in the Chapter 14 of the GSL reference manual or even better, in a linear algebra book.

Just an example, if you had symmetric positive definite matrices and wanted to find its inverse, you could use the Cholesky factorization, which is optimized for that kind of matrices. Then you could use gsl_linalg_cholesky_decomp() and gsl_linalg_cholesky_invert() functions in the GSL to make it efficient.

I hope it helps!

Related questions

What are the most widely used C++ vector/matrix math/linear algebra libraries, and their cost and benefit tradeoffs?

It seems that many projects slowly come upon a need to do matrix math, and fall into the trap of first building some vector classes and slowly adding in functionality until they get caught building a half-assed custom linear algebra …

Read More
How to check if m n-sized vectors are linearly independent?

Disclaimer This is not strictly a programming question, but most programmers soon or later have to deal with math (especially algebra), so I think that the answer could turn out to be useful to someone else in the future. Now …

Read More
Efficient 4x4 matrix inverse (affine transform)

I was hoping someone can point out an efficient formula for 4x4 affine matrix transform. Currently my code uses cofactor expansion and it allocates a temporary array for each cofactor. It's easy to read, but it's slower than it should …

Read More