Search code examples
matlabmatrixoctavematrix-inverse

Why inv(matrix)*matrix is not exact identity matrix in Octave?


Why inv(A)*A is not exact identity matrix? All the diagonal elements are correct but rest are not. I learnt that this is residual error, then how to deal with it?

CODE:

A = [1,2,0;0,5,6;7,0,9]
A_inv = inv(A)
A_invA = inv(A)*A

OUTPUT:

code


Solution

  • An exploration of the documentation of inv leads you down the following path, which answers your question nicely (emphases mine):

    octave:1> help inv
    'inv' is a built-in function from the file libinterp/corefcn/inv.cc
    
     -- X = inv (A)
     -- [X, RCOND] = inv (A)
     -- [...] = inverse (...)
         Compute the inverse of the square matrix A.
    
         Return an estimate of the reciprocal condition number if requested,
         otherwise warn of an ill-conditioned matrix if the reciprocal
         condition number is small.
    
         In general it is best to avoid calculating the inverse of a matrix
         directly.  For example, it is both faster and more accurate to
         solve systems of equations (A*x = b) with 'Y = A \ b', rather than
         'Y = inv (A) * b'.
    
    

    In your particular case, you will see that:

    A = [1,2,0;0,5,6;7,0,9];
    [X, RCOND] = inv(A);
    RCOND
    % RCOND = 0.070492
    

    So, what does this value mean? You can find the answer in the relevant function rcond, which calculates this value directly:

    octave:2> help rcond
    'rcond' is a built-in function from the file libinterp/corefcn/rcond.cc
    
     -- C = rcond (A)
         Compute the 1-norm estimate of the reciprocal condition number as
         returned by LAPACK.
    
         If the matrix is well-conditioned then C will be near 1 and if the
         matrix is poorly conditioned it will be close to 0.
    
         [...]
    
         See also: cond, condest.
    

    Your value is 0.07, which is quite close to 0, therefore your A matrix is rather poorly conditioned.

    To learn a bit more what "poorly conditioned" means exactly, we can have a look at the cond function:

    octave:26> help cond
    'cond' is a function from the file /opt/octave-6.2.0/share/octave/6.2.0/m/linear-algebra/cond.m
    
     -- cond (A)
     -- cond (A, P)
         Compute the P-norm condition number of a matrix with respect to
         inversion.
    
         'cond (A)' is defined as 'norm (A, P) * norm (inv (A), P)'.
    
         [...]
    
         The condition number of a matrix quantifies the sensitivity of the
         matrix inversion operation when small changes are made to matrix
         elements.  Ideally the condition number will be close to 1.  When
         the number is large this indicates small changes (such as underflow
         or round-off error) will produce large changes in the resulting
         output.  In such cases the solution results from numerical
         computing are not likely to be accurate.
    

    In your case:

    cond(A,2)
    % ans = 7.080943875445246
    

    So there you have it. Your matrix is relatively poorly-conditioned, which means that its inversion is more susceptible to precision errors. You may get better results if you use the mldivide (i.e. \ operator) instead.