bugGNU Octave - Bugs: bug #51246, pinv(0) different from Matlab

 
 

bug #51246: pinv(0) different from Matlab

Submitted by:  Guillaume <gyom>
Submitted on:  Fri 16 Jun 2017 12:26:04 PM UTC  
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Regression
Status:  Fixed Assigned to:  None
Originator Name:  Guillaume Open/Closed:  Closed
Release:  dev Operating System:  GNU/Linux

Add a New Comment (Rich Markup)
   

You are not logged in

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

 

( Jump to the original submission)

Fri 06 Apr 2018 05:47:45 PM UTC, comment #19: 

It looks like this has been resolved, but nevertheless I'd like to chime in.

In reference to comment 4:

"I would still open a discussion on Octave Maintainers list just to see if there is a good mathematical reason for returning either 0 or Inf. "

There is a good mathematical reason for returning pinv(0) = 0. Although there are different mathematical notions of generalized inverse, the most widely used is the Moore-Penrose inverse (https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse). For any matrix A, there exists a unique matrix B with the following four properties:

  1. A*B*A = A
  2. B*A*B = B
  3. A*B is Hermitian
  4. B*A is Hermitian

B is the Moore-Penrose inverse of A. It is easy to verify that the Moore-Penrose inverse of a zero matrix is a zero matrix.

If you look at the documentation of pinv in Matlab, this is precisely what is returned there. This is what Octave should return as well.

Travis Arnold <teerav42>
Fri 03 Nov 2017 05:10:35 PM UTC, comment #18: 

I checked in a change along with BIST tests here (http://hg.savannah.gnu.org/hgweb/octave/rev/d3dc76efb38b).

I also found that we were not calculating the tolerance correctly for single matrices.  Apparently there was a copy and paste error from the code for double matrices and we were using

std::numeric_limits<double>::epsilon ()

in calculations instead of the "float" version of epsilon.

Rik <rik5>
Project Administrator
Fri 03 Nov 2017 04:15:23 PM UTC, comment #17: 

>> pinv ([0 1e-15; 0 0], 1e-14)


ans =

     0     0
     0     0

>> pinv (zeros (2,2))


ans =

     0     0
     0     0

>> version


ans =

    '9.3.0.713579 (R2017b)'

>>

Looks OK.

Michael Godfrey <godfrey>
Project Member
Fri 03 Nov 2017 04:09:16 PM UTC, comment #16: 

Hmm, okay.  What about the following:

pinv ([0 1e-15; 0 0], 1e-14)

I think numbers below the tolerance should be rounded to zero and this will be equivalent to

pinv (zeros (2,2))

which should be a 2x2 zero matrix according to the results so far.

Rik <rik5>
Project Administrator
Fri 03 Nov 2017 03:55:32 PM UTC, comment #15: 

Last pinv (D) needs fixing:
D = diag ([0 0 0])

D =

     0     0     0
     0     0     0
     0     0     0

pinv (D)

ans =

     0     0     0
     0     0     0
     0     0     0

>> version


ans =

    '9.3.0.713579 (R2017b)'

>>
>> D = full (diag ([0 0 0]))


D =

     0     0     0
     0     0     0
     0     0     0

>> pinv (D)


ans =

     0     0     0
     0     0     0
     0     0     0

>>

Michael Godfrey <godfrey>
Project Member
Fri 03 Nov 2017 03:29:27 PM UTC, comment #14: 

@Michael: Thanks.  I've got a cset that addresses this issue.  I know it is a corner case, but I'm actually not pleased that Octave has been returning the wrong answer for zero vectors and scalars for many years.

A final test to run through Matlab for diagonal matrices

D = diag ([0 0 0])
pinv (D)

In Octave, this creates a special Diagonal matrix and the result is

octave:4> D = diag ([0 0 0])
D =

Diagonal Matrix

   0   0   0
   0   0   0
   0   0   0

octave:5> pinv (D)
ans =

Diagonal Matrix

   0   0   0
   0   0   0
   0   0   0

But if I use full matrices, I get NaN

octave:6> D = full (diag ([0 0 0]))
D =

   0   0   0
   0   0   0
   0   0   0

octave:7> pinv (D)
ans =

   NaN   NaN   NaN
   NaN   NaN   NaN
   NaN   NaN   NaN

So what does Matlab do?

Rik <rik5>
Project Administrator
Fri 03 Nov 2017 10:36:29 AM UTC, comment #13: 

>> pinv (complex (0,0))

pinv ([complex(0,0), complex(0,0), complex(0,0)])

ans =

     0

ans =

     0
     0
     0

>> version


ans =

    '9.3.0.713579 (R2017b)'

>>

Michael Godfrey <godfrey>
Project Member
Fri 03 Nov 2017 02:46:19 AM UTC, comment #12: 

I've been looking at resolving the issues with inv(0) and am now  re-visiting this bug.  What does Matlab return for

pinv (complex (0,0))
pinv ([complex(0,0), complex(0,0), complex(0,0)])

Rik <rik5>
Project Administrator
Tue 20 Jun 2017 09:37:18 AM UTC, comment #11: 

So Matlab and Octave agree on the 2x2 matrices. Is the patch of comment #9 ready to be pushed or are there other issues to take care of?

Guillaume <gyom>
Sat 17 Jun 2017 10:11:53 AM UTC, comment #10: 

>> inv (zeros (2,2))
Warning: Matrix is singular to working precision.

ans =

   Inf   Inf
   Inf   Inf

>> inv (repmat (inf, [2,2]))
Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate.
RCOND = NaN.

ans =

   NaN   NaN
   NaN   NaN

>> version

ans =

8.4.0.150421 (R2014b)

Marco Caliari <caliari>
Project Member
Fri 16 Jun 2017 06:23:08 PM UTC, comment #9: 

Attached is a quick hack (inverse.diff) which applies only to double matrices.  It explicitly checks for the case of a scalar which is either 0 or Inf and returns the correct output value.

diff -r 5e0acf41a9e1 -r 26ead836cbea liboctave/array/dMatrix.cc
--- a/liboctave/array/dMatrix.cc        Fri Jun 16 09:03:10 2017 -0700
+++ b/liboctave/array/dMatrix.cc        Fri Jun 16 11:18:36 2017 -0700
@@ -639,10 +639,19 @@ Matrix::inverse (MatrixType& mattype, oc
       if (! mattype.ishermitian ())
         ret = finverse (mattype, info, rcon, force, calc_cond);

-      if ((calc_cond || mattype.ishermitian ()) && rcon == 0.
-          && ! (numel () == 1))
-        ret = Matrix (rows (), columns (),
-                      octave::numeric_limits<double>::Inf ());
+      if ((calc_cond || mattype.ishermitian ()) && rcon == 0.)
+        {
+          if (numel () == 1)
+            {
+              if (octave::math::isinf (this->xelem (0)))
+                ret = Matrix (1, 1, 0.0);
+              else
+                ret = Matrix (1, 1, octave::numeric_limits<double>::Inf ());
+            }
+          else
+            ret = Matrix (rows (), columns (),
+                          octave::numeric_limits<double>::Inf ());
+        }
     }

   return ret;

What does Matlab return for these two cases

inv (zeros (2,2))
inv (repmat (inf, [2,2]))

Rik <rik5>
Project Administrator
Fri 16 Jun 2017 05:36:24 PM UTC, comment #8: 

Any fix also needs to respect bug #49690 where the inverse of Inf should be 0 for Matlab compatibility.  My guess is that is where the problem in inverse for the input 0 has crept in.  See cset 22154.

changeset:   22154:1d242ae72240
user:        Rik <rik@octave.org>
date:        Tue Nov 29 16:46:35 2016 -0800
files:       libinterp/corefcn/inv.cc liboctave/array/CMatrix.cc liboctave/array/dMatrix.cc liboctave/array/fCMatrix.cc liboctave/array/fMatrix.cc
description:
Return 0 for special case of scalar Inf input to inverse (bug #49690).

* inv.cc (Finv): Add BIST tests for new behavior.

* CMatrix.cc, fCMatrix.cc, dMatrix.cc, fMatrix.cc (::inverse):
If function has had to calculate the reciprocal condition number, and
the rcond value is zero, then return 0 if the input was a scalar, or
a matrix of all Inf values.

Rik <rik5>
Project Administrator
Fri 16 Jun 2017 05:00:48 PM UTC, comment #7: 

The code for the pseudo_inverse of a double matrix is in liboctave/array/dMatrix.cc.  It looks like it calculates the svd, and then also calculates an inverse.  My guess is the inverse code is where the problem creeps in.

Rik <rik5>
Project Administrator
Fri 16 Jun 2017 04:47:32 PM UTC, comment #6: 

The pseudoinverse of 0 is 0 (wikipedia, for instance). And the pseudoinverse of a null vector is the transposed null vector (this answers comment #5).
inv(0) was Inf on 4.2.0 and this is correct.

Marco Caliari <caliari>
Project Member
Fri 16 Jun 2017 04:39:21 PM UTC, comment #5: 

pinv is a weird function. 

 pinv ([0 0])
ans =

   Inf
   NaN

What does Matlab do for this test case?

Rik <rik5>
Project Administrator
Fri 16 Jun 2017 04:25:30 PM UTC, comment #4: 

I would still open a discussion on Octave Maintainers list just to see if there is a good mathematical reason for returning either 0 or Inf.

And from the results of your test it looks like there are two Matlab incompatibilities.

Matlab
--------
inv (0) = Inf
pinv (0) = 0

Octave
--------
inv (0) = 0
pinv (0) = Inf

Rik <rik5>
Project Administrator
Fri 16 Jun 2017 04:09:37 PM UTC, comment #3: 

In Matlab:

>> inv(0)
Warning: Matrix is singular to working precision.
ans =
   Inf

I agree this is a design choice but I would say that 100% compatible with Matlab is required for functions as fundamental as inv & pinv.

Guillaume <gyom>
Fri 16 Jun 2017 03:50:23 PM UTC, comment #2: 

Leaving aside what Matlab does for the moment, what would you like to see for the inverse of a scalar?

In general, the inverse would just be 1/x for a scalar.  So,

pinv (5)
ans =  0.20000

If I follow this analogy, pinv(0) should return 1/0

octave:2> 1/0
warning: division by zero
ans = Inf

which is what the function seems to be doing.

Of course, maybe it is better to be Matlab compatible for this case.  Or maybe pinv should follow inv. 

octave:5> inv (0)
warning: matrix singular to machine precision
ans = 0

There is definitely room for interpretation here, rather than being absolutely the wrong answer.  Maybe you should post on the Octave Maintainer's list to gather some other opinions.

Rik <rik5>
Project Administrator
Fri 16 Jun 2017 02:20:48 PM UTC, comment #1: 

I confirm and changed to regression, since it was 0 on 4.0.2. Moreover

octave:1> pinv ([1, 0])
ans =

   1
   0

is correct on dev.

Marco Caliari <caliari>
Project Member
Fri 16 Jun 2017 12:26:04 PM UTC, original submission:  

pinv(0) returns Inf in Octave and 0 in Matlab.

Guillaume <gyom>

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #40928:  inverse.diff added by rik5 (1KiB - text/x-patch)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by teerav42 (Posted a comment)
  • -email is unavailable- added by godfrey (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by caliari (Posted a comment)
  • -email is unavailable- added by gyom (Submitted the item)
  •  

    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 6 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2017-11-03 rik5 StatusIn Progress => Fixed
        Open/ClosedOpen => Closed
    2017-06-16 rik5 Attached File- => Added inverse.diff, #40928
    2017-06-16 rik5 StatusNone => In Progress
        Summarypinv(0) => pinv(0) different from Matlab
    2017-06-16 caliari Item GroupMatlab Compatibility => Regression

    Back to the top


    Powered by Savane 3.5