Efficient multiplication by a diagonal matrix

classic Classic list List threaded Threaded
9 messages Options
Reply | Threaded
Open this post in threaded view
|

Efficient multiplication by a diagonal matrix

Mario Storti-3

I found  myself repeatedly with the following  problem. Given a matrix
A(n,m)  and a vector  v(n), I  have to  multiply   each row  A(j,:) by
v(j). This is equivalent to compute:

B = diag(v) * A                     (1)

Now, for large    n, (1) is  very  inefficient,   because it  requires
constructing the square matrix diag(v) which requires storage and many
redundant operations since most elements  of diag(v) are null. If n>>m
then:

B= kron(v,ones(1,m)).*A             (2)

does the job  and is better.  But the more  efficient way is computing
row by row if m>>n and column by column  if n>>m. However, I repeat, I
find this problem so many times and in so many  areas that it seems to
me that some system call should do it.

I wrote some code  of my own  to  do this  task,  but I wonder  if I'm
redeveloping the wheel. Does anyone have a betetr solution?

Mario

%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%
Mario Alberto Storti               | Fax: (54)(42) 55.09.44               |
Grupo de Tecnologia Mecanica       | Tel: (54)(42) 55.91.75               |
INTEC, Guemes 3450 - 3000 Santa Fe | http://venus.unl.edu.ar/gtm-eng.html |
Argentina                          | Home: Gob. Vera 3161                 |
Reply: [hidden email]  |       (54)(42) 55.00.23              |

Reply | Threaded
Open this post in threaded view
|

Re: Efficient multiplication by a diagonal matrix

Heber Farnsworth

There is a function called dmult written by KH <[hidden email]>
which is available at

ftp://ftp.tsc.uvigo.es//pub/octave/contrib/m-files/lalgebra/dmult.m

I use it all the time.  In my work I often have the need to do something
which is little more complicated  but related.  That is a row-wise
kronecker product of two matrices.  If you need to do this let me know and
I'll send you an m-file I wrote which does it pretty fast.

  Heber Farnsworth                               | Department of Finance
  Univerity of Washington                        | Box 353200
  tele:  (206) 528-0793 home                     | Seattle, WA 98195-3200
  tele:  (206) 543-4773 finance     web: http://weber.u.washington.edu/~heberf
  fax:   (206) 685-9392             email: [hidden email]

On Tue, 12 Nov 1996, Mario Storti wrote:

>
> I found  myself repeatedly with the following  problem. Given a matrix
> A(n,m)  and a vector  v(n), I  have to  multiply   each row  A(j,:) by
> v(j). This is equivalent to compute:
>
> B = diag(v) * A                     (1)
>
> Now, for large    n, (1) is  very  inefficient,   because it  requires
> constructing the square matrix diag(v) which requires storage and many
> redundant operations since most elements  of diag(v) are null. If n>>m
> then:
>
> B= kron(v,ones(1,m)).*A             (2)
>
> does the job  and is better.  But the more  efficient way is computing
> row by row if m>>n and column by column  if n>>m. However, I repeat, I
> find this problem so many times and in so many  areas that it seems to
> me that some system call should do it.
>
> I wrote some code  of my own  to  do this  task,  but I wonder  if I'm
> redeveloping the wheel. Does anyone have a betetr solution?
>
> Mario
>
> %%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%
> Mario Alberto Storti               | Fax: (54)(42) 55.09.44               |
> Grupo de Tecnologia Mecanica       | Tel: (54)(42) 55.91.75               |
> INTEC, Guemes 3450 - 3000 Santa Fe | http://venus.unl.edu.ar/gtm-eng.html |
> Argentina                          | Home: Gob. Vera 3161                 |
> Reply: [hidden email]  |       (54)(42) 55.00.23              |
>


Reply | Threaded
Open this post in threaded view
|

Re: Efficient multiplication by a diagonal matrix

Ted.Harding
In reply to this post by Mario Storti-3
( Re Message From: Mario Storti )

>
>
> I found  myself repeatedly with the following  problem. Given a matrix
> A(n,m)  and a vector  v(n), I  have to  multiply   each row  A(j,:) by
> v(j). This is equivalent to compute:
>
> B = diag(v) * A                     (1)
>
> Now, for large    n, (1) is  very  inefficient,   because it  requires
> constructing the square matrix diag(v) which requires storage and many
> redundant operations since most elements  of diag(v) are null.

A good while ago I wrote to John Eaton suggesting an extension of the
octave ".*" multiplication operator so that, if A (m x n) is a matrix,
u (m x 1) and v (1 x n) are vectors, then the following could be written:

        u .* A  =  A .* u
        =  rows of A times corresponding elements of u (the above case)

        v .* A  =  A .* v
        =  columns of A times corresponding elements of v

If this were implemented by built-in code it could be done very
efficiently and fast. The decision between row-wise or column-wise
multiplication would be taken according to the dimensions of u (or v).
Any other dimensional relationships (except u or v scalar) would be an
error.

John reacted favourably at the time, but I have heard no more since.

I certainly have frequent call for such operations: in Statistics,
applying the same weighting factor along a whole row (or column) is a
daily requirement; I have tended to adopt the diag(u)*A method, but as
Mario points out this can be inefficient, especially in computations
involving iterative re-weighting and large matrices (a matrix with
several thousand rows -- "cases" -- is not uncommon).

Best wishes to all,
Ted.                                    ([hidden email])

Reply | Threaded
Open this post in threaded view
|

Re: Efficient multiplication by a diagonal matrix

John W. Eaton-6
On 13-Nov-1996, Ted Harding <[hidden email]> wrote:

: A good while ago I wrote to John Eaton suggesting an extension of the
: octave ".*" multiplication operator so that, if A (m x n) is a matrix,
: u (m x 1) and v (1 x n) are vectors, then the following could be written:

  [...]

: John reacted favourably at the time, but I have heard no more since.

I do want to add this feature.  Someone has even offered patches.  The
problem is mostly lack of time on my part, though I expect that
eventually Octave will have this feature.

: I certainly have frequent call for such operations: in Statistics,

I also found just yesterday that I could have used this feature, but
for a different operator.  What I wanted to do was subtract the mean
of each column of a matrix from the elements of the corresponding
column.  What I ended up doing was

  x - ones (rows (x), :) * mean (x)

but it would have been much simpler to be able to just write

  x .- mean (x)

jwe

Reply | Threaded
Open this post in threaded view
|

Re: Efficient multiplication by a diagonal matrix

SANDS-3
In reply to this post by Mario Storti-3

Mario Storti writes:
>I found  myself repeatedly with the following  problem. Given a matrix
>A(n,m)  and a vector  v(n), I  have to  multiply   each row  A(j,:) by
>v(j). This is equivalent to compute:
>
>B = diag(v) * A                     (1)
>
>Now, for large    n, (1) is  very  inefficient ...

Right, try using "Tony's Trick", it's been around for years.

Assuming v is a column vector (if not do v=v(:);) try

B=v(:,ones(length(v),1)) .* A;

This is a kind-of MATLAB specific thing but I think
it works in octave o.k.

Scott Sands

Reply | Threaded
Open this post in threaded view
|

Re: Efficient multiplication by a diagonal matrix

SANDS-3
In reply to this post by Mario Storti-3

Sorry about the previous post, the right answer is

B=v(:,ones(1,m)) .*A;

ss

Reply | Threaded
Open this post in threaded view
|

Re: Efficient multiplication by a diagonal matrix

Mario Storti-3

>>>>> On Thu, 14 Nov 1996 17:43:56 -0400 (EDT),
                       [hidden email] said:

> Mario Storti writes:
> >I found  myself repeatedly with the following  problem. Given a matrix
> >A(n,m)  and a vector  v(n), I  have to  multiply   each row  A(j,:) by
> >v(j). This is equivalent to compute:
> >
> >B = diag(v) * A                     (1)
> >
> >Now, for large    n, (1) is  very  inefficient ...
>
> Right, try using "Tony's Trick", it's been around for years.
>
> Assuming v is a column vector (if not do v=v(:);) try
>
> B=v(:,ones(length(v),1)) .* A;
>
> This is a kind-of MATLAB specific thing but I think
> it works in octave o.k.
>
> Scott Sands

(...)

> Sorry about the previous post, the right answer is
>
> B=v(:,ones(1,m)) .*A;
>
> ss


I tried, and it works also in Octave. Curiously enough, both work:

>>    octave:3> v=rand(5,1)
>>    v =
>>    
>>      0.397911
>>      0.517248
>>      0.876939
>>      0.025083
>>      0.895795
>>    
>>    octave:4> v(:,ones(5,1))
>>    ans =
>>    
>>      0.397911  0.397911  0.397911  0.397911  0.397911
>>      0.517248  0.517248  0.517248  0.517248  0.517248
>>      0.876939  0.876939  0.876939  0.876939  0.876939
>>      0.025083  0.025083  0.025083  0.025083  0.025083
>>      0.895795  0.895795  0.895795  0.895795  0.895795
>>    
>>    octave:5> v(:,ones(1,5))
>>    ans =
>>    
>>      0.397911  0.397911  0.397911  0.397911  0.397911
>>      0.517248  0.517248  0.517248  0.517248  0.517248
>>      0.876939  0.876939  0.876939  0.876939  0.876939
>>      0.025083  0.025083  0.025083  0.025083  0.025083
>>      0.895795  0.895795  0.895795  0.895795  0.895795
>>    

But it is equivalent to kron(v,ones(1,m)) and it has the same problems
of inefficiency as mentioned in the original post:

> Now, for large    n, (1) is  very  inefficient,   because it  requires
> constructing the square matrix diag(v) which requires storage and many
> redundant operations since most elements  of diag(v) are null. If n>>m
                                                                 ^^^^^^^
> then:
>
> B= kron(v,ones(1,m)).*A             (2)
>
> does the job  and is better.  But the more  efficient way is computing
  ^^^^^^^^^^^^^^^^^^^^^^^^^^^
> row by row if m>>n and column by column  if n>>m. However, I repeat, I
> find this problem so many times and in so many  areas that it seems to
(...)

Mario

%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%
Mario Alberto Storti               | Fax: (54)(42) 55.09.44               |
Grupo de Tecnologia Mecanica       | Tel: (54)(42) 55.91.75               |
INTEC, Guemes 3450 - 3000 Santa Fe | http://venus.unl.edu.ar/gtm-eng.html |
Argentina                          | Home: Gob. Vera 3161                 |
Reply: [hidden email]  |       (54)(42) 55.00.23              |

Reply | Threaded
Open this post in threaded view
|

Re: Efficient multiplication by a diagonal matrix

SANDS-3
In reply to this post by Mario Storti-3
Mario Storti writes:
[stuff deleted]
>But it is equivalent to kron(v,ones(1,m)) and it has the same problems
>of inefficiency as mentioned in the original post:

Not so! In
B=kron(v,ones(1,m)).*A;
the first factor (kron(v,ones(1,m))) is calculated
by actually carrying out a number of scalar multiplications
in the computer (if you look at the source for kron
you'll see that you have the additional overhead of
performing loops in the high level MATLABish language--
very costly!) However by using "Tony's Trick" the first
factor in
B=v(:,ones(1,m)) .* A;
is determined by simply copying the elements of v into
memory--a potentially *much* faster operation (depending
on your hardware). I'm not sure about how Tony's Trick
is carried out in octave but I'm fairly certain that
the MATLAB interpreter is smart about how it executes
such indexing operations. I've never performed actual
tests to compare the two (because I never thought of
using kron to accomplish this operation) so I can't
say with absolute certainty which is faster.

Scott Sands

Reply | Threaded
Open this post in threaded view
|

Re: Efficient multiplication by a diagonal matrix

John W. Eaton-6
On 14-Nov-1996, [hidden email] <[hidden email]> wrote:

: B=kron(v,ones(1,m)).*A;

Yes, kron will definitely be slow because it uses interpreted loops.

: B=v(:,ones(1,m)) .* A;

This will be fairly fast.  In the test I ran, with v = rand (1000,1)
and m = 10, it was nearly 600 times faster than kron.  For some cases
though,

  v * ones(1,m) .* A

will actually be faster, at least in Octave.

If the .* operator is overloaded to do the job of row and column
scaling, then I would expect that

  v .* A

will be faster than any of the other methods.  It will also use less
memory, which might turn out to be important if length(v) columns(A)
are both large.  It might even be easier to read, too.

So, why am I writing this instead of implementing that?  :-)

jwe