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

 
 

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

Submitter:  Ceral Paquet <octavebugs>
Submitted:  Tue 11 Feb 2014 10:06:20 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  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

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>
Group 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>
Group administrator
Thu 13 Feb 2014 05:54:19 PM UTC, comment #2: 

Confirmed on the default branch as well.

Mike Miller <mtmiller>
Group 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 Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by mtmiller (Posted a comment)
  • -email is unavailable- added by octavebugs (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
    2014-02-15 rik5 StatusConfirmed Fixed
        Open/ClosedOpen Closed
    2014-02-13 mtmiller CategoryNone Octave Function
        Item GroupNone Incorrect Result
        StatusNone Confirmed
        Release3.2.4 dev

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code