# bug #51246: pinv(0) different from Matlab

 Submitter: Guillaume 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

## Discussion

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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
Fri 16 Jun 2017 12:26:04 PM UTC, original submission:

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

Guillaume <gyom>

## Attached Files

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

## Dependencies

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.