使用 scipy 优化求解具有逐元素约束的参数矩阵

Question

我有一个相对简单的优化问题，但我不是特别精通 scipy 并且不知道如何应用必要的约束。我的 objective 函数是最小化 10 元素向量 y 和 x 之间的绝对距离，其中 x 是 10x3 参数矩阵的加权行和 p

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

我需要在参数矩阵p中添加以下约束：

p的每个元素，但>=0且<=1
p每一列的总和必须为1
每行p的总和不得超过1

我曾尝试根据其他 SO 问题实施约束，但我承认我并不完全理解它，也不完全理解它产生的错误。

# 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)

根据运行此代码，我收到以下错误消息：

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

非常感谢任何帮助。

Answer 1

如 @sascha 所述，scipy.minimize() 仅求解一维数组。因此，您需要展平矩阵（将其转换为向量），并在完成最小化后重建矩阵。这可以通过多种方式完成。在下面的代码中，我使用 numpy.vstack() 和 numpy.concatenate() 进行转换。我也实施并加入了你的约束。请参阅以下代码：

# 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)

不幸的是，虽然代码可以正确编译和运行，但是 returns 输出如下：

     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])

因此，您可能需要进一步调试代码（我怀疑约束部分），或尝试 scipy.minimize_docs 中提供的不同方法。

Answer 2

此答案基于 SuperKogito 的回答：

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)

给我

     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])

我们可以进一步检查我们的解决方案是否满足所有约束：

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

使用 scipy 优化求解具有逐元素约束的参数矩阵

Solve parameter matrix with element-wise constraints using scipy optimize

optimization

scipy

scipy-optimize

scipy-optimize-minimize