bugGNU Octave - Bugs: bug #58215, sparse multiplication with Inf

 
 

bug #58215: sparse multiplication with Inf

Submitted by:  Guillaume <gyom>
Submitted on:  Mon 20 Apr 2020 03:40:47 PM UTC  
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Matlab Compatibility
Status:  Confirmed Assigned to:  None
Originator Name:  Guillaume Open/Closed:  Open
Release:  dev Operating System:  Any

Add a New Comment (Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

( Jump to the original submission)

Fri 24 Apr 2020 09:12:56 PM UTC, comment #7: 

Matlab special case with scalar operands is indeed a bit odd - I wonder if that's an historical decision they have to live with for compatibility reasons. It is also a situation where it is 'easy' to implement the IEEE rules as one only need to check whether said operand is finite or not before branching to the special case.

I was mentioning diagonal matrices because it creates a situation where it might be difficult to write code that will behave similarly in both Matlab and Octave; given that something as ubiquitous as eye(X) follows the 'full' rules in Matlab and the 'sparse' rules in Octave.

Finally, I was curious to find out how the decisions were made for the Julia language and there are some interesting threads, e.g.:
https://github.com/JuliaLang/julia/issues/22733

Guillaume <gyom>
Fri 24 Apr 2020 04:16:22 PM UTC, comment #6: 

I agree that the fact that the end result depends on the representation is not a good result.  The user shouldn't have to understand that a full, diagonal, or sparse matrix representation is being used when they contemplate multiplying two variables.

For the moment we can set aside diagonal matrices which Matlab doesn't implement (it just uses full matrices even for diag matrices).

Here is a table of possibilities showing the various outcomes:

full scalar * full matrix   : IEEE rules
full matrix * full matrix   : IEEE rules
full scalar * sparse matrix : IEEE rules (Matlab only)
sparse scalar * sparse matrix : IEEE rules (Matlab only)
sparse matrix * sparse matrix : ! IEEE rules

I'm not a big user of sparse matrices so I don't understand why you would want to have an exception if one of the arguments is a scalar.  I can understand why you would want to have different rules for sparse matrix by sparse matrix multiplication because it does simplify the algorithm.  I guess what I would want to know is how many practitioners who regularly use sparse matrices would find strict implementation of IEEE rules useful.  My guess is that most algorithms have been developed assuming the current behavior and so it wouldn't be particularly useful.

Rik <rik5>
Project Administrator
Fri 24 Apr 2020 01:22:00 PM UTC, comment #5: 

Thanks for the feedback and extra readings. So, in Matlab, one departure of the behavior you describe is when one of the operands is a scalar:

Inf * sparse (2, 2, 0)
Inf (2, 2) * sparse (2, 2, 0)

where Matlab returns this:

>> Inf * sparse (2, 2, 0)
ans =
   (1,1)      NaN
   (2,1)      NaN
   (1,2)      NaN
   (2,2)      NaN
>> Inf (2, 2) * sparse (2, 2, 0)
ans =
     0     0
     0     0

I am starting to feel a little uneasy when adding to this the Diagonal matrices that Octave has but not Matlab:

eye (2) * [Inf; 1]
speye (2) * [Inf; 1]
[1, 0; 0, 1] * [Inf; 1]
disable_diagonal_matrix (true)
eye (2) * [Inf; 1]

and if sparse matrices can automatically mutate to full matrices with sparse_auto_mutate then it becomes very difficult to know which operation will take place.

Guillaume <gyom>
Mon 20 Apr 2020 05:36:01 PM UTC, comment #4: 

See bug #36562 from 2012 which involves sparse storage interacting with NaN.  In general, the IEEE special values (-Inf, +Inf, NaN) behave differently in Octave depending on whether the matrix is full or sparse.

I wrote up some comments there a long time ago.  This is fixable, and probably should be, but it will significantly change the way sparse matrix multiplication gets handled.  Currently, multiplication only takes place between the set of matrix elements which are both non-zero, because the code assumes that 0 multiplied by any other value is still zero.

0 * x == 0 for all x

This is a valid statement in pure math, but not in the IEEE system used to represent numbers.  In this case, the following holds

0 * x == 0 for all normal numerical values of x
0 * Inf == NaN
0 * -Inf == NaN
0 * NaN == NaN

To see what I mean

octave:1> x = sparse ([1, 0]);
octave:2> y = sparse ([0, NaN]);
octave:3> x + y
ans =

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

  (1, 1) -> 1
  (1, 2) -> NaN

octave:4> x .* y
ans =

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


octave:5> z = sparse ([0, Inf]);
octave:6> x .* z
ans =

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

Roughly, the operation in Octave is

z = x .* y =>
idx1 = find (x);   # Find non-zeros in x
idx2 = find (y);   # Find non-zeros in y
idx3 = intersect (idx1, idx2);
z = zeros (size (x));
z(idx) = x(idx) .* y(idx);

Rik <rik5>
Project Administrator
Mon 20 Apr 2020 04:11:27 PM UTC, comment #3: 

I don't know that it is noise.  I also came across this difference recently.

John W. Eaton <jwe>
Project Administrator
Mon 20 Apr 2020 04:08:21 PM UTC, comment #2: 

And, at least, Julia matches Octave.

julia> Inf * 0
NaN

julia> Inf * sparse([0])
1-element SparseVector{Float64,Int64} with 0 stored entries

julia> Inf * sparse([1;0])
2-element SparseVector{Float64,Int64} with 1 stored entry:
  [1]  =  Inf

Sorry about the noise, you can close this report. I will probably open a similar one every few years...

Guillaume <gyom>
Mon 20 Apr 2020 04:02:07 PM UTC, comment #1: 

I should have found and read this before:
https://octave.org/doc/v5.2.0/Zeros-Treatment.html

Guillaume <gyom>
Mon 20 Apr 2020 03:40:47 PM UTC, original submission:  

After being surprised that output would differ because of matrix storage (but where Matlab and Octave agree), I found a difference with Inf*sparse(0): Matlab says sparse(NaN) while Octave says sparse(0).

octave:1> diag ([1, Inf]) * sparse ([1; 0])
ans =
Compressed Column Sparse (rows = 2, cols = 1, nnz = 1 [50%])
  (1, 1) -> 1

octave:2> diag ([1, Inf]) * full   ([1; 0])
ans =
     1
   NaN

octave:3> Inf * 0
ans = NaN

octave:4> Inf * sparse (0)
ans = Compressed Column Sparse (rows = 1, cols = 1, nnz = 0 [0%])

octave:5> Inf * sparse ([1; 0])
ans =
Compressed Column Sparse (rows = 2, cols = 1, nnz = 1 [50%])
  (1, 1) -> Inf

In Matlab:

>> diag ([1, Inf]) * sparse ([1; 0])
ans =
     1
     0

>> diag ([1, Inf]) * full   ([1; 0])
ans =
     1
   NaN

>> Inf * 0
ans =
   NaN

>> Inf * sparse (0)
ans =
   (1,1)      NaN

>> Inf * sparse ([1; 0])
ans =
   (1,1)      Inf
   (2,1)      NaN

Another difference highlighted with the first test is that full x sparse is sparse in Octave and full in Matlab.

Guillaume <gyom>

 

(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 jwe (Posted a comment)
  • -email is unavailable- added by gyom (Submitted the item)
  •  

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

     

     

    Follows 1 latest change.

    Date Changed by Updated Field Previous Value => Replaced by
    2020-04-20 rik5 StatusNone => Confirmed

    Back to the top


    Powered by Savane 3.7