Re: diagonal & permutation matrices

[Jaroslav Hajek]

On Sun, Nov 16, 2008 at 3:53 PM, David Bateman <address@hidden> wrote:
> Jaroslav Hajek wrote:
>>
>> hi all,
>>
>> I'm working on patches allowing special treatment of diagonal and
>> permutation matrices in the Octave interpreter.
>>
>> The motivation is that the natural syntax using diagonal & permutation
>> matrices, e.g.
>> D = norm (A, 'cols');
>> A /= D;
>> [Q,R,P] = qr (A); X = P' * R \ (Q' * B);
>> is currently formidable if performance is of any concern.
>> Multiplication & division by a permutation or diagonal matrix
>> is orders of magnitude slower that the "direct" counterpart
>> (subscripting by the permutation vector or using dmult).
>> And even though Matlab does the same, I think it's a bit shameful.
>>
>> So my aim is to supply five octave_base_value subclasses:
>> octave_diag_matrix, octave_float_diag_matrix,
>> octave_complex_diag_matrix, octave_float_complex_diag_matrix,
>> octave_permutation_matrix
>> (it seems that a single permutation matrix subclass will be enough,
>> provided it will remember its numeric class)
>> and implement the necessary operations. Diagonal matrix support is
>> already present in liboctave, permutation matrix support needs to be
>> implemented (but I think it's no big deal, really).
>>
>> I've run into two issues:
>>
>> 1. What to do when printing/saving? To honor Matlab compatibility, we
>> should display & save them as full matrices. However, it is certainly
>> much more informational to print them specially, and seems a little
>> bit weird for Octave to cripple its own data when saving variables.
>>
>> 2. If D is a diagonal matrix, the even D(i,i) = x with scalar i will
>> promote D to a full matrix. There does not seem to be an easy solution
>> avoiding this. (Just like even emtpy assignment to a range promotes it
>> to a matrix)
>>
>> comments?
>>
>
>
> Isn't this why Matlab added the 'vector' option to functions like QR?

Yes, apparently. That was the "quick and lazy" solution.

> Then
> you replace the matrix multiplication with a matrix index operation. If you
> prefer to keep the matrix multiply syntax, then why not just make P a sparse
> matrix, whose matrix multiplication should be relatively fast.

Well, it will be still slower than indexing, though at least O(N^2)
and not O(N^3).
The inconvenience is that instead of the natural syntax `[l,u,p] =
lu(a)', you need to do [l,u,p] = lu(a,"vector"); p = speye(n)(p,:)`,
i.e. the default, well readable syntax is useless if you don't ignore
performance.
IMHO, functions like diag, eye, and the 3rd output args of lu and qr
should have always returned sparse matrices (because diagonal and
permutation matrices *are* quite sparse). It's another place where
Matlab compatibility is a pain in the ass.
I always disliked the fact that expressions like diag(x)*y are, in
fact, to be avoided in library functions. Using spdiag(x)*y is almost
always much better.

>
> D.
>
>
> --
> David Bateman                                address@hidden
> 35 rue Gambetta                              +33 1 46 04 02 18 (Home)
> 92100 Boulogne-Billancourt FRANCE            +33 6 72 01 06 33 (Mob)
>



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