MATLAB Central Newsreader - Vectorization of matrix operations--need guru's help!

Vectorization of matrix operations--need guru's help!

Jerry — Mon, 24 Nov 2008 17:59:03 -0500

is there any way to improve the speed of the following code? I need to call it 10000 times. thanks!

function m4=calm4(mwk)
[n0,n1]=size(mwk);
m4=0;
for i=1:n0
    tmp1=mwk(:,i);
    for ii=1:n0
        if ii~=i
            tmp2=mwk(:,ii)*tmp1';
            m4=m4+sum(sum(tmp2))-sum(diag(tmp2));
        end
    end
end
m4=m4/(n0*n1*(n0-1)*(n1-1));

Re: Vectorization of matrix operations--need guru's help!

Roger Stafford — Mon, 24 Nov 2008 19:08:02 -0500

"Jerry" <mricad@yahoo.no000spppam.com> wrote in message <ggeq17$hj$1@fred.mathworks.com>...
> is there any way to improve the speed of the following code? I need to call it 10000 times. thanks!
>
> function m4=calm4(mwk)
> [n0,n1]=size(mwk);
> m4=0;
> for i=1:n0
> tmp1=mwk(:,i);
> for ii=1:n0
> if ii~=i
> tmp2=mwk(:,ii)*tmp1';
> m4=m4+sum(sum(tmp2))-sum(diag(tmp2));
> end
> end
> end
> m4=m4/(n0*n1*(n0-1)*(n1-1));

  Jerry, notice that 'tmp2' consists of every possible product between elements of the i_th and ii_th columns of 'mwk'. When you do sum(sum(tmp2)) you are taking the sum of all these possible products and by the distributive law of arithmetic that is the same as the product of the two column sums, which is a lot easier to compute.

  Notice further that by similar reasoning, if we temporarily remove the restriction that ii~=i and take the sum of all possible values of these column-sum products above, this is the same as the product of the sum of all column-sums by itself, namely the square of the sum of all elements of 'mwk', which again would be an enormous simplification.

  In the ii==i cases the quantity 'sum(sum(tmp2))' is simply the square of the sum of the ii_th (i_th) column, and the sum of all these would be the sum of the squares of all the column sums, again a simply quantity to compute in compensating for their erroneous inclusion in the previous paragraph.

  The quantity 'sum(diag(tmp2))' is simply the dot product of the ii_th column by the i_th column. If we sum for every possible combination of ii and i, again ignoring the ii~=i constraint, this gives the dot product of the row-sums of 'mwk' by itself, which again is much easier to find.

  Finally, for the ii==i cases 'sum(diag(tmp2))' is the dot product of the ii_th (i_th) column by itself. The sum of these must then be subtracted to compensate for their erroneous inclusion above.

  I leave it to you to implement these notions.

  There seems to be some confusion in your code as to whether n0 is the size of the row or the column of 'mwk'.

Roger Stafford

Re: Vectorization of matrix operations--need guru's help!

Roger Stafford — Mon, 24 Nov 2008 20:19:02 -0500

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <ggeu2i$1l4$1@fred.mathworks.com>...
> ......
> I leave it to you to implement these notions.
> ......

Despite my intention to leave the implementation to you, curiosity got the better of me and I decided to see what the resulting code would look like. Here it is. See if you agree with it.

m4 = (sum(sum(mwk,1))^2-sum(sum(mwk,1).^2)-sum(sum(mwk,2).^2) ...
+sum(sum(mwk.^2,1)))/(n0*n1*(n0-1)*(n1-1));

Roger Stafford

Re: Vectorization of matrix operations--need guru's help!

Jerry — Mon, 24 Nov 2008 22:34:02 -0500

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message
>
> Despite my intention to leave the implementation to you, curiosity got the better of me and I decided to see what the resulting code would look like. Here it is. See if you agree with it.
>
> m4 = (sum(sum(mwk,1))^2-sum(sum(mwk,1).^2)-sum(sum(mwk,2).^2) ...
> +sum(sum(mwk.^2,1)))/(n0*n1*(n0-1)*(n1-1));
>
> Roger Stafford

YES! And it's over 100 times faster! Thanks a million!
I was trying to implement your ideas but made a mistake in the last term.

Re: Vectorization of matrix operations--need guru's help!

Jerry — Mon, 24 Nov 2008 22:38:02 -0500

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message
> There seems to be some confusion in your code as to whether n0 is the size of the row or the column of 'mwk'.
>
> Roger Stafford

You are right. should be
[n1,n0]=size(mwk);