[Octave-bug-tracker] [bug #41454] Inaccurate/incorrect result for inverse of symmetric matrix

[Mattias Linde]

URL:
  <http://savannah.gnu.org/bugs/?41454>

                 Summary: Inaccurate/incorrect result for inverse of symmetric
matrix
                 Project: GNU Octave
            Submitted by: linde
            Submitted on: Tue 04 Feb 2014 11:32:42 AM GMT
                Category: Octave Function
                Severity: 3 - Normal
                Priority: 5 - Normal
              Item Group: Incorrect Result
                  Status: None
             Assigned to: None
         Originator Name: 
        Originator Email: 
             Open/Closed: Open
         Discussion Lock: Any
                 Release: dev
        Operating System: GNU/Linux

    _______________________________________________________

Details:


Double precision calculations, machine epsilon at 2.2204e-16.
Input is a symmetric matrix with condition number 4.9337e+07
so at best I cant expect results with more precision than e-09.

AB = BA = Identity, where B = inverse(A)

Since A is symmetric, B should also be symmetric.

Problems I see is that:
* B is not close to symmetric.
* norm(AB-I) is 5e-1 where as norm(BA-I) is 5e-8

Problem seen in all the versions I've tested: 3.2.4, 
3.6.0 and build from dev sources as of today.

Script:

load('hs.mat')

matrixSize = size(Hs)
matrixRank = rank(Hs)
matrixCond = cond(Hs)

machineEps = eps

I = eye( size(Hs) );

myInv  = Hs \ I;
octInv = inv(Hs);

% AB = BA = I
withBackslashEye = [ norm(full( Hs*myInv-I))  norm(full(myInv*Hs-I)) ]
withInv  = [ norm(full( Hs*octInv-I))  norm(full(octInv*Hs-I)) ]

symmetryHs     = norm(full(Hs-Hs') )
symmetryMyInv  = norm(full(myInv-myInv'))
symmetryOctInv = norm(full(octInv-octInv'))


Results from octave 

matrixSize =

   239   239

matrixRank =  239
matrixCond =    4.9337e+07
machineEps =    2.2204e-16
withBackslashEye =

   3.6139e-09   5.4825e-08

withInv =

   5.0890e-01   5.1131e-08

symmetryHs = 0
symmetryMyInv =    2.9754e-08
symmetryOctInv =  0.17683


Results from matlab 8.0:

withBackslashEye =
   1.0e-07 *
    0.0332    0.4951
withInv =
   1.0e-07 *
    0.0332    0.4951
symmetryHs =
     0
symmetryMyInv =
   2.1488e-08
symmetryOctInv =
   2.1488e-08



I also tested changing from using sparse matrices, Hs = full(Hs) and that made
things far worse.




    _______________________________________________________

File Attachments:


-------------------------------------------------------
Date: Tue 04 Feb 2014 11:32:42 AM GMT  Name: invprecision.m  Size: 434B   By:
linde

<http://savannah.gnu.org/bugs/download.php?file_id=30454>
-------------------------------------------------------
Date: Tue 04 Feb 2014 11:32:42 AM GMT  Name: hs.mat  Size: 7kB   By: linde

<http://savannah.gnu.org/bugs/download.php?file_id=30455>

    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?41454>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/