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
>

OK, I've pushed the patches.

Rationale: as in the first mail.

Summary of changes:

Rectangular diagonal matrices are now special objects. They can be
constructed using the "diag" builtin function, as well as returned by
other builtin functions (svd, eig). Row/column scaling & unscaling of
matrices can easily be done using the *, /, \ operators. The division
operators return a minimum-norm solution.
If D is a diagonal matrix, D(:,:) can be used to convert it to a full matrix.

Permutation matrices are also special objects. They can be
conveniently constructed using the syntax
`eye(n)(p,:)', `eye(n)(:,p)', `eye(n)(p,q)' as well as returned by
some builtins (lu, qr). Permuting / unpermuting is easily done using
*, /, \ or '*, *'.
If P is a permutation matrix, P(:,:) can be used to convert it to a full matrix.

The *,/,\ operations with diagonal & permutation matrices have changed
their complexity from O(N^3) to O(N^2), which significantly improves
the speed and allows to use them normally in matrix algebra without
performance penalty. A simplistic benchmark follows (set n to a
suitable value):

n = 3000;
A = rand(n);
D = diag(rand(n,1));
D1 = diag(D);
tic; C1 = dmult (dmult (D1, A), D1); toc
tic; C = D * A * D; toc
norm (C1 - C, 1)

p = randperm (n);
P = eye(n)(p,:);

tic; C1 = A(p,p); toc
tic; C = P * A * P'; toc
norm (C1 - C, 1)


before the changes, I get (Intel Core 2 Duo @ 2.8 GHz)

Elapsed time is 0.328404 seconds.
Elapsed time is 10.1832 seconds.
ans = 0
Elapsed time is 0.06 seconds.
Elapsed time is 1e+01 seconds.
ans = 0

with my changes, I get:

Elapsed time is 0.329819 seconds.
Elapsed time is 0.267672 seconds.
ans = 0
Elapsed time is 0.07 seconds.
Elapsed time is 0.1 seconds.
ans = 0


i.e. a 38x, resp. 100x speed-up (will increase with n).

TODO:
1. special printing/saving of diagonal & permutation matrices
2. make "balance" benefit from the changes
3. I'm sure there's a lot of stuff I forgot.
4. Sparse code is mostly untouched. When dealing with diagonal &
permutation matrices as sparse, the performance hit is much smaller
(not by an order, but a constant factor). As I don't work with sparse
matrices much, I'd appreciate advice about what interactions would be
beneficial to implement (e.g., diagonal-sparse multiplication etc).

I'll attend to 1,2,3 shortly.

cheers

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