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

 Submitted by: Ceral Paquet Submitted on: 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

You are not logged in

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

Confirmed on the default branch as well.

Mike Miller <mtmiller>
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.

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

Do you think this task is very important?
This task has 0 encouragements so far.

Only project members can vote.