3.1.10. horton.quadprog – A light-weight quadratic programming solver

Problems of the following type can be solved:

\[\begin{split}min_{x} \frac{1}{2} x^T A x - b^T x \\ R x = s x \ge 0\end{split}\]

This is not the most general type of quadratic programming but it is sufficient for the needs in HORTON.

The equality constraints are optional. When A is positive definite, a polynomial-scaling interior point algorithm can be used: QPSolver.find_local. If A has negative eigenvalues, only the brute force solver is reliable but it’s cost scales exponentially with problem size.

No a priori feasibility tests are carried out and it is assumed that in QPSolver.find_local one is able to provide a feasible initial guess. The routine QPSolver.find_brute will complain a posteriori if the problem is infeasible.

When A is negative definite, the solution may be unbounded. In the QPSolver.find_local algorithm, this will lead to a convergence failure in a limited number of steps. In QPSolver.find_brute`, this situation is properly detected.

In summary, it is assumed that one only works with well-behaved problems. Issues with feasibilty and boundedness are not always handled smoothly.

3.1.10.1. Notation

The following naming conventions are used in the quadprog code:

free
Refers to the components of the solution that is not constrained to zero. The variable free is always an array with booleans indicating which components are free (True) or frozen (False).
frozen
Opposite of free.
x
A solution vector
nx
The dimension of the solution vector
nl
The number of constraints
class horton.quadprog.QPSolver(a, b, r=None, s=None, eps=1e-10)

Bases: object

The problem is defined as follows;

\[\begin{split}\min_x \frac{1}{2} x^T A x - b^T x \\ R x = s \\ x \ge 0\end{split}\]
Parameters:
  • a (np.ndarray) – The symmetric matrix \(A \in \mathbb{R}^{n \times n}\)
  • b (np.ndarray) – The column vector \(b \in \mathbb{R}^n\).
  • r (np.ndarray or None) – The matrix with constraint coefficients, \(R \in \mathbb{R}^{l \times n}\).
  • s (np.ndarray or None) – The matrix with constraint targets, \(s \in \mathbb{R}^l\).
  • eps (float) – A general threshold used for several purposes, e.g. the validity of a solution.
__init__(a, b, r=None, s=None, eps=1e-10)

The problem is defined as follows;

\[\begin{split}\min_x \frac{1}{2} x^T A x - b^T x \\ R x = s \\ x \ge 0\end{split}\]
Parameters:
  • a (np.ndarray) – The symmetric matrix \(A \in \mathbb{R}^{n \times n}\)
  • b (np.ndarray) – The column vector \(b \in \mathbb{R}^n\).
  • r (np.ndarray or None) – The matrix with constraint coefficients, \(R \in \mathbb{R}^{l \times n}\).
  • s (np.ndarray or None) – The matrix with constraint targets, \(s \in \mathbb{R}^l\).
  • eps (float) – A general threshold used for several purposes, e.g. the validity of a solution.
check_feasible(x, free=None)

Check if a solution, x, matches the constraints

Arguments:

x
The solution to test

Optional arguments:

free
When given, it is assumed that the fixed variables must be zero and that the rest can be anything (also negative). When not given, all components of x must be zero or positive.
check_solution(x)

Check if a solution, x, is a valid local minimum

compute_cn(free)

Compute the condition number of the problem.

Arguments:

free
Boolean array with the free components.
compute_cost(x)

Compute the function to be minimized for the given x

find_brute()

Brute force solution of the quadratic programming problem

Returns:

cost
The value of the cost function at the solution
x
The solution vector
find_local(x, trust_radius, maxiter=None)

A local solver for the quadratic programming problem

Arguments:

x
The initial guess for the solution. This must be a feasible point.
trust_radius
The maximum step size. This can be fairly large.

Optional arguments:

maxiter
The maximum number of iterations, this is 10*nx by default.

Returns:

cost
The value of the cost function at the solution
x
The solution vector

Note that this local search will always lead to the global optimum if the matrix \(A\) is positive definite.

get_free_problem(free)

Return the matrix a, b, r and s without the frozen columns/rows

Arguments:

free
Boolean array with the free components

Returns:

a2
The coefficients of the free system, including rows and colums for the lagrange multipliers.
b2
The right-hand side of the equations, including the targets of the equations.
nfree
The number of equations due to free components of x.
get_rmsds(x)

Quantify how far x deviates from a local solution

Returns:

gradient
The gradient projected onto the equality constraints, if any.
rmsd_free
The gradient RMSD for the non-zero components
rmsd_frozen
The gradient RMSD for the zero components. Only contributions from negative components are included.
rmsd_neg
The rmsd of the negative components of the solution.
log(x=None)

Print out the qps for debugging purposes

Optional arguments:

x
An solution vector (that causes problems).
solve(free)

Solve the problem, keeping certain components fixed at zero

The equality constraints are honored. However, free components may become negative. Fixed components are internally excluded from the equations and their result is manually set to zero.

Arguments:

free
Boolean array with the free components
solve_radius(free, center, radius)

Solve the equations, keeping certain components fixed at zero, with maximum radius.

The equality constraints are honored. However, free components may become negative. Fixed components are internally excluded from the equations and their result is manually set to zero.

Arguments:

free
Boolean array with the free components
center
The center of the restraint, a vector with length nx. (\(x_0\))
radius
The maximum Euclidean distance of the returned solution x from the given center
nl

The number of constraints

nx

The number of unkowns