bugGNU Octave - Bugs: bug #46103, inversion of a singular diagonal...

 
 

bug #46103: inversion of a singular diagonal matrix A produces inv(A)==A

Submitter:  jan <pfa>
Submitted:  Fri 02 Oct 2015 09:51:36 AM UTC
   
 
Category:  Libraries Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * 4.0.0 Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Wed 07 Oct 2015 02:46:26 PM UTC, comment #17: 

I changed the behavior of Octave to return Inf where the diagonal element is 0 and to preserve the return class as diagonal matrix.  See this changeset (http://hg.savannah.gnu.org/hgweb/octave/rev/f61c67865d9f).  To access the fix the original reporter will need to build from source, or wait until release 4.2.

Rik <rik5>
Group administrator
Sat 03 Oct 2015 01:46:27 PM UTC, comment #16: 

Hi,

I teach my students that the inverse of a diagonal matrix is simply a diagonal matrix with the inverse diagonal. So, I prefer solution a) in #6, I consider it the most natural and the least surprising. I agree with #5 that Octave should give the same result as with a regular matrix, but the challenge should be have Octave to recognize a regular diagonal matrix A as a diagonal matrix and compute the inverse by diag(1./diag(A)). My 2 cents...

Marco

Marco Caliari <caliari>
Group Member
Sat 03 Oct 2015 09:09:29 AM UTC, comment #15: 


A sidenote on the questions of sparsity:

blkdiag() returns sparse if one of the arguments is sparse. I believe diagonal matrices - at least here - should be considered sparse rather that full? Is there a general rule on which side diagonal matrix falls? I wonder if issparse() should return true for diagonal matrix?

jan <pfa>
Sat 03 Oct 2015 06:42:08 AM UTC, comment #14: 

Thank you everybody for quick replies. Still some comments:

Indeed, diagonal (and e.g. permutation) matrices are one of the few, but important, classes that (mathematically) preserve the very sparsity they posess under inversion. In fact diagonal matrices form a group, so all algebraic operations preserve this structure. I think this should be respected in the code. 


Now, I, too, tried this in matlab with the following finndings

1) diag() produces a full matrix, i.e. no respective diagonal class exists. 

2) Under the  inversion of the sparse matrix a sparse result is produced.

3)  the inversion of a spoarse matrix seems to be some Gaussian type of elimimination with potential divisions by zero producing inf and divisions by inf nan, without any warnings. This makes sense and I don't think this is a bad solution.


I would be more forgiving to matlab.  As a matter of principle: Producing infs and nans as a natural product of (elementwise) operations (when feasible) follows the principle of "least loss of information" and indeed makes the debugging of operations in the code easiest. So I would still keep my vote for alternative (a) in post #6.

 



8.5.0.197613 (R2015a)

>> diag(0:2)

ans =

     0     0     0
     0     1     0
     0     0     2

>> inv(sparse(diag(0:2)))
Warning: Matrix is singular to working precision.

ans =

   (1,1)          Inf
   (1,2)          NaN
   (2,2)       1.0000
   (1,3)          NaN
   (3,3)       0.5000

>> full(ans)

ans =

       Inf       NaN       NaN
         0    1.0000         0
         0         0    0.5000



jan <pfa>
Fri 02 Oct 2015 11:46:27 PM UTC, comment #13: 

WTF Matlab?  That is a weird result for the sparse case.  I don't think we should follow that.  I'll look at modifying inv.cc.

Rik <rik5>
Group administrator
Fri 02 Oct 2015 11:34:40 PM UTC, comment #12: 

Some Matlab 2015a results:


>> x = sparse (diag ([1 0]))

x =

   (1,1)        1

>> y = inv (x)
Warning: Matrix is singular to working precision.

y =

   (1,1)        1
   (2,1)      NaN
   (2,2)      Inf

>> inv (diag (0:2))
Warning: Matrix is singular to working precision.

ans =

   Inf   Inf   Inf
   Inf   Inf   Inf
   Inf   Inf   Inf

>>


Michael Godfrey <godfrey>
Group Member
Fri 02 Oct 2015 11:25:53 PM UTC, comment #11: 

I'll ask on the Maintainer's list for someone to run the Matlab code.  I knew the mutation would be hard so I think we should avoid it in dDiagMatrix.cc.  Maybe, as I said, we can figure it out in inv.cc and do the mutation there.

But, I don't think we should always promote to full.  Lots, even most, of the diagonal matrices will not be singular and so we don't want to go from linear memory consumption (N) to quadratic consumption (N^2) every time.

Rik <rik5>
Group administrator
Fri 02 Oct 2015 11:18:45 PM UTC, comment #10: 

Sorry, that's it for my Matlab hookup for now, someone else will have to test the sparse case.

I started working on a changeset for this, but changing the functions so that they return either a diagonal matrix or a full matrix proved a bit more involved than I thought.

If the inverse of a diagonal always returns a full matrix, singular or not, do we think that is also an acceptable solution?

Mike Miller <mtmiller>
Group Member
Fri 02 Oct 2015 09:56:50 PM UTC, comment #9: 

@Mike: Could you try the sparse example under Matlab?


x = sparse (diag ([1 0]))
y = inv (x)


Under Octave it simply refuses to calculate an answer


octave:6> x = sparse (diag ([1 0]))
x =

Compressed Column Sparse (rows = 2, cols = 2, nnz = 1 [25%])

  (1, 1) ->  1

octave:7> y = inv (x)
error: inv: division by zero


This could be another option to use for singular diag matrices.  It would certainly be easier to implement.  The C++ inverse function defined in dDiagMatrix.cc returns a DiagMatrix, not an octave_value.  Thus we can't easily mutate the return value from the special DiagMatrix class to the ordinary Matrix class.  I suppose one could do so at a higher level in inv.cc by checking for any element equal to 0.  That's pretty easy check on a diagonal matrix.


Rik <rik5>
Group administrator
Fri 02 Oct 2015 09:14:29 PM UTC, comment #8: 


Ok, the argument of compatibility (with matlab) is valid. Does matlab nowadays implement diagonal class? If not matlab's  behaviour might change in future.

What is behaviour of matlab for the case of the sparse matrixes.

Apart from the current behaviour, all mentioned solutions are acceptable for a mathematician (at least this one). Still the case (a) in post #6 however is the most natural one to me.

jan <pfa>
Fri 02 Oct 2015 09:10:44 PM UTC, comment #7: 

Indeed in Matlab it returns a full matrix of inf. Despite the possibility of memory overflow if the user is expecting the operation to return a diag of the same size as the input, I think that's the best fix to make here.

Mike Miller <mtmiller>
Group Member
Fri 02 Oct 2015 09:07:09 PM UTC, comment #6: 


a correction to my previous post:

My suggestion is either a) or b) below, due to the potential memory blow up.

A=diag([1,0])
a) inv(A) -> diag([1,inf])
b) inv(A) -> diag(inf,inf])
c) inv(A) -> [inf,inf;inf,inf];






jan <pfa>
Fri 02 Oct 2015 09:04:01 PM UTC, comment #5: 

For something simple like 2+2 there is a very well-defined output.  For something more complex like the inverse of a singular matrix there are several possibilities and arguments can be made for each one.  I am advocating following the Principle of Least Surprise.  The implementation of diagonal matrices is optimized in Octave, but users don't need to know or understand that.  To them, a diagonal matrix is the same as a regular matrix and the same calculation should yield the same result irrespective of what Octave has used for the underlying object type.

One useful test is to have the following code run under Matlab


inv (diag (0:2))


If this produces a 3x3 matrix of Inf values then that is what Octave should also produce.

I can't tell from comment #3 which is considered the odd behavior: 1) returning the original matrix, or 2) returning an all Inf matrix.  I certainly want to stop using method 1.  And for the next question, I would prefer to return Inf (3,3) rather than diag (Inf (3, 1))

Rik <rik5>
Group administrator
Fri 02 Oct 2015 08:59:00 PM UTC, comment #4: 

Indeed, there is already difference in the behaviour of the inverse operation for the sparse and respectively the full matrix classes. The sparse matrix only gives an error and no result as shown below.  

One aspect of using the elementwise inversion (with 1/0=inf) is that it preserves the sparsity structure and does not produce memory blow up. If one wants to follow the approach similar to the full matrix case, perhaps a diagonal matrix with the entire column inf's is a better choice:

inv(diag([1,0])) -> diag([1,inf])



octave:158> inv(sparse([1,0;0,0]))
error: inv: division by zero

octave:159> inv(diag([1,0]))
warning: inverse: matrix singular to machine precision, rcond = 1
ans =

Diagonal Matrix

   1   0
   0   0

octave:160> inv(full(diag([1,0])))
warning: inverse: matrix singular to machine precision, rcond = 0
ans =

   Inf   Inf
   Inf   Inf


jan <pfa>
Fri 02 Oct 2015 08:39:55 PM UTC, comment #3: 



With all due respect, inverse is mathematically a different concept to pseudo-inverse. And no function should knowlingly give such a misleading result that  can give very odd behaviour.

And since the diagonal matrix is now produced by e.g. diag(), already the bw-compatibility requires this.

It seems that simply the following patch would do, but as this is my first peak into octave source, I might be wrong.


jvpf@kuo:/opt/octave/latest/octave$ hg diff
diff -r 6256f6e366ac liboctave/array/fCDiagMatrix.cc
--- a/liboctave/array/fCDiagMatrix.cc        Fri Oct 02 05:50:43 2015 +0200
+++ b/liboctave/array/fCDiagMatrix.cc        Fri Oct 02 23:35:33 2015 +0300
@@ -375,10 +375,7 @@
   for (octave_idx_type i = 0; i < length (); i++)
     {
       if (elem (i, i) == 0.0f)
-        {
-          info = -1;
-          return *this;
-        }
+        retval.elem (i, i) = std::numeric_limits<double>::infinity ();
       else
         retval.elem (i, i) = 1.0f / elem (i, i);
     }
diff -r 6256f6e366ac liboctave/array/fDiagMatrix.cc
--- a/liboctave/array/fDiagMatrix.cc        Fri Oct 02 05:50:43 2015 +0200
+++ b/liboctave/array/fDiagMatrix.cc        Fri Oct 02 23:35:33 2015 +0300
@@ -280,10 +280,7 @@
   for (octave_idx_type i = 0; i < len; i++)
     {
       if (elem (i, i) == 0.0)
-        {
-          info = -1;
-          return *this;
-        }
+        retval.elem (i, i) = std::numeric_limits<double>::infinity ();
       else
         retval.elem (i, i) = 1.0 / elem (i, i);
     }


jan <pfa>
Fri 02 Oct 2015 03:49:52 PM UTC, comment #2: 

The diagonal matrix class is supposed to be a memory saving optimization.  All calculations on diagonal matrices should result in exactly the same behavior as when a full matrix is used.  My suggestion would be to change the diagonal class to return an Inf array for a singular matrix the same as is done fore a full matrix.

Rik <rik5>
Group administrator
Fri 02 Oct 2015 01:12:58 PM UTC, comment #1: 

Thanks for your bug report. This does not appear to be documented, but when inverting a singular diagonal matrix, the original matrix is returned.

As you've noticed, and as the documentation for inv() clearly states, there are often better ways to solve the problem using mldivide, or pinv.

Also, if you like this answer better,


octave:1> inv(diag(0:2))
warning: matrix singular to machine precision, rcond = 1
ans =

Diagonal Matrix

   0   0   0
   0   1   0
   0   0   2

octave:2> inv(full(diag(0:2)))
warning: matrix singular to machine precision
ans =

   Inf   Inf   Inf
   Inf   Inf   Inf
   Inf   Inf   Inf


Mike Miller <mtmiller>
Group Member
Fri 02 Oct 2015 09:51:36 AM UTC, original submission:  


Inverting (with inv) a singular diagonal matrix results the
matrix itself, but mrdivide and mldivide act correctly.

This bug affects also at least version 3.8.2. (from ubuntu 15.04. packaging.)


>> version
ans = 4.0.0
>> w=diag(0:2); id=eye(3); inv(w), w\id, id/w
warning: matrix singular to machine precision, rcond = 1
ans =
Diagonal Matrix
   0   0   0
   0   1   0
   0   0   2
ans =
Diagonal Matrix
   0.00000         0         0
         0   1.00000         0
         0         0   0.50000
ans =
Diagonal Matrix
   0.00000         0         0
         0   1.00000         0
         0         0   0.50000
>>



jan <pfa>

 

(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 caliari (Posted a comment)
  • -email is unavailable- added by godfrey (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by mtmiller (Posted a comment)
  • -email is unavailable- added by pfa (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 3 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2015-10-07 rik5 StatusNeed Info Fixed
        Open/ClosedOpen Closed
    2015-10-02 rik5 StatusNone Need Info

    Back to the top

    Powered by Savane 3.13-d3ae.
    Corresponding source code