Inverting a matrix of any size

Question

I'm using the GNU Scientific Library in implementing a calculator which needs to be able to raise matrices to powers. Unfortunately, there doesn't seem to be such a function available in the GSL for even multiplying matrices (the gsl_matrix_mul_elements() function only multiplies using the addition process), and by extension, no raising to powers.

I want to be able to raise to negative powers, which requires the ability to take the inverse. From my searching around, I have not been able to find any sound code for calculating the inverses of arbitrary matrices (only ones with defined dimensions), and the guides I found for doing it by hand use clever "on-paper tricks" that don't really work in code.

Is there a common algorithm that can be used to compute the inverse of a matrix of any size (failing of course when the inverse cannot be calculated)?

Only square matrixes (and then only some of them) can be inverted. — pmg, May 31 '15 at 15:55
That's fine. However, I want one function that can do this for any size square matrix. — 2mac, May 31 '15 at 16:16
Once you calculate the inverse you can just write a simple recursive function to quickly calculate the negative power. — Intrepid Traveller, May 31 '15 at 18:48
Assuming you're not interested in getting into the weeds of calculating an LU decomposition and inverting that using `GSL`, the C Apophenia open statistical library (http://apophenia.info) might be what you're looking for. There's a bit of a learning curve to wrap your head around the data structures and what not, but once you understand how it works, a *ton* of stuff in GSL (such as matrix inversion) is well abstracted. Once you've read through the gentle introduction, see http://apophenia.info/group__linear__algebra.html for more info on matrix inversion specifically. — Byte Lab, May 31 '15 at 23:12

score 0 · Answer 1 · answered Jun 01 '15 at 13:00

0

As mentioned in the comments, power of matrices can be computed for square matrices for integer exponents. The n power of A is A^n = A*A*...A where A appears n times. If B is the inverse of A, then the -n power of A is A^(-n) = (A^-1)^n = B^n = B*B*...B.

So in order to compute the n power of A I can suggest the following algorithm using GSL:

gsl_matrix_set_identity();         // initialize An as I
for(i=0;i<n;i++) gsl_blas_dgemm(); // compute recursive product of A

For computing B matrix you can use the following routine

gsl_linalg_LU_decomp();       // compute A decomposition
gsl_linalg_complex_LU_invert  // comput inverse from decomposition

answered Jun 01 '15 at 13:00

ztik

3,482
16
34

This explanation does not explain the method for obtaining the `gsl_permutation_t` required by `gsl_linalg_LU_decomp()` as a parameter. The data in it need to already be initialized. – 2mac Jun 10 '15 at 04:02
@2mac `gsl_permutation_t` is not always needed. It depends on the way the data are stored. The question does not give any information on it. – ztik Jun 10 '15 at 06:57
But, when looking at the source of `gsl_linalg_LU_decomp()`, the second thing it does is check that the permutation's size is equal to the matrix's `size1` value. If I do not initialize it properly, this will return false (well, really undefined), and the call will fail. – 2mac Jun 10 '15 at 11:55
@2mac I do not get the source of the problem, but note that only square matrices can be invertible. – ztik Jun 10 '15 at 12:13
`gsl_linalg_LU_decomp()` requires a `gsl_permutation` as its second argument. This `gsl_permutation` is used within the function, and it needs to have certain values based on the matrix it's given as its first argument. If I do not properly initialize the permutation, the call will fail, because there's no telling what data I give it, and just initializing the size is not an option, since the permutation value is used later on in the decomp. This has nothing to do with whether the matrix is square. It has to do with memory and how data need to be initialized, which is elementary in C. – 2mac Jun 10 '15 at 12:23

score 0 · Answer 2 · answered Jun 02 '15 at 13:50

I recommend reading up about the SVD, which the gsl implements. If your matrix is invertible, then computing it via the SVD is a not bad, though somewhat slow, way to go. If your matrix is not invertible, the SVD allows you to compute the next best thing, the generalised inverse.

In matrix calculations the errors inherent in floating point arithmetic can accumulate alarmingly. One example is the Hilbert matrix, an innocent looking thing with a remarkably large condition number, even for quite moderate dimension. A good test of an inversion routine is to see how big a Hilbert matrix it can invert, and how close the computed inverse times the matrix is to the identity.

Inverting a matrix of any size

2 Answers2