Linear Algebra Essentials

Solving systems, computing eigenvalues, and the linear-algebra patterns that hide inside almost every scientific algorithm

If numerical computing has a secret, it is this: almost every scientific problem reduces to a linear algebra problem. A nonlinear least squares fit? A sequence of linear systems. A PDE discretization? A sparse linear system. A neural network training step? A chain of matrix multiplies and a linear solve. Eigenvalue problems? They are everywhere — vibration modes, principal component analysis, page ranking, quantum mechanics.

This page is the working tour of the linear algebra primitives in NumPy and SciPy. It assumes you know what a matrix and a vector are. It does not assume you remember how to find an eigenvalue by hand.

The four problem types

Most of scientific computing falls into one of these four buckets. The next chapter will cover the factorizations behind each. This chapter is about using them.

Solving linear systems

The single most-used routine in scientific Python is numpy.linalg.solve(A, b). It solves $A x = b$ where $A$ is square and (numerically) non-singular. It calls LAPACK's dgesv, which uses LU factorization with partial pivoting — the backward-stable algorithm we mentioned in the stability chapter.

The residual is at machine epsilon — that is the right answer to seven sig figs. Never compute np.linalg.inv(A) @ b as a substitute. It is slower, less stable, and rarely what you want.

Special structure: trust SciPy

A general solve works on any non-singular matrix. But many real problems have structure — symmetric, positive-definite, banded, triangular, sparse — and dedicated solvers can be ten or a hundred times faster. Use them.

Structure	NumPy / SciPy call	Speedup vs general
Triangular	`scipy.linalg.solve_triangular`	~ $3\times$
Symmetric positive definite	`scipy.linalg.cho_solve` (Cholesky)	~ $2\times$
Tridiagonal	`scipy.linalg.solve_banded`	~ $100\times$ on large $n$
Sparse	`scipy.sparse.linalg.spsolve`	$\propto$ sparsity

When you know your matrix is SPD, choosing Cholesky over LU buys you a factor of $\sim 2$ in time and a factor of $\sim 2$ in stability. Cholesky requires no pivoting.

Eigenvalues and eigenvectors

An eigenvector of $A$ is a vector $v$ that $A$ stretches but does not rotate: $A v = \lambda v$ . The scalar $\lambda$ is the eigenvalue.

Eigenvalues are how we diagnose the structure of a linear transformation. They tell us:

Whether a dynamical system is stable (all eigenvalues in the left half plane $\Leftrightarrow$ exponential decay)
Which directions a dataset varies the most along (PCA $\Leftrightarrow$ top eigenvectors of the covariance)
The natural frequencies of a vibrating system
The convergence rate of an iterative solver

A small classic application: principal component analysis (PCA), which finds the directions of largest variance in a dataset.

Two arrows, two eigenvalues, one of the most-used data-analysis methods in science.

Norms, ranks, and condition numbers

Norms quantify "size" of a vector or matrix; ranks count the number of independent directions; the condition number we met before is the ratio of largest to smallest singular value.

When cond(A) is large and your inputs only have a few correct digits, no algorithm can recover a high-precision answer for $x = A^{-1} b$ .

Least squares

When $A$ is not square — say it has more rows than columns — $A x = b$ generally has no exact solution. The least-squares solution minimizes $\|A x - b\|^2$ . This is how you fit any linear model to data.

The same call generalizes effortlessly: polynomial fits, multi- feature regressions, even non-linear-in-parameters fits via a nonlinear least squares routine in scipy.optimize we will meet later.

Sparse matrices

Many scientific problems produce matrices that are mostly zero: a 1-D heat equation matrix has three non-zeros per row, an image- processing matrix is block-tridiagonal, etc. Storing all those zeros is wasteful and slow. SciPy's sparse module provides CSR, CSC, COO, and LIL formats and dedicated solvers.

Sparse solvers let you handle systems with millions of unknowns on a laptop. PDE discretizations, graph Laplacians, and finite- element matrices are almost always sparse.

Practice: a tiny linear-algebra package

Implement the Thomas algorithm for solving a tridiagonal system $T x = d$ where $T$ has subdiagonal a (length $n-1$), diagonal b (length $n$), and superdiagonal c (length $n-1$).

The algorithm is two passes — a forward elimination and a back substitution — and runs in $O(n)$ time, vs $O(n^3)$ for a dense solve. It is the basis of every implicit 1-D PDE solver.

Return the solution x as a NumPy array of length $n$. Do not build the dense matrix.

Check your understanding

QuestionSelect one

When you call np.linalg.solve(A, b) for a $1000 \times 1000$ general matrix, what numerical algorithm is being executed under the hood?

Cramer's rule

Direct matrix inversion

LU factorization with partial pivoting, then forward and back substitution

A neural-network-based estimator

QuestionSelect one

A symmetric positive-definite matrix admits a Cholesky factorization $A = L L^T$ . Why is Cholesky preferred over general LU when applicable?

It is the only factorization that exists for SPD matrices

It requires roughly half the operations of LU, needs no pivoting (because it is unconditionally stable for SPD inputs), and uses half the storage

It allows complex eigenvalues

It is required to compute inverses

Linear Algebra Essentials

On this page