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

Submitter:  None
Submitted:  Fri 18 Jun 2021 10:32:04 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Ready For Test Assigned to:  None
Originator Name:  Originator Email:  -email is unavailable-
Open/Closed:  * Open Release:  * dev
Operating System:  * Any Fixed Release:  None
Planned Release:  10.1.0 (current default)
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Wed 17 Jan 2024 08:18:47 AM UTC, comment #9: 

Thanks for reminding me of this old report and the additional test result.

We are early in the release cycle of Octave 10, so I opted for pushing the patch from comment #3 to the default branch:
https://hg.savannah.gnu.org/hgweb/octave/rev/d8b5addca4ea

If issues occur, we could probably fix them during the next year or revert this change before Octave 10 will be released.

Marking as ready for test.

If it is an issue that different (but valid) solutions are returned for singular matrices on different platforms, it might be better to open a new report specifically for that.

Markus Mützel <mmuetzel>
Group administrator
Wed 17 Jan 2024 02:47:01 AM UTC, comment #8: 

Just to add some data point (ARM aarch64):

octave:1> sqrtm(ones(4))
warning: sqrtm: matrix is singular, may not have a square root
ans =

   0.5000 +      0i   0.5000 + 0.0000i   0.5000 + 0.9692i   0.5000 - 0.9692i
   0.5000 +      0i   0.5000 - 0.0000i   0.5000 - 0.3231i   0.5000 + 0.3231i
   0.5000 +      0i   0.5000 - 0.0000i   0.5000 - 0.3231i   0.5000 + 0.3231i
   0.5000 +      0i   0.5000 - 0.0000i   0.5000 - 0.3231i   0.5000 + 0.3231i

octave:2> ans^2
ans =

   1.0000 - 0.0000i   1.0000 - 0.0000i   1.0000 + 0.0000i   1.0000 - 0.0000i
   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i
   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i
   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i

octave:3> test sqrtm.cc-tst
warning: sqrtm: matrix is singular, may not have a square root
warning: called from
    __test__ at line 2 column 1
    test at line 682 column 11

PASSES 4 out of 4 tests
octave:4> version -hgid
ans = 038fd58e257e


Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Wed 17 Jan 2024 01:31:56 AM UTC, comment #7: 

@Markus: You did the heavy work on this bug report.  Is it ready to close?  It seems that the first solution in comment #1 was pushed, but that possible the patch you made in comment #3 could also be pushed.

Rik <rik5>
Group administrator
Tue 07 Dec 2021 07:58:45 AM UTC, comment #6: 

I pushed a change to stable here that modifies the test as outlined in comment #5 and also marks the test with this bug number:
https://hg.savannah.gnu.org/hgweb/octave/rev/3ab3feae06cb

Markus Mützel <mmuetzel>
Group administrator
Mon 06 Dec 2021 06:17:10 PM UTC, comment #5: 

Installed the version from octave.space here.
It fails like this:

>> sqrtm(ones(4))
warning: sqrtm: matrix is singular, may not have a square root
ans =

   0.500000   0.500000   0.024457   0.975543
   0.500000   0.500000   0.658514   0.341486
   0.500000   0.500000   0.658514   0.341486
   0.500000   0.500000   0.658514   0.341486


Seems to be a pretty good solution though:

>> ans^2
ans =

   1   1   1   1
   1   1   1   1
   1   1   1   1
   1   1   1   1


The 64-bit version calculates the expected solution:

>> sqrtm(ones(4))
warning: sqrtm: matrix is singular, may not have a square root
ans =

   0.5000 + 0.0000i   0.5000 - 0.0000i   0.5000 - 0.0000i   0.5000 - 0.0000i
   0.5000 - 0.0000i   0.5000 + 0.0000i   0.5000 + 0.0000i   0.5000 + 0.0000i
   0.5000 - 0.0000i   0.5000 + 0.0000i   0.5000 + 0.0000i   0.5000 + 0.0000i
   0.5000 - 0.0000i   0.5000 + 0.0000i   0.5000 + 0.0000i   0.5000 + 0.0000i


Do we expect that this function returns the same value on all platforms for this test?
Or should we just change the test to something like:

 assert (sqrtm (ones (4))^2, ones (4), 4*eps)


Markus Mützel <mmuetzel>
Group administrator
Mon 06 Dec 2021 05:55:41 PM UTC, comment #4: 

The 32bit Windows runner is failing this test:
https://github.com/gnu-octave/octave-buildbot/runs/4433396804?check_suite_focus=true#step:10:87936

>>>>> processing D:\a\OCTAVE~1\OCTAVE~1\OCTAVE~1\mingw32\share\octave\7.0.1\etc\tests\libinterp\corefcn\sqrtm.cc-tst
***** assert (real (sqrtm (ones (4))), 0.5*ones (4), 4*eps)
!!!!! test failed
ASSERT errors for:  assert (real (sqrtm (ones (4))),0.5 * ones (4),4 * eps)

  Location  |  Observed  |  Expected  |  Reason
   (1,3)       0.024457       0.5        Abs err 0.47554 exceeds tol 8.8818e-16 by 0.5
   (2,3)       0.65851        0.5        Abs err 0.15851 exceeds tol 8.8818e-16 by 0.2
   (3,3)       0.65851        0.5        Abs err 0.15851 exceeds tol 8.8818e-16 by 0.2
   (4,3)       0.65851        0.5        Abs err 0.15851 exceeds tol 8.8818e-16 by 0.2
   (1,4)       0.97554        0.5        Abs err 0.47554 exceeds tol 8.8818e-16 by 0.5
   (2,4)       0.34149        0.5        Abs err 0.15851 exceeds tol 8.8818e-16 by 0.2
   (3,4)       0.34149        0.5        Abs err 0.15851 exceeds tol 8.8818e-16 by 0.2
   (4,4)       0.34149        0.5        Abs err 0.15851 exceeds tol 8.8818e-16 by 0.2


I'm not sure why only the 32bit runner seems to be affected.

Markus Mützel <mmuetzel>
Group administrator
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>
Group administrator
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>
Group administrator
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 dasergatskov (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by mmuetzel (Posted a comment)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 6 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2024-01-17 mmuetzel Planned ReleaseNone 10.1.0 (current default)
    2024-01-17 mmuetzel StatusPatch Submitted Ready For Test
    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.13-4448.
    Corresponding source code