Solution Sets and the Homogeneous Structure

Describing the Solution Set

When a linear system has infinitely many solutions, we need a compact description. The standard form is parametric: express each pivot variable in terms of the free variables, which become the parameters.

Example

Consider the reduced system:

$$\left[\begin{array}{rrrr|r} 1 & 2 & 0 & 1 & 3 \\ 0 & 0 & 1 & -1 & 2 \\ 0 & 0 & 0 & 0 & 0 \end{array}\right]$$

Pivot variables: $x_1$ (column 1) and $x_3$ (column 3). Free variables: $x_2$ and $x_4$.

Reading off the equations:

• $x_1 = 3 - 2x_2 - x_4$
• $x_3 = 2 + x_4$
• $x_2, x_4$ free (any value)

Parametric form (let $x_2 = s$, $x_4 = t$):

$$\vec{x} = \begin{pmatrix} 3 \\ 0 \\ 2 \\ 0 \end{pmatrix} + s\begin{pmatrix} -2 \\ 1 \\ 0 \\ 0 \end{pmatrix} + t\begin{pmatrix} -1 \\ 0 \\ 1 \\ 1 \end{pmatrix}, \quad s,t \in \mathbb{R}$$

This is the solution set expressed as an affine combination of vectors. Its geometric shape is a 2-dimensional flat (plane) in $\mathbb{R}^4$.

General = Particular + Homogeneous

The parametric form reveals a fundamental structural theorem. Let $\vec{p}$ be any particular solution (any single solution to $A\vec{x} = \vec{b}$) and let $\vec{h}$ range over all solutions of the homogeneous system $A\vec{x} = \vec{0}$.

Then the complete solution set of $A\vec{x} = \vec{b}$ is:

$$\{\vec{x} : A\vec{x} = \vec{b}\} = \{\vec{p} + \vec{h} : A\vec{h} = \vec{0}\}$$

Why: if $A\vec{p} = \vec{b}$ and $A\vec{h} = \vec{0}$, then $A(\vec{p} + \vec{h}) = A\vec{p} + A\vec{h} = \vec{b} + \vec{0} = \vec{b}$. And if $A\vec{x} = \vec{b}$ then $A(\vec{x} - \vec{p}) = \vec{0}$, so $\vec{x} - \vec{p}$ is a homogeneous solution.

Geometric picture: the solution set of $A\vec{x} = \vec{b}$ is the solution set of the homogeneous system translated by $\vec{p}$. It is a flat (affine subspace) passing through $\vec{p}$, parallel to the null space $\mathcal{N}(A)$.

The Homogeneous System

The system $A\vec{x} = \vec{0}$ always has at least one solution: $\vec{x} = \vec{0}$. This is the trivial solution.

The solution set $\mathcal{N}(A) = \{\vec{x} : A\vec{x} = \vec{0}\}$ is the null space of $A$. It is always a subspace (it contains $\vec{0}$ and is closed under addition and scaling).

When is the trivial solution the only one? When there are no free variables — i.e., every column of $A$ contains a pivot. For an $n \times n$ square matrix, this is equivalent to $A$ being invertible.

The Linear Combination Lemma

This core lemma formalizes why row operations preserve the solution set:

If $\vec{x}$ solves both $a_1 x_1 + \cdots + a_n x_n = b$ and $c_1 x_1 + \cdots + c_n x_n = d$, then it also solves every linear combination $k(a_1 x_1 + \cdots + a_n x_n) + m(c_1 x_1 + \cdots + c_n x_n) = kb + md$.

This means: any equation derivable from the original equations by linear combination is automatically satisfied by every solution of the original system. Row operations are just a systematic way of creating useful linear combinations.

Existence and Uniqueness

Existence ($A\vec{x} = \vec{b}$ has a solution): $\vec{b}$ must lie in the column space of $A$ — it must be expressible as a linear combination of the columns of $A$.

Uniqueness (the solution, if it exists, is unique): $\mathcal{N}(A) = \{\vec{0}\}$ — the only homogeneous solution is trivial.

For a square $n \times n$ matrix $A$:

This interplay between existence, uniqueness, rank, and determinant is the central theorem of linear systems — you will encounter it repeatedly as you study optimization and regression.

Computing the Decomposition

The parametric form derived above is not just pencil-and-paper algebra — you can read it directly into code. Using the same matrix:

import numpy as np

A  = np.array([[1., 2., 0.,  1.],
               [0., 0., 1., -1.],
               [0., 0., 0.,  0.]])
b  = np.array([3., 2., 0.])

# Particular solution: set free variables (x2, x4) to zero
p  = np.array([ 3., 0.,  2.,  0.])

# Null-space basis read directly from the parametric form
h1 = np.array([-2., 1.,  0.,  0.])  # coefficient of free variable x2
h2 = np.array([-1., 0.,  1.,  1.])  # coefficient of free variable x4

# Any x = p + s*h1 + t*h2 satisfies A @ x = b  (General = Particular + Homogeneous)
s, t = 3., -1.
x = p + s*h1 + t*h2
print('A @ x:', np.einsum('ij,j->i', A, x))        # [3. 2. 0.]
print('residual:', np.einsum('ij,j->i', A, x) - b)  # [0. 0. 0.]

# Programmatic null space via SVD (when you cannot read it off by hand)
_, S, Vh = np.linalg.svd(A, full_matrices=True)
rank       = int(np.sum(S > 1e-9))
null_basis = Vh[rank:]   # rows are null vectors; shape (2, 4)
print('null basis (rows):')
print(np.round(null_basis, 4))

The einsum string 'ij,j->i' is identical in NumPy, PyTorch (torch.einsum), and TensorFlow (tf.einsum). The SVD-based null space is in an orthonormal basis rather than the free-variable basis, but it spans the same subspace as the parametric vectors above. For large systems, np.linalg.lstsq (or its PyTorch/TF equivalents) finds a minimum-norm particular solution automatically, without needing to identify free variables by hand.