Solving system using linalg with constraints

Question

I want to solve some system in the form of matrices using linalg, but the resulting solutions should sum up to 1. For example, suppose there are 3 unknowns, x, y, z. After solving the system their values should sum up to 1, like .3, .5, .2. Can anyone please tell me how I can do that?

Currently, I am using something like result = linalg.solve(A, B), where A and B are matrices. But this doesn't return solutions in the range [0, 1].

Answer 1

Per the docs,

linalg.solve is used to compute the "exact" solution, x, of the well-determined, i.e., full rank, linear matrix equation ax = b.

Being linear, there can be at most one solution. If the solution you found does not sum up to 1, then adding the extra constraint would yield no solution.

However, you could use scipy.optimize.minimize to find the point on the constraint plane which minimizes the quantity ||Ax-b||^2:

def f(x):
    y = np.dot(A, x) - b
    return np.dot(y, y)

cons = ({'type': 'eq', 'fun': lambda x: x.sum() - 1})
res = optimize.minimize(f, [0, 0, 0], method='SLSQP', constraints=cons, 
                        options={'disp': False})

---

For example, given this system of equations

import numpy as np
import numpy.linalg as LA
import scipy.optimize as optimize

A = np.array([[1, 3, 4], [5, 6, 9], [1, 2, 3]])
b = np.array([1, 2, 1])
x = LA.solve(A, b)

The solution does not add up to 1:

print(x)
# [-0.5 -1.5  1.5]

But you could try to minimize f:

def f(x):
    y = np.dot(A, x) - b
    return np.dot(y, y)

subject to the constraint cons:

cons = ({'type': 'eq', 'fun': lambda x: x.sum() - 1})
res = optimize.minimize(f, [0, 0, 0], method='SLSQP', constraints=cons, 
                        options={'disp': False})
xbest = res['x']
# array([ 0.30000717,  1.89998823, -1.1999954 ])

xbest sums to 1:

print(xbest.sum())
1

The difference A·xbest - b is:

print(np.dot(A, xbest) - b)
# [ 0.19999026  0.10000663 -0.50000257]

and the sum of the squares of the difference, (also computable as f(xbest)) is :

print(res['fun'])
0.30000000014542572

No other value of x minimizes this quantity more while satisfying the constraint.

Answer 2

You could add a row consisting of ones to A and add one to B. After that use result = linalg.lstsq(A, B)[0]

Or you can replace one of A's rows to row consisting of ones, also replace value in B to one in the same row. Then use result = linalg.solve(A, B)

Answer 3

As an alternative you could use numpy.linalg.lstsq combined with the new row to your A matrix: [1, 1, 1] and B matrix [1].

import numpy as np
A = np.array([[1, 0, 1],
              [0, 1, 1],
              [1, 1, 0],
              [1, 1, 1]]
B = np.array([[0.5, 0.7, 0.8, 1.0]]).T
x, _, _, _ = np.linalg.lstsq(A, B)

This gives x as [[0.3], [0.5], [0.2]].