Scipy is a python library for scientific computing. It contains modules for linear algebra, root finding, optimization, etc.This tutorial will cover basic uses of scipy for:
The scipy.linalg module contains several useful linear algebra operation. The example below shows how you can use scipy to invert a matrix.
import numpy as np
import scipy.linalg as sp
a = np.array([[1,2],[3,4]])
ainv = sp.inv(a)
print("ainv = ", ainv)
print(a.dot(ainv))
If you want to find the determinant, then use the det() function.
print(sp.det(a))
If you want to solve, a systems of linear equations, use the solve function.
$$ x + 2y = 10 \\ 3x + 4y = 20 $$
b = np.array([[10,20]])
sp.solve(a, b.T)
Scipy has implemented several root finding and optimization algorithms. This website contains all the optimization and root finding methods implemented in scipy and their syntaxes.
For example if want to find the root of $x^3 - 10$ using bisection.
import numpy as np
import scipy.optimize as sp
g = lambda x: (x**3) - 10
x = sp.bisect(g, 1.5, 3.0)
print(x)
If you change the convergence criterion.
# Changing convergence criterion
x = sp.bisect(g, 1.5, 3.0, xtol = 0.1)
print(x)
Scipy also contains other root finding techniques such as false position, Brent, and Newton's method.
Scipy has several optimization routines also. Let us say you want to minimize $f(x) = (x-1)$^2 + e^x$ in the interval [0,2] using Golden Section Method. This can be done using the minimize_scalar() function.
import scipy.optimize as sp
import numpy as np
g = lambda x: (x-1)**2 + np.exp(x)
res = sp.minimize_scalar(g, bracket=(-1,1,2), method='Golden')
print(res.x)
According to the help pages, "bracket defines the bracketing interval and can either have three items (a, b, c) so that a < b < c and fun(b) < fun(a), fun(c)". So pick the bracket interval carefully. It is also relatively easy to code some of the popular optimization algorithms in Python. For example, the golden section method can be coded as shown below.
import numpy as np
g = lambda x: (x-1)**2 + np.exp(x)
def goldenopt(f,x_l, x_u, tol, max_iter):
iter_ctr = 0
while (x_u-x_l)>tol and iter_ctr<max_iter:
x_r1 = x_u - 0.618*(x_u - x_l)
x_r2 = x_l + 0.618*(x_u - x_l)
print(iter_ctr, x_l, x_r1, x_r2, x_u)
if f(x_r1) <= f(x_r2):
x_u = x_r2
x_r2 = x_r1
x_r1 = x_u - 0.618*(x_u - x_l)
else:
x_l = x_r1
x_r1 = x_r2
x_r2 = x_l + 0.618*(x_u - x_l)
iter_ctr = iter_ctr + 1
return (x_r1 + x_r2)/2.0, iter_ctr
x_l_init = 0
x_u_init = 2.0
tolerance = 0.0001
max_iter = 100
x_f, iter_ctr = goldenopt(g, x_l_init, x_u_init, tolerance, max_iter)
print("Final Solution is: ", x_f)
print("Iteration counter: ", iter_ctr)
Similarly, the bisection can be coded as shown below. There are plenty of examples on the web of other good optimization codes.
import numpy as np
g = lambda x: (x-1)**2 + np.exp(x)
g1 = lambda x: 2*(x-1) + np.exp(x)
def bisectopt(f,x_l, x_u, tol, max_iter):
iter_ctr = 0
while (x_u-x_l)>tol and iter_ctr<max_iter:
x_r = (x_l + x_u)/2.0
print(iter_ctr, x_l, x_u, x_r)
if f(x_r) == 0:
return x_r, iter_ctr
elif f(x_r) > 0:
x_u = x_r
else:
x_l = x_r
iter_ctr = iter_ctr + 1
return x_r,iter_ctr
x_l_init = 0
x_u_init = 2.0
tolerance = 0.0001
max_iter = 100
x_f, iter_ctr = bisectopt(g1, x_l_init, x_u_init, tolerance, max_iter)
print("Final Solution is: ", x_f)
print("Iteration counter: ", iter_ctr)
Some of the gradient based methods such as quasi newton methods or newton conjugate are also available in python through the minimize() function. For example, the quasi newton method for minimizing $f(x,y) = (x-1)^4 + 5(y-1)^4 + xy$ using quasi newton is shown below.
import scipy.optimize as sp
import numpy as np
def g(x):
return (x[0]-1)**4 + 5*(x[1]-2)**4 + x[0]*x[1]
x_init = np.array([4.0,4.0])
res = sp.minimize(g, x_init, method='BFGS')
print(res.x)
If you want to use Conjugate Gradient, then you can just change the method parameter.
res = sp.minimize(g, x_init, method='CG')
print(res.x)
Writing your own code for some of these optimization algorithms is not that difficult and should be possible with some practise. For example, the following code implements Newton method for minimizing $f(x)=x^2 - 10x + 25$.
import numpy as np
g = lambda x: x**2 - 10*x + 25
g1 = lambda x: 2*x - 10
def grad_desc(f, f1, x_init, step_size, tol, max_iter):
x_new = x_init
x_old = 0
iter_ctr = 0
while abs(x_new - x_old) > tol and iter_ctr < max_iter:
x_old = x_new
gradient = f1(x_old)
move = gradient * step_size
x_new = x_old - move
iter_ctr = iter_ctr + 1
print(iter_ctr, x_new, f(x_new))
return x_new, iter_ctr
x_init = 15.0
step_size = 0.1
tolerance = 0.01
max_iter = 100
x_f, iter_ctr = grad_desc(g, g1, x_init, step_size, tolerance, max_iter)
print("Final Solution is: ", x_f)
print("Iteration counter: ", iter_ctr)
The following implements multi-dimensional gradient descent for minimizing $f(x,y) = (x-1)^4 + 5(y-1)^4 + xy$.
import numpy as np
def g(x):
return (x[0]-1)**4 + 5*(x[1]-2)**4 + x[0]*x[1]
def g1(x):
grad = np.array([0.0,0.0])
grad[0] = 4*(x[0]-1)**3 + x[1]
grad[1] = 20*(x[1]-2)**3 + x[0]
return grad
def back_track_ls(x):
eta = 1.0
beta = 0.1
gamma = 0.1
while g(x - eta*g1(x)) > g(x) - gamma*eta*np.dot(g1(x),g1(x)):
eta = beta*eta
return eta
def grad_desc_2D(f, f1, x_init, step_size, tol, max_iter):
x_new = np.array([0.0,0.0])
x_old = np.array([0.0,0.0])
gradient = np.array([10.0,10.0])
move = np.array([0.0,0.0])
x_new = x_init
iter_ctr = 0
while np.linalg.norm(x_new - x_old) > tol and iter_ctr < max_iter:
x_old = x_new
gradient = f1(x_old)
step_size = back_track_ls(x_old)
move = gradient * step_size
x_new = x_old - move
iter_ctr = iter_ctr + 1
return x_new, iter_ctr
x_init = np.array([4.0,4.0])
x_f = np.array([0.0,0.0])
step_size = 0.0025
tolerance = 0.01
max_iter = 100
x_f, iter_ctr = grad_desc_2D(g, g1, x_init, step_size, tolerance, max_iter)
print("Final Solution is: ", x_f)
print("Functional Evaluation is: ", g(x_f))
print("Iteration counter: ", iter_ctr)
The following implements multi-dimensional Newton for the same minimization problem.
import numpy as np
def g(x):
return (x[0]-1)**4 + 5*(x[1]-2)**4 + x[0]*x[1]
def g1(x):
grad = np.array([0.0,0.0])
grad[0] = 4*(x[0]-1)**3 + x[1]
grad[1] = 20*(x[1]-2)**3 + x[0]
return grad
def hessian(x):
h = np.array([[0.0,0.0],[0.0,0.0]])
h[0,0] = 12*(x[0] - 1)**2
h[0,1] = 1
h[1,0] = 1
h[1,1] = 60*(x[1] - 2)**2
return h
def newt_2D(f, f1, x_init, step_size, tol, max_iter):
x_new = np.array([0.0,0.0])
x_old = np.array([0.0,0.0])
gradient = np.array([10.0,10.0])
move = np.array([0.0,0.0])
x_new = x_init
iter_ctr = 0
while np.linalg.norm(x_new - x_old) > tol and iter_ctr < max_iter:
x_old = x_new
gradient = f1(x_old)
move = step_size*np.dot(np.linalg.inv(hessian(x_old)),gradient)
x_new = x_old - move
iter_ctr = iter_ctr + 1
return x_new, iter_ctr
x_init = np.array([4.0,4.0])
x_f = np.array([0.0,0.0])
step_size = 1.00
tolerance = 0.01
max_iter = 100
x_f, iter_ctr = newt_2D(g, g1, x_init, step_size, tolerance, max_iter)
print("Final Solution is: ", x_f)
print("Functional Evaluation is: ", g(x_f))
print("Iteration counter: ", iter_ctr)
Linear Programs can be solved using Linprog. The sample code for solving the following LP is shown below:
$$ \text{max}\,\, 90x + 120y \\ 30x + 40y \leq 2000\\ 50x + 20y \leq 2500 \\ 60x + 90y \leq 4000 \\ x,y \geq 0 $$
import numpy as np
import scipy.optimize as sp
Aub = np.array([[30 , 40], [50 , 20] , [60 , 90]])
bub = np.array([2000 , 2500 , 4000])
c = np.array([-90 , -120])
res = sp.linprog(c, A_ub =Aub , b_ub = bub)
if(res. status == 0):
print("Status : Optimal" )
print("Optimal value :", res.fun)
print("X:", res.x)
else :
print("Optimal solution not found")
The objective function should be in minimization form. Equality constraints can be modeled by providing another parameter $A_eq$ in the linprog function. The default assumption is non-negativity of all variables. If you want to model variables with specific bounds as shown below.
$$ \text{max}\,\, 10x + 26y \\ x + y \leq 9\\ x \geq 9.5 \\ 0 \leq y \leq 5 $$
Aub = np.array([[1, 1]])
bub = np.array([9])
c = np.array([10, 26])
b = [(9.5, None), (0, 5)]
res = sp.linprog(c, A_ub=Aub, b_ub=bub, bounds = b)
print(res.status)
The status of 2 indicates the problem is infeasible.
The minimize function can also be used to solves constrained nonlinear optimization problems. Let us look at solving the following constrained optimization problem.
$$ \text{min} \,\, (x_1 - 1)^2 + (x_2 -2.5)^2 \\ x_1 - 2x_2 + 2 \geq 0 \\ -x_1 - 2x_2 + 6 \geq 0 \\ -x_1 + 2x_2 + 2 \geq 0 \\ x_1, x_2 \geq 0 $$
This is example 16.4 from Jorge Nocedal, Stephen Wright. Numerical Optimization, 1999. This is also the same example solved in the scipy docs available here.
Note that for constrained nonlinear optimization problems:
import numpy as np
import scipy.optimize as sp
objfun = lambda x: (x[0] - 1)**2 + (x[1] - 2.5)**2
con1 = lambda x: x[0] - 2 * x[1] + 2
con2 = lambda x: -x[0] - 2 * x[1] + 6
con3 = lambda x: -x[0] + 2 * x[1] + 2
cons = ({'type': 'ineq', 'fun': con1}, {'type': 'ineq', 'fun': con2}, {'type': 'ineq', 'fun': con3})
x_init = (2,0)
res = sp.minimize(objfun, x_init, method='SLSQP', constraints=cons)
print(res.x)
In general for linear and integer programs, I recommend other specialized solvers like CPLEX and Gurobi which all have python interfaces. For other open source solvers which interface with python see here.