Simultaneous equations run the world

Why bother

At school you solved two equations in two unknowns, then three. Most people file it under things that were easy and never return to it.

Here is what that same method is doing today.

Every chip is simulated before it is built, by solving a system of a few million linear equations, once per step of simulated time. The two biggest companies selling that software, Synopsys and Cadence, earn about twenty billion dollars a year between them. No chip ships without it.
Google began as one system of linear equations with one unknown per web page. The answer was the order of the search results.
A neural network is a stack of matrices. Training one is a curve-fitting problem, which is a system of linear equations underneath. Running one is a matrix multiplied by a column of numbers, repeated a few hundred times.
Weather forecasts, bridge and aircraft design, MRI scanners and the power grid come down to the same problem with more unknowns.

The method is not the hard part. You already know it. The hard part is doing it with a million unknowns instead of two. This essay covers what the method is, where these huge systems come from, why size is the problem, and what was built to deal with it. Every term is defined as it comes up. You need nothing beyond school algebra.

The method you already know

Two equations, two unknowns:

\begin{aligned} 2x - y &= 1 \\ -x + 2y &= 0 \end{aligned}

The second says x = 2y⁠. Put that into the first: 3y = 1⁠, so y = 1/3 and x = 2/3⁠.

You made one kind of move: use one equation to remove an unknown from another. With three unknowns you make the same move more times. With ten you need paper. The method does not get harder as the system grows. It gets longer.

The figure below runs the method on a system of three to six unknowns. Some words first, so the figure makes sense.

Each row is one equation. The numbers in the row are the coefficients, the numbers multiplying each unknown, and the number after the vertical bar is the right-hand side. To remove the first unknown from row two, you subtract a multiple of row one from row two. Row one is then called the pivot row, and the coefficient you divide by to pick the multiple is the pivot. When every row below the first has a zero in the first column, you move on to row two as the new pivot row and clear the second column below it. When the bottom row has one unknown left, you read its value, then work upwards.

Elimination, one move at a time. The shaded row is the pivot row. The outlined row is the one being cleared. Zeros fill in below the diagonal, the line of cells running from top left to bottom right. When each row has one unknown left, the answers are read from the bottom up. Try six unknowns and watch the counter.

Two things to take from the figure.

Every step is the same step: one row minus a multiple of another. A machine that can do that move can do all of them, given a list of which row, which multiple, and when to stop. What made it tedious for you is bookkeeping, and bookkeeping is what computers are for.

The counter grows fast. Going from three unknowns to six does not double the work. It multiplies it by about eight. Each extra unknown means one more row to clear, one more cell in each row, and one more pivot to clear under. Three factors of n⁠, so the work grows like n^3⁠. For a million unknowns that is about 3 \times 10^{17} multiply-and-add operations, plus a table of a million million numbers to hold the coefficients. This is the wall. The schoolroom method hits it at a few thousand unknowns.

A grid of numbers

The letters x and y are labels. They mean “the first unknown” and “the second unknown”. Everything that matters is the coefficients and the right-hand sides. Keep the numbers, drop the letters, and keep the numbers in their places:

\begin{bmatrix} 2 & -1 \\ -1 & 2 \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \end{bmatrix}

The grid of coefficients is a matrix. A single column of numbers is a vector. Call the matrix A⁠, the vector of unknowns \mathbf{x}⁠, and the right-hand side \mathbf{b}⁠. The whole system is written A\mathbf{x} = \mathbf{b}⁠, and it is written the same way whether there are two unknowns or a million.

What does “A times \mathbf{x}” mean? Take a row of A⁠, lay it against the column \mathbf{x}⁠, multiply the pairs, add them up. That gives one number. Do it for each row and you get a new column. So the matrix is a machine: a column goes in, a column comes out. The figure runs that machine. Change the input until the output is (1, 0, 0)⁠.

A matrix multiplying a vector. Each row of A meets the column \mathbf{x}⁠, the pairs are multiplied and summed, and one number comes out. Solving A\mathbf{x} = \mathbf{b} means finding the input column that produces the output \mathbf{b}⁠.

Three consequences of writing the system this way.

The size no longer matters to the method. Elimination is “for each row, for each row below it, subtract a multiple.” Written once as a loop, it works for any n⁠.

The matrix and the right-hand side are separate. A records how the unknowns depend on each other. \mathbf{b} records what is being asked this time. If only \mathbf{b} changes, the work done on A can be reused.

You can multiply without solving. Feed in a guess, compare the output with \mathbf{b}⁠, adjust the guess, repeat. If one multiplication is cheap, this reaches the answer without doing elimination at all. This becomes important later.

Nobody writes a million equations by hand. They come from a simple rule applied at every point of some physical or logical structure. Here are three such structures.

Where the equations come from: a circuit

A circuit is wires meeting at junctions, called nodes. Charge cannot pile up at a node, so the current flowing in equals the current flowing out. That is Kirchhoff’s current law.

A resistor lets current through in proportion to the voltage difference across it. The constant of proportion is the conductance, G⁠, measured in siemens. A larger G means an easier path. The current from node i to node j through a resistor is G(v_i - v_j)⁠, where v_i is the voltage at node i⁠.

Write “in equals out” at every node, using that formula for each resistor touching it. You get one linear equation per node. The unknowns are the node voltages.

The figure has four resistors in a chain, with both ends connected to ground (zero volts) and one ampere of current pushed into the first node.

A circuit and its matrix. Moving dashes show current. The matrix is rebuilt from the sliders every time you move one. Changing the source changes only \mathbf{b}⁠. Changing a resistor changes only the matrix cells for the two nodes it touches.

With every conductance at 1 siemens, the three equations are

\begin{bmatrix} 2 & -1 & 0 \\ -1 & 2 & -1 \\ 0 & -1 & 2 \end{bmatrix} \begin{bmatrix} v_1 \\ v_2 \\ v_3 \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}

and the solution is v = (3/4,\ 1/2,\ 1/4) volts. Three quarters of the ampere goes straight to ground through the first resistor. One quarter goes the long way round.

The matrix follows a pattern. The cell on the diagonal for a node holds the total conductance touching that node. The cell in row i⁠, column j holds minus the conductance of the resistor between nodes i and j⁠. If there is no resistor between them, the cell is zero. Nothing had to be derived. Each resistor adds its conductance to four cells, and that is the whole matrix.

This is how chip simulation works. A chip has billions of transistors. Before the design is sent to a factory, where a production run costs hundreds of millions of dollars, the whole circuit is simulated. The simulator, descended from a 1973 Berkeley program called SPICE, writes one equation per node, treats each transistor as a resistor for the current instant of time, solves, moves time forward a fraction of a nanosecond, and solves again. Millions of nodes, millions of time steps.

One more thing to notice. A resistor touches two nodes, so a row of the matrix has an entry only for the nodes that share a wire with it. In a chip, a node has about five neighbours. A million nodes gives a matrix of a million rows and a million columns, a million million cells, of which only about five million are not zero. Keep that in mind.

Where the equations come from: Google

Picture a large crowd of people on the web, each on some page. Every tick of a clock, each person picks one of the links on their page at random and follows it. Let this run for a long time. Pages where many people end up are, by this measure, the important ones. That is PageRank, and it was the idea Google was built on.

The rule as an equation: the share of the crowd on page i after the next tick is the sum, over every page j that links to i⁠, of the share on j divided by the number of links on j⁠. Once the crowd has settled, the share after the tick equals the share before. That gives one linear equation per page.

Three pages. A links to B and C. B links to C. C links to A.

\begin{aligned} r_A &= r_C \\ r_B &= \tfrac{1}{2} r_A \\ r_C &= \tfrac{1}{2} r_A + r_B \end{aligned}

Ninety people following links at random. The dots are people. The bars are the exact shares after each tick, with the final answer marked. A and C end with two fifths each and B with one fifth. Check against the equations: r_B = \tfrac{1}{2}(0.4) = 0.2 and r_C = 0.2 + 0.2 = 0.4⁠.

As a matrix, with the cell in row i⁠, column j equal to 1/\text{links}(j) when j links to i and zero otherwise:

\mathbf{r} = M\mathbf{r}, \qquad M = \begin{bmatrix} 0 & 0 & 1 \\ \tfrac{1}{2} & 0 & 0 \\ \tfrac{1}{2} & 1 & 0 \end{bmatrix}

Moving everything to one side gives (I - M)\mathbf{r} = \mathbf{0}⁠, where I is the matrix with ones on the diagonal and zeros elsewhere. That is A\mathbf{x} = \mathbf{b} again. Google adds a small chance that a person gets bored and jumps to a random page, which makes the right-hand side nonzero and guarantees the crowd settles.

The crowd never solved an equation. It applied the rule over and over until the shares stopped changing. Each tick is one multiplication of the matrix by the current shares, the third consequence from earlier. And one multiplication is cheap, because a page has about ten links, so a row has about ten nonzero cells. One tick for the whole web is one pass down the list of links. Each row can be computed on a different computer, because no row needs another row’s result.

Eighteen pages and the matrix their links make. Each filled cell is one link. Press Play and the shares settle within a few dozen ticks. The matrix was never stored as a full grid or eliminated. It was only multiplied by, and multiplying touches only the filled cells.

In 1998 this was run on a few hundred million pages, on hardware two graduate students could afford. The results were better than every other search engine’s.

Where the equations come from: fitting a curve

You have measured points (x_k, y_k) and want a parabola through them: y = c_0 + c_1 x + c_2 x^2⁠. A parabola is curved, so where are the linear equations?

The unknowns are not x⁠. You measured every x_k⁠. The unknowns are the three coefficients c_0, c_1, c_2⁠, and the parabola is linear in those. Each measured point gives one equation:

c_0 + c_1 x_k + c_2 x_k^2 = y_k

The row of the matrix for that point is (1,\ x_k,\ x_k^2)⁠: the powers of a number you already have.

Three points give three equations, and the parabola passes through all three. Real data gives you a thousand points, and no parabola passes through all of them. So you ask instead for the parabola that misses them by the least. The miss at each point is called the residual. “Least” is measured as the sum of the squared residuals, which is why the method is called least squares.

Least squares. Drag a point. The dashed sticks are the residuals. The curve is the one with the smallest sum of squared residuals. Below it is the small square system that determines the coefficients. Switch to a line or a cubic: same method, different row.

With more equations than unknowns, there is no exact solution. Multiplying both sides by A^{\mathsf T}⁠, which is A with its rows and columns swapped, gives a square system with as many equations as coefficients:

A^{\mathsf T} A\, \mathbf{c} = A^{\mathsf T} \mathbf{y}

Solving that gives the least-squares coefficients. For a million data points and three coefficients, A has a million rows and three columns, but A^{\mathsf T} A is three by three. The data is read once and summarised into nine numbers.

Replace the powers of x with any quantities you like and this is linear regression, the workhorse of statistics. Then do one more thing: take the output, apply a simple bend (replace negative values with zero, for example), and feed the result into another matrix. Repeat several times. That is a neural network. Every layer in every current model is a matrix times a vector, plus a bend. Training is the fitting problem above with billions of coefficients instead of three. Answering a question is the machine from the second figure, run a few hundred times in a row.

Why size is the problem

In all three cases, a row came from a local rule. A node only knows its neighbours. A page only knows its links. A measured point only knows its own coordinates. Local rules produce rows that are mostly zeros. A matrix that is mostly zeros is called sparse. One with most cells filled is called dense.

Sparsity is the first way past the wall. If you only store the nonzero cells, memory grows with the number of connections rather than with the square of the number of unknowns.

Dense versus sparse. Grow the grid of nodes. The full matrix grows as the square of the node count, but the filled cells grow only in proportion to it: about five per row, whatever the size. The faded bars show what a million nodes would cost each way.

The next figure puts the pieces together. Choose a number of unknowns and a machine. It shows what elimination on the full matrix would cost, and what the sparse alternative costs: storing only the nonzero cells and multiplying by the matrix a few hundred times, as the crowd did.

The wall. Elimination on a dense matrix needs about n^3/3 operations and 8n^2 bytes of memory. The sparse method needs a few hundred passes over about five cells per row. As n grows, the dense method runs out of memory before it runs out of time, on every machine. A GPU does not make the dense method workable at a million unknowns. It makes the sparse method instant, and it makes dense work practical at around a hundred thousand, which is the size of a neural network layer's matrix.

Three facts from the figure.

Store only what exists. At a million unknowns, the full table is eight terabytes. No single machine holds that. The nonzero cells fit in tens of megabytes.

Repeat a cheap step instead of doing one expensive one. The crowd never eliminated anything. It multiplied by the matrix until the answer stopped moving. For large sparse systems this is the standard approach, with refinements that choose each new guess more cleverly than the plain rule does. The cost is the number of passes times the number of nonzero cells.

Rows do not wait for each other. When you multiply a matrix by a vector, row one does not need row two’s result. A million rows can be computed by a million workers at the same time.

What a GPU is

In the 1990s a graphics card had one job. For every corner of every triangle in a game scene, multiply its position by a four-by-four matrix to rotate it and project it onto the screen. Millions of corners, the same small matrix, sixty times a second. So the cards were built as hundreds of small arithmetic units doing the same multiply-and-add on different data at the same time. That is the third fact above, built into hardware.

People doing scientific computing noticed that this is exactly what multiplying a large matrix by a vector requires. In 2007 NVIDIA released CUDA, which let programmers use the card for any calculation, and the practice was named general-purpose GPU computing, or GPGPU. A current data-centre GPU performs on the order of sixty trillion multiply-and-add operations a second at the precision scientists use, and over a thousand trillion at the lower precision neural networks use. Parts of the chip, the tensor cores, do nothing except multiply matrices.

This is why a company that made game graphics is now worth more than the largest oil companies, banks and car makers. Chip simulation, search, weather and neural network training all have matrix arithmetic as their inner loop, and the GPU is the fastest machine ever built for it.

What to learn

The method is simple and you already know it. The notation A\mathbf{x} = \mathbf{b} makes it the same at every size. The systems write themselves from local rules. The businesses built on solving them are among the largest there are. The only real difficulty is scale, and the answers to scale are sparsity, repetition of a cheap step, and parallel hardware. None of this is deeper than substitution.

If you want to go further, this is the order that makes sense: what a matrix does to a vector; when A\mathbf{x} = \mathbf{b} has a solution and when it has exactly one; how elimination works and why it costs n^3⁠; how sparse matrices are stored; how repeated multiplication reaches the answer and what it costs; least squares and why A^{\mathsf T} A appears; and the condition number, which measures how much the answer moves when the inputs move slightly. If you tried the cubic in the curve-fitting figure and saw the coefficients jump, that was the condition number at work.

Once a problem has been written as rows of a matrix, the physics or the business logic is finished. What remains is linear algebra.