Решите матрицу параметров с поэлементными ограничениями, используя scipy optimize

Question

У меня относительно простая проблема оптимизации, но я не особо разбираюсь в scipy и не могу понять, как применить необходимые ограничения.Моя цель - минимизировать абсолютное расстояние между вектором из 10 элементов y и x , где x - взвешенная сумма строк матрицы параметров 10x3 p

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

Мне нужно добавить следующие ограничения в матрицу параметров p :

1. Каждый элемент p , но be> = 0 и <= 1 </li>
2. Сумма каждого столбца p должна составлять 1
3. Сумма каждой строки 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: ось 1 выходит за границы массива измерения 1

Любая помощь очень ценится.

Answer 1

Этот ответ основан на ответе 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

Answer 2

Как упоминалось @ sascha , scipy.minimize() разрешается только для 1-го массива.Следовательно, вам нужно выровнять свою матрицу (преобразовать ее в вектор), и после того, как минимизация будет выполнена, вы восстановите свою матрицу.Это можно сделать разными способами.В следующем коде я использовал 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)

К сожалению, хотя код компилируется и работает правильно, но возвращает следующий вывод:

     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 .