Inner Product Spaces, Modified Gram-Schmidt & QR Matrix Factorization

In multi-dimensional data analysis, geometric concepts such as length, distance, angle, and perpendicularity provide the structural bedrock for optimization algorithms, signal processing filters, and regression estimation. An inner product space enriches an abstract vector space with a geometric metric, enabling rigorous projections and orthogonal factorizations.

1. Architectural Foundations of Inner Product Spaces#

In standard Euclidean space R^n, the familiar dot product u · v = sum(ui · vi) provides geometric measures of magnitude and angular separation. An axiomatic inner product space generalizes these mechanics to arbitrary vector spaces V over the real numbers R (or complex numbers C).

The Four Axioms of an Inner Product#

An inner product on a real vector space V is a function that associates with each pair of vectors u, v in V a real scalar <u, v> in R, satisfying four fundamental axioms for all vectors u, v, w in V and every scalar c in R:

1. Symmetry (Commutativity):

<u, v> = <v, u>

The ordering of the vector operands does not alter the resulting scalar value.

1. Linearity in the First Argument:

- Additivity: <u + v, w> = <u, w> + <v, w>

- Homogeneity: <c · u, v> = c · <u, v>

Together, these establish that the inner product is a linear functional in its first argument:

<a · u + b · v, w> = a · <u, w> + b · <v, w>.

1. Positive-Definiteness:

<v, v> >= 0 for all v in V, with <v, v> = 0 if and only if v = 0 (the zero vector).

The inner product of any non-zero vector with itself is strictly positive.

1. Induced Norm (Metric Length):

Every inner product naturally induces an absolute metric norm on V:

||v|| = sqrt(<v, v>).

The induced distance between two vectors u and v is given by d(u, v) = ||u - v|| = sqrt(<u - v, u - v>).

The Cauchy-Schwarz Inequality: The Cornerstone of Geometric Analysis#

Theorem: For all vectors u and v in an inner product space V:

|<u, v>| <= ||u|| · ||v||

The equality holds if and only if u and v are linearly dependent (one is a scalar multiple of the other).

Rigorous Analytical Derivation:

Consider an arbitrary real parameter t in R and evaluate the norm squared of the vector combination u + t · v. By positive-definiteness, this quantity is always non-negative:

p(t) = ||u + t · v||^2 = <u + t · v, u + t · v> >= 0

Expanding using linearity and symmetry:

p(t) = <u, u> + 2 · t · <u, v> + t^2 · <v, v>

p(t) = ||u||^2 + 2 · t · <u, v> + t^2 · ||v||^2 >= 0

This is a quadratic polynomial in t of the form a · t^2 + b · t + c >= 0, where:

• a = ||v||^2
• b = 2 · <u, v>
• c = ||u||^2

For a quadratic polynomial to remain non-negative for all real values of t, it can have at most one real root. Consequently, its discriminant Δ = b^2 - 4 · a · c must be non-positive:

b^2 - 4 · a · c <= 0

(2 · <u, v>)^2 - 4 · ||v||^2 · ||u||^2 <= 0

4 · <u, v>^2 <= 4 · ||u||^2 · ||v||^2

<u, v>^2 <= ||u||^2 · ||v||^2

Taking the positive square root on both sides proves the inequality:

|<u, v>| <= ||u|| · ||v||.

Generalized Angles and Function Spaces#

The Cauchy-Schwarz inequality guarantees that the ratio <u, v> / (||u|| · ||v||) always resides within the closed interval [-1, 1]. This justifies defining the generalized angle θ between non-zero vectors in arbitrary high-dimensional spaces:

cos(θ) = <u, v> / (||u|| · ||v||), where θ in [0, π].

Inner products extend beyond discrete finite tuples into continuous function spaces. In the Hilbert space L^2[a, b] of square-integrable functions, the inner product of two functions f(x) and g(x) is defined as:

<f, g> = integral from a to b of (f(x) · g(x)) dx.

Under this continuous inner product, the Legendre polynomials P_n(x) on [-1, 1] and trigonometric Fourier harmonics { 1, cos(n x), sin(n x) } on [-π, π] form mutually orthogonal sets of functions.

2. Orthogonality, Orthogonal Complements, and Subspace Geometry#

Orthogonality provides the ultimate geometric decoupling in multi-variable mathematics, eliminating cross-talk between coordinate components.

Orthogonal and Orthonormal Sets#

• Two vectors u and v in an inner product space V are orthogonal (denoted u ⊥ v) if and only if their inner product vanishes:

<u, v> = 0.

• A set of vectors S = { q1, q2, ..., qk } is an orthogonal set if every pair of distinct vectors is orthogonal:

<qi, qj> = 0 for all i ≠ j.

• If, in addition, every vector in the orthogonal set has unit length (||qi|| = 1 for all i), the set is termed orthonormal:

<qi, qj> = 1 if i = j, and 0 if i ≠ j (the Kronecker delta δ_ij).

The Generalized Pythagorean Theorem#

Theorem: If u and v are orthogonal vectors in an inner product space (u ⊥ v), then:

||u + v||^2 = ||u||^2 + ||v||^2

Proof:

||u + v||^2 = <u + v, u + v> = <u, u> + 2 · <u, v> + <v, v>.

Since u ⊥ v, <u, v> = 0, leaving ||u||^2 + ||v||^2.

Orthogonal Complements and Direct Sum Decomposition#

Let W be a subspace of an inner product space V. The orthogonal complement of W, denoted W^⊥ (pronounced "W-perp"), is the set of all vectors in V that are orthogonal to every vector in W:

W^⊥ = { v in V : <v, w> = 0 for all w in W }.

Theorem: If W is a finite-dimensional subspace of V, then:

1. W^⊥ is a valid linear subspace of V.
2. The only vector belonging to both W and W^⊥ is the zero vector: W ∩ W^⊥ = {0}.
3. The space V decomposes into the direct sum of W and its orthogonal complement:

V = W ⊕ W^⊥.

Every vector y in V can be uniquely represented as:

y = w + z, where w in W and z in W^⊥.

3. Orthogonal Projections and Projection Matrices#

When an exact solution to a linear system does not exist, finding the closest possible approximation within a given subspace requires computing an orthogonal projection.

Orthogonal Projection onto a 1-Dimensional Subspace#

The orthogonal projection of a vector y onto a 1-dimensional line spanned by a non-zero vector u is given by:

proj_u(y) = (<y, u> / <u, u>) · u = (<y, u> / ||u||^2) · u.

The vector z = y - proj_u(y) represents the perpendicular error vector, satisfying <z, u> = 0.

Projection onto an Arbitrary Subspace#

Let W be a k-dimensional subspace of R^m. The orthogonal projection of y onto W is the unique vector y_hat in W that minimizes the Euclidean distance ||y - y_hat|| across all candidate vectors in W:

y_hat = proj_W(y) = argmin_{w in W} ||y - w||.

Case 1: When W is spanned by an Orthonormal Basis { q1, q2, ..., qk }:

The projection decomposes into the uncoupled sum of projections onto each basis axis:

proj_W(y) = sum(i=1 to k, <y, qi> · qi).

If we collect the orthonormal basis vectors into the columns of an m×k matrix Q = [ q1 | q2 | ... | qk ], where Q^T · Q = I_k, the projection operator becomes a matrix multiplication:

proj_W(y) = Q · (Q^T · y) = (Q · Q^T) · y.

Hence, the orthogonal projection matrix onto subspace W is:

P_W = Q · Q^T.

Case 2: When W is spanned by an Arbitrary (Non-Orthogonal) Basis A:

If A = [ a1 | a2 | ... | ak ] has linearly independent columns that are not orthonormal, the normal equations give:

P_W = A · (A^T · A)^(-1) · A^T.

       Vector y
          ^
          |           |  \  Perpendicular Error Vector:
          |   \   z = y - y_hat  (z in W^perp)
          |              +-----> y_hat = P_W · y  (y_hat in Subspace W)
        [ Subspace W ]

Algebraic Invariants of Projection Matrices#

Every orthogonal projection matrix P satisfies two fundamental algebraic properties:

1. Idempotence: P^2 = P. Projecting an already-projected vector onto the subspace leaves it unaltered.
2. Symmetry: P^T = P. Orthogonal projection operators are self-adjoint.
3. Spectrum: The eigenvalues of any projection matrix P can only be 0 or 1. Eigenvectors with eigenvalue 1 span the subspace W; eigenvectors with eigenvalue 0 span the orthogonal complement W^⊥.

4. The Gram-Schmidt Process: Classical vs. Modified Formulations#

The Gram-Schmidt algorithm is the foundational numerical engine for transforming an arbitrary set of linearly independent vectors into an orthonormal basis spanning the exact same nested subspaces.

The Classical Gram-Schmidt (CGS) Formulation#

Given linearly independent vectors { a1, a2, ..., an }, CGS constructs orthonormal vectors { q1, q2, ..., qn } iteratively:

For step k = 1 to n:

1. Compute the unnormalized orthogonal vector v_k by subtracting projections onto all previously computed unit vectors:

v_k = a_k - sum(j=1 to k-1, (a_k · q_j) · q_j).

1. Normalize the vector to unit length:

q_k = v_k / ||v_k||.

The Catastrophic Failure Mode of CGS: Loss of Orthogonality#

In pure mathematical theory with infinite precision, CGS produces a perfectly orthonormal set. However, in finite-precision floating-point arithmetic (IEEE-754 double precision), CGS suffers from catastrophic loss of orthogonality.

When column vectors a_k are nearly collinear (ill-conditioned matrices), subtracting the inner product terms involves catastrophic cancellation. Small floating-point rounding errors in calculating (a_k · q_j) fail to completely remove the directional component along q_j. As subsequent vectors are processed, these residual non-zero components accumulate exponentially, causing computed vectors q_i and q_j to lose perpendicularity entirely:

||Q^T · Q - I|| >> ε_machine.

The Modified Gram-Schmidt (MGS) Formulation#

To resolve this instability without resorting to computationally expensive Householder reflections, the Modified Gram-Schmidt (MGS) algorithm updates the working vectors sequentially. Instead of subtracting projections of the original vector a_k against previously finalized vectors, MGS deflates all remaining candidate vectors immediately after each new unit vector is finalized.

+-----------------------------------------------------------------------------+
|                     MODIFIED GRAM-SCHMIDT (MGS) ALGORITHM                   |
+-----------------------------------------------------------------------------+
| Input:  Matrix A in R^(m x n) with linearly independent columns             |
| Output: Matrix Q in R^(m x n) (orthonormal), Matrix R in R^(n x n) (triang) |
+-----------------------------------------------------------------------------+
| 1. Initialize working copy: V = copy(A), R = zeros(n, n)                    |
| 2. For i = 1 to n:                                                          |
|      R[i, i] = norm(V[:, i])                                                |
|      Q[:, i] = V[:, i] / R[i, i]                                            |
|      For j = i + 1 to n:                                                    |
|        R[i, j] = dot_product(Q[:, i], V[:, j])                              |
|        V[:, j] = V[:, j] - R[i, j] * Q[:, i]  <-- Immediate in-place update  |
+-----------------------------------------------------------------------------+

By immediately purging the component along q_i from all remaining vectors V[:, j], MGS drastically suppresses roundoff error amplification, ensuring that computed vectors maintain orthogonality to within O(ε_machine · κ(A)).

Comparative Matrix Orthogonalization Taxonomy#

5. QR Factorization Mechanics & Overdetermined Least-Squares Systems#

The Gram-Schmidt process can be expressed elegantly as a matrix factorization: decomposing a matrix A into the product of an orthogonal matrix Q and an upper triangular matrix R.

Thin (Reduced) vs. Full QR Factorization#

For any real matrix A in R^(m×n) with m >= n:

1. Reduced (Thin) QR Factorization:

A = Q · R

- Q in R^(m×n) has orthonormal columns: Q^T · Q = I_n.

- R in R^(n×n) is an upper triangular, invertible matrix with strictly positive diagonal entries (R_ii > 0).

1. Full QR Factorization:

A = Q_full · R_full

- Q_full in R^(m×m) is a square orthogonal matrix: Q_full^T · Q_full = Q_full · Q_full^T = I_m.

- R_full in R^(m×n) consists of an upper triangular n×n block stacked on top of an (m - n)×n block of zeros.

Solving Overdetermined Least-Squares Systems via QR#

In data science, machine learning, and parameter estimation, overdetermined linear systems A · x ≈ b arise whenever the number of experimental observations m exceeds the number of unknown parameters n.

The classical method solves the Normal Equations:

A^T · A · x = A^T · b

The Fatal Flaw of Normal Equations:

Computing the cross-product matrix A^T · A squares the matrix condition number:

κ(A^T · A) = (κ(A))^2.

If a data matrix has a modest condition number κ(A) = 10^5, the normal equations matrix has κ(A^T · A) = 10^10. In IEEE-754 double precision (approx. 16 decimal digits), this results in the catastrophic loss of 10 digits of precision, rendering the computed parameters unreliable.

The QR Solution Method:

Substitute the thin QR factorization A = Q · R into the normal equations:

(Q · R)^T · (Q · R) · x = (Q · R)^T · b

R^T · (Q^T · Q) · R · x = R^T · Q^T · b

Because Q has orthonormal columns, Q^T · Q = I_n:

R^T · R · x = R^T · Q^T · b

Because A has full column rank, R is invertible. Multiplying by (R^T)^(-1) yields:

R · x = Q^T · b

This is an upper triangular system! It is solved directly by backward substitution in O(n^2) operations. Most importantly, backward substitution on R solves the problem using the original condition number κ(A), completely avoiding condition number squaring.

6. Systems, Robotics & Machine Learning Applications#

Orthogonal projections and QR decompositions are indispensable tools across multiple engineering domains.

1. Embedded Robotics & Kalman Filtering: Square Root Filtering#

In autonomous drones, self-driving vehicles, and spacecraft navigation, state estimation is governed by the Extended Kalman Filter (EKF). The filter continuously updates the state covariance matrix P. Under standard matrix subtraction P = (I - K · H) · P_prior, floating-point roundoff errors can cause P to lose its positive-definiteness, leading to negative variances and filter divergence.

Production aerospace systems deploy Square Root Information Filters (SRIF) using QR factorization. Instead of propagating the full covariance matrix P, the algorithm propagates its Cholesky factor S where P = S · S^T. At each measurement step, a QR decomposition updates S directly:

[ R_measurement, H · S_prior ]^T = Q · [ S_new, 0 ]^T

Because S_new is guaranteed to be upper triangular with positive diagonal, the reconstructed covariance matrix P = S_new · S_new^T is mathematically guaranteed to remain symmetric and positive-definite, preventing filter failure in flight.

2. Computer Vision: Camera Extrinsic & Intrinsic Calibration#

In photogrammetry, SLAM (Simultaneous Localization and Mapping), and augmented reality, a camera's projection matrix P (size 3×4) maps 3D world coordinates to 2D image coordinates:

P = K · [ R | t ]

where:

• K (3×3) is the upper triangular intrinsic calibration matrix containing focal lengths, principal point offsets, and skew.
• R (3×3) is the orthogonal rotation matrix describing camera orientation.
• t (3×1) is the camera translation vector.

To separate the camera's internal optics from its spatial orientation, computer vision libraries (OpenCV) perform an RQ decomposition on the leading 3×3 submatrix of P. Because RQ factorizes a matrix into an upper triangular factor (K) and an orthogonal factor (R), the camera parameters are decoupled cleanly without numerical distortion.

3. Machine Learning: Ordinary Least Squares & Ridge Regression#

In high-throughput machine learning engines (Scikit-learn, PyTorch, TensorFlow), solving linear regression models y = X · β + ε relies on QR factorization rather than explicit matrix inversion. When fitting millions of samples with correlated features, QR factorization provides stable coefficients and allows fast residual norm calculation without computing the predicted values explicitly:

||y - X · β_hat||_2 = ||(I - Q · Q^T) · y||_2.

7. Common Architectural Pitfalls & Misconceptions#

• Pitfall 1: Assuming Q · Q^T = I in Thin QR Factorization.

When A is an m×n rectangular matrix (m > n), the matrix Q in thin QR has dimension m×n. The product Q^T · Q = I_n (an n×n identity matrix). However, the reverse product Q · Q^T is an m×m matrix representing the orthogonal projection operator P onto Col(A). It equals I_m if and only if m = n (i.e., when A is square and invertible).

• Pitfall 2: Implementing Classical Gram-Schmidt in Production Floating-Point Software.

Due to severe loss of orthogonality, CGS should never be used in numerical production code. Always use Modified Gram-Schmidt (MGS) or Householder QR reflectors.

• Pitfall 3: Explicitly Inverting Matrix R to Solve R · x = Q^T · b.

Never call an explicit matrix inversion routine inv(R) to solve a triangular system. Inverting R requires O(n^3) operations and introduces unnecessary rounding errors. Always use backward substitution, which takes O(n^2) operations and is backward stable.

• Pitfall 4: Neglecting Column Rank Deficiency in QR Factorization.

If the columns of A are linearly dependent (rank(A) < n), the standard Gram-Schmidt process will encounter a zero norm (||v_k|| = 0), triggering a division-by-zero error. Production QR solvers apply column pivoting (QR with Pivoting: A · P = Q · R), dynamically swapping the column with the largest remaining norm to ensure robust diagonal decay.

8. Summary & Checkpoint Exercises#

Solidify your comprehension of inner product spaces, orthogonal projections, Modified Gram-Schmidt, and QR factorization with these four checkpoint problems. Complete step-by-step solutions are provided.

Problem 1 (Verifying Custom Weighted Inner Product & Cauchy-Schwarz):

Consider the vector space R^2. Define a candidate inner product with weighting coefficients:

<u, v>_W = 3 · u1 · v1 + 2 · u2 · v2, where u = (u1, u2) and v = (v1, v2).

1. Prove that <u, v>_W satisfies all four inner product axioms.
2. Given vectors u = (2, -1) and v = (1, 3), compute <u, v>_W, ||u||_W, ||v||_W, and verify that the Cauchy-Schwarz inequality holds strictly.

Solution:

Part 1: Axiom Verification.

1. Symmetry:

<u, v>_W = 3 · u1 · v1 + 2 · u2 · v2 = 3 · v1 · u1 + 2 · v2 · u2 = <v, u>_W. Holds.

1. Linearity in First Argument:

Let w = (w1, w2) and c in R.

<c · u + w, v>_W = 3 · (c · u1 + w1) · v1 + 2 · (c · u2 + w2) · v2

= c · (3 · u1 · v1 + 2 · u2 · v2) + (3 · w1 · v1 + 2 · w2 · v2)

= c · <u, v>_W + <w, v>_W. Holds.

1. Positive-Definiteness:

<v, v>_W = 3 · v1^2 + 2 · v2^2.

Since squares of real numbers are non-negative and the weights (3 and 2) are strictly positive, <v, v>_W >= 0.

Furthermore, 3 · v1^2 + 2 · v2^2 = 0 if and only if v1 = 0 and v2 = 0, meaning v = (0, 0). Holds.

1. Induced Norm:

||v||_W = sqrt(3 · v1^2 + 2 · v2^2). Valid norm induced.

Conclusion: <u, v>_W is a valid inner product on R^2.

Part 2: Calculation and Cauchy-Schwarz Verification.

For u = (2, -1) and v = (1, 3):

• Inner Product:

<u, v>_W = 3 · (2)(1) + 2 · (-1)(3) = 3(2) + 2(-3) = 6 - 6 = 0.

• Induced Norms:

||u||_W = sqrt(3 · (2)^2 + 2 · (-1)^2) = sqrt(3 · 4 + 2 · 1) = sqrt(12 + 2) = sqrt(14).

||v||_W = sqrt(3 · (1)^2 + 2 · (3)^2) = sqrt(3 · 1 + 2 · 9) = sqrt(3 + 18) = sqrt(21).

• Cauchy-Schwarz Check:

|<u, v>_W| = |0| = 0.

||u||_W · ||v||_W = sqrt(14) · sqrt(21) = sqrt(294) ≈ 17.146.

Since 0 < 17.146, the inequality |<u, v>_W| <= ||u||_W · ||v||_W holds strictly. (In fact, because <u, v>_W = 0, u and v are orthogonal under this weighted metric!).

Problem 2 (Step-by-Step Modified Gram-Schmidt Orthogonalization):

Apply the Modified Gram-Schmidt (MGS) algorithm to construct a thin QR factorization (A = Q · R) for the 3×2 matrix:

A = [

[ 1, 2 ],

[ 1, 1 ],

[ 0, 1 ]

]

Solution:

Let column 1 be a1 = [ 1, 1, 0 ]^T and column 2 be a2 = [ 2, 1, 1 ]^T.

Initialize working vectors: v1 = [ 1, 1, 0 ]^T, v2 = [ 2, 1, 1 ]^T.

Iteration 1 (Processing Column 1):

1. Compute the norm of v1:

R_11 = ||v1|| = sqrt((1)^2 + (1)^2 + (0)^2) = sqrt(1 + 1 + 0) = sqrt(2).

1. Compute the first orthonormal unit vector q1:

q1 = v1 / R_11 = [ 1/sqrt(2), 1/sqrt(2), 0 ]^T.

1. Deflate remaining candidate vectors (v2):

Project v2 onto q1:

R_12 = <q1, v2> = (1/sqrt(2))(2) + (1/sqrt(2))(1) + (0)(1) = 2/sqrt(2) + 1/sqrt(2) = 3/sqrt(2).

Subtract projection from v2 in-place:

v2_new = v2 - R_12 · q1

v2_new = [ 2, 1, 1 ]^T - (3/sqrt(2)) · [ 1/sqrt(2), 1/sqrt(2), 0 ]^T

v2_new = [ 2, 1, 1 ]^T - [ 3/2, 3/2, 0 ]^T

v2_new = [ 2 - 1.5, 1 - 1.5, 1 - 0 ]^T = [ 1/2, -1/2, 1 ]^T.

Iteration 2 (Processing Column 2):

1. Compute the norm of the deflated vector v2_new:

||v2_new||^2 = (1/2)^2 + (-1/2)^2 + (1)^2 = 1/4 + 1/4 + 1 = 2/4 + 1 = 3/2.

R_22 = sqrt(3/2) = sqrt(3) / sqrt(2).

1. Compute the second orthonormal unit vector q2:

q2 = v2_new / R_22 = (sqrt(2)/sqrt(3)) · [ 1/2, -1/2, 1 ]^T

q2 = [ 1 / (2 · sqrt(3/2)), -1 / (2 · sqrt(3/2)), 1 / sqrt(3/2) ]^T

q2 = [ 1/sqrt(6), -1/sqrt(6), sqrt(2/3) ]^T = [ 1/sqrt(6), -1/sqrt(6), 2/sqrt(6) ]^T.

Assembling Matrices Q and R:

Q = [

[ 1/sqrt(2), 1/sqrt(6) ],

[ 1/sqrt(2), -1/sqrt(6) ],

[ 0, 2/sqrt(6) ]

]

R = [

[ sqrt(2), 3/sqrt(2) ],

[ 0, sqrt(3/2) ]

]

Verification:

• Q^T · Q:

||q1||^2 = (1/sqrt(2))^2 + (1/sqrt(2))^2 + 0 = 1/2 + 1/2 = 1.

||q2||^2 = (1/sqrt(6))^2 + (-1/sqrt(6))^2 + (2/sqrt(6))^2 = 1/6 + 1/6 + 4/6 = 6/6 = 1.

<q1, q2> = (1/sqrt(2))(1/sqrt(6)) + (1/sqrt(2))(-1/sqrt(6)) + (0)(2/sqrt(6)) = 1/sqrt(12) - 1/sqrt(12) + 0 = 0.

• Q · R reconstruction:

Entry (1, 1): (1/sqrt(2)) · sqrt(2) + 0 = 1.

Entry (2, 1): (1/sqrt(2)) · sqrt(2) + 0 = 1.

Entry (3, 1): 0 · sqrt(2) + 0 = 0.

Entry (1, 2): (1/sqrt(2)) · (3/sqrt(2)) + (1/sqrt(6)) · (sqrt(3/2)) = 3/2 + sqrt(3)/(sqrt(12)) = 3/2 + 1/2 = 2.

Entry (2, 2): (1/sqrt(2)) · (3/sqrt(2)) + (-1/sqrt(6)) · (sqrt(3/2)) = 3/2 - 1/2 = 1.

Entry (3, 2): 0 + (2/sqrt(6)) · (sqrt(3/2)) = 2 · sqrt(3) / sqrt(12) = 2 · (1/2) = 1.

Matrix product Q · R exactly reproduces matrix A.

Problem 3 (Subspace Projection Matrix and Error Vector Derivation):

Let W be the 2-dimensional subspace of R^3 spanned by vectors a1 = [ 1, 0, 1 ]^T and a2 = [ 0, 1, 1 ]^T.

Given the target vector y = [ 2, 4, 3 ]^T:

1. Compute the orthogonal projection matrix P_W using the formula P_W = A · (A^T · A)^(-1) · A^T.
2. Find the projected vector y_hat = proj_W(y).
3. Find the orthogonal error vector z = y - y_hat, and verify that z is orthogonal to both a1 and a2.

Solution:

Step 1: Construct A and compute A^T · A.

A = [

[ 1, 0 ],

[ 0, 1 ],

[ 1, 1 ]

]

A^T · A = [

[ 1, 0, 1 ],

[ 0, 1, 1 ]

] · [

[ 1, 0 ],

[ 0, 1 ],

[ 1, 1 ]

]

• Entry (1, 1): 1(1) + 0(0) + 1(1) = 2
• Entry (1, 2): 1(0) + 0(1) + 1(1) = 1
• Entry (2, 1): 0(1) + 1(0) + 1(1) = 1
• Entry (2, 2): 0(0) + 1(1) + 1(1) = 2

A^T · A = [

[ 2, 1 ],

[ 1, 2 ]

]

Compute (A^T · A)^(-1):

det(A^T · A) = (2)(2) - (1)(1) = 4 - 1 = 3.

(A^T · A)^(-1) = (1/3) · [

[ 2, -1 ],

[ -1, 2 ]

]

Compute A · (A^T · A)^(-1):

A · (A^T · A)^(-1) = (1/3) · [

[ 1, 0 ],

[ 0, 1 ],

[ 1, 1 ]

] · [

[ 2, -1 ],

[ -1, 2 ]

]

= (1/3) · [

[ 2, -1 ],

[ -1, 2 ],

[ 1, 1 ]

]

Now multiply by A^T:

P_W = (1/3) · [

[ 2, -1 ],

[ -1, 2 ],

[ 1, 1 ]

] · [

[ 1, 0, 1 ],

[ 0, 1, 1 ]

]

• Row 1: [ 2(1) + (-1)(0), 2(0) + (-1)(1), 2(1) + (-1)(1) ] = [ 2, -1, 1 ]
• Row 2: [ -1(1) + 2(0), -1(0) + 2(1), -1(1) + 2(1) ] = [ -1, 2, 1 ]
• Row 3: [ 1(1) + 1(0), 1(0) + 1(1), 1(1) + 1(1) ] = [ 1, 1, 2 ]

Thus, the projection matrix is:

P_W = (1/3) · [

[ 2, -1, 1 ],

[ -1, 2, 1 ],

[ 1, 1, 2 ]

]

Step 2: Compute y_hat = P_W · y.

y = [ 2, 4, 3 ]^T.

P_W · y = (1/3) · [

[ 2, -1, 1 ],

[ -1, 2, 1 ],

[ 1, 1, 2 ]

] · [

2,

4,

3

]

• Row 1: 2(2) - 1(4) + 1(3) = 4 - 4 + 3 = 3 => (1/3)(3) = 1
• Row 2: -1(2) + 2(4) + 1(3) = -2 + 8 + 3 = 9 => (1/3)(9) = 3
• Row 3: 1(2) + 1(4) + 2(3) = 2 + 4 + 6 = 12 => (1/3)(12) = 4

Projected vector:

y_hat = [ 1, 3, 4 ]^T.

Notice that y_hat = 1 · a1 + 3 · a2 = 1[ 1, 0, 1 ]^T + 3[ 0, 1, 1 ]^T = [ 1, 3, 4 ]^T, confirming y_hat in W.

Step 3: Compute and Verify Error Vector z.

z = y - y_hat = [ 2, 4, 3 ]^T - [ 1, 3, 4 ]^T = [ 1, 1, -1 ]^T.

Verify orthogonality to basis vectors:

• <z, a1> = 1(1) + 1(0) + (-1)(1) = 1 + 0 - 1 = 0.
• <z, a2> = 1(0) + 1(1) + (-1)(1) = 0 + 1 - 1 = 0.

Both inner products vanish. The error vector z is strictly orthogonal to subspace W.

Problem 4 (Solving Least Squares via Thin QR and Back-Substitution):

Given an overdetermined system A · x ≈ b with:

A = [

[ 1, 1 ],

[ 1, 2 ],

[ 1, 3 ]

], b = [ 6, 5, 7 ]^T

Given that the thin QR factorization of A is:

Q = [

[ 1/sqrt(3), -1/sqrt(2) ],

[ 1/sqrt(3), 0 ],

[ 1/sqrt(3), 1/sqrt(2) ]

]

R = [

[ sqrt(3), 2 · sqrt(3) ],

[ 0, 1/sqrt(2) ]

]

1. Compute the projected vector c = Q^T · b.
2. Solve the upper triangular system R · x = c using backward substitution.
3. Compute the least-squares residual error ||A · x_hat - b||_2.

Solution:

Step 1: Compute c = Q^T · b.

Q^T = [

[ 1/sqrt(3), 1/sqrt(3), 1/sqrt(3) ],

[ -1/sqrt(2), 0, 1/sqrt(2) ]

]

b = [ 6, 5, 7 ]^T.

• Entry c1:

c1 = (1/sqrt(3))(6) + (1/sqrt(3))(5) + (1/sqrt(3))(7) = (6 + 5 + 7) / sqrt(3) = 18 / sqrt(3) = 6 · sqrt(3).

• Entry c2:

c2 = (-1/sqrt(2))(6) + (0)(5) + (1/sqrt(2))(7) = (-6 + 0 + 7) / sqrt(2) = 1 / sqrt(2).

Thus:

c = [ 6 · sqrt(3), 1/sqrt(2) ]^T.

Step 2: Backward Substitution on R · x = c.

[

[ sqrt(3), 2 · sqrt(3) ],

[ 0, 1/sqrt(2) ]

] · [ x1, x2 ]^T = [ 6 · sqrt(3), 1/sqrt(2) ]^T

From row 2:

(1/sqrt(2)) · x2 = 1/sqrt(2) => x2 = 1.

From row 1:

sqrt(3) · x1 + 2 · sqrt(3) · x2 = 6 · sqrt(3)

Divide through by sqrt(3):

x1 + 2 · x2 = 6

x1 + 2(1) = 6 => x1 = 4.

Least-squares parameter solution:

x_hat = [ 4, 1 ]^T.

Step 3: Compute Residual Error.

Calculate predicted vector y_pred = A · x_hat:

y_pred = [

[ 1, 1 ],

[ 1, 2 ],

[ 1, 3 ]

] · [ 4, 1 ]^T

• Row 1: 1(4) + 1(1) = 5
• Row 2: 1(4) + 2(1) = 6
• Row 3: 1(4) + 3(1) = 7

y_pred = [ 5, 6, 7 ]^T.

Calculate residual error vector:

e = b - y_pred = [ 6, 5, 7 ]^T - [ 5, 6, 7 ]^T = [ 1, -1, 0 ]^T.

Calculate Euclidean norm:

||e||_2 = sqrt((1)^2 + (-1)^2 + (0)^2) = sqrt(1 + 1 + 0) = sqrt(2) ≈ 1.4142.

The optimal least-squares parameter vector is x_hat = [ 4, 1 ]^T with residual norm sqrt(2).