Summing anonymous functions in a for loop: adding penalty terms for optimization

I have the following piece of code for passing into lsqnonlin:

    baseFun = @(p) ((x-p(1))*cos(p(3))+(y-p(2))*sin(p(3))).^2/a^2+...
        ((x-p(1))*sin(p(3))+(y-p(2))*cos(p(3))).^2/b^2-1;
    
    totalFun = baseFun;
    % Penalty stiffness
    lam = 1e5;
    j = 1;
    for k = whichVis
        totalFun = @(p) totalFun(p) + lam*(x(j) - p(1) - xeye(1,k)*cos(p(3)) + xeye(2,k)*sin(p(3)))^2+lam*(y(j)-p(2)-xeye(1,k)*sin(p(3)) -xeye(2,k)*cos(p(3)))^2;
        j = j + 1;
    end

I am trying to rigidly fit an ellipse of known size in the plane to a set of markers. I know how the markers should relate to the local x-y coordinates of the ellipse (x aligned along semi-major, y aligned along semi-minor). Sometimes, these markers are invisible. Hence, I only add penalty terms for markers indexed in whichVis.

p is a 3X1 array of ellipse center x, ellipse center y, and alpha, the inclination of the ellipse. These are the parameters to be optimized.

baseFun is an anonymous function, the equation of an ellipse in the plane.

totalFun is the entire anonymous function to be passed into lsqnonlin, with penalty terms added.

lam is a penalty stiffness.

xeye is a 2X8 array of known marker positions in local x-y coordinates of the ellipse.

x and y are marker positions of the markers which are visible.

The indices in x and y are ordered appropriately to correspond to the indicies in whichVis.

Am I constructing my objective function properly? When I print its value the j and k indices are unevaluated. I expect p to be the only thing to remain symbolic. I would have expected MATLAB to find the values in the arrays xeye, x, and y and fill them in as numbers.

Sample arrays to run code:

  xeye = [31.6016  -15.3975   24.3983   42.1441  -35.7844    0.0596  -27.5122    5.8206
   20.2474   23.6097  -22.9670   -3.5785   17.3939  -26.7238  -20.1340   26.0083]
  
  whichVis = [1     2     3     4     5     6     7     8]

x = [-237.4388
 -283.9550
 -219.5598
 -218.5919
 -291.4497
 -238.4623
 -265.7021
 -263.3386]

y = [-141.6660
 -160.0548
 -181.3641
 -158.6064
 -178.2309
 -196.2995
 -200.7612
 -147.0627]

a = 44.5499
b = 26.5382

The code executes but I am skeptical that it's set up properly.

Thank you.

Answer 1

[Torsten]

basefun returns a vector of size (8x1) . Now to each component of this vector you add the sum of the following scalar values

lam*(x(1) - p(1) - xeye(1,1)*cos(p(3)) + xeye(2,1)*sin(p(3)))^2+lam*(y(1)-p(2)-xeye(1,1)*sin(p(3)) -xeye(2,1)*cos(p(3)))^2;
lam*(x(2) - p(1) - xeye(1,2)*cos(p(3)) + xeye(2,2)*sin(p(3)))^2+lam*(y(2)-p(2)-xeye(1,2)*sin(p(3)) -xeye(2,2)*cos(p(3)))^2;
...
lam*(x(8) - p(1) - xeye(1,8)*cos(p(3)) + xeye(2,8)*sin(p(3)))^2+lam*(y(8)-p(2)-xeye(1,8)*sin(p(3)) -xeye(2,8)*cos(p(3)))^2;

Of course, x, y and xeye have to be defined before setting up the function handle.

If this is what you want, your code is correct.

[Evan Hemingway]

Thanks, I've changed the code to the following:

    baseFun = @(p) ((x-p(1))*cos(p(3))+(y-p(2))*sin(p(3))).^2/a^2+...
        ((x-p(1))*sin(p(3))+(y-p(2))*cos(p(3))).^2/b^2-1;
    
    % Penalty stiffness
    lam = 1e5;
    totalFun = @(p) baseFun(p) + lam*(x - p(1) - xeye(1,whichVis)'*cos(p(3)) + xeye(2,whichVis)'*sin(p(3))).^2+lam*(y-p(2)-xeye(1,whichVis)'*sin(p(3)) -xeye(2,whichVis)'*cos(p(3))).^2;

Now I do not have an unevaluated k and j in totalFun.

[Torsten]

I tested a recursive update of a function handle in a loop as you did and you are right: it doesn't work this way.

But the way you did it above with one line is alright.

[Evan Hemingway]

Ironically, the results were better with my original code. I think what may have been happening before was that maybe my objective function included only one spring associated with the final k and j value. Now it has a tough time converging properly. Perhaps I should lower my penalty stiffness. Anyway, thanks for your help.