1

I am using Octave 3.8.1, a Matlab-like program. I'd like to generalize 1/x to the case where x may be a scalar or a matrix. Replacing 1/x with inv(x) or pinv(x) works for most x, except:

octave:1> 1/inf
ans = 0

octave:2> pinv([inf])
ans = NaN

octave:3> inv([inf])
warning: inverse: matrix singular to machine precision, rcond = 0
ans = Inf

Should I convert NaN to 0 afterwards to get this to work? Or have I missed something? Thanks!

Racing Tadpole
  • 4,270
  • 6
  • 37
  • 56
  • `pinv` calculates the **pseudo-inverse** of a matrix. It does not have the same implications as it does for scalar values. – rayryeng Feb 11 '15 at 23:37
  • 1
    MATLAB results: `1/inf -> 0`, `pinv([inv]) -> error`, `inv([inf]) -> 0 + warning`. Your best bet will be to not rely on those `inf`-shenanigans at all... – knedlsepp Feb 11 '15 at 23:39
  • interesting - thanks knedlsepp. There's a typo in your `pinv` call though - the argument should be `inf` not `inv` - does this change the result in Matlab? – Racing Tadpole Feb 12 '15 at 00:43
  • @RacingTadpole: Sorry for the typo. `pinv([inf])` also gives an error. – knedlsepp Feb 16 '15 at 14:01

1 Answers1

1

The Moore–Penrose pseudo inverse, which is the basis for Matab and octave's pinv, is implemented via completely different algorithm than the inv function. More specifically, singular value decomposition is used, which require's finite-valued matrices (they also can't be sparse). You didn't say if your matrices are square or not. The real use of pinv is for solving non-square systems (over- or underdetermined).

However, you shouldn't be using pinv or inv for your application, no matter the dimension of your matrices. Instead you should use mldivide (octave, Matlab), i.e., the backslash operator, \. This is much more efficient and numerically robust.

A1 = 3;
A2 = [1 2 1;2 4 6;1 1 3];
A1inv = A1\1
A2inv = A2\eye(size(A2))

The mldivide function handles rectangular matrices too, but you will get different answers for underdetermined systems compared to pinv because the two use different methods to choose the solution.

A3 = [1 2 1;2 4 6]; % Underdetermined
A4 = [1 2;2 4;1 1]; % Overdetermined
A3inv = A3\eye(min(size(A3))) % Compare to pinv(A3), different answer
A4inv = A4\eye(max(size(A4))) % Compare to pinv(A4), same answer

If you run the code above, you'll see that you get a slightly different result for A3inv as compared to what is returned by pinv(A3). However, both are valid solutions.

Community
  • 1
  • 1
horchler
  • 18,384
  • 4
  • 37
  • 73
  • The real benefit of the backslash operator as opposed to `inv` is for directly solving linear systems (especially if there is only one right hand side it is stable+faster than inversion). There is no numerical benefit in doing `Ainv = A\eye(size(A)); Ainv*V`, as opposed to `inv(A)*V`! The way I read this answer it suggests to use `Ainv = A\eye(size(A));` instead of `Ainv = inv(A)`, which is of no use. – knedlsepp Feb 15 '15 at 17:47
  • @knedlsepp: Yes, this answer absolutely suggests using something like `A\eye(size(A))` instead of `inv(A)` if one just wants the inverse! Of course, [rarely does one need to calculate an explicit matrix inverse](http://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/) (but it [does occur](http://blogs.mathworks.com/loren/2007/05/16/purpose-of-inv/)). `mldivide` is much more stable and has less numerical error (especially as the matrix becomes more singular). For explicit inverses, in some cases `inv(A)` is faster, but these days `A\eye(size(A))` is often faster or comparable. – horchler Feb 15 '15 at 19:17
  • Do you have any sources that would back the assumption that using `mldivide` will be more stable? I highly doubt this. – knedlsepp Feb 15 '15 at 19:21
  • @knedlsepp: This is right in the documentation for [`inv`](http://mathworks.com/help/matlab/ref/inv.html). It's fundamental to the algorithm used. It still doesn't mean that one should use `Ainv=A\eye(size(A))` to calculate an explicit inverse and use the result to solve a linear system, e.g., `x=Ainv*b`, instead of doing `x=A\b` even if this is slightly less bad than `x=inv(A)*b`. It only means that if one needs/wants the inverse, `A\eye(size(A)` is more robust. If you know your matrix isn't singular and has a low condition number then `inv` can be perfectly fine and return identical results. – horchler Feb 15 '15 at 19:40
  • I think you may misinterpret the docs there. It states that solving a linear system `A\b` is numerically more stable than doing `inv(A)*b`, which I perfectly agree upon. But generating the inverse using a different approach (as for example using `A\eye(size(A))`) will have the same numerical problems as going with `inv` in the first place. The instability stems from a possibly ill conditioned inverse and the matrix multiplication. `mldivide` is just more robust when you use `A\b` because there are iterative methods to calculate the result. – knedlsepp Feb 15 '15 at 19:51
  • @knedlsepp: No, as I said, it is fundamental to the algorithm. There is no such thing an "ill-conditioned inverse," only an ill-conditioned matrix, i.e., one that is nearly singular. If a matrix is singular, or nearly so, both approaches will have problems and issue warnings, and find inverse and solving a linear system will both be problematic. Using `A\eye(size(A))` is equivalent to solving n linear systems via Gaussian elimination. A benefit of `inv`, however, is that it uses less memory, though one can use `speye`. – horchler Feb 15 '15 at 21:03
  • Beyond this we're arguing into the wind. You might be better off asking a question at [MatlabCentral](http://mathworks.com/matlabcentral/) or possibly [SciComp.StackExchange](http://scicomp.stackexchange.com/questions). And/or see if you can come up with a concrete demonstration that `inv` is a better or equal choice for a wide set of test cases. – horchler Feb 15 '15 at 21:06
  • Yes, there is such a thing as an "ill-conditioned inverse", as the condition number of the inverse will be simply the same as the original matrix (up to floating point precision and invertibility of course). And the condition of this inverse matrix is what's in our interest, as it is the multiplication that makes `inv(A)*x` unstable. It is neither true that `inv` uses less memory, nor that it is less stable. I'm sorry if I might sound too harsh, it's just I've heard this rumour quite often and I don't like it to be spread even further. – knedlsepp Feb 15 '15 at 21:33
  • http://stackoverflow.com/questions/1419580/why-is-matlabs-inv-slow-and-inaccurate – knedlsepp Feb 15 '15 at 21:46
  • @knedlsepp: I think you're the one that needs to back up your assumptions then. Condition numbers are for the original matrix, `A`, *before* it is inverted. If `A` is singular then `A` is ill-conditioned and it won't matter what inverting `A` does if one wants to solve AX=b. Allocating `eye(A)` for large `A` can use much more memory, which is why I said that `inv` can use less, so I don't know why you say otherwise. JIT compilation may be able to do something about it, but not always and not for older versions. – horchler Feb 15 '15 at 21:59
  • I agree with the link you posted –the point of Loren's blog post on `inv` is more about avoiding calculating a explicit inverse as one should normally try to do. But it does not address the issue of `inv(A)` vs. `A\eye(size(A))`. – horchler Feb 15 '15 at 22:10
  • Didn't Loren's post clarify enough? In addition to me mentioning that inversion with theoretical arbitrary precision having the same condition number, Loren tells us that in practice using any numerical inverse will result in the square of the original condition number of the overall system. If we are talking about inversion there is obviously no place for singular matrices. The base LAPACK routines incorporate many in-place algorithms, so it is no surprise that MATLAB will take advantage of those! Have a look at Edric's remarks in the linked question (he also works at Mathworks)! – knedlsepp Feb 15 '15 at 22:11
  • I never said `inv` was bad. Reread what I've written. I've tried to be very careful with my words and semantics so please don't trivialize this conversation with such trite statements. I've even argued elsewhere on this site that `inv` is acceptable when appropriate and the faster solution. I've also suggested reasons why MathWorks would leave the existing `inv(A)` and not use `A\eye(size(A))`: memory for one, and possibly speed in other cases. Another would be backwards compatibility as the two do not produce the same numeric results. – horchler Feb 15 '15 at 22:28
  • *this answer absolutely suggests using something like `A\eye(size(A))` instead of `inv(A)` if one just wants the inverse!*. There is absolutely no reason to do so. [Try it yourself](http://ideone.com/rAbn9S)! – knedlsepp Feb 15 '15 at 22:55
  • Sorry for my verbal lapse in the surge of emotion by the way! My sole intention was to make some clarifications on this issue! – knedlsepp Feb 15 '15 at 23:14
  • I updated the [test code](http://ideone.com/rAbn9S) to include the matrices from the [`gallery`](http://de.mathworks.com/help/matlab/ref/gallery.html) that were mentioned to be 'ill-conditioned' to make for a more conclusive comparison. Running R2014a, `inv` gives better results in every case except for the Toeplitz matrix `'prolate'`. – knedlsepp Feb 16 '15 at 00:51
  • Here is a test for all the matrices of `gallery`, each for varying sizes: http://ideone.com/ZQu2Gj 261 times `inv` gives better results than `mldivide`-`eye`. 111 times `mldivide`-`eye` gives better results than `inv`. To me this a compelling case *for* `inv` not *against*! – knedlsepp Feb 16 '15 at 14:10