Vectorizing for loops

I have a function with 3 for loops and wondered if:

1. Can these for loops be vectorized?
2. If yes, could I get some guidance vectorizing them?
3. If no, why not?
4. As a newbie to this forum and MatLab, does the forum consider these appropriate questions to ask?

Function follows:

function[E, L, G] = SetPrincipalComponentBasis(S, A)

% S : covariance matrix
% A : conditioning matrix

if isempty (A)
    N = size (S, 1);
    K = 0;
    [E_, L_] = eig (S);
    E = E_;

      for n = 1 : N
         E (: , n) = E_ (: , N - n + 1);
         L (n) = L_ (N - n + 1, N - n + 1);
      end
else
   [K,N]=size(A);
    E=[];
    B=A;

      for n=1:N-K
          if ~isempty(E) 
              B=[A
                  E'*S];
          end
          e=SetFirstEigenVectors(S,B);
          E=[E e];
      end
      for n=N-K+1:N
          B=E'*S;
          e=SetFirstEigenVectors(S,B);
          E=[E e];
      end
      E=[E(:,N-K+1:N) E(:,1:N-K)];
  end
  L=diag(E'*S*E);
  G=diag(sqrt(L))*inv(E);
  G=G(K+1:N,:);

Some additional info as requested by Walter and Matt...

Sizes of S and A typically less than 30 x 30. Sometimes only 5 X 5. As you discerned they are regular matrices. S & A will always have the same dimensions.

SetFirstEigenVectors code follows:

function e = SetFirstEigenVectors(S,A) N=size(S,1); P=eye(N); if rank(A)>0 P=eye(N)-A'*inv(A*A')*A; end [E_,L_]=eig(P*S*P'); [m,I]=max(diag(L_)); e=E_(:,I);

[Doug Hull]

Moved Andreas' elaboration to the question here (instead of as an answer) and Matt's comment here.

But if A is 5x5, then N-K is 0 and your the first FOR loop in the ELSE clause will not run (for n = 1:0), leaving E empty. This then causes an error in the inner loop when you try to multiply E'*S.

[Matt Fig]

Doug, please format Andreas' second function and put quotes around my comment you moved, or somehow set it off from yours. Thanks.

[Andreas]

Doug's comments led me do some digging through my sources. The code came to me from a friend, who just advised me he got it from: http://www.mathworks.com/matlabcentral/fileexchange/23271

So anyone can look at the original and the zip file has a data.mat file to run an example.

I agree with Doug. The N-K looks funny. Note that the original version has slightly different names for the functions I supplied, GenPCBasis and GenFirstEigVect. You guys go well beyond the call of duty on this forum. Many thanks for any guidance. -- A

[Andreas]

One more quick comment. At this stage of my learning MatLab, I'm more interested in getting a clearer idea of how to go about converting for loops into more vectorized solutions than making this particular set of code work (that may come latter ;-)

Answer 1

[Walter Roberson]

Your loop

      for n = 1 : N
         E (: , n) = E_ (: , N - n + 1);
         L (n) = L_ (N - n + 1, N - n + 1);
      end

is equivalent to

E = fliplr(E_);
L = flipud(fliplr(L_));

provided that your covariance matrix is square. (The square constraint comes out of the fact that you defined N as size(S,1) but use N in the second index rather than the first; defining N as size(S,2) would remove this constraint.)

We cannot determine whether the loops in your "else" condition can be vectorized without seeing the code for SetFirstEigenVectors()

Answer 2

[Matt Fig]

In addition to Walter's comment, this:

diag(sqrt(L))*inv(E)

is equivalent to the preferred:

diag(sqrt(L))/E

Also, it would help if you gave some typical sizes for inputs S and A, as well as any special characteristics they may possess.