import numpy as np5 Linear algebra and optimization
In this chapter we will see how Python can be used for performing linear algebra tasks, like solving systems of linear equations, and (integer) linear optimization.
Linear algebra forms the basis of many statistical tasks like principle component analysis (PCA), which can be used to identify important properties/dimensions of large data sets. This typically makes, e.g., machine learning algorithms more efficient as this avoid the so-called curse of dimensionality.
Applications of (integer) linear optimization range from scheduling, routing, transportation and food aid problems to solving Sudoku puzzles.
Let us first again import NumPy. Later on we will also use its linear algebra package (or module) linalg.
5.1 Linear algebra
In this section we will often refer to one-dimensional arrays as row or column vectors, and to two-dimensional arrays as matrices.
5.1.1 Matrix multiplications
Recall from Section 4.1.1 that for an m \times n matrix A and m \times 1 column vector x, the command A*x gives a matrix in which every column of A gets pointwise multiplied by the column vector x, that is,
\begin{align*} A*x = & \left[ \begin{matrix} a_{00} & a_{01} & \dots & a_{0(n-1)} \\ a_{10} & a_{11} & \dots & a_{1(n-1)} \\ \vdots & \vdots & \ddots & \vdots \\ a_{(m-1)0} & a_{(m-1)1} & \dots & a_{(m-1)(n-1)} \end{matrix} \right]*\left[ \begin{matrix} x_{0} \\ x_{1} \\ \vdots \\ x_{(m-1)} \\ \end{matrix} \right] \\ = & \left[ \begin{matrix} a_{00} \cdot x_0 & a_{01} \cdot x_0& \dots & a_{0(n-1)} \cdot x_0\\ a_{10} \cdot x_1 & a_{11} \cdot x_1& \dots & a_{1(n-1)} \cdot x_1\\ \vdots & \vdots& \ddots & \vdots \\ a_{(m-1)0} \cdot x_{m-1} & a_{(m-1)1} \cdot x_{m-1} & \dots & a_{(m-1)(n-1)} \cdot x_{m-1} \end{matrix} \right] \end{align*}
A = np.array([
[1,2],
[3,4]])
x = np.array([1,5])
x = x[:,None] # Turn x into (2,1) shaped column vector
print('A = \n', A)
print('x = \n', x)
print('A*x = \n', A*x)A =
[[1 2]
[3 4]]
x =
[[1]
[5]]
A*x =
[[ 1 2]
[15 20]]In the linear algebra sense, the multiplication Ax gives an m \times 1 column vector defined by
\begin{align*} Ax = & \left[ \begin{matrix} a_{00} & a_{01} & \dots & a_{0(n-1)} \\ a_{10} & a_{11} & \dots & a_{1(n-1)} \\ \vdots & \vdots & \ddots & \vdots\\ a_{(m-1)0} & a_{(m-1)1} & \dots & a_{(m-1)(n-1)} \end{matrix} \right]\left[ \begin{matrix} x_{0} \\ x_{1} \\ \vdots \\ x_{(m-1)} \\ \end{matrix} \right] \\ = & \left[ \begin{matrix} a_{00}x_0 + a_{01}x_1 + \dots + a_{0(n-1)}x_{n-1}\\ a_{10}x_0 + a_{11}x_1 + \dots + a_{1(n-1)}x_{n-1}\\ \vdots \\ a_{(m-1)0}x_0 + a_{(m-1)1}x_1 + \dots + a_{(m-1)(n-1)}x_{n-1}\\ \end{matrix} \right] \end{align*}
This can be achieved in Python with the @ (‘at’) operator, that is, by writing A @ x (the white spaces are not needed; they are included here for readability).
print('A = \n', A)
print('x = \n', x)
print('Ax = \n', A @ x)A =
[[1 2]
[3 4]]
x =
[[1]
[5]]
Ax =
[[11]
[23]]If we talk about matrix multiplication, it will always be clear from context whether we mean A*x or A @ x.
We can also use @ to multiply a k \times m and m \times n matrix with each other in the linear algebra sense.
A = np.array([
[1,2],
[3,4],
[3,7]])
B = np.array([
[2,5,1,7],
[3,4,-1,8]])
# Note that k = 3, m = 2, n = 4
print('AB = \n',A @ B)AB =
[[ 8 13 -1 23]
[18 31 -1 53]
[27 43 -4 77]]A special case of this is if we multiple a column vector
x = \left[\begin{matrix} x_0\\ x_1 \\ \vdots \\x_{n-1} \end{matrix} \right]
with its transpose (being a row vector) resulting in a matrix that contains x_ix_j on position (i,j).
x = np.array([
[1],
[2],
[4]])
print('xx^T = \n', x @ x.T)xx^T =
[[ 1 2 4]
[ 2 4 8]
[ 4 8 16]]5.1.2 Matrix properties
Using the linalg subpackage from NumPy, we can compute various well-known properties of a matrix (that you should recall from your Linear Algebra course).
linalg.matrix_rank(A): Computes rank of matrix A.linalg.det(A): Computes determinant of square A.linalg.eig(A):Computes eigenvalues and (right) eigenvectors of matrix A, i.e., values \lambda that satisfy Av = \lambda v for some (eigen)vector v.linalg.inv(A): Computes the inverse matrix A^{-1} of A, that is, the matrix that satisfies A^{-1}A = AA^{-1} = I, where I is the identity matrix.linalg.norm(A): Computes the norm of matrix (or vector).
We can access these function with the syntac np.linalg.function_name(A).
Finally, there is also the trace function trace(A) in NumPy that computes the trace of A, i.e., the sum of the elements on the diagonal. This function is implemented directly in NumPy, so you don’t have to go to the subpackage linalg first.
Let us look at some examples of these commands.
A = np.array([
[1,2],
[3,4]])# Rank of A
print('Rank of A: ', np.linalg.matrix_rank(A))Rank of A: 2# Determinant of A
print('Determinant of A: ', np.linalg.det(A))Determinant of A: -2.0000000000000004# Trace of A
print('Trace of A: ', np.trace(A))Trace of A: 5# Inverse of A
print('Inverse of A: \n', np.linalg.inv(A))Inverse of A:
[[-2. 1. ]
[ 1.5 -0.5]]The function linalg.eig() is special in that it does not output one, but two arrays.
The first array is one-dimensional and contains the eigenvalues [\lambda_0,\dots,\lambda_{n-1}] of A; the second array is two-dimensional, and contains the corresponding eigenvectors v^0,\dots,v^{n-1} as columns of the matrix. The eigenvalues and eigenvectors are paired in the sense that Av^i = \lambda_iv^i for i = 0,\dots,n-1.
# Eigenvalues and eigenvectors of A
lambdas, V = np.linalg.eig(A)
print('Eigenvalues of A: \n', lambdas)
print('Matrix with eigenvectors of A: \n', V)Eigenvalues of A:
[-0.37228132 5.37228132]
Matrix with eigenvectors of A:
[[-0.82456484 -0.41597356]
[ 0.56576746 -0.90937671]]You can output only the eigenvalues or the eigenvectors by suppressing the other output argment with _.
# Only eigenvalues of A
lambdas, _ = np.linalg.eig(A)
# Only eigenvectors of A
_, V = np.linalg.eig(A)If you want to recall how you can handle multiple outputs, and suppress the ones that your are not interested, have a look at Appendix B.
Let us check that the eigenvalues and eigenvectors indeed satsify Av^i = \lambda_iv^i for i = 0,1.
v0 = V[:,0]
v1 = V[:,1]
print('Verify first eigenvector:', A @ v0 - lambdas[0]*v0)
print('Verify second eigenvector:', A @ v1 - lambdas[1]*v1)Verify first eigenvector: [0.00000000e+00 5.55111512e-17]
Verify second eigenvector: [0. 0.]The function linalg.norm(A) computes the L^2-norm of a matrix A = (a_{ij})_{i = 1,\dots,m, j = 1,\dots,n}, also known as the Frobenius norm. It is defined as
||A||_{\text{F}} = \sqrt{\sum_{i=0}^{m-1} \sum_{j = 0}^{n-1} a_{ij}^2 }.
If you apply the function two a row or column vector x = [x_0,\dots,x_{n-1}], you get the usual L^2-norm defined as
||x||_{2} = \sqrt{\sum_{i=0}^{n-1} x_i^2 }.
You can also use this function in a vectorized manner. For example, if you want to compute the L^2-norm of every row in a matrix, you can use the axis keyword argument.
A = np.array([
[1,2],
[3,4],
[5,6]])
row_norms = np.linalg.norm(A,axis=1)
print(row_norms)[2.23606798 5. 7.81024968]5.1.3 Equation solving
For a given m \times n matrix A = (a_{ij}) and m \times 1 column vector b = [b_0,\dots,b_{m-1}]^T, perhaps the most important question in linear algebra is to compute an x = [x_0,\dots,x_{n-1}]^T that satisfies Ax = b, i.e., the linear systems of equalities \begin{array}{rrrrcrrrl} a_{00}x_0 &+& a_{01}x_1 &+& \dots &+& a_{0(n-1)}x_{n-1} &=& b_0 \\ a_{10}x_0 &+& a_{11}x_1 &+& \dots &+& a_{1(n-1)}x_{n-1} &=& b_1 \\ &&&& \vdots &&&& \\ a_{(m-1)0}x_0 &+& a_{(m-1)1}x_1 &+& \dots &+& a_{(m-1)(n-1)}x_{n-1} &=& b_{m-1} \end{array}.
If A is a square, invertible matrix then this can be done with linalg.solve(A,b). This gives the unique solution x = A^{-1}b.
A = np.array([
[1,2],
[3,4]])
b = np.array([[1],[4]])
print('Solution to Ax = b: \n', np.linalg.solve(A,b))Solution to Ax = b:
[[ 2. ]
[-0.5]]In fact, you can compute x directly as x = A^{-1}b using the inverse of A that we have seen earlier.
print('Solution to Ax = b: \n', np.linalg.inv(A) @ b)Solution to Ax = b:
[[ 2. ]
[-0.5]]If the system is overdetermined, then there is not necessarily a unique solution. This can happen when there are more constraints than variables, i.e., when m > n. In this case, we can find a solution x that is approximately optimal with the least squares method. It finds a solution x that solves the (non-linear) problem
\min_x ||Ax - b||_2.
The least squares method can be executed with linalg.lstsq(A,b); we will see this method in the next chapter, when we consider non-linear optimization problems. This function can also be used to solve underdetermined systems, where n > m.
5.2 Linear optimization
Recall that in a linear optimization problem the goal is to compute a vector x = [x_0,\dots,x_{n-1}] of decision variables that maximizes, or minimizes, a linear objective function subject to a collection of linear constraints, which can either be a \leq, \geq or = constraint.
In this section we will consider a general linear optimization problem of the form \begin{array}{ll} \min & c^Tx \\ {\rm s.t.} & Ax = b \\ & Wx \leq z \\ & \ell \leq x \leq u \end{array} where c, \ell, u \in \mathbb{R}^n, A \in \mathbb{R}^{m \times n}, b \in \mathbb{R}^m, U \in \mathbb{R}^{k \times n}, and z \in \mathbb{R}^k are the input data, or input parameters.
Written out this means we have a system with m equality constraints, k inequality constraints and upper and lower bounds on the decision variables.
\begin{array}{rrrcrrlll} \min & c_1 x_1 &+& \cdots& + & c_nx_n& & & \\ {\rm s.t.} & a_{01}x_0& +& \cdots& +& a_{0(n-1)} x_{n-1} &= &b_0 &\\ & & & \vdots & & & & \\ & a_{(m-1)1}x_0& +& \cdots& +& a_{(m-1)(n-1)} x_{n-1} &= &b_{m-1}& \\ & w_{01}x_0& +& \cdots& +& w_{0(n-1)} x_{n-1} &\leq &z_0 &\\ & & & \vdots & & & & &\\ & w_{(k-1)1}x_0& +& \cdots& +& w_{(k-1)(n-1)} x_{n-1} &\leq &z_{k-1} & \\ \ell_0 \leq & x_0 & && & & \leq & u_0& \\ & & & \ddots & & & & & \\ \ell_{n-1} \leq & & & & & x_{n-1} & \leq & u_{n-1}& \\ \end{array}
You can model a \geq constraint by multiplying it with -1, so that it turns into a \leq constraint.
If, in addition, we require that x_i \in \mathbb{N} for i = 1,\dots,n, i.e., that the variables are integral, then we call the above problem an integer linear optimization problem.
We will look at two packages that can be used to solve (integer) linear optimization problems. The difference is that the first function linprog allows us to explicitly input the data c, A, B, \ell and u, whereas the second package pulp allows us to define constraints one-by-one. The latter is more convenient if the input data his not given explicitly.
5.2.1 Explicit input data
If the input data is given explicitly, we can solve linear optimization problems quickly with the linprog function, which is part of the optimize package of SciPy. We will see more functionality of SciPy in later chapters of this book.
We can import linprog as follows:
from scipy.optimize import linprogThis package is suitable for solving problem of the general form above. We will implement the following example:
\begin{array}{lrrrrrrl} \text{max } & z = & 15x_{1} &+ &20x_{2} & & & \\ \text{s.t.} & & 2x_1 &+ & 2x_2& \leq & 8\frac{1}{2} & \\ & & x_1 &+ & 2x_2& \leq & 6 & \\ && 12x_1 &+ & 17x_2& = & 51 & \\ & &x_1 & & & \geq & 0 & \\ & && & x_2 & \geq & 0 & \end{array}
Because linprog works with explicitly input data, we first define these. We define all the input data as NumPy arrays, but you can actually also use plain lists for this.
Note that we multiply the objective function c with -1, so that our maximization problem turns into a minimization problem.
## Define input data
# Coefficients of the objective function
c = np.array([-15, -20]) # Minimize -15x1 - 20x2
# Coefficients of the equality constraints (Ax = b)
A = np.array([[12, 17]]) # 12x1 + 17x2 = 51
b = np.array([51])
# Coefficients of the inequality constraints (Ux <= z)
U = np.array([
[2, 2], # 2x1 + 2x2 <= 8.5
[1, 2]]) # x1 + 2x2 <= 6
z = np.array([[8.5], [6]])The bounds on x_i have to be inputted as tuple (\ell_i,u_i). If the lower or upper bound is not present, meaning either \ell_i = -\infty or u_i = + \infty, you can replace the entry with None.
# Bounds for the variables x1 and x2 (l = 0, u = infinity)
x1_bounds = (0, None) # x1 >= 0
x2_bounds = (0, None) # x2 >= 0We can solve the problem with linprogr() It takes input arguments
c: The vector with objective function coefficientsA_eq: The constraint matrix for the equality constraintsb_eq: The right hand side vector of the equality constraintsA_ub: The constraint matrix for the inequality constraintsb_ub: The right hand side vector of the inequality constraintsbounds: List of tuples with each tuple having the upper and lower bound for the corresponding variable
Note that A_eq and b_eq can be left out as input arguments if there are no equality constraints; the same holds for A_ub and b_ub.
# Solve the linear programming problem
# (You can use \ to continue command on next line)
result = linprog(c, A_eq=A, b_eq= b, \
A_ub=U, b_ub=z, \
bounds=[x1_bounds, x2_bounds])The result variable contains various aspects of the solution. The two most important for us are
result.x: The optimal solution xresult.fun: The function value attained by the the optimal solution.
# Output the results
print('The problem is solved by', result.x, \
'with objective function value', result.fun)The problem is solved by [4.25 0. ] with objective function value -63.75If we want to solution to be integral, we add the input argument integrality=1. See the documentation to read about this and to find more options.
# Solve the integer linear optimization problem
# (You can use \ to continue command on next line)
result_integral = linprog(c, A_eq=A, b_eq= b, \
A_ub=U, b_ub=z, \
bounds=[x1_bounds, x2_bounds], \
integrality=1)
# Output the results
print('The problem is solved by', result_integral.x, \
'with objective function value', result_integral.fun)The problem is solved by [0. 3.] with objective function value -60.0The complete code of this section can be found below.
Show complete code for linear optimization with linprog
import numpy as np
from scipy.optimize import linprog
## Define input data
# Coefficients of the objective function
c = np.array([-15, -20]) # Minimize -15x1 - 20x2
# Coefficients of the equality constraints (Ax = b)
A = np.array([[12, 17]]) # 12x1 + 17x2 = 51
b = np.array([51])
# Coefficients of the inequality constraints (Ux <= z)
U = np.array([
[2, 2], # 2x1 + 2x2 <= 8.5
[1, 2]]) # x1 + 2x2 <= 6
z = np.array([[8.5], [6]])
## Bounds for the variables x1 and x2 (l = 0, u = infinity)
x1_bounds = (0, None) # x1 >= 0
x2_bounds = (0, None) # x2 >= 0
## Solve the linear optimization problem
## (You can use \ to continue command on next line)
result = linprog(c, A_eq=A, b_eq= b, \
A_ub=U, b_ub=z, \
bounds=[x1_bounds, x2_bounds])
# Output the results
print('The linear optimization problem is solved by', result.x, \
'with objective function value', result.fun)
##################################
## With integrality constraints ##
##################################
# Solve the problem
result_integral = linprog(c, A_eq=A, b_eq= b, \
A_ub=U, b_ub=z, \
bounds=[x1_bounds, x2_bounds], \
integrality=1)
# Output the results
print('The integer linear optimization problem is solved by', \
result_integral.x, 'with objective function value', result_integral.fun)5.2.2 Implicit input data
If the constraints are not given explicitly, creating the input data matrices like A and U can be quite tedious. In this case, it works better to define the constraints directly, based on the given problem description.
In this section we will consider the maximum weight bipartite matching problem. We are given a (complete) bipartite graph G = (V,W,E) with node sets V = \{0,\dots,n-1\} and W = \{0,\dots,n-1\} and edges ij \in E = \{ab : a \in V, b \in W\} connecting the nodes i \in V with j \in W. Every such edge has a weight w_{ij}. We want to match up every node in V with a node in W so that the total summed weight of the selected edges is maximal.
We introduce (binary) decision variables x_{ij} with the interpretation that x_{ij} = 1 if i and j get matched, and x_{ij} = 0 otherwise. The constraints then mean that every node i \in V should get matched up with precisely one j \in W and vice versa.
\begin{array}{lllll} \max & \displaystyle \sum_{i=0}^{n-1} \sum_{j=0}^{n-1} w_{ij}x_{ij} & & \\ \rm{s.t. } & \displaystyle \sum_{i=0}^{n-1} x_{ij} & =& 1 & \text{for } j = 0,\dots,n-1 \\ &\displaystyle \sum_{j=0}^{n-1} x_{ij} & = &1 & \text{for } i = 0,\dots,n-1\\ & x_{ij} &\geq & 0 & \text{for all } i,j = 0,\dots,n-1 \\ & x_{ij} & \in & \{0,1\} & \text{for all } i,j = 0,\dots,n-1 \end{array}
In fact, the last conditions defining the decisions variables to be binary, are redundant: If we leave them out then the optimal solution that is found will satisfy them regardless (you might have seen this in a Combinatorial Optimization course). For now, we will ignore this mathematical fact, and define the decision variables as binary variables.
Defining the constraint matrix explicitly is quite a hassle, and given that many of its entries are zero, it is not very efficient for Python to store it explicitly, so using linprog here is not very convenient.
Instead, in this section we will use the pulp package to solve this problem.
import pulpWe create some input data (weights) for the problem that we are going to build.
# Size of the node sets
n = 4
# Weights of the edges ij
w = np.arange(1,n**2+1).reshape(n,n)
print(w) [[ 1 2 3 4]
[ 5 6 7 8]
[ 9 10 11 12]
[13 14 15 16]]Problem instantiation. We first initialize the problem object that we want to solve with LpProblem() that takes as input:
- Problem name: String variable
- Problem type: Maximization (
pulp.LpMaximize), or minimization (pulp.LpMinimize).
prob = pulp.LpProblem("Weighted_Bipartite_Matching", pulp.LpMaximize)As opposed to linprog we need to explicitly instantiate the decision variables, so that we can use them later to add constraints to our problem.
Decision variables. A decision variable y can be instantiated with LpVariable() that takes as input
- Variable name: String being name of variable
- Variable category: Keyword argument
catbeing Binary, Integer or Continuous (default ifcatis not specified).
pulp.LpVariable('y',cat='Binary')To initiate our decision variables we need to loop over the indices i and j, which can be done as follows. Note that the syntax f"y_{i}_{j}" allows us to incorporate/format the loop indices i and j into the variable name (which is a string).
x = np.zeros((n,n), dtype=object) # Data type of variable is 'object'
for i in range(n):
for j in range(n):
x[i,j] = pulp.LpVariable(f"x_{i}_{j}", cat='Binary') This can be done more compactly using list comprehension, which is a quick alternative for doing the for-loops explicitly. Avoiding for-loops altogether here is not possible, because variables cannot be generated in a vectorized manner if we want to give them all a separate name.
# With list comprehension
x = np.array([[pulp.LpVariable(f"x_{i}_{j}", cat='Binary') \
for j in range(n)] for i in range(n)])In both cases we have created an n \times n NumPy array whose elements are decision variables x_{ij}.
# Print names of decision variables
print(x)[[x_0_0 x_0_1 x_0_2 x_0_3]
[x_1_0 x_1_1 x_1_2 x_1_3]
[x_2_0 x_2_1 x_2_2 x_2_3]
[x_3_0 x_3_1 x_3_2 x_3_3]]We will use these variables to define the objective function and constraints.
Objective function. We next construct the objective function
\begin{array}{lllll}
\max & \displaystyle \sum_{i=0}^{n-1} \sum_{j=0}^{n-1} w_{ij}x_{ij} & &
\end{array}
Because both the weight matrix w and the decision variables x are stored in NumPy arrays, we can do a pointwise multiplication of these arrays, and sum up the elements of the resulting array, to get the objective.
Summing up decision variables in arrays, in our case the array w*x, can be done using lpSum() within PuLP. To add the objective to the problem, and later constraints, we use the syntax prob += [objective] where [objective] contains the objective function expresssion in terms of the decision variables that we instantiated. You should think of the syntax += as adding something to the instantiated problem prob.
prob += pulp.lpSum(w*x) Note that we write lpSum() and not LpSum() with a capital L. This has to do with the fact that lpSum() is a function, whereas LpProblem and LpVariable are classes (that have different naming conventions in Python modules).
# Print objective function
print(prob.objective)x_0_0 + 2*x_0_1 + 3*x_0_2 + 4*x_0_3 + 5*x_1_0 + 6*x_1_1 + 7*x_1_2 + 8*x_1_3 + 9*x_2_0 + 10*x_2_1 + 11*x_2_2 + 12*x_2_3 + 13*x_3_0 + 14*x_3_1 + 15*x_3_2 + 16*x_3_3Constraints. Finally, we add the constraints
\begin{array}{lllll}
& \displaystyle \sum_{i=0}^{n-1} x_{ij} & =& 1 & \text{for } j = 0,\dots,n-1 \\
&\displaystyle \sum_{j=0}^{n-1} x_{ij} & = &1 & \text{for } i = 0,\dots,n-1
\end{array}
to the problem, which can also be done with the prob += [constraint] syntax, where [constraint] is the constraint we want to add. We can add our constraints as follows.
for j in range(n):
prob += pulp.lpSum([x[i,j] for i in range(n)]) == 1
for i in range(n):
prob += pulp.lpSum([x[i,j] for j in range(n)]) == 1To print the constraints, we need to realize that they are stored in an (ordered) dictionary. We can print the keys and values of this dictipnary as follows.
# Print constraints
for key, value in prob.constraints.items():
print(f"{key}: {value}")_C1: x_0_0 + x_1_0 + x_2_0 + x_3_0 = 1
_C2: x_0_1 + x_1_1 + x_2_1 + x_3_1 = 1
_C3: x_0_2 + x_1_2 + x_2_2 + x_3_2 = 1
_C4: x_0_3 + x_1_3 + x_2_3 + x_3_3 = 1
_C5: x_0_0 + x_0_1 + x_0_2 + x_0_3 = 1
_C6: x_1_0 + x_1_1 + x_1_2 + x_1_3 = 1
_C7: x_2_0 + x_2_1 + x_2_2 + x_2_3 = 1
_C8: x_3_0 + x_3_1 + x_3_2 + x_3_3 = 1 Solving the problem. We solve the problem with the solve() function.
We can access the values of the objective function and the variables using pulp.value(). The objective function of the problem is stored in prob.objective, and so its value in the optimized model can be accessed with pulp.value(prob.objective).
We want to represent the optial matching nicely in a binary matrix. Unfortunately, printing the values of the optimized model is difficult to do without for-loops since the function pulp.value(x[i,j]) that returns the value of variable x_{ij} only works for numbers and not lists or arrays.
# Solve the problem
prob.solve()
# Store the results in matrix 'matching'
matching = np.zeros((n,n), dtype=int)
for i in range(n):
for j in range(n):
matching[i,j] = pulp.value(x[i,j])
# Print solution
print("The optimal matching is: \n", matching)
print(f"Optimal value of the objective function: {pulp.value(prob.objective)}")The optimal matching is:
[[1 0 0 0]
[0 1 0 0]
[0 0 0 1]
[0 0 1 0]]
Optimal value of the objective function: 34.0The complete code of this section can be found below.
Show complete code for linear optimization with PuLp
import pulp
# Size of the node sets (input data)
n = 4
# Weights of the edges ij (input data)
w = np.arange(1,n**2+1).reshape(n,n)
# Instantiate problem
prob = pulp.LpProblem("Weighted_Bipartite_Matching", pulp.LpMaximize)
# Instantiate decision variables and store in NumPy array
x = np.array([[pulp.LpVariable(f"x_{i}_{j}", cat='Binary') \
for i in range(n)] for j in range(n)])
# Set objective function
prob += pulp.lpSum(w*x)
# Set constraints
for j in range(n):
prob += pulp.lpSum([x[i,j] for i in range(n)]) == 1
for i in range(n):
prob += pulp.lpSum([x[i,j] for j in range(n)]) == 1
# Store the results in matrix 'matching'
matching = np.zeros((n,n), dtype=int)
for i in range(n):
for j in range(n):
matching[i,j] = pulp.value(x[i,j])
# Print solution and objective value
print("The optimal matching is: \n", matching)
print(f"Optimal value of the objective function: {pulp.value(prob.objective)}")5.2.3 Remarks
The linprog function is typically convenient for small linear optimization problems, especially for problems whose input data (such as the constraint matrices) you can define explicitly.
The pulp package is more useful when the input data is defined implicitly, e.g., using different sets of inequalities defined throug indices. Furthermore, when problems become of larger scale, as you will typically encounter in practice, it is better to use state-of-the-art solvers such as CPLEX, GUROBI, or MOSEK. PuLP is able to be coupled to such solvers, that is, you can define a problem in PuLP and then have it solved with the external solver software.
Many of these solvers can also handle other problems such as mixed integer linear optimization problems, in which there is a mixture of continuous, integer and binary variables, and various non-linear optimization problems. We will see more of those in the next chapter.