vectorizing my code

can anyone give an idea on how to vectorize my code.

for loop=1:1:10
          %X=1*loop*ones(n,4);
          %Y=1*loop*ones(n,4)';
          X=rand(n,4);
          Y=rand(4,n);
          U=zeros(n,4);
          V=zeros(4,n);
          while(norm(U-X)+norm(V-Y)>0.01)
              U=X;
              V=Y;
              for i=1:1:n
                  f=@(r)norm((D(i,:)-[r(1),r(2),r(3),r(4)]*Y));
                  X(i,:)=fminsearch(f,[1,1,1,1]);
              end

                for j=1:1:n
                    g=@(s)norm((D(:,j)-X*[s(1);s(2);s(3);s(4)]));
                    Y(:,j)=fminsearch(g, [1;1;1;1]);
                end
                if(max(max(X*Y))>3)
                    X=rand(n,4)/10;
                    Y=rand(4,n)/10;
                    U=zeros(n,4);
                    V=zeros(4,n);
                end
            end

end

Thanks

[Roger Stafford]

A comment:

This is a linear least squares problem. You should be using matlab's matrix left and right divisions rather than 'fminsearch'. This is the kind of problem they are designed to do. Do:

X(i,:) = (D(i,:)*Y')/(Y*Y');

and

Y(:,j) = (X'*X)\(X'*D(:,j));

instead of

f=@(r)norm((D(i,:)-[r(1),r(2),r(3),r(4)]*Y));

X(i,:)=fminsearch(f,[1,1,1,1]);

and

g=@(s)norm((D(:,j)-X*[s(1);s(2);s(3);s(4)]));

Y(:,j)=fminsearch(g,[1;1;1;1]);

That should save quite a lot of cpu time.

(Note: I believe you should be saving copies of X to use for the purpose of subsequently changing Y before you begin the process of changing X. Otherwise you get strange results.)

Roger Stafford

[Roger Stafford]

I should have said;

X = (D*Y')/(Y*Y');

and

Y = (X'*X)\(X'*D);

in place of the two entire inner for-loops. No need to separately compute the least squares rows or columns in X or Y, respectively.

[Roger Stafford]

Third (embarrassed) comment:

My brain is not up to par today, I'm afraid. What I should have said was the very simple:

X = D/Y;

and

Y = X\D;

in place of the two inner for loops.

As I wrote these earlier I would have been doing some of the matrix divisions' work for them. My apologies.

[srinadh]

let me try what you said, and I will let you know the results.

Thanks

Answer 1

[Walter Roberson]

I think you had better explain what you are attempting to accomplish.

[srinadh]

I am trying to rewrite the code without the outer for loop

[bym]

I don't see an _outer_ for loop

[Mark Shore]

(Walter's comment)^3

Answer 2

[Walter Roberson]

The code as it is now does not use the value of the index variable within the body of the "for" loop. If we examine the body of the "for" loop we see that it has no memory, no reliance on previous iterations: no matter what has gone before, it re-initializes all of the variables. The only memory that that loop has is that it uses different random number values in each iteration. Is it the case that it is important to "advance" the random number generator by a particular amount, or is one set of random numbers as good as another for your purposes? If one set is as good as another, then the way to eliminate the "for" loop is simply to perform what is inside the loop exactly once: that will generate the same answer (statistically speaking) as running it with 10 different sets of random numbers and keeping only the information about the last set.

[srinadh]

Thanks for your replies. Sorry for not giving more details before. Let me explain a bit more what I am trying to do.

I have to run the code in the for loop for different starting values of X and Y. I am not looking for averages. I am checking convergence of the minimization from different starting points.

Main issue in vectorizing the code is that the while loop termination condition for different runs of X and Y. I am not sure how to write it.

[srinadh]

Also I figured out using profiler that most of the time is taken b y fminsearch function. So i guess vectorizing wont help much.

Thanks guys for your help.