Systems of Linear Equations

Introduction

Systems of linear equations arise whenever several unknown quantities are constrained by simultaneous linear relationships. In artificial intelligence and machine learning, they appear in linear regression, feature transformations, optimization, signal reconstruction, graph algorithms, state estimation, computer vision, and neural-network computations.

A linear system may have one solution, infinitely many solutions, or no solution. Determining which case applies is as important as calculating numerical values. This chapter develops the theory and algorithms needed to formulate systems in matrix form, solve them by row reduction and Gaussian elimination, interpret rank and nullity, and assess computational cost and numerical reliability.

Learning Objectives

After completing this chapter, students should be able to:

• formulate simultaneous linear equations as a matrix equation;
• perform elementary row operations and construct echelon forms;
• solve systems using Gaussian and Gauss–Jordan elimination;
• classify systems as consistent, inconsistent, unique, or underdetermined;
• use rank and nullity to characterize solution sets;
• express infinitely many solutions in parametric vector form;
• estimate the computational cost of dense elimination;
• explain pivoting, conditioning, residuals, and numerical stability;
• select appropriate methods for dense, sparse, or least-squares problems.

1. Linear Systems and Matrix Form

A system of $m$ linear equations in $n$ unknowns has the form

$$

\begin{aligned}

a_{11}x_1+\cdots+a_{1n}x_n&=b_1,\\

a_{21}x_1+\cdots+a_{2n}x_n&=b_2,\\

&\vdots\\

a_{m1}x_1+\cdots+a_{mn}x_n&=b_m.

\end{aligned}

$$

It can be written compactly as

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

where $A\in\mathbb{R}^{m\times n}$ is the coefficient matrix, $\mathbf{x}\in\mathbb{R}^n$ is the unknown vector, and $\mathbf{b}\in\mathbb{R}^m$ is the right-hand-side vector.

The augmented matrix is

$$[A\mid \mathbf{b}],$$

which places the coefficients and right-hand side in one array. For example,

$$

\begin{aligned}

2x-y&=3,\\

x+4y&=7

\end{aligned}

$$

is represented by

$$

\left[

\begin{array}{cc|c}

2&-1&3\\

1&4&7

\end{array}

\right].

$$

A vector $\mathbf{x}$ is a solution if substituting it into every equation produces the corresponding component of $\mathbf{b}$.

Geometric interpretation

Each equation in two variables represents a line, while each equation in three variables represents a plane. A solution is a common intersection. Two distinct lines may intersect once, be parallel, or coincide, corresponding respectively to a unique solution, no solution, or infinitely many solutions.

2. Elementary Row Operations

Three operations may be performed on an augmented matrix without changing its solution set:

1. interchange two rows;

2. multiply a row by a nonzero scalar;

3. add a scalar multiple of one row to another row.

Each operation corresponds to replacing equations by equivalent equations. Algebraically, a row operation is multiplication by an invertible elementary matrix, so it preserves equivalence of the system.

Row-echelon form

A matrix is in row-echelon form (REF) when:

• all zero rows appear below nonzero rows;
• each leading nonzero entry, called a pivot, lies to the right of the pivot above it;
• entries below each pivot are zero.

A matrix is in reduced row-echelon form (RREF) when it is in REF and, additionally:

• every pivot equals one;
• each pivot is the only nonzero entry in its column.

REF is sufficient for solving by back substitution. RREF directly reveals pivot and free variables. Although REF is not unique, every matrix has a unique RREF.

3. Gaussian and Gauss–Jordan Elimination

Gaussian elimination transforms a system into upper-triangular or row-echelon form. The resulting system is solved by back substitution.

For pivot column $k$, the basic elimination update is

$$m_{ik}=\frac{a_{ik}}{a_{kk}},\qquad

R_i\leftarrow R_i-m_{ik}R_k,$$

for rows $i>k$, provided the pivot $a_{kk}$ is nonzero. If it is zero or numerically unsuitable, rows should be interchanged.

Gauss–Jordan elimination continues above each pivot until RREF is obtained. It simplifies interpretation but generally requires more arithmetic than Gaussian elimination when only a solution is needed.

Algorithmic outline

1. Select a pivot in the current column.

2. If necessary, interchange rows.

3. Eliminate entries below the pivot.

4. Repeat for successive columns.

5. Use back substitution, or continue to RREF for Gauss–Jordan elimination.

4. Worked Example: A Unique Solution

Solve

$$

\begin{aligned}

3x+y+z&=8,\\

x+2y-z&=3,\\

3x+y-z&=0.

\end{aligned}

$$

The augmented matrix is

$$

\left[

\begin{array}{ccc|c}

3&1&1&8\\

1&2&-1&3\\

3&1&-1&0

\end{array}

\right].

$$

Interchange the first two rows:

$$

\left[

\begin{array}{ccc|c}

1&2&-1&3\\

3&1&1&8\\

3&1&-1&0

\end{array}

\right].

$$

Apply $R_2\leftarrow R_2-3R_1$ and $R_3\leftarrow R_3-3R_1$:

$$

\left[

\begin{array}{ccc|c}

1&2&-1&3\\

0&-5&4&-1\\

0&-5&2&-9

\end{array}

\right].

$$

Then apply $R_3\leftarrow R_3-R_2$:

$$

\left[

\begin{array}{ccc|c}

1&2&-1&3\\

0&-5&4&-1\\

0&0&-2&-8

\end{array}

\right].

$$

Back substitution gives

$$z=4,$$

$$-5y+4(4)=-1\quad\Rightarrow\quad y=\frac{17}{5},$$

and

$$x+2\left(\frac{17}{5}\right)-4=3

\quad\Rightarrow\quad x=\frac15.$$

Therefore,

$$\boxed{(x,y,z)=\left(\frac15,\frac{17}{5},4\right)}.$$

Verification in the original equations gives

$$3\left(\frac15\right)+\frac{17}{5}+4=8,$$

$$\frac15+2\left(\frac{17}{5}\right)-4=3,$$

$$3\left(\frac15\right)+\frac{17}{5}-4=0.$$

All equations are satisfied.

5. Consistency and Rank

The rank of a matrix is the number of pivot positions in its echelon form. Equivalently, it is the dimension of its column space or row space.

A system is consistent if it has at least one solution. The rank criterion, sometimes called the Rouché–Capelli theorem, states that

$$A\mathbf{x}=\mathbf{b}\text{ is consistent}

\iff \operatorname{rank}(A)=\operatorname{rank}([A\mid\mathbf{b}]).$$

For a system with $n$ unknowns:

• it has a unique solution when

$$\operatorname{rank}(A)=\operatorname{rank}([A\mid\mathbf{b}])=n;$$

• it has infinitely many solutions when

$$\operatorname{rank}(A)=\operatorname{rank}([A\mid\mathbf{b}])<n;$$

• it is inconsistent when

$$\operatorname{rank}(A)<\operatorname{rank}([A\mid\mathbf{b}]).$$

For a square matrix $A\in\mathbb{R}^{n\times n}$, a unique solution exists for every $\mathbf{b}$ exactly when $A$ is invertible, or equivalently when $\operatorname{rank}(A)=n$.

Example: An inconsistent system

Consider

$$x+y=2,\qquad 2x+2y=5.$$

Elimination gives

$$

\left[

\begin{array}{cc|c}

1&1&2\\

2&2&5

\end{array}

\right]

\longrightarrow

\left[

\begin{array}{cc|c}

1&1&2\\

0&0&1

\end{array}

\right].

$$

The second row represents $0=1$, which is impossible. Here $\operatorname{rank}(A)=1$ but $\operatorname{rank}([A\mid\mathbf{b}])=2$.

6. Underdetermined Systems, Free Variables, and Nullity

A system is underdetermined when it has fewer independent equations than unknowns. If it is consistent, at least one variable is free and the system has infinitely many solutions over $\mathbb{R}$.

Consider

$$

\begin{aligned}

x+2y-z&=3,\\

2x+4y-2z&=6.

\end{aligned}

$$

The second equation is twice the first, so the rank is one. Let $y=s$ and $z=t$. Then

$$x=3-2s+t.$$

Thus

$$

\begin{bmatrix}x\\y\\z\end{bmatrix}

=

\begin{bmatrix}3\\0\\0\end{bmatrix}

+s\begin{bmatrix}-2\\1\\0\end{bmatrix}

+t\begin{bmatrix}1\\0\\1\end{bmatrix},

\qquad s,t\in\mathbb{R}.

$$

The first vector is a particular solution. The remaining vectors form a basis for the null space

$$\mathcal{N}(A)=\{\mathbf{v}:A\mathbf{v}=\mathbf{0}\}.$$

For $A\in\mathbb{R}^{m\times n}$, the rank–nullity theorem states

$$\operatorname{rank}(A)+\operatorname{nullity}(A)=n.$$

Therefore, this example has nullity $3-1=2$, matching its two free parameters.

More generally, if $A\mathbf{x}=\mathbf{b}$ is consistent and $\mathbf{x}_p$ is one particular solution, then every solution has the form

$$\mathbf{x}=\mathbf{x}_p+\mathbf{z},\qquad \mathbf{z}\in\mathcal{N}(A).$$

The solution set is an affine translation of the null space. It is not generally a subspace unless $\mathbf{b}=\mathbf{0}$.

Selecting one solution

Applications may require a single solution from an infinite set. A common choice is the minimum-Euclidean-norm solution. If $A$ has full row rank and $m<n$, it is

$$\mathbf{x}_{\min}=A^T(AA^T)^{-1}\mathbf{b}.$$

In general, the Moore–Penrose pseudoinverse gives $\mathbf{x}_{\min}=A^+\mathbf{b}$. Additional criteria such as nonnegativity, sparsity, or physical bounds require constrained optimization rather than arbitrary assignment of free variables.

7. LU Factorization and Repeated Systems

Gaussian elimination can be organized as a factorization. Without row interchanges, it produces

$$A=LU,$$

where $L$ is lower triangular and contains elimination multipliers, while $U$ is upper triangular. With partial pivoting, the standard form is

$$PA=LU,$$

where $P$ is a permutation matrix recording row exchanges.

To solve $A\mathbf{x}=\mathbf{b}$:

1. solve $L\mathbf{y}=P\mathbf{b}$ by forward substitution;

2. solve $U\mathbf{x}=\mathbf{y}$ by back substitution.

LU factorization is particularly valuable when many systems share the same coefficient matrix but have different right-hand sides. The expensive factorization is performed once, after which each new right-hand side requires only triangular solves.

8. Computational Cost

For a dense $n\times n$ matrix, elimination updates the trailing submatrix at every pivot. Counting one multiplication or one addition as one floating-point operation, the dominant work is

$$2\sum_{j=1}^{n-1}j^2

=\frac{(n-1)n(2n-1)}{3}

=\frac{2}{3}n^3+O(n^2).$$

Thus, dense Gaussian or LU factorization requires approximately

$$\boxed{\frac{2}{3}n^3\text{ floating-point operations}}.$$

Back substitution alone requires approximately $n^2$ operations. After LU factorization, one forward and one backward triangular solve require about $2n^2$ operations per right-hand side. Dense storage requires $O(n^2)$ memory.

For a general dense $m\times n$ matrix, elimination has complexity on the order of

$$O\bigl(mn\min(m,n)\bigr).$$

The cubic growth for square matrices becomes expensive at large scale. Sparse matrices should therefore be processed using sparse factorizations that avoid unnecessary operations on zeros. However, elimination may create new nonzero entries, called fill-in. Reordering rows and columns can substantially reduce fill-in and memory usage.

Parallel processors and GPUs accelerate matrix updates because much of factorization can be expressed through block matrix operations. Mixed-precision solvers may perform much of the work at lower precision and then improve the answer through iterative refinement, provided conditioning and hardware precision permit reliable correction.

9. Numerical Stability

Exact arithmetic and floating-point arithmetic behave differently. Computers store rounded approximations, so subtraction, division, and elimination introduce rounding errors.

Small pivots and partial pivoting

Dividing by a very small pivot produces large multipliers and may amplify rounding errors. Partial pivoting selects the entry with the largest absolute value in the current column among the available rows and swaps it into the pivot position.

For example, consider

$$

\begin{aligned}

0.0001x+y&=1,\\

x+y&=2.

\end{aligned}

$$

Using $0.0001$ as the first pivot creates the multiplier $10{,}000$. Interchanging the rows uses the safer pivot $1$. The exact solution is

$$x=\frac{1}{0.9999}\approx1.00010001,

\qquad

y\approx0.99989999.$$

Partial pivoting is effective for most practical dense systems, although exceptional matrices can still exhibit substantial element growth. Complete pivoting searches both rows and columns but costs more and requires tracking column permutations.

Conditioning versus stability

Conditioning describes sensitivity inherent in the mathematical problem. For an invertible matrix under a chosen induced norm,

$$\kappa(A)=\lVert A\rVert\lVert A^{-1}\rVert$$

is its condition number. A large $\kappa(A)$ means that small input perturbations may cause large changes in the solution. If $A$ is fixed and only $\mathbf{b}$ changes, then

$$

\frac{\lVert\delta\mathbf{x}\rVert}{\lVert\mathbf{x}\rVert}

\leq

\kappa(A)

\frac{\lVert\delta\mathbf{b}\rVert}{\lVert\mathbf{b}\rVert},

$$

for nonzero quantities in compatible norms.

Stability, by contrast, concerns the algorithm. A backward-stable method returns the exact solution of a nearby problem. A stable algorithm cannot remove sensitivity caused by severe ill-conditioning.

Residual and error

For a computed solution $\widehat{\mathbf{x}}$, the residual is

$$\mathbf{r}=\mathbf{b}-A\widehat{\mathbf{x}}.$$

A small residual indicates that the computed vector nearly satisfies the equations, but it does not necessarily imply a small solution error when $A$ is ill-conditioned. A useful scale-independent measure is the relative backward error

$$

\eta=

\frac{\lVert\mathbf{r}\rVert}

{\lVert A\rVert\lVert\widehat{\mathbf{x}}\rVert+\lVert\mathbf{b}\rVert}.

$$

Practical validation should consider residuals, conditioning, data uncertainty, and domain constraints together.

10. Overdetermined Systems and Least Squares

When $m>n$, exact consistency is often impossible because noisy observations do not intersect at a common solution. Instead, one may minimize the residual norm:

$$\min_{\mathbf{x}}\lVert A\mathbf{x}-\mathbf{b}\rVert_2.$$

If $A$ has full column rank, the least-squares solution satisfies the normal equations

$$A^TA\mathbf{x}=A^T\mathbf{b}.$$

Although mathematically valid, explicitly forming $A^TA$ approximately squares the condition number. QR factorization is therefore generally preferred for numerical reliability, while singular value decomposition is especially useful for rank-deficient or ill-conditioned problems.

Least squares is fundamental to linear regression. However, solving an inconsistent system approximately must not be confused with finding an exact solution by elimination.

11. Illustrative Case Study: Calibrated Sensor Fusion

Background

This fictional case illustrates a linear sensor-reconstruction problem. A monitoring device receives signals from three latent sources. Calibration has produced dimensionless linear response coefficients. Some coefficients are negative because the readings are baseline-centered differences rather than raw physical measurements.

Problem

Let $x$, $y$, and $z$ denote source intensities in calibrated intensity units. The three sensor outputs are

$$

A=

\begin{bmatrix}

3&1&1\\

1&2&-1\\

3&1&-1

\end{bmatrix},

\qquad

\mathbf{b}=

\begin{bmatrix}8\\3\\0\end{bmatrix}.

$$

The reconstruction model is $A\mathbf{x}=\mathbf{b}$.

Analysis and approach

Row reduction gives an upper-triangular system with three pivots. Therefore,

$$\operatorname{rank}(A)=3,$$

so the square coefficient matrix is invertible and the model has a unique solution. Gaussian elimination gives

$$

\mathbf{x}=

\begin{bmatrix}

1/5\\17/5\\4

\end{bmatrix}

=

\begin{bmatrix}

0.2\\3.4\\4.0

\end{bmatrix}.

$$

Direct verification yields

$$3(0.2)+3.4+4=8,$$

$$0.2+2(3.4)-4=3,$$

$$3(0.2)+3.4-4=0.$$

The recovered intensities are also nonnegative, as required by their interpretation.

Outcome

The system reconstructs all three source intensities exactly under the calibrated linear model. In an operational implementation, the organization should also report residuals and monitor the condition number. A unique algebraic solution is not necessarily trustworthy if calibration coefficients are inaccurate or the system is poorly conditioned.

Learning points

• Matrix formulation separates calibration coefficients from measured outputs.
• Rank determines whether source intensities are identifiable.
• Substitution into the original equations verifies arithmetic correctness.
• Physical constraints provide an additional validation layer.
• Uncertainty in sensors and calibration must be considered alongside algebraic solvability.

12. Applications in AI and Machine Learning

Linear systems support many AI workflows:

• Linear regression: least-squares fitting estimates model coefficients.
• Computer vision: geometric reconstruction and image alignment produce structured systems.
• Graph learning: ranking and diffusion methods can require sparse linear solves.
• Optimization: Newton and interior-point methods repeatedly solve systems involving Hessian or constraint matrices.
• Gaussian models: inference may involve covariance or precision-matrix systems.
• Signal processing: source separation, filtering, and inverse problems rely on linear models.

For large AI problems, practitioners typically avoid explicitly computing $A^{-1}$. Solving through factorization is faster and usually more stable. Iterative methods such as conjugate gradients may be preferable for large sparse symmetric positive-definite systems, often with a preconditioner.

13. Limitations, Governance, and Best Practices

A linear model assumes that output variables are linear combinations of inputs. Real systems may include nonlinear responses, interactions, drift, missing variables, and measurement bias. A precise numerical solution cannot compensate for an invalid model.

When solutions affect people or regulated processes, practitioners should document data origins, calibration procedures, precision, constraints, and uncertainty. Sensitive applications should include audit trails and independent validation. A low residual does not establish fairness, causality, or ethical acceptability.

Recommended practices include:

1. Check dimensions and units before solving.

2. Use pivoting rather than unpivoted elimination in floating-point computation.

3. Reuse LU factors when solving for multiple right-hand sides.

4. Use QR or SVD for difficult least-squares problems.

5. Exploit sparsity and matrix structure.

6. Avoid forming a matrix inverse solely to solve a system.

7. Compute and interpret residuals.

8. Estimate conditioning when reliability matters.

9. Apply scaling when rows or columns have highly different magnitudes.

10. Validate solutions against physical, statistical, and governance constraints.

Chapter Summary

A linear system is represented as $A\mathbf{x}=\mathbf{b}$ and solved by equivalent elementary row operations. Gaussian elimination produces echelon form followed by back substitution, while Gauss–Jordan elimination produces reduced row-echelon form. Rank determines consistency and uniqueness: equal coefficient and augmented ranks imply consistency, while full column rank gives uniqueness. Consistent rank-deficient systems have free variables, and their solutions form an affine translation of the null space.

Dense Gaussian or LU factorization requires approximately $(2/3)n^3$ floating-point operations under the conventional flop count. Partial pivoting reduces the risk of error amplification, but reliability also depends on the condition number of the problem. Residuals, backward error, model validity, and domain constraints should all be evaluated. In AI and machine learning, these principles support regression, optimization, sensor reconstruction, graph methods, and many other computational tasks.

Key Terms

Augmented matrix: A matrix formed by appending the right-hand side to the coefficient matrix.

Back substitution: Solving an upper-triangular system from the last equation upward.

Condition number: A measure of the sensitivity of a problem's solution to input perturbations.

Consistency: The property that a system has at least one solution.

Elementary row operation: A reversible row interchange, scaling, or row replacement.

Free variable: A non-pivot variable that can be assigned parameter values.

Gaussian elimination: Reduction to echelon or triangular form followed by back substitution.

Gauss–Jordan elimination: Row reduction continued until reduced row-echelon form is obtained.

LU factorization: Representation of a matrix through lower- and upper-triangular factors.

Nullity: The dimension of a matrix's null space.

Null space: The set of vectors mapped to the zero vector by a matrix.

Partial pivoting: Selecting the largest-magnitude available entry in the pivot column.

Pivot: A leading nonzero entry used during elimination.

Rank: The number of pivots, equal to the dimension of the row or column space.

Reduced row-echelon form: Echelon form in which each pivot is one and is alone in its column.

Residual: The vector $\mathbf{b}-A\widehat{\mathbf{x}}$ for a computed solution.

Row-echelon form: A staircase matrix form with zeros below each pivot.

Underdetermined system: A system with fewer independent constraints than unknowns.

Review Questions

Short Answer Questions

1. What three elementary row operations preserve the solution set of a linear system?

2. How does row-echelon form differ from reduced row-echelon form?

3. State the rank condition for consistency of $A\mathbf{x}=\mathbf{b}$.

4. What is the relationship between rank and nullity?

5. Why is partial pivoting used in floating-point Gaussian elimination?

Long Answer Questions

1. Explain Gaussian elimination and Gauss–Jordan elimination, comparing their procedures and computational purposes.

2. Derive the classification of a linear system into unique, infinite, and inconsistent cases using matrix ranks.

3. Explain how LU factorization is obtained from elimination and why it is useful for multiple right-hand sides.

4. Discuss the distinction between numerical conditioning and algorithmic stability, including the roles of residual and backward error.

5. Explain how underdetermined systems are parameterized and why their consistent solution sets are affine spaces.

Analytical / Application Questions

1. Solve by Gaussian elimination: $2x+y-z=3$, $x-y+2z=2$, and $3x+2y+z=7$. Verify the result in all three equations.

2. For a $4\times7$ matrix of rank three, determine its nullity and describe the number of free variables in any consistent system $A\mathbf{x}=\mathbf{b}$.

3. An AI pipeline repeatedly solves systems with the same dense coefficient matrix and 1,000 different right-hand sides. Compare refactoring the matrix for every system with computing one pivoted LU factorization and reusing it. Discuss computational cost and numerical checks.

Multiple Choice Questions with Answers

Question 1

For a system A x = b to be consistent, which rank condition must hold?

A. rank(A) = rank([A | b])

B. rank(A) must equal the number of rows

C. rank(A) must be zero

D. rank([A | b]) must exceed rank(A)

Correct Answer: A. rank(A) = rank([A | b])

Question 2

A matrix has 5 columns and rank 3. What is its nullity?

A. 1

B. 2

C. 3

D. 5

Correct Answer: B. 2

Question 3

What is the solution of x + y = 3 and x - y = 1?

A. x = 1, y = 2

B. x = 2, y = 1

C. x = 3, y = 0

D. x = 1, y = 1

Correct Answer: B. x = 2, y = 1

Question 4

Which operation is not an elementary row operation?

A. Interchanging two rows

B. Multiplying a row by a nonzero scalar

C. Adding a multiple of one row to another

D. Multiplying two rows componentwise

Correct Answer: D. Multiplying two rows componentwise

Question 5

A consistent system has 6 unknowns and coefficient-matrix rank 4. How many free variables does it have?

A. 1

B. 2

C. 4

D. 6

Correct Answer: B. 2

Question 6

Using the conventional flop count, what is the leading computational cost of dense n by n Gaussian or LU factorization?

A. Approximately n squared flops

B. Approximately one-third n cubed flops

C. Approximately two-thirds n cubed flops

D. Approximately 2n cubed flops

Correct Answer: C. Approximately two-thirds n cubed flops

Question 7

What does partial pivoting select at each elimination stage?

A. The smallest entry in the entire matrix

B. The largest-magnitude available entry in the current pivot column

C. The largest right-hand-side value

D. The first zero entry in the current row

Correct Answer: B. The largest-magnitude available entry in the current pivot column

Question 8

If x_p is a particular solution of a consistent system A x = b, what is its general solution?

A. x = x_p + z, where z belongs to the null space of A

B. x = A x_p

C. x = x_p + b

D. x = z for every vector z

Correct Answer: A. x = x_p + z, where z belongs to the null space of A

Question 9

For a computed solution x-hat, which expression is the residual?

A. A x-hat + b

B. x-hat - b

C. b - A x-hat

D. A inverse minus x-hat

Correct Answer: C. b - A x-hat

Question 10

In the calibrated sensor case, what is the recovered value of x?

A. -0.2

B. 0.2

C. 3.4

D. 4.0

Correct Answer: B. 0.2