Re: About diagonal matrices

[Jaroslav Hajek]

On Sun, Feb 22, 2009 at 11:16 AM, Søren Hauberg <address@hidden> wrote:
> søn, 22 02 2009 kl. 10:14 +0100, skrev Jaroslav Hajek:
>> > 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
>> >
>>
>> Why? This behaviour has been around for ages for sparse matrices, and
>> nobody complained.
>> It's just the distinction between an assumed zero and numerical zero.
>> It's *standard* in numerical software.
>> It is usually both numerically (or practically) superior and more
>> effective, you just need to be aware of it in
>> certain cases.
>
> When the diagonal matrix type was introduced I was under the impression
> that it would only be an optimization (i.e. my programs would run faster
> but produce the same results).

It is not really so. They are data structures in their own right, and
they employ different algorithms to get results of operations. The
incompatibility you see is just a consequence of different (smarter)
algorithm being employed.
It still numerically approximates the same mathematical operation, it
just does so in a better way.
Look at sparse matrix multiplication, look at triangular system
solution (an example shown in reply to John) - they're an instance of
the same issue - you employ a different algorithm, you get slightly
different results.

Let me show you another example:

beta(Inf,1)
ans = NaN

in both Matlab and Octave.

This is "mathematically wrong" in D. Bateman's sense, because the
limit of beta(x,1) for x->Inf is zero.
This result is just an inherent consequence of the formula used for
the calculation.

Now, if you really want to eliminate all stuff like this, it means
that you virtually can't implement functions in m-files without
cluttering them by various NaN and Inf checks and suffering a
significant performance loss.
I don't regard the current Octave's beta function as useless. Fixing
even the Inf cases would be a superior implementation, but it is
ineffective to do so within an m-file implementation.
Now, however, suppose I write such an implementation and you complain
that it is inconsistent because it no longer holds that beta(x,y) =
gamma(x)*gamma(y)/gamma(x+y).
OK, this is not an extremely good analogy, but I hope you get the point.

cheers

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