Algorithms¶
Generic optimization algorithms.
biogeme.algorithms module¶
Optimization algorithms.
 author
Michel Bierlaire
 date
Sun Apr 5 16:48:54 2020

biogeme.algorithms.
bfgs
(H, d, y)[source]¶ Update the BFGS matrix. Formula (13.12) of Bierlaire (2015) where the method proposed by Powell (1977) is applied
 Parameters
H (numpy.array (2D)) – current approximation of the inverse of the Hessian
d (numpy.array (1D)) – difference between two consecutive iterates.
y (numpy.array (1D)) – difference between two consecutive gradients.
 Returns
updated approximation of the inverse of the Hessian.
 Return type
numpy.array (2D)

biogeme.algorithms.
bfgsLineSearch
(fct, x0, initBfgs=None, eps=6.06273418136464e06, maxiter=1000)[source]¶ BFGS method with inexact line search (Wolfe conditions)
 Parameters
fct (optimization.functionToMinimize) – object to calculate the objective function and its derivatives.
x0 (numpy.array) – starting point
initBfgs (numpy.array) – matrix used to initialize BFGS. If None, the identity matrix is used. Default: None.
eps (float) – the algorithm stops when this precision is reached. Default: \(\varepsilon^{\frac{1}{3}}\)
maxiter (int) – the algorithm stops if this number of iterations is reached. Default: 1000
 Returns
tuple x, messages, where
x is the solution found,
messages is a dictionary reporting various aspects related to the run of the algorithm.
 Return type
numpy.array, dict(str:object)
 Raises
biogeme.exceptions.biogemeError – if the dimensions of the matrix initBfgs do not match the length of x0.

biogeme.algorithms.
bfgsTrustRegion
(fct, x0, initBfgs=None, delta0=1.0, eps=6.06273418136464e06, dl=False, maxiter=1000, eta1=0.01, eta2=0.9)[source]¶ BFGS method with trust region
 Parameters
fct (optimization.functionToMinimize) – object to calculate the objective function and its derivatives.
x0 (numpy.array) – starting point
initBfgs (numpy.array) – matrix used to initialize BFGS. If None, the identity matrix is used. Default: None.
delta0 (float) – initial radius of the trust region. Default: 100.
eps (float) – the algorithm stops when this precision is reached. Default: \(\varepsilon^{\frac{1}{3}}\)
dl (bool) – If True, the Dogleg method is used to solve the trut region subproblem. If False, the truncated conjugate gradient is used. Default: False.
maxiter (int) – the algorithm stops if this number of iterations is reached. Default: 1000.
eta1 (float) – threshold for failed iterations. Default: 0.01.
eta2 (float) – threshold for very successful iterations. Default 0.9.
 Returns
tuple x, messages, where
x is the solution found,
messages is a dictionary reporting various aspects related to the run of the algorithm.
 Return type
numpy.array, dict(str:object)
 Raises
biogeme.exceptions.biogemeError – if the dimensions of the matrix initBfgs do not match the length of x0.

class
biogeme.algorithms.
bioBounds
(b)[source]¶ Bases:
object
This class is designed for the management of simple bound constraints

__init__
(b)[source]¶  Parameters
b (list(tuple)) – list of tuples (ell,u) containing the lower and upper bounds for each free parameter.
 Raises
biogeme.exceptions.biogemeError – if the bounds are incompatible

activity
(x, epsilon=2.220446049250313e16)[source]¶ Determines the activity status of each variable.
 Parameters
x (numpy.array) – point for which the activity must be determined.
epsilon (float) – a bound is considered active if the distance to it is less rhan epsilon.
 Returns
a vector, same length as x, where each entry reports the activity of the corresponding variable:
0 if no bound is active
1 if the lower bound is active
1 if the upper bound is active
 Return type
numpy.array
 Raises
biogeme.exceptions.biogemeError – if the vector x is not feasible
biogeme.exceptions.biogemeError – if the dimensions of x and bounds do not match.

bounds
¶ list of tuples (ell,u) containing the lower and upper bounds for each free parameter.

breakpoints
(x, d)[source]¶ Projects the direction d, starting from x, on the intersection of the bound constraints
 Parameters
x (numpy.array) – current point
d (numpy.array) – search direction
 Returns
list of tuple (index, value), where index is the index of the variable, and value the value of the corresponding breakpoint.
 Return type
list(tuple(int,float))
 Raises
biogeme.exceptions.biogemeError – if the dimensions are inconsistent
biogeme.exceptions.biogemeError – if x is infeasible

feasible
(x)[source]¶ Check if point verifies the bound constraints
 Parameters
x (numpy.array) – point to project
 Returns
True if x is feasible, False otherwise.
 Return type
bool
 Raises
biogeme.exceptions.biogemeError – if the dimensions are inconsistent

generalizedCauchyPoint
(xk, gk, H, direction)[source]¶ Implementation of Step 2 of the Specific Algorithm by Conn et al. (1988).
The quadratic model is defined as
\[m(x) = f(x_k) + (x  x_k)^T g_k + \frac{1}{2} (xx_k)^T H (xx_k).\] Parameters
xk (numpy.array. Dimension n.) – current point
gk (numpy.array. Dimension n.) – vector g involved in the quadratic model definition.
H (numpy.array. Dimension n x n.) – matrix H involved in the quadratic model definition.
direction (numpy.array. Dimension n.) – direction along which the GCP is searched.
 Returns
generalized Cauchy point based on inexact line search.
 Return type
numpy.array. Dimension n.
 Raises
biogeme.exceptions.biogemeError – if the dimensions are inconsistent
biogeme.exceptions.biogemeError – if xk is infeasible

intersect
(otherBounds)[source]¶ Create a bounds object representing the intersection of two regions.
 Parameters
otherBounds (class bioBounds) – other bound object that must be intersected.
 Returns
bound object, intersection of the two.
 Return type
class bioBounds
 Raises
biogeme.exceptions.biogemeError – if the dimensions are inconsistent

intersectionWithTrustRegion
(x, delta)[source]¶ Create a bioBounds object representing the intersection between the feasible domain and the trust region.
 Parameters
x (numpy.array) – center of the trust region
delta (float) – radius of the tust region (infinity norm)
 Returns
intersection between the feasible region and the trus region
 Return type
class bioBounds
 Raises
biogeme.exceptions.biogemeError – if the dimensions are inconsistent

lowerBounds
¶ List of lower bounds

maximumStep
(x, d)[source]¶ Calculates the maximum step that can be performed along a direction while staying feasible.
 Parameters
x (numpy.array) – reference point
d (numpy.array) – direction
 Returns
the largest alpha such that x + alpha * d is feasible and the list of indices achieving this value.
 Return type
float, int
 Raises
biogeme.exceptions.biogemeError – if the point is infeasible

n
¶ number of optimization variables

project
(x)[source]¶ Project a point onto the feasible domain defined by the bounds.
 Parameters
x (numpy.array) – point to project
 Returns
projected point
 Return type
numpy.array
 Raises
biogeme.exceptions.biogemeError – if the dimensions are inconsistent

subspace
(selectedVariables)[source]¶ Generate a bioBounds object for selected variables
 Parameters
selectedVariables (numpy.array(bool)) – boolean vector. If an entry is True, the corresponding variables is considered.
 Returns
bound object
 Return type
class bioBounds
 Raises
biogeme.exceptions.biogemeError – if the dimensions are inconsistent

upperBounds
¶ List of upper bounds


biogeme.algorithms.
cauchyNewtonDogleg
(g, H)[source]¶ Calculate the Cauchy, the Newton and the dogleg points.
The Cauchy point is defined as
\[d_c =  \frac{\nabla f(x)^T \nabla f(x)} {\nabla f(x)^T \nabla^2 f(x) \nabla f(x)} \nabla f(x)\]The Newton point \(d_n\) verifies Newton equation:
\[H_s d_n =  \nabla f(x)\]where \(H_s\) is a positive definite matrix generated with the method by Schnabel and Eskow (1999).
The Dogleg point is
\[d_d = \eta d_n\]where
\[\eta = 0.2 + 0.8 \frac{\alpha^2}{\beta \nabla f(x)^T d_n}\]and \(\alpha= \nabla f(x)^T \nabla f(x)\), \(\beta=\nabla f(x)^T \nabla^2 f(x)\nabla f(x).\)
 Parameters
g (numpy.array) – gradient \(\nabla f(x)\)
H (numpy.array) – hessian \(\nabla^2 f(x)\)
 Returns
tuple with Cauchy point, Newton point, Dogleg point
 Return type
numpy.array, numpy.array, numpy.array
 Raises
biogeme.exceptions.biogemeError – if the quadratic model is not convex.

biogeme.algorithms.
dogleg
(g, H, delta)[source]¶ Find an approximation of the trust region subproblem using the dogleg method
 Parameters
g (numpy.array) – gradient of the quadratic model.
H (numpy.array) – hessian of the quadratic model.
delta (float) – radius of the trust region.
 Returns
d, diagnostic where
d is an approximate solution of the trust region subproblem
diagnostic is the nature of the solution:
2 if negative curvature along Newton direction
1 if negative curvature along Cauchy direction (i.e. along the gradient)
1 if partial Cauchy step
2 if Newton step
3 if partial Newton step
4 if Dogleg
 Return type
numpy.array, int

class
biogeme.algorithms.
functionToMinimize
[source]¶ Bases:
object
This is an abstract class. The actual function to minimize must be implemented in a concrete class deriving from this one.

abstract
f
(batch=None)[source]¶ Calculate the value of the function
 Parameters
batch (float) – for data driven functions (such as a log likelikood function), it is possible to approximate the value of the function using a sample of the data called a batch. This argument is a value between 0 and 1 representing the percentage of the data that should be used for thre random batch. If None, the full data set is used. Default: None pass
 Returns
value of the function
 Return type
float

abstract
f_g
(batch=None)[source]¶ Calculate the value of the function and the gradient
 Parameters
batch (float) – for data driven functions (such as a log likelikood function), it is possible to approximate the value of the function using a sample of the data called a batch. This argument is a value between 0 and 1 representing the percentage of the data that should be used for the random batch. If None, the full data set is used. Default: None pass
 Returns
value of the function and the gradient
 Return type
tuple float, numpy.array

abstract
f_g_bhhh
(batch=None)[source]¶ Calculate the value of the function, the gradient and the BHHH matrix
 Parameters
batch (float) – for data driven functions (such as a log likelikood function), it is possible to approximate the value of the function using a sample of the data called a batch. This argument is a value between 0 and 1 representing the percentage of the data that should be used for the random batch. If None, the full data set is used. Default: None pass
 Returns
value of the function, the gradient and the BHHH
 Return type
tuple float, numpy.array, numpy.array

abstract
f_g_h
(batch=None)[source]¶ Calculate the value of the function, the gradient and the Hessian
 Parameters
batch (float) – for data driven functions (such as a log likelikood function), it is possible to approximate the value of the function using a sample of the data called a batch. This argument is a value between 0 and 1 representing the percentage of the data that should be used for the random batch. If None, the full data set is used. Default: None pass
 Returns
value of the function, the gradient and the Hessian
 Return type
tuple float, numpy.array, numpy.array

abstract

biogeme.algorithms.
inverseBfgs
(Hinv, d, y)[source]¶ Update the inverse BFGS matrix. Formula (13.13) of Bierlaire (2015)
 Parameters
Hinv (numpy.array (2D)) – current approximation of the inverse of the Hessian
d (numpy.array (1D)) – difference between two consecutive iterates.
y (numpy.array (1D)) – difference between two consecutive gradients.
 Returns
updated approximation of the inverse of the Hessian.
 Return type
numpy.array (2D)

biogeme.algorithms.
lineSearch
(fct, x, f, g, d, alpha0=1.0, beta1=0.0001, beta2=0.99, lbd=2.0)[source]¶ Calculate a step along a direction that satisfies both Wolfe conditions
 Parameters
fct (optimization.functionToMinimize) – object to calculate the objective function and its derivatives.
x (numpy.array) – current iterate.
f (float) – value of the objective function at the current iterate.
g (numpy.array) – value of the gradient at the current iterate.
d (numpy.array) – descent direction.
alpha0 (float) – first step to test.
beta1 (float) – parameter of the first Wolfe condition.
beta2 (float) – parameter of the second Wolfe condition.
lbd (float) – expansion factor for a short step.
 Returns
a step verifing both Wolfe conditions
 Return type
float
 Raises
biogeme.exceptions.biogemeError – if
lbd
\(\leq\) 1biogeme.exceptions.biogemeError – if
alpha0
\(\leq\) 0biogeme.exceptions.biogemeError – if
beta1
\(\geq\) beta2biogeme.exceptions.biogemeError – if
d
is not a descent direction

biogeme.algorithms.
newtonLineSearch
(fct, x0, eps=6.06273418136464e06, maxiter=100)[source]¶ Newton method with inexact line search (Wolfe conditions)
 Parameters
fct (optimization.functionToMinimize) – object to calculate the objective function and its derivatives.
x0 (numpy.array) – starting point
eps (float) – the algorithm stops when this precision is reached. Default: \(\varepsilon^{\frac{1}{3}}\)
maxiter (int) – the algorithm stops if this number of iterations is reached. Defaut: 100
 Returns
x, messages
x is the solution generated by the algorithm,
messages is a dictionary describing information about the lagorithm
 Return type
numpay.array, dict(str:object)

biogeme.algorithms.
newtonTrustRegion
(fct, x0, delta0=1.0, eps=6.06273418136464e06, dl=False, maxiter=1000, eta1=0.01, eta2=0.9)[source]¶ Newton method with trust region
 Parameters
fct (optimization.functionToMinimize) – object to calculate the objective function and its derivatives.
x0 (numpy.array) – starting point
delta0 (float) – initial radius of the trust region. Default: 100.
eps (float) – the algorithm stops when this precision is reached. Default: \(\varepsilon^{\frac{1}{3}}\)
dl (bool) – If True, the Dogleg method is used to solve the trut region subproblem. If False, the truncated conjugate gradient is used. Default: False.
maxiter (int) – the algorithm stops if this number of iterations is reached. Default: 1000.
eta1 (float) – threshold for failed iterations. Default: 0.01.
eta2 (float) – threshold for very successful iterations. Default 0.9.
 Returns
tuple x, messages, where
x is the solution found,
messages is a dictionary reporting various aspects related to the run of the algorithm.
 Return type
numpy.array, dict(str:object)

biogeme.algorithms.
relativeChange
(x, xpred, typx)[source]¶ Calculates the relative change.
It is typically used as stopping criterion.
 Parameters
x (numpy.array) – current iterate.
xpred (numpy.array) – previous iterate.
typx (numpy.array) – typical value for x.
 Returns
relative change:
\[\max_i \frac{(x_k)_i  (x_{k1})_i}{\max((x_k)_i, \text{typx}_i)}\] Return type
float

biogeme.algorithms.
relativeGradient
(x, f, g, typx, typf)[source]¶ Calculates the relative gradients.
It is typically used as stopping criterion.
 Parameters
x (numpy.array) – current iterate.
f (float) – value of f(x)
g (numpy.array) – \(\nabla f(x)\), gradient of f at x
typx (numpy.array) – typical value for x.
typf (float) – typical value for f.
 Returns
relative gradient
\[\max_{i=1,\ldots,n}\frac{(\nabla f(x))_i \max(x_i,\text{typx}_i)} {\max(f(x), \text{typf})}\] Return type
float

biogeme.algorithms.
schnabelEskow
(A, tau=6.06273418136464e06, taubar=3.6756745753887175e11, mu=0.1)[source]¶ Modified Cholesky factorization by Schnabel and Eskow (1999).
If the matrix is ‘safely’ positive definite, the output is the classical Cholesky factor. If not, the diagonal elements are inflated in order to make it positive definite. The factor \(L\) is such that \(A + E = PLL^TP^T\), where \(E\) is a diagonal matrix contaninig the terms added to the diagonal, \(P\) is a permutation matrix, and \(L\) is w lower triangular matrix.
 Parameters
A (numpy.array) – matrix to factorize. Must be square and symmetric.
tau (float) – tolerance factor. Default: \(\varepsilon^{\frac{1}{3}}\). See Schnabel and Eskow (1999)
taubar (float) – tolerance factor. Default: \(\varepsilon^{\frac{2}{3}}\). See Schnabel and Eskow (1999)
mu (float) – tolerance factor. Default: 0.1. See Schnabel and Eskow (1999)
 Returns
tuple \(L\), \(E\), \(P\), where \(A + E = PLL^TP^T\).
 Return type
numpy.array, numpy.array, numpy.array
 Raises
biogeme.exceptions.biogemeError – if the matrix A is not square.
biogeme.exceptions.biogemeError – if the matrix A is not symmetric.

biogeme.algorithms.
simpleBoundsNewtonAlgorithm
(fct, bounds, x0, proportionTrueHessian=1.0, infeasibleConjugateGradient=False, delta0=1.0, tol=6.06273418136464e06, steptol=1e05, cgtol=6.06273418136464e06, maxiter=1000, eta1=0.01, eta2=0.9, enlargingFactor=10)[source]¶ Trust region algorithm for problems with simple bounds
 Parameters
fct (optimization.functionToMinimize) – object to calculate the objective function and its derivatives.
bounds (class bounds) – bounds on the variables
x0 (numpy.array) – starting point
proportionTrueHessian (float) – proportion of the iterations where the true hessian is calculated. When not, the BFGS update is used. If 1.0, it is used for all iterations. If 0.0, it is not used at all.
infeasibleConjugateGradient (bool) – if True, the conjugate gradient algorithm may generate until termination. The result will then be projected on the feasible domain. If False, the algorithm stops as soon as an infeasible iterate is generated. Default: False.
delta0 (float) – initial radius of the trust region. Default: 100.
tol (float) – the algorithm stops when this precision is reached. Default: \(\varepsilon^{\frac{1}{3}}\)
steptol (float) – the algorithm stops when the relative change in x is below this threshold. Basically, if p significant digits of x are needed, steptol should be set to 1.0ep. Default: \(10^{5}\)
cgtol (float) – the conjugate gradient algorithm stops when this precision is reached. Default: \(\varepsilon^{\frac{1}{3}}\)
maxiter (int) – the algorithm stops if this number of iterations is reached. Default: 1000.
eta1 (float) – threshold for failed iterations. Default: 0.01.
eta2 (float) – threshold for very successful iterations. Default 0.9.
enlargingFactor (float) – if an iteration is very successful, the radius of the trust region is multiplied by this factor. Default 10.
 Returns
x, messages
x is the solution generated by the algorithm,
messages is a dictionary describing information about the lagorithm
 Return type
numpay.array, dict(str:object)
 Raises
biogeme.exceptions.biogemeError – if the dimensions of the matrix initBfgs do not match the length of x0.

biogeme.algorithms.
truncatedConjugateGradient
(g, H, delta)[source]¶ Find an approximation of the trust region subproblem using the truncated conjugate gradient method
 Parameters
g (numpy.array) – gradient of the quadratic model.
H (numpy.array) – hessian of the quadrartic model.
delta (float) – radius of the trust region.
 Returns
d, diagnostic, where
d is the approximate solution of the trust region subproblem,
diagnostic is the nature of the solution:
1 for convergence,
2 if out of the trust region,
3 if negative curvature detected.
4 if a numerical problem has been encountered
 Return type
numpy.array, int

biogeme.algorithms.
truncatedConjugateGradientSubspace
(xk, gk, Hk, delta, bounds, infeasibleIterate=False, tol=6.06273418136464e06)[source]¶ Find an approximation of the solution of the trust region subproblem using the truncated conjugate gradient method within the subspace of free variables. Free variables are those corresponding to inactive constraints at the generalized Cauchy point.
 Parameters
xk (numpy.array) – current iterate.
gk (numpy.array) – gradient of the quadratic model.
Hk (numpy.array) – hessian of the quadrartic model.
delta (float) – radius of the trust region.
bounds (class bioBounds) – bounds on the variables.
infeasibleIterate (bool) – if True, the algorithm may continue until termination. The result will then be projected on the feasible domain. If False, the algorithm stops as soon as an infeasible iterate is generated. Default: False.
tol (float) – tolerance used for stopping criterion.
 Returns
d, diagnostic, where
d is the approximate solution of the trust region subproblem,
diagnostic is the nature of the solution:
1 for convergence,
2 if out of the trust region,
3 if negative curvature detected.
4 if a numerical problem has been encountered
 Return type
numpy.array, int
 Raises
biogeme.exceptions.biogemeError – if the dimensions are inconsistent

biogeme.algorithms.
trustRegionIntersection
(dc, d, delta)[source]¶ Calculates the intersection with the boundary of the trust region.
Consider a trust region of radius \(\delta\), centered at \(\hat{x}\). Let \(x_c\) be in the trust region, and \(d_c = x_c  \hat{x}\), so that \(\d_c\ \leq \delta\). Let \(x_d\) be out of the trust region, and \(d_d = x_d  \hat{x}\), so that \(\d_d\ \geq \delta\). We calculate \(\lambda\) such that
\[\ d_c + \lambda (d_d  d_c)\ = \delta\] Parameters
dc (numpy.array) – xc  xhat.
d (numpy.array) – dd  dc.
delta (float) – radius of the trust region.
 Returns
\(\lambda\) such that \(\ d_c + \lambda (d_d  d_c)\ = \delta\)
 Return type
float