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