Numpy: Линейная система с особыми условиями. Нет отрицательных решений

Question

Я пишу код Python, используя NumPy. В моем коде я использую "linalg.solve" для решения линейной системы из n уравнений с n переменными. Конечно, решения могут быть как положительными, так и отрицательными. Мне нужно иметь всегда положительные решения или, по крайней мере, равные 0. Для этого я сначала хочу, чтобы программное обеспечение решало мою линейную систему уравнений в таком виде

x=np.linalg.solve(A,b)

в которой x является массивом с n переменными в определенном порядке (x1, x2, x3.....xn), A является размерной квадратной матрицей, а b является n-мерным массивом. Теперь я подумал сделать это:

решить систему уравнений

-проверьте, если каждый х положительный

-Если нет, каждый минус х я хочу, чтобы они были =0 (например, х2=-2 ----> х2 = 0)

-с помощью общего xn=0 хотите исключить n-строку и n-столбец в n-мерной квадратной матрице A (я получу другую квадратную матрицу A1) и исключить n-элемент в b, получив b1.

- Решите систему снова с матрицами A1 и b1.

повторять до тех пор, пока каждый х не будет положительным или нулевым

-при последней сборке окончательный массив из n элементов, в который я поместил бы последние итерационные решения и каждую переменную, равную нулю (ИМ НУЖНО, ЧТОБЫ БЫ НЕ БЫЛО Итераций, так что если во время итераций это было x2=0 -----> xfinal=[x1, 0, x3,.....,xn]

Думаю, это сработает, но не знаю, как это сделать на python.

Надеюсь, мне было ясно. Не могу понять это!

Answer 1

У вас есть проблема минимизации, т.е.

min ||Ax - b||
s.t. x_i >= 0 for all i  in [0, n-1]

Вы можете использовать модуль Optimize от Scipy

import numpy as np
from scipy.optimize import minimize

A = np.array([[1., 2., 3.],[4., 5., 6.],[7., 8., 10.]], order='C')
b = np.array([6., 12., 21.])
n = len(b)

# Ax = b --> x = [1., -2., 3.]

fun = lambda x: np.linalg.norm(np.dot(A,x)-b)
# xo = np.linalg.solve(A,b)
# sol = minimize(fun, xo, method='SLSQP', constraints={'type': 'ineq', 'fun': lambda x:  x})
sol = minimize(fun, np.zeros(n), method='L-BFGS-B', bounds=[(0.,None) for x in xrange(n)])

x = sol['x'] # [2.79149722e-01, 1.02818379e-15, 1.88222298e+00]

С твоим методом я получаю x = [ 0.27272727, 0., 1.90909091],

Если вы все еще хотите использовать свой алгоритм, он ниже

n = len(b)
x = np.linalg.solve(A,b)
pos = np.where(x>=0.)[0]

while len(pos) < n:
    Ap = A[pos][:,pos]
    bp = b[pos]
    xp = np.linalg.solve(Ap, bp)
    x = np.zeros(len(b))
    x[pos] = xp
    pos = np.where(x>=0.)[0]

Но я не рекомендую вам использовать его, вы должны использовать опцию минимизации.

Answer 2

Еще быстрее и надежнее также с использованием минимизации является методscipy.optimize.lsq_linearкоторый специально предназначен для линейной оптимизации.

Границы транспонируются и используются np.infвместо Noneдля верхней границы.

Рабочий пример:

      from scipy.optimize import lsq_linear
n = A.shape[1]
res = lsq_linear(A, b, bounds=np.array([(0.,np.inf) for i in range(n)]).T, lsmr_tol='auto', verbose=1)
y = res.x

Вы предоставляете матрицу Aи вектор bкоторые имеют одинаковое количество строк (= размерность диапазона A).