Linear Algebra

Question

Is NumPy a good choice for linear algebra?

Answer

Yes, it is. NumPy offers a rich set of routines comparable to Matlab’s one. Matlab is a numeric computing environment that is particularly popular in academia and industry. Currently, there are five families of routines concerning the field of linear algebra at large:

products
decomposition
eigenvalues
norms
equations and inversions

For an overview of NumPy’s linear algebra routines, see the tables below.

Matrix and Vector Products

Routine	Synopsis
`np.dot(a, b[, out])`	Dot product of two arrays
`np.linalg.multi_dot(arrays, *[, out])`	Compute the dot product of two or more arrays in a single function call, while automatically selecting the fastest evaluation order
`np.vdot(a, b, /)`	Return the dot product of two vectors
`np.inner(a, b, /)`	Inner product of two arrays
`np.outer(a, b[, out])`	Compute the outer product of two vectors
`np.matmul(x1, x2, /[, out, casting, order, ...])`	Matrix product of two arrays
`np.tensordot(a, b[, axes])`	Compute tensor dot product along specified axes
`np.einsum(subscripts, *operands[, out, dtype, ...])`	Evaluates the Einstein summation convention on the operands
`np.einsum_path(subscripts, *operands[, optimize])`	Evaluates the lowest cost contraction order for an einsum expression by considering the creation of intermediate arrays
`np.linalg.matrix_power(a, n)`	Raise a square matrix to the (integer) to the power n
`np.kron(a, b)`	Kronecker product of two arrays

Note: the statements included in the ‘Routine’ column assume NumPy is loaded with the np alias.

Decompositions

Routine	Synopsis
`np.linalg.cholesky(a)`	Cholesky decomposition
`np.linalg.qr(a[, mode])`	Compute the qr factorization of a matrix
`np.linalg.svd(a[, full_matrices, compute_uv, ...])`	Singular Value Decomposition

Note: the statements included in the ‘Routine’ column assume NumPy is loaded with the np alias.

Matrix Eigenvalues

Routine	Synopsis
`np.linalg.eig(a)`	Compute the eigenvalues and right eigenvectors of a square array
`np.linalg.eigh(a[, UPLO])`	Return the eigenvalues and eigenvectors of a complex Hermitian (conjugate symmetric) or a real symmetric matrix
`np.linalg.eigvals(a)`	Compute the eigenvalues of a general matrix
`np.linalg.eigvalsh(a[, UPLO])`	Compute the eigenvalues of a complex Hermitian or real symmetric matrix

Note: the statements included in the ‘Routine’ column assume NumPy is loaded with the np alias.

Norms and Other Numbers

Routine	Synopsis
`np.linalg.norm(x[, ord, axis, keepdims])`	Matrix or vector norm
`np.linalg.cond(x[, p])`	Compute the condition number of a matrix
`np.linalg.det(a)`	Compute the determinant of an array
`np.linalg.matrix_rank(A[, tol, hermitian])`	Return matrix rank of array using SVD method
`np.linalg.slogdet(a)`	Compute the sign and (natural) logarithm of the determinant of an array
`np.trace(a[, offset, axis1, axis2, dtype, out])`	Return the sum along diagonals of the array

Note: the statements included in the ‘Routine’ column assume NumPy is loaded with the np alias.

Solving Equations and Inverting Matrices

Routine	Synopsis
`np.linalg.solve(a, b)`	Solve a linear matrix equation, or system of linear scalar equations
`np.linalg.tensorsolve(a, b[, axes])`	Solve the tensor equation a x = b for x
`np.linalg.lstsq(a, b[, rcond])`	Return the least-squares solution to a linear matrix equation
`np.linalg.inv(a)`	Compute the (multiplicative) inverse of a matrix
`np.linalg.pinv(a[, rcond, hermitian])`	Compute the (Moore-Penrose) pseudo-inverse of a matrix
`np.linalg.tensorinv(a[, ind])`	Compute the ‘inverse’ of an N-dimensional array

Note: the statements included in the ‘Routine’ column assume NumPy is loaded with the np alias.

Matrix and Vector Products

Question

Matrix and vector products

Answer

The code snippet below shows how to use NumPy to get dot, inner, and outer products — whose definitions are:

Dot product: \(x \cdot y = \sum_{i=1}^n x_{i} y_{i}\) where \(x\) and \(y\) are vectors of length \(n\).

Inner product: \(\langle x, y\rangle = x^{T} y\) where \(x^{T}\) is the transpose of \(x\).

Outer product: \(x \otimes y = x y^{T}\) where \(y^{T}\) is the transpose of \(y\).

# import numpy with the socially accepted alias 'np'
>>> import numpy as np

# the arrays
>>> X = np.arange(1, 11, 1)
>>> X
array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10])

>>> Y = np.arange(10, 0, -1)
>>> Y
array([10,  9,  8,  7,  6,  5,  4,  3,  2,  1])

# dot product
>>> np.dot(X, Y)
220

# inner product
# --+ let's reshape the arrays
>>> X = X.reshape(5, 2)
>>> X
array([[ 1,  2],
       [ 3,  4],
       [ 5,  6],
       [ 7,  8],
       [ 9, 10]])
>>> Y = Y.reshape(5, 2)
>>> Y
array([[10,  9],
       [ 8,  7],
       [ 6,  5],
       [ 4,  3],
       [ 2,  1]])
>>> np.inner(X, Y)
array([[ 28,  22,  16,  10,   4],
       [ 66,  52,  38,  24,  10],
       [104,  82,  60,  38,  16],
       [142, 112,  82,  52,  22],
       [180, 142, 104,  66,  28]])

# outer product
>>> np.outer([0, 1, 2], [4, 5, 6, 7])
array([[ 0,  0,  0,  0],
       [ 4,  5,  6,  7],
       [ 8, 10, 12, 14]])

Solving a System of Equations

Question

Solving a system of equations

Answer

The code snippet below shows how to solve the following system of equations:

\[\left\{ \begin{array}{lr} 4x + 2y = 2000\\ 7x + 13y = 3000 \end{array} \right.\]

# import numpy with the socially accepted alias 'np'
>>> import numpy as np

# the system of equations
# --+ left-hand side
>>> a = np.array([[4, 2], [7, 13]])
# --+ right-hand side
>>> b = np.array([2000, 3000])

# solve the system of equations
>>> np.linalg.solve(a, b)
array([526.31578947, -52.63157895])

Least-Square Estimation

Question

How do I carry out a least-square estimation in NumPy?

Answer

The least-square estimator is a popular choice for estimating the relationship between two variables, one of the most common tasks in statistics and econometrics. Students who want to familiarize with the field of econometrics are warmly encouraged to read Angrist, Joshua D., and Jörn-Steffen Pischke. Mostly harmless econometrics: An empiricist’s companion. Princeton university press, 2009.

Specifically, the least-square estimator is defined as follows:

\[\hat{\beta} = (X^T X)^{-1} X^T y\]

Where \(\hat{\beta}\) is the estimate of the regression coefficients, \(X\) is the matrix of regressors, \(y\) is the outcome vector, and \(X^T X\).

In the code snippet below, we create two arrays, x — playing the role of the regressor — and y — which we pretend to be the outcome variable. Then, we arrange the data in an array containing x and a vector of ones. In this way, we add an intercept to the model, namely a scalar that does not change across the observations in the data (see line 11). Finally, we estimate the two regression coefficients of interest:

The regression coefficient for x, capturing the expected change in y for a unitary increase in x
The regression coefficient for the intercept, namely, the expected value of y when the association between x and y is partialled-out

Let me stress that .linalg.lstsq returns four pieces of information:

An array with the estimated regression coefficients
The sums of squared residuals (i.e., the differences between the observed y and the values predicted by the model)
The rank of the matrix containing the regressors
The singular values of the matrix containing the regressors

Thus, in line 28, we fetch the first element of the returned array, i.e., the regression coefficients.

# import numpy with the socially accepted alias 'np'
>>> import numpy as np

# the arrays
>>> x = np.arange(0, 15, 1)
>>> y = [3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31]

# data preparation / adding a constant to the model
>>> A = np.vstack([x, np.ones(len(x))]).T
>>> A
array([[ 0.,  1.],
       [ 1.,  1.],
       [ 2.,  1.],
       [ 3.,  1.],
       [ 4.,  1.],
       [ 5.,  1.],
       [ 6.,  1.],
       [ 7.,  1.],
       [ 8.,  1.],
       [ 9.,  1.],
       [10.,  1.],
       [11.,  1.],
       [12.,  1.],
       [13.,  1.],
       [14.,  1.]])

# estimate the regression coefficients (a.k.a., the regression slopes) of
# the linear model
>>> b = np.linalg.lstsq(A, y)[0]
>>> b
array([2., 3.])