bugGNU Octave - Bugs: bug #55256, Wrong result raise a...

 
 

bug #55256: Wrong result raise a non-diagonalizable matrix by a fractional exponent

Submitter:  None
Submitted:  Thu 20 Dec 2018 07:54:01 PM UTC
   
 
Category:  Octave Function Severity:  2 - Minor
Priority:  3 - Low Item Group:  Missed Error or Warning
Status:  None Assigned to:  None
Originator Name:  Volker Becker Originator Email:  -email is unavailable-
Open/Closed:  * Open Release:  * dev
Operating System:  * Any Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Fri 21 Dec 2018 07:44:57 AM UTC, comment #2: 

Here is what matlab does:


>> version
ans = '9.4.0.813654 (R2018a)'

>> A=[1 1 1; 0 1 1; 0 0 1];

>> A^0.5
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.930381e-32.
> In matlab.internal.math.mpower.viaEig (line 10)

ans =

     1     0     0
     0     1     0
     0     0     1


A.R. Burgers <arb>
Fri 21 Dec 2018 01:13:08 AM UTC, comment #1: 

Confirmed.

If you calculate the condition number with respect to eig, you can see that this is going to be a problem.


A=[1 1 1; 0 1 1; 0 0 1];
condeig (A)
warning: matrix singular to machine precision, rcond = 2.46519e-32
warning: called from
    condeig at line 94 column 8
ans =

   6.3691e+15
   2.0282e+31
   2.0282e+31


Also, for grins, in this case one could use sqrtm rather than X^0.5 and that functions succeeds.


sqrtm (A)
ans =

   1.00000   0.50000   0.37500
   0.00000   1.00000   0.50000
   0.00000   0.00000   1.00000


So, it looks like at a minimum Octave should warn when the matrix is ill-conditioned and the results are unlikely to be accurate.

Rik <rik5>
Group administrator
Thu 20 Dec 2018 07:54:01 PM UTC, original submission:  

If raise a non-diagonalizable matrix to a fractional power, the result is wrong.

For example:



>> A=[1 1 1; 0 1 1; 0 0 1]
A =

   1   1   1
   0   1   1
   0   0   1

>> B = A^0.5
B =

   1.00000   0.00000   0.00000
   0.00000   1.00000   0.00000
   0.00000   0.00000   1.00000

So B*B is obiously not equal to A. probaly the result will be calculated by diagonalistion of the matrix which fails when the matrix is  non-diagonalizable.

The correct result (due to wolfram aplha) would be

C=[1 0.5 0.375; 0 1 0.5; 0 0 1]
C =

   1.00000   0.50000   0.37500
   0.00000   1.00000   0.50000
   0.00000   0.00000   1.00000

>> C*C
ans =

   1   1   1
   0   1   1
   0   0   1



Anonymous

 

(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 arb (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by None (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 4 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2019-02-27 mtmiller Severity3 - Normal 2 - Minor
        Priority5 - Normal 3 - Low
        Item GroupNone Missed Error or Warning
        Release4.4.1 dev

    Back to the top

    Powered by Savane 3.13-d3ae.
    Corresponding source code