5. Linear Algebra and Matrix Computation

Introduction and Context

Linear algebra provides the language and computational machinery for representing structure, relationships, and transformations in computer science. A digital image can be modeled as an array of intensity values; a three-dimensional object as a collection of vectors and transformation matrices; a dataset as a matrix whose rows represent observations and columns represent features; and the parameters of a machine-learning model as vectors or higher-order tensors. In each case, linear algebra makes it possible to express complex operations compactly and to implement them efficiently.

The subject is not limited to solving simultaneous equations. It addresses vector spaces, linear independence, bases, linear mappings, eigenvalues, orthogonality, matrix factorizations, numerical stability, and tensor operations. These ideas support computer graphics, information retrieval, signal processing, scientific simulation, optimization, data analytics, and artificial intelligence. Matrix computation translates the abstract theory into algorithms that operate on finite-precision computers.

A central theme of this chapter is the relationship between mathematical structure and computational behavior. An exact algebraic method may be unsuitable when a matrix is very large, sparse, or ill-conditioned. Conversely, a numerical approximation may be highly useful when its error is controlled and its computational cost is acceptable. Graduate-level computer science therefore requires both theoretical understanding and judgment about representation, algorithms, accuracy, and scalability.

Learning Objectives

After studying this chapter, students should be able to:

1. Explain vector spaces, subspaces, bases, dimension, rank, and linear independence.

2. Represent and analyze linear transformations using matrices.

3. Solve linear systems and interpret consistency, rank, and null spaces.

4. Explain eigenvalues, eigenvectors, diagonalization, and spectral interpretations.

5. Apply orthogonality, least squares, QR factorization, singular value decomposition, and related decompositions.

6. Distinguish direct and iterative numerical methods and assess their stability and computational cost.

7. Explain the role of tensors and matrix computation in graphics, analytics, machine learning, and scientific computing.

8. Select appropriate representations and algorithms for practical computational problems.

Vector Spaces: Structure Before Coordinates

A vector space over a field such as the real numbers \(\mathbb{R}\) or complex numbers \(\mathbb{C}\) is a set whose elements can be added and multiplied by scalars while satisfying closure, associativity, distributivity, and identity properties. Although vectors in \(\mathbb{R}^n\) are familiar ordered lists, vectors may also be polynomials, functions, signals, matrices, or probability distributions.

A subset \(W\) of a vector space \(V\) is a subspace if it contains the zero vector and is closed under linear combinations. Thus, if \(u,v\in W\) and \(\alpha,\beta\) are scalars, then \(\alpha u+\beta v\in W\). The set of all linear combinations of vectors \(v_1,\ldots,v_k\) is their span. A vector belongs to the span precisely when it can be constructed from those vectors.

Vectors are linearly independent when the equation

\[c_1v_1+c_2v_2+\cdots+c_kv_k=0\]

has only the solution \(c_1=\cdots=c_k=0\). Dependence means that at least one vector is redundant. A basis is a linearly independent spanning set. The number of vectors in every basis of a finite-dimensional vector space is its dimension. Coordinates are not intrinsic to a vector; they depend on the selected basis. This observation is important in computer graphics and numerical analysis, where changing coordinates can simplify computation.

For a matrix \(A\in\mathbb{R}^{m\times n}\), four fundamental subspaces are especially significant: the column space, row space, null space \(\{x:Ax=0\}\), and left null space \(\{y:A^Ty=0\}\). The rank of \(A\) is the dimension of its column space, equivalently its row space. The rank–nullity theorem states

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

This theorem connects the number of independent constraints with the number of remaining degrees of freedom.

Linear Transformations and Matrix Representation

A function \(T:V\rightarrow W\) is linear if \(T(\alpha u+\beta v)=\alpha T(u)+\beta T(v)\). Its kernel, or null space, consists of vectors mapped to zero, while its image or range consists of all attainable outputs. Once bases are chosen for the domain and codomain, every linear transformation can be represented by a matrix. Matrix multiplication then represents composition: if \(S\) follows \(T\), the composite transformation is represented by \([S][T]\).

In two-dimensional graphics, a rotation through angle \(\theta\) is represented by

\[R=\begin{bmatrix}\cos\theta&-\sin\theta\\\sin\theta&\cos\theta\end{bmatrix}.\]

A vector is rotated by computing \(Rv\). Scaling, reflection, shear, and projection have similarly simple matrix representations. Homogeneous coordinates extend this framework so that translations can be represented by matrix multiplication in three dimensions. This permits a graphics pipeline to combine modeling, viewing, and projection transformations into a single composite matrix.

Matrix multiplication is associative but generally not commutative. The order of transformations therefore matters. A transformation that scales and then rotates an object need not produce the same result as one that rotates and then scales it. This is both a mathematical fact and a practical source of errors in graphics software.

Linear Systems and Their Interpretation

A linear system is written as \(Ax=b\), where \(A\) is the coefficient matrix, \(x\) is the unknown vector, and \(b\) is the observation or output vector. A system may have no solution, one solution, or infinitely many solutions. It is consistent exactly when \(b\) belongs to the column space of \(A\). If \(A\) is square and nonsingular, the solution is unique and \(x=A^{-1}b\), although explicitly computing the inverse is usually less efficient and less stable than solving the system through factorization.

Gaussian elimination transforms a system into upper triangular form using elementary row operations. Back substitution then obtains the unknowns. Pivoting is commonly used to avoid division by very small elements and to reduce numerical error. For multiple right-hand sides, factorizing \(A\) once and solving repeatedly is more efficient than repeating elimination.

The determinant is a scalar associated with a square matrix. It is zero exactly when the matrix is singular. Geometrically, its absolute value gives the volume scaling factor of the transformation, while its sign indicates orientation reversal. Although determinants are conceptually useful, determinant-based formulas are generally poor numerical methods for large systems.

Inner Products, Norms, and Orthogonality

An inner product generalizes the dot product. In \(\mathbb{R}^n\), \(\langle x,y\rangle=x^Ty\). It induces the Euclidean norm \(\|x\|_2=\sqrt{x^Tx}\), measuring length. Other norms, such as the Manhattan norm \(\|x\|_1\) and maximum norm \(\|x\|_\infty\), are useful for different modeling and error criteria.

Two vectors are orthogonal when their inner product is zero. Orthogonal sets of nonzero vectors are linearly independent, and orthonormal sets additionally have unit length. Orthogonal transformations preserve lengths and angles, making them numerically attractive. The projection of a vector \(b\) onto a nonzero vector \(a\) is

\[\operatorname{proj}_a(b)=\frac{a^Tb}{a^Ta}a.\]

The Gram–Schmidt process constructs an orthonormal basis from independent vectors. In practical numerical software, modified Gram–Schmidt or Householder transformations are often preferred because the classical procedure can lose orthogonality through round-off error.

Least Squares and Approximation

When \(Ax=b\) is inconsistent, an exact solution does not exist. A least-squares solution minimizes the residual norm:

\[\hat{x}=\arg\min_x\|Ax-b\|_2^2.\]

The normal equations are \(A^TA\hat{x}=A^Tb\). Geometrically, \(A\hat{x}\) is the orthogonal projection of \(b\) onto the column space of \(A\). In linear regression, the columns of \(A\) contain feature values, \(b\) contains observed responses, and \(\hat{x}\) contains fitted coefficients.

Although normal equations are simple, forming \(A^TA\) squares the condition number and may amplify numerical error. QR factorization, \(A=QR\), provides a more stable approach because the least-squares problem becomes \(R\hat{x}=Q^Tb\). Singular value decomposition is even more informative when the matrix is rank-deficient or nearly so.

Eigenvalues and Eigenvectors

For a square matrix \(A\), a nonzero vector \(v\) is an eigenvector if \(Av=\lambda v\), where \(\lambda\) is its eigenvalue. The transformation changes the magnitude, and possibly direction, of an eigenvector without rotating it into a different independent direction. Eigenvalues satisfy the characteristic equation

\[\det(A-\lambda I)=0.\]

If a matrix has enough linearly independent eigenvectors, it can be diagonalized as \(A=PDP^{-1}\), where the columns of \(P\) are eigenvectors and \(D\) is diagonal. Then \(A^k=PD^kP^{-1}\), simplifying repeated transformations. Diagonalization is not always possible; defective matrices lack a complete eigenvector basis.

For symmetric real matrices, the spectral theorem guarantees an orthonormal eigenbasis: \(A=Q\Lambda Q^T\). Such matrices have real eigenvalues and particularly well-behaved spectral interpretations. Eigenvalues can describe stability in dynamical systems, principal directions in data, vibration modes in engineering, and dominant structures in networks. For a discrete system \(x_{t+1}=Ax_t\), eigenvalues with magnitude less than one generally correspond to decaying modes, while eigenvalues with magnitude greater than one indicate growth.

Matrix Decompositions as Computational Models

A decomposition expresses a matrix as a product of simpler matrices. The choice of decomposition should reflect the matrix structure and the intended task.

LU factorization. A matrix is factored as \(PA=LU\), where \(P\) represents row permutations, \(L\) is lower triangular, and \(U\) is upper triangular. LU is effective for solving many systems with the same coefficient matrix.

Cholesky factorization. If \(A\) is symmetric positive definite, then \(A=LL^T\). Cholesky is approximately twice as efficient as general LU and is widely used in covariance computations and optimization.

QR factorization. The decomposition \(A=QR\), with orthogonal \(Q\) and upper triangular \(R\), supports stable least-squares computation and eigenvalue algorithms.

Singular value decomposition. Every real matrix can be written

\[A=U\Sigma V^T,\]

where \(U\) and \(V\) are orthogonal and \(\Sigma\) contains nonnegative singular values in descending order. The right singular vectors identify important directions in the input space; the singular values indicate their scaling; and the left singular vectors describe corresponding output directions. The rank equals the number of nonzero singular values. A truncated SVD provides a best low-rank approximation in the spectral and Frobenius norms, making it useful for compression, denoising, and recommender systems.

Principal component analysis (PCA) can be performed through the eigendecomposition of a centered covariance matrix or through the SVD of a centered data matrix. The first principal component captures the greatest variance subject to unit-length constraints, and subsequent components capture maximum remaining variance while maintaining orthogonality. PCA is not automatically a causal model; it is a representation technique whose results depend on scaling, centering, and the data distribution.

Numerical Matrix Computation

Exact arithmetic is rare in computer systems. Floating-point numbers have finite precision, so algorithms must be evaluated in terms of conditioning and stability. The condition number of an invertible matrix in a selected norm is \(\kappa(A)=\|A\|\|A^{-1}\|\). A large condition number means that small perturbations in data or round-off may produce large changes in the solution. An ill-conditioned problem is intrinsically sensitive; a numerically unstable algorithm can introduce additional error even when the problem itself is well-conditioned.

A stable algorithm produces a result that is the exact solution to a nearby problem, or whose error is controlled relative to the conditioning of the original problem. Residuals are often more informative than raw differences: after computing \(\hat{x}\), evaluate \(r=b-A\hat{x}\). A small residual does not always imply a small solution error if \(A\) is ill-conditioned, but it is an essential diagnostic.

For large sparse systems, direct factorization may require excessive memory because of fill-in. Iterative methods generate a sequence of approximations. The Jacobi and Gauss–Seidel methods are simple stationary iterations; convergence depends on properties such as diagonal dominance. The conjugate-gradient method is designed for symmetric positive-definite systems. GMRES is suitable for more general nonsymmetric systems, though memory management and restart strategies may be required.

Preconditioning transforms a system into an equivalent one that is easier for an iterative method to solve. A good preconditioner reduces the effective condition number without being too expensive to apply. In scientific computing, sparse storage, parallel matrix-vector multiplication, domain decomposition, and hardware accelerators are central to scalable computation.

Tensors and Higher-Order Computation

A scalar is a zero-order tensor, a vector a first-order tensor, and a matrix a second-order tensor. A tensor of order three or higher can represent data indexed by several dimensions, such as video frames indexed by height, width, and time, or a neural-network weight array indexed by layer, input channel, and output channel.

Tensor operations include contraction, reshaping, permutation, outer products, and multilinear transformations. Matrix multiplication is itself a contraction over one shared index. In machine learning, tensor libraries execute these operations on CPUs, GPUs, and specialized accelerators. Broadcasting rules allow lower-dimensional arrays to participate in operations without explicitly copying data, but careless broadcasting can create large hidden memory costs.

Tensor decompositions generalize matrix factorizations. CP decomposition represents a tensor as a sum of rank-one tensors, while Tucker decomposition separates a smaller core tensor from factor matrices. These methods support compression, feature extraction, and multiway data analysis, but tensor rank and decomposition uniqueness are more complicated than their matrix counterparts.

Applications Across Computing

In computer graphics, vectors represent positions, directions, colors, and normals. Matrices implement transformations, camera projections, and coordinate changes. Homogeneous coordinates combine translation with linear transformations, while quaternions can provide efficient and stable rotation interpolation.

In data analytics, matrices hold observations and variables. Linear systems support regression, covariance estimation, and network analysis. SVD and PCA reduce dimensionality, while sparse matrices represent user–item interactions and graph connections.

In machine learning, matrix and tensor operations dominate forward propagation, gradient computation, attention mechanisms, and parameter updates. Automatic differentiation applies the chain rule through computational graphs. The quality of numerical linear algebra affects convergence, reproducibility, and model reliability. Normalization, regularization, and low-rank approximation often have direct linear-algebraic interpretations.

In scientific computing, discretized differential equations produce large sparse systems. Eigenvalue problems describe physical modes, while factorizations and iterative solvers enable simulation. In information retrieval, documents and queries can be represented as vectors, with inner products or cosine similarity measuring relevance. In cryptography and coding theory, finite-field linear algebra supports error-correcting codes and structured algebraic constructions.

Illustrative Case Study: Low-Rank Recommendation at Scale

Background. A digital platform wants to recommend products to users. Its interaction matrix \(R\) has one row per user and one column per product. Most entries are missing because each user interacts with only a small fraction of products.

Situation and problem. Directly storing and factorizing the dense matrix is wasteful. Moreover, a missing interaction does not necessarily mean dislike. The platform needs a compact representation that captures latent preferences while handling sparse observations.

Analysis. A low-rank model approximates observed ratings by \(R\approx UV^T\), where rows of \(U\) represent user factors and rows of \(V\) represent product factors. The rank controls the number of latent dimensions. SVD gives a useful conceptual model, but practical recommender systems often optimize factors using regularized matrix factorization over observed entries rather than computing a full dense SVD.

Approach. The system minimizes a regularized loss such as

\[\sum_{(i,j)\in\Omega}(R_{ij}-u_i^Tv_j)^2+\lambda(\|U\|_F^2+\|V\|_F^2),\]

where \(\Omega\) is the set of observed interactions and \(\lambda\) controls overfitting. Alternating optimization or stochastic gradient methods update user and product factors. Sparse data structures ensure that computation is proportional mainly to observed interactions. Evaluation uses a time-based holdout and metrics appropriate to ranking, not only prediction error.

Outcome. The factorized representation reduces storage and enables recommendations through inexpensive vector similarities. Performance depends on rank, regularization, data freshness, popularity bias, and the quality of the evaluation design.

Learning points. Low-rank structure can make a large problem computationally feasible. Missing data must be interpreted carefully. Numerical compression does not eliminate modeling bias, and an accurate algebraic approximation may still produce unfair or commercially undesirable recommendations.

Challenges, Risks, and Common Mistakes

A frequent mistake is to confuse an invertible matrix with a numerically safe matrix. A matrix may be theoretically nonsingular but so ill-conditioned that computed solutions are unreliable. Explicitly forming an inverse, using normal equations without considering conditioning, and ignoring pivoting are avoidable errors.

Other common problems include mixing row and column conventions, applying transformations in the wrong order, treating correlation as causation in PCA, and assuming that a low-rank approximation preserves every feature of the data. In iterative methods, failure to define a stopping criterion can waste computation or produce an inaccurate result. Stopping criteria should consider residual size, relative change, maximum iterations, and problem scale.

Large matrices create memory and communication challenges. Dense storage may be impossible, while sparse storage depends on the pattern of nonzero entries. Parallel implementations must address synchronization, load imbalance, and reproducibility. GPU acceleration is valuable for regular, highly parallel operations but may be inefficient for irregular sparse workloads or data transfers dominated by communication.

Professional and Ethical Considerations

Linear algebra systems influence decisions in credit, employment, health, education, and public services. A model based on matrix factorization may reproduce historical discrimination even when no protected attribute is explicitly included. Practitioners should examine data provenance, missingness, representation, fairness, privacy, and the consequences of errors.

Dimensionality reduction can obscure sensitive information, but it can also make auditing more difficult. Numerical reproducibility requires documenting software versions, precision, initialization, preprocessing, and hardware. In safety-critical scientific applications, residual checks, uncertainty analysis, validation against independent evidence, and human review are essential. Technical correctness is not sufficient when the modeled objective is inappropriate or the deployment context is harmful.

Emerging Directions and Best Practices

Modern linear algebra increasingly combines mathematical structure with hardware-aware computation. Randomized SVD and sketching methods approximate large matrix decompositions using fewer passes over data. Differentiable linear algebra integrates solvers and decompositions into machine-learning pipelines. Mixed-precision computation uses lower precision for speed while retaining higher precision where accuracy is sensitive. Quantum algorithms and tensor-network methods offer active research directions, although their practical advantages are problem-dependent.

Effective practice follows several principles:

• Choose representations that exploit sparsity, symmetry, block structure, or low rank.
• Use QR or SVD when stability matters, rather than relying automatically on normal equations.
• Inspect condition numbers, residuals, convergence behavior, and sensitivity.
• Scale features and variables when appropriate, while recording the transformation.
• Separate training, validation, and test data in analytics applications.
• Match solver choice to matrix properties such as symmetry, definiteness, and sparsity.
• Use libraries with tested kernels, but understand their assumptions and failure modes.
• Report approximation error, precision, stopping criteria, and reproducibility information.

Chapter Summary

Linear algebra models vectors, spaces, transformations, and relationships that occur throughout computer science. Bases and subspaces reveal structure; rank and nullity describe information and degrees of freedom; and matrices provide coordinate representations of linear transformations. Linear systems can be solved directly or approximately, while least squares handles inconsistent observations through projection.

Eigenvalues identify invariant directions and modes of behavior. LU, Cholesky, QR, and SVD decompositions convert difficult matrix operations into structured computations, with SVD offering particularly useful insight into rank and low-dimensional representation. Numerical computation requires attention to conditioning, stability, residuals, finite precision, sparsity, and convergence. Tensors extend these ideas to multidimensional data and underpin modern machine-learning systems.

The practical value of linear algebra depends on matching mathematical methods to computational and organizational contexts. Graphics, analytics, recommendation, simulation, and AI all benefit from efficient matrix computation, but ethical analysis, validation, and awareness of bias remain necessary whenever the resulting systems affect people or important decisions.

Key Terms

Vector space: A set closed under vector addition and scalar multiplication.

Subspace: A subset of a vector space that is itself a vector space under the inherited operations.

Span: The set of all linear combinations of a specified collection of vectors.

Linear independence: A property in which only the trivial linear combination equals the zero vector.

Basis: A linearly independent spanning set.

Rank: The dimension of a matrix’s column space, equivalently its row space.

Null space: The set of vectors mapped to zero by a matrix or linear transformation.

Linear transformation: A mapping that preserves vector addition and scalar multiplication.

Eigenvalue: A scalar satisfying \(Av=\lambda v\) for some nonzero vector \(v\).

Eigenvector: A nonzero vector whose direction is preserved by a linear transformation.

Orthogonality: The condition that two vectors have zero inner product.

Norm: A function measuring the size or length of a vector or matrix.

Least squares: An optimization method minimizing the squared residual norm.

Condition number: A measure of a problem’s sensitivity to perturbations.

LU factorization: A decomposition involving lower- and upper-triangular matrices.

QR factorization: A decomposition into an orthogonal matrix and an upper-triangular matrix.

Singular value decomposition: The factorization \(A=U\Sigma V^T\) using orthogonal factors and singular values.

Sparse matrix: A matrix containing mostly zero entries.

Tensor: A multidimensional generalization of scalars, vectors, and matrices.

Preconditioner: A transformation intended to improve iterative solver convergence.

Review Questions

Short Answer Questions

1. What is the difference between a span and a basis?

2. State the rank–nullity theorem and explain its significance.

3. Why is matrix multiplication generally order-dependent?

4. What geometric interpretation is associated with a least-squares solution?

5. What does a large condition number indicate?

Long Answer Questions

1. Explain the four fundamental subspaces of a matrix and describe their relationships.

2. Compare LU, Cholesky, QR, and SVD in terms of structure, purpose, and numerical behavior.

3. Derive the normal equations for least squares and explain why QR may be preferred computationally.

4. Explain eigenvalues, eigenvectors, diagonalization, and the spectral theorem for real symmetric matrices.

5. Discuss the role of tensors and matrix operations in modern machine-learning systems.

Analytical / Application Questions

1. A data matrix has highly correlated columns and a large condition number. Which computational and modeling choices would you make for regression, and why?

2. A graphics program produces different results when scaling and rotation operations are reordered. Use matrix composition to explain the behavior and propose a debugging strategy.

3. A recommendation matrix is extremely sparse and too large for dense SVD. Design a low-rank computational approach, including objective function, regularization, evaluation, and ethical checks.

22. Multiple Choice Questions with Answers

Question 1

Which statement correctly expresses the rank–nullity theorem for A in R^(m×n)?

A. rank(A) − nullity(A) = n

B. rank(A) + nullity(A) = n

C. rank(A) × nullity(A) = n

D. rank(A) + nullity(A) = m

Correct Answer: B. rank(A) + nullity(A) = n

Question 2

Why is QR factorization often preferred to normal equations for least-squares problems?

A. QR always requires forming A inverse

B. QR avoids the conditioning effects of forming A^T A

C. QR works only for square matrices

D. QR guarantees an exact solution to every inconsistent system

Correct Answer: B. QR avoids the conditioning effects of forming A^T A

Question 3

For a real symmetric matrix, which property is guaranteed by the spectral theorem?

A. It has an orthonormal eigenbasis and real eigenvalues

B. It is always singular

C. All eigenvalues are zero

D. It cannot be diagonalized

Correct Answer: A. It has an orthonormal eigenbasis and real eigenvalues

Question 4

In the decomposition A = UΣV^T, what do the singular values primarily describe?

A. The row labels of A

B. Scaling along paired input and output directions

C. The determinant only

D. The number of columns regardless of rank

Correct Answer: B. Scaling along paired input and output directions

Question 5

A system Ax = b is consistent precisely when:

A. b belongs to the column space of A

B. A has more rows than columns

C. det(A) is always equal to one

D. the null space contains only b

Correct Answer: A. b belongs to the column space of A

Question 6

What is the main purpose of a preconditioner in an iterative solver?

A. To increase the dimension of the system

B. To make an equivalent system easier to solve iteratively

C. To convert every matrix into a diagonal matrix exactly

D. To eliminate the need for a stopping criterion

Correct Answer: B. To make an equivalent system easier to solve iteratively

Question 7

In a discrete system x_(t+1) = Ax_t, eigenvalues with magnitude greater than one generally indicate:

A. A growing mode

B. A guaranteed zero solution

C. A preserved norm in every direction

D. A matrix with no eigenvectors

Correct Answer: A. A growing mode

Question 8

Which representation is most appropriate for a matrix containing mostly zeros?

A. Dense storage without exception

B. Sparse storage that records nonzero structure

C. A scalar representation

D. A determinant table

Correct Answer: B. Sparse storage that records nonzero structure

Question 9

In low-rank recommendation, regularization is primarily used to:

A. Increase overfitting to observed entries

B. Control model complexity and reduce overfitting

C. Replace all missing interactions with zeros as facts

D. Guarantee causal recommendations

Correct Answer: B. Control model complexity and reduce overfitting

Question 10

Why can explicitly computing A^−1 be a poor way to solve Ax=b?

A. It is often less efficient and less numerically stable than factorization-based solving

B. An inverse exists only for rectangular matrices

C. It cannot be used with floating-point numbers

D. It always changes the dimensions of b

Correct Answer: A. It is often less efficient and less numerically stable than factorization-based solving