Solving Linear Equations for Data Science

Introduction

A linear equation is an equation where every variable appears only to the first power and is not multiplied by another variable. For example, 2x + 3y = 12 is linear, while x² + y = 4 and xy = 6 are not.

A system of linear equations is a set of linear equations that must all be true at the same time. Solving a system means finding values for the unknowns that satisfy every equation.

2x + 3y = 12
 x −  y = 1

Systems of linear equations appear everywhere in machine learning. Linear regression is the problem of finding the best-fitting solution to a (usually inconsistent) system of equations, and many optimisation and inference algorithms solve such systems internally.

Methods for solving a system

1. Substitution

Solve one equation for one variable, then substitute into the other.

x − y = 1  →  x = y + 1
2(y + 1) + 3y = 12  →  5y = 10  →  y = 2,  x = 3

2. Elimination

Add or subtract multiples of the equations to cancel a variable.

3. Graphical method

Each equation is a line in the plane. The solution is where the lines intersect.

Three possible outcomes

Case	Description	Geometry (2 variables)
Unique solution	Exactly one solution	Lines cross at one point
No solution (inconsistent)	Equations contradict each other	Parallel lines
Infinitely many solutions	Equations are redundant	The lines coincide

Examples:

x + y = 2 and x + y = 5: no solution (contradiction).
x + y = 2 and 2x + 2y = 4: infinitely many solutions (the second is just twice the first).

Matrix form: Ax = b

Any system of linear equations can be written compactly as

A · x = b

A is the coefficient matrix,
x is the vector of unknowns,
b is the vector of constants (right-hand sides).

2x + 3y = 12         ⎡ 2  3 ⎤ ⎡ x ⎤   ⎡ 12 ⎤
 x −  y = 1     →    ⎣ 1 −1 ⎦ ⎣ y ⎦ = ⎣  1 ⎦

The augmented matrix [A | b] places the constants next to the coefficients and is convenient for row operations.

Gaussian elimination

Gaussian elimination solves a system by applying row operations to turn the augmented matrix into a triangular form, then using back-substitution. The three allowed row operations are:

Swap two rows.
Multiply a row by a non-zero number.
Add a multiple of one row to another row.

Example. Solve

 x +  y +  z =  6
      2y + 5z = −4
2x + 5y −  z = 27

Eliminate x from the third equation: R₃ ← R₃ − 2R₁ gives 3y − 3z = 15.
Eliminate y from the new third equation: R₃ ← R₃ − 1.5·R₂ gives −3z − 7.5z = 15 + 6, so −10.5z = 21 and z = −2.
Back-substitute: 2y + 5(−2) = −4 gives y = 3. Then x + 3 − 2 = 6 gives x = 5.

The solution is (x, y, z) = (5, 3, −2).

Solving with the inverse

If A is square and invertible, multiply both sides by A⁻¹:

x = A⁻¹ b

Example. [[2, 1], [5, 3]] · x = [4, 11]. The determinant is 6 − 5 = 1 and A⁻¹ = [[3, −1], [−5, 2]]. So x = [3×4 − 11, −5×4 + 2×11] = [1, 2].

In practice, software uses factorisation-based solvers (LU, QR) rather than explicitly computing A⁻¹, because they are faster and more numerically stable.

Cramer's rule (2 × 2)

For a₁x + b₁y = c₁ and a₂x + b₂y = c₂, with D = a₁b₂ − a₂b₁ ≠ 0:

x = (c₁b₂ − c₂b₁) / D ,    y = (a₁c₂ − a₂c₁) / D

Cramer's rule is elegant for small systems but far too slow for large ones.

When does a unique solution exist?

For a square system Ax = b, all of these are equivalent:

A is invertible.
det A ≠ 0.
A has full rank (its columns are linearly independent).
Ax = 0 has only the solution x = 0.

If det A = 0 the system has either no solution or infinitely many.

Overdetermined and underdetermined systems

Type	Equations vs unknowns	Typical situation
Square	equal	Exact solution possible
Overdetermined	more equations than unknowns	Real data: many samples, few parameters. Usually no exact solution
Underdetermined	fewer equations than unknowns	Infinitely many solutions (common in high-dimensional problems)

Least squares and the normal equation

For an overdetermined system Xw = y (many samples, few weights), an exact solution usually does not exist because of noise. Instead we find the w that minimises the total squared error ‖Xw − y‖². Setting the gradient to zero gives the normal equation:

XᵀX w = Xᵀy      →      w = (XᵀX)⁻¹ Xᵀy

This is the closed-form solution of linear regression.

Ridge regression adds a penalty: w = (XᵀX + λI)⁻¹ Xᵀy. The extra λI makes the matrix invertible and stabilises the solution when features are correlated.

Conditioning

If two features are nearly identical, XᵀX is nearly singular: tiny changes in the data cause huge changes in w. Such a system is ill-conditioned. Remedies include removing redundant features, regularisation, and using gradient-based optimisation.

Real-life applications

Engineering: solving for currents in electrical circuits and forces in structures.
Economics: equilibrium prices and input-output models.
Chemistry: balancing chemical equations.
Traffic and networks: flow conservation equations.
Computer graphics and robotics: solving for positions and joint angles.

Where it appears in machine learning

Linear regression (normal equation and least squares).
Ridge regression and other regularised solutions.
Gaussian processes, Kalman filters and other models solve linear systems repeatedly.
Second-order optimisation (Newton's method) solves a linear system at every step.
PageRank is defined as the solution to a system of linear equations.

Python implementation

import numpy as np

# 2x + 3y = 12 and x - y = 1
A = np.array([[2, 3],
              [1, -1]])
b = np.array([12, 1])

x = np.linalg.solve(A, b)
print(x)                 # [3. 2.]

# The 3x3 system from the Gaussian elimination example
A3 = np.array([[1, 1, 1],
               [0, 2, 5],
               [2, 5, -1]])
b3 = np.array([6, -4, 27])
print(np.linalg.solve(A3, b3))      # [ 5.  3. -2.]

# Check singularity before solving
M = np.array([[1, 2], [2, 4]])
print(np.linalg.det(M))             # 0.0  -> no unique solution

Least squares (linear regression)

import numpy as np

# Points (1,1), (2,2), (3,2): find the best line y = w0 + w1 * x
x = np.array([1, 2, 3])
y = np.array([1, 2, 2])

X = np.column_stack([np.ones_like(x), x])          # design matrix with intercept column

# Method 1: normal equation
w = np.linalg.inv(X.T @ X) @ X.T @ y
print(w)                                            # [0.6667 0.5]

# Method 2: numerically stable least squares (preferred)
w2, *_ = np.linalg.lstsq(X, y, rcond=None)
print(w2)                                           # [0.6667 0.5]

Solved examples

Example 1: Elimination

Solve 2x + 3y = 12 and x − y = 1.

Solution. From the second equation, x = y + 1. Substitute: 2(y + 1) + 3y = 12, so 5y = 10, y = 2, and x = 3.

Check: 2(3) + 3(2) = 12 ✓ and 3 − 2 = 1 ✓.

Example 2: Using the inverse

Solve [[2, 1], [5, 3]] · [x, y] = [4, 11] with the inverse.

Solution. det = 2×3 − 1×5 = 1, so A⁻¹ = [[3, −1], [−5, 2]].

x = A⁻¹b = [3×4 − 1×11, −5×4 + 2×11] = [1, 2].

Check: 2(1) + 1(2) = 4 ✓ and 5(1) + 3(2) = 11 ✓.

Example 3: Gaussian elimination

Solve x + y + z = 6, 2y + 5z = −4, 2x + 5y − z = 27.

Solution. As worked above: R₃ − 2R₁ gives 3y − 3z = 15; eliminating y gives z = −2; back-substitution gives y = 3 and x = 5. The answer is (5, 3, −2).

Check: 5 + 3 − 2 = 6 ✓; 2(3) + 5(−2) = −4 ✓; 2(5) + 5(3) − (−2) = 27 ✓.

Example 4: Classifying systems

Classify (a) x + y = 2, x + y = 5 and (b) x + y = 2, 2x + 2y = 4.

Solution. (a) The left sides are identical but the right sides differ, so it is contradictory: no solution. (b) The second equation is 2 × the first, so they describe the same line: infinitely many solutions.

Example 5: Cramer's rule

Solve 3x + 2y = 7 and x − y = 1.

Solution. D = 3(−1) − 2(1) = −5.

x = (7(−1) − 1·2) / D = −9/−5 = 1.8. y = (3·1 − 1·7) / D = −4/−5 = 0.8.

Check: 3(1.8) + 2(0.8) = 7 ✓ and 1.8 − 0.8 = 1 ✓.

Example 6: Least squares

Find the best-fit line y = w₀ + w₁x through the points (1, 1), (2, 2), (3, 2).

Solution. The design matrix is X = [[1, 1], [1, 2], [1, 3]] and y = [1, 2, 2].

XᵀX = [[3, 6], [6, 14]] and Xᵀy = [5, 11].

det(XᵀX) = 42 − 36 = 6, so (XᵀX)⁻¹ = (1/6)[[14, −6], [−6, 3]].

w = (1/6)[14×5 − 6×11, −6×5 + 3×11] = (1/6)[4, 3] = [0.667, 0.5].

The best-fit line is y = 0.667 + 0.5x.

Common mistakes

Trying to solve a system when det A = 0 as if it had a unique solution.
Computing A⁻¹ explicitly for large systems instead of using solve or lstsq.
Forgetting the column of ones for the intercept in the design matrix.
Ignoring ill-conditioning when features are nearly duplicates.

Practice problems

Solve x + y = 10 and x − y = 2.
Solve 3x + y = 9 and x + 2y = 8.
Write 2x − y = 3 and x + 4y = 6 in matrix form Ax = b, then solve.
Classify 2x + 4y = 6 and x + 2y = 3.
Solve x + 2y = 5 and 3x − y = 1 with Cramer's rule.
Find the least-squares line through (0, 1), (1, 3), (2, 5).

Answers

Adding gives 2x = 12, so x = 6, y = 4
From the first, y = 9 − 3x. Then x + 18 − 6x = 8, so x = 2, y = 3
[[2, −1], [1, 4]] · [x, y] = [3, 6]; y = 2x − 3 gives x + 8x − 12 = 6, so x = 2, y = 1
The second equation times 2 gives the first, so infinitely many solutions
D = −7; x = (5(−1) − 2(1))/−7 = 1; y = (1·1 − 5·3)/−7 = 2, so (1, 2)
The points lie exactly on a line: y = 1 + 2x

Key takeaways

A system of linear equations can be written Ax = b and has zero, one or infinitely many solutions.
Gaussian elimination is the workhorse; use np.linalg.solve in practice.
A unique solution exists exactly when A is square and invertible (det A ≠ 0).
With more equations than unknowns (real data), we use least squares and the normal equation, which is exactly linear regression.
Correlated features make the system ill-conditioned; regularisation helps.

Linear Equations