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 11 Feb 2014 10:06:20 AM UTC  
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 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 18 Feb 2014 09:51:50 PM UTC, 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 18 Feb 2014 08:23:29 AM UTC, 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 17 Feb 2014 06:28:47 PM UTC, 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 17 Feb 2014 10:25:42 AM UTC, 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 15 Feb 2014 10:43:30 PM UTC, 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 13 Feb 2014 05:54:19 PM UTC, comment #2:

Confirmed on the default branch as well.

Mike Miller <mtmiller>
Project Administrator
Thu 13 Feb 2014 01:47:02 PM UTC, 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 11 Feb 2014 10:06:20 AM UTC, 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 15 Feb 2014 10:43:30 PM UTCrik5StatusConfirmed=>Fixed
    Thu 13 Feb 2014 05:54:19 PM UTCmtmillerCategoryNone=>Octave Function
      Item GroupNone=>Incorrect Result

    Back to the top

    Powered by Savane 3.1-cleanup1