Re: compound operators - what next?

[Jaroslav Hajek]

On Mon, May 19, 2008 at 9:15 AM, dbateman <address@hidden> wrote:
>
>
>
> Jaroslav Hajek-2 wrote:
>>
>> My top candidate is optimizing expressions like diag(v)*A because row
>> and column scaling tend to be quite common in my work (what about
>> yours?)
>> Since this will be more like a syntactic sugar (the above
>> optimizations really use different algorithms, or at least differently
>> optimized), I do not consider it that important; however, diag(v)*A is
>> IMHO more readable than dmult (v, A) (btw. what you do in Matlab?)
>> I see 3 ways open here:
>>
>> 1. implement compound operators diag_mul, mul_diag, diag_ldiv,
>> div_diag, so that the following expressions would get optimized:
>> diag(v)*A, A*diag(v), diag(v)\A, A/diag(v). Perhaps some checking will
>> be needed if `diag' refers to the intrinsic function.
>>
>
> But then if the user wrote d=diag(v); d*A instead your optimization would be
> lost.
>

Precisely - that is what compound operators are about. Do you use
expressions like the above (i.e. store a diagonal matrix for later
scaling)? Or better; would you like to use it?

>
>> 2. this is suggested in Octave's old TODO: make v.*A work like
>> diag(v)*A. I don't like it much; but it's an option.
>>
>
> Sometime around 2.1.43 Paul Kienzle implemented this. The patch for it
> stayed in octave-forge for years and as NDArrays were added afterwards
> became unusable. I suspect John doesn't really like the idea either.
>

OK, that solves the matter, I think.

>
>
>> 3. make `diag' and `eye' actually output special objects (i.e.
>> specialize octave_value for holding DiagMatrix). This would be the
>> most complex, since it would allow doing things like:
>> D = diag (v); % create scaling matrix
>> S1 = D*S;  % scale
>> S2 = S1 / (D+lambda*eye (n))
>> It would also partially eliminate the need for spdiag and speye,
>> because things like S+eye(n) or sparse(eye(n)) would work efficiently.
>> The question here is sensible diagonal matrix arithmetic vs. matlab
>> compatibility. For instance, should exp(D) for D diagonal only
>> exponentiate the diagonal? (this is sensible from a practical point of
>> view, but breaks Matlab compatibility).
>> Another question is whether this syntactic sugar is worth creating
>> another two octave_value classes (and thus increasing code bloat),
>> especially when most of the code will have to ensure compatibility.
>>
>
> Is the above just a diagonal sparse matrix stored in a different form. In
> that case what does diags(v)*A give you? Improvements I'd suggest are
>
> * Make diags flag the matrix type as diagonal if it is
> * For the spare-dense multiply functions special case the diagonal scaling
> case
>

Yes! I've considered this, too. If you intend diag and eye to return a
sparse matrix, then that would break backwards & Matlab compatibility.
I would like it, however, because I bet that were sparse matrices in
Matlab from the beginning, these functions would return sparse
matrices. They *are* sparse, after all.

Other possibility is to flag even dense matrix with MatrixType
"diagonal" and make diag and eye return such a matrix; in such case,
diag(v)*A would be time-efficient but not memory-efficient (because
the matrix would still get constructed).


>
> This would pretty much implement the third option with little work and only
> require the user to use sparse diagonal matrices instead for the scaling.
> Its also not using tricks in the parser to get the wanted behavior and so
> can't be fooled by d=diag(v);d*A as your other technique will be.
>

exactly, exactly ...
I hope I did not misunderstood you; it seems to me that you are OK
with diag and eye
returning sparse matrices, right?

What do others think?
I think there is good chance that Mathworks will do this anyway some
day, because (I think) it makes more sense.


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