Linear Algebra and Matrix Operations

As soon as we start working with data in multiple dimensions, we step into the world of linear algebra. Linear algebra is the branch of mathematics that studies vectors, matrices, and linear transformations — in other words, how sets of numbers relate to one another through operations like addition, scaling, and multiplication.

It forms the mathematical backbone of nearly every modern data science method — from machine learning and image analysis to network mapping and text similarity measurement.

What Is a Matrix?

A matrix is simply a rectangular grid of numbers arranged in rows and columns. In NumPy, we represent matrices using two-dimensional arrays. Each row can represent an observation (like a book or an author), and each column can represent a feature (like word frequency, sentiment, or publication year).

Let’s create two small matrices to work with:

import numpy as np

A = np.array([[1, 2],
              [3, 4]])

B = np.array([[5, 6],
              [7, 8]])

print("Matrix A:\n", A)
print("Matrix B:\n", B)

Matrix A:
 [[1 2]
 [3 4]]
Matrix B:
 [[5 6]
 [7 8]]

Here, A and B each have two rows and two columns. They could represent, for example, two time periods of cultural activity across two regions, where each number records a frequency or measurement.

Matrix Multiplication

In linear algebra, matrix multiplication is one of the most fundamental operations. It allows us to combine information from two matrices to form a new one.

Unlike adding or subtracting matrices (which works element by element), multiplication involves combining rows from the first matrix with columns from the second.

In NumPy, we perform this using np.dot(A, B) or the shorthand A @ B.

C = np.dot(A, B)
print("Matrix C (A multiplied by B):\n", C)

Matrix C (A multiplied by B):
 [[19 22]
 [43 50]]

What happened?

To calculate the element at row 1, column 1 (value 19):

$$1\cdot5+2\cdot7=19$$

To calculate row 2, column 1 (value 43):

$$3\cdot5+4\cdot7=43$$

Matrix multiplication is almost everywhere in data science. It’s used to combine datasets, transform coordinates in image processing, or project data into lower-dimensional spaces during, e.g., Principal Component Analysis (PCA) - a statistical method for simplifying complex datasets. PCA reduces large amounts of numerical data into just a few key dimensions (called principal components) that capture the most important patterns or variations. In humanities contexts, you might use it to compute document similarities, model cultural relationships, or apply transformations to text or metadata features.

Matrix Determinant

The determinant is a single number that summarizes key properties of a matrix. It tells us whether a matrix can be inverted (that is, whether it has a unique reverse operation), and in geometric terms, it represents how the matrix scales or transforms space.

In NumPy, we can calculate it using np.linalg.det():

det_A = np.linalg.det(A)
print(det_A)

-2.0000000000000004

A determinant of zero means the matrix is singular, meaning it cannot be inverted. A nonzero determinant (like -2 here) means it can be. For a two-by-two matrix, the absolute determinant is the area of the parallelogram spanned by its row vectors. Its sign records orientation. If that area is zero, the vectors lie on the same line — meaning there’s no “volume” or space between them, so the system of equations they represent is dependent and can’t be uniquely solved.

For a square matrix, a nonzero determinant means its columns are linearly independent in exact arithmetic. Linear independence is different from statistical independence, and it does not establish predictive usefulness. Floating-point rank and conditioning require additional care; a tiny rounding residual or a nonzero computed determinant alone is not a quality certificate.

Matrix Inversion

Matrix inversion is the process of finding another matrix that “undoes” the effect of the first one. Just as dividing by a number reverses multiplication (for example, 8 ÷ 2 = 4), multiplying by a matrix’s inverse reverses the transformation that matrix performs.

If A is a matrix, its inverse is written as A⁻¹, and when we multiply them together, we get the identity matrix — a matrix that behaves like multiplying by 1:

$$A \times A^{-1} = I$$

Let’s see this in action:

inv_A = np.linalg.inv(A)
print("Inverse of A:\n", inv_A)

Inverse of A:
 [[-2.   1. ]
 [ 1.5 -0.5]]

Now, if we multiply A by its inverse:

I = np.dot(A, inv_A)
print("A multiplied by its inverse:\n", I)

A multiplied by its inverse:
 [[1.0000000e+00 0.0000000e+00]
 [8.8817842e-16 1.0000000e+00]]

The result is approximately the identity matrix. The tiny off-diagonal value is floating-point roundoff, not an extra relationship in the data. Check a suitable tolerance with np.allclose(A @ inv_A, np.eye(2)) rather than requiring exact equality.

Why is this important?

Solving linear systems is a common analytical task. An inverse explains the algebra, but a numerical solver can obtain a solution without explicitly forming that inverse. For instance, in regression analysis, we often solve equations of the form:

$$Y=X\beta+\varepsilon$$

Here, Y represents the output (such as observed data), X represents the input variables, and β (beta) contains the coefficients we want to estimate. Here $\varepsilon$ represents residual differences between observations and the chosen linear representation. Least squares chooses coefficients to minimize $\lVert Y-X\beta\rVert_2^2$. The inverse formula below explains the full-column-rank case; the implementation uses np.linalg.lstsq, as documented in the NumPy least-squares reference.

Simple Example

Consider three invented novels used to demonstrate the arithmetic of pages, characters and reading time. These are teaching values, not a measured reading study. We fit a model without an intercept for this example; in real work, that choice needs a substantive reason.
This relationship can be modeled as:

$$Y = X\beta$$

where:

• Y represents the observed reading times for several novels,
• X contains the input features (pages and characters),
• β contains coefficients describing the fitted linear association, not causal effects.

Using a Small Dataset

We can express this in matrix form:

$$X = \begin{bmatrix} 200 & 50 \\ 150 & 40 \\ 300 & 80 \end{bmatrix}, \quad Y = \begin{bmatrix} 5 \\ 4 \\ 8 \end{bmatrix}$$

We want to find β (the coefficients):

$$\beta = \begin{bmatrix} b_1 \\ b_2 \end{bmatrix}$$

so that

$$Y = X\beta$$

Solving for β Using Matrix Algebra

Consider the simple linear model again:

$$Y = X\beta$$

where: - y is a column vector (the observed outcomes), - X is a matrix containing the input data (the independent variables), - β is a column vector of coefficients we want to find.

Our goal is to solve for β.

Step 1: Express the equation

$$Y = X\beta$$

To isolate β, we need to “undo” the effect of X.
However, X is not always a square matrix, so we can’t directly divide by it.
Instead, we multiply both sides by the transpose of X, written $X^T$.

Step 2: Multiply both sides by $X^T$

$$X^T Y = X^T X \beta$$

Here:

• $X^T$ means the transpose of X, obtained by flipping rows into columns.
  For example, if $X$ is a $3 \times 2$ matrix, then $X^T$ is a $2 \times 3$ matrix.

Step 3: Multiply by the inverse of $X^T X$

If $X^T X$ is invertible, we can multiply both sides by its inverse to isolate $β$:

$$(X^T X)^{-1} X^T Y = (X^T X)^{-1} X^T X \beta$$

Since $(X^T X)^{-1} X^T X$ equals the identity matrix $I$, which leaves vectors unchanged, we simplify to:

$$\beta = (X^T X)^{-1} X^T Y$$

Step 4: Interpretation of notation

Step 5: The Final Formula

$$\boxed{\beta = (X^T X)^{-1} X^T Y}$$

This is the ordinary least-squares coefficient formula when $X$ has full column rank. More generally, noisy observations need not satisfy $Y=X\beta$ exactly: minimizing the squared residuals yields the normal equations $X^T X\hat\beta=X^T Y$. The inverse formula is not a recommendation to compute that inverse explicitly.

Example Meaning

If $X$ contains features such as the number of pages and the number of characters in a novel,
and $Y$ contains the reading times, then solving for $β$ gives the coefficients (weights)
that describe how much each feature contributes to the reading time.

Python Example

import numpy as np

# Input data (pages and characters)
X = np.array([[200, 50],
              [150, 40],
              [300, 80]])

# Observed reading times
Y = np.array([[5],
              [4],
              [8]])

# Compute least-squares coefficients without forming an explicit inverse
beta = np.linalg.lstsq(X, Y, rcond=None)[0]
print("Estimated coefficients (β):\n", beta)

Estimated coefficients (β):
 [[-5.00366212e-17]
 [ 1.00000000e-01]]

The computed coefficients are approximately 0 hours per page and 0.1 hours per thousand characters in this fitted, no-intercept example. Indeed, every supplied reading time equals 0.1 times its character count. The tiny negative page coefficient is numerical roundoff. These three invented observations establish neither a causal reading law nor that pages are unimportant in a different dataset.

In Summary

Matrix inversion lets us “reverse” a system of relationships to find the unknowns — just as dividing by a number undoes multiplication. For computation, distinguish an exactly solvable square system (np.linalg.solve) from a least-squares fit (np.linalg.lstsq). In either case, inspect the residual and the assumptions before interpreting the coefficients.

Check the fitted result before interpreting it

Predict the fitted reading times from the character column alone. Then compare them with the matrix calculation.

predicted = X @ beta
print(predicted.ravel())
print(np.linalg.matrix_rank(X))
assert np.allclose(beta.ravel(), [0.0, 0.1])
assert np.allclose(predicted, Y)
assert np.linalg.matrix_rank(X) == 2

[5. 4. 8.]
2

The predicted values are 5, 4 and 8 hours. A perfect fit here follows from the numbers we constructed; it is not evidence of generalization. Change one reading time, refit, and explain why exact reconstruction may disappear. Retain the units when describing each coefficient.

Chapter exercisesWorked answers after an attempt