Discrepancy in sparse matrix math, when NaN's present

I expect result1 and result2 below to be identical, but they aren't. The discrepancy must be a bug, right? I'm working in R2024b, but the Run output below shows the issue exists as well in whatever the Matlab online engine is now running.

n = 5; m = 3;

S = sparse(n, m);
v=nan(n, 1); 
D=sparse(diag(v));

result1=full(v.*S)  %correct

result1 = 5×3

NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

result2=full(D*S) %incorrect

result2 = 5×3

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Matt J]

A possibly related observation. Below, I get different results depending only on whether S is in full or sparse form:

n = 5;
m = 3;
S = sparse(n, m);

v=nan(n, 1); 
D=sparse(diag(v));

result3 = full(D*full(S))

result3 = 5×3

NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

result4 = full(D*S)

result4 = 5×3

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Walter Roberson]

Are you sure that you want I and J to be empty, making S all-zero ?

randperrm(something, something_bigger) gives the empty array.

[Matt J]

Yes, I am sure I want S to be all zero. I have now simplified the code to reflect that.

[Walter Roberson]

It looks like the key is that S is sparse in D*S

n = 2; m = 1;

S = sparse(n, m);
v = zeros(n,1);
v(1) = nan;
D=sparse(diag(v))

D = 2×2 sparse double matrix (1 nonzero)
   (1,1)      NaN

result1=full(v.*S)  %correct

result1 = 2×1

NaN 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

result2a=full(D)*(S) %incorrect

result2a = 2×1

0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

result2 = full(result2a)

result2 = 2×1

0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Matt J]

Yes, it seems that when S is sparse, all NaNs in the left operand D of D*S are treated as finite, nonzeros.

n = 5;
S = sparse(n, n);

v=nan(n, 1); 
D=sparse(diag(v));
 
testOps(D,S)

ans = "D=full  S=full"

DmtimesS = 5×5

NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

ans = "D=full  S=sparse"

DmtimesS = 5×5

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

ans = "D=sparse  S=full"

DmtimesS = 5×5

NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

ans = "D=sparse  S=sparse"

DmtimesS = 5×5

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

function testOps(D,S)

  ops={@full,@sparse};

  for i=1:2

      L=ops{i};

      for j=1:2

          R=ops{j};
          
          "D="+func2str(L)+"  S="+func2str(R)

          DmtimesS=full( L(D)*R(S) ),

      end
  end


end

[dpb]

As an observation point, the same behavior exists at least as far back as R2017b; one might guess it would go back much further.

[Matt J]

@dpb Yes, I think there's little doubt that it has always been there. It's clear, from a performance optimization point of view, why you mgiht want sparse mtimes() to ignore NaNs. The whole idea of sparse type is to exploit the fact that 0*somethingFinite=0.If you have to worry about "somethingFinite" being non-finite, I wonder how much of the intended performance can be salvaged.

[dpb]

There was a similar thread on sparse performance here just a week or so ago where it ended up being full for similar reasons. I looked but didn't find it in a cursory search, iirc, it was @Steven Lord who pointed out the issue about timing.

Answer 1

[Matt J]

The response I seem to be getting from Tech Support is that this is expected behavior. NaNs are not guaranteed to propagate the same in sparse matrix multiplication as in full. In particular, sparse mtimes does not try to detect 0*NaN scenarios when the zero is sparse/implicit.

[Paul]

Is this expected behavior documented anywhere? I poked around a bit but didn't see anything; admittedly didn't try very hard (shouldn't have to).

[Matt J]

No, Tech Support seems to acknowledge that it isn't documented:

"I apologize for any confusion caused by the lack of documentation regarding NaN propagation in these operations.

To clarify, the time differences you are seeing between “.*” and “*” are due to underlying algorithms which handle NaN propagation differently. The same concept applies to the result of the second example that you provided. Specifically, the propagation rules between a sparse-sparse matrix multiply are different from a sparse-full matrix multiply.

So, although theoretically the results should be the same for both cases, the underlying code goes through a different algorithm and propagates NaN differently."