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

 
 

bug #58215: sparse multiplication with Inf

Submitter:  Guillaume <gyom>
Submitted:  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
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

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

    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.

     

    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.13-f8d8.
    Corresponding source code