Numerical Linear Algebra

The Central Idea

A basic linear system has the form

where \(A\) is a matrix, \(\mathbf{x}\) is the unknown vector, and \(\mathbf{b}\) is the right-hand side.

In exact mathematics, the question is whether a solution exists and how it can be characterized. In numerical computation, we must also ask:

• How accurately can the solution be computed?
• How sensitive is it to perturbations?
• Which algorithm is stable?
• How much computation and storage are required?

Why Linear Algebra Is Central to Numerical Computing

Many nonlinear or continuous problems become linear algebra problems after approximation.

Matrix systems arise in:

• finite-difference and finite-element methods;
• least-squares fitting;
• optimization;
• differential equations;
• network models;
• statistics and machine learning;
• eigenvalue and stability analysis;
• discretized physical fields.

Exact Linear Algebra and Numerical Linear Algebra

Exact linear algebra studies identities such as

Numerical linear algebra studies what happens when matrix entries and arithmetic are represented approximately.

A computed result may therefore satisfy

rather than reproduce the exact identity bit for bit.

Solving a \(2\times2\) System

Consider

In matrix form,

The numerical methods in this lesson generalize this familiar process to much larger systems.

Direct and Iterative Methods

Numerical methods for linear systems fall broadly into two families.

1. Direct methods aim to obtain the solution after a finite sequence of algebraic operations.
2. Iterative methods produce a sequence of approximations that ideally converges to the solution.

Direct Methods

Important direct methods include:

• Gaussian elimination;
• LU factorization;
• Cholesky factorization;
• QR factorization.

In finite precision, even a direct method produces an approximation rather than an exact symbolic result.

Iterative Methods

Iterative methods begin with an approximation

and construct

with the goal

They become especially important for large sparse systems.

Gaussian Elimination

Gaussian elimination systematically removes variables from a linear system.

The matrix is transformed into upper-triangular form:

where \(U\) has zeros below the diagonal.

The unknowns are then recovered by back substitution.

Elimination as Row Operations

The legal elementary row operations are:

1. interchange two rows;
2. multiply a row by a nonzero scalar;
3. add a multiple of one row to another.

In numerical algorithms, these operations must also be evaluated for stability and rounding behavior.

A First Elimination Step

Consider

To eliminate the \(4\) below the first pivot, use the multiplier

Replace the second row by

The result is

Pivots

A pivot is the matrix entry used to eliminate entries beneath it.

Dividing by a pivot that is zero is impossible.

Dividing by a pivot that is extremely small can also create large multipliers and numerical instability.

Why Small Pivots Are Dangerous

Suppose the pivot is

Dividing an ordinary-sized entry by this pivot can produce a multiplier on the order of

Large intermediate values can magnify rounding error and produce severe cancellation.

Pivoting

Pivoting reorders equations so that a more suitable pivot is used.

The most common strategy is partial pivoting, which selects a large-magnitude entry from the current column.

Partial Pivoting

At elimination step \(k\), partial pivoting selects a row \(p\) satisfying

Rows \(p\) and \(k\) are then exchanged before elimination continues.

This helps control the size of elimination multipliers.

Complete Pivoting

Complete pivoting searches the remaining submatrix for a large pivot and may interchange both rows and columns.

It can provide additional protection but requires more searching and bookkeeping.

Partial pivoting is the standard practical compromise for many dense systems.

LU Factorization

Gaussian elimination can often be organized as a matrix factorization

where \(L\) is lower triangular and \(U\) is upper triangular.

With row pivoting, the relationship is commonly written

where \(P\) is a permutation matrix.

Why LU Factorization Is Useful

Suppose many systems share the same matrix:

Factor \(A\) once:

Each new right-hand side can then be solved using inexpensive triangular solves.

Forward Substitution

For a lower-triangular system

the first equation determines \(y_1\), the second determines \(y_2\), and so on.

This sequential process is called forward substitution.

Back Substitution

For an upper-triangular system

the final equation determines the last unknown first.

The solution then proceeds upward through the system.

This is back substitution.

Do Not Form the Inverse Just to Solve a System

Mathematically,

But explicitly computing \(A^{-1}\) is usually unnecessary when the goal is only to solve

Factorization and triangular solves are generally more efficient and numerically preferable.

Cholesky Factorization

If \(A\) is symmetric positive definite, it has a special factorization

where \(L\) is lower triangular.

Cholesky factorization exploits matrix structure and can be more efficient than a general LU factorization.

Positive Definite Matrices

A real symmetric matrix \(A\) is positive definite if

for every nonzero vector \(\mathbf{x}\).

Positive-definite matrices arise naturally in optimization, energy models, least-squares systems, and differential-equation discretizations.

QR Factorization

A matrix \(A\) may be factored as

where \(Q\) has orthonormal columns and \(R\) is upper triangular.

QR factorization is especially important in least squares and eigenvalue computation.

Orthogonal Matrices Are Numerically Valuable

For an orthogonal matrix \(Q\),

Orthogonal transformations preserve Euclidean length:

This norm-preserving behavior is one reason orthogonal methods are important in stable numerical algorithms.

Householder Transformations

A Householder transformation has the form

It reflects vectors across a hyperplane.

Householder transformations can systematically introduce zeros and are widely used to compute QR factorizations stably.

Givens Rotations

A Givens rotation acts in a selected coordinate plane to eliminate one matrix entry.

These rotations are particularly useful when matrices are sparse or when a factorization must be updated incrementally.

Least Squares

An overdetermined system

may have more equations than unknowns and no exact solution.

Least squares chooses \(\mathbf{x}\) to minimize

The Residual Vector

For an approximate solution \(\widehat{\mathbf{x}}\), define

The residual measures how well the computed vector satisfies the original equations.

Residual Is Not the Same as Solution Error

The solution error is

The residual is

A small residual does not necessarily imply a small solution error when the matrix is ill-conditioned.

Matrix Conditioning

For a nonsingular matrix, a standard condition number is

The exact numerical value depends on the norm.

A large condition number signals sensitivity of the solution to perturbations.

Condition Number in the \(2\)-Norm

In the Euclidean operator norm,

where \(\sigma_{\max}\) and \(\sigma_{\min}\) are the largest and smallest singular values.

A very small smallest singular value indicates near singularity.

Near-Singular Systems

If a matrix is close to singular, small perturbations in its entries or in \(\mathbf{b}\) can produce large changes in the solution.

This is an issue of conditioning, not necessarily an algorithmic failure.

Perturbation in the Right-Hand Side

Suppose

and the right-hand side changes to

The solution changes to

Condition numbers provide bounds relating relative changes in \(\mathbf{b}\) to relative changes in \(\mathbf{x}\).

Backward Error for Linear Systems

A computed solution \(\widehat{\mathbf{x}}\) can be interpreted through a nearby problem:

If the required perturbations are small, the algorithm has small backward error.

Gaussian Elimination with Pivoting and Stability

Gaussian elimination with partial pivoting is widely used because it is generally effective at controlling numerical growth.

Its practical success comes from keeping elimination multipliers moderate and avoiding unnecessary division by small pivots.

Growth Factor

During elimination, intermediate entries may become larger than the original matrix entries.

The growth factor measures this amplification.

Large growth can increase the effect of floating-point rounding.

Scaling a Linear System

Rows or columns may differ greatly in magnitude.

Scaling can bring quantities to more comparable numerical ranges before factorization or iteration.

Proper scaling can improve numerical behavior and make tolerances easier to interpret.

Sparse Matrices

A sparse matrix contains mostly zeros.

Large systems from differential equations and networks are often sparse.

Sparse numerical linear algebra stores and computes primarily with the nonzero entries rather than treating the matrix as dense.

Fill-In

During Gaussian elimination, entries that were originally zero can become nonzero.

This phenomenon is called fill-in.

Sparse direct solvers often reorder rows and columns to reduce fill-in and memory usage.

Jacobi Iteration

Split the matrix as

where \(D\) is diagonal and \(L\) and \(U\) contain the lower and upper off-diagonal parts.

Jacobi iteration can be written

Component Form of Jacobi Iteration

For the system

Jacobi updates each component using values from the previous iteration:

Gauss-Seidel Iteration

Gauss-Seidel improves on the Jacobi idea by using newly updated components as soon as they become available.

Its matrix form is

Why Iterative Methods Can Be Attractive

For very large sparse systems, direct factorization may require substantial memory and computation.

Iterative methods can exploit:

• sparsity;
• matrix-vector products;
• good initial guesses;
• problem-specific structure;
• parallel computation.

Convergence of Stationary Iterations

Many linear iterations have the form

The error satisfies

Convergence is therefore controlled by powers of \(T\).

Spectral Radius

The spectral radius of \(T\) is

For a stationary linear iteration, the condition

is the fundamental convergence criterion.

Diagonal Dominance

A matrix is strictly diagonally dominant by rows when

for every row.

Diagonal dominance provides useful sufficient conditions for convergence of common stationary iterative methods.

Conjugate Gradient

For large symmetric positive-definite systems, the conjugate-gradient method is a major iterative solver.

Rather than correcting one coordinate at a time, it builds search directions adapted to the geometry induced by \(A\).

It can converge efficiently without forming a matrix factorization.

Krylov Subspaces

Many modern iterative methods build approximations in a Krylov subspace:

Conjugate gradient and GMRES are important examples.

Preconditioning

A preconditioner transforms a linear system into one that is easier for an iterative method to solve.

Ideally, choose \(M\) so that

has more favorable spectral or conditioning properties while systems involving \(M\) remain inexpensive to solve.

The Goal of Preconditioning

A good preconditioner balances two competing goals:

1. make the transformed system easier to solve;
2. keep application of the preconditioner inexpensive.

A perfect preconditioner that is as expensive as solving the original problem offers little practical advantage.

Stopping Criteria for Iterative Solvers

A common stopping test uses the residual:

A relative criterion may require

Residual Reduction Is Not the Whole Story

An ill-conditioned matrix can convert a small residual into a larger solution error.

Therefore, iterative stopping criteria should be interpreted together with matrix conditioning and application requirements.

Eigenvalue Problems

An eigenvalue problem seeks

where \(\mathbf{v}\neq\mathbf{0}\).

Numerical eigenvalue methods are essential in stability analysis, vibration, quantum models, principal-component methods, and dynamical systems.

The Power Method

Starting from a nonzero vector \(\mathbf{x}_0\), repeatedly compute

and normalize:

Under suitable conditions, the direction approaches an eigenvector associated with the dominant eigenvalue.

Dominant Eigenvalue

A dominant eigenvalue has magnitude greater than the magnitudes of the competing eigenvalues:

The ratio

influences the convergence rate of the power method.

The Rayleigh Quotient

For a nonzero vector \(\mathbf{x}\), define

When \(\mathbf{x}\) approximates an eigenvector, the Rayleigh quotient provides an associated eigenvalue estimate.

The QR Algorithm

QR factorization also supports a fundamental eigenvalue iteration.

Starting with \(A_0=A\), factor

and form

With appropriate refinements, this process forms the basis of modern dense eigenvalue computation.

Singular Value Decomposition

Any real \(m\times n\) matrix has a singular value decomposition

The diagonal entries of \(\Sigma\) are the singular values.

The SVD reveals rank, conditioning, dominant directions, and low-dimensional structure.

Numerical Rank

In exact arithmetic, rank depends on whether singular values are exactly zero.

In finite precision, a tiny singular value may be indistinguishable from zero relative to the scale and accuracy of the problem.

This leads to the idea of numerical rank.

Least Squares and Normal Equations

The least-squares problem

leads formally to the normal equations

However, forming \(A^TA\) can worsen conditioning.

Why QR Is Often Preferred for Least Squares

If

then

Orthogonal transformations preserve Euclidean norms and avoid explicitly squaring the condition number through \(A^TA\).

Conditioning of the Normal Equations

In the \(2\)-norm,

for full-column-rank \(A\).

A moderately ill-conditioned matrix can therefore become much more difficult when converted to normal equations.

Backward Stability and Matrix Algorithms

A numerically stable matrix algorithm often produces the exact result for a nearby matrix problem.

For example, a computed factorization may satisfy

for a small perturbation \(\Delta A\).

This is the matrix analogue of backward stability developed earlier.

Floating-Point Effects in Dot Products

A dot product is

Each multiplication and addition can introduce rounding.

Since dot products appear throughout matrix multiplication and factorizations, their error behavior is fundamental to numerical linear algebra.

Matrix Multiplication and Accumulated Error

Each entry of

is a dot product:

Large matrix computations therefore involve many opportunities for finite- precision error to accumulate.

Computational Complexity

Numerical algorithms must be judged not only by accuracy but also by cost.

Dense Gaussian elimination for an \(n\times n\) system requires on the order of

arithmetic operations.

Triangular solves require on the order of \(n^2\) operations.

Why Factorization Reuse Matters

If one matrix \(A\) is used with many right-hand sides, paying the \(O(n^3)\) factorization cost once can be efficient.

Each later solve then costs roughly

Algorithm design should therefore consider the full workload, not only a single solve.

Memory Complexity Matters Too

A dense \(n\times n\) matrix contains

entries.

For very large \(n\), memory can become the dominant constraint even before arithmetic cost does.

Sparse storage and matrix-free methods address this issue.

Matrix-Free Computation

Some iterative methods require only the ability to compute

for arbitrary vectors \(\mathbf{x}\).

The full matrix need not always be stored explicitly.

This can be valuable for extremely large discretized systems.

Verification by Residual

After computing \(\widehat{\mathbf{x}}\), evaluate

A small residual is an essential consistency check.

It should be interpreted relative to the scales of \(A\), \(\widehat{\mathbf{x}}\), and \(\mathbf{b}\).

A Scaled Residual

One useful dimensionless diagnostic is

This places the raw residual in context and is closely related to backward-error reasoning.

Iterative Refinement

After obtaining \(\widehat{\mathbf{x}}\), compute the residual

Solve the correction equation

and update

This process is called iterative refinement.

Why Iterative Refinement Can Help

A computed factorization may already contain rounding error, yet residual information can reveal how the approximate solution fails to satisfy the original equations.

Correcting that residual can improve the solution, especially when residuals are evaluated accurately enough.

Choosing a Solver by Matrix Structure

Matrix structure should guide method selection.

• general dense matrix: LU with pivoting;
• symmetric positive definite matrix: Cholesky or conjugate gradient;
• least-squares problem: QR or SVD;
• large sparse system: sparse direct or iterative Krylov methods;
• rank-deficient or highly ill-conditioned problem: SVD may be especially informative.

A Reliable Direct-Solve Workflow

1. Inspect matrix dimensions and structure.
2. Check scaling and obvious singularity issues.
3. Select a structure-appropriate factorization.
4. Use pivoting when required.
5. Solve the resulting triangular systems.
6. Compute a residual.
7. Estimate conditioning when accuracy matters.
8. Apply iterative refinement if useful.

A Reliable Iterative-Solve Workflow

1. Identify matrix structure and sparsity.
2. Select an appropriate iterative method.
3. Choose an initial guess.
4. Apply a useful preconditioner when available.
5. Monitor residual reduction.
6. Watch for stagnation or divergence.
7. Stop using absolute and relative criteria appropriate to the problem.
8. Verify the final residual and interpret it through conditioning.

A Reliable Least-Squares Workflow

1. Form the design matrix \(A\) and data vector \(\mathbf{b}\).
2. Inspect scaling and rank.
3. Avoid normal equations when conditioning is a concern.
4. Use QR for a stable general solve.
5. Use SVD when rank deficiency or sensitivity must be diagnosed.
6. Compute and inspect the residual.
7. Interpret parameter uncertainty in the context of conditioning.

Common Error: Forming \(A^{-1}\) Unnecessarily

The formula

is mathematically correct but does not mean the inverse should be explicitly formed in numerical software.

Factorization-based solves are normally the better computational strategy.

Common Error: Eliminating Without Pivoting

A tiny pivot can create enormous multipliers and severe rounding amplification.

Row exchanges are not merely cosmetic; they can be essential to numerical stability.

Common Error: Treating a Small Residual as Proof of Accuracy

A small residual shows that the computed vector nearly satisfies the equations.

An ill-conditioned matrix can still produce a large solution error.

Residual and condition number should be interpreted together.

Common Error: Ignoring Matrix Structure

Applying a generic dense solver to every problem can waste computation and sacrifice numerical advantages.

Symmetry, positive definiteness, sparsity, band structure, and rank information can all guide a better method.

Common Error: Using Normal Equations Automatically

Normal equations are mathematically valid for least squares, but

can make the numerical problem substantially more sensitive.

QR or SVD is often preferable when accuracy matters.

Common Error: Ignoring Scaling

Large disparities in row or column magnitudes can make numerical behavior harder to control.

Rescaling can improve algorithmic performance without changing the underlying solution in an essential way.

Common Error: Iterating Without a Convergence Test

An iterative method should not be assumed to converge simply because it produces a sequence of vectors.

Analyze the iteration matrix, spectral properties, residual trend, and structural convergence conditions.

Common Error: Stopping at a Fixed Iteration Count

A fixed iteration budget may stop too early on a difficult system or waste work on an easy one.

Error- or residual-based stopping criteria are usually more informative.

Example: One Jacobi Step

Consider

Rewrite:

Starting from

Jacobi gives

and

Example: One Gauss-Seidel Step

Using the same system and starting point, compute

Gauss-Seidel immediately uses this new value:

The distinction is that Gauss-Seidel uses newly available information within the same iteration.

Example: Matrix Condition Number

Suppose

and

Then

This indicates substantial sensitivity compared with a condition number near \(1\).

Example: Residual Check

Suppose

Then

If

the residual is

Practice 1: Matrix Form

Write the system

in the form

Practice 2: Elimination Multiplier

In

what multiplier eliminates the \(15\) below the first pivot?

Practice 3: Why Pivot?

Explain why a row exchange may be desirable when the current pivot has very small magnitude.

Practice 4: Condition Number

If

and

compute \(\kappa(A)\).

Practice 5: Residual

Let

Compute

Practice 6: Choose the Factorization

Which factorization is especially natural for a symmetric positive-definite matrix?

Practice 7: Least Squares

What quantity is minimized in the standard linear least-squares problem?

Practice 8: Jacobi Concept

In Jacobi iteration, are newly computed components used immediately within the same iteration or only on the next iteration?

Practice 9: Spectral Radius

A stationary iterative method has iteration matrix \(T\) with

What does this suggest about convergence?

Practice 10: Normal Equations

If

what is

for a full-column-rank matrix?

Practice Answers

Practice 1

Practice 2

The multiplier is

Practice 3

Dividing by a tiny pivot can create very large elimination multipliers, which magnify rounding error and produce large intermediate values. A row exchange can provide a safer pivot.

Practice 4

Practice 5

First,

Therefore,

Practice 6

The natural factorization is

the Cholesky factorization.

Practice 7

Standard least squares minimizes

Practice 8

Jacobi uses values from the previous iteration. Newly computed components are not used until the next iteration.

Practice 9

Since

the stationary iteration satisfies the fundamental spectral-radius condition for convergence.

Practice 10

What to Remember

Numerical linear algebra studies how matrix computations behave under finite precision.

• linear systems are written \(A\mathbf{x}=\mathbf{b}\);
• Gaussian elimination transforms systems to triangular form;
• pivoting protects elimination from dangerous small pivots;
• LU factorization organizes elimination for efficient repeated solves;
• Cholesky exploits symmetric positive-definite structure;
• QR uses orthogonal transformations and is central to least squares;
• condition numbers measure sensitivity to perturbations;
• residual size and solution error are different quantities;
• sparse and iterative methods become essential for very large systems;
• spectral radius controls convergence of stationary iterations;
• preconditioning can accelerate iterative convergence;
• eigenvalue methods and the SVD reveal important matrix structure;
• stable algorithms, appropriate scaling, and structure-aware methods are fundamental to reliable matrix computation.

Where We Go Next

Numerical linear algebra provides machinery for discrete matrix systems. The next lesson turns to another central approximation task: constructing functions that pass through or closely represent known data.

The next lesson, Interpolation, develops polynomial interpolation, Lagrange and Newton forms, divided differences, interpolation error, piecewise methods, and the numerical issues that arise when discrete data are used to construct continuous approximations.