

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/gtmeng.html 
Argentina  Home: Gob. Vera 3161 
Reply: [hidden email]  (54)(42) 55.00.23 


There is a function called dmult written by KH < [hidden email]>
which is available at
ftp://ftp.tsc.uvigo.es//pub/octave/contrib/mfiles/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 rowwise
kronecker product of two matrices. If you need to do this let me know and
I'll send you an mfile I wrote which does it pretty fast.
Heber Farnsworth  Department of Finance
Univerity of Washington  Box 353200
tele: (206) 5280793 home  Seattle, WA 981953200
tele: (206) 5434773 finance web: http://weber.u.washington.edu/~heberf fax: (206) 6859392 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/gtmeng.html 
> Argentina  Home: Gob. Vera 3161 
> Reply: [hidden email]  (54)(42) 55.00.23 
>


( 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 builtin code it could be done very
efficiently and fast. The decision between rowwise or columnwise
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 reweighting and large matrices (a matrix with
several thousand rows  "cases"  is not uncommon).
Best wishes to all,
Ted. ( [hidden email])


On 13Nov1996, 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


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


>>>>> 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 kindof 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/gtmeng.html 
Argentina  Home: Gob. Vera 3161 
Reply: [hidden email]  (54)(42) 55.00.23 


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


On 14Nov1996, [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

