How to Solve Systems of Linear Equations (The Complete Guide)

How to Truly Solve Linear Equations (The Math & The Code)

A system of linear equations isn't just a problem from high school algebra; it's the hidden engine that powers a huge part of machine learning. The goal is simple: find a unique set of values that satisfy all constraints simultaneously. Visually, it's about finding the exact point where the lines cross.

But how do we solve this rigorously? This post provides a more detailed breakdown of the two primary methods: the procedural, step-by-step approach of Gaussian Elimination, and the high-level, conceptual solution using the Matrix Inverse.

Watch the video for the full visual explanation, then scroll down for the detailed examples and code.

The Mathematician's Way: Gaussian Elimination to RREF

To solve a system systematically, we first represent it as an Augmented Matrix. This structure combines the coefficients (A) and the constant outcomes (b) into a single matrix: `[A | b]`.

The goal of Gaussian Elimination is to transform this matrix into Reduced Row Echelon Form (RREF), where the solution becomes obvious. The RREF form looks like `[I | x]`, where `I` is the Identity Matrix and `x` is our solution vector.

We achieve this using three Elementary Row Operations:

1. Row Swapping (Rᵢ ↔ Rⱼ): Swap any two rows.
2. Row Scaling (k * Rᵢ → Rᵢ): Multiply any row by a non-zero scalar.
3. Row Addition (Rᵢ + k * Rⱼ → Rᵢ): Add a multiple of one row to another row.

Example 1: A Unique Solution

Let's solve the system from the video: `2x + y = 4` and `x - y = -1`.

Initial Augmented Matrix:
[ 2  1 |  4 ]
[ 1 -1 | -1 ]

Step 1: R1 ↔ R2
[ 1 -1 | -1 ]
[ 2  1 |  4 ]

Step 2: R2 → R2 - 2*R1
[ 1 -1 | -1 ]
[ 0  3 |  6 ]

Step 3: R2 → (1/3)*R2
[ 1 -1 | -1 ]
[ 0  1 |  2 ]

Step 4: R1 → R1 + R2
[ 1  0 | 1 ]
[ 0  1 | 2 ]  --> RREF Achieved.
    

The final matrix tells us that x = 1 and y = 2. We have found the unique solution.

The Engineer's Insight: When Solutions Don't Exist (or Aren't Unique)

The Matrix Inverse method (`x = A⁻¹b`) is elegant, but it fails if the matrix `A` is singular (i.e., its determinant is zero). This happens when the system does not have a unique solution. Gaussian Elimination, however, handles these cases gracefully and reveals the nature of the problem.

Example 2: No Solution (Inconsistent System)

Consider the system: `x + y = 2` and `x + y = 4`. Visually, these are parallel lines that never intersect.

Initial Augmented Matrix:
[ 1  1 | 2 ]
[ 1  1 | 4 ]

Step 1: R2 → R2 - R1
[ 1  1 | 2 ]
[ 0  0 | 2 ]
    

Look at the last row. It translates to the equation `0x + 0y = 2`, which simplifies to `0 = 2`. This is a logical contradiction. When Gaussian Elimination results in a row of this form, it proves that the system is inconsistent and has no solution.

Example 3: Infinite Solutions (Dependent System)

Now consider: `x + y = 2` and `2x + 2y = 4`. These are two different equations for the exact same line.

Initial Augmented Matrix:
[ 1  1 | 2 ]
[ 2  2 | 4 ]

Step 1: R2 → R2 - 2*R1
[ 1  1 | 2 ]
[ 0  0 | 0 ]
    

The last row now translates to `0x + 0y = 0`, which simplifies to `0 = 0`. This is always true, but it provides no new information. This indicates that one of our equations was redundant. The remaining equation, `x + y = 2`, has an infinite number of solutions (e.g., if x=0, y=2; if x=1, y=1, etc.). This system has infinite solutions.

The Engineer's Way: Python & NumPy

In practice, we use libraries like NumPy to solve these systems. The key is to use the right tool for the job.

import numpy as np

# Our solvable system from Example 1
A = np.array([[2, 1], [1, -1]])
b = np.array([4, -1])

# Method 1: Direct Inverse (The "Sledgehammer")
# Fails on singular matrices.
try:
  x_inv = np.linalg.inv(A) @ b
  print(f"Inverse Method Solution: {x_inv}")
except np.linalg.LinAlgError:
  print("Matrix is singular. Inverse method fails.")

# Method 2: The Preferred Solver (More Robust)
# This is the function you should use in production code.
try:
  x_solve = np.linalg.solve(A, b)
  print(f"Solve Method Solution: {x_solve}")
except np.linalg.LinAlgError:
  print("Matrix is singular. System has no unique solution.")
    

While both methods give the same answer for a solvable system, np.linalg.solve() is the professional choice. It is faster, more numerically stable for ill-conditioned matrices, and its underlying algorithms (often based on LU Decomposition) are more robust than manually computing the inverse.

Conclusion

Understanding both the procedural "how" of Gaussian Elimination and the conceptual "what" of the Matrix Inverse gives you the complete picture. You can now appreciate the elegance of the theory and the robust engineering required to apply it in practice. This foundation is essential for truly understanding how models like Linear Regression work under the hood.

Your Turn...

Which method was more intuitive for you—the step-by-step row operations or the high-level inverse abstraction? Let me know your thoughts in the comments!

This post is part of the "Linear Algebra: The Core of Machine Learning" series. For the previous part, check out: Matrix Inverse Explained - For Developers & ML Beginners.