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 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 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.

E.g.

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 Member
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.

pinv([q],2)

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.

http://octave.sourceforge.net/doxygen/html/de/d18/fCDiagMatrix_8cc_source.html#l00390

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])
pinv(q,2)

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):
   
   
Comment:
   

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
      Open/ClosedOpen=>Closed
    Thu 13 Feb 2014 05:54:19 PM UTCmtmillerCategoryNone=>Octave Function
      Item GroupNone=>Incorrect Result
      StatusNone=>Confirmed
      Release3.2.4=>dev

    Back to the top


    Powered by Savane 3.1-cleanup