Re: About diagonal matrices

[Daniel J Sebald]

John W. Eaton 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.

It depends on how one wants to interpret 1/0, 1/Inf I suppose.  (And perhaps 
the solution is to have some type of option specifying how to treat it.)

If I'm a user who wants to work with just rational numbers (notice I didn't say "real"), 
then 1/0 is NaN, i.e., undefined.  Also, the "number" Inf then means the CPU has reached 
the limits of its number system, i.e., 1/Inf is meaningless (or NaN) in this context as well.

If I'm a user who wants to think in terms of the extended real number system, 
then limits are of interest, i.e., Inf is the limit as series x_i becomes 
positively unbounded (but not negatively unbounded).  1/0 is fraught with 
problems though because it doesn't indicate how the limit is approached.  That 
is, how do we know that 1/0 is +Inf and not -Inf?  1/+0 could mean approach 
zero from the left, so 1/+0 = +Inf.  Likewise 1/-0 = -Inf.  Octave agrees with 
that train of thought, e.g.,

octave:20> 1/(+0)
warning: division by zero
ans = Inf
octave:21> 1/(-0)
warning: division by zero
ans = -Inf

But here is where the trouble shows up:

octave:22> 1/(+0-0)
warning: division by zero
ans = Inf
octave:23> 1/(-0+0)
warning: division by zero
ans = Inf

The problem is trying to treat 0 as both a number in the number system and the 
notion of a one-sided limit.  (Perhaps what some mathematician should have done 
long ago was to define the extended real number system as R union {+Inf, -Inf, 
+Zed, -Zed}.  Anyway...)


And how about this one, irrespective of the previous comments?

octave:42> sqrt(-0)
ans = -0

Shouldn't that be complex 0 + 0i?  That is,

octave:44> sqrt(-2)
ans =  0.00000 + 1.41421i
octave:45> sqrt(-1)
ans =  0 + 1i
octave:46> sqrt(-0.5)
ans =  0.00000 + 0.70711i
...
sqrt(lim x-> 0 from left) = 0 + 0i.


Dan