bugGNU Octave - Bugs: bug #49904, Inversion of a Hermitian matrix...

 
 

bug #49904: Inversion of a Hermitian matrix produces incorrect result

Submitter:  Jonathan Lai <jlai319>
Submitted:  Wed 21 Dec 2016 04:06:58 PM UTC
   
 
Category:  Libraries Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Regression
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * 4.2.0 Operating System:  * Microsoft Windows
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Wed 21 Dec 2016 08:42:19 PM UTC, comment #3: 

I think this is the changeset that introduced the bug:

http://hg.savannah.gnu.org/hgweb/octave/rev/5ce959c55cc0

It looks like some cutting and pasting wiped out calls to std::conj...

I checked in the following changeset on stable and merged it with default:

http://hg.savannah.gnu.org/hgweb/octave/rev/db2a531717b3

John W. Eaton <jwe>
Group administrator
Wed 21 Dec 2016 07:45:46 PM UTC, comment #2: 

Here is what I see with 4.2.0:


>> Rxx = [ 1.00003 + 0.00000i  -0.68292 + 0.26406i  -0.57670 + 0.25291i   0.34986 +  0.03335i
  -0.68292 - 0.26406i   0.94539 + 0.00000i   0.69006 - 0.03967i  -0.43847 +  0.38619i
  -0.57670 - 0.25291i   0.69006 + 0.03967i   0.95260 + 0.00000i  -0.34734 +  0.40728i
   0.34986 - 0.03335i  -0.43847 - 0.38619i  -0.34734 - 0.40728i   0.98356 + 0.00000i]
Rxx =

   1.00003 + 0.00000i  -0.68292 + 0.26406i  -0.57670 + 0.25291i   0.34986 + 0.03335i
  -0.68292 - 0.26406i   0.94539 + 0.00000i   0.69006 - 0.03967i  -0.43847 + 0.38619i
  -0.57670 - 0.25291i   0.69006 + 0.03967i   0.95260 + 0.00000i  -0.34734 + 0.40728i
   0.34986 - 0.03335i  -0.43847 - 0.38619i  -0.34734 - 0.40728i   0.98356 + 0.00000i

>> Rxx * inv (Rxx)
ans =

   3.8447e+03 + 4.1605e+03i  -4.2505e+02 - 5.3359e+03i   2.1636e+02 - 2.2698e+03i  -2.2737e-13 - 1.1369e-12i
   4.8212e+03 - 6.0472e+03i   4.0808e+03 + 6.5882e+03i   2.5055e+03 + 2.8447e+03i   0.0000e+00 + 1.3642e-12i
   5.6085e+03 - 5.2994e+03i   4.1757e+03 + 6.4500e+03i   2.6433e+03 + 2.2534e+03i  -6.8212e-13 + 0.0000e+00i
  -5.1015e+03 - 7.0805e+03i   1.2351e+03 - 1.1137e+04i   4.2633e-13 - 6.3810e+03i   1.0000e+00 - 2.2737e-13i

>> inv (Rxx)
ans =

   7918.16 +    0.00i   4705.22 - 5155.62i   3582.48 - 1374.29i  -2047.24 - 6383.13i
   4705.22 - 5155.62i   6155.56 +    0.00i   3022.13 + 1516.30i   2940.40 - 5126.34i
   3582.48 - 1374.29i   3022.13 + 1516.30i   1861.79 +    0.00i    181.66 - 3243.85i
  -2047.24 - 6383.13i   2940.40 - 5126.34i    181.66 - 3243.85i   5676.72 +    0.00i


and with 4.0.3:


>> Rxx = [ 1.00003 + 0.00000i  -0.68292 + 0.26406i  -0.57670 + 0.25291i   0.34986 +  0.03335i
  -0.68292 - 0.26406i   0.94539 + 0.00000i   0.69006 - 0.03967i  -0.43847 +  0.38619i
  -0.57670 - 0.25291i   0.69006 + 0.03967i   0.95260 + 0.00000i  -0.34734 +  0.40728i
   0.34986 - 0.03335i  -0.43847 - 0.38619i  -0.34734 - 0.40728i   0.98356 + 0.00000i]
Rxx =

   1.00003 + 0.00000i  -0.68292 + 0.26406i  -0.57670 + 0.25291i   0.34986 + 0.03335i
  -0.68292 - 0.26406i   0.94539 + 0.00000i   0.69006 - 0.03967i  -0.43847 + 0.38619i
  -0.57670 - 0.25291i   0.69006 + 0.03967i   0.95260 + 0.00000i  -0.34734 + 0.40728i
   0.34986 - 0.03335i  -0.43847 - 0.38619i  -0.34734 - 0.40728i   0.98356 + 0.00000i

>> Rxx * inv (Rxx)
ans =

   1.00000 + 0.00000i  -0.00000 - 0.00000i   0.00000 - 0.00000i  -0.00000 - 0.00000i
  -0.00000 - 0.00000i   1.00000 + 0.00000i  -0.00000 + 0.00000i  -0.00000 + 0.00000i
  -0.00000 - 0.00000i  -0.00000 + 0.00000i   1.00000 - 0.00000i  -0.00000 + 0.00000i
   0.00000 + 0.00000i   0.00000 + 0.00000i   0.00000 + 0.00000i   1.00000 - 0.00000i

>> inv (Rxx)
ans =

   7918.16 +    0.00i   4705.22 - 5155.62i   3582.48 - 1374.29i  -2047.24 - 6383.13i
   4705.22 + 5155.62i   6155.56 +    0.00i   3022.13 + 1516.30i   2940.40 - 5126.34i
   3582.48 + 1374.29i   3022.13 - 1516.30i   1861.79 +    0.00i    181.66 - 3243.85i
  -2047.24 + 6383.13i   2940.40 + 5126.34i    181.66 + 3243.85i   5676.72 +    0.00i


Something seems messed up with the signs here.  I'll take a look and see if I can spot what changed.

John W. Eaton <jwe>
Group administrator
Wed 21 Dec 2016 06:30:20 PM UTC, comment #1: 

The matrix Rxx is Hermitian, but is also ill-conditioned. The condition number is on the order of 10^4, so large numerical errors in calculating its inverse are to be expected.


>> ishermitian (Rxx)
ans = 1
>> cond (Rxx)
ans =    5.7645e+04
>> [a, rcond] = inv (Rxx);
>> rcond
rcond =    1.3312e-05
>> ishermitian (a)
ans = 0
>> issymmetric (a)
ans = 1
>> abs (a - a')
ans =

   0.0000e+00   1.0311e+04   2.7486e+03   1.2766e+04
   1.0311e+04   0.0000e+00   3.0326e+03   1.0253e+04
   2.7486e+03   3.0326e+03   0.0000e+00   6.4877e+03
   1.2766e+04   1.0253e+04   6.4877e+03   0.0000e+00


Do you expect this to behave differently? Does Matlab behave differently?

Mike Miller <mtmiller>
Group Member
Wed 21 Dec 2016 04:06:58 PM UTC, original submission:  


>> Rxx

Rxx =

   1.00003 + 0.00000i  -0.68292 + 0.26406i  -0.57670 + 0.25291i   0.34986 + 0.03335i
  -0.68292 - 0.26406i   0.94539 + 0.00000i   0.69006 - 0.03967i  -0.43847 + 0.38619i
  -0.57670 - 0.25291i   0.69006 + 0.03967i   0.95260 + 0.00000i  -0.34734 + 0.40728i
   0.34986 - 0.03335i  -0.43847 - 0.38619i  -0.34734 - 0.40728i   0.98356 + 0.00000i

>> inv(Rxx)

ans =

   7412.81 +    0.00i   4404.95 - 4826.59i   3353.84 - 1286.54i  -1916.52 - 5975.72i
   4404.95 - 4826.59i   5762.91 +    0.00i   2829.13 + 1419.57i   2752.81 - 4799.16i
   3353.84 - 1286.54i   2829.13 + 1419.57i   1743.11 +    0.00i    170.06 - 3036.83i
  -1916.52 - 5975.72i   2752.81 - 4799.16i    170.06 - 3036.83i   5314.48 +    0.00i

Comment: Rxx is a Hermitian, but inv(Rxx) is not

Jonathan Lai <jlai319>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by jlai319 (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 5 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2019-02-26 mtmiller Carbon-CopyRemoved 80942 -
    2016-12-21 jwe StatusNeed Info Fixed
        Open/ClosedOpen Closed
    2016-12-21 jwe Item GroupIncorrect Result Regression
    2016-12-21 mtmiller StatusNone Need Info

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code