LINEAR ALGEBRA / 5. LINEAR SYSTEMS

Systems of Linear Equations

Ax = b — solving the core problem of linear algebra

EXPLANATION

A system of linear equations is represented as Ax = b:
• A → coefficient matrix (m×n)
• x → unknown vector (n×1)
• b → right-hand side vector (m×1)

Three possible situations:
1. Unique solution → det(A) ≠ 0, A is invertible. x = A⁻¹b
2. No solution → inconsistent (parallel lines in 2D)
3. Infinite solutions → underdetermined (a line, plane of solutions)

Gaussian Elimination: the standard algorithm.
Row operations that don't change the solution:
• Swap two rows
• Multiply a row by a scalar
• Add a multiple of one row to another

Augmented matrix [A|b] → row reduce to [I|x]

The inverse A⁻¹:
• AA⁻¹ = I
• Exists only when det(A) ≠ 0
• For linear systems: x = A⁻¹b — but computing A⁻¹ is expensive!
• In practice: use np.linalg.solve(A, b) — more stable than A⁻¹ @ b

In ML: normal equation for linear regression:
θ = (XᵀX)⁻¹Xᵀy — the closed-form solution to least squares.

DIAGRAM

System:  2x + y = 5
           x - y  = 1

  Matrix form: [[2, 1],  [x]   [5]
                [1,-1]] ×[y] = [1]

  Augmented matrix → row reduce:
  [2  1 | 5]       [1  0 | 2]
  [1 -1 | 1]  →→→  [0  1 | 1]

  Solution: x=2, y=1

  Geometric interpretation (2D):
  Each equation = a line
  Solution = intersection of lines

  2x+y=5:  ─────╲
  x-y=1:   ─────╱
  Intersection: (2,1)  ← unique solution

CODE

PYTHON
1import numpy as np
2from scipy import linalg
3
4# ── Solve Ax = b ──────────────────────────────────────────────────
5A = np.array([[2., 1.],
6              [1., -1.]])
7b = np.array([5., 1.])
8
9# Method 1: np.linalg.solve (preferred — numerically stable)
10x = np.linalg.solve(A, b)
11print(f"Solution x = {x}")
12print(f"Verify Ax = {A @ x} (should be {b})")
13
14# Method 2: inverse (slower, less stable — avoid for large systems)
15x_inv = np.linalg.inv(A) @ b
16print(f"Via inverse: {x_inv}")
17
18# ── Gaussian elimination (manual demonstration) ───────────────────
19def gaussian_elimination(A, b):
20    n = len(b)
21    # Augmented matrix
22    M = np.hstack([A.astype(float), b.reshape(-1,1)])
23    print("
24Augmented matrix:")
25    print(M)
26
27    for col in range(n):
28        # Find pivot
29        pivot = np.argmax(abs(M[col:, col])) + col
30        M[[col, pivot]] = M[[pivot, col]]   # swap rows
31
32        # Eliminate below
33        for row in range(col+1, n):
34            factor = M[row, col] / M[col, col]
35            M[row] -= factor * M[col]
36
37    print("After forward elimination:")
38    print(M)
39
40    # Back substitution
41    x = np.zeros(n)
42    for i in range(n-1, -1, -1):
43        x[i] = (M[i, -1] - M[i, i+1:n] @ x[i+1:n]) / M[i, i]
44    return x
45
46solution = gaussian_elimination(A, b)
47print(f"
48Gaussian elimination solution: {solution}")
49
50# ── Normal equation: linear regression closed form ────────────────
51# θ = (XᵀX)⁻¹Xᵀy
52np.random.seed(42)
53n, d = 100, 3
54X = np.hstack([np.ones((n,1)), np.random.randn(n, d)])  # add bias
55true_theta = np.array([1.0, 2.0, -1.0, 0.5])
56y = X @ true_theta + 0.1 * np.random.randn(n)
57
58# Normal equation
59theta_normal = np.linalg.solve(X.T @ X, X.T @ y)
60print(f"
61Linear regression (normal equation):")
62print(f"  True theta:      {true_theta}")
63print(f"  Estimated theta: {theta_normal.round(3)}")
64
65# sklearn comparison
66from sklearn.linear_model import LinearRegression
67lr = LinearRegression(fit_intercept=False).fit(X, y)
68print(f"  sklearn theta:   {lr.coef_.round(3)}")
← PREV4. DeterminantsNEXT →6. Eigenvalues & Eigenvectors