Linear Algebra

import numpy as np
import numpy.linalg as la

# 1) Make some matrices
A = np.array([[3., 1.], [1., 2.]])
b = np.array([9., 8.])

# 2) Solve a linear system Ax = b — use solve, not inv()
x = la.solve(A, b)
# array([2.        , 3.        ])

# DON'T do this — slower, less numerically stable
# x = la.inv(A) @ b

# 3) Least squares — overdetermined systems
X = np.random.randn(100, 3)
y = X @ np.array([1.0, 2.0, -1.0]) + 0.1 * np.random.randn(100)
beta, residuals, rank, sv = la.lstsq(X, y, rcond=None)

# 4) Determinant + inverse — small matrices only
print(la.det(A))            # 5.0
print(la.inv(A))            # ok for 2x2, avoid for huge matrices

# 5) Eigenvalues + eigenvectors — symmetric uses eigh (faster, stable)
M = np.array([[4., 1.], [1., 3.]])
w, v = la.eigh(M)
print(w)                    # eigenvalues
print(v)                    # eigenvectors as columns
# Recover: M ≈ v @ np.diag(w) @ v.T

# 6) SVD — the workhorse decomposition
U, S, Vt = la.svd(X, full_matrices=False)
rank = (S > 1e-10).sum()

# Low-rank reconstruction (e.g. for PCA or compression)
k = 2
X_k = (U[:, :k] * S[:k]) @ Vt[:k]

# 7) QR — orthogonal factorisation
Q, R = la.qr(X)

# 8) Cholesky — for positive-definite matrices (e.g. covariance)
C = np.cov(X.T)
L = la.cholesky(C)             # L @ L.T == C

# 9) Norms
la.norm(b)                    # L2 norm — default
la.norm(b, ord=1)             # L1 (Manhattan)
la.norm(b, ord=np.inf)        # max absolute value
la.norm(A, ord='fro')         # Frobenius

# 10) Solving sparse systems — switch libraries
# from scipy.sparse import csr_matrix
# from scipy.sparse.linalg import spsolve
# x = spsolve(csr_matrix(A_big), b)

# 11) Beware: BLAS does the heavy lifting
# np.show_config()  → look for openblas/mkl. On Apple silicon, vecLib is great.
# Set OMP_NUM_THREADS to limit cores in multi-process settings.

import numpy as np
A = np.array([[1, 2], [3, 4]])
b = np.array([5, 6])
x = np.linalg.solve(A, b)  # Ax = b
print(x, np.linalg.det(A))