Re: About diagonal matrices

[John W. Eaton-3]

On 20-Feb-2009, Jaroslav Hajek wrote:

| On Fri, Feb 20, 2009 at 5:45 PM, John W. Eaton <jwe@...> 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.

For diagonal matrices, we can't fill in the zero values and keep the
diagonal matrix storage scheme.  But we could convert back to diagonal
if the off diagonal elements become zero.  Should we be doing that?
It doesn't appear to be implemented now.

jwe

[Jaroslav Hajek-2]

On Fri, Feb 20, 2009 at 6:04 PM, John W. Eaton <jwe@...> wrote:

> On 20-Feb-2009, Jaroslav Hajek wrote:
>
> | On Fri, Feb 20, 2009 at 5:45 PM, John W. Eaton <jwe@...> 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

[John W. Eaton-3]

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

[Jaroslav Hajek-2]

On Fri, Feb 20, 2009 at 6:30 PM, John W. Eaton <jwe@...> 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).

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

[Jaroslav Hajek-2]

On Fri, Feb 20, 2009 at 6:40 PM, Jaroslav Hajek <highegg@...> wrote:

> On Fri, Feb 20, 2009 at 6:30 PM, John W. Eaton <jwe@...> 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

[dbateman]

Jaroslav Hajek-2 wrote:

On Fri, Feb 20, 2009 at 6:40 PM, Jaroslav Hajek <[e-mail]> wrote:
> On Fri, Feb 20, 2009 at 6:30 PM, John W. Eaton <[e-mail]> 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.
 

[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 <jwe@...> 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

[Daniel J Sebald]

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).

A conversation about this came up at the "OctaConf" a few years back, one of the ideas being that a sparse matrix not necessarily have 0 as the "sparse default value".  That is, the "sparse default value" could be any number, 0, NaN, Inf, 2.33, etc.

I'll say more about this in a follow up to David's email...

Dan

[Daniel J Sebald]

dbateman wrote:

> 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

What does 'nnz' mean?  "something-not-zero"?  If so, then the above is wrong because any number divided by zero is not zero.  Hence, "nnz" should be 9, or 100%.

> 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.

Maybe it would make more sense to define sparse matrices in terms of index sets (or I'll offer another REALLY creative alternative below).  Say we have non-intersecting sets IJ1, IJ2, IJ3, ..., IJN which are assigned numbers V0, V1, V2, V3, ..., VN.  Say Omega is the universe, or set of all indeces for an M x N matrix.  Then the assignment for sparse matrix S would be:

(i,j) in IJ1 -> S(i,j) = V1
(i,j) in IJ2 -> S(i,j) = V2
(i,j) in IJ3 -> S(i,j) = V3
...
(i,j) in IJN -> S(i,j) = VN

and finally

(i,j) in Omega-(IJ1+IJ2+IJ3+...+IJN) -> S(i,j) = V0

That last statement means that the index in question is not in any of the other "sparse subsets".

What that would do is allow some "default sparse value", as well as perhaps give a bit more compaction in the case of several repeating

Dan


PS:  One could maybe dream up a "sub-block" scheme for sparse matrices, but that might be a programming nightmare.

[Jaroslav Hajek-2]

On Sat, Feb 21, 2009 at 1:27 AM, dbateman <dbateman@...> wrote:

>
>
>
> Jaroslav Hajek-2 wrote:
>>
>> On Fri, Feb 20, 2009 at 6:40 PM, Jaroslav Hajek <highegg@...> wrote:
>>> On Fri, Feb 20, 2009 at 6:30 PM, John W. Eaton <jwe@...> 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

[Jaroslav Hajek-2]

On Sat, Feb 21, 2009 at 1:27 AM, dbateman <dbateman@...> wrote:

>
>
>
> Jaroslav Hajek-2 wrote:
>>
>> On Fri, Feb 20, 2009 at 6:40 PM, Jaroslav Hajek <highegg@...> wrote:
>>> On Fri, Feb 20, 2009 at 6:30 PM, John W. Eaton <jwe@...> 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

[Søren Hauberg]

lør, 21 02 2009 kl. 19:48 +0100, skrev Jaroslav Hajek:
> 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).

I don't like the idea of having some variable that changes the behaviour
of programs. Assume you provide me with a function that depends the
value of this variable. Now, I have to incorporate this function in a
program of mine that depends on the variable having the opposite value
as yours. Then my combined code would have to look like this

  special_diagonal_treatment (true);
  some_function ();

  special_diagonal_treatment (false);
  some_other_function ();

Now assume we have couple of variables like this one, and what you have
is a collaboration (and debugging) nightmare.

Søren

[Jaroslav Hajek-2]

On Sat, Feb 21, 2009 at 9:08 PM, Søren Hauberg <soren@...> wrote:

> lør, 21 02 2009 kl. 19:48 +0100, skrev Jaroslav Hajek:
>> 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).
>
> I don't like the idea of having some variable that changes the behaviour
> of programs. Assume you provide me with a function that depends the
> value of this variable. Now, I have to incorporate this function in a
> program of mine that depends on the variable having the opposite value
> as yours. Then my combined code would have to look like this
>
>  special_diagonal_treatment (true);
>  some_function ();
>
>  special_diagonal_treatment (false);
>  some_other_function ();
>
> Now assume we have couple of variables like this one, and what you have
> is a collaboration (and debugging) nightmare.
>
> Søren
>
>

There are already such variables in Octave. Think about save_precision
or string_fill_char (OK, I dunno whether there are more). Sometimes it
seems the best solution. Usually it controls only corner cases (like
string_fill_char), so the risk of messing up is fairly low. The
important thing is the default; whoever sets such a global variable
should be aware of the possible risk.
Besides, I can hardly imagine a function depending on the
incompatibilities discussed.



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

[John W. Eaton-3]

On 21-Feb-2009, Jaroslav Hajek wrote:

| There are already such variables in Octave. Think about save_precision
| or string_fill_char (OK, I dunno whether there are more).

We used to have many more, but removed all the ones that could cause
obvious trouble.  The things that remain (like print_empty_dimensions,
or history_timestamp_format_string) shouldn't cause trouble except in
really strange cases.  I think save_precision falls into this category
as well.  Maybe string_fill_char should be removed.  It is probably
rarely used anyway.  I would not want to add a new option to control
the way operations on diagnoal matrices work.

jwe

[John W. Eaton-3]

On 20-Feb-2009, dbateman wrote:

| 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.

With the current sources, which I think includes your patch, I now see

  octave:2> speye (3)/0
  warning: division by zero
  ans =

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

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

  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)*NaN
  ans =

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

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


Is that the behavior you want?

jwe

[John W. Eaton-3]

On 21-Feb-2009, Jaroslav Hajek wrote:

| Mathematically, there's no such thing as a NaN.

Substitute "numerically" for "mathematically" then, where
"numerically" means according to accepted conventions of the IEEE
standard for floating point arithmetic.

| 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.

Yes, I understand that it is convenient for many uses for diagonal
and sparse matrices to have the properties you want.  But I'm also
don't like having things like

  full_matrix == diag_matrix

yet

  diag_matrix * scalar != full_matrix * scalar

for some values of scalar.

jwe

[John W. Eaton-3]

On 20-Feb-2009, Daniel J Sebald wrote:

| 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.

We try to do what is right for IEEE floating point.  Maybe sqrt should
be fixed.

There are a number of other problems with the way we (and every other
language that I know of) handle complex numbers.  For example, since
we don't have pure imaginary values, I think things like i/0 or i*Inf
fail to do the right thing.

jwe

[Søren Hauberg]

lør, 21 02 2009 kl. 16:46 -0500, skrev John W. Eaton:

> Yes, I understand that it is convenient for many uses for diagonal
> and sparse matrices to have the properties you want.  But I'm also
> don't like having things like
>
>   full_matrix == diag_matrix
>
> yet
>
>   diag_matrix * scalar != full_matrix * scalar
>
> for some values of scalar.

I don't like this either. But I must say that I really like the reduced
memory usage that comes with the diagonal matrix type (this has already
saved me a couple of times). I guess I'm missing the obvious here, but
couldn't the diagonal matrix class be extended to have a variable that
holds the value of the non-diagonal elements of the matrix? Usually this
would be zero, but when the matrix is multiplied with NaN or divided
with zero, this value would change.

Søren

[Daniel J Sebald]

Søren Hauberg wrote:

> couldn't the diagonal matrix class be extended to have a variable that
> holds the value of the non-diagonal elements of the matrix? Usually this
> would be zero, but when the matrix is multiplied with NaN or divided
> with zero, this value would change.

That is what I meant by "default sparse value", or value V0 in the set of

Omega-(IJ1+IJ2+IJ3+...+IJN) -> V0

i.e., the most prevalent value in the matrix.

Dan

[Daniel J Sebald]

John W. Eaton wrote:

> On 20-Feb-2009, Daniel J Sebald wrote:
>
> | 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.
>
> We try to do what is right for IEEE floating point.  Maybe sqrt should
> be fixed.
>
> There are a number of other problems with the way we (and every other
> language that I know of) handle complex numbers.  For example, since
> we don't have pure imaginary values, I think things like i/0 or i*Inf
> fail to do the right thing.

octave:2> i/0
warning: division by zero
ans = NaN - NaNi
octave:3> i*Inf
ans = NaN + Infi

That is a bit messy, isn't it?  As a slight fix for consistency, I'd think these two should match, i.e., i/0 = NaN + Infi.

Dan