Solving Ax = b

"Two numbers add up to 3 and differ by 1. What are they?" That is a system of equations, the kind you solved in middle school. Solving problems like this is where linear algebra gets used the most. The only difference in real life is that there are not two equations but thousands, or millions.

Developers run into systems of equations all the time. Finding where two rays meet in a game, balancing forces in a circuit or physics simulation, undoing a filter in image processing, finding the weights in machine learning: each of these solves "several unknowns, several equations" all at once.

This lesson looks at a system of equations through two pictures. The row picture sees each equation as one line, and the column picture sees it as a linear combination from Lesson 4. You will also see that the determinant from Lesson 7 is the signal that tells you "is there exactly one solution?", and finally you will solve systems yourself with NumPy's np.linalg.solve.

8-1. One equation = one line (10 min)

Plot every point (x, y) that satisfies x + y = 3 and they line up along a single line: (3, 0), (2, 1), (0, 3), (−1, 4), and so on. In general, an equation of the form a x + b y = c is one line on the plane. x and y are never squared or multiplied together, only scaled by numbers and added, and that is where the name "linear" comes from.

Two equations means two lines. A point that satisfies both equations at the same time lies on both lines: it is where they cross. Seen as a picture, solving a system of equations means finding the point where two lines meet. This way of picturing each equation (each row of the matrix) on its own is called the row picture.

The lab below starts with x + y = 3 (teal) and x − y = 1 (crimson). The black point where the two lines cross is (2, 1), and indeed 2 + 1 = 3 and 2 − 1 = 1, so it satisfies both equations.

The sliders a, b and c on the right of the lab are the numbers in each equation a x + b y = c. Check a few things for yourself.

• Change only c and the line moves parallel to itself. How it tilts is decided by a and b.
• Change a or b and the line turns. If b is 0 you get a x = c, a vertical line; if a is 0 you get a horizontal line.
• The coordinates written next to the crossing point are the solution. The sliders move in steps of 0.5, so sometimes the solution shows up rounded, like 1.67.

The slider named b means "the number in front of y". It is not the same as the right-hand-side vector b you will meet in a moment. Books often use the same letter for both too, so tell them apart from context.

8-2. One solution, none, or infinitely many (7 min)

There are only three ways two lines can sit. The lab's 'Ready-made equations' include all three.

The two equations in 'No solution' have the same left side but different right sides, 2 and −2. "The same thing is 2 and −2 at once" cannot happen, so there is no solution. In 'Infinitely many', the second equation is just the whole first equation multiplied by 2. It carries no new information, so in effect there is only one equation. With only one equation, every point on that one line is a solution.

Also remember that you never get "exactly two solutions" or "exactly three". Two lines either meet at one point, never meet, or are the same line. Even with more unknowns, these three (one, none, infinitely many) are the only possibilities.

The two cases without exactly one solution have something in common: the two lines have the same slope. In the lab, pick 'No solution' and move only the c slider of either equation. The two lines move while staying parallel, and the moment c of the second equation reaches 2 they land on top of each other and the status changes to 'Infinitely many solutions: the two lines are the same line'.

8-3. The same equations as a matrix: Ax = b (6 min)

With many equations, writing out x and y every time gets tedious. So we pull out just the numbers and write them as a matrix and vectors. The numbers on the left side form the matrix $A$, the unknowns form the vector $\mathbf{x}$, and the right side is the vector $\mathbf{b}$.

Each row of the matrix is one equation. Compute $A\mathbf{x}$ with the "dot product with each row" view from Lesson 5: the dot product of the first row $(1, 1)$ with $(x, y)$ is $x + y$, and the dot product with the second row $(1, -1)$ is $x - y$, which gives back the original equations.

Read $A\mathbf{x} = \mathbf{b}$ as a transformation and it becomes this question: "Which vector $\mathbf{x}$ lands on $\mathbf{b}$ after going through the transformation $A$?" In Lesson 7, $A^{-1}$ was the matrix that undoes the transformation, so undoing $\mathbf{b}$ brings you back to the starting point $\mathbf{x}$. That thinking is right, but in real calculations you do not build the inverse. You will see why in 8-5.

8-4. The column picture: building b from the columns (8 min)

Now read the same equations down the columns. In Lesson 5, A @ v was "x·column 1 + y·column 2". So A x = b asks this:

This is exactly the linear combination problem from Lesson 4. When you moved the a and b sliders in the Lesson 4 lab to land on target T, you were solving a system of equations. This way of seeing the columns of a matrix as arrows is called the column picture.

In the lab above, pick 'One solution' and turn on 'Show the column picture (x·col 1 + y·col 2 = b)'. The teal arrow is column 1 (1, 1), the crimson arrow is column 2 (1, −1), and the thick black arrow is b = (3, 1). The dashed lines show "take column 1 twice, then attach column 2 once, and you arrive at b". The solution (x, y) = (2, 1) is exactly the set of coefficients of the linear combination.

Now read the three cases again with the column picture, and you can see why they happen.

• One solution: if the two columns point in different directions, they span the whole plane (Lesson 4). You can build any b, and there is only one way to do it.
• No solution: pick 'No solution' and look at the column picture. The two columns (1, 1) and (2, 2) lie on one line. With these two you can only move along that line, but b = (2, −2) is off the line.
• Infinitely many: in 'Infinitely many', the two columns (1, 2) and (−1, −2) also lie on one line, but this time b = (1, 2) is on that line. You can get there, but since the two columns point along the same line, there are endlessly many ways to mix them.

Put together, this ties back to the determinant from Lesson 7 like so.

When det = 0, whether you get "none" or "infinitely many" depends on whether b lies on the line the plane got squashed onto. In the next lesson this distinction gets its proper name: the column space.

8-5. Solving with NumPy: np.linalg.solve (12 min)

Follow the code in the "Check it with NumPy" panel below the lab, one line at a time. First make A and b. Each equation is one row, so write the equations from top to bottom, one inner list each.

>>> import numpy as np
>>> A = np.array([[1, 1], [1, -1]])
>>> b = np.array([3, 1])
>>> A.shape, b.shape
((2, 2), (2,))

A has 2 rows and 2 columns, and b is a vector with two components. The number of equations (rows) must match the number of components in b. Now check with the determinant that there is exactly one solution, and solve with np.linalg.solve(A, b). The arguments go in the order A first, then b.

>>> np.linalg.det(A)
np.float64(-2.0)
>>> x = np.linalg.solve(A, b)
>>> x
array([2., 1.])

You put in only integers, but the answer comes out as real numbers (float64) like 2. and 1.. Solutions to systems of equations are usually fractions, so solve always computes in floating point. Once you have an answer, always check it. If A @ x equals b, you solved it correctly.

>>> A @ x
array([3., 1.])
>>> np.allclose(A @ x, b)
True

Floating-point arithmetic can sneak in tiny errors (Lesson 3), so check with np.allclose instead of ==. You can confirm the row picture and the column picture in NumPy too. A[0] is the first row (the first equation), and A[:, 0] is the first column.

>>> A[0] @ x
np.float64(3.0)
>>> x[0] * A[:, 0] + x[1] * A[:, 1]
array([3., 1.])

The first line says "plug the solution into the first equation and you get 3" (row picture), and the second says "2 × column 1 + 1 × column 2 = b" (column picture).

To keep the screen tidy, the lab panel shows the determinant as np.round(np.linalg.det(A), 10), rounded to 10 decimal places. If you compute it yourself, you may get a value whose last digits are slightly off, like -1.9999999999999996 (Lesson 7).

Three unknowns work the same way

With three equations and three unknowns (x, y, z), the code does not change. As a picture it is three planes meeting at a point, but to NumPy it is just a 3×3 matrix.

 2x +  y -  z =   8
-3x -  y + 2z = -11
-2x +  y + 2z =  -3

>>> A3 = np.array([[2, 1, -1],
...                [-3, -1, 2],
...                [-2, 1, 2]])
>>> b3 = np.array([8, -11, -3])
>>> x3 = np.linalg.solve(A3, b3)
>>> x3
array([ 2.,  3., -1.])
>>> np.allclose(A3 @ x3, b3)
True

The solution is x = 2, y = 3, z = −1. By hand, adding and subtracting equations, this takes a few minutes. Even with 1000 unknowns, np.linalg.solve does it in the same single line.

No unique solution means an error

Pick 'No solution' or 'Infinitely many' and an error appears in the lab panel. NumPy calls a matrix with det 0 a singular matrix, and since solve cannot settle on a single answer, it raises LinAlgError.

>>> A_none = np.array([[1, 2], [1, 2]])
>>> b_none = np.array([2, -2])
>>> np.linalg.det(A_none)
np.float64(0.0)
>>> np.linalg.solve(A_none, b_none)
Traceback (most recent call last):
  ...
numpy.linalg.LinAlgError: Singular matrix

Put in the 'Infinitely many' case, A = [[1, -1], [2, -2]] and b = [1, 2], and you get exactly the same error. solve only handles the "exactly one solution" case, so it does not tell "none" apart from "infinitely many". You will learn how to tell them apart in Lesson 9, and how to find the "closest answer" when there is no solution in Lesson 13.

In a program, wrap the spot that can fail in try.

>>> try:
...     np.linalg.solve(A_none, b_none)
... except np.linalg.LinAlgError:
...     print("No unique solution")
No unique solution

Why solve instead of inv(A) @ b

In 8-3 we said the solution is A⁻¹b, so np.linalg.inv(A) @ b gives the same answer too.

>>> np.linalg.inv(A) @ b
array([2., 1.])

In practice, though, you almost always use solve. There are two reasons.

1. It is slower. Building the inverse is roughly like solving the system once per column (n times for an n×n matrix). You only need the answer for one b, so you would be computing information you never use.
2. Errors grow more easily. The rounding errors from building the inverse pile up once more in the multiplication. solve factors A directly and solves, so there are fewer steps.

If you want to compare them yourself, run the code below. Timings differ from computer to computer, so the results are not written here. Usually inv takes longer, and the leftover error (the length of A @ x - b) is smaller or about the same for solve.

import time
import numpy as np

rng = np.random.default_rng(0)
A = rng.standard_normal((2000, 2000))   # a random 2000×2000 matrix
b = rng.standard_normal(2000)

t = time.perf_counter(); x1 = np.linalg.solve(A, b); print("solve:", time.perf_counter() - t)
t = time.perf_counter(); x2 = np.linalg.inv(A) @ b;  print("inv  :", time.perf_counter() - t)

print(np.allclose(x1, x2))                                   # the answers are almost the same
print(np.linalg.norm(A @ x1 - b), np.linalg.norm(A @ x2 - b))  # compare the leftover errors

Unless you need the inverse itself (say, the matrix is small and you have to undo the same transformation many times), remember "solving A x = b means solve".

8-6. Common mistakes (3 min)

Swapping the arguments. Swap the order as in np.linalg.solve(b, A) and you get an error saying the first argument is one-dimensional, so it cannot be treated as a matrix.

>>> np.linalg.solve(b, A)
Traceback (most recent call last):
  ...
numpy.linalg.LinAlgError: 1-dimensional array given. Array must be at least two-dimensional

The number of equations and the length of b differ. With two equations but three numbers on the right side, you get an error saying the shapes do not match. The message is long, but you only need to read the last parenthesis, (size 3 is different from 2).

>>> np.linalg.solve(A, np.array([3, 1, 0]))
Traceback (most recent call last):
  ...
ValueError: solve1: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (m,m),(m)->(m) (size 3 is different from 2)

Copying numbers before tidying the equation. An equation written as x = 3 − y should first be rewritten as x + y = 3, with "unknowns on the left, constants on the right", before you copy it into the matrix. Write a missing unknown as 0. For example, 2x = 4 becomes the row [2, 0] with 4 on the right.

Exercises (5 min)

1. Build the system 2x + y = 4, x − y = 0.5 with the lab sliders and read off the crossing point. Then solve it with NumPy and check it.
2. In the lab, pick 'No solution' and move only one slider so that there is exactly one solution. Which slider did you move, and how did det change?
3. For x + 2y = 2, 2x + 4y = c, which c gives infinitely many solutions? Explain with the column picture.

Solutions

>>> A1 = np.array([[2, 1], [1, -1]])
>>> b1 = np.array([4, 0.5])
>>> x1 = np.linalg.solve(A1, b1)
>>> x1
array([1.5, 1. ])
>>> A1 @ x1
array([4. , 0.5])

For question 2, you have to move the a or b slider, not c. c only shifts a line parallel to itself, so the slopes stay the same. For example, change b of the second equation from 2 to 1.5 and det = 1×1.5 − 2×1 = −0.5, which is no longer 0, and the two lines meet at one point.

For question 3, c = 4. The two columns (1, 2) and (2, 4) lie on one line, and the only place you can reach is the line in the direction of (1, 2). For b = (2, c) to be on that line, it must be (2, 4). Then the second equation is twice the first and carries no new information.

In the next lesson you will sort out, when det = 0 and the plane gets squashed, where things go and what disappears, under the names column space and null space. That picture also explains why there are "infinitely many" solutions.