spore234 - 10 months ago 101

Python Question

I have a dataset like this:

`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])

and want to fit a constrained linear regression

y = b1*a + b2*b + b3*c

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

A similar problem in R is specified here:

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

How can I do this in python?

Answer Source

```
import numpy as np
from scipy.optimize import minimize
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])
def loss(x):
return np.sum(np.square((np.dot(x, m) - y)))
cons = ({'type': 'eq',
'fun' : lambda x: np.sum(x) - 1.0})
x0 = np.zeros(m.shape[1])
res = minimize(loss, x0, method='SLSQP', constraints=cons,
bounds=[(0, np.inf) for i in range(m.shape[1])], options={'disp': True})
print(res.x)
print(np.dot(res.x, m.T))
print(np.sum(np.square(np.dot(res.x, m) - y)))
```

```
Optimization terminated successfully. (Exit mode 0)
Current function value: 18.817792344
Iterations: 5
Function evaluations: 26
Gradient evaluations: 5
[ 0.7760881 0. 0.2239119]
[ 1.87173571 2.11955951 4.61630834]
18.817792344
```

- It's obvious that the model-capabilitites / model-complexity is not enough to obtain a good performance (high loss!)

Advantages:

- cvxpy prooves, that the problem is convex
- convergence for convex optimization problems is guaranteed (might be also true for the above)
- in general: better accuracy
- in general: more robust in regards to numerical instabilities (solver is only able to solve SOCPs; not non-convex models like the SLSQP-approach above!)

```
import numpy as np
from cvxpy import *
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])
X = Variable(m.shape[1])
constraints = [X >= 0, sum_entries(X) == 1.0]
product = m.T * diag(X)
diff = sum_entries(product, axis=1) - y
problem = Problem(Minimize(norm(diff)), constraints)
problem.solve(verbose=True)
print(problem.value)
print(X.value)
```

```
ECOS 2.0.4 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS
It pcost dcost gap pres dres k/t mu step sigma IR | BT
0 +0.000e+00 -0.000e+00 +2e+01 5e-01 1e-01 1e+00 4e+00 --- --- 1 1 - | - -
1 +2.451e+00 +2.539e+00 +4e+00 1e-01 2e-02 2e-01 8e-01 0.8419 4e-02 2 2 2 | 0 0
2 +4.301e+00 +4.306e+00 +2e-01 5e-03 7e-04 1e-02 4e-02 0.9619 1e-02 2 2 2 | 0 0
3 +4.333e+00 +4.334e+00 +2e-02 4e-04 6e-05 1e-03 4e-03 0.9326 2e-02 2 1 2 | 0 0
4 +4.338e+00 +4.338e+00 +5e-04 1e-05 2e-06 4e-05 1e-04 0.9698 1e-04 2 1 1 | 0 0
5 +4.338e+00 +4.338e+00 +3e-05 8e-07 1e-07 3e-06 7e-06 0.9402 7e-03 2 1 1 | 0 0
6 +4.338e+00 +4.338e+00 +7e-07 2e-08 2e-09 6e-08 2e-07 0.9796 1e-03 2 1 1 | 0 0
7 +4.338e+00 +4.338e+00 +1e-07 3e-09 4e-10 1e-08 3e-08 0.8458 2e-02 2 1 1 | 0 0
8 +4.338e+00 +4.338e+00 +7e-09 2e-10 2e-11 9e-10 2e-09 0.9839 5e-02 1 1 1 | 0 0
OPTIMAL (within feastol=1.7e-10, reltol=1.5e-09, abstol=6.5e-09).
Runtime: 0.000555 seconds.
4.337947939 # needs to be squared to be compared to scipy's output!
# as we are using l2-norm (outer sqrt) instead of sum-of-squares
# which is nicely converted to SOCP-form and easier to
# tackle by SOCP-based solvers like ECOS
# -> does not change the solution-vector x, only the obj-value
[[ 7.76094262e-01]
[ 7.39698388e-10]
[ 2.23905737e-01]]
```