3.1 Numerical Linear Algebra

: 20 minutes

NumPy offers functions for matrix operations such as addition, multiplication, dot product, decompositions, etc. These operations also apply to higher-dimensional arrays.

Basic Operations

The dot product of two vectors can be computed using numpy.dot().

Recall that * computes the element-wise product of two arrays. To perform matrix multiplication, one can either use @ or numpy.dot.

numpy.linalg

Matrix decompositions such as LU, SVD, and operations like inverse and determinant are offered through the numpy.linalg module.

Commonly used numpy.linalg functions (McKinney 2017)

Method	Description
diag	Return the diagonal (or off-diagonal) elements of a square matrix as a 1D array, or convert a 1D array into a square matrix with zeros on the off-diagonal
dot	Matrix multiplication
trace	Compute the sum of the diagonal elements
det	Compute the matrix determinant
eig	Compute the eigenvalues and eigenvectors of a square matrix
inv	Compute the inverse of a square matrix
pinv	Compute the Moore-Penrose pseudoinverse of a matrix
qr	Compute the QR decomposition
svd	Compute the singular value decomposition (SVD)
solve	Solve the linear system Ax = b for x, where A is a square matrix

McKinney, Wes. 2017. Python for Data Analysis. 2nd ed. O’Reilly Media. https://www.oreilly.com/library/view/python-for-data/9781491957653/.

Solving A\mathbf{x} = \mathbf{b}

It is tempting to solve a linear system by computing an inverse and multiplying. Do not.

The two agree here, but solve is faster and numerically better behaved. It works by factoring A into triangular pieces—an LU decomposition—and then solving two easy triangular systems. NumPy does this internally; if you ever need the factors themselves they live in scipy.linalg.lu.

Forming inv(A) does strictly more work than the problem requires: it solves n systems (one per column of the identity) in order to answer a question about one right-hand side, and every one of those solves contributes its own rounding error. The rule of thumb worth carrying: if you find yourself writing inv(A) @ b, you wanted solve(A, b).

When A is singular, solve refuses rather than guessing.

That refusal is a feature. A singular A means the columns do not span enough of \mathbb{R}^m to reach every \mathbf{b}, so there is either no solution or infinitely many. Section 3.2 shows what to ask for instead.

The QR decomposition

Any matrix A can be written A = QR, where Q has orthonormal columns (Q^\top Q = I) and R is upper triangular.

Orthonormal columns are the point. As we will see in 3.2, they turn projection—and therefore least squares—into a matter of multiplication rather than inversion.

Note solve vs inv

You need to solve A\mathbf{x} = \mathbf{b} once, for a single dense 2000\times2000 matrix A and a single right-hand side \mathbf{b}. What is the best reason to prefer np.linalg.solve(A, b) over np.linalg.inv(A) @ b?

inv returns a wrong answer whenever A is not symmetric.solve factors A once and then runs two cheap triangular solves, while inv effectively solves 2000 systems and then does one more multiplication—more arithmetic, and more accumulated rounding error, for the same answer.solve succeeds when A is singular, whereas inv does not.The two cost the same; solve is just shorter to type.

solve performs an LU factorization of A and then two triangular solves, at a cost of roughly \tfrac{2}{3}n^3 operations. Computing inv(A) amounts to solving A X = I, i.e. n right-hand sides, and then you still pay an n^2 matrix-vector product—so you do several times the work in order to answer a question about one right-hand side, and every extra operation is another opportunity to lose precision.

The other options are wrong for specific reasons: inv is perfectly correct for non-symmetric invertible matrices; neither function tolerates a singular A (solve raises LinAlgError, and inv does too); and the two are certainly not equal in cost.

Eigenvalues and eigenvectors

A vector \mathbf{v} \neq \mathbf{0} is an eigenvector of a square matrix A when A merely stretches it:

A\mathbf{v} = \lambda\mathbf{v}

The scalar \lambda is the corresponding eigenvalue. Eigenvectors are the directions the matrix does not rotate.

For symmetric matrices—and covariance matrices are symmetric—the eigenvalues are real and the eigenvectors can be chosen orthogonal. Use eigh rather than eig in that case; it is faster and it will not hand you complex numbers from rounding error.

Those eigenvectors are the principal components, which is where 3.3 picks the story up.

For a symmetric A, eigh returns eigenvectors that are not merely orthogonal but orthonormal, collected as the columns of a matrix V with V^\top V = I. That single fact turns the eigendecomposition into something you can write down without an inverse:

A = V\,\Lambda\,V^{-1} = V\,\Lambda\,V^{\top}, \qquad \Lambda = \operatorname{diag}(\lambda_1,\dots,\lambda_n).

This is the spectral decomposition, and it is the symmetric special case of the SVD you will meet in 3.3.

Note Exercise: rebuild a matrix from its spectrum

A is symmetric and positive definite; vals and vecs come from eigh, with the eigenvectors stored as the columns of vecs. Reassemble A from its spectrum. Do not use inv.

answer = vecs @ np.diag(vals) @ vecs.T

np.diag(vals) turns the 1D array of eigenvalues into \Lambda. Scaling the columns directly, (vecs * vals) @ vecs.T, gives the same matrix without building the diagonal.

Condition number

np.linalg.cond measures how much a matrix amplifies error. Formally, if A\mathbf{x} = \mathbf{b} and we perturb the right-hand side to \mathbf{b} + \delta\mathbf{b}, the solution moves by \delta\mathbf{x} with

\frac{\|\delta\mathbf{x}\|}{\|\mathbf{x}\|} \;\le\; \kappa(A)\,\frac{\|\delta\mathbf{b}\|}{\|\mathbf{b}\|}.

The condition number \kappa(A) is the worst-case amplification factor for relative error. A large \kappa means the columns are nearly linearly dependent, and that solving with the matrix will be unreliable no matter how good the algorithm is—it is a property of the problem, not of NumPy.

For a symmetric positive definite matrix there is a clean formula: \kappa(A) = \lambda_{\max}/\lambda_{\min}. In general \kappa(A) = \sigma_{\max}/\sigma_{\min}, a ratio of singular values, which is one more reason 3.3 matters.

Note Exercise: the condition number is a ratio

vals holds the eigenvalues of the symmetric positive definite matrix A, in ascending order. Write the condition number of A using vals alone—no call to cond.

answer = vals.max() / vals.min()

Since eigvalsh returns the eigenvalues in ascending order, vals[-1] / vals[0] works too.

Note Exercise: watch the error get amplified

A is badly conditioned. b and b2 differ in the last entry by 10^{-4}. Compute the amplification factor: the relative change in the solution divided by the relative change in the right-hand side. Predict its size before you run it, then compare it to \kappa(A).

answer = (np.linalg.norm(x2 - x1) / np.linalg.norm(x1)) / \
         (np.linalg.norm(b2 - b) / np.linalg.norm(b))

The solution jumps from (2, 0) to (1, 1). Nothing went wrong in the arithmetic: the two columns of A are nearly identical, so many different \mathbf{x} produce nearly the same A\mathbf{x}, and the data cannot tell them apart.

Note Conditioning

You compute np.linalg.cond(A) for a design matrix A and get roughly 10^{8}. Which statement is true?

A is singular, so no least-squares solution exists.The condition number depends on which NumPy routine you call, so switching from inv to solve will bring it down.A small change in \mathbf{b} can produce a large change in the coefficients \hat{\boldsymbol\beta}, even though the fitted values A\hat{\boldsymbol\beta} barely move.A condition number that large just means the entries of A are large.

The condition number \kappa(A) = \sigma_{\max}/\sigma_{\min} is a property of the matrix, not of the algorithm, and it bounds how much a relative perturbation of the data can be amplified in the solution:

\frac{\|\delta\mathbf{x}\|}{\|\mathbf{x}\|} \le \kappa(A)\,\frac{\|\delta\mathbf{b}\|}{\|\mathbf{b}\|}.

A large \kappa says the columns of A are nearly linearly dependent, so many different coefficient vectors give nearly the same predictions. The coefficients are therefore poorly determined while the fit itself is fine—which is precisely why multicollinearity is so easy to miss by looking at R^2 alone.

\kappa is also scale-free in the sense that matters: multiplying A by 1000 multiplies every singular value by 1000 and leaves \kappa unchanged, so large entries are not the issue. And \kappa = 10^8 is large but finite, so A is not singular; a least-squares solution exists, it is simply untrustworthy.

Keep that amplification in mind. When we reach regression on Nov 3, a design matrix with near-duplicate predictors produces exactly this: the coefficients become enormous and unstable while the fitted values barely move. That is multicollinearity, and ridge regression is the standard response.