Principal Component Analysis, Part 2

From mathematical foundations to practical implementation

Principal Component Analysis (PCA) is a fundamental dimensionality reduction technique that identifies an optimal $r$-dimensional subspace capturing the maximum variance in the data. Its mathematical foundation rests on finding orthogonal directions — called principal components — that sequentially maximize the projected variance.

Given a centered data matrix $\mathbf{X}$, we seek unit vectors $\mathbf{u}$ that maximize the variance of projected data $\mathbf{X}\mathbf{u}$:

This optimization leads to the key eigenvalue problem:

where $\mathbf{C}$ is a covariance matrix.

In the following sections, we will systematically derive the solution to this problem and demonstrate its geometric interpretation.

Linear Algebra Foundations¶

Data Matrix Representation¶

A data matrix $\mathbf{X}$ with $n$ rows (samples) and $d$ columns (features) can be written as:

Each row $\mathbf{x}_i = (x_{i1}, x_{i2}, \ldots, x_{id})^T$ represents one sample in $d$-dimensional space.

Vector Operations¶

Dot Product¶

Vector Length (Norm)¶

Unit Vector¶

Distance between Vectors¶

Angle between Vectors¶

Vector Projections (Derivation & Geometric Intuition)¶

Derivation of Projection Formula¶

We begin with fundamental vector operations to derive the projection formula naturally:

1. Dot Product and Angle Relationship:
   From our previous definitions, we know:

   $$\mathbf{a}^T \mathbf{b} = \|\mathbf{a}\| \|\mathbf{b}\| \cos \theta$$

   (9)

   This measures alignment between vectors.

2. Scalar Projection (Length of Shadow):
   The length of $\mathbf{b}$'s shadow along $\mathbf{a}$ is:

   $$\text{Scalar projection} = \|\mathbf{b}\| \cos \theta = \frac{\mathbf{a}^T \mathbf{b}}{\|\mathbf{a}\|}$$

   (10)

   For unit vectors, this simplifies to $\mathbf{a}^T \mathbf{b}$.

3. Vector Projection:
   To get the vector form, we multiply the scalar projection by $\mathbf{a}$'s unit vector:

   $$\mathbf{b}_{\parallel} = \left( \frac{\mathbf{a}^T \mathbf{b}}{\|\mathbf{a}\|} \right) \left( \frac{\mathbf{a}}{\|\mathbf{a}\|} \right) = \frac{\mathbf{a}^T \mathbf{b}}{\|\mathbf{a}\|^2} \mathbf{a}$$

   (11)

Final Projection Formula¶

Key Insights:

1. Dot Product Interpretation:

   • The numerator $\mathbf{a}^T \mathbf{b}$ measures alignment

   • The denominator $\|\mathbf{a}\|^2$ normalizes for vector length

2. Special Case - Unit Vector: When $\|\mathbf{a}\|=1$:

   $$\mathbf{b}_{\parallel} = (\mathbf{a}^T \mathbf{b}) \mathbf{a}$$

   (13)

3. Geometric Meaning:

   • $\mathbf{b}_{\parallel}$ is $\mathbf{b}$'s shadow on $\mathbf{a}$

   • The residual $\mathbf{b}_{\perp} = \mathbf{b} - \mathbf{b}_{\parallel}$ is orthogonal to $\mathbf{a}$

Example

Let $\mathbf{a} = \begin{bmatrix} 1 \\ 0 \end{bmatrix}$ and $\mathbf{b} = \begin{bmatrix} 3 \\ 2 \end{bmatrix}$.

• Step 1: Compute $\mathbf{a}^T \mathbf{b} = 1 \times 3 + 0 \times 2 = 3$.

• Step 2: Compute $\|\mathbf{a}\|^2 = 1^2 + 0^2 = 1$.

• Step 3: Projection $\mathbf{b}_{\parallel} = 3\begin{bmatrix} 1 \\ 0 \end{bmatrix} = \begin{bmatrix} 3 \\ 0 \end{bmatrix}$.

• Perpendicular Component: $\mathbf{b}_{\perp} = \begin{bmatrix} 3 \\ 2 \end{bmatrix} - \begin{bmatrix} 3 \\ 0 \end{bmatrix} = \begin{bmatrix} 0 \\ 2 \end{bmatrix}$.

Interpretation:
Since $\mathbf{a}$ points along the x-axis, $\mathbf{b}_{\parallel}$ retains only the x-component of $\mathbf{b}$, and $\mathbf{b}_{\perp}$ captures the remaining y-component, which is orthogonal to $\mathbf{a}$.

Connection to PCA:
This projection concept is the geometric heart of PCA. In the next section we extend it to multiple dimensions, leading to the PCA solution via eigenvalue decomposition.

Basis Transformation¶

Change of Basis¶

Any vector $\mathbf{x}$ can be expressed in terms of a new orthonormal basis $\{\mathbf{u}_1, \mathbf{u}_2, \ldots, \mathbf{u}_d\}$:

where $\mathbf{U} = [\mathbf{u}_1 \mid \mathbf{u}_2 \mid \cdots \mid \mathbf{u}_d]$ has the new basis vectors as its columns, and $\mathbf{a}$ holds the coordinates of $\mathbf{x}$ in the new basis:

This works because orthonormality gives $\mathbf{U}^T\mathbf{U} = \mathbf{I}$, so the coordinate in direction $\mathbf{u}_i$ is simply the dot product $a_i = \mathbf{u}_i^T \mathbf{x}$.

Rotation Matrix Mathematical Foundation¶

2D Rotation Fundamentals¶

A rotation matrix transforms vectors by rotating them through a specified angle θ while preserving their length. For a 45° rotation in the xy-plane (θ = π/4 radians), the fundamental rotation matrix is:

For θ = 45°:

Thus, the 45° rotation matrix becomes:

Basis Transformation¶

In the previous code, this rotation matrix defines new basis vectors:

u1 = np.array([1/np.sqrt(2), 1/np.sqrt(2), 0])  # Rotated x-axis
u2 = np.array([-1/np.sqrt(2), 1/np.sqrt(2), 0]) # Rotated y-axis 
u3 = np.array([0, 0, 1])                        # Unchanged z-axis

These form an orthonormal basis:

1. $\mathbf{u}_1$: Original x-axis rotated 45° counterclockwise

2. $\mathbf{u}_2$: Original y-axis rotated 45° counterclockwise

3. $\mathbf{u}_3$: Original z-axis (unchanged)

Key Properties¶

1. Orthonormality:

   • $\mathbf{u}_i^T\mathbf{u}_j = \begin{cases} 1 & \text{if } i=j \\ 0 & \text{otherwise} \end{cases}$

2. Determinant: $\det(\mathbf{R}) = 1$ (preserves orientation)

3. Inverse: $\mathbf{R}^{-1} = \mathbf{R}^T$ (transpose reverses rotation)

Vector Transformation¶

For any vector $\mathbf{x} = [x,y,z]^T$, the rotated version is:

This matches our basis transformation approach where:

Why This Works¶

The columns of $\mathbf{R}$ are exactly the images of the standard basis vectors under the rotation. When we multiply:

So $\mathbf{R}$ tells us precisely where each axis goes — which is why its columns form the new basis.

PCA Mathematical Foundation¶

Variance Maximization Problem¶

PCA seeks an orthonormal basis $\{\mathbf{u}_1, \dots, \mathbf{u}_r\}$ that captures maximal variance in the data through sequential optimization:

1. First Principal Component:
   Find unit vector $\mathbf{u}_1$ maximizing projected variance:

   $$\max_{\|\mathbf{u}_1\|=1} \text{Var}(\mathbf{X}\mathbf{u}_1)$$

   (22)

2. Subsequent Components:
   Each $\mathbf{u}_k$ is found by maximizing variance orthogonal to previous components:

   $$\max_{\|\mathbf{u}_k\|=1, \, \mathbf{u}_k \perp \{\mathbf{u}_1, \dots, \mathbf{u}_{k-1}\}} \text{Var}(\mathbf{X}\mathbf{u}_k)$$

   (23)

In short, the first principal component is the direction of greatest spread in the data. Each subsequent component must be orthogonal to all previous ones, ensuring it captures genuinely new, independent information.

Finding the first principal component¶

We seek a unit vector $\mathbf{u}$ that maximizes the variance of the projected data:

where $\mu_{\mathbf{u}} = \frac{1}{n} \sum_{i=1}^n \mathbf{u}^T \mathbf{x}_i = \mathbf{u}^T \bar{\mathbf{x}}$ is the mean of the projected values (using linearity of the dot product).

Expanding the variance:

1. Let $a_i = \mathbf{u}^T \mathbf{x}_i$ denote the scalar projection of the $i$-th sample.

2. The variance becomes:

   $$\sigma^2_{\mathbf{u}} = \frac{1}{n} \sum_{i=1}^n (a_i - \mu_{\mathbf{u}})^2 = \frac{1}{n} \sum_{i=1}^n \left( \mathbf{u}^T \mathbf{x}_i - \mathbf{u}^T \bar{\mathbf{x}} \right)^2 = \frac{1}{n} \sum_{i=1}^n \left[ \mathbf{u}^T (\mathbf{x}_i - \bar{\mathbf{x}}) \right]^2$$

   (25)

Let $\mathbf{z}_i = \mathbf{x}_i - \bar{\mathbf{x}}$ denote the centered data points. Expanding step by step:

1. Use the identity $(A^TB)^T = B^TA$, noting that $(\mathbf{u}^T \mathbf{z}_i)^2 = (\mathbf{u}^T \mathbf{z}_i)(\mathbf{z}_i^T \mathbf{u})$:

2. Rearrange, grouping $\mathbf{z}_i\mathbf{z}_i^T$ between $\mathbf{u}^T$ and $\mathbf{u}$:

3. Pull $\mathbf{u}^T$ and $\mathbf{u}$ outside the sum:

4. Recognize the sample covariance matrix $\mathbf{C} = \frac{1}{n} \sum_{i=1}^n \mathbf{z}_i \mathbf{z}_i^T$:

5. Equivalently, using the centered data matrix $\mathbf{Z} = [\mathbf{z}_1 \cdots \mathbf{z}_n]^T \in \mathbb{R}^{n \times d}$:

Hence the variance of the projected data is:

Constrained Optimization via Lagrange Multipliers¶

We maximize the projected variance $\sigma^2_\mathbf{u} = \mathbf{u}^T \mathbf{C} \mathbf{u}$ subject to $\|\mathbf{u}\|=1$. Introducing Lagrange multiplier $\lambda$, the Lagrangian is:

Solution¶

Setting the gradient with respect to $\mathbf{u}$ to zero:

This is precisely the eigenvalue equation: $\mathbf{u}$ must be an eigenvector of the covariance matrix $\mathbf{C}$, with $\lambda$ the corresponding eigenvalue.

Choosing the Right Eigenvalue¶

Substituting $\mathbf{C}\mathbf{u} = \lambda\mathbf{u}$ back into the variance expression:

The projected variance equals the eigenvalue $\lambda$. To maximize variance we therefore pick $\lambda_1$, the largest eigenvalue of $\mathbf{C}$. The corresponding eigenvector $\mathbf{u}_1$ defines the direction of maximum variance — the first principal component.

PCA Algorithm¶

Given a data matrix $\mathbf{X} \in \mathbb{R}^{n \times d}$:

1. Center the data: $\mathbf{Z} = \mathbf{X} - \mathbf{1}\cdot\mathbf{\mu}^T$

2. Compute the covariance matrix: $\mathbf{C} = \frac{1}{n}\mathbf{Z}^T\mathbf{Z}$

3. Compute eigenvalues and eigenvectors of $\mathbf{C}$

4. Sort eigenvalues in descending order

5. Select the top $r$ eigenvectors as the principal components

6. Project the data: $\mathbf{Z}' = \mathbf{Z}\mathbf{U}_r$

PCA Workflow: From Single Sample to Full Matrix¶

Projection & Reconstruction (Single Data Point)¶

Step 1: Mean-Centering¶

• Given a raw data vector: $\mathbf{x} = [x_1, x_2, \dots, x_d]^T \quad \text{(original space)}$

• Subtract the mean vector: $\mathbf{z} = \mathbf{x} - \mathbf{\mu} \quad \text{(centered)}$

  where $\mathbf{\mu}$ is the mean of all samples.

Step 2: Projection (Dimensionality Reduction)¶

• Project $\mathbf{z}$ onto the first $r$ principal components (PCs) $\mathbf{u}_1, \dots, \mathbf{u}_r$: $a_i = \mathbf{u}_i^T \mathbf{z}$ (scalar, projection of $\mathbf{z}$ onto $\mathbf{u}_i$)

• The reduced representation (in PCA space) is: $\mathbf{a} = [a_1, a_2, \dots, a_r]^T$

Step 3: Reconstruction (Approximation)¶

• Reconstruct the centered data using the PCs: $\mathbf{z}' = \sum_{i=1}^r a_i \mathbf{u}_i = a_1 \mathbf{u}_1 + a_2 \mathbf{u}_2 + \dots + a_r \mathbf{u}_r$

• Add back the mean to return to original space: $\mathbf{x}_{\text{recon}} = \mathbf{z}' + \mathbf{\mu}$

Geometric Interpretation¶

• The reconstruction $\mathbf{z}'$ is the closest approximation of $\mathbf{z}$ using only $r$ directions.

• The error $\|\mathbf{z} - \mathbf{z}'\|$ comes from discarding the remaining PCs.

Generalization to Full Data Matrix¶

Step 1: Mean-Centering (Matrix Form)¶

• Original data matrix:

  $$\mathbf{X} = [\mathbf{x}_1 | \mathbf{x}_2 | \dots | \mathbf{x}_n]^T \quad \text{($n \times d$)}$$

  (37)

• Mean-centered data:

  $$\mathbf{Z} = \mathbf{X} - \mathbf{1} \mathbf{\mu}^T$$

  (38)

  where $\mathbf{1}$ is a column vector of ones.

Step 2: Projection (PCA Transformation)¶

• Project all samples onto the first $r$ PCs ($\mathbf{U}_r = [\mathbf{u}_1 | \dots | \mathbf{u}_r]$):

  $$\mathbf{A} = \mathbf{Z} \mathbf{U}_r \quad \text{($n \times r$ matrix of PCA coordinates)}$$

  (39)

  • Each row of $\mathbf{A}$ contains the reduced representation of a sample.

Step 3: Reconstruction (Matrix Form)¶

• Reconstruct the centered data:

  $$\mathbf{Z}' = \mathbf{A} \mathbf{U}_r^T = \sum_{i=1}^r \mathbf{a}_i \mathbf{u}_i^T$$

  (40)

  • $\mathbf{a}_i$ = $i$-th column of $\mathbf{A}$ (scores for PC $i$).

  • Each term $\mathbf{a}_i \mathbf{u}_i^T$ is a rank-1 matrix.

• Return to original space:

  $$\mathbf{X}_{\text{recon}} = \mathbf{Z}' + \mathbf{1} \mathbf{\mu}^T$$

  (41)

Connection to Python Implementation¶

Projection (Dimensionality Reduction)¶

X_transformed = Z @ U_r  # ≡ Z U_r (n × r matrix)

Reconstruction¶

Z_recon = X_transformed @ U_r.T  # ≡ A U_r^T (n × d)
X_recon = Z_recon + mean        # Add back mean

Key Takeaways¶

1. Single Sample $\rightarrow$ Full Matrix:

   • Start with a single vector to understand projection/reconstruction.

   • Generalize to the full dataset using matrix operations.

2. Lossy Compression:

   • Keeping only $r < d$ PCs means the reconstruction is approximate.

   • The error depends on the discarded eigenvalues.

3. PCA as a Rotation:

   • Projection = Rotate data to align with maximum variance directions.

   • Reconstruction = Rotate back, using only the most important directions.

Equivalence of PCA’s Two Optimization Perspectives and Error Analysis (Further Reading)¶

PCA can be derived from two equivalent perspectives:

1. Variance Maximization: Find directions that maximize projected variance.

2. Reconstruction Error Minimization: Find directions that minimize the mean squared reconstruction error.

We’ll prove their equivalence mathematically and analyze the error.

1. Variance Maximization Formulation¶

Goal: Find unit vector $\mathbf{u}$ that maximizes the variance of projected data:

• The solution is the eigenvector of $\mathbf{Z}^T \mathbf{Z}$ with largest eigenvalue.

2. Reconstruction Error Minimization¶

Goal: Find $\mathbf{u}$ that minimizes the MSE between original data $\mathbf{z}_j$ and its reconstruction $\mathbf{z}_j' = (\mathbf{u}^T \mathbf{z}_j)\mathbf{u}$:

3. Proof of Equivalence¶

Expand the reconstruction error:

Thus, the MSE becomes:

Key Observations:

1. The term $\frac{1}{n} \sum \|\mathbf{z}_j\|^2$ is constant (total variance in data).

2. Minimizing MSE is equivalent to maximizing $\frac{1}{n} \sum (\mathbf{u}^T \mathbf{z}_j)^2$ (projected variance).

4. Error Analysis in PCA¶

Projection and Residual Components¶

For a data point $\mathbf{z}_j$ (centered) and its reconstruction $\mathbf{z}_j'$ using $r$ PCs:

• Error vector: $\mathbf{\epsilon}_j = \mathbf{z}_j - \mathbf{z}_j' = \sum_{i=r+1}^d (\mathbf{u}_i^T \mathbf{z}_j)\mathbf{u}_i$.

Mean Squared Error (MSE) Derivation¶

The total MSE across all samples is:

Link to Eigenvalues¶

By definition, the variance along $\mathbf{u}_i$ is:

Thus:

5. Conclusion and Summary¶

• Dual Formulations: Variance maximization and error minimization are equivalent in PCA.

• Eigenvalue Interpretation: The MSE equals the sum of discarded eigenvalues because:

  • Each $\lambda_i$ represents variance along PC $\mathbf{u}_i$.

  • Discarding PCs $r+1$ to $d$ loses exactly $\sum_{i=r+1}^d \lambda_i$ variance.

• Optimality: PCA’s solution simultaneously maximizes preserved variance and minimizes reconstruction error.

Key Relationships:

Examples and Applications¶

Iris Dataset Example¶

Reconstruction Error Analysis¶

High-Dimensional Data Example¶

Comparison with Scikit-learn¶

Some of my previous papers are about PCA, including Amintoosi & Farbiz (2022)Fathy et al. (2008)Amintoosi et al. (2007)Amintoosi & Poursadeghi (1394)

Summary¶

This tutorial covered:

1. Linear Algebra Foundations: Vector operations, norms, and angles

2. Vector Projections: Mathematical foundation for dimensionality reduction

3. Basis Transformation: How to represent data in different coordinate systems

4. PCA Theory: Variance maximization and eigenvalue decomposition

5. Implementation: Step-by-step PCA algorithm from scratch

6. Applications: Real examples with Iris dataset and high-dimensional data

Key Takeaways:¶

• PCA finds directions of maximum variance in the data

• Principal components are eigenvectors of the covariance matrix

• The eigenvalues represent the amount of variance captured by each component

• PCA provides optimal linear dimensionality reduction in terms of reconstruction error

• The number of components can be chosen based on desired explained variance threshold