Singular Value Decomposition, Moore-Penrose Pseudoinverse & Low-Rank Approximation

While eigendecomposition requires square, non-defective matrices, the Singular Value Decomposition (SVD) represents the master theorem of linear algebra because it applies unconditionally to every matrix—regardless of dimensions, rank deficiency, or symmetry. SVD provides complete insight into the geometric transformation of space, establishes explicit orthonormal bases for all four fundamental subspaces, delivers optimal low-rank matrix compression, and forms the mathematical backbone of modern machine learning, computer vision, and Large Language Model optimization.

1. Geometric Intuition & The SVD Master Theorem#

In Euclidean geometry, every linear transformation A : R^n → R^m can be visualized as mapping geometric shapes in the domain into geometric shapes in the codomain. A unit sphere in R^n (the set of all vectors x with ||x|| = 1) is invariably transformed by matrix A into a hyper-ellipsoid in R^m.

The Singular Value Decomposition discovers the perpendicular axes of this hyper-ellipsoid:

• The directions of the principal semi-axes are given by the left singular vectors u_i.
• The lengths of the principal semi-axes are given by the singular values σ_i.
• The original perpendicular directions on the unit sphere that mapped to these semi-axes are the right singular vectors v_i.

The Fundamental SVD Theorem#

Theorem: Every real rectangular matrix A in R^(m×n) of rank r <= min(m, n) can be factored into the product of three matrices:

A = U · Σ · V^T

where:

1. U in R^(m×m) is a square orthogonal matrix:

U^T · U = U · U^T = I_m.

The column vectors of U, { u_1, u_2, ..., u_m }, are the left singular vectors. They form an orthonormal basis for R^m and are the eigenvectors of the symmetric matrix A · A^T.

1. Σ in R^(m×n) is a rectangular diagonal matrix:

The main diagonal contains r strictly positive real numbers arranged in non-increasing order:

σ_1 >= σ_2 >= ... >= σ_r > 0

All remaining entries are zero:

σ_{r+1} = σ_{r+2} = ... = σ_{min(m, n)} = 0.

These scalars σ_i are the singular values of matrix A. They are uniquely determined as the positive square roots of the non-zero eigenvalues of A^T · A (or A · A^T).

1. V in R^(n×n) is a square orthogonal matrix:

V^T · V = V · V^T = I_n.

The column vectors of V, { v_1, v_2, ..., v_n }, are the right singular vectors. They form an orthonormal basis for R^n and are the eigenvectors of the symmetric matrix A^T · A.

+-----------------------------------------------------------------------------+
|                     GEOMETRIC FACTORIZATION OF SVD                          |
+-----------------------------------------------------------------------------+
|                                                                             |
|  Unit Sphere in R^n                                 Hyper-Ellipsoid in R^m  |
|       { v_1, v_2 }        --- A = U * Sigma * V^T --->  { sigma_1 * u_1,    |
|                                                           sigma_2 * u_2 }   |
|                                                                             |
|  Three-Step Geometric Transformation:                                       |
|  1. Rotate / Reflect in R^n via V^T (align axes with right singular vectors)|
|  2. Scale axes by singular values sigma_i & change dimension via Sigma       |
|  3. Rotate / Reflect in R^m via U (align with left singular vectors)        |
+-----------------------------------------------------------------------------+

Derivation from the Spectral Theorem#

Consider the symmetric Gram matrix S = A^T · A in R^(n×n). By the Spectral Theorem for real symmetric matrices:

1. All eigenvalues of A^T · A are real and non-negative:

For any x in R^n, x^T · (A^T · A) · x = (A · x)^T · (A · x) = ||A · x||^2 >= 0.

Hence, λ_i >= 0.

1. There exists an orthonormal eigenbasis { v_1, v_2, ..., v_n } such that:

(A^T · A) · v_i = λ_i · v_i.

1. Define the singular values as σ_i = sqrt(λ_i) for i = 1, ..., r.
2. For i = 1, ..., r, define the left singular vectors:

u_i = (1 / σ_i) · A · v_i.

Check their orthonormality:

<u_i, u_j> = u_i^T · u_j = (1 / (σ_i · σ_j)) · (A · v_i)^T · (A · v_j)

= (1 / (σ_i · σ_j)) · v_i^T · (A^T · A) · v_j = (1 / (σ_i · σ_j)) · v_i^T · (λ_j · v_j)

= (λ_j / (σ_i · σ_j)) · <v_i, v_j>.

Because { v_i } are orthonormal:

- If i ≠ j, <v_i, v_j> = 0, so <u_i, u_j> = 0.

- If i = j, λ_i / σ_i^2 = λ_i / λ_i = 1, so <u_i, u_i> = 1.

Hence, { u_1, ..., u_r } form an orthonormal set in R^m, which can be extended via Gram-Schmidt to an orthonormal basis for R^m. In matrix form, A · v_i = σ_i · u_i, which yields A · V = U · Σ, or A = U · Σ · V^T.

2. Complete Unification of the Four Fundamental Subspaces#

SVD achieves what no other matrix decomposition can: it provides an explicit, decoupled orthonormal basis for all four fundamental subspaces established by the Fundamental Theorem of Linear Algebra.

The Four Subspaces Mapped by SVD#

Let A in R^(m×n) have rank r and SVD factorization A = U · Σ · V^T:

Because A · v_i = σ_i · u_i for i = 1, ..., r, matrix A acts as a diagonal scaling between the row space and the column space!

For i > r, A · v_i = 0 · u_i = 0, meaning { v_{r+1}, ..., v_n } span the exact kernel that collapses to zero.

SVD Variants: Full, Thin, and Compact#

Depending on computational and memory requirements, SVD is deployed in three standard configurations:

1. Full SVD:

A = U_(m×m) · Σ_(m×n) · V^T_(n×n)

Preserves the complete orthogonal coordinate frames for both entire ambient spaces R^m and R^n.

1. Thin (Economy) SVD:

For m >= n, discards the last m - n columns of U and rows of Σ that correspond to the left null space:

A = U_(m×n) · Σ_(n×n) · V^T_(n×n)

where U^T · U = I_n.

1. Compact SVD:

Retains strictly the r non-zero singular values and their corresponding singular vectors:

A = U_r · Σ_r · V_r^T

where U_r is m×r, Σ_r is r×r, and V_r is n×r.

Numerical Rank & Condition Number#

In digital computation, rounding errors prevent exact zeros from appearing in matrices. Counting exact non-zero pivot elements is mathematically unsound. SVD provides the gold standard for determining numerical rank:

A singular value σ_i is treated as zero if:

σ_i <= τ = ε_machine · max(m, n) · σ_1.

The matrix condition number under the Euclidean L2 operator norm is computed directly from SVD:

κ(A) = ||A||_2 · ||A^+||_2 = σ_max / σ_min = σ_1 / σ_r.

3. The Eckart-Young-Mirsky Low-Rank Approximation Theorem#

In dyadic (outer product) form, matrix A can be expressed as a linear combination of r rank-one matrices:

A = sum(i=1 to r, σ_i · (u_i · v_i^T))

Each term E_i = u_i · v_i^T is a rank-one matrix with unit operator norm (||u_i · v_i^T||_2 = 1) and unit Frobenius norm (||u_i · v_i^T||_F = 1). The singular values σ_i act as decreasing weight coefficients prioritizing each structural mode.

The Truncated SVD Formulation#

The rank-k truncated SVD approximation A_k (for k < r) is formed by retaining only the first k dominant terms:

A_k = sum(i=1 to k, σ_i · (u_i · v_i^T)) = U_k · Σ_k · V_k^T

where U_k in R^(m×k), Σ_k in R^(k×k), and V_k in R^(n×k).

The Eckart-Young-Mirsky Theorem#

Theorem (Eckart & Young 1936, Mirsky 1960):

Let A be an m×n real matrix of rank r, and let k < r. The truncated SVD matrix A_k is the strictly optimal rank-k approximation of A under both the Spectral (L2) Norm and the Frobenius Norm:

1. Optimality under Spectral Norm:

min_{B : rank(B) <= k} ||A - B||_2 = ||A - A_k||_2 = σ_{k+1}

1. Optimality under Frobenius Norm:

min_{B : rank(B) <= k} ||A - B||_F = ||A - A_k||_F = sqrt(sum(i=k+1 to r, σ_i^2))

The theorem proves that no matrix of rank k or less can approximate A better than A_k, and the exact error is dictated exclusively by the discarded singular values!

Massive Data Compression Economics#

Consider a 1000×1000 dense image matrix A (e.g., 1 megapixel grayscale). Storing A requires:

1000 × 1000 = 1,000,000 floating-point numbers.

Using a rank-20 approximation A_20:

• U_20: 1000 × 20 = 20,000 floats
• Σ_20: 20 diagonal floats
• V_20: 1000 × 20 = 20,000 floats

Total storage = 20,000 + 20 + 20,000 = 40,020 floats.

Compression Ratio: 40,020 / 1,000,000 ≈ 4% of original size (a 96% data reduction!) while preserving over 95% of the visual fidelity.

4. The Moore-Penrose Pseudoinverse & Regularization#

When a matrix A in R^(m×n) is rectangular (m ≠ n) or rank-deficient (rank(A) < min(m, n)), the classical inverse A^(-1) does not exist. The Moore-Penrose Pseudoinverse A^+ in R^(n×m) generalizes matrix inversion unconditionally.

Definition via SVD#

Using the SVD factorization A = U · Σ · V^T, the pseudoinverse is defined as:

A^+ = V · Σ^+ · U^T

where Σ^+ in R^(n×m) is formed by transposing Σ and taking the reciprocal of every strictly non-zero singular value:

(Σ^+)_ii = 1 / σ_i if σ_i > 0, and 0 if σ_i = 0.

The Four Penrose Conditions#

Theorem: The Moore-Penrose pseudoinverse A^+ is the unique matrix in R^(n×m) satisfying all four Penrose algebraic equations:

1. A · A^+ · A = A
2. A^+ · A · A^+ = A^+
3. (A · A^+)^T = A · A^+ (A · A^+ is the orthogonal projection matrix onto Col(A))
4. (A^+ · A)^T = A^+ · A (A^+ · A is the orthogonal projection matrix onto Row(A))

The Minimum-Norm Least-Squares Solution#

For any general linear system A · x = b:

• If m > n (overdetermined): No exact solution may exist.
• If m < n (underdetermined): Infinitely many solutions exist.
• If rank(A) < min(m, n) (rank-deficient): Both inconsistencies and null space ambiguities arise.

Theorem: The vector x_hat = A^+ · b is the unique optimal solution satisfying two simultaneous criteria:

1. Least-Squares Minimization: x_hat minimizes the residual error norm:

||A · x_hat - b||_2 = min_{x} ||A · x - b||_2.

1. Minimum Euclidean Norm: Among all vectors that achieve this minimal residual error, x_hat has the smallest Euclidean length:

||x_hat||_2 = min { ||x||_2 : ||A · x - b||_2 is minimal }.

Truncated SVD (TSVD) & Tikhonov Regularization#

In physical inverse problems (e.g., medical CT scans, seismology, deblurring), matrix A is frequently ill-conditioned, possessing tiny singular values σ_i ≈ 10^(-8).

In computing A^+ · b, these tiny values become massive multipliers: 1 / σ_i ≈ 10^8. Any microscopic sensor noise in b is magnified by 10^8, completely corrupting x_hat.

Two standard regularization techniques resolve this instability:

1. Truncated SVD (TSVD):

Set 1 / σ_i = 0 for all σ_i < τ (a chosen noise threshold). Small singular values are completely suppressed.

1. Tikhonov Regularization (Ridge Regression):

Instead of hard thresholding, apply smooth filter factors:

f_i = σ_i^2 / (σ_i^2 + α)

where α > 0 is the regularization parameter.

The pseudoinverse diagonal entries become:

(Σ^+)_ii = σ_i / (σ_i^2 + α).

When σ_i >> sqrt(α), the entry approaches 1/σ_i. When σ_i << sqrt(α), the entry smoothly decays toward 0, preventing noise amplification.

5. Machine Learning, NLP & Generative AI Applications#

SVD is one of the most widely deployed mathematical frameworks in modern Artificial Intelligence.

1. Principal Component Analysis (PCA)#

PCA is the standard unsupervised technique for dimensionality reduction, feature extraction, and high-dimensional data visualization.

Given a dataset of N observations with D features stored in matrix X in R^(N×D):

1. Mean-center the data: Subtract the column mean from every feature:

X_centered = X - 1 · μ^T.

1. The sample covariance matrix is C = (1 / (N - 1)) · X_centered^T · X_centered.
2. Compute the SVD of the centered data matrix:

X_centered = U · Σ · V^T.

Relationships Revealed:

• The right singular vectors (columns of V) are the principal component directions (eigenvectors of covariance matrix C).
• The singular values σ_i are related to the variances along principal directions:

Variance_i = σ_i^2 / (N - 1).

• The projected low-dimensional coordinates are given by:

Z_k = X_centered · V_k = U_k · Σ_k in R^(N×k).

SVD computes PCA directly on X_centered in O(N · D · min(N, D)) flops, avoiding the formation and squaring of the covariance matrix C!

2. Latent Semantic Analysis (LSA) & Recommender Systems#

In natural language processing, a text corpus is represented as a Term-Document matrix A in R^(words × documents), where A_ij is the frequency (or TF-IDF score) of word i in document j.

• Problem: Synonymy (different words with identical meaning) and polysemy (one word with multiple meanings) create noise and high dimensionality.
• Truncated SVD Solution: A_k = U_k · Σ_k · V_k^T (typically k = 100 to 300).

- Columns of U_k represent latent semantic concepts.

- Columns of V_k represent document embeddings in concept space.

- Documents discussing similar concepts cluster together even if they share zero overlapping vocabulary!

In collaborative filtering (e.g., Netflix Prize, Amazon recommendations), user ratings on items form an immense, highly sparse matrix R in R^(users × items). Low-rank matrix completion via SVD decomposes R ≈ P · Q^T, discovering latent user preferences and latent item characteristics to predict unobserved ratings.

3. Low-Rank Adaptation (LoRA) for Large Language Models#

Modern Generative AI models (GPT-4, Claude, Llama 3) contain tens to hundreds of billions of parameters. Full parameter fine-tuning of a 70B model requires updating and storing dense gradient matrices W in R^(d×k) (where d, k ≈ 4096 to 8192), requiring hundreds of gigabytes of GPU VRAM.

The LoRA Paradigm:

Empirical studies grounded in SVD demonstrate that weight update matrices ΔW have a very low intrinsic rank.

Instead of updating the full dense matrix W = W_0 + ΔW, LoRA freezes the pretrained weights W_0 and parameterizes the update as a low-rank factorization:

ΔW = B · A

where:

• B in R^(d×r) is initialized to all zeros.
• A in R^(r×k) is initialized with Gaussian random weights.
• The rank r is chosen to be extremely small (e.g., r = 8 or 16).

Parameter Economics:

• Dense update parameters: d × k = 4096 × 4096 = 16,777,216 parameters per layer.
• LoRA parameters (r = 8): (4096 × 8) + (8 × 4096) = 32,768 + 32,768 = 65,536 parameters.
• Parameter reduction: 65,536 / 16,777,216 ≈ 0.39% (a 99.6% reduction in trainable parameters!).

Because SVD guarantees that a low-rank subspace captures the dominant update geometry, LoRA achieves performance identical to full fine-tuning while running on commodity consumer GPUs.

6. Common Architectural Pitfalls & Misconceptions#

• Pitfall 1: Confusing Singular Values with Eigenvalues.

Singular values are ALWAYS non-negative real numbers (σ_i >= 0), representing the lengths of geometric ellipsoid axes. Eigenvalues can be negative, zero, or complex conjugate pairs. For an arbitrary matrix A, singular values are the square roots of the eigenvalues of A^T · A, NOT the eigenvalues of A itself.

• Pitfall 2: Omitting Mean-Centering Prior to SVD in PCA.

If you compute SVD on raw uncentered data X, the first right singular vector v_1 will point toward the global mean vector of the data, rather than the direction of maximum variance. Always mean-center features before applying SVD for dimensionality reduction.

• Pitfall 3: Unfiltered Pseudoinversion of Ill-Conditioned Matrices.

Blindly executing x = pinv(A) @ b in numerical code without setting a threshold parameter rcond allows microscopic singular values (e.g., σ = 10^(-15)) to be inverted into 10^(15), destroying calculation stability. Always configure an explicit noise cutoff threshold τ.

• Pitfall 4: Assuming Left and Right Singular Vectors Have Unique Signs.

The signs of singular vectors are not uniquely determined. For any singular triplet (σ_i, u_i, v_i), simultaneously multiplying u_i and v_i by -1 produces the exact same rank-one matrix:

(-u_i) · (-v_i)^T = u_i · v_i^T.

Numerical software libraries across different programming languages (NumPy, MATLAB, LAPACK) may choose opposite sign conventions for singular vector columns.

7. Summary & Checkpoint Exercises#

Solidify your comprehension of Singular Value Decomposition, the Eckart-Young-Mirsky Theorem, and Moore-Penrose pseudoinverses with these four checkpoint problems. Full step-by-step mathematical solutions are provided.

Problem 1 (Full Manual SVD Calculation for a 3×2 Matrix):

Compute the complete thin SVD factorization A = U · Σ · V^T for the 3×2 matrix:

A = [

[ 3, 2 ],

[ 2, 3 ],

[ 2, -2 ]

]

1. Compute the Gram matrix S = A^T · A.
2. Find the eigenvalues and orthonormal right singular vectors v1, v2 (forming matrix V).
3. Determine the singular values σ1, σ2 (forming matrix Σ).
4. Compute the orthonormal left singular vectors u1, u2 (forming matrix U).
5. Verify that U · Σ · V^T exactly equals matrix A.

Solution:

Step 1: Compute S = A^T · A.

A^T = [

[ 3, 2, 2 ],

[ 2, 3, -2 ]

]

A^T · A = [

[ 3, 2, 2 ],

[ 2, 3, -2 ]

] · [

[ 3, 2 ],

[ 2, 3 ],

[ 2, -2 ]

]

• Entry (1, 1): (3)(3) + (2)(2) + (2)(2) = 9 + 4 + 4 = 17
• Entry (1, 2): (3)(2) + (2)(3) + (2)(-2) = 6 + 6 - 4 = 8
• Entry (2, 1): (2)(3) + (3)(2) + (-2)(2) = 6 + 6 - 4 = 8
• Entry (2, 2): (2)(2) + (3)(3) + (-2)(-2) = 4 + 9 + 4 = 17

A^T · A = [

[ 17, 8 ],

[ 8, 17 ]

]

Step 2: Eigenvalues and Right Singular Vectors (Matrix V).

Characteristic equation:

det(A^T · A - λ · I) = (17 - λ)^2 - 8^2 = 0

(17 - λ)^2 = 64 => 17 - λ = ± 8

λ1 = 17 + 8 = 25

λ2 = 17 - 8 = 9

Eigenvectors of A^T · A:

• For λ1 = 25:

[ [ 17 - 25, 8 ], [ 8, 17 - 25 ] ] · v = [ [ -8, 8 ], [ 8, -8 ] ] · v = 0 => -8 v1 + 8 v2 = 0 => v1 = v2.

Unit vector: v1 = [ 1/sqrt(2), 1/sqrt(2) ]^T.

• For λ2 = 9:

[ [ 17 - 9, 8 ], [ 8, 17 - 9 ] ] · v = [ [ 8, 8 ], [ 8, 8 ] ] · v = 0 => 8 v1 + 8 v2 = 0 => v1 = -v2.

Unit vector: v2 = [ -1/sqrt(2), 1/sqrt(2) ]^T.

Matrix V:

V = [

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

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

]

Step 3: Singular Values (Matrix Σ).

σ1 = sqrt(λ1) = sqrt(25) = 5

σ2 = sqrt(λ2) = sqrt(9) = 3

Matrix Σ (thin 2×2):

Σ = [

[ 5, 0 ],

[ 0, 3 ]

]

Step 4: Left Singular Vectors (Matrix U).

Use the relationship u_i = (1 / σ_i) · A · v_i:

• Compute u1:

A · v1 = [

[ 3, 2 ],

[ 2, 3 ],

[ 2, -2 ]

] · [

1/sqrt(2),

1/sqrt(2)

] = [

(3 + 2)/sqrt(2),

(2 + 3)/sqrt(2),

(2 - 2)/sqrt(2)

] = [

5/sqrt(2),

5/sqrt(2),

0

]

u1 = (1 / 5) · [ 5/sqrt(2), 5/sqrt(2), 0 ]^T = [ 1/sqrt(2), 1/sqrt(2), 0 ]^T.

• Compute u2:

A · v2 = [

[ 3, 2 ],

[ 2, 3 ],

[ 2, -2 ]

] · [

-1/sqrt(2),

1/sqrt(2)

] = [

(-3 + 2)/sqrt(2),

(-2 + 3)/sqrt(2),

(-2 - 2)/sqrt(2)

] = [

-1/sqrt(2),

1/sqrt(2),

-4/sqrt(2)

]

Wait, let's re-verify row 3:

v2 = [ -1/sqrt(2), 1/sqrt(2) ]^T.

Row 3: (2)(-1/sqrt(2)) + (-2)(1/sqrt(2)) = (-2 - 2)/sqrt(2) = -4/sqrt(2).

Wait, let's check norm of A · v2:

||A · v2||^2 = (-1/sqrt(2))^2 + (1/sqrt(2))^2 + (-4/sqrt(2))^2 = 1/2 + 1/2 + 16/2 = 18/2 = 9.

Norm = sqrt(9) = 3 = σ2. Exactly matches!

u2 = (1 / 3) · [ -1/sqrt(2), 1/sqrt(2), -4/sqrt(2) ]^T = [ -1/(3·sqrt(2)), 1/(3·sqrt(2)), -4/(3·sqrt(2)) ]^T.

Matrix U (thin 3×2):

U = [

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

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

[ 0, -4/(3·sqrt(2)) ]

]

Step 5: Verification of U · Σ · V^T = A.

Compute U · Σ:

U · Σ = [

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

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

[ 0, -4/(3·sqrt(2)) ]

] · [

[ 5, 0 ],

[ 0, 3 ]

] = [

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

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

[ 0, -4/sqrt(2) ]

]

Now multiply by V^T = [ [ 1/sqrt(2), 1/sqrt(2) ], [ -1/sqrt(2), 1/sqrt(2) ] ]:

• Row 1:

Col 1: (5/sqrt(2))(1/sqrt(2)) + (-1/sqrt(2))(-1/sqrt(2)) = 5/2 + 1/2 = 6/2 = 3.

Col 2: (5/sqrt(2))(1/sqrt(2)) + (-1/sqrt(2))(1/sqrt(2)) = 5/2 - 1/2 = 4/2 = 2.

• Row 2:

Col 1: (5/sqrt(2))(1/sqrt(2)) + (1/sqrt(2))(-1/sqrt(2)) = 5/2 - 1/2 = 4/2 = 2.

Col 2: (5/sqrt(2))(1/sqrt(2)) + (1/sqrt(2))(1/sqrt(2)) = 5/2 + 1/2 = 6/2 = 3.

• Row 3:

Col 1: (0)(1/sqrt(2)) + (-4/sqrt(2))(-1/sqrt(2)) = 0 + 4/2 = 2.

Col 2: (0)(1/sqrt(2)) + (-4/sqrt(2))(1/sqrt(2)) = 0 - 4/2 = -2.

The product is:

[

[ 3, 2 ],

[ 2, 3 ],

[ 2, -2 ]

] = A.

The thin SVD factorization is mathematically exact.

Problem 2 (Eckart-Young-Mirsky Rank-1 Approximation and Error Analysis):

Using the SVD factorization from Problem 1:

σ1 = 5, u1 = [ 1/sqrt(2), 1/sqrt(2), 0 ]^T, v1 = [ 1/sqrt(2), 1/sqrt(2) ]^T

σ2 = 3, u2 = [ -1/(3·sqrt(2)), 1/(3·sqrt(2)), -4/(3·sqrt(2)) ]^T, v2 = [ -1/sqrt(2), 1/sqrt(2) ]^T

1. Construct the optimal rank-1 approximation matrix A_1 = σ1 · u1 · v1^T.
2. Calculate the difference matrix E = A - A_1.
3. Verify the Eckart-Young-Mirsky Theorem by calculating the Spectral norm error ||A - A_1||_2 and Frobenius norm error ||A - A_1||_F.

Solution:

Part 1: Construct A_1.

u1 · v1^T = [

1/sqrt(2),

1/sqrt(2),

0

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

= [

[ 1/2, 1/2 ],

[ 1/2, 1/2 ],

[ 0, 0 ]

]

Scale by σ1 = 5:

A_1 = 5 · [

[ 1/2, 1/2 ],

[ 1/2, 1/2 ],

[ 0, 0 ]

] = [

[ 2.5, 2.5 ],

[ 2.5, 2.5 ],

[ 0.0, 0.0 ]

]

Part 2: Error Matrix E = A - A_1.

E = [

[ 3, 2 ],

[ 2, 3 ],

[ 2, -2 ]

] - [

[ 2.5, 2.5 ],

[ 2.5, 2.5 ],

[ 0.0, 0.0 ]

]

E = [

[ 0.5, -0.5 ],

[ -0.5, 0.5 ],

[ 2.0, -2.0 ]

]

Part 3: Verify Norm Errors.

By SVD expansion: E = A - A_1 = σ2 · u2 · v2^T.

• Spectral Norm Error ||A - A_1||_2:

||E||_2 = ||σ2 · u2 · v2^T||_2 = σ2 · ||u2||_2 · ||v2||_2 = 3 · (1) · (1) = 3.

Theoretical error: σ_{k+1} = σ_2 = 3.

The spectral norm matches the theoretical bound exactly.

• Frobenius Norm Error ||A - A_1||_F:

Sum of squared entries of E:

(0.5)^2 + (-0.5)^2 + (-0.5)^2 + (0.5)^2 + (2.0)^2 + (-2.0)^2

= 0.25 + 0.25 + 0.25 + 0.25 + 4.0 + 4.0 = 1.0 + 8.0 = 9.0.

||E||_F = sqrt(9.0) = 3.0.

Theoretical error: sqrt(sum(i=k+1 to r, σ_i^2)) = sqrt(σ_2^2) = sqrt(3^2) = 3.

The Frobenius norm error matches the theoretical bound perfectly.

Problem 3 (Moore-Penrose Pseudoinverse and Minimum-Norm Underdetermined Solution):

Consider the underdetermined linear equation with 1 constraint and 2 variables:

[ 1, 2 ] · [ x1, x2 ]^T = [ 5 ]

Here matrix A = [ 1, 2 ] has dimension 1×2 and rank r = 1.

1. Compute the Moore-Penrose pseudoinverse A^+ using the full-rank row formula A^+ = A^T · (A · A^T)^(-1).
2. Compute the pseudoinverse solution x_hat = A^+ · [ 5 ].
3. Prove that among all possible solutions to x1 + 2 x2 = 5, x_hat achieves the minimum Euclidean norm ||x||_2.

Solution:

Part 1: Compute A^+.

Matrix A = [ 1, 2 ].

A^T = [ 1, 2 ]^T.

Compute A · A^T:

A · A^T = [ 1, 2 ] · [ 1, 2 ]^T = (1)(1) + (2)(2) = 1 + 4 = 5 (a 1×1 scalar).

(A · A^T)^(-1) = 1/5.

Compute A^+:

A^+ = A^T · (A · A^T)^(-1) = [

1,

2

] · (1/5) = [

1/5,

2/5

] = [

0.2,

0.4

]

Part 2: Compute x_hat = A^+ · b.

b = [ 5 ].

x_hat = [ 0.2, 0.4 ]^T · [ 5 ] = [ 0.2 · 5, 0.4 · 5 ]^T = [ 1, 2 ]^T.

Verify constraint: x1 + 2 x2 = 1 + 2(2) = 1 + 4 = 5. Satisfied!

Part 3: Minimum Norm Proof.

Any general solution to x1 + 2 x2 = 5 can be parameterized by setting x2 = t:

x1 = 5 - 2t.

General solution vector:

x(t) = [ 5 - 2t, t ]^T.

Evaluate the Euclidean norm squared:

f(t) = ||x(t)||^2 = (5 - 2t)^2 + t^2 = 25 - 20t + 4t^2 + t^2 = 5t^2 - 20t + 25.

To find the minimum, set derivative df/dt = 0:

df/dt = 10t - 20 = 0 => 10t = 20 => t = 2.

At t = 2:

x1 = 5 - 2(2) = 1.

x2 = 2.

x_optimal = [ 1, 2 ]^T.

The minimum Euclidean norm is ||x_optimal|| = sqrt(1^2 + 2^2) = sqrt(5) ≈ 2.236.

The Moore-Penrose pseudoinverse solution x_hat = [ 1, 2 ]^T exactly identifies the unique point on the constraint line closest to the origin.

Problem 4 (Truncated SVD Denoising for an Ill-Conditioned Sensor System):

A satellite sensor measures two physical parameters x = [ x1, x2 ]^T through a coupled linear system A · x = b_noisy.

The matrix A has SVD:

U = [ [ 1, 0 ], [ 0, 1 ] ] = I_2

V = [ [ 1, 0 ], [ 0, 1 ] ] = I_2

Σ = [ [ 10.0, 0.0 ], [ 0.0, 0.001 ] ]

The true state is x_true = [ 2.0, 1.0 ]^T, yielding clean measurement b_clean = A · x_true = [ 20.0, 0.001 ]^T.

Suppose sensor noise perturbs the measurement to b_noisy = [ 20.0, 0.005 ]^T (noise = [ 0, 0.004 ]^T).

1. Compute the unregularized pseudoinverse solution x_unreg = A^(-1) · b_noisy. Show how noise is amplified.
2. Apply Truncated SVD (TSVD) with threshold τ = 0.01 to compute x_tsvd.
3. Compare the estimation error ||x - x_true||_2 for both solutions.

Solution:

Part 1: Unregularized Solution.

Since U = V = I, A = Σ = [ [ 10.0, 0 ], [ 0, 0.001 ] ].

The inverse is:

A^(-1) = [ [ 1/10.0, 0 ], [ 0, 1/0.001 ] ] = [ [ 0.1, 0 ], [ 0, 1000.0 ] ].

Calculate x_unreg = A^(-1) · b_noisy:

x_unreg = [

[ 0.1, 0 ],

[ 0, 1000.0 ]

] · [

20.0,

0.005

] = [

0.1 · 20.0,

1000.0 · 0.005

] = [

2.0,

5.0

]

Analysis of Noise Amplification:

The true x2 was 1.0. Because σ2 = 0.001, the measurement noise of 0.004 was magnified by 1/0.001 = 1000, creating an error of 4.0 in parameter x2!

Estimation error:

||x_unreg - x_true||_2 = || [ 2.0 - 2.0, 5.0 - 1.0 ]^T ||_2 = || [ 0, 4.0 ]^T ||_2 = 4.0.

Part 2: Truncated SVD Solution (TSVD).

Threshold τ = 0.01.

Check singular values against threshold:

• σ1 = 10.0 >= 0.01: Retain. Invert: 1 / 10.0 = 0.1.
• σ2 = 0.001 < 0.01: Discard. Set 1 / σ2 = 0.

The regularized diagonal matrix is:

Σ_tsvd^+ = [

[ 0.1, 0.0 ],

[ 0.0, 0.0 ]

]

Compute x_tsvd = V · Σ_tsvd^+ · U^T · b_noisy:

x_tsvd = [

[ 0.1, 0.0 ],

[ 0.0, 0.0 ]

] · [

20.0,

0.005

] = [

2.0,

0.0

]

Part 3: Error Comparison.

For TSVD:

||x_tsvd - x_true||_2 = || [ 2.0 - 2.0, 0.0 - 1.0 ]^T ||_2 = || [ 0, -1.0 ]^T ||_2 = 1.0.

Comparison:

• Unregularized error: 4.0
• TSVD regularized error: 1.0

Truncated SVD prevented the noise from blowing up by a factor of 1000, reducing the total estimation error by 75% and stabilizing the physical parameters.