1. Welcome to the NSR Physiome Project!
  2. Physiome JSim
  3. The JSim Documentation Set
  4. JSim Optimization Algorithms

JSim Optimization Algorithms

Introduction

There is no perfect optimization algorithm that is best for all problems. Optimization algorithms vary in their approach, efficiency, robustness and applicability to particular problem domains. This document describes the optimization algorithms currently supported by JSim in some detail so users can make intelligent use of them. JSim's currently available optimizers are listed below. Other algorithms are in development.

Prerequisites:

Contents:

  • Overview
  • Simplex
  • GGopt
  • GridSearch
  • Nelder-Mead
  • Nl2sol
  • Sensop
  • Simulated Annealing
  • Genetic Algorithm (version 2.01 and above)
  • Particle Swarm Algorithm (version 2.19 and above)
  • Principle Axis Algorithm (version 2.20 and above)
  • Comments or Questions?
    Give feedback

Overview

Some terminology is useful when discussing the merits of optimization algorithms:

  • Bounded algorithms are those that require upper and lower bounds for each parameter varied. Unbounded algorithms require no such bounds.
  • Parallel algorithms are those that can take advantage of multiple system processors for faster processing. See Using JSim Multiprocessing for further information on JSim multiprocessing.

All of JSim's currently available optimization algorithms share the following control parameters:

  • max # runs: The optimizer will stop if it has run the model this many times.
  • min RMS error: The optimizer will stop if the mean RMS error between reference data and model output is less than this value.

Algorithm-specific control parameters are listed with each algorithm description.

Simplex

JSim's simplex is a bounded, non-linear steepest-descent algorithm. This algorithm does not currently support parallel processing. (description needs work)

Algorithm-specific control parameters:

  • parameter initial step: default=0.01
  • minimum par step: The optimizer will stop if it is considering a parameter step value that vary less than this value.

Reference: Dantzig GB, Orden A, Wolfe P: The generalized simplex method for minimizing a linear form under linear inequality restraints. Pacific J Math. 1955;5(2): 183–195.

GGopt

GGopt is an unbounded non-linear algorithm originally written by Glad and Goldstein. This algorithm does not currently support parallel processing. (description needs work)

Algorithm-specific control parameters:

  • minimum par step: The optimizer will stop if it is considering a parameter step value that vary less than this value.
  • maximum # iterations: default=10
  • minimum gradient: default=1e-6
  • relative error: default=1e-6

Reference: Bassingthwaighte JB, Chan IS, Goldstein AA, et al.: GGOPT: an unconstrained non-linear optimizer. Comput Methods Programs Biomed. 1988;26(3): 275–81. PubMed link

Nelder-Mead

Nelder-Mead is an unbounded, steepest descent algorithm by Nelder and Mead. It is also called non-linear Simplex. This algorithm supports multiprocessing (MP).

During each iteration in a P-parameter optimization, this algorithm performs a P or P+1 parmeter queries (model runs). Several additional single queries are also performed. Ideal MP speedup on an N-processor system on be roughly of order P (if P is a factor of N), or order N (if N is a factor of P).

This algorithm differs from the previously available JSim "Simplex" (above) in that:

  • it is unbounded, while "Simplex" is bounded;
  • it supports MP, while "Simplex" does not;
  • it is a newer implementation of the algorithm (Java vs. Fortran).

Algorithm-specific control parameters:

  • parameter initial step: default=0.01
  • minimum par step: The optimizer will stop if it is considering a parameter step value that vary less than this value.

Reference: Nelder JA, Mead R: A simplex method for function minimization. The Computer Journal. 1965;7(4): 308–313.

GridSearch

GridSearch is a bounded, parallel algorithm. The algorithm operates via progressively restricted search of parameter space on a regularly spaced grid of npoints per dimension. Each iteration, npoints^nparm points are searched for the minimum residual. Each parameter dimension is then restricted to one grid delta around that minimum and the search repeats until stopping criteria are met.

Search bounds in each dimension narrow by a factor of at least 2/(npoints-1) each iteration. Thus, npoints must be at least 4. Each iteration requires up to npoints^nparm residual calculations. Residual calculations are reused when possible, and this reuse is most efficient when npoints is 1 + 2^N for some N. Therefore, npoints defaults to 5, which is the smallest "efficient" value.

This algorithm is very not efficient for very smooth residual functions in high-dimensional space. It works well on noisy functions when low accuraccy situations (e.g. 3 significant digits required). With npoints large, it copes well with multiple local minima. An effective strategy may be to use several interations of GridSearch to estimate a global minimum, and then use a steepest-descent algorithm to fine tune the answer.

The number of points searched during each iteration is typically large compared to the number of available processors. Typical MP speedup on an N-processor system is therefore on the order of N.

Algorithm-specific control parameters:

  • minimum par step: The optimizer will stop if it is considering a parameter step value that vary less than this value.
  • max # iterations: stop after this many iterations, default=10.
  • # grid points: npoints above, default=5.

Reference: Kolda TG, Lewis RM, Torczon V: Optimization by direct search: New perspectives on some classical and modern methods. Siam Review. 2003;45(3): 385–482.

NL2SOL

This version of Nl2sol is derived from NL2SOL, a library of FORTRAN routines which implement an adaptive nonlinear least-squares algorithm. It has been modified to perform all calculations in double precision. It is an unbounded optimizer. This optimizer does not support multi-processing.

Algorithm-specific control parameters:

  • maximum # runs: default=50
  • relative error: default=1e-6

References:

  • Dennis JE, Gay DM, Welsch RE: NL2SOL: An adaptive nonlinear least-squares algorithm. ACM Trans Math Softw. 1981;7(3): 348–368.
  • Dennis JE, Schnabel RB: Numerical methods for unconstrained optimization and nonlinear equation. N. Y.: Prentice-Hall, 1983; 378

 

Sensop

Sensop is a variant of Levenberg-Marquardt algorithm that utilized the maximum parameter sensitivity to determine step size. It is a bounded optimizer, supporting multiprocessing.

Algorithm-specific control parameters:

  • minimum gradient: default=1e-6
  • minimum par step: The optimizer will stop if it is considering a parameter step value that vary less than this value.
  • (JSim v2.16 or greater) Max stat iter: (default =0) When this parameter is non-zero, Sensop will terminate if it makes more than the maximum number of iterations specified by 'Max stat iter' without improving the best parameter estimates. This control parameter is useful for users using JSim w/out user intervention or if the optimizer has found a good fit but none of the other stopping criteria have been met.
    For example, if you set 'maximum # runs' to 5000 and Sensop found a best fit at run #160 but continued to try 4840 times to improve the fit. But if 'Max stat iter' is set to 100 then it will stop when all of the parameter estimates have not changed for 100 consecutive optimizer iterations. The correlation between optimizer iterations and model evaluations (runs) is dependent on the number of parameters to be optimized.

The JSim java implimented Sensop algorithm can be found here (Github).

Reference: Chan IS, Goldstein AA, Bassingthwaighte JB: SENSOP: a derivative-free solver for non-linear least squares with sensitivity scaling. Ann Biomed Eng. 1993;21(6): 621–31. PubMed link.

Simulated Annealing

Simulated annealing is an algorithm inspired by the annealing process in metullurgy. As the problem "cools" from its initial "temperature", random fluctuations of model parameters are reduced. JSim implements a bounded version of this algorithm that supports multiprocessing.

Algorithm-specific control parameters:

  • initial temperature:default=100. This parameter has no physical meaning (e.g. Kelvin), but must be positive.

The JSim java implimented Simulated Annealing algorithm can be found here (Github).

Reference: Kirkpatrick S, Gelatt CD Jr, Vecchi MP: Optimization by simulated annealing. Science. 1983; 220(4598): 671–680

Genetic Algorithms

This algorithm is available only in JSim version 2.01 and above.

Genetic algorithms are a family of algorithms that generate a population of candidate solutions and then select the best solutions in each iteration. A new population of solutions is created during each iteration. There are different ways of specifying how a new population is generated from the existing population. The error calculation is used to score and rank the candidate solution. The "fit" individuals in the existing population are selected using one of three methods: (1) roulette-wheel, (2) tournament selection, or (3) elitism. In the roulette-wheel method, the probability of a solution being selected is inversely proportional to the error. In the tournament selection, two random solutions are selected and the one with the lower error is placed in the new population. In elitism, all the solutions with the lowest errors (with a cutoff fraction) are selected. New solutions are selected by "mutating" and "crossing over" existing solutions.

Algorithm-specific control parameters:

  • Population:default= 25. Number of trials in each generation. If max # runs < Population then Population is defaulted to the value of max # runs and only the parent generation is calculated. Suggested values are max # runs = 400, Population=100 for four generations.
  • Mutation rate:default = 0.1. The probability of any single parameter getting changed.
  • Crossover rate:default = 0.5. The probability of creating a new solution by combining existing solutions.
  • Select Method: default= 1. Acceptable values are (1) roulette-wheel, (2) tournament selection, and (3) elitism. See above description identifying these terms.
  • Mutation step:default=0.05. The amount by which a mutation changes a parameter, specified as a fraction of the range.
  • Elite cutoff:default=0.5. Fraction of the population considered fit for the elitism selection method (Select Method=3.)

 

The JSim java implimented Genetic algorithm can be found here (Github). The JSim genetic algorithm supports multiprocessors (parallel processing).

Reference: Holland JH: Adaptation in natural and artificial Systems: an introductory analysis with applications to biology, control, and artificial intelligence. MIT Press. 1992; 183.

Particle Swarm Algorithm (PSO)

This algorithm is available only in JSim version 2.19 and above.

The PSO algorithm works by having swarm candidate solutions (called particles) that move around in the search-space according to a few simple rules. The movements of the particles are guided by their own best known position in the search-space as well as the entire swarm's best known position. When improved positions are being discovered these will then come to guide the movements of the swarm. The process is repeated and by doing so it is hoped, but not guaranteed, to converge to a minimum value. This optimization can be useful when fitting a large number of parameters (>10) when estimates of their values are unknown or have a large range of values. The PSO algorithm supports multiprocessors (parallel processing) in JSim release v 2.20 and later.

Algorithm-specific control parameters:

  • # of particles: default=25. Number of particles in the swarm. May need larger values (50 or greater) when optimizing 20, 30 or more parameters, or the range of possible values for a parameter is large.
  • velocity coeff: default = 0.25 (range: 0 - 1.0). Velocity coefficient used to calculate initial particle velocity. Calculation for initial indiviual particle velocity (parameter,n, and particle,i,) is: 
    v_particle[n][i] = vel_coeff*(xmax[n]-xmin[n])*(random number[0 < # <1.0])
    Where xmax[n], xmin[n] are a parameter's maximum and minimum values as specified by the user.
  • Min inertia: default=0.3 (range: 0 - 1.0). Coefficient that reflects a particle's minimum resistance to changing velocity.
  • Max inertia: default=1.0 (range: min inertia - 1.0). Coefficient that reflects a particle's maximum resistance to changing velocity. The Inertia calculation for each particle is: 
     w = w_min +it*dw,
    where w is particle inertia, w_min is minimum inertia, 'it' is iteration number (max_iter is the maximum amount of iterations specified), and dw is change in inertia with each iteration and is given by: 
     dw = (w_max-w_min)/max_iter
  • Cog learn: default=1.05 (range:0 - 3.0). Determines how much of the particle's personal best fit influences its trajectory.
  • Soc learn: default=1.05 (range:0 - 3.0). Determines how much of the swarm's best fit influences the individual particle's trajectory.
    The updated position of each particle (parameter, n, and particle, i,) is given by: 
    particle_x[n][i] = particle_x[n][i] + v_particle[n][i]
    The updated velocity of each particle (parameter,n, and particle,i,) is given by: 
    v_particle[n][i] = w*v_particle[n][i] + 
c1*(random number[0 < # <1.0]) *(x_pbest[n][i] - particle_x[n][i]) +
c2*(random number[0 < # <1.0]) * (x_gbest[n] - particle_x[n][i])
    where c1 is cog learn factor, c2 is soc learn factor, x_pbest is particle's personal best position, x_gbest is the swarm's best position so far.

 

Further work/issues:

  • Add control of randomization: Allow user to specify the seed for random number generator so that user can repeat optimization. Currently, each optimization run is a little different due to seed being different with each run.

 

The JSim java implimented Particle Swarm algorithm can be found here (Github).

Reference: James Kennedy and Russell C. Eberhart. Particle swarm optimization. In Proceedings of the 1995 IEEE International Conference on Neural Networks, pages 1942–1948, Piscataway, New Jersey, 1995. IEEE Service Center. IEEE paper.

Principle Axis Algorithm (PRAXIS)

This algorithm is available only in JSim version 2.20 and above.

PRAXIS seeks an M-dimensional point X (parameters to fit) which minimizes a given fitness, or cost, function. PRAXIS minimizes along directions that are conjugates of each other to locate a minima and is a derivative free optimizer. The code is a refinement of Powell's method of conjugate search directions and the fitness function need not be smoothly differentiable.

Algorithm-specific control parameters:

  • Bounded?: default = true (checked). This algorithm is originally unbounded but is modified to accept bounds. The decision to bound the parameter space depends on the type of equations in the model. Often ODEs and PDEs can be very stiff causing the solvers to stop. Using boundaries allows the user to keep the parameter values inside realistic limits. To use the unbounded optimizer just uncheck box. Note: the the JSim optimizer configuration page will still check that the initial parameter value is within the listed min/max values so you will have to adjust them to some large value.
  • Max dist to min: default = 1 (range 0.01 - 1.0). Relative distance from initial parameter value to actual value at minimum. This value is applied to the average range of starting values for all parameters. This value does not affect convergence but can increase or decrease the amount of iterations needed to reach the minimum. This setting is less useful when parameters' starting values are very different (ie a=0.01, b=2500).
  • Iter no improve: default = 1 (range: 1-4). Number of iterations with no improvement in minimum before it stops. '1' is adequate, '4' is conservative.
  • Fit tolerence: default = 1E-10 (range: number greater than zero). Max Tolerence (T0) between actual minimum and current minimum. Based on the precision of the computer's arithmetic. From comments in code: 
    PRAXIS attempts to return praxis = f(x) such that 
if X0 is the true local minimum near X, 
then norm ( x - x0 ) < T0 + sqrt ( EPSILON ( X ) ) * norm ( X ),
    where EPSILON ( X ) is the machine precision for X.

 

 

The JSim java implimented PRAXIS algorithm can be found here (Github).

Reference:

  • Brent, Richard P., Algorithms for minimization without derivatives, chapter 7, 1973, Prentice-Hall, NJ, USA
  • Code reference: C++ PRAXIS code

 

Comments or Questions?

Give feedback

 

Model development and archiving support at https://www.imagwiki.nibib.nih.gov/physiome provided by the following grants: NIH U01HL122199 Analyzing the Cardiac Power Grid, 09/15/2015 - 05/31/2020, NIH/NIBIB BE08407 Software Integration, JSim and SBW 6/1/09-5/31/13; NIH/NHLBI T15 HL88516-01 Modeling for Heart, Lung and Blood: From Cell to Organ, 4/1/07-3/31/11; NSF BES-0506477 Adaptive Multi-Scale Model Simulation, 8/15/05-7/31/08; NIH/NHLBI R01 HL073598 Core 3: 3D Imaging and Computer Modeling of the Respiratory Tract, 9/1/04-8/31/09; as well as prior support from NIH/NCRR P41 RR01243 Simulation Resource in Circulatory Mass Transport and Exchange, 12/1/1980-11/30/01 and NIH/NIBIB R01 EB001973 JSim: A Simulation Analysis Platform, 3/1/02-2/28/07.