bugGNU Octave - Bugs: bug #60786, mpower: negative scalar to power...

 
 

bug #60786: mpower: negative scalar to power of matrix containing non-integers returns nan

Submitter:  None
Submitted:  Wed 16 Jun 2021 01:48:43 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Need Info Assigned to:  None
Originator Name:  Originator Email:  -email is unavailable-
Open/Closed:  * Open Release:  * dev
Operating System:  * Any Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Sat 03 Jul 2021 07:04:16 PM UTC, comment #21: 

@jwe: shall I back out the one changeset?  Or does it matter given that a bigger overhaul is required?

Rik <rik5>
Group administrator
Sat 03 Jul 2021 07:01:53 PM UTC, comment #20: 

Testing code from comment #17 produces different results on my machine.  Instead of jwe's results of


1
(nan,nan)
(-nan,-nan)
(nan,nan)


I get


1
(1,0)
(-nan,-nan)
(1,0)


This is with a C++ stdlib from Kubuntu 18.04.


ii  libstdc++6:amd64                                8.4.0-1ubuntu1~18.04                                amd64        GNU Standard C++ Library v3


The fact that results seem to change with changing library versions is not a good sign, and points to implementing our own wrappers to get the results we want.

Rik <rik5>
Group administrator
Wed 30 Jun 2021 06:15:03 PM UTC, comment #19: 

@rik: Yes, definitely a bummer.  I'm tempted to use our own "oct-math.h" and "oct-complex.h" that would provide all the math functions that we need with the behavior we want.  Then we would also not allow the standard header files to be included directly in Octave's source files and we could at least ensure consistency in Octave functions.  But of course that doesn't "fix" the behavior of other dependencies that we have.

If we don't override the complex pow function to do what we want we will have to check all places where it is used.  I've been looking at the element-wise operations in sparse-xpow.cc and it first checks whether conversion to complex is needed for any of the calculations, and if so just uses complex pow for all.  But then that can produce these "bad" results:


octave> x = sparse ([-2, 0, -2]);
octave> y = sparse ([-2, 0, 0.25]);
octave> x .^ y
ans =

Compressed Column Sparse (rows = 1, cols = 3, nnz = 3 [100%])

  (1, 1) ->  2.5000e-01 + 6.1232e-17i
  (1, 2) ->  1 + 0i
  (1, 3) ->  0.8409 + 0.8409i


John W. Eaton <jwe>
Group administrator
Wed 30 Jun 2021 04:40:01 PM UTC, comment #18: 

@jwe: Bummer.  This is why I made many small changesets so each one could be reversed as necessary.  I'll take a look later today or tomorrow.  If the std::pow template really doesn't work well with Complex arguments with imaginary parts then I think we will need to duplicate the testing logic in many more of the templates in xpow.cc.  In other words, reversing the changeset is probably not enough.  We will need to examine every instance where a Complex value is used and potentially add a test.

Rik <rik5>
Group administrator
Wed 30 Jun 2021 02:49:03 AM UTC, comment #17: 

RE: comment #14 and the changeset http://hg.savannah.gnu.org/hgweb/octave/rev/004e63a7a50, I don't think we get the expected result if we just use std::pow.  For example, with the following program


#include <iostream>
#include <complex>

int main (void)
{
  std::complex<double> c (0.0, 0.0);
  double d = 0.0;

  std::cerr << std::pow (d, d) << std::endl;
  std::cerr << std::pow (d, c) << std::endl;
  std::cerr << std::pow (c, d) << std::endl;
  std::cerr << std::pow (c, c) << std::endl;

  return 0;
}


I see


1
(nan,nan)
(-nan,-nan)
(nan,nan)


so it doesn't appear that std::pow for complex values with zero imaginary part gives the same result as converting to double and then calling std::pow (double, double).

The differences aren't limited to just the complex (0, 0) case.  If I change the program above to set c = (-2, 0) and d = -2, I see:


0.25
(0.25,6.12323e-17)
(0.25,6.12323e-17)
(0.25,6.12323e-17)


I've been looking at bug #60846 and noticed that complex (0, 0) ^ complex (0, 0) produced NaNs instead of 1 as 0 ^ 0 does.  We have this problem for the scalar operators and for element-wise operators for both full and sparse matrices, so this is nothing new but it does appear to be incompatible behavior with Matlab.

John W. Eaton <jwe>
Group administrator
Mon 28 Jun 2021 07:03:15 PM UTC, comment #16: 

I checked in a series of changesets that generally improve things.

See:

1. http://hg.savannah.gnu.org/hgweb/octave/rev/212322c0e44e
2. http://hg.savannah.gnu.org/hgweb/octave/rev/0474c3a27d16
3. http://hg.savannah.gnu.org/hgweb/octave/rev/ef138252af51

I think this is good enough now for this bug report so I am going to mark the issue as fixed.

There are still opportunities for optimization:

1. Use templates to more concisely cover Matrix, FloatMatrix, NDArray, FloatNDArray cases rather than having lots of code duplication.

2. Don't call octave_quit() for every iteration of elem_xpow().  For example, in some of the library routines in liboctave/Array we check for interrupts only after 4 executions of the innermost loop.  That same syntax could be copied to xpow.cc.  In addition, even 4 loop iterations is probably less than a millisecond execution time so the number 4 could probably be increased.

3. The code uses operator indexing which is the equivalent of calling elem().  However, the code knows the exact size of matrices so these could be re-written to use xelem() for a possible performance benefit.

Rik <rik5>
Group administrator
Mon 28 Jun 2021 05:17:13 PM UTC, comment #15: 

Finished step #2, better error reporting, by moving call to EIG outside of try/catch block in this cset: http://hg.savannah.gnu.org/hgweb/octave/rev/e4c65b8276a6.

Interestingly, there was only one case of this that needed to be fixed.  It looks like an accident that xpow (double, Matrix) was the only case where this code pattern was used.

Rik <rik5>
Group administrator
Mon 28 Jun 2021 05:07:46 PM UTC, comment #14: 

I completed point #3 (eliminating check for imag==0.0 and just use C++ library to figure out correct template) in this cset: http://hg.savannah.gnu.org/hgweb/octave/rev/004e63a7a501.

Will now complete step #2.

Rik <rik5>
Group administrator
Sat 26 Jun 2021 08:01:43 AM UTC, comment #13: 

I second your proposal to move the EIG call into the try/catch block, your example looks convincing.

A.R. Burgers <arb>
Thu 24 Jun 2021 09:07:47 PM UTC, comment #12: 

Next question is step 2 and moving the eigenvalue decomposition into OR out of the try/catch block.  The error messages will be different.  Here is a sample session to show those differences.


octave:14> mpower (-2, x)
error: Failed to diagonalize matrix while calculating matrix exponential
octave:15> mpower (-2, complex (x))
error: EIG: matrix contains Inf or NaN values


In the first case, the eig call is within the try/catch block.  In the second case, it is outside the try/catch block.

Does anyone have a preference?  I sort of appreciate the extra specificity of the error message from eig() when it fails.  In this case, I purposefully did have a NaN value in the matrix 'x'.

Rik <rik5>
Group administrator
Thu 24 Jun 2021 09:00:38 PM UTC, comment #11: 

I finished step 1 and checked in a changeset here: http://hg.savannah.gnu.org/hgweb/octave/rev/f54cfd60725e.

Rik <rik5>
Group administrator
Wed 23 Jun 2021 04:52:10 PM UTC, comment #10: 

In comment #7 I wrote:

> Additionally, would it make sense to have early exits when a=0.


That is not a good idea as for example 0^0=1 and 0^1=0, so you need to get the eigenvalues to get the answer correct for a=0

A.R. Burgers <arb>
Wed 23 Jun 2021 04:45:46 PM UTC, comment #9: 

@jwe, @arb: My plan right now is to proceed in a series of small steps/changesets.  The idea being that any one change can be quickly reversed by backing out the appropriate changeset if it is somehow problematic.

Steps

1. Fix output of NaN for mixed scalar/complex matrix combinations
2. Use common try/catch block with eig calculation outside block
3. Eliminate test for complex value with imag=0.0 and rely on C++ library templates to correctly process this case.

4. Templating.  This would be the final changeset, but I have other coding projects that are higher priority for me.  If someone else is interested in tackling this that would be great.

Rik <rik5>
Group administrator
Wed 23 Jun 2021 04:02:19 PM UTC, comment #8: 

Rik: This code has origins even earlier than I thought.  The first commit in the hg archive dated August 1993 has this code for the scalar x matrix case:


tree_constant
xpow (double a, Matrix& b)
{
  tree_constant retval;

  int nr = b.rows ();
  int nc = b.columns ();

  if (nr == 0 || nc == 0 || nr != nc)
    error ("for x^A, A must be square");
  else
    {
      EIG b_eig (b);
      ComplexColumnVector lambda (b_eig.eigenvalues ());
      ComplexMatrix Q (b_eig.eigenvectors ());

      for (int i = 0; i < nr; i++)
        {
          Complex elt = lambda.elem (i);
          if (imag (elt) == 0.0)
            lambda.elem (i) = pow (a, real (elt));
          else
            lambda.elem (i) = pow (a, elt);
        }
      ComplexDiagMatrix D (lambda);

      ComplexMatrix result = Q * D * Q.inverse ();
      retval = tree_constant (result);
    }

  return retval;
}


so the check for real eigenvalues and the conversion to real has been in the code from the beginning.  If this isn't the part of the code you were asking about yesterday, then please let me know which function it was because I'm curious now to know if there was a later change that introduced the complex to real conversion that you were asking about.

And, if there is a better way to do this now using templates, I'm all in favor of fixing it, but only if someone has time and is interested in working on it, of course.

BTW, the hg archive apparently doesn't have ALL changes since it begins with a date of 1993.  I guess that's the earliest change that was in the CVS archive that we converted to hg.  I'm fairly sure I used some form of version control from the very beginning though (probably RCS), so it is a little surprising that there are no earlier changes recorded since I started working on Octave in February of 1992.

John W. Eaton <jwe>
Group administrator
Wed 23 Jun 2021 06:54:41 AM UTC, comment #7: 

Here is test script:


a_vals = [2, -2, 0];
a_vals = [0];
a_vals = [-2];
for a_s = a_vals
                       a_fs = single(a_s);
  b_m  = 0.5*eye(2);   b_fm = single(b_m);
  b_cm = complex(b_m); b_fcm = complex(b_fm);
  b_sm = sparse(b_m);  b_scm = sparse(b_cm);

  % powers of float/double matrices
  pow_s_m = a_s ^ b_m
  pow_s_fm = a_s ^ b_fm
  pow_fs_m = a_fs ^ b_m
  pow_fs_fm = a_fs ^ b_fm

  % powers of float/double complex matrices
  pow_s_cm = a_s ^ b_cm
  pow_s_fcm = a_s ^ b_fcm
  pow_fs_cm = a_fs ^ b_cm
  pow_fs_cfm = a_fs ^ b_fcm

  % powers of sparse  matrices
  pow_s_sm = a_s ^ b_sm
  pow_s_scm = a_s ^ b_scm
  if 0
    % not implemented
    pow_fs_sm = a_fs ^ b_sm
  end
end


which produces this after your patch, showing the sparse and real matrices have been fixed, but not the complex ones.


pow_s_m =
   0.0000 + 1.4142i        0 +      0i
        0 +      0i   0.0000 + 1.4142i

pow_s_fm =
  -0.0000 + 1.4142i        0 +      0i
        0 +      0i  -0.0000 + 1.4142i

pow_fs_m =
  -0.0000 + 1.4142i        0 +      0i
        0 +      0i  -0.0000 + 1.4142i

pow_fs_fm =
  -0.0000 + 1.4142i        0 +      0i
        0 +      0i  -0.0000 + 1.4142i

pow_s_cm =
   NaN - NaNi   NaN - NaNi
   NaN - NaNi   NaN - NaNi

pow_s_fcm =
   NaN - NaNi   NaN - NaNi
   NaN - NaNi   NaN - NaNi

pow_fs_cm =
   NaN - NaNi   NaN - NaNi
   NaN - NaNi   NaN - NaNi

pow_fs_cfm =
   NaN - NaNi   NaN - NaNi
   NaN - NaNi   NaN - NaNi

pow_s_sm =
   0.0000 + 1.4142i        0 +      0i
        0 +      0i   0.0000 + 1.4142i

pow_s_scm =
   0.0000 + 1.4142i        0 +      0i
        0 +      0i   0.0000 + 1.4142i


also there was already an inconsistency in xpow whether EIG is called in- or outside the try block. In xpow (double a, const Matrix& b) it is called inside, otherwise outside. On first glance outside seems OK, since what fails and needs to be caught is the Q.inverse.

I wonder whether the separate code path for a>0 is necessary for the real matrix case, as later on for a>0 the real part of the result is taken anyway before returning.

Additionally, would it make sense to have early exits when a=0.


A.R. Burgers <arb>
Tue 22 Jun 2021 08:57:33 PM UTC, comment #6: 

Do you have a test case for scalars raised to complex powers that I can use for testing?

I asked about the complicated structure of xpow.cc at an Octave Developer's Meeting today.  The reason is that this was coded back in the late '90s before templating was widespread.  Simplification is possible, but the easiest thing right now is just to fix the code, rather than rewrite this file.  Maybe sometime in the future that can happen.

I pushed a changeset here (http://hg.savannah.gnu.org/hgweb/octave/rev/540f25090412) that fixes the example cited in comment #0.

Rik <rik5>
Group administrator
Mon 21 Jun 2021 08:55:51 PM UTC, comment #5: 

scalars raised to complex matrix powers have the issue too.
There seems to be a bit of code duplication in xpow.cc Would be nice if some form of templating could be used?

Like ComplexMatrix <float> or EIG <float>/EIG <double> instead of FloatComplexMatrix and FloatEIG

A.R. Burgers <arb>
Mon 21 Jun 2021 03:44:25 AM UTC, comment #4: 

I should have posted sooner.  I used gdb to resolve the issue down to the same call to std::pow.  I've coded a patch, but haven't pushed it yet because I wanted to ask some questions on the Octave Discourse site about potentially simplifying the code.

Rik <rik5>
Group administrator
Sun 20 Jun 2021 08:36:08 PM UTC, comment #3: 

The float-float version must be changed accordingly


diff -r c75127492941 libinterp/corefcn/xpow.cc
--- a/libinterp/corefcn/xpow.cc Fri Jun 18 08:26:54 2021 +0200
+++ b/libinterp/corefcn/xpow.cc Sun Jun 20 22:25:20 2021 +0200
@@ -124,7 +124,7 @@

   if (nr != nc)
     err_nonsquare_matrix ();
-
+  std::complex<double> pow_arg(a, 0);
   try
     {
       EIG b_eig (b);
@@ -136,9 +136,9 @@
         {
           Complex elt = lambda(i);
           if (std::imag (elt) == 0.0)
-            lambda(i) = std::pow (a, std::real (elt));
+            lambda(i) = std::pow (pow_arg, std::real (elt));
           else
-            lambda(i) = std::pow (a, elt);
+            lambda(i) = std::pow (pow_arg, elt);
         }
       ComplexDiagMatrix D (lambda);

@@ -1558,10 +1558,13 @@
   if (nr != nc)
     err_nonsquare_matrix ();

-  FloatEIG b_eig (b);
+  std::complex<float> pow_arg(a, 0);

   try
     {
+
+      FloatEIG b_eig (b);
+
       FloatComplexColumnVector lambda (b_eig.eigenvalues ());
       FloatComplexMatrix Q (b_eig.right_eigenvectors ());

@@ -1569,9 +1572,9 @@
         {
           FloatComplex elt = lambda(i);
           if (std::imag (elt) == 0.0)
-            lambda(i) = std::pow (a, std::real (elt));
+            lambda(i) = std::pow (pow_arg, std::real (elt));
           else
-            lambda(i) = std::pow (a, elt);
+            lambda(i) = std::pow (pow_arg, elt);
         }
       FloatComplexDiagMatrix D (lambda);


A.R. Burgers <arb>
Sun 20 Jun 2021 05:51:13 PM UTC, comment #2: 

this patch fixes the case. Not sure this is the best fix, but isolates the issue. The issue is std::pow(a, ...) produces NaN when a <= 0.


diff -r c75127492941 libinterp/corefcn/xpow.cc
--- a/libinterp/corefcn/xpow.cc Fri Jun 18 08:26:54 2021 +0200
+++ b/libinterp/corefcn/xpow.cc Sun Jun 20 19:46:16 2021 +0200
@@ -124,7 +124,7 @@

   if (nr != nc)
     err_nonsquare_matrix ();
-
+  std::complex<double> pow_arg(a, 0);
   try
     {
       EIG b_eig (b);
@@ -136,9 +136,9 @@
         {
           Complex elt = lambda(i);
           if (std::imag (elt) == 0.0)
-            lambda(i) = std::pow (a, std::real (elt));
+            lambda(i) = std::pow (pow_arg, std::real (elt));
           else
-            lambda(i) = std::pow (a, elt);
+            lambda(i) = std::pow (pow_arg, elt);
         }
       ComplexDiagMatrix D (lambda);


A.R. Burgers <arb>
Wed 16 Jun 2021 08:40:26 PM UTC, comment #1: 

Confirmed.  The result should be


[0+sqrt(2)*i, 0+0i;
 0+0i       , 0+sqrt(2)*i]


Setting "Release" field to "dev" since the fix is likely to be involved and will have to be made on the development branch.

This is also a compatibility issue as Matlab gets the correct answer.

Rik <rik5>
Group administrator
Wed 16 Jun 2021 01:48:43 PM UTC, original submission:  

If a negative scalar A is raised to a matrix power A^P that contains decimals (tried integers for a simple example first and it worked) mpower returns a matrix of NaN. A Taylor series approximation of A^P and expm(log(A)*P) return the expected result.


A=-2;
P=[.5 0;0 .5];
#Taylor expected result
v=zeros(2);
for k=0:10
  v+=(P*log(A))^k/factorial(k);
end
v
#explog expected result
expm(log(A)*P)
#mpower incorrect all NaN
A^P
mpower(A,P)


Anonymous

 

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

Attach Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by siko1056 (Updated the item)
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by arb (Posted a comment)
  • -email is unavailable- added by rik5 (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 8 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2021-06-30 siko1056 StatusFixed Need Info
        Open/ClosedClosed Open
    2021-06-28 rik5 StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2021-06-22 rik5 StatusConfirmed Ready For Test
    2021-06-16 rik5 Summarympower: negative scalar to power of matrix containing decimals returns nan mpower: negative scalar to power of matrix containing non-integers returns nan
    2021-06-16 rik5 StatusNone Confirmed
        Release6.2.0 dev

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code