Supplementary material for“Optima: a model for HIV epidemic analysis,
program prioritization, and resource optimization”

Contents

HIV epidemic model structure .................................................................................................................................. 2
Optimization algorithm............................................................................................................................................... 51

HIV epidemic model structure

The aim of the HIV transmission and progression model that underlies Optima is to estimate the movementof people between population groups and between health states over time. Each combination of populationgroup and health state is a single compartment in the model, and these evolve in time according to a set ofdifferential equations. Seven types of movement are possible: people can enter the model by being born orleave the model by dying, or they can move between compartments by becoming infected, having their CD4
count decline, becoming diagnosed, or by going on or off treatment. Each compartment is described by asingle ordinary differential equation, and each of the seven types of movement described above is a term in adifferential equation. Positive terms result in an increase in population in a given compartment, while negative
terms reduce the population of a given compartment. Movement between population groups (e.g., fromfemale sex workers to low-risk females and vice versa; from male children to male adolescents, etc.) is also
possible.
The model dynamics are determined by the initial conditions (i.e., the number of people in each compartment
at the initial time point) and the model parameter values. Model parameters can depend on time 𝑡 (e.g., 2002
vs. 2007), population group 𝑔 (e.g., female sex workers vs. male injecting drug users), health state ℎ (e.g.
infected, CD4>500 vs. infected, CD4<200), and interaction type 𝑖 (i.e., casual, regular, commercial, or
injecting partnerships). Optima supports up to 14 distinct population groups; since there are 21 health states(susceptible plus undiagnosed, diagnosed, first-line treatment, treatment failure, and second-line treatment foreach of four CD4 stratifications), the model may have up to 294 compartments.
Transmission dynamics are determined by the rate at which susceptible (i.e., uninfected) individuals become
infected. The probability that a susceptible individual in population group 𝑔1 will be infected by an individualof population group 𝑔2 at time 𝑡 via sexual intercourse of interaction type 𝑖 is defined as the force-ofinfection:𝑚(𝑡,𝑔1 ,𝑔2 ,𝑖)𝜆(𝑡, 𝑔1 , 𝑔2 , 𝑖) = 1 − (1 − 𝜍(𝑡, 𝑔1 )𝜓(𝑡, 𝑔1 )𝜌(𝑡, 𝑔1 )𝜚(𝑡, 𝑔1 )𝐼(𝑡, 𝑔2 , 𝑖)),where 𝜆 is the force-of-infection; 𝜍 is the effect of circumcision; 𝜓 is the effect of coinfection with
discharging and ulcerative sexually transmitted infections (STIs); 𝜌 is the effect of pre-exposure prophylaxis(PrEP); 𝜚 is the effect of post-exposure prophylaxis (PEP); 𝐼 is the infectiousness, which is determined byprevalence and transmissibility; and 𝑚 is the number of unprotected acts.
The effect of circumcision is defined as 𝜍(𝑡, 𝑔) = 1 − 𝑒𝜍 𝑝𝜍 (𝑡, 𝑔), where 𝑒𝜍 is the fractional reduction intransmission probability per act due to circumcision (i.e., its efficacy), and 𝑝𝜍 is the probability that anindividual is circumcised. This probability is zero for female populations, so circumcision only affects malepopulations. The effect of STI coinfection is 𝜓(𝑡, 𝑔) = 𝑒𝜓 𝑝𝜓 (𝑡, 𝑔), where 𝑒𝜓 is the fractional increase intransmissibility caused by ulcerative STIs (such as HSV-2 and syphilis), and 𝑝𝜓 is the probability of infection(i.e., prevalence) of these STIs. The effects of PrEP and PEP are defined as 𝜌(𝑡, 𝑔) = 1 − 𝑒𝜌 𝑝𝜌 (𝑡, 𝑔) and𝜚(𝑡, 𝑔) = 1 − 𝑒𝜚 𝑝𝜚 (𝑡, 𝑔), which is similar to circumcision, except that changes in transmissibility due to
PrEP and PEP apply to both male and female populations.
The infectiousness of a population is𝐼(𝑡, 𝑔) = ∑ 𝑃(𝑡, 𝑔, ℎ)𝛽(ℎ) ,ℎwhere 𝐼 is the infectiousness, 𝛽 is the biological transmission probability, which depends on whether theintercourse is insertive, vaginal receptive, or anal receptive; 𝑃 is the prevalence of each health state (such that
𝑃(𝑡, 𝑔) = ∑ℎ 𝑃(𝑡, 𝑔, ℎ) is the prevalence of HIV in that population); and 𝛽(ℎ) is the biological transmission
2

probability, which depends on the viral load, which in turn can be approximated from the individual’s CD4
count category ℎ.
The number of unprotected acts is 𝑚(𝑡, 𝑔1 , 𝑔2 , 𝑖) = 𝑛(𝑡, 𝑔1 , 𝑔2 , 𝑖)(1 − 𝑒𝜅 𝑝𝜅 (𝑡, 𝑔1 )), where 𝑛 is the numberof acts between individuals in two populations, 𝑒𝜅 is the fractional reduction in transmission probability peract due to condom use, and 𝑝𝜅 is the probability of condom use. While it may seem odd that the effect ofcondom use appears in the exponent while the effect of circumcision appears in the base term, circumcisionaffects the biological transmission probability and thus affects the risk per act, whereas condom use can beconsidered as a reduction in the number of acts. The two approaches are equivalent to first order, and since
the transmission probability is typically <10−3, it makes negligible difference whether these terms appear in
the base or exponent.
The force-of-infection due to injecting drug use is
𝜔(𝑡)𝑚𝑑 (𝑡)𝜆𝑑 (𝑡) = 1 − (1 − 𝜒(𝑡)𝐼𝑑 (𝑡)),where 𝜆𝑑 is the force-of-infection for injection interactions between injecting drug users, 𝜒 is the effect ofsyringe cleaning, 𝐼𝑑 is the infectiousness of injecting drug users, 𝜔 is the effect of methadone treatment, and𝑚𝑑 is the number of receptively shared injections.
The effect of syringe cleaning is 𝜒(𝑡) = 1 − 𝑒𝜒 𝑝𝜒 (𝑡), where 𝑒𝜒 is the fractional reduction in transmission
probability per act, and 𝑝𝜒 is the probability that a syringe is cleaned before use. The infectiousness 𝐼𝑑 has thesame form as before, but now includes only injecting drug users. The effect of methadone is 𝜔(𝑡) = 1 −𝑒𝜔 𝑝𝜔 (𝑡), where 𝑒𝜔 is the fractional reduction in the number of injections, and 𝑝𝜔 is the probability that agiven injecting drug user is on methadone. The number of receptively-shared injections is given by 𝑚𝑑 (𝑡) =𝑛𝑑 (𝑡)𝑠(𝑡), where 𝑛𝑑 is the total number of injections, and 𝑠 is the probability that a given syringe has been
receptively shared prior to use.
Finally, the force-of-infection for a given population group 𝑔1 at a given point in time 𝑡 is the product of all
individual forces-of-infection (including injecting interactions, if applicable):Λ(𝑡, 𝑔1 ) = 1 − ∏ ∏(1 − 𝜆(𝑡, 𝑔1 , 𝑔2 , 𝑖))𝑔2𝑖This quantity is the instantaneous risk of an individual in population group 𝑔1 becoming infected. The indices𝑔2 and 𝑖 include all interactions an individual has. For example, direct female sex workers have commercial
interactions with male clients, as well as regular and casual relationships with all four heterosexual malepopulations, for a total of nine interaction types.
In contrast to the force-of-infection, most types of movement between population compartments aredescribed by a single rate, and thus the remaining equations are relatively straightforward. Labeling the healthstates ℎ such that 0=acute, 1=CD4>500, 2=350<CD4<500, 3=200<CD4<350, 4=CD4<200, 5=CD4<50,
the full set of differential equations is as follows.
The change in the number of susceptible individuals is𝑑𝑆(𝑡, 𝑔)= ε(𝑡, 𝑔) − (Λ(𝑡, 𝑔) + 𝜇0𝑔 )𝑆(𝑡, 𝑔),𝑑𝑡
3

where 𝑆 is the number of susceptible individuals, 𝜖 is the entry rate into the population, 𝜆 is the force-ofinfection, and 𝜇0𝑔 is the background mortality rate for population group 𝑔. Note that the term Λ𝑆 in thisequation is the number of new infections, while 𝜆𝑆 would be the incidence.
The change in the number of infected but undiagnosed individuals is
𝑑𝑈ℎ=0 (𝑡, 𝑔)= Λ(𝑡, 𝑔) 𝑆(𝑡, 𝑔) − (𝜇ℎ𝑔 + 𝜁ℎ (𝑡, 𝑔) + 𝜋ℎ )𝑈ℎ (𝑡, 𝑔),𝑑𝑡
𝑑𝑈ℎ>0 (𝑡, 𝑔)= 𝜋ℎ−1 𝑈ℎ−1 (𝑡, 𝑔) − (𝜇ℎ𝑔 + 𝜁ℎ (𝑡, 𝑔) + 𝜋ℎ )𝑈ℎ (𝑡, 𝑔),𝑑𝑡
where 𝑈ℎ is the number of undiagnosed individuals in health state ℎ (note that newly infected individuals Λ𝑆
only enter 𝑈ℎ=0 , the CD4>500 compartment), 𝜇ℎ𝑔 is the mortality rate for the given health state andpopulation group, 𝜁ℎ is the testing rate, and 𝜋ℎ is the disease progression rate (𝜋ℎ=5 = 0, since noprogression occurs after CD4<50).
The change in the number of diagnosed individuals is
𝑑𝐷ℎ (𝑡, 𝑔)= 𝜁ℎ (𝑡, 𝑔)𝑈ℎ (𝑡, 𝑔) + 𝜋ℎ−1 𝐷ℎ−1 (𝑡, 𝑔) − (𝜇ℎ𝑔 + 𝜏ℎ (𝑡, 𝑔) + 𝜋ℎ )𝐷ℎ (𝑡, 𝑔),𝑑𝑡
where 𝐷ℎ is the number of infected and diagnosed individuals in a particular health state ℎ, and 𝜏ℎ is the HIVtreatment rate.
The change in the number of individuals on treatment or with treatment failure is
𝑑𝑇1 (𝑡, 𝑔)= ∑ 𝜎ℎ (𝑡, 𝑔)𝐷ℎ (𝑡, 𝑔) − (𝜇𝐴 + 𝜙1 )𝑇1 (𝑡, 𝑔),𝑑𝑡ℎ𝑑𝑇2 (𝑡, 𝑔)= 𝜎𝐹 (𝑡, 𝑔)𝑇𝐹 (𝑡, 𝑔) − (𝜇𝐴 + 𝜙2 )𝑇2 (𝑡, 𝑔),𝑑𝑡
𝑑𝑇𝐹 (𝑡, 𝑔)= 𝜙1 𝑇1 (𝑡, 𝑔) + 𝜙2 𝑇2 (𝑡, 𝑔) − (𝜇𝐹 + 𝜎𝐹 (𝑡, 𝑔))𝑇𝐹 (𝑡, 𝑔),𝑑𝑡
where 𝑇1 , 𝑇2 , and 𝑇𝐹 are the numbers of people on first- and second-line antiretroviral therapy (ART) andtreatment failure, respectively; 𝜇𝐴 and 𝜇𝐹 are the mortality rates for individuals on ART and with treatment
failure, respectively; and 𝜙1 and 𝜙2 are the failure rates for first- and subsequent-lines of ART, respectively.
This set of coupled differential equations can be solved numerically; in practice, maximum computational
efficiency is achieved by converting these to difference equations with a suitable time step (usually 0.1–1
years). Running the MATLAB implementation of this model for 10 population groups for a 50-year timespan with an 0.2-year time step on an Intel Xeon 2.4 GHz CPU takes approximately 2 s.

4

Optimization algorithm
The following manuscript has been submitted to Nature Scientific Reports and is currently under review.

5

1

Optimization by Bayesian adaptive locally linear
stochastic descent
Cliff C. Kerr1,2,3 *, Tomasz G. Smolinski4 , Salvador Dura-Bernal2 , David P. Wilson1
1

Kirby Institute for Infection and Immunity in Society, University of New South Wales, Sydney, NSW, Australia
Department of Physiology and Pharmacology, SUNY Downstate Medical Center, Brooklyn, NY
3
Complex Systems Group, School of Physics, University of Sydney, Sydney, NSW, Australia
4
Department of Computer and Information Sciences, Delaware State University, Dover, DE
2

Abstract
When standard optimization methods fail to find a satisfactory solution for a parameter fitting problem, a tempting recourse
is to adjust parameters manually. While tedious, this approach can be surprisingly powerful in terms of achieving optimal ornear-optimal solutions. This paper outlines an optimization algorithm, Bayesian Adaptive Locally Linear Stochastic Descent(BALLSD), that has been designed to replicate the essential aspects of manual parameter fitting in an automated way. Specifically,BALLSD uses simple Bayesian principles to form probabilistic assumptions about (a) which parameters have the greatest effecton the objective function, and (b) optimal step sizes for each parameter. We show that for a certain class of optimization problems(namely, those with a moderate to large number of scalar parameter dimensions, especially if some dimensions are more important
than others), BALLSD is capable of minimizing the objective function with far fewer function evaluations than classic optimization
methods, such as the Nelder-Mead nonlinear simplex, Levenberg-Marquardt gradient descent, simulated annealing, and genetic
algorithms.

1 Introduction
Consider a human H who is attempting to minimize a nonlinear objective function, E = f (x), by manually adjusting
parameters in the vector x. H typically begins with a uniform prior regarding which parameters to vary, and chooses step
sizes that are a fixed fraction (e.g., 10%) of the initial parameter values. H will then pseudorandomly choose one or moreparameters to adjust. Every time a parameter xi is found to reduce E, the probability that H will select xi in the future
increases; conversely, if changes in xi are not found to improve E, the probability that H will select xi decreases (formally,
H forms “hunches” about which parameters are “good”). H also adaptively adjusts the step size based on the information
H obtains about the curvature of parameter space with respect to each dimension (e.g., if ∆E/∆xi ≈ const. over multiple
iterations, H will increase the step size). Despite its drawbacks, the adaptive nature of manual parameter fitting makes it aremarkably powerful method.
Thus, despite the smörgåsbord of available automated optimization algorithms, manual fitting of parameters remains a familiar
bane of researchers (e.g., Castillo et al. (2005); Brom et al. (2012)), especially in cases where evaluations of the objectivefunction are computationally intensive, such as climate models (Wilby, 2005), neuronal network models (Prinz et al., 2003;
Baker et al., 2011; Kerr et al., 2013), or detailed epidemiological models (Kwon et al., 2012). However, it is difficult toestimate how commonly manual parameter fitting is performed, since authors often do not explicitly mention its use (e.g.,Song et al. (2013)).
In many types of optimization problems, it is more important to need only a small number of function evaluations to find areasonable local minimum than it is to find the global minimum (Goodner et al., 2012). Indeed, the latter may be ill-defined
given the large uncertainties that are often present when models of complex systems are fitted to empirical data, as in thecitations listed above.
With the increasing availability of high-performance computers and clusters (Hilbert and López, 2011), easily parallelizable
optimization methods such as evolutionary algorithms (where different individuals can be run on different cores) and Monte
Carlo methods (where different initializations can be run on different cores) have a notable advantage for certain types ofproblems. The common theme in these algorithms is the ability to use a different random seed for each parallel instance.
However, as the size of parameter space increases, the advantage of this approach is diluted: whereas a 3- or even 5-dimensionalparameter space may be reasonably densely sampled by a Monte Carlo initialization, a 20- or 100-dimensional space cannot.
This is because parameter space grows exponentially with an increasing number of dimensions, whereas parallelization increases
sampling rates linearly.
In high-dimensional parameter spaces, it is unlikely that all parameters contribute equally to the objective function. Identifying
those that contribute more, thereby allowing computational resources to be focused on them, has the potential to significantly* Corresponding author:
Address: 450 Clarkson Ave., Brooklyn, NY 11203, USA
Phone: +1 (347) 721-7328
Email: cliff@thekerrlab.com

2

reduce the total number of function evaluations required. Despite humans’ limited capacity to implement Bayesian-optimal
strategies (Charness et al., 2007; Steyvers et al., 2009), we speculate that this adaptive approach to both parameter selection
and step size is the key reason why manual parameter fitting can be highly effective.
In this paper, we present an algorithm, Bayesian Adaptive Locally Linear Stochastic Descent (BALLSD), that was inspired bymanual parameter fitting. This approach is most applicable to optimization problems with more than approximately 5 dimensions– i.e., large enough so that performing function evaluations across all dimensions is inefficient. The paper is organized asfollows: Sec. 2 describes the basic BALLSD algorithm; Sec. 3.1 quantifies the performance of BALLSD against four common
nonlinear optimization methods (Nelder-Mead nonlinear simplex, Levenberg-Marquardt gradient descent, simulated annealing,
and a genetic algorithm) for several classic optimization problems; Sec. 3.2 provides a real-world example of optimizing HIV
resource allocations using BALLSD; Sec. 4.1 describes possible modifications to the algorithm that may make it more suitablefor particular types of problems; and Sec. 4.2 summarizes the advantages and limitations of BALLSD.

2 Methods
Consider an objective function E = f (x), where E is the scalar error (or other quantity) to be minimized (or maximized)
and x = [x1 , x2 , ..., xn ] is an n-element vector of parameters. There are 2n possible directions j to step in: an increase or−
decrease in the value of each parameter. Associated with each parameter xi are (a) two initial step sizes: sj = s+i or si , which+−
define the step size in the directions of increasing or decreasing xi , respectively (i.e., si > 0 and si < 0); and (b) two initial−probabilities: pj = p+
which define the likelihood of selecting direction j (for a uniform prior, pj = 1/2n – satisfying
iPor pi ,P2n
p = j=1 pj = 1). Thus, the vectors s and p have length 2n.
the requirement that
At each step k, the algorithm maps a random variable α ∈ (0, 1) onto p, thereby choosing a direction j ∈ (1 . . . 2n) and acorresponding parameter i = dj/2e ∈ (1 . . . n), where d·e denotes the ceiling operator. The algorithm then evaluates
Ek± = f (x + δ(i)),(1)

where δ(i) is an n-element vector such that δi = sj and 0 otherwise. Then:
1) If Ek± < Ek−1 :a) The new parameter value is adopted: xi → xi + sj ;b) The error is updated: Ek → Ek± ;c) sj is increased: sj → sj · sinc (sinc > 1);Pd) pj is increased: pj → pj · pinc (pinc > 1), and p is renormalized such thatp = 1.
2) Otherwise:a) The parameter vector x and error E are not changed;b) sj is decreased: sj → sj /sdec (sdec > 1);c) pj is decreased: pj → pj /pdec (pdec > 1), and p is renormalized as above.
The algorithm thus has four meta-parameters: sinc , sdec , pinc , and pdec . In practice, a reasonable choice is sinc = sdec =pinc = pdec = 2, although in our experience any value from approximately 1.2 to 3 appears to work reasonably well for thetest cases used here; in general, the smoother and more linear the objective function is, the larger the learning rates shouldbe. In addition to these meta-parameters, four initial value vectors need to be specified: the initial parameter vector x0 , step
sizes s (which in general can be initialized as a fixed fraction of the corresponding initial parameter value), and probabilities
p (where typically pj = 1/2n suffices for an n-parameter problem).
By modifying s and p after each iteration, the algorithm learns which directions are most effective to step in and by how
much, using a loosely Bayesian approach (in the sense that it estimates the posterior distributions of s and p by modifying thepriors according to accumulated evidence). This, combined with the stochastic choice of which parameters to modify on each
iteration, resembles the way in which humans (imperfectly) perform Bayesian decision-making in situations such as N -armed
bandit problems (Steyvers et al., 2009).

3 Results
3.1 Comparison to other optimization methods
Here we compare BALLSD to four standard optimization methods: the Nelder-Mead nonlinear simplex algorithm (Nelder
and Mead, 1965), Levenberg-Marquardt gradient descent (Marquardt, 1963), simulated annealing (Kirkpatrick et al., 1983), anda genetic algorithm (Bethke, 1978). All methods were implemented in Matlab 2012b (The Mathworks) via the Optimization
Toolbox functions “fminsearch”, “lsqnonlin”, “simulannealbnd”, and “ga”, respectively.
For BALLSD, we used meta-parameters sinc = pinc = sdec = pdec = 2, initial step sizes sj of 20% of the parametervalues in x0 (which are given below; the step size for any parameter with an initial value of 0 is the mean of the otherstep sizes), and initial probabilities pj of 1/2n for an n-dimensional problem. Matlab’s default meta-parameters were used
for the other four algorithms, except that the initial temperature of the simulated annealing algorithm was set to be equal to3

Rosenbrock's valley (2D)

10

-2

10

-4

10

-6

10-8

0

2000

4000

6000

Powell's quartic (4D)100

Relative error

Relative error

100

8000

10

-2

10

-4

10

-6

10-8

10000

0

2000

Function evaluations
Powell's quartic (12D)

0

10

10-2

Relative error

Relative error

10

10-4
10

-6

10-8

0

2000

4000

6000

Relative error

10000

8000

10000

8000

10000

10-4
10

-6

10-8

0

2000

4000

6000

Function evaluations

Powell's quartic (100D)

0

BALLSD
Simplex-4

Levenberg-Marquardt
Simulated annealing

10-6
10

8000

10-2

10-2
10

6000

Powell's quartic (20D)0Function evaluations
10

4000

Function evaluations

Genetic algorithm

-8

0

2000

4000

6000

8000

10000

Function evaluations

F IGURE 1: Performance of BALLSD compared to four standard nonlinear optimization algorithms: Nelder-Mead nonlinear simplex,
Levenberg-Marquardt gradient descent, simulated annealing, and a genetic algorithm. While standard methods – especially the simplex method– are most efficient for low-dimensional problems (e.g., Rosenbrock’s valley), in many cases BALLSD is the most efficient algorithm forhigh-dimensional parameter spaces (e.g., the 100-dimensional version of Powell’s quartic function). For the stochastic methods (BALLSD,
simulated annealing, and the genetic algorithm), the shaded regions show the interquartile range for 40 different random seeds.

10 · h|x0 |i following manual exploration of meta-parameter space, since the default choice of 100 did not generalize well acrossproblems of different scales. Indeed, one of the major disadvantages of this type of algorithm is its sensitivity to the values of
its meta-parameters (Ben-Ameur, 2004).
To test this suite of algorithms, we used original and modified versions of the two classic optimization problems given in
Nelder and Mead (1965):
1) Rosenbrock’s parabolic valley (two-dimensional):
E = 100(x2 − x21 )2 + (1 − x1 )2 ,(2)with the starting point at x = (−1.2, 1). The optimum is at x = (1, 1).
2) A modified 10-dimensional version of Rosenbrock’s valley, with the functional form as given in Eq. 2, but with a 10element parameter vector x; the remaining 8 parameters do not contribute to the objective function. The starting point
is at x = (1.5, −1.5, 0, 0 . . . 0). The optimum is at x = (0, 0 . . . 0).
3) A modified Powell’s quartic function (N -dimensional):XE=(x1 + 10x2 )2 + 5(x3 − x4 )2 + (x2 − 2x3 )4 + 10(x1 − x4 )4 ,(3)

4

F IGURE 2: Optimization of the 10-dimensional version of Rosenbrock’s valley starting from the point (1.5, −1.5). (A) Trajectories of each
optimization method up to 300 function evaluations; each iteration is shown with a square, but note that multiple function evaluations may
occur at each iteration. Note the locally linear steps of BALLSD that rapidly adapt in size. (B) Relative error of each method for the first
100 function evaluations, showing the initial stage of the algorithms. (C) Relative error for the first 300 function evaluations, showing theasymptotic stage of the algorithms.

where xq is a vector of length N/4. The starting point is at x10 = (3), x20 = (−1), x30 = (0), and x40 = (1),where each component is repeated N/4 times (e.g., if N = 8, x0 = (3, 3, −1, −, 1, 0, 0, 1, 1)). The optimum is at
x = (x1 , x2 , x3 , x4 ) = (0, 0, 0, . . . 0). Here, we used 4, 12, 20, and 100-dimensional versions of Powell’s function.
As shown in Fig. 1, for the two-dimensional optimization problem, the nonlinear simplex method is most efficient, with allother algorithms requiring considerably more function evaluations to obtain the same error. Notably, BALLSD was especiallyinefficient, since its assumption of local linearity is violated by the shallow, curved valley. However, with the modified version of
the problem described above, BALLSD is the most efficient algorithm over most of the first several hundred function evaluations,as shown in Fig. 2 for a single random seed. For small numbers of iterations (<30), for this particular seed, simulated annealing
was by far the most efficient algorithm, reducing the error by a remarkable 98% after just 4 function evaluations. However, thisalgorithm became mired near the point (1.5, 2.4), far from the true minimum of (1, 1), and did not significantly reduce the error
beyond the first 20 function evaluations. After 50 function evaluations, BALLSD had reduced the error by 99.9%, compared to
99.7% for simulated annealing, 96% for the Levenberg-Marquardt method, 82% for the nonlinear simplex method, and 0% for5Simplex
4

3

3
Parameter value

Parameter value

BALLSD
4

2
1
0
−1
−2

2
1
0
−1

0

500

1000
1500
Function evaluations

−2

2000

0

500

4

3

3

2
1
0
−1
−2

2
1
0
−1

0

500

1000
1500
Function evaluations

−2

2000

Genetic algorithm

0

500

1000
1500
Function evaluations

2000

Objective function

0

4

10

3

−2

10
Relative error

Parameter value

2000

Simulated annealing

4

Parameter value

Parameter value

Levenberg−Marquardt

1000
1500
Function evaluations

2
1
0

−4

10

BALLSD
SimplexLevenberg−Marquardt
Simulated annealing
Genetic algorithm

−6

10
−1
−2

−8

0

500

1000
1500
Function evaluations

2000

10

0

500

1000
1500
Function evaluations

2000

F IGURE 3: Parameter update strategies for each algorithm applied to a 20-dimensional Powell’s quartic function. Each line is a separate
parameter; the optimum is at (0, 0, 0, . . . 0). The error relative to the starting point is shown in the bottom right panel. For small numbers of
iterations (the adaptive phase of BALLSD), the Levenberg-Marquardt method reduces error most quickly; for larger numbers of iterations,BALLSD achieves 1–4 orders of magnitude smaller error for a given number of iterations than the other methods. (Note: since the genetic
algorithm does not use a single initial point, individuals were instead initialized using a uniform random distribution in the range [−1, 3]. The
Levenberg-Marquardt algorithm operates on the 20-dimension Powell’s function identically to the 4-dimensional version, with the exceptionthat each iteration requires 5 times as many function evaluations.)

the genetic algorithm. Similarly, BALLSD reduced the error by 99.99% after 70 function evaluations; in comparison, the next
best algorithm (the simplex method) required 220 function evaluations to reach the same error level. During the descent into
the shallow curved valley (comprising ∼99.9% of the total error), the most efficient algorithms were BALLSD and simulated
annealing; within the valley (the remaining ∼0.1% of the total error), the simplex algorithm was by far the most efficient.
Hence, this example illustrates that in optimization problems where some parameters are significantly more important thanothers, BALLSD has significant advantages. In contrast, for problems in which all parameters have equal importance, as in
the original Rosenbrock’s valley problem, other algorithms have superior performance.

6

For the 4-dimensional Powell’s quartic function, the nonlinear simplex method was again the most efficient, followed byBALLSD. For the 12- and 20-dimensional version, BALLSD was most efficient for 60–1700 and 250–4400 function evaluations
respectively (corresponding to roughly 99.9999% of the total error at the upper limit in each case), after which the simplex
method was most efficient. For the 100-dimensional version, the Levenberg-Marquardt method was most efficient for the first
1000 function evaluations (corresponding to 97% of the total error), but BALLSD was the most efficient algorithm for larger
numbers of function evaluations.
The five optimization methods discussed here employ very different parameter update strategies, as shown strikingly in
Fig. 3. The approach used in BALLSD is most similar to the Levenberg-Marquardt method, with the exception that the rate ofconvergence of the former increases over time (due to its adaptive step size), whereas for the latter it decreases. In the example
shown here (a 20-dimensional Powell’s quartic function), the Levenberg-Marquardt method has the lowest error for 250 or
fewer iterations; for large numbers of iterations, BALLSD has by far the lowest error – indeed, for 2000 or more iterations, it
has nearly 2 orders of magnitude less error than the Levenberg-Marquardt method, and 4 orders of magnitude less error thannonlinear simplex, simulated annealing, and genetic algorithms. The superior performance of BALLSD compared to the othermethods is surprising since, unlike in Fig. 2, in this problem all parameters are of roughly equal importance, so the adaptive
probability p is unlikely to significantly contribute to the efficiency of the optimization. Thus, even in cases where BALLSD’s
only advantage is its adaptive step size, it is still capable of outperforming traditional algorithms.
3.2 Optimizing HIV resource allocations
In contrast to the foregoing theoretical discussion of error minimization for analytical functions, here we describe the practicalapplication that BALLSD was designed for, a topic that has recently garnered considerable interest (Anderson et al., 2014):
finding the allocation of resources across different HIV prevention programs that minimizes new infections. Full details of
the HIV model are presented in Supplementary Methods S1. In brief, the model describes HIV transmission and progression
in a number of interacting subpopulations (14 in this case), including female sex workers, men who have sex with men, andgeneral males and females in different age groups. The model incorporates roughly 200 parameters, which describe the sexual
behavior, injecting behavior, HIV testing and treatment rates, and sexual and injecting partnerships of each population, as wellas basic clinical parameters such as HIV transmissibility and disease progression rates. In addition to empirical estimates ofthese parameters, the model is calibrated to match surveillance data on HIV prevalence, diagnoses, and numbers of people ontreatment. In this example, the model was based on and calibrated to behavioral and surveillance data from Swaziland. Since
the model is relatively computationally intensive, requiring 1–2 s per function evaluation on a standard laptop, large numbers(>102 ) of evaluations are wearisome.
To optimize the allocation of Swaziland’s HIV budget, we assumed that spending on particular HIV programs produceschanges in corresponding behavioral parameters or testing and treatment rates (for example, programs targeting female sex
workers increase their probability of condom use). In this case, the objective we chose to minimize was the number of new
infections over the period 2015–2020, subject to the constraint that total funding was held constant at current (2014) levels.
Allocations across 9 different HIV prevention, testing, and treatment programs were considered.
As shown in Fig. 4A, under current conditions, the model predicts an average of approximately 2500 new infections per
year in Swaziland. However, if funding is optimally allocated, as shown in Fig. 4B (which consists largely of shifting funds
from programs for orphans and vulnerable children towards treatment and male circumcision programs), this can be reduced
to approximately 1260 new infections per year. BALLSD found this allocation after 65 function evaluations. The second bestalgorithm, the Levenberg-Marquardt method, found a nearly identical allocation after 830 function evaluations. None of theother methods reached this level of optimization within 2000 function evaluations; by that point, the genetic algorithm had
achieved 99.3% of the reduction in new infections found by BALLSD and the Levenberg-Marquardt method, the nonlinear
simplex algorithm 95%, and the simulated annealing algorithm 90%.

4 Discussion
4.1 Extensions of the algorithm
This section describes several modifications to the basic algorithm that may make it more suitable for a broader range of
optimization problems.
First, the assumptions of uniform priors p and uniform initial step sizes s can easily be relaxed, allowing assumptions about
the scale or relative importance of parameters to be incorporated. However, due to the adaptive nature of the algorithm, evensilly initial choices of p and s will be corrected, as long as all pj and sj are nonzero.
Second, rather than updating the probability pj by a fixed amount after each successful iteration, the change in pj (∆pj )could be proportional to the change in the error function E (∆E), with a larger ∆E resulting in larger ∆pj . Since the expectedchange in E at step k is proportional to both Ek − min(E) and the ratio of the step size to the characteristic scale of eachparameter, and since in general neither of these quantities are known, the constant of proportionality between ∆pj and ∆E
cannot typically be estimated a priori. One can partially circumvent this problem by comparing the current ∆E to its previous

7

A

2500

New infections per year

BALLSD
SimplexLevenberg-Marquardt
Simulated annealing
Genetic algorithm
2000

1500

1000

B

0

200

400

600

800

1000
1200
Function evaluations

Current allocation

1400

1600

1800

2000

Optimal allocation

MSM
Circumcision
FSW
CondomBehavior change
Testing
PMTCT
OVC
ART

F IGURE 4: Comparison of optimization methods for a real-world example of HIV resource allocation. (A) The minimum number of new
infections calculated by each method for the first 2000 function evaluations. As above, the shaded regions show the interquartile ranges over
40 different random seeds. (B) Resource allocations in Swaziland for each HIV program, showing the current budget allocation (left) and theallocation that minimizes new infections (right). MSM = men who have sex with men; FSW = female sex workers; PMTCT = preventionof mother-to-child transmission; OVC = orphans and vulnerable children; ART = antiretroviral therapy.

values; however, more weight would need to be given to more recent values, since ∆E tends to decrease as the algorithm
converges on a solution.
Third, it is possible to relax the assumption of local linearity by varying multiple parameters on a single iteration. However,assuming a separate probability is stored for each parameter combination, this reduces the learning rate; for an n-parameter
problem, modifying a single parameter at each iteration results in a learning rate of 1/2n on average for each parameter; inthe limit where all possible combinations of parameters are considered, the learning rate would be 1/22n . While manageable
for small numbers of parameters (e.g., ≤4), this quickly becomes intractable as the number of parameters grows. Conversely,
if multiple parameters are modified simultaneously, the probabilities of all modified parameters could be updated individually;
this approach is likely to be most effective in very high-dimensional systems where the function E is nearly flat with respect

8

to many of the dimensions, in which case varying parameters one by one may be time-consuming. The superior performanceof simulated annealing compared to BALLSD for small numbers of function evaluations in the 10-parameter Rosenbrock’svalley problem shown in Fig. 2 is likely due to this effect.
Fourth, to circumvent the problem of local minima, the method may be used with either Monte Carlo initialization (Metropolis
and Ulam, 1949) or a Metropolis-Hastings approach (Metropolis et al., 1953). In the former, the procedure described in Sec. 2
is repeated multiple times (typically, 102 − 103 ) with pseudorandom choices of x0 . In the latter, instead of always performing
step 2 of the BALLSD algorithm if the new iteration does not reduce error, step 1 is performed with nonzero acceptance ratio
ρ, where ρ is a function of the change in error; e.g., ρ ∝ Ek± /Ek−1 . Although the parameter set resulting from each iteration
can be kept, as in a standard Metropolis-Hastings algorithm, the value of doing so is limited since the asymptotic distributionof parameter sets is not guaranteed to reach a stationary distribution, due to the adaptive method for choosing which parameters
to vary. Instead, it would suffice to keep two parameter sets, the current one and the best one. As a simpler alternative toimplementing a Metropolis-Hastings approach, rather than always reducing the step size if the new iteration does not reduce
the error, step size could have a nonzero probability of increasing, potentially allowing the algorithm to escape local minima.
Finally, the BALLSD algorithm is only loosely Bayesian. While a more formal Bayesian approach may be desirable in certain
situations, in general it is difficult to determine whether new information should be used to update the existing distribution,or whether the system is in a sufficiently dissimilar part of the parameter space that information from much earlier iterations
is no longer relevant. Nonetheless, for certain problems, additional capacity for adaptation may be beneficial. For example,
the basic implementation of BALLSD described above performs poorly in the classic version of Rosenbrock’s valley (Fig. 1);
for this particular problem, an algorithm that was capable of learning nonlinear parameter combinations would be far more
efficient.
4.2 Conclusions
This paper has presented a simple optimization method inspired by the process of manual parameter fitting that is capable ofoutperforming traditional algorithms for certain classes of problems. The algorithm is most effective for problems with moderate
to large dimensionality (>5 dimensions), which corresponds to the case in which there are enough parameters that different
parameters are likely to have very different overall contributions to the objective function. Indeed, the relative uniformityof parameters in the test functions used here (in terms of both scale and effectiveness) does not necessarily reflect certain
real-world situations in which some – or even most – of the objective function’s parameters may have little influence on its
value. In such situations, BALLSD is especially effective, as it is able to adapt to those parameters (and those scales) that
produce the greatest improvements in the objective function. An example of this is provided in Fig. 4, where BALLSD finds
what appears to be the globally optimal solution more than 10 times faster than any other algorithm.
This study has two main limitations. First, Matlab’s default values of the meta-parameters were used for the other algorithms(except the initial temperature of the simulated annealing, as noted above). Meta-parameter tuning would likely increase theperformance of these algorithms more than it would for BALLSD, since these algorithms are not adaptive – but conversely, anadvantage of BALLSD is that it typically does not require any meta-parameter tuning, so in that sense the comparison is fair.
Second, we chose the four algorithms to compare against BALLSD based on their popularity, as evidenced by their inclusion
in Matlab’s Optimization Toolbox. However, many other optimization algorithms exist, some of which have also been shownto significantly outperform these more traditional methods (e.g., Rios and Sahinidis (2013)).
As noted in Sec. 3.2, BALLSD has already been used successfully in the real-world application of optimizing the allocationof HIV budgets. For this problem, standard optimization methods, including the four compared against BALLSD in this paper,
were found to require an excessively large number of function evaluations to obtain acceptable solutions. This led the authors
to resort to manual parameter fitting until BALLSD was developed. It is our hope that this algorithm may be able to free otherresearchers from similar unpleasantries. Finally, we note that Python and Matlab implementations of BALLSD, as well as thesuite of tests used in this paper, are available for download at http://thekerrlab.com/ballsd.

Acknowledgements
CCK and DPW were supported by World Bank Assignment 1045478. CCK and SDB were supported by DARPA Contract
N66001-10-C-2008. TGS was supported by NIH NCRR 5P20RR016472-12, NIH NIGMS 8P20GM103446-12, NSF EPSCoR
0814251, and NSF HRD 1242067. The authors wish to thank W. W. Lytton, D. Pokrajac, J. Francis, and Z. McGrath for theirhelpful comments on the manuscript.

References
Anderson SJ, Cherutich P, Kilonzo N, Cremin I, Fecht D, Kimanga D, Harper M, Masha RL, Ngongo PB, Maina W et al.
(2014) Maximising the effect of combination HIV prevention through prioritisation of the people and places in greatest
need: a modelling study. The Lancet 384:249–256.

9

Baker JL, Perez-Rosello T, Migliore M, Barrionuevo G, Ascoli GA (2011) A computer model of unitary responses fromassociational/commissural and perforant path synapses in hippocampal CA3 pyramidal cells. Journal of Computational
Neuroscience 31:137–158.
Ben-Ameur W (2004) Computing the initial temperature of simulated annealing. Computational Optimization andApplications 29:369–385.
Bethke AD (1978) Genetic algorithms as function optimizers. University of Michigan.
Brom C, Vyhnnek J, Lukavsk J, Waller D, Kadlec R (2012) A computational model of the allocentric and egocentric spatial
memory by means of virtual agents, or how simple virtual agents can help to build complex computational models. Cognitive
Systems Research 17–18:1–24.
Castillo P, Lozano R, Dzul A (2005) Stabilization of a mini rotorcraft with four rotors. IEEE Control Systems
Magazine 25:45–55.
Charness G, Karni E, Levin D (2007) Individual and group decision making under risk: An experimental study of Bayesian
updating and violations of first-order stochastic dominance. Journal of Risk and Uncertainty 35:129–148.
Goodner J, Tsianos GA, Li Y, Loeb GE (2012) BioSearch: A physiologically plausible learning model for the sensorimotor
system. In Proceedings of the Society for Neuroscience Annual Meeting 275.22/LL11.
Hilbert M, López P (2011) The world’s technological capacity to store, communicate, and compute information. Science 332:60–65.
Kerr CC, Van Albada SJ, Neymotin SA, Chadderdon GL, Robinson P, Lytton WW (2013) Cortical information flow in
Parkinson’s disease: a composite network/field model. Frontiers in Computational Neuroscience 7:1–14.
Kirkpatrick S, Gelatt CD, Vecchi MP (1983) Optimization by simmulated annealing. Science 220:671–680.
Kwon JA, Anderson J, Kerr CC, Thein HH, Zhang L, Iversen J, Dore GJ, Kaldor JM, Law MG, Maher L, Wilson DP (2012)
Estimating the cost-effectiveness of needle-syringe programs in Australia. AIDS 26:2201–2210.
Marquardt DW (1963) An algorithm for least-squares estimation of nonlinear parameters. Journal of the Society for Industrial& Applied Mathematics 11:431–441.
Metropolis N, Rosenbluth AW, Rosenbluth MN, Teller AH, Teller E (1953) Equation of state calculations by fast computing
machines. Journal of Chemical Physics 21:1087–1092.
Metropolis N, Ulam S (1949) The Monte Carlo method. Journal of the American Statistical Association 44:335–341.
Nelder JA, Mead R (1965) A simplex method for function minimization. Computer Journal 7:308–313.
Prinz AA, Billimoria CP, Marder E (2003) Alternative to hand-tuning conductance-based models: construction and analysis of
databases of model neurons. Journal of Neurophysiology 90:3998–4015.
Rios LM, Sahinidis NV (2013) Derivative-free optimization: A review of algorithms and comparison of software implementations. Journal of Global Optimization 56:1247–1293.
Song W, Kerr CC, Lytton WW, Francis JT (2013) Cortical plasticity induced by spike-triggered microstimulation in primate
somatosensory cortex. PLOS ONE 8:e57453.
Steyvers M, Lee MD, Wagenmakers EJ (2009) A Bayesian analysis of human decision-making on bandit problems. Journalof Mathematical Psychology 53:168–179.
Wilby RL (2005) Uncertainty in water resource model parameters used for climate change impact assessment. Hydrological
Processes 19:3201–3219.