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

You are not logged in

Discussion

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.

Markus Mützel <mmuetzel>
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

Markus Mützel <mmuetzel>
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'

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

Attached Files

(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)

Dependencies

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?
This task has 0 encouragements so far.

Only project members can vote.