Re: About diagonal matrices

[Jaroslav Hajek]

On Fri, Feb 20, 2009 at 6:04 PM, John W. Eaton <address@hidden> wrote:
> On 20-Feb-2009, Jaroslav Hajek wrote:
>
> | On Fri, Feb 20, 2009 at 5:45 PM, John W. Eaton <address@hidden> wrote:
> |
> | > But I think there are also some bugs.  Shouldn't the following return
> | > full matrices with the zero elements replaced by NaN (or -0, in the
> | > case of dividing by -Inf)?
> | >
> | >  diag ([1,2,3]) / 0
> | >  diag ([1,2,3]) / NaN
> | >  diag ([1,2,3]) / -Inf
> | >  diag ([1,2,3]) * NaN
> | >  diag ([1,2,3]) * Inf
> | >
> |
> | I think it's better in the current manner. I don't like the idea that
> | the memory can suddenly explode just because the multiplier happened
> | to be Inf. Right now, scalar * diag gives invariantly diag. This is
> | somewhat analogous to how sparse matrices behave.
>
> I see
>
>  octave:3> speye (3)*Inf
>  ans =
>
>  Compressed Column Sparse (rows = 3, cols = 3, nnz = 3 [33%])
>  )
>
>    (1, 1) -> Inf
>    (2, 2) -> Inf
>    (3, 3) -> Inf
>
>  octave:4> speye (3)/0
>  warning: division by zero
>  ans =
>
>     Inf   NaN   NaN
>     NaN   Inf   NaN
>     NaN   NaN   Inf
>
>
> So there seems to be some disagreement here, even within Octave.
>
> I'm not sure what is best here, but I don't like it that we produce
> numerically incorrect or incompatible answers just because of a
> storage scheme.
>
> BTW, for the two examples above, Matlab appears to keep the sparse
> storage scheme, but does now fill with NaN values.  If I recall
> correctly, some previous versions did not do the operations on the
> zero values at all.
>
> For sparse matrices, I don't mind switching the storage scheme, though
> I can imagine some cases where that would not be desirable either.
>

I think that what Octave does now for sparse * scalar is certainly
better than what Matlab does; I'd keep it that way. Otherwise, when
you do scalar * sparse, and just by coincidence scalar happens to be
Inf or NaN, you fill up the memory; bang, you're dead (or your
computation is).

Note that Matlab is not strictly numerically consistent either; try
"speye (3) * [NaN; 1; 1]" vs. "eye (3) * [NaN; 1; 1]".
There is simply a difference between a "numeric zero" and "assumed
zero" (or "defined zero"). The second one just always nullifies, which
is what the user actually expects.
I'd say document this somewhere, but keep the assumed zeros. IMHO
these are really corner cases and we do not need to copy the less
intelligent behavior of Matlab just for compatibility.

regards

-- 
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz