13

I have lots of large (around 5000 x 5000) matrices that I need to invert in Matlab. I actually need the inverse, so I can't use mldivide instead, which is a lot faster for solving Ax=b for just one b.

My matrices are coming from a problem that means they have some nice properties. First off, their determinant is 1 so they're definitely invertible. They aren't diagonalizable, though, or I would try to diagonlize them, invert them, and then put them back. Their entries are all real numbers (actually rational).

I'm using Matlab for getting these matrices and for this stuff I need to do with their inverses, so I would prefer a way to speed Matlab up. But if there is another language I can use that'll be faster, then please let me know. I don't know a lot of other languages (a little but of C and a little but of Java), so if it's really complicated in some other language, then I might not be able to use it. Please go ahead and suggest it, though, in case.

Daniel
  • 944
  • 1
  • 7
  • 24
  • 3
    Matlab's linear algebra routines are usually pretty optimized, so you may not get any speed improvement from re-implementing anything by hand in another language. You may want to look at Eigen (c++), though. For Matlab, you may want to try the parallel processing toolbox. – Jonas Jun 28 '11 at 16:01
  • 2
    Also, it might help if you told us _what properties_ your matrix has, or how you get them. Perhaps some linear algebra manipulations before plugging it into MATLAB might help you more than any library can. – abcd Jun 28 '11 at 16:22
  • 2
    You may like to elaborate more on your particular case. Are you 100% sure you need the inverse? Thanks – eat Jun 28 '11 at 16:22
  • I'm guessing they are sparse since most real world problems don't come up with full 5000 x 5000 matrices. Assuming they are, can you exploit the sparse structure somehow? Have you looked into this? – Chris A. Jun 28 '11 at 16:36
  • How sparse is the matrix? Does the matrix have any special shape, or relationship between entries? – nibot Jun 28 '11 at 17:51
  • What about L-U decomposition since you said that you need each matrix multiple time for solving equation? http://www.mathworks.com/help/techdoc/ref/lu.html – Luka Rahne Jun 28 '11 at 18:36
  • Do you really need the inverse? [Don’t invert that matrix](http://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/), I would say there are probably other ways to solve the problem you really want to solve. – HelloGoodbye May 25 '16 at 15:19

3 Answers3

19

I actually need the inverse, so I can't use mldivide instead,...

That's not true, because you can still use mldivide to get the inverse. Note that A-1 = A-1 * I. In MATLAB, this is equivalent to 

invA = A\speye(size(A));

On my machine, this takes about 10.5 seconds for a 5000x5000 matrix. Note that MATLAB does have an inv function to compute the inverse of a matrix. Although this will take about the same amount of time, it is less efficient in terms of numerical accuracy (more info in the link).

First off, their determinant is 1 so they're definitely invertible

Rather than det(A)=1, it is the condition number of your matrix that dictates how accurate or stable the inverse will be. Note that det(A)=∏i=1:n λi. So just setting λ1=M, λn=1/M and λi≠1,n=1 will give you det(A)=1. However, as M → ∞, cond(A) = M2 → ∞ and λn → 0, meaning your matrix is approaching singularity and there will be large numerical errors in computing the inverse.

My matrices are coming from a problem that means they have some nice properties.

Of course, there are other more efficient algorithms that can be employed if your matrix is sparse or has other favorable properties. But without any additional info on your specific problem, there is nothing more that can be said.

I would prefer a way to speed Matlab up

MATLAB uses Gauss elimination to compute the inverse of a general matrix (full rank, non-sparse, without any special properties) using mldivide and this is Θ(n3), where n is the size of the matrix. So, in your case, n=5000 and there are 1.25 x 1011 floating point operations. So on a reasonable machine with about 10 Gflops of computational power, you're going to require at least 12.5 seconds to compute the inverse and there is no way out of this, unless you exploit the "special properties" (if they're exploitable)

abcd
  • 41,765
  • 7
  • 81
  • 98
  • 1
    Does this actually provide any benefit? – nibot Jun 28 '11 at 17:50
  • @yoda: You state your method takes 10.5 seconds. I'm curious to know how long `inv` takes on the same input. – PengOne Jun 28 '11 at 18:24
  • @PengOne: `inv` roughly takes about the same time (~10 seconds), but is less accurate. See [`inv`](http://www.mathworks.com/help/techdoc/ref/inv.html), where they talk about this. – abcd Jun 28 '11 at 18:29
  • @yoda: Your answer (which I like, btw) implies that using `mldivide` is faster, when really you mean that it's more accurate and not less fast. That might be worth mentioning. – PengOne Jun 28 '11 at 18:32
  • @PengOne: I just edited my answer and made sure to add that in. – abcd Jun 28 '11 at 18:38
  • 1
    The general guidance is to use `A\b` instead of `inv(A)*b` to solve `Ax=b` **when `b` is a vector**. When `b` is a matrix— _in particular, the identity matrix, as suggested by yoda_ —I think the `mldivide` form has no benefit. – nibot Jun 28 '11 at 18:48
  • @nibot: I haven't done extensive testing, but `mldivide` still does have a _slightly_ better numerical accuracy than `inv` in this case too, although it is not as drastic as when solving a linear system of equations. We can go back and forth on how good is good enough for an engineer, but I shall accept your point for now (although I wouldn't put it as _no benefit_). Besides, my original intent was to correct OP's claim that it cannot be done using `mldivide`. The use of `inv` in MATLAB is almost universally prohibited and its presence is solely for backwards compatibility and for toy examples – abcd Jun 28 '11 at 20:24
  • 1
    @nibot, @yoda: solving routines usually perform 1) a decomposition of the matrix (usually some form of LU or QR, or Cholesky for symmetric matrices) 2) backsubstitution (solving triangular or orthogonal systems) 3) solution polish. Canned matrix inversion routines usually skip step 3 (or do full gaussian elimination). Since solving systems is way more frequent than inverting matrices, matrix inversion routines tend to be less optimized, and it can very well be the case than 1 + 2 + 3 with n right hand sides (1 has to be done once) is faster than a general purpose inversion routine nobody uses. – Alexandre C. Jun 28 '11 at 21:31
  • I doubt this is better than calculating the inverse of `A` directly. Otherwise that would mean that you had come up with a better way to calculate the matrix inverse than Mathworks. – HelloGoodbye May 25 '16 at 15:21
7

Inverting an arbitrary 5000 x 5000 matrix is not computationally easy no matter what language you are using. I would recommend looking into approximations. If your matrices are low rank, you might want to try a low-rank approximation M = USV'

Here are some more ideas from math-overflow:

https://mathoverflow.net/search?q=matrix+inversion+approximation

Community
  • 1
  • 1
BlessedKey
  • 1,615
  • 1
  • 10
  • 16
  • 2
    shmeh. Just because a matrix has det(A)==1 doesn't mean that a low rank approximation isn't good enough for engineering purposes. See e.g. A = diag([1e5,1,1e-5]). The determinent of A is still 1, but A might as well be diag([1e5 0 0]) if you are using it for some purposes. – BlessedKey Jun 28 '11 at 18:13
  • The fact that det A = 1 **does not** mean that A is numerically invertible. For large matrices, the odds that the condition number is huge are big enough that most matrices are not invertible (unless they come eg. from a discretization of an invertible operator). The SVD approach is probably the only sensible one, unless OP provides us with background on the problem. – Alexandre C. Jun 28 '11 at 18:36
  • @AlexandreC.: What is the SVD approach? What does SVD stand for? – HelloGoodbye May 25 '16 at 15:23
  • @hellogoodbye : this is the M=USV' approach from the post (Singular Value Decomposition) – Alexandre C. May 25 '16 at 15:25
1

First suppose the eigen values are all 1. Let A be the Jordan canonical form of your matrix. Then you can compute A^{-1} using only matrix multiplication and addition by

A^{-1} = I + (I-A) + (I-A)^2 + ... + (I-A)^k

where k < dim(A). Why does this work? Because generating functions are awesome. Recall the expansion

(1-x)^{-1} = 1/(1-x) = 1 + x + x^2 + ...

This means that we can invert (1-x) using an infinite sum. You want to invert a matrix A, so you want to take

A = I - X

Solving for X gives X = I-A. Therefore by substitution, we have

A^{-1} = (I - (I-A))^{-1} = 1 + (I-A) + (I-A)^2 + ...

Here I've just used the identity matrix I in place of the number 1. Now we have the problem of convergence to deal with, but this isn't actually a problem. By the assumption that A is in Jordan form and has all eigen values equal to 1, we know that A is upper triangular with all 1s on the diagonal. Therefore I-A is upper triangular with all 0s on the diagonal. Therefore all eigen values of I-A are 0, so its characteristic polynomial is x^dim(A) and its minimal polynomial is x^{k+1} for some k < dim(A). Since a matrix satisfies its minimal (and characteristic) polynomial, this means that (I-A)^{k+1} = 0. Therefore the above series is finite, with the largest nonzero term being (I-A)^k. So it converges.

Now, for the general case, put your matrix into Jordan form, so that you have a block triangular matrix, e.g.:

A 0 0
0 B 0
0 0 C

Where each block has a single value along the diagonal. If that value is a for A, then use the above trick to invert 1/a * A, and then multiply the a back through. Since the full matrix is block triangular the inverse will be

A^{-1} 0      0
0      B^{-1} 0
0      0      C^{-1}

There is nothing special about having three blocks, so this works no matter how many you have.

Note that this trick works whenever you have a matrix in Jordan form. The computation of the inverse in this case will be very fast in Matlab because it only involves matrix multiplication, and you can even use tricks to speed that up since you only need powers of a single matrix. This may not help you, though, if it's really costly to get the matrix into Jordan form.

PengOne
  • 48,188
  • 17
  • 130
  • 149
  • 1
    Why should the eigenvalues be either `±1` if the determinant is `1` and the entries are real? Try `A=[1,3;2,7]`. The characteristic polynomial is `1 - 8*x + x^2`. Its roots are `4±sqrt(15)` and not `±1`. – abcd Jun 28 '11 at 16:59
  • @Yoda: Somehow I was thinking integers... good point. I'll edit my answer. – PengOne Jun 28 '11 at 17:01
  • @Alexandre C: What part of it? It's possible I was hasty with the general case, but when all eigen values are `1`, this solution does work. – PengOne Jun 28 '11 at 17:08
  • @PengOne: you edited the answer in the meantime. You said "since det(A) = 1, the eigenvalues are +-1". Apart from that, computing the Jordan form is way slower than computing the inverse matrix. I don't know about the stability of the Jordan form, but I suspect that it heavily depends on the condition number of A, which will likely be very high since the matrix is large. But if all the eigenvalues are one, then A = I + N with N nilpotent, and the series trick works (albeit with 5000 terms in the worst case, which gives the same O(N^3) complexity as plain inversion). – Alexandre C. Jun 28 '11 at 18:19
  • @Alexandre C: Your comment came after my edit, and states that the entire answer ("this") is false. Yes, it can be costly to compute the Jordan form, but based on the line "They aren't diagonalizable, though, or I would try to diagonlize them, invert them, and then put them back." in the question, I figured this approach was worth mentioning. – PengOne Jun 28 '11 at 18:22
  • @PengOne: Fair enough. I doubt reduction will beat plain inversion here, and there may be stability concerns (for large matrices, having det A = 1 **does not** guarantee that A is numerically invertible, since this depends on the condition number). Imho, the right approach here is SVD, but this depends on what the OP wants to do. – Alexandre C. Jun 28 '11 at 18:34