[Octave-bug-tracker] [bug #60797] sqrtm: returns nan for matrix of ones with rows and columns >=4

[anonymous]

Follow-up Comment #2, bug #60797 (project octave):

I think there could still be a case where the previous solution may still have
difficulty. A large Hermitian matrix with several eigenvalues found to be
exactly zero but a small numerical error above the main diagonal. The problem
would be a tiny number divided by zero equals + or - inf. We cant remove the
divide by zero because some matrices don't have a square root and nan inf
would be appropriate to return. However, the Schur decomposition is diagonal
for Hermitian matrixes so we can skip the upper triangular portion of the
matrix square root entirely. To implement this I would:

* loop over the upper triangular portion of the matrix
** If any value is greater than the tolerance break the loop and do the
existing algorithm.
** If no value exceeds the tolerance square root the diagonal and finish

The following Octave code is what I described above

A=ones(4);
n = rows (A);
[u T]=schur(A,'complex');
tol = n * eps (max (abs (diag (T))));
if (max (abs (triu (T,1))(:)) < tol)
  T=diag(sqrt(diag(T)));
else
  for j = 1:n
    T(j,j) = sqrt (T(j,j));
    for i = j-1:-1:1
      if T(i,j)~=0
        T(i,j) /= (T(i,i) + T(j,j));
      end
      k = 1:i-1;
      T(k,j) -= T(k,i) * T(i,j);
    endfor
  endfor
endif
sqrtAA=u*T*u'
norm(sqrtAA^2-A,"fro")

As an added advantage the square roots of Hermitian will be faster (no
dividing (by non-zero) and subtracting with lots of zeros) and no slower for
upper triangular matrices as it is likely the first element will be larger
than the tolerance. This suggested fix will solve what the previous fix solved
more robustly.

I'm not sure if the previous fix should be kept or removed. I think the
difference will be inf appearing in the matrix rather than nan but as the inf
will be subtracted it will quickly turn into nan. Also, is the check of zero
going to change how fast it runs?

    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?60797>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/