Complex number on optimization

3 次查看（过去 30 天）
Camelot Abeits
Camelot Abeits 2017-12-9
编辑： Matt J 2017-12-9
Hi! I have a code for the simmulated annealing method, to maximize a function. My function is the log likelihood. The problem is that the annealing reports a maximizer that when you evaluate on the log likelihood, it gets you an complex number! How can i put in the code a line that says something like this:
"If the loglikelihood is complex, then try again"
Thanks a lot! Here's the code:
fu
nction [minimum,fval] = anneal(loss, parent, options)
% ANNEAL Minimizes a function with the method of simulated annealing
% (Kirkpatrick et al., 1983)
%
% ANNEAL takes three input parameters, in this order:
%
% LOSS is a function handle (anonymous function or inline) with a loss
% function, which may be of any type, and needn't be continuous. It does,
% however, need to return a single value.
%
% PARENT is a vector with initial guess parameters. You must input an
% initial guess.
%
% OPTIONS is a structure with settings for the simulated annealing. If no
% OPTIONS structure is provided, ANNEAL uses a default structure. OPTIONS
% can contain any or all of the following fields (missing fields are
% filled with default values):
%
% Verbosity: Controls output to the screen.
% 0 suppresses all output
% 1 gives final report only [default]
% 2 gives temperature changes and final report
% Generator: Generates a new solution from an old one.
% Any function handle that takes a solution as input and
% gives a valid solution (i.e. some point in the solution
% space) as output.
% The default function generates a row vector which
% slightly differs from the input vector in one element:
% @(x) (x+(randperm(length(x))==length(x))*randn/100)
% Other examples of possible solution generators:
% @(x) (rand(3,1)): Picks a random point in the unit cube
% @(x) (ceil([9 5].*rand(2,1))): Picks a point in a 9-by-5
% discrete grid
% InitTemp: The initial temperature, can be any positive number.
% Default is 1.
% StopTemp: Temperature at which to stop, can be any positive number
% smaller than InitTemp.
% Default is 1e-8.
% StopVal: Value at which to stop immediately, can be any output of
% LOSS that is sufficiently low for you.
% Default is -Inf.
% CoolSched: Generates a new temperature from the previous one.
% Any function handle that takes a scalar as input and
% returns a smaller but positive scalar as output.
% Default is @(T) (.8*T)
% MaxConsRej: Maximum number of consecutive rejections, can be any
% positive number.
% Default is 1000.
% MaxTries: Maximum number of tries within one temperature, can be
% any positive number.
% Default is 300.
% MaxSuccess: Maximum number of successes within one temperature, can
% be any positive number.
% Default is 20.
%
%
% Usage:
% [MINIMUM,FVAL] = ANNEAL(LOSS,NEWSOL,[OPTIONS]);
% MINIMUM is the solution which generated the smallest encountered
% value when input into LOSS.
% FVAL is the value of the LOSS function evaluated at MINIMUM.
% OPTIONS = ANNEAL();
% OPTIONS is the default options structure.
%
%
% Example:
% The so-called "six-hump camelback" function has several local minima
% in the range -3<=x<=3 and -2<=y<=2. It has two global minima, namely
% f(-0.0898,0.7126) = f(0.0898,-0.7126) = -1.0316. We can define and
% minimise it as follows:
% camel = @(x,y)(4-2.1*x.^2+x.^4/3).*x.^2+x.*y+4*(y.^2-1).*y.^2;
% loss = @(p)camel(p(1),p(2));
% [x f] = ANNEAL(loss,[0 0])
% We get output:
% Initial temperature: 1
% Final temperature: 3.21388e-007
% Consecutive rejections: 1027
% Number of function calls: 6220
% Total final loss: -1.03163
% x =
% -0.0899 0.7127
% f =
% -1.0316
% Which reasonably approximates the analytical global minimum (note
% that due to randomness, your results will likely not be exactly the
% same).
% Reference:
% Kirkpatrick, S., Gelatt, C.D., & Vecchi, M.P. (1983). Optimization by
% Simulated Annealing. _Science, 220_, 671-680.
% joachim.vandekerckhove@psy.kuleuven.be
% $Revision: v5 $ $Date: 2006/04/26 12:54:04 $
def = struct(...
'CoolSched',@(T) (.9*T),...
'Generator',@(x) (x+(randperm(length(x))'==length(x))*randn/100),...
'InitTemp',1,...
'MaxConsRej',50,... %1000 default
'MaxSuccess',20,... %20 default
'MaxTries',120,... %300 default
'StopTemp',1e-8,...
'StopVal',-Inf,...
'Verbosity',1);
% Check input
if ~nargin %user wants default options, give it and stop
minimum = def;
return
elseif nargin<2, %user gave only objective function, throw error
error('MATLAB:anneal:noParent','You need to input a first guess.');
elseif nargin<3, %user gave no options structure, use default
options=def;
else %user gave all input, check if options structure is complete
if ~isstruct(options)
error('MATLAB:anneal:badOptions',...
'Input argument ''options'' is not a structure')
end
fs = {'CoolSched','Generator','InitTemp','MaxConsRej',...
'MaxSuccess','MaxTries','StopTemp','StopVal','Verbosity'};
for nm=1:length(fs)
if ~isfield(options,fs{nm}), options.(fs{nm}) = def.(fs{nm}); end
end
end
% main settings
newsol = options.Generator; % neighborhood space function
Tinit = options.InitTemp; % initial temp
minT = options.StopTemp; % stopping temp
cool = options.CoolSched; % annealing schedule
minF = options.StopVal;
max_consec_rejections = options.MaxConsRej;
max_try = options.MaxTries;
max_success = options.MaxSuccess;
report = options.Verbosity;
k = 1; % boltzmann constant
% counters etc
itry = 0;
success = 0;
finished = 0;
consec = 0;
T = Tinit;
initenergy = loss(parent);
oldenergy = initenergy;
total = 0;
if report==2, fprintf(1,'\n T = %9.7f, loss = %10.5f\n',T,oldenergy); end
while ~finished;
itry = itry+1; % just an iteration counter
current = parent;
% % Stop / decrement T criteria
if itry >= max_try || success >= max_success;
if T < minT || consec >= max_consec_rejections;
finished = 1;
total = total + itry;
break;
else
T = cool(T); % decrease T according to cooling schedule
if report==2, % output
fprintf(1,' T = %9.7f, loss = %10.5f\n',T,oldenergy);
end
total = total + itry;
itry = 1;
success = 1;
end
end
newparam = newsol(current);
newenergy = loss(newparam);
if (newenergy < minF),
parent = newparam;
oldenergy = newenergy;
break
end
if (oldenergy-newenergy > 1e-6) %default: 1e-6
parent = newparam;
oldenergy = newenergy;
success = success+1;
consec = 0;
else
if (rand < exp( (oldenergy-newenergy)/(k*T) ));
parent = newparam;
oldenergy = newenergy;
success = success+1;
else
consec = consec+1;
end
end
end
minimum = parent;
fval = oldenergy;
if report;
fprintf(1, '\n Initial temperature: \t%g\n', Tinit);
fprintf(1, ' Final temperature: \t%g\n', T);
fprintf(1, ' Consecutive rejections: \t%i\n', consec);
fprintf(1, ' Number of function calls:\t%i\n', total);
fprintf(1, ' Total final loss: \t%g\n', fval);
end

请先登录，再回答此问题。

回答（2 个）

John D'Errico
John D'Errico 2017-12-9
编辑：John D'Errico 2017-12-9
It is just simulated annealing. It won't be harmed if you just slap a penalty on the objective. Return inf to the annealing code when a complex result is generated.
Star Strider
Star Strider 2017-12-9
I did not even attempt to go through all your code.
"If the loglikelihood is complex, then try again"
Test to see if the imaginary part is zero. If it not, the test for a complex result is true.
Example —
X = 1
IsComplex = imag(X) ~= 0
X = 1 + 2i
IsComplex = imag(X) ~= 0

类别

在 Help Center 和 File Exchange 中查找有关 Simulated Annealing 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by