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

[anonymous]

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

                 Summary: sqrtm: returns nan for matrix of ones with rows and
columns >=4
                 Project: GNU Octave
            Submitted by: None
            Submitted on: Fri 18 Jun 2021 10:32:04 AM UTC
                Category: Octave Function
                Severity: 3 - Normal
                Priority: 5 - Normal
              Item Group: Incorrect Result
                  Status: None
             Assigned to: None
         Originator Name: 
        Originator Email: 
             Open/Closed: Open
                 Release: 6.2.0
         Discussion Lock: Any
        Operating System: Any

    _______________________________________________________

Details:

sqrtm returns nan for the matrix of ones with rows and columns >=4


A=ones(4);
sqrtm(A)
#Symmetric matrix so Schur factorization diagional
[u s]=schur(A)
sqrtmA=u*diag(diag(s).^.5)*u'
#answer matrix ones(4)*.5


The reason for the issue is all but one eigenvalue is zero and 0/0 results in
nan.

The algorithm according to the c++ comments is

A=ones(4);
n = rows (A);
[u T]=schur(A,'complex');
for j = 1:n
  T(j,j) = sqrt (T(j,j));
  for i = j-1:-1:1
    T(i,j) /= (T(i,i) + T(j,j));
    k = 1:i-1;
    T(k,j) -= T(k,i) * T(i,j);
  endfor
endfor
sqrtAA=u*T*u'
norm(sqrtAA^2-A,"fro")

The following modification fixes the issue and I don't think it will cause
other issues.

A=ones(4);
n = rows (A);
[u T]=schur(A,'complex');
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
sqrtAA=u*T*u'
norm(sqrtAA^2-A,"fro")

I think changing line 79 of the c++ file from

const element_type colji = colj[i] /= (coli[i] + colj[j]);

to

element_type colji=colj[i];
if (colj[i] != zero){
  colji = colj[i] /= (coli[i] + colj[j]);
}

will fix it but I have not tested it and am not completely sure.




    _______________________________________________________

Reply to this item at:

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

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