%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% Global Optimization Sketch -------------------------------------------
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
%% Issue with fmin
% Consider the function y = e^{-.1x}*sin(x) defined on x \in [0,\infty).
% Over its support [0,\infty), this function has a global minimum at
% x=4.613. Additionally, it has many local minima. For example, x=10.896 is
% a local minimizer.
% First, use fmincon to find the global minimum of this function.
% Initial guess %
param0 = 4;
% Upper & Lower bounds of search space %
lb = 0;
ub = inf;
fun=@(param)y(param);
[solution,~] = fmincon(fun,param0,lb,ub);
solution
% The solution is 4.613, so it found the global minimizer. But I sort of
% cheated. I picked an initial guess very close to the global minimizer,
% which I knew beforehand.
% Let me instead pick a point close to the local minimizer x=10.896
param0 = 10;
[solution,~] = fmincon(fun,param0,lb,ub);
solution
% Now, it reports a solution of 10.896. Fmincon, or any of the "fmins,"
% will always report back the first minima it finds, whether it is global
% or local.
%% Global Minimization Algorithms
% To avoid local minima, we can use a global optimizer. Simulated annealing
% and particle swarm are two covenient choices. Below is some quick code
% for running both. For both, I'll use initial guess param0=10.896 just to
% show that even with a poor initial guess, these algorithms will still
% find the global minimizer (at least most of the time).
param0 = 10.896;
fun=@(param)y(param);
% Upper & Lower bounds of search space %
lb = 0;
ub = inf;
%%% Simulated Annealing %%%
options = optimoptions(@simulannealbnd,'Display','diagnose','InitialTemperature',400);
[solution,~] = simulannealbnd(fun,param0,lb,ub);
solution
%%% Particle Swarm %%%
x0 = repmat(param0,100,1);
options = optimoptions(@particleswarm,'SwarmSize',100,'Display','iter','InitialSwarmMatrix',x0);
[solution,~] = particleswarm(fun,1,lb,ub);
solution
% Run this section over and over again and you'll see that most of the
% time, both algorithms find the global minimizer despite me initializing
% it right at a local minimizer. Occasionally it gets stuck in local
% minima. There are a lot of settings that can be tweaked to prevent this
% from happening.
%% Parallelizing Particle Swarm
% For computationally-intensive models, particle swarm is preferable to
% simulated annealing. The reason is that PS can be parallelized while SA
% cannot. Below I'll show how to run PS in parallel.
% # of cpus to be used %
cpu=4; % I set # of cpus = # of function evaluations in each iteration.
% Note that this number cannot be larger than the # of cpus on your
% computer, unless you're running on longleaf. When
% running on longleaf, it is good to have the # of cpus be a multiple of
% 12.
% Upper & Lower bounds of search space %
lb = 0;
ub = inf;
% Initial guess %
% PS doesn't need an initial guess. If the swarm size (# of
% particles on each iteration) is N, then PS will by default use N times 1
% points uniformly spaced over the support of the function to be minimized
% as its initial guess.
% If you would like to manually pick an initial guess, you can first pick a point:
param0 = 7;
% Turn this point into a vector, where each element of the vector will be
% evaluated by a separate particle:
x0=repmat(param0,[cpu,1]); % matrix of initial guesses (one guess per particle)
% Then slightly perturb each guess (so they aren't all evaluating at the
% same point):
x0=x0.*(1+.2*(randi(3,size(x0))-2));
% The above line randomly scales up or down (or leaves unchanged) all of
% the elements of x0.
%parpool('local', cpu); %manually create parallel pool (not usually needed)
fun=@(param)y(param);
warning('off','all');
options = optimoptions(@particleswarm,'UseParallel',true,'SwarmSize',cpu,'Display','iter','InitialSwarmMatrix',x0);
[solution,~] = particleswarm(fun,size(param0,2),lb,ub,options);
solution
% If you run this over and over again, you'll see that most of the time it
% successfully finds the global minimizer x=4.613. Again, it occasionally
% settles in local minima. To avoid this, it is good to make the swarm size
% as large as possible. I only used 4 particles in this example, so the coverage
% of the search space is pretty horrible. In practice you should use as many
% particles as is computationally feasible.
%% Functions
% y = e^{-.1x}*sin(x)
function obj = y(param)
obj = exp(-.1*param).*sin(param);
end