bugGNU Octave - Bugs: bug #60797, sqrtm: returns nan for matrix of...

 
 

bug #60797: sqrtm: returns nan for matrix of ones with rows and columns >=4

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:  Patch Submitted Assigned to:  None
Originator Name:  Originator Email:  -email is unavailable-
Open/Closed:  Open Release:  dev
Operating System:  Any

Add a New Comment (Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

Sun 01 Aug 2021 03:01:38 PM UTC, comment #3: 

The attached patch adds the shortcut for Hermitian matrices along the lines described in comment #2.
For performance reasons, I opted to take the square root of the diagonal elements in-place and don't touch the off-diagonal elements of the Schur matrix. Those off-diagonal elements are checked to be below the tolerance. But they might be non-zero.

Seems to be working correctly for me.

(file #51725)

Markus Mützel <mmuetzel>
Project Member
Fri 18 Jun 2021 11:26:57 PM UTC, comment #2: 

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?

Anonymous
Fri 18 Jun 2021 02:29:26 PM UTC, comment #1: 

Thanks for the report and the description of a fix.

I pushed a change to the default branch here:
https://hg.savannah.gnu.org/hgweb/octave/rev/90d7137b7dc6

Marking as ready for test.

Markus Mützel <mmuetzel>
Project Member
Fri 18 Jun 2021 10:32:04 AM UTC, original submission:  

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.

Anonymous

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #51725:  bug60797-sqrtm-hermitian.patch added by mmuetzel (5KiB - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by mmuetzel (Posted a comment)
  •  

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

     

     

    Follow 4 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2021-08-01 mmuetzel Attached File- => Added bug60797-sqrtm-hermitian.patch, #51725
        StatusReady For Test => Patch Submitted
        Release6.2.0 => dev
    2021-06-18 mmuetzel StatusNone => Ready For Test

    Back to the top


    Powered by Savane 3.6