Least-squares optimization in Python with single equal constraint

user_13037272 0
Viewed 376

I am looking to solve a least-squares problem in python such that I minimize 0.5*||np.dot(A,x) - b||^2 subject to the constraint np.dot(C,x) = d, with the bounds 0<x<np.inf, where

A : shape=[m, n]
C : shape=[m, n]
b : shape=[m]
d : shape=[m]

are all known matrices.

Unfortunately, it looks like the scipy.optimize.lsq_linear() function only works for an upper/lower bound constraint:

minimize 0.5 * ||A x - b||**2
subject to lb <= x <= ub

not an equality constraint. In addition, I would like to add bounds for the solution such that it is positive ONLY. Is there an easy or clean way to do this using scipy or NumPy built-in functions?

least-squares scipy numpy optimization python
1 Answers

If you stick to scipy.optimize.lsq_linear, you could add a penalty term to the objective function, i.e.

minimize 0.5 * ||A x - b||**2 + beta*||Cx-d||**2
subject to lb <= x <= ub

and choose the scalar beta sufficiently large in order to transform your constrained optimization problem into a problem with simple box-constraints. However, it is more convenient to use scipy.optimize.minimize for constrained optimization problems:

import numpy as np
from scipy.optimize import minimize

# ... Define your np.ndarrays A, C, b, d here

# constraint
cons = [{'type': 'eq', 'fun': lambda x: C @ x - d}]

# variable bounds
bnds = [(0, None) for _ in range(A.shape[1])]

# initial point
x0 = np.ones(A.shape[1])

# res.x contains your solution
res = minimize(lambda x: np.linalg.norm(A@x-b)**2, x0=x0, bounds=bnds, constraints=cons)

where I used the euclidian norm. Timing this for a toy example with m = 2000 and n = 1000 yields 3min and 10s on my machine.

Note that both the objective gradient and the constraint Jacobian are approximated by finite differences, which requires many evaluations of the objective and the constraint function. Consequently, the code might be too slow for large problems. In this case, it's worth providing the exact gradient, jacobian and hessian. Furthermore, since

# f(x) = ||Ax-b||^2 = x'A'Ax - 2b'Ax + b'b
# grad f(x) = 2A'Ax - 2A'b
# hess f(x) = 2A'A

we can significantly reduce the time needed to evaluate the functions by precalculating A'A, b'A and b'b. Additionally, we only need to calculate A'Ax once:

from scipy.optimize import NonlinearConstraint

# Precalculate the matrix-vector products
AtA = A.T @ A
btA = b.T @ A
btb = b.T @ b
Atb = btA.T

# return the objective value and the objective gradient
def obj_and_grad(x):
    AtAx = AtA @ x
    obj = x.T @ AtAx - 2*btA @ x + btb
    grad = 2*AtAx - 2*Atb
    return obj, grad

#constraint d <= C@x <= d (including the jacobian and hessian)
con = NonlinearConstraint(lambda x: C @ x, d, d, jac=lambda x: C, hess=lambda x, v: v[0] * np.zeros(x.size))

# res.x contains your solution
res = minimize(obj_and_grad, jac=True, hess=lambda x: AtA2, x0=x0, bounds=bnds, constraints=(con,), method="trust-constr")

Here, the option jac=True tells minimize that our function obj_and_grad returns a tuple containing the objective function and the gradient. We further choose the 'trust-constr' method, since it is the only one available that supports hessians. Timing this again for the same aforementioned toy example, yields 7s on my machine.

user_47455290
Related