bugGNU Octave - Bugs: bug #50846, Unexpected results with ordered qz...

 
 

bug #50846: Unexpected results with ordered qz decomposition

Submitter:  Johannes Pfeifer <jpfeifer>
Submitted:  Fri 21 Apr 2017 12:50:38 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  Johannes Pfeifer Open/Closed:  * Closed
Release:  * 4.2.1 Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Mon 24 Apr 2017 04:27:47 AM UTC, comment #7: 

I overhauled the qz function on the development branch in this changeset (http://hg.savannah.gnu.org/hgweb/octave/rev/32d532b8e7d0).  It's still not perfect, but it is a lot better.  Closing report.

Rik <rik5>
Group administrator
Sun 23 Apr 2017 02:02:16 PM UTC, comment #6: 

I re-did the documentation for qz in this cset (http://hg.savannah.gnu.org/hgweb/octave/rev/1bc0e610e293).  In addition to fixing facual errors, I used @tex blocks more to make equations appear better in formats like pdf.

Rik <rik5>
Group administrator
Sun 23 Apr 2017 09:48:54 AM UTC, comment #5: 

I got the part with the machine precision. But it did not explain why 1 was suddenly bigger than 1. But if the documentation is incorrect and the returned block has lambda strictly smaller than 1, then this explains the behavior. Thanks.

Johannes Pfeifer <jpfeifer>
Sat 22 Apr 2017 02:51:12 PM UTC, comment #4: 

I was trying to point out that the pure math value 1 has a representation on a computer which includes all adjacent values +/-eps.  And it's not clear what eps is in this case because it is a calculated value:


double eps = std::numeric_limits<double>::epsilon () * inf_norm * nn


If this calculation has eps > 1.1e-16 you will get the results shown.

It's not unexpected that eig and qz should disagree because they are not necesarilly performing the same calculation.  If you check the bottom of the documentation for qz it says


 Note: 'qz' performs permutation balancing, but not scaling (*note
     XREFbalance::).


You get better numerical results if you both scale and permute, and qz is only doing one of them.  However, you can manually do the scale operation and see the results improve.


load octave_qz.mat
format long
dd = balance (d, 'noperm');
ee = balance (e, 'noperm');
[ss,tt,w,eigval] = qz(ee,dd,'S');
abs(eigval)-1
ans =

  -1.000000000000000
  -0.808043615874273
  -0.342232732515283
  -0.100000000000000
   0.262626262626262
   0.000000000000000
   0.122499899205143
   0.010101010101010
   0.600000000000000


Now, the 6th eigenvalue is exactly 1.

I have to update the documentation for qz because the 'S' flag should indicate that the leading block contains eigenvalues |lambda| < 1, not <= 1.


Rik <rik5>
Group administrator
Sat 22 Apr 2017 08:14:42 AM UTC, comment #3: 

Dear Rik,
thanks for looking into this. I am not convinced it is entirely due to numerical precision. The eigenvalue is actually 1 and therefore should be indistinguishable from 1. Thus, it should be contained in the stable block that stores everything <=1. But qz(,'S') counts it as being strictly bigger than 1.

For this reason, one cannot rely on the output of qz() with the 'S' option to contain the stable sub-block in the upper corner for the purpose of continuing computations in Octave.

The fifth eigenvalue is actually a unit eigenvalue that is not properly recognized by qz(,'S') in Octave.
If you run

[ss,tt,w,eigvalS] = qz(e,d,'S');
res1=diag(ss)./diag(tt)
[ss,tt,w,eigvalS] = qz(e,d,'N');
res2=diag(ss)./diag(tt)
(res2(2)-1)*1e10
ans = 0
(res1(6)-1)*1e10
ans =  -1.11022302462516e-006

The eigenvalue is exactly 1 for the non-sorted case, while it becomes slightly smaller than 1 for the sorted case. But it is counted as being strictly larger than 1. Numerical precision alone cannot explain why something smaller equal to 1 is suddenly bigger than 1.

The result from qz() also contradicts the one returned by eig() in Octave


lambda = sort(abs(eig(e,d)));
abs(lambda)>1

returns

  0
  0
  0
  0
  0
  1
  1
  1
  1
  1

Here, the fifth eigenvalue is properly computed to be smaller equal to 1.

Johannes Pfeifer <jpfeifer>
Fri 21 Apr 2017 09:25:41 PM UTC, comment #2: 

This looks like it is probably roundoff error, rather than a bug.  The sixth eigenvalue is only missing being equal to 1 by eps/2.  Once you are below machine precision, odd things can happen.

Try this


load octave_qz.mat
format bank
[ss,tt,w,eigvalS] = qz(e,d,'S');
[ss,tt,w,eigvalB] = qz(e,d,'B');
abs (eigvalS)
ans =

                   0.00
                   0.43
                   0.80
                   0.90
                   1.26
                   1.00
                   1.12
   11088015243066408.00
                   2.36
                   1.01

size (eigvalS)
ans =

   10.00    1.00

abs (eigvalB)
ans =

   1.26
   1.00
   1.12
   2.36
   1.01
   0.00
   0.43
   0.80
   0.90

size (eigvalB)
ans =

   9.00   1.00


The bank format is just a quick hack to force Octave not to use scientific notation.  The sixth eigenvalue really looks like it is >=1 and is being included in the larger block.  Also note that if you reverse the sorting by using 'B' then there are only 9 eigenvalues returned.  The one that is essentially infinite was discarded.

The code for this is in libinterp/corefcn/qz.cc if you want to take a look.  Around line 698 it eventually calls a Fortran function


      F77_XFCN (dsubsp, DSUBSP,
                (nn, nn, aa.fortran_vec (), bb.fortran_vec (),
                 ZZ.fortran_vec (), sort_test, eps, ndim, fail,
                 ind.fortran_vec ()));


Importantly, the precision to solve for (the variable eps) is not necessarily the same eps as the Octave interpreter uses.  It was calculated earlier as


      double inf_norm;

      F77_XFCN (xdlange, XDLANGE,
                (F77_CONST_CHAR_ARG2 ("I", 1),
                 nn, nn, aa.data (), nn, work.fortran_vec (), inf_norm
                 F77_CHAR_ARG_LEN (1)));

      double eps = std::numeric_limits<double>::epsilon () * inf_norm * nn;




Rik <rik5>
Group administrator
Fri 21 Apr 2017 12:56:03 PM UTC, comment #1: 

Sorry, it seems the markup got mangled. The report should continue:
The code returns


ans =
  -1.00000000000000e+000
  -5.71920035902214e-001
  -2.00000000000000e-001
  -1.00000000000000e-001
  2.62626262626263e-001
  -1.11022302462516e-016
  1.22334455667789e-001
  1.10880152430664e+016
  1.35960823868476e+000
  1.01010101010102e-002

Thus, the fifth eigenvalue has |lambda|>1 while the sixth one again has |lambda|<1 and should therefore still be in the first block.

Johannes Pfeifer <jpfeifer>
Fri 21 Apr 2017 12:50:38 PM UTC, original submission:  

According to  https://octave.sourceforge.io/octave/function/qz.html
the 'S' flag of the qz-decomposition should return a leading block has all eigenvalues |lambda|≤ 1

But running

load octave_qz.mat
[ss,tt,w,eigval] = qz(e,d,'S');
abs(eigval)-1

returns


ans =

  -1.00000000000000e+000
  -5.71920035902214e-001
  -2.00000000000000e-001
  -1.00000000000000e-001
  2.62626262626263e-001
  -1.11022302462516e-016
  1.22334455667789e-001
  1.10880152430664e+016
  1.35960823868476e+000
  1.01010101010102e-002
+verbatim+
Thus, the fifth eigenvalue has |lambda|>1 while the sixth one again has |lambda|<1 and should therefore still be in the first block.

Johannes Pfeifer <jpfeifer>

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #40465:  octave_qz.mat added by jpfeifer (2KiB - application/octet-stream - data file with matrices for QZ)

 

Depends on the following items: None found

Items that depend on this one: None found

 

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

    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
    2017-04-24 rik5 StatusNone Fixed
        Open/ClosedOpen Closed
    2017-04-21 rik5 Operating SystemMicrosoft Windows Any
    2017-04-21 jpfeifer Attached File- Added octave_qz.mat, #40465
        Carbon-Copy- Added -email is unavailable-

    Back to the top

    Powered by Savane 3.13-02a9.
    Corresponding source code