Linear Algebra with NumPy: dot, matmul, solve, eig & norms

Linear algebra forms the mathematical foundation of machine learning, physics simulations, and scientific computing. Utilizing optimized C-level BLAS and LAPACK libraries under the hood, NumPy executes matrix operations significantly faster than pure Python implementations. This comprehensive guide covers matrix dot products, the @ operator, stable linear system solvers, eigenvalue decompositions, least-squares fitting, and numerical conditioning. For related element-wise mathematical ufuncs, review Math & ufuncs and Broadcasting.

Overview of Linear Algebra in NumPy 🔝

Operations in numpy.linalg are designed to mirror standard mathematical notation while maintaining high performance. Arrays are treated as linear transformations acting on vector spaces. Ensuring proper shape alignment (e.g., matching inner dimensions for matrix multiplication) is essential for avoiding shape mismatch errors.

Code Block Modernization Comparison 🔝

The table below compares legacy NumPy linear algebra snippets with the updated, robust, and self-verifying implementation:

1. Matrix Multiplication: dot vs @ (matmul) 🔝

np.dot and the @ operator perform matrix multiplication for 2D arrays. For 1D arrays, dot computes the inner product. The @ operator (PEP 465) is preferred for readability and supports batched matrix multiplication across higher dimensions.

import numpy as np

A = np.array([
    [1, 2],
    [3, 4]
])    # shape (2,2)
B = np.array([
    [5, 6],
    [7, 8]
])    # shape (2,2)
x = np.array([10, 20])     # shape (2,)

# Matrix @ Matrix
print("Matrix product A @ B:
", A @ B)
# Matrix @ Vector
print("Matrix-vector product A @ x:", A @ x)

# Dot product for 1D vectors (inner product)
u = np.array([1, 2, 3])
v = np.array([4, 5, 6])
print("Inner product np.dot(u, v):", np.dot(u, v))

assert np.array_equal(A @ x, np.array([50, 110]))
assert np.dot(u, v) == 32

Expected output

Matrix product A @ B:
 [[19 22]
 [43 50]]
Matrix-vector product A @ x: [ 50 110]
Inner product np.dot(u, v): 32

2. Solving Linear Systems: np.linalg.solve 🔝

Solve linear matrix equations of the form Ax = b using np.linalg.solve(). This approach is numerically stable and computationally faster than computing the matrix inverse directly.

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("Solution vector x:", x)

# Verify identity Ax = b
print("Verification A @ x == b:", np.allclose(A @ x, b))

assert np.allclose(x, np.array([2.0, 3.0]))

Expected output

Solution vector x: [2. 3.]
Verification A @ x == b: True

3. Eigenvalues & Eigenvectors: eig and eigh 🔝

np.linalg.eig() computes eigenvalues and eigenvectors for general square matrices. For real symmetric or Hermitian matrices, np.linalg.eigh() is specialized, offering faster execution and superior numerical stability.

S = np.array([
    [2.0, 1.0],
    [1.0, 2.0]
])
w, v = np.linalg.eigh(S)

print("Eigenvalues (w):", w)
print("Eigenvectors (v):
", v)

# Verify eigen-decomposition identity Sv = vw
print("Eigen decomposition verified:", np.allclose(S @ v, v * w))
assert np.allclose(S @ v, v * w)

Expected output

Eigenvalues (w): [1. 3.]
Eigenvectors (v):
 [[-0.70710678  0.70710678]
 [ 0.70710678  0.70710678]]
Eigen decomposition verified: True

4. Norms & Numerical Conditioning 🔝

Vector and matrix norms quantify magnitude and scale, while condition numbers estimate how sensitive a matrix inversion or linear system solver is to rounding errors.

M = np.array([
    [1.0, 2.0],
    [3.0, 4.0]
])
x = np.array([3.0, 4.0])

print("L2 Vector Norm:", np.linalg.norm(x))
print("Spectral Matrix Norm:", np.linalg.norm(M, ord=2))
print("Condition Number:", np.linalg.cond(M))

assert np.linalg.norm(x) == 5.0

Expected output

L2 Vector Norm: 5.0
Spectral Matrix Norm: 5.464985704219043
Condition Number: 14.933034373659688

5. Least Squares & Rank-Deficient Cases 🔝

When solving overdetermined linear systems where equations exceed unknowns, np.linalg.lstsq() computes the least-squares best-fit solution minimizing squared residuals:

A = np.array([
    [1.0, 1.0],
    [1.0, 2.0],
    [1.0, 3.0]
])  # (3,2) design matrix
b = np.array([1.0, 2.0, 2.5])

x_ls, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None)
print("Least squares coefficients (intercept, slope):", x_ls)

assert x_ls.shape == (2,)

Expected output

Least squares coefficients (intercept, slope): [0.66666667 0.75      ]

6. Determinant, Inverse & Pseudo-Inverse 🔝

Determinants, inverses, and Moore-Penrose pseudo-inverses are supported, though solvers should be preferred over direct matrix inversion for numerical stability:

A = np.array([
    [2.0, 1.0],
    [5.0, 3.0]
])

print("Determinant:     ", np.linalg.det(A))
print("Matrix Inverse:
", np.linalg.inv(A))
print("Pseudo-Inverse:
", np.linalg.pinv(A))

assert round(np.linalg.det(A), 4) == 1.0

Expected output

Determinant:      1.0
Matrix Inverse:
 [[ 3. -1.]
 [-5.  2.]]
Pseudo-Inverse:
 [[ 3. -1.]
 [-5.  2.]]

6. Advanced Matrix Decompositions: SVD & QR 🔝

Beyond eigenvalues, np.linalg supports core matrix factorizations utilized in dimensionality reduction (PCA) and least-squares regression:

import numpy as np

M = np.array([
    [1.0, 2.0],
    [3.0, 4.0],
    [5.0, 6.0]
])

# 1. Singular Value Decomposition (SVD)
U, s, Vt = np.linalg.svd(M)
print("Singular values (s):", s)
assert len(s) == 2

# 2. QR Decomposition (A = QR)
Q, R = np.linalg.qr(M)
print("Q orthogonal matrix shape:", Q.shape)
print("R upper triangular shape:", R.shape)

assert np.allclose(Q @ R, M)

Expected output

Singular values (s): [9.52551809 0.51430058]
Q orthogonal matrix shape: (3, 2)
R upper triangular shape: (2, 2)

7. Comprehensive Troubleshooting & Numerical Stability Matrix 🔝

Common Pitfalls & Troubleshooting 🔝

• Direct Inversion Instability: Using np.linalg.inv(A) @ b accumulates rounding errors. Always prefer np.linalg.solve(A, b).
• Shape Mismatches in Matmul: Inner dimensions must match (e.g., multiplying (2, 3) by (3, 4) yields (2, 4)).
• Ill-Conditioned Systems: High condition numbers amplify tiny data perturbations into massive output errors. Use regularization or least-squares.

Frequently Asked Questions 🔝

Practice in Google Colab 🔝

Test matrix multiplication, linear system solvers, eigenvalue decomposition, least-squares regression, and run verification test suites in our interactive Google Colab companion notebook:

Interactive Practice: Colab & GitHub

Launch the notebook directly in Google Colab to run live code cells or inspect the source repository on GitHub:

View on GitHub

Continue exploring related NumPy features: Math & ufuncs, Vectorization, File I/O, Aggregations, and Views & Strides.

Subhendu Mohapatra

Author & Instructor at plus2net

I write and maintain practical tutorials on Python, PHP, SQL, JavaScript, HTML, jQuery, and web development at plus2net. The tutorials focus on clear explanations, working examples, and code that readers can test and adapt while learning.

About the author   Watch on YouTube

Connect: