Re: diagonal matrices specializations

[Jaroslav Hajek]

On Tue, Dec 2, 2008 at 8:23 AM, Jaroslav Hajek <address@hidden> wrote:
> hi all,
>
> this is another proposed improvement in the interpreter - diagonal
> matrix objects.
>
> Rationale:
> Diagonal matrices are common objects in linear algebra. They are
> frequently used to perform scaling and unscaling, and are also
> naturally yielded by the svd and eig factorizations. In Octave and
> Matlab alike, however, direct use of these matrices suffers a
> significant performance penalty, because they're treated as full.
> So, while expressions like `diag(D) *A' are common in manuals and
> examples, they're avoided in practice by using something ugly like `A
> .* repmat (D, 1, columns(A))'. Octave has introduced `dmult' for this
> purpose, however, this function is not much used because it is not
> Matlab compatible and it is not much known even to Octave users. And
> anyway, the whole point of having a language like Octave is to allow
> user to type equations as directly as possible, with the interpreter
> taking care of performance.
>
> Summary of changes:
> There is existing support for rectangular diagonal matrices in
> liboctave. (Yeah, rectangular. Search me why, but after a while I
> started to like the idea, it's not much additional work and it fits
> more uses nicely, SVD for instance). Some modifications were made to
> improve performance (e.g., inline element access) and extend
> functionality.
> In the interpreter code, 4 new octave_value classes are introduced (
> real/complex diagonal matrix + float counterparts), with appropriate
> operators registered and specialized in OPERATORS. Several builtins
> were modified to take advantage of the new functionality (diag, eye,
> inv) while others, due to their fine coding, benefit automatically
> (svd, eig). As a result, some tests need to be modified.
>
> Two issues are left open:
>
> 1. Any assignment into a diagonal matrix converts it to a full matrix,
> even an assigment to a single diagonal element. This is an inherent
> consequence of current Octave's type system, where you can specialize
> assignment based on types, but at the time when indices are resolved,
> the octave_base_value polymorphic instance is already fixed (see
> ops.h, octave_base_value::numeric_assign etc). I don't think this is
> much of an issue, as diagonal matrices are seldom manipulated by
> elements, but still it's a defect.
>
> 2. There was a prior debate how the special diagonal matrices should
> be displayed & saved. Saving to external formats is out of question,
> of course, but displaying and Octave's native saving is debatable.
> Displaying a diagonal matrix by printing out just the diagonal (after
> a notification of diagonality) is no doubt more informative than
> displaying it as full. Similarly, when Octave saves a matrix in its
> native format, it can be useful to preserve the diagonality
> information rather than forgetting it. OTOH, there's compatibility.
> In the patch presented, neither option is exploited, i.e. the matrix
> is always converted to full.
> My personal preference is to alter the printing and leave the saving,
> but I'd like to hear more opinions about this.
>
> In very near future, I intend to contribute a similar improvement
> concerning permutation matrices, another common type of object in
> linear algebra.
>
> The patch (>100K) is available for download here:
>
> http://artax.karlin.mff.cuni.cz/~hajej2am/ulozna/octave/diag-mat.diff
>
> cheers
>
> --
> RNDr. Jaroslav Hajek
> computing expert
> Aeronautical Research and Test Institute (VZLU)
> Prague, Czech Republic
> url: www.highegg.matfyz.cz
>

Oh yes, and as usual, here's a simple benchmark for your convenience
(set n to a suitable value):

n = 3000;
A = rand(n);
D = diag(rand(n,1));
D1 = diag(D);
# now do a double scaling
tic; C1 = dmult (dmult (D1, A), D1); toc
tic; C = D * A * D; toc
norm (C1 - C, 1)

with a recent tip, I get on my Intel Core 2 Duo @ 2.83 GHz:

Elapsed time is 0.320868 seconds.
Elapsed time is 10.1951 seconds.
ans = 0

whereas with the current patch, I get:

Elapsed time is 0.325157 seconds.
Elapsed time is 0.274062 seconds.
ans = 0

i.e. a 28x speed-up (will increase with n).
Note that division is slightly slower than multiplication, because
zeros are checked (minimum-norm solution) and because true division is
used rather than precomputing inverses due to slightly better
numerical properties (this can be revised, but I dont' think it's an
issue).

cheers

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