EQP#

purpose#

The eqp package uses an iterative method to solve a given equality-constrained quadratic program. The aim is to minimize the quadratic objective function

\[q(x) = f + g^T x + \frac{1}{2} x^T H x,\]

or the shifted-least-distance objective function

\[s(x) = f + g^T x + \frac{1}{2} \sum_{j=1}^n w_j^2 (x_j - x_j^0)^2,\]

subject to the general linear equality constraints

\[A x + c = 0,\]

where $H$ and $A$ are, respectively, given $n$ by $n$ symmetric and $m$ by $n$ general matrices, $g$, $w$, $x^0$ and $c$ are vectors, and $f$ is a scalar. The method is most suitable for problems involving a large number of unknowns $x$.

See Section 4 of $GALAHAD/doc/eqp.pdf for additional details.

terminology#

Any required solution $x$ necessarily satisfies the primal optimality conditions

\[A x + c = 0 \;\;\mbox{(1)}\]

and the dual optimality conditions

\[H x + g = A^T y,\]

where the vector $y$ is known as the Lagrange multipliers for the general linear constraints.

In the shifted-least-distance case, $g$ is shifted by $-W^2 x^0$, and $H = W^2$, where $W$ is the diagonal matrix whose entries are the $w_j$.

method#

A solution to the problem is found in two phases. In the first, a point $x_F$ satisfying (1) is found. In the second, the required solution $x = x_F + s$ is determined by finding $s$ to minimize $q(s) = \frac{1}{2} s^T H s + g_F^T s + f_F^{}$ subject to the homogeneous constraints $A s = 0$, where $g_F^{} = H x_F^{} + g$ and $f_F^{} = \frac{1}{2} x_F^T H x_F^{} + g^T x_F^{} + f$. The required constrained minimizer of $q(s)$ is obtained by implictly applying the preconditioned conjugate-gradient method in the null space of $A$. Any preconditioner of the form

\[\begin{split}K_{G} = \left(\begin{array}{cc} G & A^T \\ A & 0 \end{array}\right)\end{split}\]

is suitable, and the package SBLS provides a number of possibilities. In order to ensure that the minimizer obtained is finite, an additional, precautionary trust-region constraint $\|s\| \leq \Delta$ for some suitable positive radius $\Delta$ is imposed, and the package GLTR is used to solve this additionally-constrained problem.

references#

The preconditioning aspcets are described in detail in

H. S. Dollar, N. I. M. Gould and A. J. Wathen. ``On implicit-factorization constraint preconditioners’’. In Large Scale Nonlinear Optimization (G. Di Pillo and M. Roma, eds.) Springer Series on Nonconvex Optimization and Its Applications, Vol. 83, Springer Verlag (2006) 61–82

and

H. S. Dollar, N. I. M. Gould, W. H. A. Schilders and A. J. Wathen ``On iterative methods and implicit-factorization preconditioners for regularized saddle-point systems’’. SIAM Journal on Matrix Analysis and Applications 28(1) (2006) 170–189,

while the constrained conjugate-gradient method is discussed in

N. I. M. Gould, S. Lucidi, M. Roma and Ph. L. Toint, ``Solving the trust-region subproblem using the Lanczos method’’. SIAM Journal on Optimization 9(2) (1999), 504–525.

matrix storage#

The unsymmetric $m$ by $n$ matrix $A$ may be presented and stored in a variety of convenient input formats.

Dense storage format: The matrix $A$ is stored as a compact dense matrix by rows, that is, the values of the entries of each row in turn are stored in order within an appropriate real one-dimensional array. In this case, component $n \ast i + j$ of the storage array A_val will hold the value $A_{ij}$ for $0 \leq i \leq m-1$, $0 \leq j \leq n-1$. The string A_type = ‘dense’ should be specified.

Dense by columns storage format: The matrix $A$ is stored as a compact dense matrix by columns, that is, the values of the entries of each column in turn are stored in order within an appropriate real one-dimensional array. In this case, component $m \ast j + i$ of the storage array A_val will hold the value $A_{ij}$ for $0 \leq i \leq m-1$, $0 \leq j \leq n-1$. The string A_type = ‘dense_by_columns’ should be specified.

Sparse co-ordinate storage format: Only the nonzero entries of the matrices are stored. For the $l$-th entry, $0 \leq l \leq ne-1$, of $A$, its row index i, column index j and value $A_{ij}$, $0 \leq i \leq m-1$, $0 \leq j \leq n-1$, are stored as the $l$-th components of the integer arrays A_row and A_col and real array A_val, respectively, while the number of nonzeros is recorded as A_ne = $ne$. The string A_type = ‘coordinate’should be specified.

Sparse row-wise storage format: Again only the nonzero entries are stored, but this time they are ordered so that those in row i appear directly before those in row i+1. For the i-th row of $A$ the i-th component of the integer array A_ptr holds the position of the first entry in this row, while A_ptr(m) holds the total number of entries. The column indices j, $0 \leq j \leq n-1$, and values $A_{ij}$ of the nonzero entries in the i-th row are stored in components l = A_ptr(i), $\ldots$, A_ptr(i+1)-1, $0 \leq i \leq m-1$, of the integer array A_col, and real array A_val, respectively. For sparse matrices, this scheme almost always requires less storage than its predecessor. The string A_type = ‘sparse_by_rows’ should be specified.

Sparse column-wise storage format: Once again only the nonzero entries are stored, but this time they are ordered so that those in column j appear directly before those in column j+1. For the j-th column of $A$ the j-th component of the integer array A_ptr holds the position of the first entry in this column, while A_ptr(n) holds the total number of entries. The row indices i, $0 \leq i \leq m-1$, and values $A_{ij}$ of the nonzero entries in the j-th columnsare stored in components l = A_ptr(j), $\ldots$, A_ptr(j+1)-1, $0 \leq j \leq n-1$, of the integer array A_row, and real array A_val, respectively. As before, for sparse matrices, this scheme almost always requires less storage than the co-ordinate format. The string A_type = ‘sparse_by_columns’ should be specified.

The symmetric $n$ by $n$ matrix $H$ may also be presented and stored in a variety of formats. But crucially symmetry is exploited by only storing values from the lower triangular part (i.e, those entries that lie on or below the leading diagonal).

Dense storage format: The matrix $H$ is stored as a compact dense matrix by rows, that is, the values of the entries of each row in turn are stored in order within an appropriate real one-dimensional array. Since $H$ is symmetric, only the lower triangular part (that is the part $H_{ij}$ for $0 \leq j \leq i \leq n-1$) need be held. In this case the lower triangle should be stored by rows, that is component $i * i / 2 + j$ of the storage array H_val will hold the value $H_{ij}$ (and, by symmetry, $H_{ji}$) for $0 \leq j \leq i \leq n-1$. The string H_type = ‘dense’ should be specified.

Sparse co-ordinate storage format: Only the nonzero entries of the matrices are stored. For the $l$-th entry, $0 \leq l \leq ne-1$, of $H$, its row index i, column index j and value $H_{ij}$, $0 \leq j \leq i \leq n-1$, are stored as the $l$-th components of the integer arrays H_row and H_col and real array H_val, respectively, while the number of nonzeros is recorded as H_ne = $ne$. Note that only the entries in the lower triangle should be stored. The string H_type = ‘coordinate’ should be specified.

Sparse row-wise storage format: Again only the nonzero entries are stored, but this time they are ordered so that those in row i appear directly before those in row i+1. For the i-th row of $H$ the i-th component of the integer array H_ptr holds the position of the first entry in this row, while H_ptr(n) holds the total number of entries. The column indices j, $0 \leq j \leq i$, and values $H_{ij}$ of the entries in the i-th row are stored in components l = H_ptr(i), …, H_ptr(i+1)-1 of the integer array H_col, and real array H_val, respectively. Note that as before only the entries in the lower triangle should be stored. For sparse matrices, this scheme almost always requires less storage than its predecessor. The string H_type = ‘sparse_by_rows’ should be specified.

Diagonal storage format: If $H$ is diagonal (i.e., $H_{ij} = 0$ for all $0 \leq i \neq j \leq n-1$) only the diagonals entries $H_{ii}$, $0 \leq i \leq n-1$ need be stored, and the first n components of the array H_val may be used for the purpose. The string H_type = ‘diagonal’ should be specified.

Multiples of the identity storage format: If $H$ is a multiple of the identity matrix, (i.e., $H = \alpha I$ where $I$ is the n by n identity matrix and $\alpha$ is a scalar), it suffices to store $\alpha$ as the first component of H_val. The string H_type = ‘scaled_identity’ should be specified.

The identity matrix format: If $H$ is the identity matrix, no values need be stored. The string H_type = ‘identity’ should be specified.

The zero matrix format: The same is true if $H$ is the zero matrix, but now the string H_type = ‘zero’ or ‘none’ should be specified.

functions#

eqp.initialize()#

Set default option values and initialize private data

Returns:

optionsdict

dictionary containing default control options:

errorint
error and warning diagnostics occur on stream error.

outint
general output occurs on stream out.

print_levelint
the level of output required is specified by print_level. Possible values are

<=0

gives no output,

1

1 gives a one-line summary for every iteration.

2

gives a summary of the inner iteration for each iteration.

>=3

gives increasingly verbose (debugging) output.

factorizationint
the factorization to be used. Possible values are

0

automatic

1

Schur-complement factorization

2

augmented-system factorization (obsolete).

max_colint
the maximum number of nonzeros in a column of A which is permitted with the Schur-complement factorization (obsolete).

indminint
an initial guess as to the integer workspace required by SBLS (obsolete).

valminint
an initial guess as to the real workspace required by SBLS (obsolete).

len_ulsminint
an initial guess as to the workspace required by ULS (obsolete).

itref_maxint
the maximum number of iterative refinements allowed (obsolete).

cg_maxitint
the maximum number of CG iterations allowed. If cg_maxit < 0, this number will be reset to the dimension of the system + 1.

preconditionerint
the preconditioner to be used for the CG. Possible values are

0

automatic.

1

no preconditioner, i.e, the identity within full factorization.

2

full factorization.

3

band within full factorization.

4

diagonal using the barrier terms within full factorization. (obsolete)

5

optionally supplied diagonal, G = D.

semi_bandwidthint
the semi-bandwidth of a band preconditioner, if appropriate (obsolete).

new_aint
how much has A changed since last problem solved: 0 = not changed, 1 = values changed, 2 = structure changed.

new_hint
how much has H changed since last problem solved: 0 = not changed, 1 = values changed, 2 = structure changed.

sif_file_deviceint
specifies the unit number to write generated SIF file describing the current problem.

pivot_tolfloat
the threshold pivot used by the matrix factorization. See the documentation for SBLS for details (obsolete).

pivot_tol_for_basisfloat
the threshold pivot used by the matrix factorization when finding the ba See the documentation for ULS for details (obsolete).

zero_pivotfloat
any pivots smaller than zero_pivot in absolute value will be regarded to zero when attempting to detect linearly dependent constraints (obsolete).

inner_fraction_optfloat
the computed solution which gives at least inner_fraction_opt times the optimal value will be found (obsolete).

radiusfloat
an upper bound on the permitted step (-ve will be reset to an appropriat large value by eqp_solve).

min_diagonalfloat
diagonal preconditioners will have diagonals no smaller than min_diagonal (obsolete).

max_infeasibility_relativefloat
if the constraints are believed to be rank defficient and the residual at a “typical” feasible point is larger than max( max_infeasibility_relative * norm A, max_infeasibility_absolute ) the problem will be marked as infeasible.

max_infeasibility_absolutefloat
see max_infeasibility_relative.

inner_stop_relativefloat
the computed solution is considered as an acceptable approximation to th minimizer of the problem if the gradient of the objective in the preconditioning(inverse) norm is less than max( inner_stop_relative * initial preconditioning(inverse) gradient norm, inner_stop_absolute ).

inner_stop_absolutefloat
see inner_stop_relative.

inner_stop_interfloat
see inner_stop_relative.

find_basis_by_transposebool
if find_basis_by_transpose is True, implicit factorization precondition will be based on a basis of A found by examining A’s transpose (obsolete).

remove_dependenciesbool
if remove_dependencies is True, the equality constraints will be preprocessed to remove any linear dependencies.

space_criticalbool
if space_critical is True, every effort will be made to use as little space as possible. This may result in longer computation time.

deallocate_error_fatalbool
if deallocate_error_fatal is True, any array/pointer deallocation error will terminate execution. Otherwise, computation will continue.

generate_sif_filebool
if generate_sif_file is True, a SIF file describing the current problem is to be generated.

sif_file_namestr
name of generated SIF file containing input problem.

prefixstr
all output lines will be prefixed by the string contained in quotes within prefix, e.g. ‘word’ (note the qutoes) will result in the prefix word.

fdc_optionsdict
default control options for FDC (see fdc.initialize).

sbls_optionsdict
default control options for SBLS (see sbls.initialize).

gltr_optionsdict
default control options for GLTR (see gltr.initialize).

eqp.load(n, m, H_type, H_ne, H_row, H_col, H_ptr, A_type, A_ne, A_row, A_col, A_ptr, options=None)#

Import problem data into internal storage prior to solution.

Parameters:

nint
holds the number of variables.

mint
holds the number of constraints.

H_typestring
specifies the symmetric storage scheme used for the Hessian $H$. It should be one of ‘coordinate’, ‘sparse_by_rows’, ‘dense’, ‘diagonal’, ‘scaled_identity’, ‘identity’, ‘zero’ or ‘none’; lower or upper case variants are allowed.

H_neint
holds the number of entries in the lower triangular part of $H$ in the sparse co-ordinate storage scheme. It need not be set for any of the other schemes.

H_rowndarray(H_ne)
holds the row indices of the lower triangular part of $H$ in the sparse co-ordinate storage scheme. It need not be set for any of the other schemes, and in this case can be None.

H_colndarray(H_ne)
holds the column indices of the lower triangular part of $H$ in either the sparse co-ordinate, or the sparse row-wise storage scheme. It need not be set when the other storage schemes are used, and in this case can be None.

H_ptrndarray(n+1)
holds the starting position of each row of the lower triangular part of $H$, as well as the total number of entries, in the sparse row-wise storage scheme. It need not be set when the other schemes are used, and in this case can be None.

A_typestring
specifies the unsymmetric storage scheme used for the constraints Jacobian $A$. It should be one of ‘coordinate’, ‘sparse_by_rows’ or ‘dense’; lower or upper case variants are allowed.

A_neint
holds the number of entries in $A$ in the sparse co-ordinate storage scheme. It need not be set for any of the other two schemes.

A_rowndarray(A_ne)
holds the row indices of $A$ in the sparse co-ordinate storage scheme. It need not be set for any of the other two schemes, and in this case can be None.

A_colndarray(A_ne)
holds the column indices of $A$ in either the sparse co-ordinate, or the sparse row-wise storage scheme. It need not be set when the dense storage scheme is used, and in this case can be None.

A_ptrndarray(m+1)
holds the starting position of each row of $A$, as well as the total number of entries, in the sparse row-wise storage scheme. It need not be set when the other schemes are used, and in this case can be None.

optionsdict, optional
dictionary of control options (see eqp.initialize).

eqp.solve_qp(n, m, f, g, H_ne, H_val, A_ne, A_val, c, x, y)#

Find a solution to the convex quadratic program involving the quadratic objective function $q(x)$.

Parameters:

nint
holds the number of variables.

mint
holds the number of residuals.

ffloat
holds the constant term $f$ in the objective function.

gndarray(n)
holds the values of the linear term $g$ in the objective function.

H_neint
holds the number of entries in the lower triangular part of the Hessian $H$.

H_valndarray(H_ne)
holds the values of the nonzeros in the lower triangle of the Hessian $H$ in the same order as specified in the sparsity pattern in eqp.load.

A_neint
holds the number of entries in the constraint Jacobian $A$.

A_valndarray(A_ne)
holds the values of the nonzeros in the constraint Jacobian $A$ in the same order as specified in the sparsity pattern in eqp.load.

cndarray(m)
holds the values of the linear term $c$ in the constraints

xndarray(n)
holds the initial estimate of the minimizer $x$, if known. This is not crucial, and if no suitable value is known, then any value, such as $x=0$, suffices and will be adjusted accordingly.

yndarray(m)
holds the initial estimate of the Lagrange multipliers $y$ associated with the general constraints, if known. This is not crucial, and if no suitable value is known, then any value, such as $y=0$, suffices and will be adjusted accordingly.

Returns:

xndarray(n)
holds the values of the approximate minimizer $x$ after a successful call.

yndarray(m)
holds the values of the Lagrange multipliers associated with the linear constraints.

eqp.solve_sldqp(n, m, f, g, w, x0, A_ne, A_val, c, x, y)#

Find a solution to the quadratic program involving the shifted least-distance objective function $s(x)$.

Parameters:

nint
holds the number of variables.

mint
holds the number of residuals.

ffloat
holds the constant term $f$ in the objective function.

gndarray(n)
holds the values of the linear term $g$ in the objective function.

wndarray(n)
holds the values of the weights $w$ in the objective function.

x0ndarray(n)
holds the values of the shifts $x^0$ in the objective function.

A_neint
holds the number of entries in the constraint Jacobian $A$.

A_valndarray(A_ne)
holds the values of the nonzeros in the constraint Jacobian $A$ in the same order as specified in the sparsity pattern in eqp.load.

cndarray(m)
holds the values of the linear term $c$ in the constraints

xndarray(n)
holds the initial estimate of the minimizer $x$, if known. This is not crucial, and if no suitable value is known, then any value, such as $x=0$, suffices and will be adjusted accordingly.

yndarray(m)
holds the initial estimate of the Lagrange multipliers $y$ associated with the general constraints, if known. This is not crucial, and if no suitable value is known, then any value, such as $y=0$, suffices and will be adjusted accordingly.

Returns:

xndarray(n)
holds the values of the approximate minimizer $x$ after a successful call.

yndarray(m)
holds the values of the Lagrange multipliers associated with the general linear constraints.

eqp.resolve_qp(n, m, f, g, c, x, y)#

Resolve the convex quadratic program when the data $f$, $g$ and/or $c$ has changed.

Parameters:

nint
holds the number of variables.

mint
holds the number of residuals.

ffloat
holds the constant term $f$ in the objective function.

gndarray(n)
holds the values of the linear term $g$ in the objective function.

cndarray(m)
holds the values of the linear term $c$ in the constraints

xndarray(n)
holds the initial estimate of the minimizer $x$, if known. This is not crucial, and if no suitable value is known, then any value, such as $x=0$, suffices and will be adjusted accordingly.

yndarray(m)
holds the initial estimate of the Lagrange multipliers $y$ associated with the general constraints, if known. This is not crucial, and if no suitable value is known, then any value, such as $y=0$, suffices and will be adjusted accordingly.

Returns:

xndarray(n)
holds the values of the approximate minimizer $x$ after a successful call.

yndarray(m)
holds the values of the Lagrange multipliers associated with the linear constraints.

[optional] eqp.information()

Provide optional output information

Returns:

informdict

dictionary containing output information:

statusint
return status. Possible values are:

0

The run was successful.

-1

An allocation error occurred. A message indicating the offending array is written on unit options[‘error’], and the returned allocation status and a string containing the name of the offending array are held in inform[‘alloc_status’] and inform[‘bad_alloc’] respectively.

-2

A deallocation error occurred. A message indicating the offending array is written on unit options[‘error’] and the returned allocation status and a string containing the name of the offending array are held in inform[‘alloc_status’] and inform[‘bad_alloc’] respectively.

-3

The restriction n > 0 or m > 0 or requirement that type contains its relevant string ‘dense’, ‘coordinate’, ‘sparse_by_rows’, ‘diagonal’, ‘scaled_identity’, ‘identity’, ‘zero’ or ‘none’ has been violated.

-5

The constraints appear to be inconsistent.

-7

The objective function appears to be unbounded from below on the feasible set.

-9

The analysis phase of a symmetric factorization failed; the return status from the factorization package is given by inform[‘factor_status’].

-10

A symmetric factorization failed; the return status from the factorization package is given by inform[‘factor_status’].

-11

The solution of a set of linear equations using factors from a symmetric factorization package failed; the return status from the factorization package is given by inform[‘factor_status’].

-12

The analysis phase of an unsymmetric factorization failed; the return status from the factorization package is given by inform[‘factor_status’].

-13

An unsymmetric factorization failed; the return status from the factorization package is given by inform[‘factor_status’].

-14

The solution of a set of linear equations using factors from an unsymmetric factorization package failed; the return status from the factorization package is given by inform[‘factor_status’].

-15

The computed preconditioner has the wrong inertia and is thus unsuitable.

-16

The residuals from the preconditioning step are large, indicating that the factorization may be unsatisfactory.

-25

eqp.resolve has been called before eqp.solve.

alloc_statusint
the status of the last attempted allocation/deallocation.

bad_allocstr
the name of the array for which an allocation/deallocation error occurred.

cg_iterint
the total number of conjugate gradient iterations required.

cg_iter_interint
see cg_iter.

factorization_integerlong
the total integer workspace required for the factorization.

factorization_reallong
the total real workspace required for the factorization.

objfloat
the value of the objective function at the best estimate of the solution determined by QPB_solve.

timedict

dictionary containing timing information:

totalfloat
the total CPU time spent in the package.

find_dependentfloat
the CPU time spent detecting linear dependencies.

factorizefloat
the CPU time spent factorizing the required matrices.

solvefloat
the CPU time spent computing the search direction.

solve_interfloat
see solve.

clock_totalfloat
the total clock time spent in the package.

clock_find_dependentfloat
the clock time spent detecting linear dependencies.

clock_factorizefloat
the clock time spent factorizing the required matrices.

clock_solvefloat
the clock time spent computing the search direction.

fdc_informdict
inform parameters for FDC (see fdc.information).

sbls_informdict
inform parameters for SBLS (see sbls.information).

gltr_informdict
inform parameters for GLTR (see gltr.information).

eqp.terminate()#

Deallocate all internal private storage.

example code#

from galahad import eqp
import numpy as np
np.set_printoptions(precision=4,suppress=True,floatmode='fixed')
print("\n** python test: eqp")

# set parameters
n = 3
m = 2

#  describe objective function

f = 1.0
g = np.array([0.0,2.0,0.0])
H_type = 'coordinate'
H_ne = 3
H_row = np.array([0,1,2])
H_col = np.array([0,1,2])
H_ptr = None
H_val = np.array([1.0,1.0,1.0])

#  describe constraints

A_type = 'coordinate'
A_ne = 4
A_row = np.array([0,0,1,1])
A_col = np.array([0,1,1,2])
A_ptr = None
A_val = np.array([2.0,1.0,1.0,1.0])
c = np.array([3.0,0.0])

# allocate internal data and set default options
options = eqp.initialize()

# set some non-default options
options['print_level'] = 0
options['fdc_options']['use_sls'] = True
options['fdc_options']['symmetric_linear_solver'] = 'sytr '
options['sbls_options']['symmetric_linear_solver'] = 'sytr '
options['sbls_options']['definite_linear_solver'] = 'sytr '
#print("options:", options)

# load data (and optionally non-default options)
eqp.load(n, m, H_type, H_ne, H_row, H_col, H_ptr,
         A_type, A_ne, A_row, A_col, A_ptr, options)

#  provide starting values (not crucial)

x = np.array([0.0,0.0,0.0])
y = np.array([0.0,0.0])

# find optimum of qp
print("\n 1st problem: solve qp")
x, y = eqp.solve_qp(n, m, f, g, H_ne, H_val, A_ne, A_val, c, x, y)
print("x:",x)
print("y:",y)

# get information
inform = eqp.information()
print("f:",inform['obj'])

# deallocate internal data

eqp.terminate()

#  describe shifted-least-distance qp

w = np.array([1.0,1.0,1.0])
x0 = np.array([0.0,0.0,0.0])
H_type = 'shifted_least_distance'

# allocate internal data
eqp.initialize()

# load data (and optionally non-default options)
eqp.load(n, m, H_type, H_ne, H_row, H_col, H_ptr,
         A_type, A_ne, A_row, A_col, A_ptr, options)

#  provide starting values (not crucial)

x = np.array([0.0,0.0,0.0])
y = np.array([0.0,0.0])

# find optimum of sldqp
print("\n 2nd problem: solve sldqp")
x, y = eqp.solve_sldqp(n, m, f, g, w, x0, A_ne, A_val, c, x, y)
print("x:",x)
print("y:",y)

# get information
inform = eqp.information()
print(" f: %.4f" % inform['obj'])
print('** eqp exit status:', inform['status'])

# deallocate internal data

eqp.terminate()

This example code is available in $GALAHAD/src/eqp/Python/test_eqp.py .