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

 
 

bug #51246: pinv(0) different from Matlab

Submitter:  Guillaume <gyom>
Submitted:  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
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

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>
Group 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>
Group 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>
Group 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>
Group 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>
Group 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>
Group 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>
Group 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>
Group 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>
Group 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>
Group 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>
Group 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>
Group 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>
Group 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>
Group 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>
Group 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>
Group 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)
  •  

    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
    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.13-f8d8.
    Corresponding source code