bugGNU Octave - Bugs: bug #41546, warning ("pinv: tol is...


bug #41546: warning ("pinv: tol is ignored for diagonal matrices");

Submitted by:  Ceral Paquet <octavebugs>
Submitted on:  Tue Feb 11 10:06:20 2014  
Category: Octave FunctionSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Incorrect Result
Status: FixedAssigned to: None
Originator Name: Open/Closed: Closed
Release: devOperating System: GNU/Linux

Add a New Comment (Rich MarkupRich Markup):

You are not logged in

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


(Jump to the original submission Jump to the original submission)

Tue Feb 18 21:51:50 2014, comment #7:

I was looking at the SVD code and realized there's another issue with diagonal pinv: no checking for Inf and NaN values.

BTW the "is_perm_matrix" has no handling of tol and no warning either. Perhaps it would be best to remove this code and let the general case handle it?

97 else if (arg.is_perm_matrix ())
98 {
99 retval = arg.perm_matrix_value ().inverse ();
100 }

Ceral Paquet <octavebugs>
Tue Feb 18 08:23:29 2014, comment #6:

Watch out for the sign. As it is stands, the new code uses retval.elem(i, i) = 1.0 / val (which is the absolute value of the diagonal).

Please consider using the same handling as for the general case, i.e.

if (tol <= 0.0)
tol = nc * sigma.elem (0) * std::numeric_limits<double>::epsilon ();

otherwise there will be slight differences for diagonal matrices.

Just for kicks, I checked to see if SVD with nargout=1 has special handling for diagonal matrices: it doesn't. Perhaps it should (it's just the sorted absolute values of the diagonal).

Ceral Paquet <octavebugs>
Mon Feb 17 18:28:47 2014, comment #5:

Given that diagonal matrices are already going to be very fast to invert, just a for loop over the elements along the diagonal and a simple division operation, I think it is okay to have an extra explicit test for the element being zero and keep compatibility with the regular pinv.

I made that change in this changeset (http://hg.savannah.gnu.org/hgweb/octave/rev/16b0cd465ecd).

Rik <rik5>
Project Administrator
Mon Feb 17 10:25:42 2014, comment #4:

I just noticed a patch got sent to fix a bug I reported last week.

Awesome - thanks. However there's a small gotcha in the new code:

+ if (std::abs (elem (i, i)) < tol)
+ retval.elem (i, i) = 0.0;
+ else
retval.elem (i, i) = 1.0 / elem (i, i);

When tol==0 and elem(i,i)==0, the code will divide by zero. Changing < to <= will remedy the problem but makes the behaviour inconsistent with non-diagonal pinv function.


q=diag([3 2 1]);
pinv([q],2) # uses non-diag code
pinv(q,2) # uses diag code

To copy from the non-diagonal pinv, we could do [url=http://octave.sourceforge.net/doxygen/html/dd/d4d/dMatrix_8cc_source.html#l00854this[/url]]:

868 if (tol <= 0.0)
869 {
873 tol = nc * sigma.elem (0) * std::numeric_limits<double>::epsilon ();
874 }

where sigma.elem (0) in this case is the max abs of the diagonal elements.

Ceral Paquet <octavebugs>
Sat Feb 15 22:43:30 2014, comment #3:

This was a pretty easy fix, although long. I made the change on the development branch in this changeset (http://hg.savannah.gnu.org/hgweb/octave/rev/c5a101de2d88). To access the fix you will need to build from Mercurial sources, or wait until the 4.2 release which is a long ways away.

Rik <rik5>
Project Administrator
Thu Feb 13 17:54:19 2014, comment #2:

Confirmed on the default branch as well.

Mike Miller <mtmiller>
Project Administrator
Thu Feb 13 13:47:02 2014, comment #1:

Interestingly pinv([q],2) gives the right answer. i.e. the extra square braces cause the diagonal property to be lost.


ans =

0.33333 0.00000 0.00000
0.00000 0.50000 0.00000
0.00000 0.00000 1.00000

I can see in the source code what needs to be done but currently I don't know how to go about it.


The pseudo_inverse(void) needs to be changed to accept tol as an argument and the loop needs to be modified a little:

400 if (abs(elem (i, i)) > tol)
401 retval.elem (i, i) = 1.0f / elem (i, i);
402 else
403 retval.elem (i, i) = 0.0f;

Ditto for the complex and double precision versions.

Regarding the loss of diagonality of [q], I have no idea what to do!

Ceral Paquet <octavebugs>
Tue Feb 11 10:06:20 2014, original submission:

This functionality differs from MATLAB, which will happily apply a singular value tolerance to diagonal matrices (as it should!).

q=diag([3 2 1])

warning: pinv: tol is ignored for diagonal matrices
ans =

Diagonal Matrix

0.33333 0 0
0 0.50000 0
0 0 1.00000

Ceral Paquet <octavebugs>


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

Attach File(s):

No files currently attached


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -unavailable- added by rik5 (Posted a comment)
  • -unavailable- added by mtmiller (Posted a comment)
  • -unavailable- added by octavebugs (Submitted the item)

    Do you think this task is very important?
    If so, you can click here to add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.


    Please enter the title of George Orwell's famous dystopian book (it's a date):



    Follow 6 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Sat Feb 15 22:43:30 2014rik5StatusConfirmed=>Fixed
    Thu Feb 13 17:54:19 2014mtmillerCategoryNone=>Octave Function
      Item GroupNone=>Incorrect Result

    Back to the top

    Powered by Savane 3.1-cleanup1