Constrained Linear Optimization problem in python

Constraints

Presuming that you are referring to a linear system of equations of the form:

ax = b

where a is a matrix of values (coefficients) in your code snippet, and b is a vector of constants (maxPoints) in your code, then you are solving for x.

Adding constraints doesn't really make sense here, as the system is soluble if:

1. coefficients is invertible (or at least, has a pseudoinverse).
2. All elements of maxPoints are finite.

The LinAlgError

In essence, you can think of numpy.linalg.solve(a,b) doing the following:

Solve: ax = b for x

One can think of this in a sense as numpy computing the inverse of the matrix a. Here, numpy raised a LinAlgError because the matrix was not square, and therefore the inverse is poorly defined. numpy actually does some other little bits behind the scenes (see LAPACK), but this is close enough for discussion here.

What I think may be the issue

You have the following problem:

• The number of columns in coefficients is 6, and the number of rows is 32.
• The number of rows in maxPoints is 32.
• This is an over-determined problem.

A first step would be to use numpy.linalg.lstsq (see here). You can use exactly the same syntax.

x_ls = np.linalg.lstsq(coefficients, maxPoints)

This solver yields an answer which minimizes the squared Euclidean 2-norm. Instead you want to minimize using the constraint that x[j] <= x[j+1] for all j < len(x).

Some steps towards a solution

I may come back and edit this later with a Github snippet, as it's quite an interesting problem.

Edit: trust-constr seemed to work much better for this than SLSQP.

I think you should have a bit of a further read of the optimization documentation in scipy, particularly the trust-constr algorithm. I think you pretty much want to do a least squares solution on part of your matrix, then iterate forward. As a plausible algorithm:

1. Solve using trust-constr on the problem with 2 columns only (corresponding to x[4] and x[5]). Bounds (choose BOUND freely):
   • x[4]: [-np.inf,BOUND]
   • x[5]: [BOUND, np.inf]
2. Solve again with 3 columns (x3, x[4], x[5]). Bounds (with _p denoting value obtained from previous step):
   • x3: [-np.inf,x_p[4]]
   • x[4]: [x_p[4],x_p[5]]
   • x[5]: [x_p[5],np.inf]
3. Repeat similarly until you solve for the full 6 columns.

It's not the best algorithm in the world, however, your constraint will be satisfied, and the iterative nature will correct a little bit if you do pick a bad value of BOUND to start with. If you wish to obtain a global solution, you can repeat for many different values of bound, and see which one returns the lowest error in the final trust-constr fit.

Edit: for an implementation:

import numpy as np
from scipy.optimize import minimize
from scipy.optimize import Bounds

# Set up bounds and make a,b objects
a,b = coefficients, maxPoints[:,0]

# Set acceptable bounds to be whatever you like here.
# I've chosen -100000 to 100000 in steps of 10000.
iter_of_acceptable_bounds = range(-100000,100000,10000)

def e2n(x,*args):
    'Compute the sum of the residuals (Euclidean 2-norm).'
    a,b = args
    return np.inner(b-a@x,b-a@x)

solutions, errors = [], []
for BOUND in iter_of_acceptable_bounds:   
    # Repeat the algorithm for each initial choice of bounds.

    # Generate an array for the bounds, based no the initial choice.
    bounds_arr = np.concatenate(([-np.inf],[BOUND,BOUND],[np.inf]))

    for i in range(2,a.shape[1]+1):
        # Initial guess of 0.
        x0 = [0 for j in range(i)]

        # Start with 2 columns, end at a.shape[1]
        # (a.shape[1]==6 for the example).

        bounds = Bounds(
                *bounds_arr.reshape(i,2).T.tolist()
                )

        # Slice up a accordingly
        A = a[:,-i:]

        # Minimize using trust-constr, and unpack the solution
        xr = minimize(e2n,x0,
                 args=(A,b),
                 bounds=bounds,
                 method='trust-constr').x

        # Generate a new array for the bounds, dynamically based on the previous choices.
        bounds_arr = np.concatenate(
                        [[-np.inf]] + [(xr[j],xr[j]) for j in range(i)] + [[np.inf]]
                        )

    # Save the solution and error from the full matrix
    solutions.append(xr)
    errors.append(np.sqrt(e2n(xr,a,b)))

This seems to perform fairly well. If we make a quick plot:

import matplotlib.pyplot as plt

# lstsq, for comparison (no constraints)
lstsq_fit, lstsq_res, *_ = np.linalg.lstsq(a,b)
lstsq_err = np.sqrt(lstsq_res)

# Plot errors
plt.plot(iter_of_acceptable_bounds,100*(errors-lstsq_err)/lstsq_err)
plt.xlabel('Starting choice of BOUND for algorithm')
plt.ylabel('Increase in error compared to numpy.linalg.lstsq')
plt.show()

In general, numpy.linalg.lstsq will solve this problem fairly well (it is just a linear regression problem), however it doesn't have the desired constraint. However, generally speaking this algorithm is okay -- the data fits only around 40% worse than numpy.linalg.lstsq for some of our choices of BOUND.

From this coarse search, the best values were given by:

>>> my_sol = solutions[errors.index(min(errors))]
>>> print(my_sol) 
array([-186017.00778511,    6680.13364168,  184099.89046232,
        220587.55247874,  257275.09670101,  257275.09670101])