Numpy: Linear system with specific conditions. No negative solutions

Question

I'm writing a Python code using numpy. In my code I use "linalg.solve" to solve a linear system of n equations in n variables. Of course the solutions could be either positive or negative. What I need to do is to have always positive solutions or at least equal to 0. To do so I first want the software to solve my linear system of equations in this form

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

in which x is an array with n variables in a specific order (x1, x2, x3.....xn), A is a n dimensional square matrix and b is a n-dimensional array. Now I thought to do this:

-solve the system of equations

-check if every x is positive

-if not, every negative x I'll want them to be =0 (for example x2=-2 ---->x2=0)

-with a generic xn=0 want to eliminate the n-row and the n-coloumn in the n dimensional square matrix A (I'll obtain another square matrix A1) and eliminate the n element in b obtaining b1.

-solve the system again with the matrix A1 and b1

-re-iterate untill every x is positive or zero

-at last build a final array of n elements in which I'll put the last iteration solutions and every variable which was equal to zero ( I NEED THEM IN ORDER AS IT WOULD HAVE BEEN NO ITERATIONS so if during the iterations it was x2=0 -----> xfinal=[x1, 0 , x3,.....,xn]

Think it 'll work but don't know how to do it in python.

Hope I was clear. Can't really figure it out!

sebacastroh sebacastroh · Accepted Answer · 2016-05-01T17:39:23

You have a minimization problem, i.e.

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

You can use the Optimize module from 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]

With your method I get x = [ 0.27272727, 0., 1.90909091].

In the case you still want to use your algorithm, it is below

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]

But I don't recommend you to use it, you should use the minimize option.

Numpy: Linear system with specific conditions. No negative solutions

1 Answers