Visighta

Section 1.4 — Infinite-Solution Systems and Free Variables

In this section we explore what happens when a linear system has infinitely many solutions. We will:

  • Recognise when a system has more than one solution.
  • Distinguish dependent (pivot) variables from free (non-pivot) variables.
  • Write solutions in parametric form and interpret them as lines or planes.
  • Use NumPy to set up, solve, and verify linear systems, including singular cases.
  • Watch Manim animations that bring the geometry to life.
import numpy as np
import sympy as sy
import matplotlib.pyplot as plt
from IPython.display import Video, display
import time

sy.init_printing()
np.set_printoptions(precision=3, suppress=True)

Setting Up a Linear System $Ax = b$

Any system of linear equations can be written in matrix form as

$$A\mathbf{x} = \mathbf{b}$$

where $A$ is the coefficient matrix, $\mathbf{x}$ is the vector of unknowns, and $\mathbf{b}$ is the right-hand side vector.

For example, the system \begin{align} x + 2y &= 6 \ 2x + 4y &= 12 \end{align}

becomes $A\mathbf{x} = \mathbf{b}$ with

$$ A = \begin{bmatrix} 1 & 2 \ 2 & 4 \end{bmatrix}, \quad \mathbf{x} = \begin{bmatrix} x \ y \end{bmatrix}, \quad \mathbf{b} = \begin{bmatrix} 6 \ 12 \end{bmatrix} $$

Notice that the second equation is just twice the first — so both equations describe the same line. This hints that the system has infinitely many solutions.

A = np.array([[1, 2],
              [2, 4]])
b = np.array([6, 12])

print("A ="); print(A)
print("\nb =", b)

We can also verify this is singular by checking its determinant:

np.linalg.det(A)  # zero (or very close) → singular matrix

Generating Random Systems for Testing

When we work with a non-singular (full-rank) square matrix, np.linalg.solve succeeds. Random matrices are almost always non-singular, which makes them handy for quick tests.

np.random.seed(42)  # for reproducibility
A_rand = np.round(10 * np.random.rand(5, 5))
b_rand = np.round(10 * np.random.rand(5))

print("A ="); print(A_rand)
print("\nb =", b_rand)

Solving the System Using np.linalg.solve

For a square, non-singular system $A\mathbf{x} = \mathbf{b}$, NumPy's np.linalg.solve computes the exact solution. Internally it performs an LU decomposition — the same row-reduction idea we use by hand, but highly optimised via LAPACK.

x_rand = np.linalg.solve(A_rand, b_rand); x_rand

Verifying the Solution by Computing the Residual

After solving, always check by computing the residual $A\mathbf{x} - \mathbf{b}$. If the answer is correct, the residual should be a vector of (near) zeros. Due to floating-point arithmetic, you may see tiny values like $-0.0$ — that is perfectly normal.

A_rand @ x_rand - b_rand  # should be (approximately) zero

The entries above are essentially zero — the tiny residuals are due to floating-point round-off, not actual errors.

Handling Singular and Near-Singular Systems

When $A$ is singular (its determinant is zero), the system either has no solution or infinitely many solutions. In either case, np.linalg.solve cannot find a unique answer and raises a LinAlgError.

Let us try to solve the singular system from earlier:

A_singular = np.array([[1, 2],
                       [2, 4]])
b_singular = np.array([6, 12])

try:
    np.linalg.solve(A_singular, b_singular)
except np.linalg.LinAlgError as e:
    print(f"LinAlgError: {e}")

NumPy raises an error because the matrix is singular — but that does not mean there is no solution! It only means there is no unique solution. For singular systems we need SymPy (row reduction) to describe the full family of solutions.

Near-singular matrices and condition numbers

A matrix can be close to singular without being exactly singular. The condition number $\kappa(A) = |A||A^{-1}|$ measures this: larger values mean the solution is more sensitive to tiny changes in $A$ or $\mathbf{b}$.

# A well-conditioned matrix
print("Condition number of A_rand:", np.linalg.cond(A_rand))

# A near-singular matrix (rows almost identical)
A_near = np.array([[1.0, 2.0],
                   [1.0, 2.0001]])
print("Condition number of A_near:", np.linalg.cond(A_near))

The near-singular matrix has a very large condition number — a warning sign that numerical solutions may be unreliable.

Solving with np.linalg.inv — and Why solve is Preferred

We could solve $A\mathbf{x} = \mathbf{b}$ by computing $\mathbf{x} = A^{-1}\mathbf{b}$. Let's try both approaches and compare.

x_solve = np.linalg.solve(A_rand, b_rand)
x_inv   = np.linalg.inv(A_rand) @ b_rand

print("Via solve:", x_solve)
print("Via inv:  ", x_inv)
print("\nDifference:", x_solve - x_inv)

Both methods give the same answer, but np.linalg.solve is preferred because:

  1. Numerical stability — solve uses LU decomposition with partial pivoting, which is more stable than explicitly computing $A^{-1}$.
  2. Efficiency — Computing $A^{-1}$ requires roughly three times more floating-point operations than LU-based solving.

Rule of thumb: never compute $A^{-1}$ just to multiply by $\mathbf{b}$; use solve instead.

Geometric Interpretation — Manim Animations

The animations below are imported from the course's Manim module. They visualise the key ideas of this section:

  • Why some systems have infinitely many solutions.
  • How to spot free versus dependent variables.
  • How to write the solution in parametric form.
  • The geometry of solution lines and planes.

Run each cell to render and display the animation inline.

from manim import *

# Import all scene classes from the section 4 animation module
from linear_algebra_course.chapter1.section4.scenes import *

config.media_width = "80%"
config.verbosity = "WARNING"
from drawsvg_renderer import DrawSVGRenderer

renderer = DrawSVGRenderer(
    width=800,
    height=600,
    frame_width=16.0,
    background_color="#000000",
    output_mode='jupyter',
    fps=30,
    show_progress=False
)

Scene 1 — A System with Many Solutions

Consider the system \begin{align} x + 2y &= 6 \ 2x + 4y &= 12 \end{align}

The second equation is just twice the first, so both equations describe the same line. Every point on that line is a solution — there are infinitely many.

renderer.display_all(Scene1_ManySolutions(), display=False)
renderer._finalize_interactivity()
renderer.display_inline()
Cell visualization output

Scene 2 — Dependent versus Free Variables

Row-reduce the augmented matrix to echelon form:

$$ \begin{bmatrix} 1 & 2 & | & 6 \ 2 & 4 & | & 12 \end{bmatrix} \xrightarrow{R_2 \to R_2 - 2R_1} \begin{bmatrix} 1 & 2 & | & 6 \ 0 & 0 & | & 0 \end{bmatrix} $$

  • Column 1 has a pivot → $x$ is a dependent variable.
  • Column 2 has no pivot → $y$ is a free variable.

The free variable $y$ can take any value; $x$ is then determined by $x = 6 - 2y$.

renderer.display_all(Scene2_DependentVsFree(), display=False)
renderer._finalize_interactivity()
renderer.display_inline()
Cell visualization output

We can verify the row reduction with SymPy:

M = sy.Matrix([[1, 2, 6],
               [2, 4, 12]])
M_rref = M.rref(); M_rref

The pivot is in column $0$ (for $x$). Column $1$ ($y$) is free — confirming what the animation showed.

Scene 3 — Free Variable as a Parameter

Replace the free variable $y$ with a parameter $t$:

$$y = t, \quad x = 6 - 2t$$

As $t$ varies, the point $(x, y) = (6-2t,\; t)$ traces out the solution line. The animation below shows a dot sliding along the line as $t$ changes, with a live readout of $t$, $x$, and $y$.

renderer.display_all(Scene3_ParameterT(), display=False)
renderer._finalize_interactivity()
renderer.display_inline()
Cell visualization output

Scene 4 — Parametric Solution as a Vector

We can rewrite the parametric solution in vector form:

$$ \begin{bmatrix} x \ y \end{bmatrix} = \underbrace{\begin{bmatrix} 6 \ 0 \end{bmatrix}}{\text{start point}} + t \underbrace{\begin{bmatrix} -2 \ 1 \end{bmatrix}} $$}

Every point on the solution line equals the start point $(6, 0)$ plus $t$ times the direction vector $(-2, 1)$. This is a fundamental pattern you will see throughout linear algebra.

renderer.display_all(Scene4_ParametricVectorForm(), display=False)
renderer._finalize_interactivity()
renderer.display_inline()
Cell visualization output

Let's verify the parametric form with a quick plot. Choose a few values of $t$ and check they all lie on the line $x + 2y = 6$.

t_vals = np.array([0, 1, 2, 3])
start = np.array([6, 0])
direction = np.array([-2, 1])

solutions = start + t_vals[:, None] * direction
print("t | x    y  | x + 2y")
print("-" * 28)
for t, (x, y) in zip(t_vals, solutions):
    print(f"{t} | {x:4.0f}  {y:3.0f}  | {x + 2*y:.0f}")

Every row satisfies $x + 2y = 6$ — confirming our parametric solution is correct.

Scene 5 — Plane of Solutions in Three Variables

Now consider a single equation in three unknowns:

$$x + y + z = 3$$

There are two non-pivot variables ($y$ and $z$), so we introduce two free parameters: $y = s$, $z = t$. Then $x = 3 - s - t$.

The solution set is a plane in $\mathbb{R}^3$. The animation shows a dot sweeping across this plane as $s$ and $t$ vary.

Key insight: one free variable → a line; two free variables → a plane.

renderer.display_all(Scene5_PlaneOfSolutions(), display=False)
renderer._finalize_interactivity()
renderer.display_inline()
Cell visualization output

We can also visualise this plane with Matplotlib:

S, T = np.mgrid[0:3:21j, 0:3:21j]
X_plane = 3 - S - T

fig = plt.figure(figsize=(9, 9))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X_plane, S, T, cmap='viridis', alpha=0.6)
ax.set_xlabel('x'); ax.set_ylabel('y (= s)'); ax.set_zlabel('z (= t)')
ax.set_title('Solution plane of $x + y + z = 3$', size=16)
plt.show()

Scene 6 — Intuition Summary

The overall workflow for infinite-solution systems:

  1. Row-reduce the augmented matrix.
  2. Identify non-pivot (free) variables.
  3. Assign parameters ($t$, $s$, ...) to free variables.
  4. Express pivot variables in terms of those parameters.
  5. Write the solution family — this gives you a line (1 free variable), a plane (2 free variables), etc.
renderer.display_all(Scene6_IntuitionSummary(), display=False)
renderer._finalize_interactivity()
renderer.display_inline()
Cell visualization output

Solving Larger Random Systems and Timing Performance

NumPy's linalg.solve relies on LAPACK's highly optimised LU factorisation, which runs in $O(n^3)$ time. Let's see how it scales on random systems of increasing size.

sizes = [10, 100, 500, 1000]

for n in sizes:
    A_big = np.random.rand(n, n)
    b_big = np.random.rand(n)
    start = time.perf_counter()
    x_big = np.linalg.solve(A_big, b_big)
    elapsed = time.perf_counter() - start
    residual = np.linalg.norm(A_big @ x_big - b_big)
    print(f"n = {n:5d}  |  time = {elapsed:.6f} s  |  residual norm = {residual:.2e}")

Even $1000 \times 1000$ systems are solved in a fraction of a second — thanks to LAPACK's compiled Fortran code running under the hood. The residual norms stay tiny, confirming numerical accuracy.

Summary of Section 1.4

Concept Key Idea
Singular matrix $\det(A) = 0$ → no unique solution
Free variables Non-pivot columns after row reduction
Parametric form Assign parameters to free variables; express pivot variables in terms of them
Geometry 1 free var → line, 2 free vars → plane
np.linalg.solve Fast, stable solver for non-singular systems
Residual check Always verify via $|A\mathbf{x} - \mathbf{b}|$