Thread Subject: Bogey

Hi,

I'm trying to avoid using loops when multiplying matrices together iteratively and was wondering if anyone has a trick to do this quickly and without loops. My algebra skills are somewhat lacking, so hopefully this is straightforward!

X = [1 2; 5 4];
Y = [1 2 1 ;2 4 1];
ksi = [1 2 ;2 3];

Here's the equation that I'm attempting to solve:
    \sum_{n=1}^N\sum_{i=1}^M \xi_{ni}X_nY_i
where a Xn would be a column (e.g. [1;5]) and a Yi would be a row (e.g. [1 2 1]). We should end up with a 2x3 matrix.

Here's the code I'm currently using:

tot = zeros(2,3);
for n = 1:2,
  for i = 1:2,
    tot = tot+ X(:,n)*Y(i,:);
  end;
end;

Does anyone have a suggestion how to do this without these loops? These matrices are quite big in the actual program, and this section takes a long time to run.

Thanks!

Bruce

I'm not sure it will help performance, but you could get rid of one of the loops with repmat.

for i = 1:2,
     tot = tot + X*repmat(Y(i,:),2,1);
end

Algebraically there are simple expressions for sums like these involving Kronecker products with identity matrices.

These perform badly in code since they involve multiplying very sparse matrices (and allocating large matrices, or using splapack rather than blas).

"b b" <bogey4@gmail.com> wrote in message <gdvjlp$kij$1@fred.mathworks.com>...
> Hi,
>
> I'm trying to avoid using loops when multiplying matrices together iteratively and was wondering if anyone has a trick to do this quickly and without loops. My algebra skills are somewhat lacking, so hopefully this is straightforward!
>
> X = [1 2; 5 4];
> Y = [1 2 1 ;2 4 1];
> ksi = [1 2 ;2 3];
>
> Here's the equation that I'm attempting to solve:
> \sum_{n=1}^N\sum_{i=1}^M \xi_{ni}X_nY_i
> where a Xn would be a column (e.g. [1;5]) and a Yi would be a row (e.g. [1 2 1]). We should end up with a 2x3 matrix.
>
> Here's the code I'm currently using:
>
> tot = zeros(2,3);
> for n = 1:2,
> for i = 1:2,
> tot = tot+ X(:,n)*Y(i,:);
> end;
> end;
>
> Does anyone have a suggestion how to do this without these loops? These matrices are quite big in the actual program, and this section takes a long time to run.
>
> Thanks!
>
> Bruce
------
 tot = sum(X,2)*sum(Y,1);

Roger Stafford

Thanks Roger, that's exactly what I was looking for! Unfortunately, I missed one variable in the code I posted. I meant to write:

tot = 0;
for n = 1:2,
  for i = 1:2,
    tot = tot+ ksi(n,i)*X(:,n)*Y(i,:);
  end;
end;

Where
ksi = [1 2 ;2 3];

Any suggestions?

Bruce


> ------
> tot = sum(X,2)*sum(Y,1);
>
> Roger Stafford

"b b" <bogey4@gmail.com> wrote in message <gdvq28$482$1@fred.mathworks.com>...
> .... I meant to write:
> ....
> for n = 1:2,
> for i = 1:2,
> tot = tot+ ksi(n,i)*X(:,n)*Y(i,:);
> end;
> end;

  In that case try this:

tot = X*ksi*Y;

It should give the correct answer. Whether it is faster you will have to determine for yourself. You can experiment with parentheses: (X*ksi)*Y as compared with X*(ksi*Y).

Roger Stafford

Brilliant! The first line ( tot = X*ksi*Y; ) gave me the correct answer. I'm glad to now know an elegant and straight forward solution to this. As for computational speed, it seems that the equation that you suggested is faster:

tot = zeros(2,3);
tic;
for n = 1: 4;
for i = 1: 4;
tot = tot+ ksi(n,i)*X(:,n)*Y(i,:);
end;
end;
toc
Elapsed time is 0.000118 seconds.

tic; X*ksi*Y; toc
Elapsed time is 0.000032 seconds.


Thanks again!
Bruce

> In that case try this:
>
> tot = X*ksi*Y;
>
> It should give the correct answer. Whether it is faster you will have to determine for yourself. You can experiment with parentheses: (X*ksi)*Y as compared with X*(ksi*Y).
>
> Roger Stafford
>

Hi Roger et al.,

I've got another question for you, and I'm not sure if this can also be done without loops.

I'm trying to sum up the squared Euclidean distance between all combinations of points in two sets (X and muEst), multiplied by a scalar value in ksi. Again, a point in X is the nth column and a point in muEst is a column in muEst.

Here's an attempt to write out the equation:
Sum(from n=1 to 4)
   Sum( from i=1 to 4)
     ksi_{ni} ||X_n - muEst_i ||^2
    end
end


and the Corresponding matlab code (which evaluates to tot=186 in this example)
X = [1 2 3 4; 5 4 3 2];
ksi = [ 1 2 2 1; 2 3 3 2; 1 3 5 1; 1 2 2 2];
muEst = [2 3 4 5; 6 5 4 4];
tot = 0;


for n=1: 4,
 for i=1:4,
  tot = tot + ksi(n,i) * sum((X(:,n) - muEst(:,i)).^2);
 end;
end;

Any ideas?
Thanks again,

Bruce

"b b" <bogey4@gmail.com> wrote in message <ge2dju$bjg$1@fred.mathworks.com>...

>
> Any ideas?

Here is one:

X = [1 2 3 4;
     5 4 3 2];
ksi = [ 1 2 2 1;
       2 3 3 2;
       1 3 5 1;
       1 2 2 2];
muEst = [2 3 4 5;
         6 5 4 4];

% Engine
[N I]=ndgrid(1:4,1:4);
tot = sum((X(:,N)-muEst(:,I)).^2,1) * ksi(:)

% Bruno

Thanks Bruno,

I'm still wrapping my head around how this works, but the results are correct. FYI, on a slightly larger dataset (N=137 vs N=4), the times are a bit faster with the loop method, although I was more interested in learning to program without loops, so this is great.


Elapsed time is 0.004989 seconds. (Loop)
Elapsed time is 0.008515 seconds.(ndgrid)


Thanks again,

Bruce
"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ge2jbm$q98$1@fred.mathworks.com>...
> "b b" <bogey4@gmail.com> wrote in message <ge2dju$bjg$1@fred.mathworks.com>...
>
> >
> > Any ideas?
>
> Here is one:
>
> X = [1 2 3 4;
> 5 4 3 2];
> ksi = [ 1 2 2 1;
> 2 3 3 2;
> 1 3 5 1;
> 1 2 2 2];
> muEst = [2 3 4 5;
> 6 5 4 4];
>
> % Engine
> [N I]=ndgrid(1:4,1:4);
> tot = sum((X(:,N)-muEst(:,I)).^2,1) * ksi(:)
>
> % Bruno

"b b" <bogey4@gmail.com> wrote in message <ge2lf8$dbb$1@fred.mathworks.com>...
> Thanks Bruno,
>
> I'm still wrapping my head around how this works, but the results are correct. FYI, on a slightly larger dataset (N=137 vs N=4), the times are a bit faster with the loop method, although I was more interested in learning to program without loops, so this is great.
>
> Elapsed time is 0.004989 seconds. (Loop)
> Elapsed time is 0.008515 seconds.(ndgrid)
>
> Bruce

Hello Bruce,

  Try the following, which expands the square operation:

tot = sum(X.^2,1)*sum(ksi,2)-2*sum(sum((X*ksi).*muEst))...
     +sum(sum(ksi,1).*sum(muEst.^2,1));

  I can't speak for its execution speed, but its intermediate matrices are smaller in size than that produced using the 'ndgrid' method.

Roger Stafford