bugGNU Octave - Bugs: bug #31773, numerical instability in the pinv...

 
 

bug #31773: numerical instability in the pinv (Moore-Penrose inverse) function

Submitter:  Sven Serneels <svens>
Submitted:  Fri 03 Dec 2010 09:04:31 PM UTC
   
 
Category:  Configuration and Build System Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Invalid / Not an Octave Bug Assigned to:  None
Originator Name:  Sven Serneels Open/Closed:  * Closed
Release:  * 3.2.4 Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Thu 16 Dec 2010 09:30:55 PM UTC, comment #4: 

Your computation involves the pseudo-inverse of a nearly singular
matrix and that's about as numerically sensitive as it gets. In
fact, if your sigma values are off in the fourteenth decimal or
so, that's sufficient to get completely different results,
whether you believe it or not. The problem is that the pinv
function has to decide which singular values are zero and which
are not. The default numerical tolerance is quite small and since
the non-zero singular values are inverted, any value which should
mathematically be zero but numerically falls just outside the
tolerance will introduce huge deviations. For a simple example,
compare the results of

octave:2> pinv(diag([1 1e-15]))
ans =

   1.0000e+00   0.0000e+00
   0.0000e+00   1.0000e+15

octave:3> pinv(diag([1 1e-16]))
ans =

   1   0
   0   0

If you want your code to be more numerically robust I propose
that you make use of the second argument to pinv and set it to
1e-10 or so. E.g.

octave:4> pinv(diag([1 1e-15]), 1e-10)
ans =

   1   0
   0   0

Incidentally, with the data you provided I see the same behavior in
Matlab as in Octave.

Gunnar Farnebäck <gf>
Mon 13 Dec 2010 06:06:23 PM UTC, comment #3: 

Please find in the attachment the data producing the error reported.

Note as well that sigmas and sigmac should be, according to the underlying mathematical equations, identic.

I do not believe that 5th and later decimal differences do lead to results differing in this order of magnitude. This is a bug and not some numerical instability acceptable within some limits.

NB - The payware alternative (Matlab) does produce identic values for all results.

(file #22196)

Anonymous
Sun 12 Dec 2010 08:30:56 PM UTC, comment #2: 

Given what you've shown, those values are only known to be identical up to the 4th decimal place. Octave simply shows you variables with limited precision (configurable using "format").

They may, of course, still yield vastly different results when used in a formula involving pseudo-inverse. Unless you supply more evidence, I assume you just didn't understand what's going on.

Jaroslav Hajek <highegg>
Sat 04 Dec 2010 04:07:53 AM UTC, comment #1: 

It would help us reproduce this if you could attach a file with the variables you used.

From within Octave, try

save -binary pinv.var wKs sigmac sigmas q Zs


and then attach the file pinv.var to this bug report.

Rik <rik5>
Group administrator
Fri 03 Dec 2010 09:04:31 PM UTC, original submission:  

In the following lines sigmas and sigmac are both singular values of the same matrix computed by different algorithms. The values these variables return in the Shell prompt are identical, i.e. for a given data set:


sigmas =
8.6534

sigmac =
8.6534


Using these values in the following lines:


kron(wKs',pinv(sigmac**2*eye(q)-Zs'*Zs)*Zs')
kron(wKs',pinv(sigmas**2*eye(q)-Zs'*Zs)*Zs')


should theoretically yield identical results. Nonetheless the first line yields values in the order of magnitude 1E-3, whereas the second yields values in the order of magnitude 1E+12.

Checked several times.
Likely due to numertical instability in the 'pinv' function.

Sven Serneels <svens>

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #22196:  Example added by None (2KiB - application/octet-stream - Data producing the error reported above.)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by gf (Posted a comment)
  • -email is unavailable- added by highegg (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by svens (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 4 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2011-01-04 rik5 Open/ClosedOpen Closed
    2010-12-13 None Attached File- Added Example, #22196
    2010-12-12 highegg StatusNeed Info Invalid / Not an Octave Bug
    2010-12-04 rik5 StatusNone Need Info

    Back to the top

    Powered by Savane 3.13-02a9.
    Corresponding source code