# Optimization in Python

You might have noticed that we didn't do anything related to sparsity with scikit-learn models. A lot of the work we covered in the machine learning class is very recent research, and as such is typically not implemented by the popular libraries.

If we want to do things like sparse regression, we're going to have to roll up our sleeves and do it ourselves. For that, we need to be able to solve optimization problems. In Julia, we did this with JuMP. In Python, we'll use a similar library called *pyomo*.

# Installing pyomo

You can run the following command to install pyomo if you haven't already.


```python
!pip install pyomo --user
```

# Intro to pyomo

Let's see how we translate a simple, 2 variable LP to pyomo code.

$$
\begin{align*}
\max_{x,y} \quad& x + 2y \\
\text{s.t.}\quad& x + y \leq 1 \\
& x, y \geq 0.
\end{align*}
$$

First thing is to import the pyomo functions:


```python
from pyomo.environ import *
from pyomo.opt import SolverFactory
```

Next, we construct a model object. This is a container for everything in our optimization problem: variables, constraints, solver options, etc.


```python
m = ConcreteModel()
```

Next, we define the two decision variables in our optimization problem. We use the ``Var`` function to create the variables. The `within` keyword is used to specify the bounds on the variables, or equivalently the `bounds` keyword. The variables are added to the model object with names `x` and `y`.


```python
m.x = Var(within=NonNegativeReals)
m.y = Var(bounds=(0, float('inf')))
```

We now add the single constraint of our problem using the ``Constraint`` function. We write it algebraically, and save the result to the model.


```python
m.con = Constraint(expr=m.x + m.y <= 1)
```

We specify the objective function with the `Objective` function:


```python
m.obj = Objective(sense=maximize, expr=m.x + 2 * m.y)
```

We solve the optimization problem by first specifying a solver using `SolverFactory` and then using this solver to solve the model:


```python
solver = SolverFactory('gurobi')
solver.solve(m)
```

We can now inspect the solution values and optimal cost.


```python
m.obj()
```


```python
m.x.value
```


```python
m.y.value
```

Let's put it all together to compare with Julia/JuMP


```python
# Create model
m = ConcreteModel()
# Add variables
m.x = Var(within=NonNegativeReals)
m.y = Var(bounds=(0, float('inf')))
# Add constraint
m.con = Constraint(expr=m.x + m.y <= 1)
# Add objective
m.obj = Objective(sense=maximize, expr=m.x + 2 * m.y)
# Solve model
solver = SolverFactory('gurobi')
solver.solve(m)
# Inspect solution
print(m.obj())
print(m.x.value)
print(m.y.value)
```

```julia
# Create model
m = Model(solver=GurobiSolver())
# Add variables
@variable(m, x >= 0)
@variable(m, y >= 0)
# Add constraint
@constraint(m, x + y <= 1)
# Add objective
@objective(m, Max, x + 2y)
# Solve model
solve(m)
# Inspect solution
@show getobjectivevalue(m)
@show getvalue(x)
@show getvalue(y)
```

### Exercise

Code and solve the following optimization problem:

$$
\begin{align*}
\min_{x,y} \quad& 3x - y \\
\text{s.t.}\quad& x + 2y \geq 1 \\
& x \geq 0 \\
& 0 \leq y \leq 1.
\end{align*}
$$


```python
# Create the model
m = ConcreteModel()
# Add the variables
m.x = Var(within=NonNegativeReals)
m.y = Var(bounds=(0, 1))
# Add the constraint
m.con = Constraint(expr=m.x + 2 * m.y >= 1)
# Add the objective
m.obj = Objective(sense=minimize, expr=3 * m.x - m.y)

solver = SolverFactory('gurobi')
solver.solve(m)

print(m.x.value, m.y.value)
```


```python
for v in m.component_data_objects(Var, active=True):
    print(v, value(v)) # doctest: +SKIP
```


```python
m.pprint()
```

# Index sets

Let's now move to a more complicated problem. We'll look at a transportation problem:

$$
\begin{align}
\min & \sum\limits_{i = 1}^{m} \sum\limits_{j = 1}^{n} c_{ij} x_{ij}\\
& \sum\limits_{j = 1}^{n} x_{ij} \leq b_i && i = 1, \ldots, m\\
& \sum\limits_{i = 1}^{m} x_{ij} = d_j && j = 1, \ldots, n\\
& x_{ij} \ge 0 && i = 1, \ldots, m, j = 1, \ldots, n
\end{align}
$$

And with some data:


```python
import numpy as np

m = 2  # Number of supply nodes
n = 5  # Number of demand nodes
# Supplies
b = np.array([1000, 4000])
# Demands
d = np.array([500, 900, 1800, 200, 700])
# Costs
c = np.array([[2, 4, 5, 2, 1], 
              [3, 1, 3, 2, 3]])
```

Now we can formulate the problem with pyomo


```python
model = ConcreteModel()
```

First step is adding variables. We can add variables with indices by passing the relevant index sets to the `Var` constructor. In this case, we need a $m$-by$n$ matrix of variables:


```python
model.x = Var(range(m), range(n), within=NonNegativeReals)
```

Now to add the constraints. We have to add one supply constraint for each factory, so we might try something like:


```python
for i in range(m):
    model.supply = Constraint(expr=sum(model.x[i, j] for j in range(n)) <= b[i])
```

Can you see the problem? We are overwriting `model.supply` in each iteration of the loop, and so only the last constraint is applied.

Luckily, pyomo has a (not-so-easy) way to add multiple constraints at a time. We first define a *rule* that takes in the model and any required indices, and then returns the expression for the constraint:


```python
def supply_rule(model, i):
    return sum(model.x[i, j] for j in range(n)) <= b[i]
```

We then add the constraint by referencing this rule along with the index set we want the constraint to be defined over:


```python
model.supply2 = Constraint(range(m), rule=supply_rule)
```

We then apply the same approach for the demand constraints


```python
def demand_rule(model, j):
    return sum(model.x[i, j] for i in range(m)) == d[j]
model.demand = Constraint(range(n), rule=demand_rule)
```

Finally, we add the objective:


```python
model.obj = Objective(sense=minimize, 
    expr=sum(c[i, j] * model.x[i, j] 
    for i in range(m) for j in range(n)))
```

Now we can solve the problem


```python
solver = SolverFactory('gurobi')
solver.solve(model)
```

It solved, so we can extract the results


```python
flows = np.array([[model.x[i, j].value for j in range(n)] for i in range(m)])
flows
```

We can also check the objective value for the cost of this flow


```python
model.obj()
```

For simplicity, here is the entire formulation and solving code together:


```python
model = ConcreteModel()
# Variables
model.x = Var(range(m), range(n), within=NonNegativeReals)
# Supply constraint
def supply_rule(model, i):
    return sum(model.x[i, j] for j in range(n)) <= b[i]
model.supply2 = Constraint(range(m), rule=supply_rule)
# Demand constraint
def demand_rule(model, j):
    return sum(model.x[i, j] for i in range(m)) == d[j]
model.demand = Constraint(range(n), rule=demand_rule)
# Objective
model.obj = Objective(sense=minimize, 
    expr=sum(c[i, j] * model.x[i, j] 
    for i in range(m) for j in range(n)))
# Solve
solver = SolverFactory('gurobi')
solver.solve(model)
# Get results
flows = np.array([[model.x[i, j].value for j in range(n)] for i in range(m)])
model.obj()
```

# Machine Learning

Now let's put our pyomo knowledge to use and implement some of the same methods we saw in the machine learning class

First, specify your solver executable location:


```python
executable='C:/Users/omars/.julia/v0.6/Ipopt/deps/usr/bin/ipopt.exe'
```

To use the version left over from Julia
### On MacOS and Linux

`executable="~/.julia/v0.6/Homebrew/deps/usr/Cellar/ipopt/3.12.4_1/bin/ipopt")`

### On Windows

The path is probably under WinRPM:

`executable='%HOME%/.julia/v0.6/WinRPM/...')")`


# Linear Regression

Let's just try a simple linear regression


```python
def linear_regression(X, y):
    n, p = X.shape

    # Create model
    m = ConcreteModel()

    # Add variables
    m.beta = Var(range(p))

    # Add constraints

    # Add objective
    m.obj = Objective(sense=minimize, expr=sum(
        pow(y[i] - sum(X[i, j] * m.beta[j] for j in range(p)), 2) 
        for i in range(n)))

    solver = SolverFactory('ipopt', executable=executable)
    
    ## tee=True enables solver output
    # results = solver.solve(m, tee=True)
    results = solver.solve(m, tee=False)

    return [m.beta[j].value for j in range(p)]
```

Let's load up some data to test it out on:


```python
from sklearn.datasets import load_boston
data = load_boston()
X = data.data
y = data.target
```

Try our linear regression function:


```python
print(linear_regression(X, y))
```

We can compare with sklearn to make sure it's right:


```python
from sklearn.linear_model import LinearRegression
m = LinearRegression(fit_intercept=False)
m.fit(X, y)
m.coef_
```

Just for reference, let's look back at how we do the same thing in JuMP!

```julia
using JuMP, Gurobi
function linear_regression(X, y)
    n, p = size(X)
    m = Model(solver=GurobiSolver())
    @variable(m, beta[1:p])
    @objective(m, Min, sum((y[i] - sum(X[i, j] * beta[j] for j = 1:p)) ^ 2 for i = 1:n))
    solve(m)
    getvalue(beta)
end
```

or even

```julia
using JuMP, Gurobi
function linear_regression(X, y)
    n, p = size(X)
    m = Model(solver=GurobiSolver())
    @variable(m, beta[1:p])
    @objective(m, Min, sum((y - X * beta) .^ 2))
    solve(m)
    getvalue(beta)
end
```

Much simpler!

### Exercise

Modify the linear regression formulation to include an intercept term, and compare to scikit-learn's LinearRegression with `fit_intercept=False` to make sure it's the same


```python
def linear_regression_intercept(X, y):
    n, p = X.shape

    # Create model
    m = ConcreteModel()

    # Add variables
    m.beta = Var(range(p))
    m.b0 = Var()

    # Add constraints

    # Add objective
    m.obj = Objective(sense=minimize, expr=sum(
        pow(y[i] - sum(X[i, j] * m.beta[j] for j in range(p)) - m.b0, 2) 
        for i in range(n)))

    solver = SolverFactory('ipopt', executable=executable)
    
    ## tee=True enables solver output
    # results = solver.solve(m, tee=True)
    results = solver.solve(m, tee=False)

    return [m.beta[j].value for j in range(p)]

linear_regression_intercept(X, y)
```


```python
m = LinearRegression(fit_intercept=True)
m.fit(X, y)
m.coef_
```

# Robust Regression

We saw in the class that both ridge and lasso regression were robust versions of linear regression. Both of these are provided by `sklearn`, but we need to know how to implement them if we want to extend regression ourselves


```python
def ridge_regression(X, y, rho):
    n, p = X.shape

    # Create model
    m = ConcreteModel()

    # Add variables
    m.beta = Var(range(p))

    # Add objective
    m.obj = Objective(sense=minimize, expr=sum(
        pow(y[i] - sum(X[i, j] * m.beta[j] for j in range(p)),2) 
        for i in range(n)) + rho * sum(pow(m.beta[j], 2) for j in range(p)))

    solver = SolverFactory('ipopt', executable=executable)
    
    ## tee=True enables solver output
    # results = solver.solve(m, tee=True)
    results = solver.solve(m, tee=False)
    return [m.beta[j].value for j in range(p)]
```


```python
ridge_regression(X, y, 100000)
```


```python
def lasso(X, y, rho):
    n, p = X.shape

    # Create model
    m = ConcreteModel()

    # Add variables
    m.beta = Var(range(p))

    # Add objective
    m.obj = Objective(sense=minimize, expr=sum(
        pow(y[i] - sum(X[i, j] * m.beta[j] for j in range(p)),2) 
        for i in range(n)) + rho * sum(pow(m.beta[j], 2) for j in range(p)))

    solver = SolverFactory('ipopt', executable=executable)
    
    ## tee=True enables solver output
    # results = solver.solve(m, tee=True)
    results = solver.solve(m, tee=False)
    return [m.beta[j].value for j in range(p)]
```

### Exercise

Implement Lasso regression


```python
def lasso_regression(X, y, rho):
    n, p = X.shape

    # Create model
    m = ConcreteModel()

    # Add variables
    m.beta = Var(range(p))
    m.absb = Var(range(p))

    # Add constraints
    def absbeta1(m, j):
        return m.beta[j] <= m.absb[j]
    m.absb1 = Constraint(range(p), rule=absbeta1)
    def absbeta2(m, j):
        return -m.beta[j] <= m.absb[j]
    m.absb2 = Constraint(range(p), rule=absbeta2)
        

    # Add objective
    m.obj = Objective(sense=minimize, expr=sum(
        pow(y[i] - sum(X[i, j] * m.beta[j] for j in range(p)), 2) 
        for i in range(n)) + rho * sum(m.absb[j] for j in range(p)))

    solver = SolverFactory('ipopt', executable=executable)
    
    ## tee=True enables solver output
    # results = solver.solve(m, tee=True)
    results = solver.solve(m, tee=False)
    return [m.beta[j].value for j in range(p)]
```


```python
lasso_regression(X, y, 1000)
```

# Sparse Regression


```python
def sparse_regression(X, y, k):
    n, p = X.shape
    M = 1000

    # Create model
    m = ConcreteModel()

    # Add variables
    m.beta = Var(range(p))
    m.z = Var(range(p), within=Binary)

    # Add constraints
    def bigm1(m, j):
        return m.beta[j] <= M * m.z[j]
    m.bigm1 = Constraint(range(p), rule=bigm1)
    def bigm2(m, j):
        return m.beta[j] >= -M * m.z[j]
    m.bigm2 = Constraint(range(p), rule=bigm2)
        
    m.sparsity = Constraint(expr=sum(m.z[j] for j in range(p)) <= k)

    # Add objective
    m.obj = Objective(sense=minimize, expr=sum(
        pow(y[i] - sum(X[i, j] * m.beta[j] for j in range(p)), 2) 
        for i in range(n)))

    solver = SolverFactory('ipopt', executable=executable)
    
    ## tee=True enables solver output
    # results = solver.solve(m, tee=True)
    results = solver.solve(m, tee=False)
    return [m.beta[j].value for j in range(p)]
```


```python
sparse_regression(X, y, 10)
```


```python
import numpy as np
l = np.array([1,2,3,4])

print(l**2)
print([sqrt(i) for i in l])
```

### Exercise

Try implementing the algorithmic framework for linear regression:
- sparsity constraints
- lasso regularization
- restrict highly correlated pairs of features
- nonlinear transformations (just $\sqrt(x)$ and $x^2$)


```python
import numpy as np
from sklearn.preprocessing import normalize

def all_regression(X_orig, y, k, rho):
    n, p_orig = X_orig.shape
    M = 10
    
    X = np.concatenate(
        [X_orig, np.sqrt(X_orig), np.square(X_orig)], axis=1
    )
    p = X.shape[1]
    
    # Normalize data
    X = normalize(X, axis=0)
    y = (y - np.mean(y)) / np.linalg.norm(y)

    # Create model
    m = ConcreteModel()

    # Add variables
    m.beta = Var(range(p))
    m.z = Var(range(p), within=Binary)
    m.absb = Var(range(p))

    # Sparsity constraints
    def bigm1(m, j):
        return m.beta[j] <= M * m.z[j]
    m.bigm1 = Constraint(range(p), rule=bigm1)
    def bigm2(m, j):
        return m.beta[j] >= -M * m.z[j]
    m.bigm2 = Constraint(range(p), rule=bigm2)
    m.sparsity = Constraint(expr=sum(m.z[j] for j in range(p)) <= k)
    
    # Lasso constraints
    def absbeta1(m, j):
        return m.beta[j] <= m.absb[j]
    m.absb1 = Constraint(range(p), rule=absbeta1)
    def absbeta2(m, j):
        return -m.beta[j] <= m.absb[j]
    m.absb2 = Constraint(range(p), rule=absbeta2)
    
    # Correlation constraints
    corX = np.corrcoef(np.transpose(X))
    def cor_rule(m, i, j):
        if i > j and abs(corX[i, j]) > 0.8:
            return (sum(m.z[k] for k in range(i, p, p_orig)) + 
                    sum(m.z[k] for k in range(j, p, p_orig)) <= 1)
        else:
            return Constraint.Skip
    m.cor = Constraint(range(p_orig), range(p_orig), rule=cor_rule)
    
    # Nonlinear constraints
    def nl_rule(m, i):
        return sum(m.z[k] for k in range(i, p, p_orig)) <= 1
    m.nl = Constraint(range(p_orig), rule=nl_rule)

    # Add objective
    m.obj = Objective(sense=minimize, expr=sum(
        pow(y[i] - sum(X[i, j] * m.beta[j] for j in range(p)), 2) 
        for i in range(n)) + rho * sum(m.absb[j] for j in range(p)))

    solver = SolverFactory('ipopt', executable=executable)
    
    ## tee=True enables solver output
#     results = solver.solve(m, tee=True)
    results = solver.solve(m, tee=False)

    return np.array([m.beta[j].value for j in range(p)]).reshape(-1, p_orig)
```


```python
all_regression(X, y, 6, 0)
```

# Logistic Regression

Like JuMP, we need to use a new solver for the nonlinear problem. We can use Ipopt as before, except we have to set it up manually. You'll need to download Ipopt and add it to the PATH. 

On Mac, you can do this with Homebrew if you have it:

The other way is to download a copy of ipopt and specify the path to it exactly when creating the solver. For example, I have a copy of Ipopt left over from JuMP, which I can use by modifying the SolverFactory line as indicated below:


```python
def logistic_regression(X, y):
    n, p = X.shape
    
    # Convert y to (-1, +1)
    assert np.min(y) == 0
    assert np.max(y) == 1
    Y = y * 2 - 1
    assert np.min(Y) == -1
    assert np.max(Y) == 1

    # Create the model
    m = ConcreteModel()

    # Add variables
    m.b = Var(range(p))
    m.b0 = Var()

    # Set nonlinear objective function
    m.obj = Objective(sense=maximize, expr=-sum(
        log(1 + exp(-Y[i] * (sum(X[i, j] * m.b[j] for j in range(p)) + m.b0)))
        for i in range(n)))

    # Solve the model and get the optimal solutions
    solver = SolverFactory('ipopt', executable=executable)
        
    solver.solve(m)
    return [m.b[j].value for j in range(p)], m.b0.value
```

Load up some data


```python
from sklearn.datasets import load_breast_cancer
data = load_breast_cancer()
X = data.data
y = data.target
```


```python
logistic_regression(X, y)
```

### Exercise

Implement the regularized versions of logistic regression that scikit-learn provides:






```python
def logistic_regression_l1(X, y, C):
    n, p = X.shape
    
    # Convert y to (-1, +1)
    assert np.min(y) == 0
    assert np.max(y) == 1
    Y = y * 2 - 1
    assert np.min(Y) == -1
    assert np.max(Y) == 1

    # Create the model
    m = ConcreteModel()

    # Add variables
    m.b = Var(range(p))
    m.b0 = Var()
    
    # Lasso constraints
    m.absb = Var(range(p))
    def absbeta1(m, j):
        return m.b[j] <= m.absb[j]
    m.absb1 = Constraint(range(p), rule=absbeta1)
    def absbeta2(m, j):
        return -m.b[j] <= m.absb[j]
    m.absb2 = Constraint(range(p), rule=absbeta2)
    
    # Set nonlinear objective function
    m.obj = Objective(sense=minimize, expr=sum(m.absb[j] for j in range(p)) + C * sum(
        log(1 + exp(-Y[i] * (sum(X[i, j] * m.b[j] for j in range(p)) + m.b0)))
        for i in range(n)))

    # Solve the model and get the optimal solutions
    solver = SolverFactory('ipopt', executable=executable)
    
    solver.solve(m)
    return [m.b[j].value for j in range(p)], m.b0.value
```


```python
logistic_regression_l1(X, y, 100)
```


```python
def logistic_regression_l2(X, y, C):
    n, p = X.shape
    
    # Convert y to (-1, +1)
    assert np.min(y) == 0
    assert np.max(y) == 1
    Y = y * 2 - 1
    assert np.min(Y) == -1
    assert np.max(Y) == 1

    # Create the model
    m = ConcreteModel()

    # Add variables
    m.b = Var(range(p))
    m.b0 = Var()
    
    # Set nonlinear objective function
    m.obj = Objective(sense=minimize, expr=0.5 * sum(pow(m.b[j], 2) for j in range(p)) + C * sum(
        log(1 + exp(-Y[i] * (sum(X[i, j] * m.b[j] for j in range(p)) + m.b0)))
        for i in range(n)))

    # Solve the model and get the optimal solutions
    solver = SolverFactory('ipopt', executable=executable)
    
    solver.solve(m)
    return [m.b[j].value for j in range(p)], m.b0.value
```


```python
logistic_regression_l2(X, y, 1000)
```


```python

```

53c0c5aa6501dc8ad0b566476cc1816921d3a560
ipynb
Jupyter Notebook
ML3 - Optimization Modeling (Complete).ipynb
oskali/mban_softwareTools
60b73c798a1f8447de22c46070d023de41d33a30
MIT
2021-03-06T21:16:13.000Z
2021-03-06T21:16:13.000Z
ML3 - Optimization Modeling (Complete).ipynb
oskali/mban_softwareTools
60b73c798a1f8447de22c46070d023de41d33a30
MIT
ML3 - Optimization Modeling (Complete).ipynb
oskali/mban_softwareTools
60b73c798a1f8447de22c46070d023de41d33a30
MIT
2019-12-03T22:35:28.000Z
2021-03-04T00:28:02.000Z
Qwen/Qwen-72B
1. YES
2. YES