[Octave-bug-tracker] [bug #50846] Unexpected results with ordered qz decomposition

[Rik]

Follow-up Comment #4, bug #50846 (project octave):

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.



    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?50846>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/