Re: fastest way to element-by-element multiply matix by vector?

[David Bateman]

Ozzy Lash wrote:
> On 6/19/07, Michael Creel <address@hidden> wrote:
>   
>> Hello all,
>> I need to element-by-element multiply an nXk matrix A by an nX1 vector b, so
>> that the result is that each row of A is multiplied by the corresponding row 
>> of
>> b. This is on the inside loop of some intense calculations, so I'd like the 
>> best
>> performance. Is C = A .* repmat(b,1,k) the best way to do this?
>> Thanks, Michael
>>     
>
> The only other way I can think of is:
>
> C = diag(b) * A
>
> Not sure which would be faster, replication and element by element
> multiply or diagonalization and matrix multiply.
>
>   
Better would be

C = spdiags(b,0,n,n) * A

as it doesn't allocate the additional memory of repmat and is as fast as
repmat.. You might also write that as

C = diag(sparse(b)) * A

which is about as  fast or alternatively

C = spdiag(b) * A;

though that is Octave specific and no guarantees that syntax would be
kept. Try the attached script, which runs in both Octave and Matlab..

D.

-- 
David Bateman                                address@hidden
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 6 72 01 06 33 (Mob) 
91193 Gif-Sur-Yvette FRANCE                  +33 1 69 35 77 01 (Fax) 

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary

1;
clear all;
N = [100,200,500,1000,2000];
K = [100,200,500,1000,2000];
tr = [];
td = [];
ts = [];

for n = N
  for k = K
    fprintf('A(%d,%d)\n', n, k);
    if exist ('OCTAVE_VERSION')
      eval('fflush(stdout);');
    end
    A = randn(n,k);
    b = randn(n,1);

    t0 = cputime();
    C = A .* repmat(b,1,k);
    tr(end+1) = cputime() - t0;

    t0 = cputime();
    C = diag(b) * A;
    td(end+1) = cputime() - t0;

    t0 = cputime();
    %% C = spdiags(b,0,n,n) * A;
    %% C = spdiag(b) * A;
    C = diag(sparse(b)) * A;
    ts(end+1) = cputime() - t0;
  end
end

fprintf('\n\n');
fprintf('(   n,   k)     A.*repmat(b,1,k)     diag(b)*A      spdiag(b)*A\n');
i = 1;
for n = N
  for k = K
    fprintf(  '(%4d,%4d)     %10.4f s     %10.4f s     %10.4f s\n', n, k, ...
           tr(i), td(i), ts(i));
    i = i + 1;
  end
end