Think of NumPy linear algebra as a calculator that works on whole tables of numbers at once. Instead of one value times another, you multiply a grid of inputs by a grid of weights in a single step. This numpy linear algebra guide covers the operations you will actually use in data science: dot products, matrix multiplication, the determinant, the inverse, solving linear systems, and eigenvalues. Plain English first, then the NumPy code, with every output run on real Python.

“Linear algebra is the mathematics of data.”

Gilbert Strang, MIT

Last Updated: July 2026 | Tested on: Python 3.14.6, NumPy 2.4.6 | Difficulty: Advanced | Reading Time: 16 minutes

Here is the plain-English version before any math. A vector is just a list of numbers, like one person’s height, weight, and age. A matrix is a stack of those lists, like a spreadsheet where every row is a person and every column is a feature. Linear algebra is the set of rules for combining these grids quickly. That is the whole idea. Everything below is detail.

Why bother? Because every machine learning model runs on this. When scikit-learn calls model.fit(X, y), under the hood it is solving a system of linear equations. When a neural network makes a prediction, it is multiplying matrices. So when you learn the NumPy linalg module here, you are not learning a side topic. You are learning the engine that .fit() and .predict() sit on top of.

Picture a coffee bill to make matrix multiplication click. You order 2 lattes and 3 espressos. A latte is 4 rupees, an espresso is 2. Your bill is 2 times 4 plus 3 times 2, which is 14. That single number, quantities times prices added up, is a dot product. A matrix multiply is just doing that bill for many customers and many price lists at once, in one step. Hold that picture. It is all of linear algebra in one sentence.

Prathamesh, a 24-year-old on our team, once spent an afternoon chasing a linear regression that spat out absurd predictions. The culprit was a singular matrix: two of his features were perfectly correlated, so the matrix could not be inverted. Knowing a little linear algebra told him exactly what had gone wrong, instead of blindly retrying .fit() and hoping. You learned vectors and matrices conceptually in the math for data science tutorial; this post makes them run fast with NumPy.

The diagram sorts NumPy’s linear algebra tools into three buckets: basic operations (dot product, matrix multiplication, transpose), decomposition (eigenvalues, Singular Value Decomposition or SVD, QR), and system solving (inverse, linear systems, least squares). Each one maps to a function in np.linalg. Treat it like a menu. Need to solve a set of equations? Look under “Solving Systems.” Need to find the directions your data spreads in? Look under “Decomposition.”

Prerequisites

Complete NumPy broadcasting tutorial and math for data science tutorial (vectors and matrices section).

Dot Product & Matrix Multiplication

Start with the coffee bill from earlier. The dot product takes two lists of numbers, multiplies them position by position, and adds it all up into one number. Quantities times prices, summed. Matrix multiplication is the same move repeated for many rows and many columns at once. In NumPy you reach for the @ operator, which reads like “times” for matrices.

📄 dot_matmul.py: the two operations you will use most

import numpy as np

# Dot product of two vectors
a = np.array([1, 2, 3])
b = np.array([4, 5, 6])
dot = np.dot(a, b)  # 1*4 + 2*5 + 3*6
print(f"Dot product: {a} · {b} = {dot}")
print(f"Using @:     {a @ b}")

# Matrix multiplication
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])

print(f"\nA @ B =\n{A @ B}")
print(f"np.matmul(A, B) =\n{np.matmul(A, B)}")

# Shape rule: (m, n) @ (n, p) -> (m, p)
X = np.random.default_rng(42).random((3, 4))  # 3 samples, 4 features
W = np.random.default_rng(42).random((4, 2))  # 4 features, 2 outputs
result = X @ W  # (3, 4) @ (4, 2) -> (3, 2)
print(f"\nNeural network forward pass:")
print(f"Input (3,4) @ Weights (4,2) = Output {result.shape}")
print(f"Output:\n{result.round(3)}")

▶ Output

Dot product: [1 2 3] · [4 5 6] = 32
Using @:     32

A @ B =
[[19 22]
 [43 50]]
np.matmul(A, B) =
[[19 22]
 [43 50]]

Neural network forward pass:
Input (3,4) @ Weights (4,2) = Output (3, 2)
Output:
[[1.587 2.032]
 [1.581 2.082]
 [1.226 1.461]]

What happened here: The dot product of [1 2 3] and [4 5 6] is 1 times 4 plus 2 times 5 plus 3 times 6, which is 32. Both np.dot(a, b) and a @ b give the same answer for vectors. For the two 2×2 grids, A @ B and np.matmul(A, B) produce the identical result.

The last block is the move every neural network makes: a batch of 3 samples, each with 4 features, multiplied by a weight grid of shape (4, 2), turns into 3 outputs of size 2. The one rule to remember is that the inner numbers must match: (m, n) @ (n, p) works and gives you (m, p). The two 4s touch in the middle and cancel out, leaving (3, 2).

Transpose, Determinant & Inverse

Three quick ideas before the code. Transpose flips a matrix on its side: rows become columns. The determinant is a single number that tells you whether a square matrix can be “undone”. If it is zero, the matrix is squashed flat and there is no way back. The inverse is that undo button: multiply a matrix by its inverse and you land back on the identity matrix, the linear algebra version of multiplying a number by its reciprocal to get 1. Think of dividing by 4 by multiplying by one quarter. The inverse is the one quarter.

📄 matrix_ops.py: transpose, determinant, and the undo button

import numpy as np

A = np.array([[1, 2, 3],
              [4, 5, 6]])

# Transpose: swap rows and columns
print(f"A (2x3):\n{A}")
print(f"A.T (3x2):\n{A.T}\n")

# Determinant: scalar measure of a square matrix
B = np.array([[3, 1],
              [2, 4]])
det = np.linalg.det(B)
print(f"B:\n{B}")
print(f"det(B) = {det:.1f}")
print(f"det != 0 means B is invertible: {det != 0}\n")

# Inverse: A @ A_inv = Identity
B_inv = np.linalg.inv(B)
print(f"B inverse:\n{B_inv.round(3)}")
print(f"B @ B_inv = Identity:\n{(B @ B_inv).round(1)}\n")

# Singular matrix (not invertible)
singular = np.array([[1, 2], [2, 4]])  # Row 2 = 2 * Row 1
print(f"Singular det: {np.linalg.det(singular):.1e}")
try:
    np.linalg.inv(singular)
except np.linalg.LinAlgError as e:
    print(f"Cannot invert: {e}")

▶ Output

A (2x3):
[[1 2 3]
 [4 5 6]]
A.T (3x2):
[[1 4]
 [2 5]
 [3 6]]

B:
[[3 1]
 [2 4]]
det(B) = 10.0
det != 0 means B is invertible: True

B inverse:
[[ 0.4 -0.1]
 [-0.2  0.3]]
B @ B_inv = Identity:
[[1. 0.]
 [0. 1.]]

Singular det: 0.0e+00
Cannot invert: Singular matrix

What happened here: For a 2×2 matrix the determinant is easy to check by hand. For B = [[3, 1], [2, 4]] it is 3 times 4 minus 1 times 2, which is 12 minus 2, so 10. NumPy agrees: det(B) = 10.0. Because that is not zero, B has an inverse, and multiplying B @ B_inv gives back the identity matrix [[1, 0], [0, 1]]. The last block is the trap Prathamesh hit.

The matrix [[1, 2], [2, 4]] has its second row as exactly twice the first, so its determinant is 1 times 4 minus 2 times 2, which is zero. A zero determinant means there is no inverse, and NumPy raises LinAlgError: Singular matrix instead of returning garbage. That error is your friend: it is telling you two of your features carry the same information.

Solving Linear Systems

Remember those school word problems? “Two coffees and one tea cost 5 rupees; one coffee and three teas cost 7. Find the price of each.” That is a linear system, and np.linalg.solve(A, b) cracks it in one line. You hand it the coefficients as a matrix A and the answers as a vector b, and it finds the unknowns. This is the exact machinery behind fitting a line to data, which is why linear regression is really just solving a system.

📄 solve_system.py: solve Ax = b, the foundation of linear regression

import numpy as np

# System of equations:
# 2x + y = 5
# x + 3y = 7
A = np.array([[2, 1],
              [1, 3]])
b = np.array([5, 7])

# Solve for x
x = np.linalg.solve(A, b)
print(f"Solution: x = {x[0]:.2f}, y = {x[1]:.2f}")
print(f"Verify: A @ x = {A @ x}")

# Real-world: Find line of best fit (y = mx + b)
# Points: (1, 3), (2, 5), (3, 7), (4, 9)
points_x = np.array([1, 2, 3, 4])
points_y = np.array([3, 5, 7, 9])

# Set up least squares: A @ [m, b] = y
A_ls = np.column_stack([points_x, np.ones(4)])
result = np.linalg.lstsq(A_ls, points_y, rcond=None)
m, b_val = result[0]
print(f"\nBest fit line: y = {m:.1f}x + {b_val:.1f}")
print(f"Predictions: {(A_ls @ result[0]).round(1)}")
print(f"Actual:      {points_y}")

▶ Output

Solution: x = 1.60, y = 1.80
Verify: A @ x = [5. 7.]

Best fit line: y = 2.0x + 1.0
Predictions: [3. 5. 7. 9.]
Actual:      [3 5 7 9]

What happened here: The first system, 2x + y = 5 and x + 3y = 7, has the answer x = 1.6, y = 1.8. You can sanity-check it: 2 times 1.6 plus 1.8 is 5, and 1.6 plus 3 times 1.8 is 7. NumPy confirms it by computing A @ x = [5. 7.], which matches the right-hand side exactly. The second block fits a straight line through four points using lstsq (least squares).

The points sit perfectly on y = 2x + 1, so the recovered slope is 2.0 and intercept is 1.0, and the predictions land right on the actual values. Notice we used solve for the exact square system and lstsq for the line fit. That choice matters, and the FAQ at the end explains why.

Eigenvalues & Eigenvectors

Eigenvalues and eigenvectors sound scary, so here is the everyday picture. Imagine spinning a globe. Most points on the surface move, but the two poles stay put on their axis: that axis is special. An eigenvector is that special direction a matrix leaves pointing the same way; it only stretches or shrinks it, and the amount of stretch is the eigenvalue. In data science, the eigenvectors of a covariance matrix point along the directions your data spreads out the most, and the biggest eigenvalue marks the direction with the most spread. That is exactly the engine behind PCA (Principal Component Analysis), which you will meet in the PCA tutorial.

📄 eigen.py: eigendecomposition for data science

import numpy as np

# Covariance matrix of some 2D data
cov_matrix = np.array([[4.0, 1.2],
                         [1.2, 2.0]])

eigenvalues, eigenvectors = np.linalg.eig(cov_matrix)
# np.linalg.eig returns a complex dtype for real matrices; the imaginary parts
# are zero here, so take the real part to work with plain floats.
eigenvalues = eigenvalues.real
eigenvectors = eigenvectors.real
print(f"Eigenvalues: {eigenvalues.round(3)}")
print(f"Eigenvectors:\n{eigenvectors.round(3)}")

# The largest eigenvalue = direction of most variance
max_idx = eigenvalues.argmax()
print(f"\nPrimary direction (PC1): {eigenvectors[:, max_idx].round(3)}")
print(f"Variance explained: {eigenvalues[max_idx]:.1f} out of {eigenvalues.sum():.1f}")
print(f"Ratio: {eigenvalues[max_idx] / eigenvalues.sum():.1%}")

# Verify: A @ v = lambda * v
v = eigenvectors[:, 0]
lam = eigenvalues[0]
print(f"\nA @ v = {(cov_matrix @ v).round(3)}")
print(f"lambda * v = {(lam * v).round(3)}")
print(f"Equal: {np.allclose(cov_matrix @ v, lam * v)}")

▶ Output

Eigenvalues: [4.562 1.438]
Eigenvectors:
[[ 0.906 -0.424]
 [ 0.424  0.906]]

Primary direction (PC1): [0.906 0.424]
Variance explained: 4.6 out of 6.0
Ratio: 76.0%

A @ v = [4.131 1.935]
lambda * v = [4.131 1.935]
Equal: True

What happened here: The covariance matrix [[4.0, 1.2], [1.2, 2.0]] gives two eigenvalues, 4.562 and 1.438. You can check those by hand for a 2×2: they add up to the diagonal sum (4 plus 2 is 6, and 4.562 plus 1.438 is 6) and multiply to the determinant (4 times 2 minus 1.2 times 1.2 is 6.56, and 4.562 times 1.438 is 6.56). The bigger eigenvalue, 4.562, owns 76.0% of the total spread, so its eigenvector [0.906, 0.424] is the main direction the data stretches, which is what PCA calls the first principal component.

The final check proves the definition: A @ v equals 4.131, 1.935, and so does lambda * v. Multiplying the matrix by its eigenvector is the same as just scaling that eigenvector by the eigenvalue. One small note: NumPy is free to flip the sign of an eigenvector (point it the opposite way), so if you see [-0.906, -0.424] on your machine, that is the same direction, just facing the other way.

Common Mistakes

📄 Mistake 1: using np.dot for matrix multiplication of 2D arrays

import numpy as np

# np.dot works but @ is clearer and recommended
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])

# Both work, but @ is the modern standard
print(np.dot(A, B))  # Works but old style
print(A @ B)          # Preferred, clearer intent

▶ Output

[[19 22]
 [43 50]]
[[19 22]
 [43 50]]

Both lines print the same grid, so why prefer @? Because np.dot quietly does different things depending on the shapes you feed it (vectors, 2D matrices, or higher dimensions), while @ says “matrix multiply” and nothing else. Clear intent beats clever overloading every time someone reads your code six months later.

📄 Mistake 2: inverting when you should solve

import numpy as np
A = np.array([[3, 1], [2, 4]])
b = np.array([5, 7])

# BAD: computing the inverse is slower and less numerically stable
x_bad = np.linalg.inv(A) @ b

# GOOD: solve directly (faster and more accurate)
x_good = np.linalg.solve(A, b)

print("x_bad: ", x_bad)
print("x_good:", x_good)

▶ Output

x_bad:  [1.3 1.1]
x_good: [1.3 1.1]

Same answer here, but on bigger or messier matrices the two part ways. Computing the full inverse and then multiplying piles up rounding error and extra work. np.linalg.solve goes straight to the answer with a more stable algorithm. The rule of thumb: if all you want is x in Ax = b, never build the inverse just to throw it away. Reach for solve.

Try It Yourself

Solve this little system by hand, then check yourself with NumPy. Say two friends stop at a tea stall. Niranjan buys 3 chai and 2 samosas for 70 rupees. His friend Viraj buys 1 chai and 4 samosas for 90 rupees. Find the price of one chai and one samosa.

1. Write the two equations as a matrix A (the quantities) and a vector b (the totals, 70 and 90).
2. Call np.linalg.solve(A, b) and read off the two prices.
3. Verify your answer by computing A @ x and confirming it returns [70, 90].

Bonus: compute np.linalg.det(A) first. If it is not zero, you know a unique solution exists before you even solve. (The chai costs 10 and the samosa costs 20.)

Conclusion

You now have the core of NumPy linear algebra in your hands. You can compute dot products and matrix multiplications with @, flip a matrix with .T, read a determinant to tell if a matrix is invertible, undo a matrix with np.linalg.inv, solve a system with np.linalg.solve, fit a line with np.linalg.lstsq, and pull out eigenvalues and eigenvectors with np.linalg.eig. More importantly, you know the plain-English meaning behind each one, so a singular matrix or a mismatched shape now reads as a message instead of a mystery.

Next up is Pandas, where these grids of numbers grow labels and become the DataFrames you will actually load your data into. The linear algebra you learned here keeps running underneath: every model you fit later leans on solve, lstsq, and eig. Continue with the Pandas introduction, or browse the full Python + AI/ML tutorial series home to see where this fits in the bigger road to machine learning.

Frequently Asked Questions

Interview Questions on NumPy Linear Algebra

Scenario questions, not trivia: this is the form this topic takes in a real interview.

Q: What does the shape rule for matrix multiplication say, and what shape does (3, 4) @ (4, 2) produce?

The inner dimensions must match: for (m, n) @ (n, p) the two n values have to be equal, and the result has shape (m, p). So (3, 4) @ (4, 2) is valid because the inner 4s match, and it produces a (3, 2) matrix. If the inner numbers differ, NumPy raises a ValueError about aligned dimensions.

Q: How do you know whether a square matrix can be inverted?

Compute its determinant with np.linalg.det. If the determinant is non-zero the matrix is invertible; if it is zero the matrix is singular and has no inverse. In practice a determinant that is extremely close to zero also signals trouble, because the inverse will be numerically unstable even if it technically exists.

Q: Why prefer np.linalg.solve over computing np.linalg.inv(A) @ b?

solve uses a direct factorization (LU decomposition) that is both faster and more numerically stable than forming the full inverse and then multiplying. Building the inverse does extra work and accumulates rounding error, especially on large or ill-conditioned matrices. If all you need is x in Ax = b, reach for solve and never build the inverse just to throw it away.

Q: What is the difference between np.linalg.solve and np.linalg.lstsq?

solve requires a square, invertible coefficient matrix and returns an exact solution. lstsq handles overdetermined systems, where you have more equations than unknowns, and returns the least-squares best fit that minimizes the squared error. Linear regression on real data uses lstsq because you almost never get a perfectly square system.

Q: Scenario: you fit a linear regression and np.linalg.inv raises LinAlgError: Singular matrix. What do you check first?

A singular matrix means the determinant is zero, which almost always points to linearly dependent features: two columns carry the same information, for example one is a scaled copy or a sum of others. Check for duplicate or perfectly correlated columns and drop the redundant one, or switch to regularization like Ridge or Lasso. Using lstsq instead of a raw inverse also handles this gracefully because it does not require invertibility.

Q: Scenario: your teammate says the eigenvector your code prints is the negative of the one in the textbook. Is that a bug?

No. An eigenvector defines a direction, and both v and -v point along the same line, so either is a valid eigenvector for the same eigenvalue. NumPy is free to return whichever sign its algorithm lands on, and it can differ across machines or library versions. As long as A @ v equals lambda * v (check with np.allclose), the result is correct.

Q: In eigendecomposition of a covariance matrix, what does the largest eigenvalue tell you?

The largest eigenvalue marks the direction of greatest variance in your data, and its eigenvector is that direction, which PCA calls the first principal component. Dividing an eigenvalue by the sum of all eigenvalues gives the fraction of total variance that component captures. A few large eigenvalues mean you can compress the data to those directions with little loss of information.

Go deeper: NumPy official documentation covers every edge case of this topic.

Related Posts

Previous: NumPy: Broadcasting, Vectorization, ufuncs

Next: Pandas Introduction: Series, DataFrame, Reading Data

Series Home: Python + AI/ML Tutorial Series