Solving linear systems in Python

solve is the right answer; inverting the matrix is not.

import numpy as np

A = np.array([[3.0, 1.0], [1.0, 2.0]])
b = np.array([9.0, 8.0])

x = np.linalg.solve(A, b)
print(x, np.allclose(A @ x, b))

print(np.linalg.det(A).round(3), np.linalg.matrix_rank(A))
print(np.linalg.norm(x).round(3))

tall = np.array([[1.0, 1.0], [1.0, 2.0], [1.0, 3.0]])
fit, *_ = np.linalg.lstsq(tall, np.array([2.0, 4.0, 6.1]), rcond=None)
print(fit.round(3))

How it works

solve is faster and better conditioned than an explicit inverse.
allclose is how you compare floating point results.
lstsq handles the over-determined case.

Keywords and builtins used here

as
print

The run, in numbers

Lines: 14
Characters to type: 382
Tokens: 169
Three-star pace: 110 tpm

At the three-star pace of 110 tokens a minute, this run takes about 92 seconds.

Type this snippet

Step 3 of 4 in Maths & memory, step 18 of 22 in NumPy.

← Previous Next →