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

[Johannes Pfeifer]

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

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.

    _______________________________________________________

Reply to this item at:

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

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