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.

Mathematically, there's no such thing as a NaN.
If by "mathematically wrong" you mean that "matrices don't behave
identically to their dense equivalents", then that is true.

Consider:
Inf * speye (3):
ans =

Compressed Column Sparse (rows = 3, cols = 3, nnz = 3 [33%])
)

  (1, 1) -> Inf
  (2, 2) -> Inf
  (3, 3) -> Inf

By your definition, this is "mathematically wrong". Now, if you'd want
to make this being filled with NaNs, then that should go along with a
change that allows "default sparse value" other than zero. Otherwise,
one will be forced to check for NaN prior to any scalar * sparse
multiplication to avoid possibly exploding memory.
This is the situation in Matlab 2007a - it's completely screwed up, because
a = sparse (100000); b = Inf; c = b * a;
makes Matlab effectively hang.

Now, in principle, that is still relatively easily possible to
implement. Consider another "mathematically wrong" behavior:
speye (3) * [NaN; 1; 1]
ans =

   NaN
     1
     1

eye (3) * [NaN; 1; 1]
ans =

   NaN
     1
     1
(thanks to diagonal matrices)

full (eye (3)) * [NaN; 1; 1]
ans =

   NaN
   NaN
   NaN

i.e. the "mathematically correct".
What does it take to implement the non-"mathematically wrong" behavior
for diagonal * vector operation?
Prior to multiplication, you need to check the vector for NaNs and
Infinities. If there is a NaN, the result will be all NaNs.
If there's an Infinity, it depends on their count. If there are at
least two, the result will be all NaNs. If there is a single infinity,
the result will have NaNs in all elements except the one corresponding
to the Infinity.
For left multiplication, this will slow down the operation, say, by a
factor of 2.
For right, it will be probably worse, due to cache problems (accessing rows).

Now, that is just for diagonal matrices. For general sparse matrices,
you solve the same problems, but it gets much more complicated.
For permutation matrices, it's analogical - it no longer transforms to
a simple indexing operation.

In my opinion, the current behaviour is the most useful and sane and
is OK if it's documented.
If you intend to work on "fixing" these "issues", the please make the
"mathematically wrong" behavior optional by a configuration variable,
because I will most definitely want it to be set.

regards


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