Re: About diagonal matrices

[Jaroslav Hajek]

On Sat, Feb 21, 2009 at 1:27 AM, dbateman <address@hidden> wrote:
>
>
>
> Jaroslav Hajek-2 wrote:
>>
>> On Fri, Feb 20, 2009 at 6:40 PM, Jaroslav Hajek <address@hidden> wrote:
>>> On Fri, Feb 20, 2009 at 6:30 PM, John W. Eaton <address@hidden> wrote:
>>>> On 20-Feb-2009, Jaroslav Hajek wrote:
>>>>
>>>> | 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.
>>>>
>>>> I agree that these cases don't come up all that often.  So why not
>>>> give the correct answer when they do?
>>>>
>>>> jwe
>>>>
>>>
>>> That depends on how you define the correct answer.  Most numerical
>>> software uses (implicitly) the "defined zero" definition shown above
>>> for sparse and structured matrices (e.g. BLAS with banded and
>>> triangular matrices).
>>> Matlab is apparently (somewhat) inconsistent with itself in this respect.
>>> Also, in most of the cases (not the scalar OP matrix), the additional
>>> checks for NaNs would also slow down the operations (and complicate
>>> their implementation).
>>>
>>> --
>>
>> I started to write a chapter in the manual documenting the diagonal
>> and permutation matrices.
>> I will also mention the difference between assumed zeros and numeric
>> zeros, and document the existing
>> behavior. Will that be OK?
>>
>> cheers
>>
>> --
>> RNDr. Jaroslav Hajek
>> computing expert
>> Aeronautical Research and Test Institute (VZLU)
>> Prague, Czech Republic
>> url: www.highegg.matfyz.cz
>>
>>
>
>
> I consider that the fact speye(3)/0 returns a full matrix a bug unless the
> sparse_auto_mutate function is set. In fact this behavior is a hang over
> from when sparse_auto_mutate(1) was not only the default but only behavior
> and certain narrowing could be donee in the operators themselves rather than
> the narrowing function of the octave_value. I pushed a patch..
>
>
>
> octave:4> a = eye(3)/0
> a =
>
>   Inf     0     0
>     0   Inf     0
>     0     0   Inf
>
> octave:5> a = speye(3)/0
> warning: division by zero
> a =
>
> Compressed Column Sparse (rows = 3, cols = 3, nnz = 3 [33%])
> )
>
>  (1, 1) -> Inf
>  (2, 2) -> Inf
>  (3, 3) -> Inf
>
> octave:14> full(a)/0
> warning: division by zero
> ans =
>
>   Inf   NaN   NaN
>   NaN   Inf   NaN
>   NaN   NaN   Inf
>
> and would expect the NaN fill in even for diagonal matrices as John
> suggested and so the NaN values. The current Octave behavior is
> mathematically wrong and saying thats the answer the user expected is no
> excuse. We should bring back the NaN fill-in.
>
> D.
>
> --
> View this message in context: 
> http://www.nabble.com/Re%3A-About-diagonal-matrices-tp22124562p22131121.html
> Sent from the Octave - Maintainers mailing list archive at Nabble.com.
>
>

Btw. if there is an agreement that this is necessary, then a
relatively quick and painless solution (that would not further delay
3.2) is to make the special treatment of diagonal & permutation
matrices optional (that would probably involve just modifying a few
methods).
Those who do not care about the NaN incompatibilities between full &
diagonal matrices would just turn on the option
(a majority, probably).

This does not solve the "mathematically wrong" behavior of sparse
matrices (that has existed for ages), but it's better than nothing.

Would you and John be happy with such a solution? What do others think?
Which one would you rather see as a default?

regards

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