Cramer’s rule vs. elimination

Cramer’s rule gives every unknown as a ratio of two determinants, x_i = det(A_i) / det(A). It is exact in theory but needs n + 1 determinants, each an O(n^3) factorization, so its cost grows like n^4 against elimination’s n^3.

import time

import matplotlib.pyplot as plt
import numpy as np

from mathematicskit.linalg import cramer_solve, lu_solve_system

A 3x3 system solved by determinants

A = np.array([[2.0, 1.0, -1.0], [-3.0, -1.0, 2.0], [-2.0, 1.0, 2.0]])
b = np.array([8.0, -11.0, -3.0])
print("det(A) =", round(np.linalg.det(A), 10))
for i in range(3):
    A_i = A.copy()
    A_i[:, i] = b
    print(f"det(A_{i}) = {np.linalg.det(A_i):+.10f}")
print("Cramer solution:", np.round(cramer_solve(A, b), 10))

det(A) = -1.0
det(A_0) = -2.0000000000
det(A_1) = -3.0000000000
det(A_2) = +1.0000000000
Cramer solution: [ 2.  3. -1.]

Cost as n grows

rng = np.random.default_rng(0)
sizes = [10, 20, 40, 80, 160]
t_cramer, t_lu = [], []
for n in sizes:
    M = rng.normal(size=(n, n))
    rhs = rng.normal(size=n)
    t0 = time.perf_counter()
    x_c = cramer_solve(M, rhs)
    t_cramer.append(time.perf_counter() - t0)
    t0 = time.perf_counter()
    x_lu = lu_solve_system(M, rhs)
    t_lu.append(time.perf_counter() - t0)
    print(f"n = {n:3d}: max |x_cramer - x_lu| = {np.max(np.abs(x_c - x_lu)):.1e}")

fig, ax = plt.subplots(figsize=(6, 4))
ax.loglog(sizes, t_cramer, "o-", label="Cramer's rule")
ax.loglog(sizes, t_lu, "s-", label="LU solve")
ax.set_xlabel("n")
ax.set_ylabel("time (s)")
ax.legend()
fig.tight_layout()

plt.show()

n =  10: max |x_cramer - x_lu| = 1.6e-14
n =  20: max |x_cramer - x_lu| = 1.6e-14
n =  40: max |x_cramer - x_lu| = 2.7e-14
n =  80: max |x_cramer - x_lu| = 4.5e-13
n = 160: max |x_cramer - x_lu| = 2.5e-13