Prev: Data Exchange Between Matlab and Third Party Software like
Next: Optimization of a discontinuous function
From: KJ on 22 Jun 2010 12:10
I am facing the following problem. Consider the following function test_func(x):

if x <= 0 y = 0;
elseif x <= 1 y = 10000*x^2;
elseif x <= sqrt(3) y = 5005*(3 - x^2);
else y = 0;
end

This function evalutes to 10000 at x=1 and 10010 at x=1+\epsilon (i.e., just above 1). If I maximize this function (basically,

use fmincon to minimize the -ve of this function), I obtain the correct maximum of 10010 and also the value 1.000000000025769.

Matlab uses the active-set (line search) method. Along with the result above, Matlab also gives me the Hessian (simply a scalar

here in this case) which evaluates to 534605.7. Note, however, that x=1 is a point of discontinuity in the function.

My first question is: How does Matlab compute the Hessian at this point of discontinuity? I have read in detail about the

active-set method from Matlab's manual and from other sources, but could not find out the details of computation at a

discontinuous point.

Now, suppose test_func(x) is a likelihood function so that x=1.000000000025769 is the maximum likelihood estimate. My second

question is: Given that x=1 is a point of discontinuity, can the formula sqrt(inv(hessian)) be used to calculate the standard

error of the estimate? If not, what is the correct way to calculate this standard error?

I will greatly appreciate any pointers to any resources (manuals, examples or theory).

Thanks a lot.

kjmatlab.


Below, I have attached the code I used.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% test_func.m %%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function y = test_func(x);

global SCALE DELTA

if x <= 0
y = 0;
elseif x <= 1
y = 10000*x^2 * SCALE;
elseif x <= sqrt(3)
y = (5000 + DELTA/2)*(3 - x^2) * SCALE;
else
y = 0;
end

y = -y;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% optimize.m %%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global SCALE DELTA

DELTA = 10;
SCALE = 1;

%full plot
argvec = -1:0.001:3;
for i=1:length(argvec)
valvec(i) = -test_func(argvec(i));
end
subplot(1, 2, 1)
scatter(argvec, valvec, 1);

%zoomed plot
argvec_fine = 0.9995:0.00001:1.0005;
for i=1:length(argvec_fine)
valvec_fine(i) = -test_func(argvec_fine(i));
end
subplot(1, 2, 2)
pos_of_1 = find(argvec_fine == 1);
plot(argvec_fine(1:pos_of_1), valvec_fine(1:pos_of_1));
hold on
plot(argvec_fine((pos_of_1+1):length(argvec_fine)), valvec_fine((pos_of_1+1):length(argvec_fine)));
diff = max(valvec_fine) - min(valvec_fine);
ylow = min(valvec_fine) - diff/2;
yup = max(valvec_fine) + diff/2;
axis([0.9995 1.0005 ylow yup])
plot(1, -test_func(1), 'ro', 'MarkerEdgeColor', 'b', 'MarkerFaceColor', 'b', 'MarkerSize', 8);
plot(1, -test_func(1.00001), 'ro', 'MarkerEdgeColor', 'b', 'MarkerFaceColor', 'w', 'MarkerSize', 8);


initial = [0.5];
lb = [-2];
ub = [ 4];

options = optimset('Display', 'final', 'MaxFunEvals', 25000);
[param, fmax, exitflag, output, lambda, grad, hessian] = fmincon(@test_func, initial, [], [], [], [], lb, ub, [], options);

format long
disp(['Maximum = ' num2str(-fmax) ' at x = ' num2str(param, '%.15f')]);
disp(['The hessian = ' num2str(hessian) ' and the std err computed as sqrt(inv(hessian)) is = ' num2str(sqrt(inv(hessian)))]);
 | 
Pages: 1
Prev: Data Exchange Between Matlab and Third Party Software like
Next: Optimization of a discontinuous function