Linear algebra

Linear algebra is the maths of scaling things and adding them up. Every layer of every AI model is that one move, repeated at enormous scale.

By Pranay Mahendrakar
18 min read
3 reading levels
Published 20 Aug 2026

Read these first

On this page 8

Why you should care
The word "linear" means something specific
The three jobs it does
A picture of the whole thing
The part that is genuinely awkward
Where you have already used it
Remember this
What to learn next

One lesson, three depths. Pick the one that fits you today — you can switch any time.

Beginner — No maths. Plain English.

Linear algebra is the maths of scaling things and adding them up.

Think about the bill at a kirana shop. Each item has a rate. The shopkeeper takes how much you bought of each thing, scales it by the rate, and adds the lot together.

You have done that arithmetic in your head your whole life. That is the entire operation linear algebra studies. The rest of it is that same move, done to millions of numbers at once.

Why you should care

Inside an AI model there is no thinking, no reasoning and no cleverness. There is scaling and adding, repeated a great many times.

A model gives each measurement an importance rating. It scales the measurement by that rating. Then it adds everything up into one number. Do it a thousand times over and you have one layer of a neural network.

So when you understand scale-and-add, you understand what a model physically does. Everything after that is bookkeeping about how many times and in what order.

The word "linear" means something specific

Linear means two promises hold, and they are promises you already expect from a shop bill.

Double the basket and the bill doubles. Buy twice as much of everything and you pay exactly twice as much. No bulk discount, no extra charge.

Two baskets billed together cost the same as billed apart. Splitting a shopping trip in half changes nothing about the total.

Anything that keeps both promises is linear. Anything that breaks either one is not.

This is a real restriction, and it is the reason models need something extra. A bulk discount breaks the first promise. So does a tax that starts only above a certain amount. Real life is full of such bends, and pure scale-and-add cannot produce a single one of them.

That is precisely why a neural network puts a small bend between its layers. Without the bend, stacking layer upon layer buys you nothing at all. There is a lesson on those bends: activation functions.

The three jobs it does

COMBINE     take many measurements, scale each one, add them up
            → one score, one prediction, one layer of a network

TRANSFORM   feed a list of numbers in, get a different list out
            → rotating an image, shrinking a big list into a small one

SOLVE       work backwards: you see the total, find the amounts that made it
            → learning the rates from past bills, which is what training does

The third job is the surprising one. Most of school teaches you the forward direction. AI runs it backwards.

You are handed thousands of finished bills. Nobody tells you the rates. Linear algebra is what lets you recover the rates from the bills. That recovery is the mathematical heart of training a model.

A picture of the whole thing

   how much you bought           the rates                the total
   ------------------           ----------              ----------
   rice, atta, dal, oil   scaled by   per-kilo price   →   one bill

   the same shape, inside a model:

   your measurements      scaled by   importance knobs  →   one score

Change the rates, and the same basket gives a different total. Training a model means hunting for the set of rates that make the totals come out right.

The part that is genuinely awkward

Working backwards does not always have an answer, and that has a plain physical reason.

Suppose a shop sells two ready-made packets. One packet is exactly double the other — same recipe, larger size. Now somebody hands you a finished mix and asks how many of each packet went in.

You cannot say. Two packets of the small one, or one of the large, would look identical in the final mix. The information was never there to recover.

Data does this constantly. Say two of your columns hold height in feet and height in inches. They carry one fact between them. The model cannot tell which of the two mattered. The maths does not fail loudly. It gives you an answer that is unstable, and that changes wildly when your data changes a little.

There is a name for how much genuinely independent information a table of numbers carries. It gets measured, and low values are a warning sign. That measure shows up in the Developer block.

Where you have already used it

The equaliser sliders in a music app. Each slider scales one band of sound, and the bands add back together.
Splitting a restaurant bill by what each person ordered.
A photo brightness slider, which scales every dot in the picture by the same amount.
Recommendations, where your ratings get scaled by how similar other people are to you, and added up.

Remember this

Linear means scale each thing and add: double the input, double the output.
AI models do nothing but scale-and-add, with a small bend placed between the layers.
Training runs the operation backwards, recovering the rates from finished totals.

What to learn next

Vectors and matrices — the two objects that hold all these numbers.
Calculus — how a model works out which way to change a rate.
Activation functions — the bend that makes stacked layers worth stacking.

Developer — Code and libraries.

Setup

bash

pip install numpy

Every array below is typed inline. Nothing downloads, and the whole page runs on a CPU in under a second.

This block takes the three jobs from the Beginner section in order. Solve backwards, fit noisy data, then stack layers and watch them collapse.

The objects themselves live in vectors and matrices. Dot products and shape rules are covered there, so this page assumes them.

Job one: solving backwards

solve.py

import numpy as np

# Two ready-made packets, each a fixed mix of rice and dal (units per packet).
#                    packet A   packet B
mix = np.array([[2.0,      1.0],      # rice
                [1.0,      3.0]])     # dal
need = np.array([8.0, 9.0])           # you want exactly 8 rice, 9 dal

counts = np.linalg.solve(mix, need)   # which combination gives exactly that?
print("packets to buy :", np.round(counts, 3))
print("check          :", mix @ counts)

# Now a pair of packets that are secretly the same recipe, doubled.
bad = np.array([[2.0, 4.0],
                [1.0, 2.0]])
print("rank of mix    :", np.linalg.matrix_rank(mix))
print("rank of bad    :", np.linalg.matrix_rank(bad))
try:
    np.linalg.solve(bad, need)
except np.linalg.LinAlgError as err:
    print("solve(bad)     :", err)

Output

packets to buy : [3. 2.]
check          : [8. 9.]
rank of mix    : 2
rank of bad    : 1
solve(bad)     : Singular matrix

Three packets of A and two of B. The check line confirms it lands exactly on the target.

The second matrix is the awkward case from the Beginner block, made concrete. Column two is column one doubled. Rank is the number of genuinely independent columns. Here it is 1, not 2. One column carries no new information, so NumPy refuses. Singular matrix is the message it uses.

Rank is worth checking on real feature tables. If a table has 40 columns and rank 37, three of your features are exact combinations of the others.

Job two: fitting when the numbers are noisy

Exact solving is a luxury. Real measurements disagree slightly, so there is usually no exact answer and you want the closest one.

fit.py

import numpy as np

# Four flats. Columns: number of fans, number of tube lights, and a 1 for the fixed charge.
usage = np.array([[2.0, 3.0, 1.0],
                  [1.0, 5.0, 1.0],
                  [4.0, 2.0, 1.0],
                  [3.0, 6.0, 1.0]])
bill = np.array([566.0, 604.0, 672.0, 858.0])   # rupees, read off a slightly noisy meter

rates, residual, rank, _ = np.linalg.lstsq(usage, bill, rcond=None)
print("recovered rates :", np.round(rates, 2))
print("squared error   :", np.round(residual, 2))
print("predicted bills :", np.round(usage @ rates, 1))

Output

recovered rates : [ 90.49  68.15 176.2 ]
squared error   : [39.02]
predicted bills : [561.6 607.4 674.4 856.5]

The true rates were 90 per fan, 70 per light and a fixed charge of 170. Each bill was then nudged by a few rupees. From four noisy readings, lstsq recovered 90.49, 68.15 and 176.2.

Look at what you have written. That is linear regression. There is no iteration, no learning rate and no training loop. lstsq computes the best answer in closed form.

That column of 1.0 values is the bias column. It gives the fit a baseline that does not depend on any input.

The residual is the total squared miss across all four flats. It is not zero, and it should not be. A model that hit noisy readings exactly would be fitting the meter's errors, which is overfitting.

Job three: why stacked linear layers collapse

collapse.py

import numpy as np

W1 = np.array([[0.5, -1.0, 2.0],
               [1.5,  0.0, -0.5]])      # 2 features in, 3 out
W2 = np.array([[1.0, -2.0],
               [0.5,  1.0],
               [-1.0, 0.5]])            # 3 in, 2 out

x = np.array([1.0, 2.0])
combined = W1 @ W2                      # one 2-in, 2-out layer that does the same job
print("through two layers :", x @ W1 @ W2)
print("through one layer  :", x @ combined)
print("combined matrix    :\n", combined)
print("rank of combined   :", np.linalg.matrix_rank(combined))

# Squeeze the middle down to a single number and the whole layer loses capacity.
narrow = np.array([[0.5], [1.5]]) @ np.array([[1.0, -2.0]])
print("rank through a 1-wide middle :", np.linalg.matrix_rank(narrow))

Output

through two layers : [ 2.  -7.5]
through one layer  : [ 2.  -7.5]
combined matrix    :
 [[-2.   -1.  ]
 [ 2.   -3.25]]
rank of combined   : 2
rank through a 1-wide middle : 1

Identical outputs. Two layers with nothing between them are one layer in a costume. The costume is a single 2-by-2 matrix, printed above.

Stack a hundred such layers and you still get one matrix. The parameter count grows, the training time grows, the expressive power does not move at all. This is the proof promised in what is a neural network, runnable in five lines.

The last line is the other half of the story. A narrow middle layer caps the rank of everything downstream. It does that however wide the neighbouring layers are.

That cap is a limit, and it is also a design tool. LoRA inserts a low-rank pair on purpose. It fine-tunes a huge model using a tiny count of trainable numbers.

Line by line, the parts that are not obvious

np.linalg.solve(A, b) rather than np.linalg.inv(A) @ b. Both give an answer for a well-behaved matrix. solve factorises and back-substitutes, which is roughly three times cheaper and noticeably more accurate. Forming an explicit inverse is almost always the wrong instinct.

rcond=None in lstsq. Without it, NumPy prints a FutureWarning and uses an old, machine-dependent cutoff for tiny singular values. Passing None selects the current, sane default. Copy it every time.

lstsq returns four things. Solution, residual, rank and singular values. The rank is free diagnostic information. If it comes back below your column count, some features are redundant. The recovered rates are then not trustworthy.

The residual array can come back empty. It is only populated when the system is overdetermined and full rank. With more columns than rows you get array([], dtype=float64), not a zero. Code that indexes residual[0] blindly will crash on exactly the data you did not test with.

x @ W1 @ W2 groups left to right, so it computes (x @ W1) @ W2. Matrix products are associative, so x @ (W1 @ W2) gives the same numbers. The cost differs enormously though — see the note on ordering below.

Common mistakes

1. Inverting a matrix out of habit. np.linalg.inv(A) @ b is slower and less numerically stable than np.linalg.solve(A, b). For a badly conditioned matrix, the inverse route can lose most of your significant digits without any warning printed.

2. Ignoring a near-singular matrix. matrix_rank uses a tolerance, so a matrix can pass as full rank while being numerically hopeless. Check the condition number too:

python

import numpy as np

nearly_bad = np.array([[2.0, 4.0], [1.0, 2.000001]])
print("rank      :", np.linalg.matrix_rank(nearly_bad))
print("condition : {:.1e}".format(np.linalg.cond(nearly_bad)))

Output

rank      : 2
condition : 1.3e+07

Rank says everything is fine. The condition number says an error in the seventh digit of your input can swamp the answer. Trust the second number.

3. Multiplying in a costly order. For matrices sized (1000, 5), (5, 1000) and (1000, 1), grouping as A @ (B @ C) costs about 10 thousand multiply-adds. Grouping as (A @ B) @ C costs about 6 million. Same answer, 600 times the work. This is why low-rank layers are written to keep the small dimension in the middle.

4. Forgetting the column of ones. Fit without it and you force the line through the origin. Here the true baseline is about 176 rupees. That one omission wrecks every coefficient, and the code still runs happily.

Try it yourself

Add a fifth flat to usage. Make its fan and light counts the sum of two existing rows, with a matching bill. Check whether rank in the lstsq output changes. Then work out why it does not.

Then take collapse.py and put np.maximum(0, ...) between the two layers. Try to find a single 2-by-2 matrix that reproduces its output for every input. You cannot, and failing to find one is the point of the exercise.

What to learn next

Vectors and matrices — shapes, dot products and similarity in detail.
Linear regression — the lstsq fit above, dressed as a model.
Derivatives and gradients — how models fit when no closed form exists.

Researcher — Mathematics and papers.

Span, basis, rank

A linear combination of vectors v_1 ... v_k with scalar coefficients c_1 ... c_k is sum_i c_i * v_i. The span is the set of all such combinations, and it is always a subspace.

A basis is a linearly independent spanning set. Every vector in the space has exactly one representation in a given basis, which is what makes coordinates meaningful.

For a matrix A of shape (m, n):

Column space C(A) is the span of the columns, a subspace of R**m. Its dimension is the rank r.
Null space N(A) is {x : A @ x = 0}, a subspace of R**n with dimension n - r.
Rank-nullity: dim C(A) + dim N(A) = n.

The system A @ x = b has a solution exactly when b is in C(A). It has a unique solution only when N(A) = {0}. The singular example in the Developer block has a one-dimensional null space. Any solution would therefore arrive as an infinite family. The requested b was not in the column space at all.

Strang's "four subspaces" adds the row space C(A.T) of dimension r and the left null space N(A.T) of dimension m - r. The row space and null space are orthogonal complements in R**n, which is the geometric statement behind least squares.

Least squares

Overdetermined systems, m > n, generally have no exact solution. The least-squares problem is

x_hat  =  argmin_x  norm( A @ x - b )**2

Setting the gradient to zero gives the normal equations:

A.T @ A @ x_hat  =  A.T @ b

Geometrically, the residual b - A @ x_hat is orthogonal to C(A). The fitted values A @ x_hat = P @ b where P = A @ inv(A.T @ A) @ A.T is the orthogonal projector onto the column space. P is symmetric, idempotent (P @ P = P), and has trace equal to r. That trace is the effective number of parameters in a linear fit. It reappears in AIC and in degrees-of-freedom corrections.

Do not solve the normal equations directly. Conditioning squares: kappa(A.T @ A) = kappa(A)**2. A matrix with kappa(A) = 1e4 is workable in float64 and hopeless after squaring in float32. np.linalg.lstsq uses an SVD-based driver (LAPACK gelsd), and scipy.linalg.lstsq exposes QR-based alternatives. QR via Householder reflections costs about 2 m n**2 - (2/3) n**3 flops and keeps kappa(A) rather than squaring it.

Ridge regression modifies the normal equations to (A.T @ A + lam * I) @ x = A.T @ b, with lam > 0. This shifts every eigenvalue of A.T @ A up by lam, bounding the condition number by (s_max**2 + lam) / lam. Regularisation is a numerical-stability device before it is a statistical one.

Eigenvalues

For a square A, a non-zero v with A @ v = lambda * v is an eigenvector, and lambda its eigenvalue. The map does not rotate such a v; it only stretches it.

When A has n linearly independent eigenvectors it diagonalises as A = V @ D @ inv(V), with D diagonal. Then A**k = V @ D**k @ inv(V), so repeated application is governed entirely by the eigenvalues.

That single fact explains a family of practical behaviours:

Vanishing and exploding gradients. Backpropagation through T timesteps of a recurrent network applies a Jacobian repeatedly. If the spectral radius rho(A) = max_i abs(lambda_i) is below 1, the product decays geometrically. Above 1, it explodes. See RNN and LSTM.
Power iteration. Repeatedly applying A and normalising converges to the dominant eigenvector at rate abs(lambda_2 / lambda_1)**k. PageRank is power iteration on a stochastic matrix. So is the spectral-norm estimate inside spectral normalisation for GANs (Miyato et al., 2018), which runs a single power-iteration step per training step.
Optimiser behaviour. Gradient descent on a quadratic with Hessian H converges at a rate governed by kappa = lambda_max / lambda_min. That is the origin of the condition-number dependence discussed in optimization.

Symmetric matrices are the well-behaved case. The spectral theorem guarantees a real orthogonal eigendecomposition A = Q @ D @ Q.T for any real symmetric A. Covariance matrices, Hessians and Gram matrices are all symmetric, which is why they are the ones people actually decompose.

A is positive semi-definite when x.T @ A @ x >= 0 for all x, equivalently when every eigenvalue is non-negative. A twice-differentiable function is convex exactly when its Hessian is PSD everywhere. PCA is the eigendecomposition of a mean-centred covariance matrix, and it is PSD by construction.

Eigenvalues are defined only for square matrices. The SVD generalises the idea to any shape. It is covered in vectors and matrices. For symmetric PSD A, the two coincide.

Cost, for the shapes you will actually meet

Operation	Cost	Note
`A @ B`, `(m,k) @ (k,n)`	`O(m k n)`	the dominant cost in every network
`solve`, LU on `(n,n)`	`O(n**3 / 3)` multiply-adds	plus `O(n**2)` per extra right-hand side
QR, Householder, `(m,n)`	`O(2 m n**2)`	stable, the workhorse for least squares
Full SVD `(m,n)`, `m>=n`	`O(m n**2)` with a large constant	avoid on tall dense data if a rank-`r` truncation will do
Symmetric eigendecomposition `(n,n)`	`O(n**3)`	`eigh`, not `eig` — faster and returns real values
Power iteration, `k` steps	`O(k * nnz(A))`	the only viable route at web scale

Two practical notes. Use np.linalg.eigh for symmetric input. eig returns complex dtypes even when the imaginary parts are zero, and downstream comparisons then behave oddly.

Randomised SVD (Halko, Martinsson and Tropp, 2011) computes a rank-r approximation in O(m n log r + (m + n) r**2). That is how PCA is done on anything large.

Where this shows up in modern systems

The forward pass of a transformer is a sequence of linear maps with non-linearities interleaved. Attention scores are Q @ K.T, a Gram-like matrix of inner products. Feed-forward blocks are two dense layers. Quantisation, pruning and low-rank adaptation are all statements about the numerical structure of those matrices.

Tensor and pipeline parallelism partition a single matrix product across devices. The partition is chosen so the communication pattern matches the interconnect. Whether a model fits your hardware is a linear-algebra question first.

References

Strang, G. Introduction to Linear Algebra, 6th ed., Wellesley-Cambridge Press, 2023. The four-subspaces framing, and MIT course 18.06.
Strang, G. "The Fundamental Theorem of Linear Algebra." American Mathematical Monthly 100(9), 848–855, 1993.
Trefethen, L. N., Bau, D. Numerical Linear Algebra. SIAM, 1997. Lectures 11 to 20 for least squares and conditioning.
Golub, G. H., Van Loan, C. F. Matrix Computations, 4th ed., Johns Hopkins University Press, 2013.
Halko, N., Martinsson, P.-G., Tropp, J. A. "Finding Structure with Randomness." SIAM Review 53(2), 217–288, 2011.
Miyato, T., Kataoka, T., Koyama, M., Yoshida, Y. "Spectral Normalization for Generative Adversarial Networks." ICLR, 2018.
Deisenroth, M. P., Faisal, A. A., Ong, C. S. Mathematics for Machine Learning. Cambridge University Press, 2020. Chapters 2 to 4.

What to learn next

Vectors and matrices — SVD, low-rank structure and numerical caveats.
Optimization — how the eigenvalues of the Hessian set your learning rate.
Statistics — what the residual from a least-squares fit is telling you.

🟢 Beginner — No maths. Plain English.

Why you should care#

The word "linear" means something specific#

The three jobs it does#

A picture of the whole thing#

The part that is genuinely awkward#

Where you have already used it#

Remember this#

What to learn next#

🟡 Developer — Code and libraries.

Setup#

Job one: solving backwards#

Job two: fitting when the numbers are noisy#

Job three: why stacked linear layers collapse#

Line by line, the parts that are not obvious#

Common mistakes#

Try it yourself#

What to learn next#

🔴 Researcher — Mathematics and papers.

Span, basis, rank#

Least squares#

Eigenvalues#

Cost, for the shapes you will actually meet#

Where this shows up in modern systems#

References#

What to learn next#

Beginner — No maths. Plain English.

Why you should care

The word "linear" means something specific

The three jobs it does

A picture of the whole thing

The part that is genuinely awkward

Where you have already used it

Remember this

What to learn next

Developer — Code and libraries.

Setup

Job one: solving backwards

Job two: fitting when the numbers are noisy

Job three: why stacked linear layers collapse

Line by line, the parts that are not obvious

Common mistakes

Try it yourself

What to learn next

Researcher — Mathematics and papers.

Span, basis, rank

Least squares

Eigenvalues

Cost, for the shapes you will actually meet

Where this shows up in modern systems

References

What to learn next