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

 
 

bug #50846: Unexpected results with ordered qz decomposition

Submitted by:  Johannes Pfeifer <jpfeifer>
Submitted on:  Fri 21 Apr 2017 12:50:38 PM UTC  
 
Category: Octave FunctionSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Incorrect Result
Status: FixedAssigned to: None
Originator Name: Johannes PfeiferOpen/Closed: Closed
Release: 4.2.1Operating System: Any

Add a New Comment(Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

(Jump to the original submission 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>
Project 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>
Project 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:

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

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.

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>
Project 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
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

returns
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

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

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

Rik <rik5>
Project 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

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
returns

Johannes Pfeifer <jpfeifer>

 

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

Attach File(s):
   
   
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
  • -unavailable- added by rik5 (Posted a comment)
  • -unavailable- added by jpfeifer (Submitted the item)
  • -unavailable- added by jpfeifer
  •  

    Do you think this task is very important?
    If so, you can click here to add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follow 5 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Mon 24 Apr 2017 04:27:47 AM UTCrik5StatusNone=>Fixed
      Open/ClosedOpen=>Closed
    Fri 21 Apr 2017 09:25:41 PM UTCrik5Operating SystemMicrosoft Windows=>Any
    Fri 21 Apr 2017 12:50:38 PM UTCjpfeiferAttached File-=>Added octave_qz.mat, #40465
      Carbon-Copy-=>Added -unavailable-

    Back to the top


    Powered by Savane 3.1-cleanup1