Are you referring to the comment in the code that says the following?
% minimum is when gradient of polynomial is zero
At each iteration of the line search you have the parameter value and
the resulting objective function. The idea is that you fit a quadratic
to the values from the previous iterations in order to figure out what
parameter value to try next. After fitting the quadratic, it is
possible to figure out what value for the parameter has the optimal
value for the objective function. Here is a little copy/paste demo...
p = [0.4 0.9 1.2]'; % Last three parameter values.
f = [-0.8 -1.2 0.9]'; % And the resulting objective function
pol = [p.^0 p.^1 p.^2]\f; % Fit a polynomial
d = -pol(2)/(2*pol(3)+eps); % This is the next value to try
pp = (0:0.01:1.4)';
plot(p,f,'o',pp',[pp.^0 pp.^1 pp.^2]*pol,'-',[d d],[-2 -1],'r')
legend('From last 3 iterations','Polynomial Fit','Next value to try');
xlabel('Parameter value');
ylabel('Negative Mutual Information')
Best regards,
-John
On Thu, 2009-08-20 at 21:40 +0100, Henrique Amaral wrote:
> Hi John
>
> Thanks for your attention.
>
> According to Human Brain Function 2Ed, Chapter 2, page 20:
>
> "Section 2.4.1 introduced a method of optimization based on the first and
> second derivatives
> of the cost function. Similar principles have been be applied to minimizing the
> VIR cost function
> [35], and also to maximizing MI [28]7. However, the most widely adopted
> scheme for maximizing
> MI is Powell's Method [see 22, , page 415], which involves a series of
> successive line searches.
> Failures occasionally arise if the voxel similarity measure does not vary
> smoothly with changes
> to the parameter estimates. This can happen because of interpolation artifact,
> or if insucient
> data contributes to the joint histogram. Alternatively, the algorithm can get
> caught within a
> local optimum, so it is important to assign starting estimates that
> approximately register the
> images."
>
> In Numerical Recipes in C, Second Edition (1992),page 415 "Powell´s
> Quadratically Convergent Method", ok.
> But in te spm_powell, the part:
>
> while t(2).f > t(3).f,
>
> % fit a polynomial to t
> tmp = cat(1,t.p)-t(2).p;
> pol = pinv([ones(3,1) tmp tmp.^2])*cat(1,t.f);
>
> % minimum is when gradient of polynomial is zero
> % sign of pol(3) (the 2nd deriv) should be +ve
> d = -pol(2)/(2*pol(3)+eps);
> u.p = t(2).p+d;
>
> ulim = t(2).p+100*(t(3).p-t(2).p);
> if (t(2).p-u.p)*(u.p-t(3).p) > 0.0,
> % u is between t(2) and t(3)
> u.f = linmineval(u.p);
> if u.f < t(3).f,
> % minimum between t(2) and t(3) - done
> t(1) = t(2);
> t(2) = u;
> return;
> elseif u.f > t(2).f,
> % minimum between t(1) and u - done
> t(3) = u;
> return;
> end;
> % try golden search instead
> u.p = t(3).p+gold*(t(3).p-t(2).p);
> u.f = linmineval(u.p);
>
> elseif (t(3).p-u.p)*(u.p-ulim) > 0.0
> % u is between t(3) and ulim
> u.f = linmineval(u.p);
> if u.f < t(3).f,
> % still no minimum as function is still decreasing
> % t(1) = t(2);
> t(2) = t(3);
> t(3) = u;
> u.p = t(3).p+gold*(t(3).p-t(2).p);
> u.f = linmineval(u.p);
> end;
>
> elseif (u.p-ulim)*(ulim-t(3).p) >= 0.0,
> % gone too far - constrain it
> u.p = ulim;
> u.f = linmineval(u.p);
>
> else,
> % try golden search instead
> u.p = t(3).p+gold*(t(3).p-t(2).p);
> u.f = linmineval(u.p);
> end;
>
> % Move all 3 points along
> t(1) = t(2);
> t(2) = t(3);
> t(3) = u;
> end;
> return;
> %_____________
>
> The gradient is calculated in the algorithm. My question is what kind of
> gradient that is.
> According to my leader (teacher), whether the gradient is calculated in
> numerical form,
> it is a descendent gradient,but if calculated analytically, not
> is a descendent gradient. After all, what is kind of gradient is calculated
> here???
>
> Thanks
>
--
John Ashburner <[log in to unmask]>
|