Solve parameter matrix with element-wise constraints using scipy optimize

Question

I have a relatively simple optimisation problem, but am not particularly well-versed in scipy and cannot figure out how to apply the necessary constraints. My objective function is to minimise the absolute distance between a 10-element vector y and x, where x is a weighted row-sum of the 10x3 parameter matrix p

x = p[:,0] + 2 * p[:,1] + 3 * p[:,2]

I need to add the following constraints to the parameter matrix p:

Each element of p, but be >=0 and <=1
The sum of each column of p must sum to 1
The sum of each row of p must not exceed 1

I have attempted to implement constraints based on other SO questions, although I confess I don't fully understand it, nor the errors it generates.

# import packages
import numpy as np
from scipy.optimize import minimize

# create objective function
def objective(p, y):
    x = p[:,0] + 2*p[:,1] + 3*p[:,2]
    return np.sum(np.abs(y-x))

# define constraints
cons = [{'type':'ineq', 'fun': lambda x: x}, 
        {'type':'ineq', 'fun': lambda x: np.sum(x, 0)},   # row sum >= 0
        {'type':'ineq', 'fun': lambda x: 1 - np.sum(x, 0)},   # row sum <= 1
        {'type':'ineq', 'fun': lambda x: np.sum(x, 1) - 1},   # row sum >= 1
        {'type':'ineq', 'fun': lambda x: 1 - np.sum(x, 1)}]  # row sum <= 1

# generate target data (y)
y = np.array([2.48, 1.75, 1.13, 0.20, 0.20, 0.20, 0.0, 0.0, 0.0, 0.0])

# initialise starting param values
p0 = np.zeros((44,3))
p0[:,:] = 1/44

# solve
sol = minimize(objective, p0, method='SLSQP', constraints=cons)

Upon running this code, I get the following error message:

AxisError: axis 1 is out of bounds for array of dimension 1

Any help is much appreciated.

If there is some error, show the error. I also highly recommend formulating it as LP and not using minimize (yes: the abs will need some reformulation; but that is then easy to google). — sascha, Sep 19 '19 at 07:46
have included error msg - thanks. However i think the implementation of the constraints is the issue. — user3725021, Sep 19 '19 at 07:55
The core problem is that minimize is completely based on `x` being flat / 1d. It's no matrix / not multidimensional. Same goes for `p0`. This will effect in your problem as flattening will occur. Of course your matrix-like constraints will get more awkward as you need to do the 1d/2d mapping extre then, — sascha, Sep 19 '19 at 08:05

score 2 · Accepted Answer · answered Sep 22 '19 at 08:18

This answer is based on SuperKogito's answer:

import numpy as np
from scipy.optimize import minimize

# create objective function
def objective(p, y):
    sep = len(y)
    x   = p[0:sep] + 2*p[sep:2*sep] + 3*p[2*sep:3*sep]
    return np.sum(np.abs(y-x))

def vec2mat(vec):
    sep = vec.size // 3
    return vec.reshape(sep, 3)

# define constraints
cons = [{'type':'eq', 'fun': lambda x: np.sum(vec2mat(x), axis=0) - 1},# sum cols == 1
        {'type':'ineq', 'fun': lambda x: 1-np.sum(vec2mat(x), axis=1)} # sum rows <= 1
        ]

# generate target data (y)
y = np.array([2.48, 1.75, 1.13, 0.20, 0.20, 0.20, 0.0, 0.0, 0.0, 0.0])

# initialise starting param values
p0      = np.zeros((10,3))
p0[:,0] = 0.001
p0[:,1] = 0.002
p0[:,2] = 0.003

# Bounds: 0 <= x_i <= 1
bnds = [(0, 1) for _ in range(np.size(p0))]

# solve
sol = minimize(objective, x0 = p0.flatten('F'), args=(y,), bounds=bnds, constraints=cons)
print(sol)

# reconstruct p
p_sol = vec2mat(sol.x)

gives me

     fun: 1.02348743648716e-05
     jac: array([-1.,  1., -1.,  1., -1.,  1.,  1.,  1.,  1.,  1., -2.,  2., -2.,
        2., -2.,  2.,  2.,  2.,  2.,  2., -3.,  3., -3.,  3., -3.,  3.,
        3.,  3.,  3.,  3.])
 message: 'Optimization terminated successfully.'
    nfev: 1443
     nit: 43
    njev: 43
  status: 0
 success: True
       x: array([3.38285914e-01, 2.39989678e-01, 1.67911460e-01, 5.98105349e-02,
       4.84049819e-02, 1.44281667e-01, 1.42377791e-15, 1.74446343e-15,
       8.10053562e-16, 1.28295919e-15, 4.76418211e-01, 2.59584012e-01,
       2.20053270e-01, 5.22055787e-02, 2.00046857e-02, 1.43742204e-02,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       3.96292063e-01, 3.30281055e-01, 1.73992026e-01, 1.19261113e-02,
       3.71950067e-02, 8.98952507e-03, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00])

and we can further check our solution if it satisfies all constraints:

In [92]: np.all((p_sol >= 0) & (p_sol <= 1))
Out[92]: True

In [93]: np.sum(p_sol, axis=0)
Out[93]: array([1., 1., 1.])

In [94]: np.all(np.sum(p_sol, axis=1) <= 1)
Out[94]: True

Thanks - this code runs without error, and gives a solution that that meets all the constraints. However it doesn't appear to do a great job of minimising the cost function. In most cases I expect the cost to be close to zero and rows where y=0, should have p[1], p[2] and p[3] =0 as well. I have tried playing with the tolerance argument but without success. Perhaps the problem isn't suitable for the scipy solver? — user3725021, Sep 23 '19 at 00:37
I think the initial guess is crucial. You could set the specific rows of p0 to 0.0 and further decrease the solver's tolerance. Then it might be worthwhile to increase the iteration limit also. With `tol = 1.0e-12, options={'maxiter' : 10000}` and `p0[-3:,:]=0.0` I get a better solution with function objective 1.46e-7. — joni, Sep 23 '19 at 06:48

score 1 · Answer 2 · answered Sep 19 '19 at 17:02

As @sascha mentioned, scipy.minimize() only solves for 1-d arrays. Therefore, you need to flatten your matrix (convert it to a vector) and after the minimization is done, you reconstruct your matrix. This can be done in various ways. In the following code I used numpy.vstack() and numpy.concatenate() to do the transformations. I also implemented and joined your constraints. See code in the following:

# import packages
import numpy as np
from scipy.optimize import minimize

# create objective function
def objective(p, y, s):
    x   = p[0:sep] + 2*p[sep:2*sep] + 3*p[2*sep:3*sep]
    return np.sum(np.abs(y-x))

def vec2mat(vec):
    sep = int(vec.size / 3) 
    return np.vstack((np.vstack((vec[0:sep], vec[sep:2*sep])), vec[2*sep:3*sep]))

def mat2vec(mat):
    return np.concatenate((np.concatenate((mat[:,0], mat[:,1]), axis=0), mat[:,2]), axis=0)

def columns_sum(mat):
    return np.array([sum(x) for x in zip(*mat)])

def rows_sum(mat):    
    return np.array([sum(x) for x in mat])
def constraint_func(x):
    c1 = np.all(0 <= x) & np.all(x <= 1)         # 0 <= x_i <= 1
    c2 = np.all(rows_sum(vec2mat(x)) <= 1)       # row sum <= 1
    c3 = np.all(1 == columns_sum(vec2mat(x)))    # columns sum == 1
    if (c1 and c2 and c3) == True: return 1
    else                         : return 0

# define constraints
cons = [{'type':'eq', 'fun': lambda x: constraint_func(x) - 1 }]  

# generate target data (y)
y = np.array([2.48, 1.75, 1.13, 0.20, 0.20, 0.20, 0.0, 0.0, 0.0, 0.0])

# initialise starting param values
p0      = np.zeros((10,3))
p0[:,0] = 0.001
p0[:,1] = 0.002 
p0[:,2] = 0.003

# flatten p
p0_flat = mat2vec(p0)
sep     = int(p0_flat.size / 3) 

# solve
sol = minimize(objective, p0_flat, args=(y,sep), method='SLSQP', constraints=cons)
print(sol)

# reconstruct p
p_sol = vec2mat(sol.x)

Unfortunately, though the code compiles and works correctly but it returns the following output:

     fun: 5.932
     jac: array([-1., -1., -1., -1., -1., -1.,  1.,  1.,  1.,  1., -2., -2., -2.,
       -2., -2., -2.,  2.,  2.,  2.,  2., -3., -3., -3., -3., -3., -3.,
        3.,  3.,  3.,  3.])
 message: 'Singular matrix C in LSQ subproblem'
    nfev: 32
     nit: 1
    njev: 1
  status: 6
 success: False
       x: array([0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001,
       0.001, 0.002, 0.002, 0.002, 0.002, 0.002, 0.002, 0.002, 0.002,
       0.002, 0.002, 0.003, 0.003, 0.003, 0.003, 0.003, 0.003, 0.003,
       0.003, 0.003, 0.003])

Therefore, you might need to further debug the code (I suspect the constraints section), or try different methods offered in scipy.minimize_docs.

Solve parameter matrix with element-wise constraints using scipy optimize

2 Answers2