bugGNU Octave - Bugs: bug #63080, Allow matrix power for sparse...

 
 

bug #63080: Allow matrix power for sparse empty matrices

Submitter:  Arun Giridhar <arungiridhar>
Submitted:  Mon 19 Sep 2022 11:04:30 AM UTC
   
 
Category:  Interpreter Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Unexpected Error or Warning
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Tue 20 Sep 2022 09:58:07 PM UTC, comment #5: 

I made a small change to sparse-xpow.cc to use the Octave coding convention of having a space between the function name and the opening parenthesis.  See http://hg.savannah.gnu.org/hgweb/octave/rev/a1318deb4584.

I checked with Matlab and it allows inverse of an all zero sparse matrix.  The result is odd with NaNs and Inf values so it isn't of much use.  I think it is okay to keep the error behavior.

Marking bug as fixed and closing report.

Rik <rik5>
Group administrator
Tue 20 Sep 2022 01:14:23 PM UTC, comment #4: 

Maybe something to consider with respect to that FIXME note:
Sparse matrices might be very large but sparsely populated (in your example not populated at all). If those operations would return a matrix of the same size as the input that is completely populated with Inf, it could take a lot of memory.

Markus Mützel <mmuetzel>
Group administrator
Tue 20 Sep 2022 01:03:11 PM UTC, comment #3: 

@Rik: Good points, both. I've made the changes you described with others and pushed it here: https://hg.savannah.gnu.org/hgweb/octave/rev/a1318deb4584

In the course of testing, I found an additional inconsistency between full and sparse, when the input matrix is all-zeros. For the full case, it returns an all-Inf matrix with a warning, but for the sparse case it gives an error.


>> A = zeros (2, 2); A ^ -1
warning: inverse: matrix singular to machine precision, rcond = 0
ans =

   Inf   Inf
   Inf   Inf

>> A = sparse (zeros (2, 2)); A ^ -1
error: inverse of the null matrix not defined


I've added a FIXME and marked that inconsistency as a known bug in the BISTs.

This is also the same inconsistency with calling inv() directly:

>> A = zeros (2, 2); inv (A)
warning: matrix singular to machine precision
ans =

   Inf   Inf
   Inf   Inf

>> A = sparse (zeros (2, 2)); inv (A)
error: inverse of the null matrix not defined


We could in principle detect an all-zeros matrix with a negative power and return an all-Inf matrix instead of taking the inverse, but the better fix for that inconsistency should probably be in the matrix inverse function, not in the matrix power function. What do you think?

Arun Giridhar <arungiridhar>
Group Member
Mon 19 Sep 2022 09:06:03 PM UTC, comment #2: 

I have one small request and one question.

For the request, could you add a BIST test marked with regression indicator for this bug report?  Something like


%!test <*63080>
%! A = sparse ([]);
%! assert (A ^ 1, A);
%! assert (A ^ 0, A);
%! assert (A ^ -1, A);


One question is how to handle empty matrices which are not 0x0.  For example,


A = zeros (2,0);
A ^ 1;


In Octave, this returns a 0x0 empty matrix (original dimensions are lost) while in Matlab this returns an error that the matrix must be square.

With this patch in place, the following code produces an error


A = sparse (zeros (2,0));
A ^ 1
error: for A^b, A must be a square matrix.  Use .^ for elementwise power.


This behavior is compatible with Matlab, but incompatible with Octave for full matrices.  I would argue that Octave should be consistent with itself and return the same result for full or sparse matrices.  One way to do that would be to modify the test condition to


-  if (nr == 0 && nc == 0)
+  if (nr == 0 || nc == 0)


This still won't be completely identical because for full matrices the return value always has dimensions 0x0.  We could make even that consistent by returning an empty SparseMatrix.

Rik <rik5>
Group administrator
Mon 19 Sep 2022 11:08:11 AM UTC, comment #1: 

Pushed here: https://hg.savannah.gnu.org/hgweb/octave/rev/a026fb2be108

Marking as ready for test

Arun Giridhar <arungiridhar>
Group Member
Mon 19 Sep 2022 11:04:30 AM UTC, original submission:  

Taking the matrix power of an empty sparse matrix fails in certain legal cases even though passing an empty full matrix works fine:

A = sparse ([]);
A ^ 1
A ^ 0
A ^ -1


Oddly, non-integer exponents work even for sparse empty inputs, like A ^ 1.3 or A ^ -1.3.

The problem turned out to be this input validation in sparse-xpow.cc:

   if (nr == 0 || nc == 0 || nr != nc)
     error ("for A^b, A must be a square matrix.  Use .^ for elementwise power.");


The following change fixes it for me:

diff -r 6646f2b5a3d1 libinterp/corefcn/sparse-xpow.cc
--- a/libinterp/corefcn/sparse-xpow.cc  Sat Sep 17 04:22:38 2022 -0400
+++ b/libinterp/corefcn/sparse-xpow.cc  Mon Sep 19 06:46:27 2022 -0400
@@ -66,9 +66,12 @@ xpow (const SparseMatrix& a, double b)
   octave_idx_type nr = a.rows ();
   octave_idx_type nc = a.cols ();

-  if (nr == 0 || nc == 0 || nr != nc)
+  if (nr != nc)
     error ("for A^b, A must be a square matrix.  Use .^ for elementwise power.");

+  if (nr == 0 && nc == 0)
+    return a;
+
   if (! xisint (b))
     error ("use full(a) ^ full(b)");


I will push this after I get a bug number.

Arun Giridhar <arungiridhar>
Group Member

 

(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 mmuetzel (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by arungiridhar (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
    2022-09-20 rik5 StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2022-09-19 rik5 CategoryNone Interpreter
    2022-09-19 arungiridhar StatusNone Ready For Test

    Back to the top

    Powered by Savane 3.13-d3ae.
    Corresponding source code