Matrices and systems

Why this matters

Real problems rarely have one unknown. A truss has a force in every member, an electrical network has a current in every branch, a chemical equation has a coefficient for every species, and a regression model has a weight for every feature. Each gives a system of linear equations, and matrices are how engineers, scientists and software solve them, often with thousands of unknowns at once.

What a matrix is

A matrix is a rectangular grid of numbers with m rows and n columns (an m × n matrix):

A = | 2   1 |
    | 1   3 |

A vector can be thought of as a matrix with one column.

Matrix times vector

To multiply A by a vector x, take the dot product of each row of A with x:

| 2  1 | | x |   | 2x + y  |
| 1  3 | | y | = | x + 3y  |

So the whole system

2x +  y = 5
 x + 3y = 10

is one equation: A x = b, with b = (5, 10). The matrix holds the coefficients, x holds the unknowns, and b holds the right-hand sides.

Matrix times matrix

To multiply A (m × n) by B (n × p), entry (i, j) of the result is row i of A dotted with column j of B. The inner sizes must match, and the result is m × p. Order matters: in general AB ≠ BA.

| 1  2 | | 5  6 |   | 1·5 + 2·7   1·6 + 2·8 |   | 19  22 |
| 3  4 | | 7  8 | = | 3·5 + 4·7   3·6 + 4·8 | = | 43  50 |

Solving A x = b

By hand, elimination works: subtract 2 × (second equation) from the first to get −5y = −15, so y = 3 and x = 1. Check: 2·1 + 3 = 5 and 1 + 3·3 = 10.

Conceptually, if A has an inverse A⁻¹ (the matrix that undoes A), then x = A⁻¹b. For a 2 × 2 matrix:

A = | a  b |      A⁻¹ = 1/(ad − bc) · |  d  −b |
    | c  d |                          | −c   a |

The number ad − bc is the determinant. If it is 0 the matrix has no inverse: the equations are either contradictory or redundant, so there is no unique solution.

In practice you never compute an inverse to solve a system. You call a solver that uses elimination (LU decomposition), which is faster and more accurate:

import numpy as np

A = np.array([[2.0, 1.0],
              [1.0, 3.0]])
b = np.array([5.0, 10.0])

x = np.linalg.solve(A, b)
print(x)                       # [1. 3.]
print(np.allclose(A @ x, b))   # True: always check the residual

@ is Python's matrix multiplication operator. np.linalg.det(A) gives the determinant (here 5).

A three-unknown example

A lab mixes three stock solutions to hit target amounts of three nutrients. Each column of A is one solution's nutrient content per litre, x is the litres of each, and b is the target. One call to np.linalg.solve(A, b) gives the recipe. If the answer contains a negative volume, the target is not reachable with those stocks, and the maths is telling you something real.

Common mistakes

• Mismatched sizes. An m × n matrix times a vector needs n components. NumPy will raise an error, which is helpful.
• Using * instead of @. In NumPy, A * B multiplies element by element; A @ B is matrix multiplication.
• Assuming AB = BA. Matrix multiplication is not commutative.
• Inverting to solve. Use np.linalg.solve, not np.linalg.inv(A) @ b.
• Ignoring a near-zero determinant. A nearly singular system amplifies measurement errors; the answer may be meaningless.

Quick recap

• A system of linear equations is A x = b: coefficients, unknowns, right-hand side.
• Matrix-vector and matrix-matrix products are built from dot products of rows and columns.
• A unique solution exists when det(A) ≠ 0.
• Solve with np.linalg.solve and check that A @ x is close to b.

Try it

In the Do tab, turn a three-equation system into a NumPy matrix and solve it.