compound operators - what next?

[Jaroslav Hajek-2]

Hello,

I've added optimized handling for expressions like A'*v and v*A' where
A is a sparse matrix and v a dense matrix or vector.
(useful, e.g., for some Krylov iterations).
I've also created a simple benchmark, that is attached.

On my laptop (Pentium M, gcc 4.2.1) I get with current Octave:

constructing sparse matrix
A'*v matrix: 1.294604 s
A'*v 100 vectors: 4.204328 s
(v'*A)' 100 vectors: 1.828778 s
v*A' matrix: 0.979149 s
v*A' 100 vectors: 4.621955 s

and with my repo:

constructing sparse matrix
A'*v matrix: 0.788072 s
A'*v 100 vectors: 0.841585 s
(v'*A)' 100 vectors: 1.388352 s
v*A' matrix: 0.305787 s
v*A' 100 vectors: 1.325876 s

on an AMD Opteron (2.6GHz, Intel C++), the results are:
constructing sparse matrix
A'*v matrix: 1.032092 s
A'*v 100 vectors: 3.318786 s
(v'*A)' 100 vectors: 1.376143 s
v*A' matrix: 0.632022 s
v*A' 100 vectors: 3.684349 s

and

A'*v matrix: 0.834610 s
A'*v 100 vectors: 0.830632 s
(v'*A)' 100 vectors: 1.122641 s
v*A' matrix: 0.603371 s
v*A' 100 vectors: 1.127208 s

thus especially the matrix-vector expressions are sped up dramatically
(which is obvious, since the transpose is avoided).
It is noteworthy that A'*v is faster than (v'*A)'.


So now, the following expressions are treated in an optimized way (A,B
dense matrices, S sparse):
A'*B A.'*B A*B' A*B.' (BLAS-GEMM)
A'*A A*A' (BLAS-SYRK, HERK)
S'*A A*S' (special code)
I'd like to ask anyone interested to try compiling, or review code,
see https://tw-math.de/highegg
What is needed to get these into mainstream? Shall I do a merge with main repo?

What next?
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.

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.

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.

I'll really appreciate anyone's comments on this issue.

Also, I'll welcome suggestions for more compound operators
optimizations, like, for instance mapping A+B*I directly to complex
(A,B)
(but is that any good?)

regards,

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

[bench_sptransmul.m]

---

disp ('constructing sparse matrix')
n = 300; % size of the grid
m = 90000; % number of points
X = (n-1)*rand (m, 1); Y = (n-1)*rand (m, 1);
IX = ceil (X); JY = ceil (Y);

A = sparse(m, n^2);
A += sparse (1:m, sub2ind ([n, n], IX  , JY  ), (IX+1-X).*(JY+1-Y), m, n^2);
A += sparse (1:m, sub2ind ([n, n], IX+1, JY  ), (X - IX).*(JY+1-Y), m, n^2);
A += sparse (1:m, sub2ind ([n, n], IX  , JY+1), (IX+1-X).*(Y - JY), m, n^2);
A += sparse (1:m, sub2ind ([n, n], IX+1, JY+1), (X - IX).*(Y - JY), m, n^2);

printf ("A'*v matrix: ")
nv = 100;
v = ones (m, nv);
tic;
u = A'*v;
printf ("%f s\n", toc)

printf ("A'*v 100 vectors: ")
v = ones (m, 1);
tic;
for i=1:nv
  u = A'*v;
endfor
printf ("%f s\n", toc)

printf ("(v'*A)' 100 vectors: ")
v = ones (m, 1);
tic;
for i=1:nv
  u = (v'*A)';
endfor
printf ("%f s\n", toc)

printf ("v*A' matrix: ")
v = ones (nv, m);
tic;
u = v*A';
printf ("%f s\n", toc)

printf ("v*A' 100 vectors: ")
v = ones (1, m);
tic;
for i=1:nv
  u = v*A';
endfor
printf ("%f s\n", toc)

---

[Thomas Weber-8]

Am Montag, den 19.05.2008, 08:13 +0200 schrieb Jaroslav Hajek:
> I'd like to ask anyone interested to try compiling, or review code,
> see https://tw-math.de/highegg

Hmm, I'm surprised the above URL actually works. I couldn't even
remember having the server configured like that :)

Anyway, http://hg.tw-math.de/highegg works as well and avoids warnings
about the certificate.


        Thomas

[dbateman]

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.

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.

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

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.

Cheers
David

[Jaroslav Hajek-2]

On Mon, May 19, 2008 at 8:47 AM, Thomas Weber
<thomas.weber.mail@...> wrote:

> Am Montag, den 19.05.2008, 08:13 +0200 schrieb Jaroslav Hajek:
>> I'd like to ask anyone interested to try compiling, or review code,
>> see https://tw-math.de/highegg
>
> Hmm, I'm surprised the above URL actually works. I couldn't even
> remember having the server configured like that :)
>
> Anyway, http://hg.tw-math.de/highegg works as well and avoids warnings
> about the certificate.
>

Thanks for the hint, Thomas! Why didn't I figure that out?


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

[Jaroslav Hajek-2]

On Mon, May 19, 2008 at 9:15 AM, dbateman <dbateman@...> 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

[Thomas Weber-8]

Am Montag, den 19.05.2008, 09:58 +0200 schrieb Jaroslav Hajek:

> On Mon, May 19, 2008 at 8:47 AM, Thomas Weber
> <thomas.weber.mail@...> wrote:
> > Am Montag, den 19.05.2008, 08:13 +0200 schrieb Jaroslav Hajek:
> >> I'd like to ask anyone interested to try compiling, or review code,
> >> see https://tw-math.de/highegg
> >
> > Hmm, I'm surprised the above URL actually works. I couldn't even
> > remember having the server configured like that :)
> >
> > Anyway, http://hg.tw-math.de/highegg works as well and avoids warnings
> > about the certificate.
> >
>
> Thanks for the hint, Thomas! Why didn't I figure that out?

Currently anyone using the repositories has an ssh account. Initially, I
wanted pull over http (so for everyone) and pushing via https, that's
why I configured https in the first place.

I thought that https would be easier for people behind firewalls, but it
turned out that Mercurial has some problems with https and proxy
servers, so we are back to ssh.

        Thomas

[Søren Hauberg]

Quoting Jaroslav Hajek <highegg@...>:
> 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.

Probably a stupid question but will such a change make 'diag' and  
'eye' useless on a system that doesn't have SuiteSparse?

Søren

[Jaroslav Hajek-2]

On Mon, May 19, 2008 at 11:00 AM,  <soren@...> wrote:

> Quoting Jaroslav Hajek <highegg@...>:
>>
>> 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.
>
> Probably a stupid question but will such a change make 'diag' and 'eye'
> useless on a system that doesn't have SuiteSparse?
>
> Søren
>
>

Absolutely not. Sparse matrices with basic arithmetic work without
suitesparse, just the factorizations don't.


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

[dbateman]

Jaroslav Hajek-2 wrote:

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).

Ok for the flag coming from diag but not that diagonal matrices are probed

>
> 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?

I said diags and not diag, but I meant spdiags. I don't thinks its a great idea to have eye and diag return sparse matrices. Rather that the the acceleration you are after might already be there or easy to achieve if the user replaced diag with spdiags.

D.

[John W. Eaton]

On 19-May-2008, Jaroslav Hajek wrote:

| What is needed to get these into mainstream? Shall I do a merge with
| main repo?

I think the change is mostly OK but I've been traveling so haven't had
much time to integrate patches.  I should be working through some of
the backlog over the next few days.

jwe

[Jaroslav Hajek-2]

On Mon, May 19, 2008 at 11:09 PM, dbateman <dbateman@...> wrote:

>
>
>
> Jaroslav Hajek-2 wrote:
>>
>> 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).
>>
>
> Ok for the flag coming from diag but not that diagonal matrices are probed
>

No, probing on every multiplication, that would be a waste.

>
>
>>>
>>> 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?
>>
>
> I said diags and not diag, but I meant spdiags. I don't thinks its a great
> idea to have eye and diag return sparse matrices. Rather that the the
> acceleration you are after might already be there or easy to achieve if the
> user replaced diag with spdiags.
>

You're right, after all. I can just use spdiags and speye for my
purposes. It's no good breaking backward compatibility just to save a
couple of letters.

I'll implement special multiplication for sparse matrices marked as
diagonal and permuted diagonal.

> D.
>
>


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

[Jaroslav Hajek-2]

On Tue, May 20, 2008 at 7:31 AM, Jaroslav Hajek <highegg@...> wrote:
> On Mon, May 19, 2008 at 11:09 PM, dbateman <dbateman@...> wrote:
>
> I'll implement special multiplication for sparse matrices marked as
> diagonal and permuted diagonal.
>

Hmm - now that I think of it, I doubt that it's really worth the
effort. Row or column scaling is seldom a bottleneck (unlike general
sparse*full multiply, which often is), so I guess I can be happy with
the general algorithm and simply use spdiags instead of diag when
coding. Funny how the issue just went away...

Since I have no more good ideas for compound operators, I would claim
the work as complete.

cheers

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