Prev: Optimization of a discontinuous function
Next: hi i need this my deadline in is in 10 hours help
From: Kinshuk on 22 Jun 2010 16:38
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)))]);
From: Matt J on 22 Jun 2010 16:59
KJ <kinshuk(a)cmu.edu> wrote in message <2019808340.10319.1277237474473.JavaMail.root(a)gallium.mathforum.org>...

>
> 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.
===============

If you supply your own function to approximate the Hessian, MATLAB will use that. Otherwise, it will use finite differences to approximate it, and will probably result in a crazy value. In either case, fmincon is not intended to handle non-smooth functions and cannot be trusted to function well near points of discontinuity. If it has worked for you empirically in the past, that's only because you were lucky enough that the algorithm did not land on/near that 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?
=================

Probably not. A formula involving a Hessian probably assumes that the Hessian exists, an assumption that is not valid in your case.
From: Matt J on 22 Jun 2010 17:04
"Kinshuk " <kinshuk(a)cmu.edu> wrote in message <hvr6vt$q1u$1(a)fred.mathworks.com>...

>
> 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 change my mind. Yes, the formula can work, provided you compute the Hessian sufficiently accurately at x=1.000000000025769

Note that a finite difference approx. is a dangerous way to go about that considering how close x=1.000000000025769 lies to the discontinuity.
From: Kinshuk on 26 Jun 2010 18:26
Thanks for your help, Matt.


"Matt J " <mattjacREMOVE(a)THISieee.spam> wrote in message <hvr8g5$6pf$1(a)fred.mathworks.com>...
> "Kinshuk " <kinshuk(a)cmu.edu> wrote in message <hvr6vt$q1u$1(a)fred.mathworks.com>...
>
> >
> > 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 change my mind. Yes, the formula can work, provided you compute the Hessian sufficiently accurately at x=1.000000000025769
>
> Note that a finite difference approx. is a dangerous way to go about that considering how close x=1.000000000025769 lies to the discontinuity.
 | 
Pages: 1
Prev: Optimization of a discontinuous function
Next: hi i need this my deadline in is in 10 hours help