paretosearch Algorithm Overview
paretosearch algorithm uses pattern search on a set of
points to search iteratively for nondominated points. See Multiobjective Terminology. The pattern search
satisfies all bounds and linear constraints at each iteration.
Theoretically, the algorithm converges to points near the true Pareto front. For a discussion and proof of convergence, see Custòdio et al. , whose proof applies to problems with Lipschitz continuous objectives and constraints.
paretosearch uses a number of intermediate quantities and
tolerances in its algorithm.
The rank of a point has an iterative definition.
Hypervolume of the set of points p in objective function space that satisfy the inequality, for every index j,
fi(j) < pi < Mi,
where fi(j) is the ith component of the jth objective function value in the Pareto set, and Mi is an upper bound for the ith component for all points in the Pareto set. In this figure, M is called the Reference Point. The shades of gray in the figure denote portions of the volume that some calculation algorithms use as part of an inclusion-exclusion calculation.
For details, see Fleischer .
Volume change is one factor in stopping the algorithm. For details, see Stopping Conditions.
Distance is a measure of the closeness of an individual to
its nearest neighbors. The
The algorithm sets the distance of
individuals at the extreme positions to
The algorithm sorts each dimension separately, so the term neighbors means neighbors in each dimension.
Individuals of the same rank with a higher distance have a higher chance of selection (higher distance is better).
Distance is one factor in the calculation of the spread, which is part of a stopping criterion. For details, see Stopping Conditions.
Spread is a measure of the movement of the Pareto set. To
calculate the spread, the
The spread is small when the extreme objective function values do not change much between iterations (that is, μ is small) and when the points on the Pareto front are spread evenly (that is, σ is small).
|Stopping condition for the search.
Minimum fraction of locations to poll during an iteration.
This option does not apply when the
To create the initial set of points,
options.ParetoSetSize points from a quasirandom sample based
on the problem bounds, by default. For details, see Bratley and Fox . When the problem
has over 500 dimensions,
paretosearch uses Latin
hypercube sampling to generate the initial points.
If a component has no bounds,
paretosearch uses an artificial
lower bound of
-10 and an artificial upper bound of
If a component has only one bound,
paretosearch uses that
bound as an endpoint of an interval of width
20 + 2*abs(bound).
For example, if there is no upper bound for a component and there is a lower bound
paretosearch uses an interval width of 20 + 2*15 = 55,
so uses an artificial upper bound of 15 + 55 = 70.
If you pass some initial points in
paretosearch uses those points as the initial points.
paretosearch generates more points, if necessary, to obtain
options.ParetoSetSize initial points.
paretosearch then checks the initial points to ensure that
they are feasible with respect to the bounds and linear constraints. If necessary,
paretosearch projects the initial points onto the linear
subspace of linearly feasible points by solving a linear programming problem. This
process can cause some points to coincide, in which case
paretosearch removes any duplicate points.
paretosearch does not alter initial points for artificial
bounds, only for specified bounds and linear constraints.
After moving the points to satisfy linear constraints, if necessary,
paretosearch checks whether the points satisfy the
paretosearch gives a penalty value of
Inf to any point that does not satisfy all nonlinear
paretosearch calculates any missing objective
function values of the remaining feasible points.
paretosearch does not support nonlinear equality
ceq(x) = 0.
Create Archive and Incumbents
paretosearch maintains two sets of points:
archive— A structure that contains nondominated points associated with a mesh size below
options.MeshToleranceand satisfying all constraints to within
archivestructure contains no more than
2*options.ParetoSetSizepoints and is initially empty. Each point in
archivecontains an associated mesh size, which is the mesh size at which the point was generated.
iterates— A structure containing nondominated points and possibly some dominated points associated with larger mesh sizes or infeasibility. Each point in
iteratescontains an associated mesh size.
iteratescontains no more than
Poll to Find Better Points
paretosearch polls points from
with the polled points inheriting the associated mesh size from the point in
uses a poll that maintains feasibility with respect to bounds and all linear
If the problem has nonlinear constraints,
computes the feasibility of each poll point.
the score of infeasible points separately from the score of feasible points. The
score of a feasible point is the vector of objective function values of that point.
The score of an infeasible point is the sum of the nonlinear infeasibilities.
paretosearch polls at least
MinPollFraction*(number of points in pattern) locations for
each point in
iterates. If the polled points give at least one
nondominated point with respect to the incumbent (original) point, the poll is
considered a success. Otherwise,
paretosearch continues to poll
until it either finds a nondominated point or runs out of points in the pattern. If
paretosearch runs out of points and does not produce a
paretosearch declares the poll unsuccessful
and halves the mesh size.
If the poll finds nondominated points,
the poll in the successful directions repeatedly, doubling the mesh size each time,
until the extension produces a dominated point. During this extension, if the mesh
options.MaxMeshSize (default value:
Inf), the poll stops. If the objective function values
paretosearch declares the
problem unbounded and stops.
After polling all the points in
iterates, the algorithm
examines the new points together with the points in the
computes the rank, or Pareto front number, of each point and then does the
Mark for removal all points that do not have rank 1 in
Mark new rank
1points for insertion into
Mark feasible points in
iterateswhose associated mesh size is less than
options.MeshTolerancefor transfer to
Mark dominated points in
iteratesfor removal only if they prevent new nondominated points from being added to
paretosearch then computes the volume and distance measures
for each point. If
archive will overflow as a result of marked
points being included, then the points with the largest volume occupy
archive, and the others leave. Similarly, the new points
marked for addition to
in order of their volumes.
iterates is full and has no dominated points, then
paretosearch adds no points to
and declares the iteration to be unsuccessful.
multiplies the mesh sizes in
iterates by 1/2.
For three or fewer objective functions,
volume and spread as stopping measures. For four or more objectives,
paretosearch uses distance and spread as stopping measures.
In the remainder of this discussion, the two measures that
paretosearch uses are denoted the
The algorithm maintains vectors of the last eight values of the applicable
measures. After eight iterations, the algorithm checks the values of the two
applicable measures at the beginning of each iteration, where
spreadConverged = abs(spread(end - 1) - spread(end)) <= tol*max(1,spread(end - 1));
volumeConverged = abs(volume(end - 1) - volume(end)) <= tol*max(1,volume(end - 1));
distanceConverged = abs(distance(end - 1) - distance(end)) <= tol*max(1,distance(end - 1));
If either applicable test is
true, the algorithm stops.
Otherwise, the algorithm computes the max of squared terms of the Fourier transforms
of the applicable measures minus the first term. The algorithm then compares the
maxima to their deleted terms (the DC components of the transforms). If either
deleted term is larger than
100*tol*(max of all other terms),
then the algorithm stops. This test essentially determines that the sequence of
measures is not fluctuating, and therefore has converged.
Additionally, a plot function or output function can stop the algorithm, or the algorithm can stop because it exceeds a time limit or function evaluation limit.
The algorithm returns the points on the Pareto front as follows.
paretosearchcombines the points in
iteratesinto one set.
When there are three or fewer objective functions,
paretosearchreturns the points from the largest volume to the smallest, up to at most
When there are four or more objective functions,
paretosearchreturns the points from the largest distance to the smallest, up to at most
Modifications for Parallel Computation and Vectorized Function Evaluation
paretosearch computes objective function values in
parallel or in a vectorized fashion (
true), there are some changes to the algorithm.
MinPollFractionoption and evaluates all poll points in the pattern.
When computing in parallel,
paretosearchsequentially examines each point in
iteratesand performs a parallel poll from each point. After returning
MinPollFractionfraction of the poll points,
paretosearchdetermines if any poll points dominate the base point. If so, the poll is deemed successful, and any other parallel evaluations halt. If not, polling continues until a dominating point appears or the poll is done.
paretosearchperforms objective function evaluations either on workers or in a vectorized fashion, but not both. If you set both
paretosearchcalculates objective function values in parallel on workers, but not in a vectorized fashion. In this case,
MinPollFractionoption and evaluates all poll points in the pattern.
The fastest way to run
paretosearch depends on several
If objective function evaluations are slow, then it is usually fastest to use parallel computing. The overhead in parallel computing can be substantial when objective function evaluations are fast, but when they are slow, it is usually best to use more computing power.
Parallel computing requires a Parallel Computing Toolbox™ license.
If objective function evaluations are not very time consuming, then it is usually fastest to use vectorized evaluation. However, this is not always the case, because vectorized computations evaluate an entire pattern, whereas serial evaluations can take just a small fraction of a pattern. In high dimensions especially, this reduction in evaluations can cause serial evaluation to be faster for some problems.
To use vectorized computing, your objective function must accept a matrix with an arbitrary number of rows. Each row represents one point to evaluate. The objective function must return a matrix of objective function values with the same number of rows as it accepts, with one column for each objective function. For a single-objective discussion, see Vectorize the Fitness Function (
ga) or Vectorized Objective Function (
 Custòdio, A. L., J. F. A. Madeira, A. I. F. Vaz, and L. N. Vicente. Direct Multisearch for Multiobjective Optimization. SIAM J. Optim., 21(3), 2011, pp. 1109–1140. Preprint available at https://www.researchgate.net/publication/220133323_Direct_Multisearch_for_Multiobjective_Optimization.
 Bratley, P., and B. L. Fox. Algorithm 659: Implementing Sobol’s quasirandom sequence generator. ACM Trans. Math. Software 14, 1988, pp. 88–100.
 Fleischer, M. The Measure of Pareto Optima: Applications to Multi-Objective Metaheuristics. In "Proceedings of the Second International Conference on Evolutionary Multi-Criterion Optimization—EMO" April 2003 in Faro, Portugal. Published by Springer-Verlag in the Lecture Notes in Computer Science series, Vol. 2632, pp. 519–533. Preprint available at https://apps.dtic.mil/dtic/tr/fulltext/u2/a441037.pdf.