Re: vector function within vector function

[Nicholas Jankowski]

On Wed, Dec 30, 2015 at 6:09 PM, ken1866 <address@hidden> wrote:

I think I understand the error message now.  Since each variable is assumed
to be a matrix, p2 will clearly be incompatible with x.  Does this mean I
have no choice but to use a for loop (like below)?

num=200;
b1=linspace(200,225,num);
b2=linspace(0.33,1.0,num);
[xx,yy]=meshgrid(b1,b2);

function val=ls(p1,p2)
  val=0;
  x=[1,2,3,5,7,10];
  y=[109,149,149,191,213,224];
  #val=sum((y-p1.*(1-exp(-p2.*x))).^2); #<== nonconformant arguments (op1 is
200x200, op2 is 1x6)
  for i=1:6
    val+=(y(i)-p1.*(1-exp(-p2.*x(i)))).^2;
  end
endfunction

z=ls(xx,yy);

Vectorizing is usually possible with straightforward algorithms. x and y have the same length, and you can consider their vector direction as being a projection orthogonal to p1 and p2.  currently, you have p1&2 as 2d matrices (the solution below would need to be adjusted if it would ever higher than 2d). so it would be convenient to project x and y along dim3.  then, you can take advantage of that broadcasting thing I mentioned. Just as a reminder:

look at

>> [1,2,3,4] .* [0.1;1;10]    %%row vector .* col vector

because it sees a size mismatch, it 'projects' or broadcasts the element wise multiplication to create an matrix including both dimensions:

ans =

    0.10000    0.20000    0.30000    0.40000
    1.00000    2.00000    3.00000    4.00000
   10.00000   20.00000   30.00000   40.00000

first row is row vector times first value in the col vectior, 2nd row is by the 2nd value, etc.

the same would happen for a 2d matrix .* a dim3 vector. try:

>> a = magic(3)

>> b = reshape([0.1,1,10],1,1,3)

>> a .* b

So, given that behavior of the . operators in Octave (matlab doesn't do that, you have to use the 'bsxfun' function instead), you can replace the for loop by calculating all of the intermediate values of val in one step, and then summing. So, just replace the three lines for the for loop with the following:

----
x = reshape(x,1,1,[]);
y = reshape(y,1,1,[]);
 
valall =   (y - p1.* (1 - exp(-p2.*x))).^2;

val = sum(valall);
----

you can actually combine those last two lines getting rid of the intermediate valall variable. 

Also, note that this is just one way of vectorizing your code. Using projection in higher dimensions often works, but it can quickly run into memory problems with larger data sets (you're basically multiplying your element count...).  But it's often straightforward for small problems, and I find it easy to keep straight. There are a lot of tips and tricks out there for vectorizing matlab code, and I'm sure there is a more elegant way to do this one.

nick j