constrained linear regression / quadratic programming python
<p>I have a dataset like this:</p>

<pre><code>import numpy as np

a = np.array([1.2, 2.3, 4.2])
b = np.array([1, 5, 6])
c = np.array([5.4, 6.2, 1.9])

m = np.vstack([a,b,c])
y = np.array([5.3, 0.9, 5.6])
</code></pre>

<p>and want to fit a constrained linear regression</p>

<p>y = b1*a + b2*b + b3*c</p>

<p>where all b's sum to one and are positive: b1+b2+b3=1</p>

<p>A similar problem in R is specified here:</p>

<p><a href="http://stats.stackexchange.com/questions/21565/how-do-i-fit-a-constrained-regression-in-r-so-that-coefficients-total-1">http://stats.stackexchange.com/questions/21565/how-do-i-fit-a-constrained-regression-in-r-so-that-coefficients-total-1</a></p>

<p>How can I do this in python?</p>

<p>You can get a good solution to this with a little bit of math and <code>scipy.optimize.nnls</code>:</p>

<p>First we do the math:</p>

<p>If </p>

<p>y = b1*a + b2*b + b3*c and b1 + b2 + b3 = 1, then b3 = 1 - b1 - b2. </p>

<p>If we substitute and simplify we end up with </p>

<p>y - c = b1(a - c) + b2(b - c)</p>

<p>Now, we don't have any equality constraints and nnls can solve directly:</p>

<pre><code>import scipy.optimize
A = np.vstack([a - c, b - c]).T
(b1, b2), norm = scipy.optimize.nnls(A, y - c)
b3 = 1 - b1 - b2
</code></pre>

<p>This recovers the same solution as obtained in the other answer using cvxpy.</p>

<pre><code>b1 = 0.77608809648662802
b2 = 0.0
b3 = 0.22391190351337198
norm = 4.337947941595865
</code></pre>

<p>This approach can be generalised to an arbitrary number of dimensions as follows. Assume that we have a matrix B constructed with a, b, c from the original question arranged in the columns. Any additional dimensions will get added to this.</p>

<p>Now, we can do</p>

<pre><code>A = B[:, :-1] - B[:, -1:]
bb, norm = scipy.optimize.nnls(A, y - B[:, -1])
bi = np.append(bb, 1 - sum(bb))
</code></pre>

python scipy linear-regression quadratic-programming