bugGNU Octave - Bugs: bug #60738, logm returning incorrect result...

 
 

bug #60738: logm returning incorrect result with some real non-symmetric matrices

Submitter:  None
Submitted:  Sun 06 Jun 2021 06:58:10 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * 6.2.0
Operating System:  * Any Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Fri 18 Jun 2021 04:05:01 PM UTC, comment #14: 

I pushed the patch from comment #10 with you as author here:
https://hg.savannah.gnu.org/hgweb/octave/rev/c44e54bb018b

Marking as ready for test.

Markus Mützel <mmuetzel>
Group administrator
Fri 18 Jun 2021 01:40:37 PM UTC, comment #13: 

The patch in comment #11 doesn't apply for me. But the gist of it looks similar to the patch in comment #10.

OT with respect to handling Mercurial repositories:
I've grown to like TortoiseHg for its easy to use GUI. With it, it is quite straight-forward to see at which revision the repository is and which changesets are applied on which parents.

I don't know which OS you are using. I'm using Ubuntu. Unfortunately, TortoiseHg is no longer packaged in that distribution. IIRC, I built it from sources with the following steps:

mkdir -p ~/Repositories/tortoisehg
cd ~/Repositories/tortoisehg/
hg clone https://foss.heptapod.net/mercurial/tortoisehg/thg -r stable
pip install pyqt5 qscintilla

cd ~/Repositories/tortoisehg/thg/
pip install mercurial --upgrade
make local
sudo python3 setup.py install


The last four lines can be used to install a newer version once in a while (after updating the repo [in thg]).

I created an icon in the dock (or whatever that bar on the left is called), by creating the file `~/.local/share/applications/tortoisehg.desktop` with the following content and adding the corresponding app to favorites:

[Desktop Entry]
Name=TortoiseHG
Exec=/usr/local/bin/thg
Comment=Launch TortoiseHG
Terminal=false
Type=Application
Icon=/usr/local/share/pixmaps/tortoisehg/icons/thg_logo.ico


Most of this was taken from here: https://askubuntu.com/questions/1182598/how-to-install-tortoisehg-on-ubuntu-19-10

Markus Mützel <mmuetzel>
Group administrator
Fri 18 Jun 2021 09:22:02 AM UTC, comment #12: 

We were typing at the same time (I took over an hour so did not see the patch you uploaded) I tried to load your patch but I get a failed message. I assume I need to do something to fix my local copy. I don't know where to start. From looking at the text it looks like what I did. Compare with the logm I uploaded in the last comment.

sqrtm is used in the original algorithm to scale the matrix before integrating.

The changes which were made yesterday did not cause any problems not sure why I said revert. The current version is not worse.

Anonymous
Fri 18 Jun 2021 08:28:09 AM UTC, comment #11: 

Attached new patch for logm incorporating the speed up and the .m file in case it does not work. The patch modifies the file after the updates you pushed earlier.

logm of ones(4) is likely to be different on each implementation as all of the eigenvalues excluding one will be zero and the log of zero is -inf. I chose -log(realmax()) for negative inf to try to avoid inf in the matrix.

I tried to test the c++ of the sqrtm file but failing. I think the line should be

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

to fix the nan issue with sqrtm(ones(4))


(file #51580, file #51581)

Anonymous
Fri 18 Jun 2021 07:04:46 AM UTC, comment #10: 

The issue with the patch not applying might have been that there were local changes before that change? At least, I cannot find hg id 9503a5411cce5cd4e11a995d02a850d6630ab58f in the public repository.

In general, patches are preferred over copies of the files with changes. One advantage of patches is that they make clearer what the author intended to change and make it easier to distinguish those changes from other ones that might have occurred in the public repository in the mean time.

It's not clear to me why we would need to revert the changes we made so far. I can't see why the current version would be any worse than what we had before.

Thanks for taking the time for explaining the motivation for the new code path you propose. I think it should be save to apply those changes on the stable branch. Please, check if the attached patch that applies to the current stable branch fits your intensions.
I set you as the author of the patch. I hope that is ok with you.

I can't comment on why the algorithm uses `sqrtm`. But there are references to the algorithm in the code:

## Reference: N. J. Higham, Functions of Matrices: Theory and Computation
##            (SIAM, 2008.)
[...]
    ## Algorithm 11.9 in "Function of matrices", by N. Higham


According to TheMathworks' documentation for that function[1], they are using a different (newer) algorithm:

> [1] Al-Mohy, A. H. and Nicholas J. Higham, “Improved inverse scaling and squaring algorithms for the matrix logarithm,” SIAM J. Sci. Comput., 34(4), pp. C153–C169, 2012
>
> [2] Al-Mohy, A. H., Higham, Nicholas J. and Samuel D. Relton, “Computing the Frechet derivative of the matrix logarithm and estimating the condition number,” SIAM J. Sci. Comput.,, 35(4), pp. C394–C410, 2013


If we really wanted identical results to Matlab (I'm not sure if this is the case), we would probably need to implement the same algorithm.

If there are issues with the algorithm used by `sqrtm`, we should probably best tackle them in a separate report. Could you please open a new one for that?
I didn't look in detail at the changes you proposed for `sqrtm`. But something seems to be odd with the scoping of `colji` in the changed C++ code.

[1]: https://de.mathworks.com/help/matlab/ref/logm.html


(file #51579)

Markus Mützel <mmuetzel>
Group administrator
Fri 18 Jun 2021 01:09:01 AM UTC, comment #9: 

Typo with the sqrtm c++ fix it should be

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

I have not tested this I only tested the Octave code from the comments of the c++ file.

Anonymous
Fri 18 Jun 2021 12:57:03 AM UTC, comment #8: 

I have attached the logm.m file. I will ask you to revert the changes just applied and follow the attached file. I am confused why the patch did not work, I updated my local copy just before starting working. The way the patch marked the differences would have been very hard to follow as I effectively made two paths with an if statement, the original code was essentially unmodified and I put in a second path for Hermitian Matrices.

The following will describe how the current algorithm works:

It performs the following integral, (a definition of the log function)

ln(A)=int_I^A y^-1 dy

It uses a quadrature with limits between 0 and 1 so if we apply a linear coordinate transformation to change the lower limit from I to 0 and the upper limit from A to 1 using y=(A-I)x+I we get the following integral

ln(A)=int_0^1 (A-I)((A-I)x+I)^-1 dx
The following lines in the function do this

s -= eye (n);
if (m > 1)
  s = logm_pade_pf (s, m);
endif

The following code is equivelent but much slower

n=3;
A=randn(n);
A*=A';#quad not a fan of complex numbers so make A symmetric and pos def
logmA=zeros(n);
for i=1:n
  for j=1:n
    logmA(i,j)=quad(@(x)((A-eye(n))/((A-eye(n)).*x+eye(n)))(i,j),0,1);
  end
end
logmA
max(abs(expm(logmA)-A)(:))

I also believe the following comment in the file is incorrect.

##LOGM_PADE_PF   Evaluate Pade approximant to matrix log by partial fractions.
##   Y = LOGM_PADE_PF(A,M) evaluates the [M/M] Pade approximation to
##   LOG(EYE(SIZE(A))+A) using a partial fraction expansion.

I have shown above we are obtaining logm by its integration definition. The Padé approximant is in the form (a0+a1*x+a2*x^2+...)/(1+b1*x+b2*x^2+...) where a and b are found by solving a system of equations formed from
log(x)=(a0+a1*x+a2*x^2+...)/(1+b1*x+b2*x^2+...)
d(log(x))dx=d((a0+a1*x+a2*x^2+...)/(1+b1*x+b2*x^2+...))dx
d^2(log(x))dx^2=d^2((a0+a1*x+a2*x^2+...)/(1+b1*x+b2*x^2+...))dx^2
...
where x is the location to form the approximant about x=I in our case.

If the method is obtaining the log with integration what is the while loop doing with sqrtm? It is rescaling the matrix as square rooting a matrix makes it closer to the identity allowing the integral to be found using fewer quadrature points. As an extra bonus the first term of the Taylor series of log ln(A+I)=sum_(k=1)^inf A^k/k is A so log(A) approx equal A-I when A is almost the identity. Therefore, there is a small chance the rescaling loop may also find the matrix log.

Why is the algorithm having issues with the log of ones(4)? There is a zero divided by zero issue in sqrtm to fix change the algorithm for sqrtm from

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

to

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

or change line 79 of sqrtm.cc from

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

to

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

Also, ones(4) is a singular matrix so the inversion required in the integration method is still going to cause issues.

Why is the alternative algorithm better for Hermitian matrixes?

The Schur decomposition A=U*S*U' is the first step of the existing algorithm. S is normally upper triangular but it is diagonal for Hermitian Matrices, therefore, when A is Hermitian (symmetric for real matrixes) the Schur decomposition is equivalent to the eigendecomposition. Matrix functions are easily applied to diagonal matrices logm(A)=U*diag(log(diag(S)))*U'. This only applies when A is Hermitian and therefore the Schur decomposition results in S being diagonal.




(file #51578)

Anonymous
Thu 17 Jun 2021 04:45:29 PM UTC, comment #7: 

Thanks for the patch. (Unfortunately, it didn't apply for me.)

I pushed the part of it that adds the tolerance to the check for real negative values in a complex vector here:
https://hg.savannah.gnu.org/hgweb/octave/rev/fe06c183cc4e

Additionally, I pushed a patch that fixes e.g. `logm (zeros (3))` here:
https://hg.savannah.gnu.org/hgweb/octave/rev/4778b21b1386

I'm not sure about the impact of the other changes in your patch. Maybe those should better go to the default branch?
I tried `logm(ones(3))` in Matlab and it gives a different result:

>> logm(ones(3))
Warning: Principal matrix logarithm is not defined for A with nonpositive real eigenvalues. A non-principal matrix logarithm is returned.
> In logm (line 63)

ans =

-24.713572198920126 + 2.094395102393196i   13.566584297695865 - 1.047197551196598i   12.245600189892365 - 1.047197551196598i
 13.566584297695865 - 1.047197551196598i  -24.713572198920119 + 2.094395102393195i   12.245600189892361 - 1.047197551196597i
 12.245600189892365 - 1.047197551196598i   12.245600189892361 - 1.047197551196597i  -23.392588091116618 + 2.094395102393195i


I don't know if one or the other is "more correct" though.

Markus Mützel <mmuetzel>
Group administrator
Tue 15 Jun 2021 01:58:47 AM UTC, comment #6: 

Patch attached, incorporated tolerance check and 10x speedup for Hermitian matrices.

For random matrices, I'm just using A=s*rand(n) and A=s*randn(n) or complex(rand(n),randn(n)) of different sizes (n) and scales (s). A*=A' for Hermitian matrices or A+=A' for symmetric not nessaraly positive definate along with A=[21 17 6;-5 -1 -6;4 4 16]; or A=[5 4 2 1;0 1 -1 -1;-1 -1 3 0;1 1 -1 2]; for repeated eigenvalues. I'm placing the matrix at the top of the file and checking with the max elementwise norm at the bottom max(abs(expm(logm(A))-A)(:))

(file #51565)

Anonymous
Mon 14 Jun 2021 07:42:11 AM UTC, comment #5: 

Basing the tolerance on the absolute value of the largest eigenvalue sound reasonable to me. (Although I'm not sure this is the "right" thing.)
Can you provide an updated patch please?

Could you also please provide the code that you used for your tests with random matrices?

Markus Mützel <mmuetzel>
Group administrator
Wed 09 Jun 2021 03:12:04 PM UTC, comment #4: 

I'm not sure either. The tolerance should be the error in the Schur factorization which should be around the machine precision of each eigenvalue. I think it would be safe to base it on the absolute value of the largest eigenvalue. I'm open to other suggestions.

As the consequence of having a tolerance that is too lax is a few results will contain small complex values that are unlikely to change the final result by much and it should be obvious to the user that they come from numerical error whereas if the tolerance is too tight the complex values that need to be present will be removed and the answer will be wrong in a non-obvious way (this bug has been present for years and I'm the first person to notice it).

I gave it the tight tolerance below because it worked for a large number of random matrices. The error I typically see for complex components of real eigenvalues from the Schur factorization is normally very small.

Anonymous
Wed 09 Jun 2021 09:34:07 AM UTC, comment #3: 

These might be stupid questions. So, please be patient with me.
I'm not exactly certain if we should use the condition in comment #2.
If the matrix has eigenvalues that differ by many orders of magnitude, should only the largest determine the tolerance (even for the (very) small eigenvalues)?
Also: Should only the magnitude of the imaginary part be taken into account for the tolerance? Or should it be the absolute values of the eigenvalues?

Markus Mützel <mmuetzel>
Group administrator
Tue 08 Jun 2021 01:27:41 AM UTC, comment #2: 

The modification you made can still have the same issue but can be fixed by changing the line


real_neg_eigv = (real (eigv) < 0) & (imag (eigv) == 0);


to allow a small tolerance

abs_imag_eigv = abs (imag (eigv));
real_neg_eigv = (real (eigv) < 0) & (abs_imag_eigv <= 10 * eps (max (abs_imag_eigv))));


Anonymous
Mon 07 Jun 2021 09:53:56 AM UTC, comment #1: 

Thanks for your report and the detective work.

I pushed a slightly different fix to the stable branch here:
https://hg.savannah.gnu.org/hgweb/octave/rev/5928dd4f7a02

Marking as ready for test.

Markus Mützel <mmuetzel>
Group administrator
Sun 06 Jun 2021 06:58:10 AM UTC, original submission:  

The following shows a case where logm returns an incorrect result


A=[0.2510   1.2808  -1.2252
   0.2015   1.0766   0.5630
  -1.9769  -1.0922  -0.5831];
max(abs(expm(logm(A))-A)(:))

It's cause is the same as the following complex array not testing true in the following:

[-2 2+i -2-i]<0

The above appears to be equivalent to abs([-2 2+i -2-i])<0

logm can be fixed by changing line 83 from

real_eig = all (eigv >= 0);

to

real_eig = all (real (eigv) >= 0);

or perhaps fixing the logical check of real numbers in complex arrays?

Anonymous

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #51580:  bug60738logm@v2.patch added by None (3KiB - application/octet-stream)
file #51581:  logm.m added by None (6KiB - text/plain)
file #51579:  bug60738-logm-hermitian.patch added by mmuetzel (3KiB - application/octet-stream)
file #51578:  logm.m added by None (6KiB - text/plain)
file #51565:  bug60738logm.patch added by None (3KiB - 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 rik5 (Updated the item)
  • -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 10 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-02-23 rik5 StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2021-06-18 mmuetzel StatusIn Progress Ready For Test
    2021-06-18 None Attached File- Added bug60738logm@v2.patch, #51580
        Attached File- Added logm.m, #51581
    2021-06-18 mmuetzel Attached File- Added bug60738-logm-hermitian.patch, #51579
    2021-06-18 None Attached File- Added logm.m, #51578
    2021-06-15 None Attached File- Added bug60738logm.patch, #51565
    2021-06-14 mmuetzel StatusReady For Test In Progress
    2021-06-07 mmuetzel StatusNone Ready For Test

    Back to the top

    Powered by Savane 3.13-3230.
    Corresponding source code