Thread Subject:

Speed up the nested loops (slower than inverse a large matrix)

I have a large matrix correlation R (n by n), a design matrix F (n by k) and I want to calculate

the sum of F(i,j)^t *R(i,j)^(-1)*F(i,j) for all possible pair i,j , where F(i,j) is a 2byk matrix with the ith, jth row of F ; R(i,j) is a 2by2matrix with entries of the intersection of i,j th row and i,jth column of R.

I use the for loops for i=1:n-1
                           for j=i+1:n
which turns out very slow compared with inv(R) which is theoretically O(n^3) operation.

Any one know how to speed up the loops in my case ? Is it possible faster than the buil-in function inv(R)? Thanks a lot .


Simon

 

Post more code, specifically what is in the FOR loop.
Also, would backslash help in this case?

"Ximing xu" <siemenxu@gmail.com> wrote in message <icpphd$cb2$1@fred.mathworks.com>...
> .......
> the sum of F(i,j)^t *R(i,j)^(-1)*F(i,j) for all possible pair i,j , where F(i,j) is a 2byk matrix with the ith, jth row of F ; R(i,j) is a 2by2matrix with entries of the intersection of i,j th row and i,jth column of R.
> .......
> I use the for loops for i=1:n-1
> for j=i+1:n
> .......
- - - - - - - - - -
  I am puzzled by your description "sum of F(i,j)^t *R(i,j)^(-1)*F(i,j) for all possible pair i,j". The range of j in F is presumably different from the range of j in R^-1, so how can we speak of using "all possible pair i,j"? Can you please clarify that statement and how can we make it compatible with the given nested for loops in which i ranges from 1 to n-1 and j ranges from i+1 to n?

  I wonder if you are actually referring to the matrix product F.'*(R^-1)*F? If so, are you asking for the sum of all its elements, and if not, what kind of a sum do you have in mind? If you answer Matt's question thoroughly, that will probably answer mine.

Roger Stafford

Thanks for your quick reply.
F is a k*n matrix , R is a n*n matrix. the loop is
A=zeros(k,k);
for i=1:n-1
for j=i+1:n
 A=A+F([i,j],:)'/ R([i,j],[i,j])*F([i,j],:) ;
end
end
A

I hope this code is clear. I hope this loop should be faster than inverse of R, but for n=100 and 1000, inverse of R is much faster than this loop.

Thanks,

Ximing

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <icpt5k$2q$1@fred.mathworks.com>...
> "Ximing xu" <siemenxu@gmail.com> wrote in message <icpphd$cb2$1@fred.mathworks.com>...
> > .......
> > the sum of F(i,j)^t *R(i,j)^(-1)*F(i,j) for all possible pair i,j , where F(i,j) is a 2byk matrix with the ith, jth row of F ; R(i,j) is a 2by2matrix with entries of the intersection of i,j th row and i,jth column of R.
> > .......
> > I use the for loops for i=1:n-1
> > for j=i+1:n
> > .......
> - - - - - - - - - -
> I am puzzled by your description "sum of F(i,j)^t *R(i,j)^(-1)*F(i,j) for all possible pair i,j". The range of j in F is presumably different from the range of j in R^-1, so how can we speak of using "all possible pair i,j"? Can you please clarify that statement and how can we make it compatible with the given nested for loops in which i ranges from 1 to n-1 and j ranges from i+1 to n?
>
> I wonder if you are actually referring to the matrix product F.'*(R^-1)*F? If so, are you asking for the sum of all its elements, and if not, what kind of a sum do you have in mind? If you answer Matt's question thoroughly, that will probably answer mine.
>
> Roger Stafford

"Matt Fig" <spamanon@yahoo.com> wrote in message <icpqj3$hu3$1@fred.mathworks.com>...
> Post more code, specifically what is in the FOR loop.
> Also, would backslash help in this case?
Thanks for your quick reply.
F is a k*n matrix , R is a n*n matrix. the loop is
A=zeros(k,k);
for i=1:n-1
for j=i+1:n
 A=A+F([i,j],:)'/ R([i,j],[i,j])*F([i,j],:) ;
end
end
A

I hope this code is clear (my original complete code is in another computer). I hope this loop should be faster than inverse of R, but for n=100 and 1000, inverse of R is much faster than this loop.

Thanks,

Ximing

"Ximing xu" <siemenxu@gmail.com> wrote in message <icq218$3ib$1@fred.mathworks.com>...
 
> F is a k*n matrix , R is a n*n matrix. the loop is
> A=zeros(k,k);
> for i=1:n-1
> for j=i+1:n
> A=A+F([i,j],:)'/ R([i,j],[i,j])*F([i,j],:) ;
> end
> end
> A
>
> I hope this code is clear (my original complete code is in another computer). I hope this loop should be faster than inverse of R, but for n=100 and 1000, inverse of R is much faster than this loop.

What you might try is to extract once all the 2x2 submatrices of R and put it on 3D array, and use this FEX to solve the system than sum them back:

http://www.mathworks.com/matlabcentral/fileexchange/27762-small-size-linear-solver

Bruno

"Ximing xu" <siemenxu@gmail.com> wrote in message <icq1vb$t6d$1@fred.mathworks.com>...
> Thanks for your quick reply.
> F is a k*n matrix , R is a n*n matrix. the loop is
> A=zeros(k,k);
> for i=1:n-1
> for j=i+1:n
> A=A+F([i,j],:)'/ R([i,j],[i,j])*F([i,j],:) ;
> end
> end
> A
>
> I hope this code is clear. I hope this loop should be faster than inverse of R, but for n=100 and 1000, inverse of R is much faster than this loop.
>
> Thanks,
>
> Ximing
- - - - - - - - - -
  I am still puzzled, this time by your statement "I hope this loop should be faster than inverse of R, but for n=100 and 1000, inverse of R is much faster than this loop." What does the inverse of the n by n matrix R have to do with this calculation? That is a very different operation from computing the inverses of individual pairs of two rows and columns in R. Why are you comparing them? For each trip within your nested for-loops you are doing an order(k^2) algorithm because of all the pairs of rows of F, which makes the over-all algorithm order(k^2*n^2), but for large k this is principally occupied with computing the k^2 elements of A rather than just finding the inverse of the 2 by 2 matrix R([i,j],[i,j]). This latter should be very fast in comparison. In fact it is quite possible that using inv(R([i,j],[i,j])) would be faster than using the backslash operator as you have done
here. After all, the inverse of a 2 x 2 correlation matrix [1,a;a,1] is simply:

 [1/(1-a^2),-a/(1-a^2);-a/(1-a^2),1/(1-a^2)]

This is a far cry from computing the inverse of an entire n x n correlation matrix, but this in turn is a far cry from doing your entire operation with the F matrix to obtain the k x k matrix A.

Roger Stafford

Ximing, yes, I agree with Roger than your calculation of sum over 2x2 matrices are very odd. May be there is some rough approximation that leads you to compute such sum. But it looks suspiciously unjustified IMO.

Bruno

I try to compare two statistical methods: the old method and my new method. For the old method I need to calculate the inverse of R , but for the new method I only need to calculate the inverse of 2by2 matrix. The code I show you is just part of my method. The motivation for my method is that the sum shown in my code is only O(n^2), which is much faster than inversion of R which is O(n^3).

The k is very small. I also tried to write the he inverse of a 2 x 2 correlation matrix [1,a;a,1] as [1/(1-a^2),-a/(1-a^2);-a/(1-a^2),1/(1-a^2)]; but the loop is still very slow.

I also tried Bruno's idea to extract all the 2x2 submatrices of R (also the F) and put it on 3D array; but it turn out the procedure to do the extraction is just as slow as doing the sum directly.

Thanks,

Ximing

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <icqhka$6og$1@fred.mathworks.com>...
> "Ximing xu" <siemenxu@gmail.com> wrote in message <icq1vb$t6d$1@fred.mathworks.com>...
> > Thanks for your quick reply.
> > F is a k*n matrix , R is a n*n matrix. the loop is
> > A=zeros(k,k);
> > for i=1:n-1
> > for j=i+1:n
> > A=A+F([i,j],:)'/ R([i,j],[i,j])*F([i,j],:) ;
> > end
> > end
> > A
> >
> > I hope this code is clear. I hope this loop should be faster than inverse of R, but for n=100 and 1000, inverse of R is much faster than this loop.
> >
> > Thanks,
> >
> > Ximing
> - - - - - - - - - -
> I am still puzzled, this time by your statement "I hope this loop should be faster than inverse of R, but for n=100 and 1000, inverse of R is much faster than this loop." What does the inverse of the n by n matrix R have to do with this calculation? That is a very different operation from computing the inverses of individual pairs of two rows and columns in R. Why are you comparing them? For each trip within your nested for-loops you are doing an order(k^2) algorithm because of all the pairs of rows of F, which makes the over-all algorithm order(k^2*n^2), but for large k this is principally occupied with computing the k^2 elements of A rather than just finding the inverse of the 2 by 2 matrix R([i,j],[i,j]). This latter should be very fast in comparison. In fact it is quite possible that using inv(R([i,j],[i,j])) would be faster than using the backslash operator as you have
done
> here. After all, the inverse of a 2 x 2 correlation matrix [1,a;a,1] is simply:
>
> [1/(1-a^2),-a/(1-a^2);-a/(1-a^2),1/(1-a^2)]
>
> This is a far cry from computing the inverse of an entire n x n correlation matrix, but this in turn is a far cry from doing your entire operation with the F matrix to obtain the k x k matrix A.
>
> Roger Stafford

"Ximing xu" <siemenxu@gmail.com> wrote in message <icretg$4gm$1@fred.mathworks.com>...

> I also tried Bruno's idea to extract all the 2x2 submatrices of R (also the F) and put it on 3D array; but it turn out the procedure to do the extraction is just as slow as doing the sum directly.
>

I guess timing depends *how* you do the extraction. A for-loop would not be optimal.

Bruno

Thanks, Bruno. I agree with you; I think using for-loop is not efficient to do the extraction or to do the sum. Is there any (possible) alternative to do the extraction in my case?

If I want to compare the speed of my method (double sums) and the traditional method (inverse R) more precisely , should I use some basic language like C to do the programming rather than use the built-in functions in Matlab?

Thanks,

Ximing


"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <icrjun$oe5$1@fred.mathworks.com>...
> "Ximing xu" <siemenxu@gmail.com> wrote in message <icretg$4gm$1@fred.mathworks.com>...
>
> > I also tried Bruno's idea to extract all the 2x2 submatrices of R (also the F) and put it on 3D array; but it turn out the procedure to do the extraction is just as slow as doing the sum directly.
> >
>
> I guess timing depends *how* you do the extraction. A for-loop would not be optimal.
>
> Bruno

"Ximing xu" <siemenxu@gmail.com> wrote in message <icuh3v$i51$1@fred.mathworks.com>...
> Thanks, Bruno. I agree with you; I think using for-loop is not efficient to do the extraction or to do the sum. Is there any (possible) alternative to do the extraction in my case?
>
> If I want to compare the speed of my method (double sums) and the traditional method (inverse R) more precisely , should I use some basic language like C to do the programming rather than use the built-in functions in Matlab?
>

I think you should, because the traditional method of inversion of R benefits a very efficient implementation of BLAS/LAPACK.

Your method should be implemented in low level language to be fair. But even that I'm not optimistic it can win on traditional processor because your method needs data to be moved a lot in the memory and that takes time. However it can be efficiently parallelized such as GPU computing.

Good luck,

Bruno