1. Question 1

Let $A=U\mathrm{\Sigma }{V}^T$ be the SVD of $A$, where $A\in {\mathbb{R}}^{m\times n},U\in {\mathbb{R}}^{m\times m}$ and $V\in {\mathbb{R}}^{n\times n}$ are orthogonal matrices, $\mathrm{\Sigma }=\mathrm{diag}({\sigma }_1,\cdots ,{\sigma }_r,0,\cdots ,0)$, and $r=\mathrm{rank}(A)$. Show that

1. The first $r$ columns of $U$ are eigenvectors of $A{A}^T$ corresponding to nonzero eigenvalues.
2. The first $r$ columns of $V$ are eigenvectors of $A^{T}A$ corresponding to nonzero eigenvalues.

1. The Background You Need to Know

To solve this, we need to understand a few fundamental linear algebra rules and definitions:

• The Transpose (${}^T$): Taking the transpose of a matrix simply means swapping its rows and columns. A crucial rule you need to know is that when you take the transpose of matrices multiplied together, you flip their order: $(ABC{)}^T=C^T{B}^T{A}^T$.

• Orthogonal Matrices ($U$ and $V$): The problem states that $U$ and $V$ are orthogonal matrices. An orthogonal matrix is a special square matrix where its columns (and rows) are perpendicular to each other and have a length of 1. The magic property of an orthogonal matrix $Q$ is that multiplying it by its transpose results in the Identity matrix $I$ (a matrix of all 1s on the diagonal and 0s elsewhere). Mathematically, this means $U^{T}U=U{U}^T=I$ and $V^{T}V=V{V}^T=I$. Multiplying by $I$ is like multiplying by 1 in regular math; it leaves the other matrix unchanged.

• Eigen-Decomposition: If you have a symmetric square matrix $S$, you can factor it into three parts: $S=Q\mathrm{\Lambda }{Q}^T$. In this specific format, the columns of the matrix $Q$ exactly represent the eigenvectors of $S$, and the diagonal elements of the middle matrix $\mathrm{\Lambda }$ represent the corresponding eigenvalues.

• The Rank ($r$): The rank $r$ of the matrix simply tells us how many of the singular values (${\sigma }_1,\dots ,{\sigma }_r$) in the diagonal matrix $\mathrm{\Sigma }$ are strictly positive (greater than 0). Any singular values after the $r$-th one are exactly zero.

  1. What exactly is the "Rank" (r)?
     • In linear algebra, the rank of a matrix is the maximum number of linearly independent rows (or equivalently the maximum number of linearly independent columns).
     • To put it simply, "linearly independent" means that a row (or column) contains truly unique information that cannot be created just by mixing or scaling the other rows. For example, if you have a dataset where column A is "Age in Years" and column B is "Age in Months", column B is completely redundant because it is just column A multiplied by 12. It does not add any new dimension of information, so it does not increase the rank.
     • Therefore, the rank (r) simply tells you exactly how many unique, non-redundant dimensions of information are hidden inside your matrix.
  2. How SVD acts as building blocks
     • When you perform SVD on a matrix $R$, you are decomposing it into three matrices: $U$, $\mathrm{\Sigma }$, and $V^T$.
     • A very intuitive way to look at this mathematical multiplication is that it breaks the original matrix down into a sequence of simpler, "rank-1" matrices summed together. The SVD formula can actually be rewritten as a sum of individual components: $R={\sigma }_1(u_1{v}_1^T)+{\sigma }_2(u_2{v}_2^T)+\dots$
     • In this view:
       • The terms like $(u_1{v}_1^T)$ are the structural building blocks (representing specific directions or patterns in your data).
       • The singular values ($\sigma$) on the diagonal of the $\mathrm{\Sigma }$ matrix act as the "scales" or "weights" for each of these building blocks.
  3. Why the first $r$ singular values are greater than zero
     • By mathematical convention, the elements in the $\mathrm{\Sigma }$ matrix are always non-negative real numbers, and we always arrange them in sorted order from largest to smallest (${\sigma }_1\ge {\sigma }_2\ge \cdots \ge 0$).
     • Here is where the two concepts perfectly connect:
       • If your original matrix has a rank of $r$, it means it only possesses $r$ unique dimensions of information.
       • Therefore, SVD only needs $r$ building blocks to perfectly reconstruct your matrix.
       • To "turn on" these $r$ necessary building blocks, the first $r$ singular values (${\sigma }_1,{\sigma }_2,\dots ,{\sigma }_r$) must be strictly greater than zero.
       • Because the matrix has no more unique information left to capture beyond those $r$ dimensions, any remaining building blocks are completely useless. To "turn them off", all subsequent singular values (${\sigma }_{r+1},{\sigma }_{r+2},\dots ,{\sigma }_n$) are set to exactly 0.
     • In fact, the relationship is so direct that the number of positive singular values directly tells you the rank of that matrix.
     • Because everything after the $r$-th singular value is multiplied by zero, they contribute absolutely nothing to the final matrix. This is why in a technique called "Compact SVD", we can literally just throw away all the rows and columns corresponding to the 0 elements without losing a single piece of data.

2. Solving Part 1: Proving $U$ contains the eigenvectors of $A{A}^T$

We are asked to look at the matrix $A{A}^T$. Since we know the SVD of $A$ is $A=U\mathrm{\Sigma }{V}^T$, we can substitute this definition into $A{A}^T$.

Step 1: Substitute $A$ with its SVD formulation. $$A A^T = (U \Sigma V^T)(U \Sigma V^T)^T$$

Step 2: Apply the transpose rule to the second part. Using our rule $(ABC{)}^T=C^T{B}^T{A}^T$, the second set of parentheses becomes $(V^T{)}^T{\mathrm{\Sigma }}^T{U}^T$. Since the transpose of a transpose returns the original matrix, $(V^T{)}^T=V$. Also, because $\mathrm{\Sigma }$ is a diagonal matrix, transposing it does nothing, so ${\mathrm{\Sigma }}^T=\mathrm{\Sigma }$. Therefore, $(U\mathrm{\Sigma }{V}^T{)}^T=V\mathrm{\Sigma }{U}^T$.

Step 3: Put it back together. $$A A^T = U \Sigma V^T V \Sigma U^T$$.

Step 4: Use the orthogonal matrix magic. Notice that $V^{T}V$ sits right in the middle. Because $V$ is an orthogonal matrix, we know that $V^{T}V=I$. $$A A^T = U \Sigma (I) \Sigma U^T$$ Since $I$ just acts like a 1, we can drop it. Then we multiply the two $\mathrm{\Sigma }$ matrices together. Since they are diagonal matrices, multiplying them just squares the values on the diagonal. $$A A^T = U \Sigma^2 U^T$$.

Step 5: The Conclusion. Look at our final equation: $A{A}^T=U{\mathrm{\Sigma }}^2{U}^T$. This perfectly matches the eigen-decomposition format $S=Q\mathrm{\Lambda }{Q}^T$. By matching the pieces, we can definitively say that the columns of $U$ are the eigenvectors of $A{A}^T$. Furthermore, the diagonal elements of ${\mathrm{\Sigma }}^2$ are its eigenvalues (which are ${\sigma }_1^2,{\sigma }_2^2,\dots$). Because the rank is $r$, only the first $r$ singular values (${\sigma }_1\dots {\sigma }_r$) are greater than zero. Therefore, the first $r$ columns of $U$ correspond to the strictly nonzero eigenvalues!

3. Solving Part 2: Proving $V$ contains the eigenvectors of $A^{T}A$

This part is almost perfectly identical to Part 1, but in reverse. We will look at $A^{T}A$ instead.

Step 1: Substitute $A$ into the equation. $$A^T A = (U \Sigma V^T)^T (U \Sigma V^T)$$

Step 2: Apply the transpose rule to the first part. Just like before, $(U\mathrm{\Sigma }{V}^T{)}^T$ becomes $V\mathrm{\Sigma }{U}^T$. $$A^T A = V \Sigma U^T U \Sigma V^T$$.

Step 3: Use the orthogonal matrix magic. This time, we have $U^{T}U$ in the middle. Because $U$ is an orthogonal matrix, $U^{T}U=I$. $$A^T A = V \Sigma (I) \Sigma V^T$$ Squaring the diagonal matrices gives us our final form: $$A^T A = V \Sigma^2 V^T$$.

Step 4: The Conclusion. Once again, $V{\mathrm{\Sigma }}^2{V}^T$ perfectly matches the eigen-decomposition format $S=Q\mathrm{\Lambda }{Q}^T$. This proves that the columns of $V$ are the orthonormal eigenvectors of $A^{T}A$. And just as before, because only the first $r$ singular values in $\mathrm{\Sigma }$ are nonzero, the first $r$ columns of $V$ correspond to the nonzero eigenvalues.

2. Question 2

Given a symmetric matrix $A\in {\mathbb{R}}^{3\times 3}$, suppose its eigen-decomposition can be written as

What is the singular value decomposition of this matrix?

This is a fantastic question that perfectly bridges the gap between Eigen-decomposition and Singular Value Decomposition (SVD). It tests a very specific and fundamental rule of linear algebra.

1. Identifying What We Have

First, let's look at the equation you were given. Because the matrix $A$ is a symmetric matrix, the equation is currently written in the exact format of an Eigen-Decomposition: $$A = Q \Lambda Q^T$$ Where:

• $Q$ is the orthogonal matrix on the left containing the eigenvectors.
• $\mathrm{\Lambda }$ (Lambda) is the diagonal matrix in the middle containing the eigenvalues: 3, -2, and 1.
• $Q^T$ is the transpose of the first matrix on the right.

2. The Core Problem: Eigen-Decomposition vs. SVD

The question asks you to convert this into a Singular Value Decomposition (SVD).

Recall that the formula for SVD is: $$A = U \Sigma V^T$$Where $U$ and $V$ are orthogonal matrices, and $\mathrm{\Sigma }$ (Sigma) is a diagonal matrix containing the singular values.

At first glance, the given eigen-decomposition $Q\mathrm{\Lambda }{Q}^T$ looks almost exactly identical to the SVD format $U\mathrm{\Sigma }{V}^T$. However, there are two strict mathematical rules for the $\mathrm{\Sigma }$ matrix in SVD that the current $\mathrm{\Lambda }$ matrix violates:

1. Non-Negative Rule: All singular values on the diagonal of $\mathrm{\Sigma }$ must be non-negative real numbers (greater than or equal to zero). Our current $\mathrm{\Lambda }$ matrix has a -2, which is illegal in SVD.
2. Sorted Rule: The singular values must be sorted in descending order from largest to smallest (${\sigma }_1\ge {\sigma }_2\ge {\sigma }_3$).

If a symmetric matrix is "positive semi-definite" (meaning all its eigenvalues are naturally positive), its eigen-decomposition and its SVD are perfectly identical. But because of that -2, we have to do a little bit of math to fix it.

3. The Mathematical Trick (Absorbing the Negative)

We need to turn that -2 into a +2 to satisfy the SVD rules. However, we cannot just erase a negative sign, because that would completely change the underlying matrix $A$!

Instead, we use a simple linear algebra trick: we factor the negative sign out of the diagonal matrix and absorb it into one of the other matrices.

In an eigen-decomposition ($Q\mathrm{\Lambda }{Q}^T$), the matrices on the left and right are exactly the same ($Q$). If we change the sign inside the right matrix, it is no longer exactly $Q^T$. This is perfectly fine for SVD, because SVD explicitly allows the left matrix ($U$) and the right matrix ($V^T$) to be two completely different matrices!

Here is the step-by-step math: We can split our diagonal matrix $\mathrm{\Lambda }$ into two separate diagonal matrices multiplied together: one containing the absolute values (which will become $\mathrm{\Sigma }$), and one containing the signs (let's call it $D$).

Now, let's substitute this back into the original equation: $$A = Q (\Sigma D) Q^T$$Because matrix multiplication is associative, we can group the sign matrix $D$ with the right-most matrix: $$A = Q \Sigma (D Q^T)$$

4. Constructing the Final SVD Answer

Now we simply map our new pieces to the SVD format $A=U\mathrm{\Sigma }{V}^T$:

• The $U$ Matrix: This remains exactly the same as the original $Q$ matrix.
• The $\mathrm{\Sigma }$ Matrix: This is our new diagonal matrix with strictly positive values. Conveniently, the values 3, 2, and 1 are already in the correct sorted descending order ($3\ge 2\ge 1\ge 0$)!
• The $V^T$ Matrix: This is the product of $(D{Q}^T)$. Multiplying a matrix by a diagonal matrix of signs simply flips the signs of the corresponding rows. Because the -1 is in the second row of $D$, it will multiply every element in the second row of $Q^T$ by -1. (Note: Changing the sign of an entire row of an orthogonal matrix keeps it completely orthogonal, so $V^T$ is still a valid orthogonal matrix!).

Answer

You can write out your final answer clearly like this:

"To find the Singular Value Decomposition $A=U\mathrm{\Sigma }{V}^T$, we must ensure that all singular values in the diagonal matrix $\mathrm{\Sigma }$ are non-negative and in descending order. The given eigen-decomposition contains a negative eigenvalue (-2). We can resolve this by factoring a -1 out of the diagonal matrix and absorbing it into the corresponding row of the right-hand orthogonal matrix. Because SVD allows $U$ and $V^T$ to be different matrices, this is a valid operation.

$A=\left(\begin{array}{lll}{u}_{11}& u_{12}& u_{13}\\ u_{21}& u_{22}& u_{23}\\ u_{31}& u_{32}& u_{33}\end{array}\right)\left(\begin{array}{lll}3& 0& 0\\ 0& 2& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{ccc}{u}_{11}& u_{21}& u_{31}\\ -u_{12}& -u_{22}& -u_{32}\\ u_{13}& u_{23}& u_{33}\end{array}\right)$

Where $\mathrm{\Sigma }$ now correctly contains only non-negative singular values sorted in descending order (${\sigma }_1\ge {\sigma }_2\ge {\sigma }_3$), and $V^T$ remains an orthogonal matrix."

3. Question 3

In this homework, you will apply the principal component analysis to a collection of handwritten digit images from the USPS dataset. The USPS dataset is in the "data" folder: USPS.mat. The starting code is in the "code" folder. The whole data has already been loaded into the matrix A. The matrix A has shape $3000\times 256$ and contains all the images. Each row in A corresponds to a handwritten digit image (between 0 and 9) with size $16\times 16$. You are expected to implement your solution based on the given codes. The only file you need to modify is the "solution.py" file. You can test your solution by running the "main.py" file.

3.1 do_pca()

1. (In Code) Complete the _do_pca() method. Your code will be tested on $\mathrm{p}=10,50,100,200$, a total of four different numbers of the principal components.

I would be happy to help you implement your PCA (Principal Component Analysis) code and explain how it works from scratch!

PCA can seem intimidating, but at its core, it is just a method for dimensionality reduction. In your dataset, each handwritten digit image is $16\times 16$ pixels, which flattens out into an array of $256$ numbers (features). Processing 256 features for every single image can be computationally heavy. PCA helps us compress these 256 features down to a much smaller number (like 10 or 50) while retaining the most important structural information (the "principal components") to successfully reconstruct or classify the image.

Here is the step-by-step breakdown of how PCA works from scratch, how the math aligns with your class notes, and exactly how to implement it in your code.

1. How PCA Works (The Theory)

Step A: Center the Data If we want to find the directions where the data varies the most, we first need to center the data around the origin $(0,0)$. We do this by calculating the mean (average) of each feature across all images, and subtracting that mean from every single image. As your notes point out, you must always center the data first if you want to achieve the minimal reconstruction error.

Step B: Find the Covariance and Eigenvectors Next, we need to find the "Principal Components." These are mathematically defined as the eigenvectors of the covariance matrix of your data.

• A covariance matrix tells us how much the pixels vary with one another.
• The eigenvectors of this matrix give us the best directional axes to project our data onto to maximize variance.

Step C: The SVD Shortcut Computing a huge covariance matrix and finding its eigenvectors directly is computationally very expensive. However, linear algebra gives us a magical shortcut: Singular Value Decomposition (SVD). If we take the SVD of our centered data matrix ($X=U\mathrm{\Sigma }{V}^T$), the columns of the $V$ matrix perfectly correspond to the orthonormal eigenvectors of the covariance matrix. Therefore, we completely skip calculating the covariance matrix and just use SVD directly on the centered data!

Step D: Truncate and Project SVD will give us 256 components, sorted from most important to least important. To reduce the dimensionality, we simply throw away the less important ones and keep the first $k$ components. Finally, we multiply our original centered data by these top components to project our data into the new, smaller dimension.

2. Implementing _do_pca(self)

In your class notes, the data matrix $X$ is written as having shape [features, samples]. However, in Python/NumPy and your assignment instructions, self.X is transposed to have shape [n_samples, n_features] (3000 x 256). We just have to adapt our NumPy code to match this shape!

Here is the code to place inside your _do_pca method:

    def _do_pca(self):
        '''
        To do PCA decomposition.
        Returns:
            Up: Principal components (transform matrix) of shape [n_features, n_components].
            Xp: The reduced data matrix after PCA of shape [n_samples, n_components].
        '''
        ### YOUR CODE HERE

        # 1. Center the Data
        # Calculate the mean of each feature (column) across all samples.
        mean_vec = np.mean(self.X, axis=0)

        # Subtract the mean from the data to center it.
        # (Numpy automatically broadcasts this subtraction across all rows)
        X_centered = self.X - mean_vec

        # 2. Compute the Singular Value Decomposition (SVD)
        # SVD decomposes X_centered into U, S, and V^T
        # full_matrices=False is standard practice to save memory
        U, S, Vt = np.linalg.svd(X_centered, full_matrices=False)

        # 3. Extract the Principal Components (Up)
        # np.linalg.svd returns V-transpose (Vt). We need to transpose it back to V.
        V = np.transpose(Vt)

        # The principal components are the first 'n_components' columns of V
        Up = V[:, :self.n_components]

        # 4. Project the Data (Xp)
        # Multiply the centered data by the principal components to reduce dimensionality
        Xp = np.dot(X_centered, Up)

        return Up, Xp

        ### END YOUR CODE

3. Detailed Code Explanation

• np.mean(self.X, axis=0): This collapses the 3,000 samples and returns an array of 256 values, representing the "average" digit image.
• self.X - mean_vec: This subtracts the average image from every single image in the dataset, effectively centering our data at zero. (Note: The README hints at using np.ones(), which is a more manual, mathematical way to broadcast the subtraction, e.g., self.X - np.dot(np.ones((n_samples, 1)), [mean_vec]). However, standard NumPy broadcasting self.X - mean_vec achieves the exact same thing much more efficiently).
• np.linalg.svd(...): This applies the SVD shortcut discussed in your notes. It returns U, S (the singular values $\mathrm{\Sigma }$), and Vt ($V^T$).
• V = np.transpose(Vt): Because SVD natively outputs the transpose of V ($V^T$), we must transpose it back to $V$ to access the eigenvectors as columns.
• V[:, :self.n_components]: This is where the dimensionality reduction actually happens. The notation : means "take all 256 feature rows", and :self.n_components means "only take the first p columns".
• np.dot(X_centered, Up): Finally, we take our 3000 x 256 dataset and matrix-multiply it by our 256 x p component matrix. The inner dimensions cancel out, leaving us with a reduced dataset Xp of shape 3000 x p!

Once you implement this and run python main.py, the test_pca() function will loop through p = 10, 50, 100, 200 to automatically test your implementation and save reconstructed images of the handwritten digits so you can visualize how much data is preserved at different compression sizes! Let me know if you need help with the reconstruction error function next!

3.2 reconstruction()

Complete the reconstruction() method to reconstruct the reduced data.

1. What is the Reconstruction Process?

In the _do_pca function, we took 256-dimensional images and compressed them down into a much smaller number of dimensions (like 10 or 50). This compressed data is your Xp matrix (which your class notes refer to as $Z$).

The reconstruction method is the "decompression" step. We want to take those 10 numbers and blow them back up into the original 256 pixels so we can actually look at the image and see how much quality we lost during the compression.

2. The Math Behind Reconstruction

According to your class notes, if your principal components matrix is $G$ and your reduced data is $Z$, the mathematical formula to reconstruct the data is $\check{X}=GZ$.

Because our data matrices in Python are transposed (the samples are rows, not columns), the order of our multiplication flips:

1. The Projection: We multiply our reduced data Xp by the transpose of our principal components self.Up. This projects our 10-dimensional data back into the 256-dimensional space.
2. Un-centering: There is one crucial final step. Remember that in _do_pca, we subtracted the "average image" (the mean) from every single image to center our data at the origin $(0,0)$. Because PCA only reconstructs this centered data, we must add that exact same mean back to shift the pixels back to their original brightness levels.

3. Implementing the Code

Here is the step-by-step code you can paste into your solution.py file:

	def reconstruction(self, Xp):
		'''
		To reconstruct reduced data given principal components Up.

		Args:
		Xp: The reduced data matrix after PCA of shape [n_samples, n_components].

		Return:
		X_re: The reconstructed matrix of shape [n_samples, n_features].
		'''
		### YOUR CODE HERE

		# 1. Transpose the principal components matrix
		# self.Up has shape [n_features, n_components].
		# We transpose it to [n_components, n_features] for multiplication.
		Up_T = np.transpose(self.Up)

		# 2. Decompress the data
		# We multiply the reduced data (Xp) by the transposed components (Up_T).
		# Shape math: [n_samples, n_components] @ [n_components, n_features] = [n_samples, n_features]
		X_re_centered = np.dot(Xp, Up_T)

		# 3. Calculate the mean of the original data
		# We need the original mean to un-center our reconstructed data
		mean_vec = np.mean(self.X, axis=0)

		# 4. Add the mean back
		# This restores the original pixel intensity values
		X_re = X_re_centered + mean_vec

		return X_re

		### END YOUR CODE

What is happening at each step?

• np.transpose(self.Up): We take our self.Up matrix (which holds our best directional axes) and transpose it so the matrix multiplication aligns properly.
• np.dot(Xp, Up_T): We matrix-multiply the compressed data by the axes. This takes the coordinates in the small space and maps them back to the coordinates in the large, 256-dimensional pixel space.
• np.mean(self.X, axis=0): We recalculate the average digit image across all 3,000 samples.
• X_re_centered + mean_vec: By adding the average image back, we reverse the centering step we performed during encoding. If we forgot this step, our reconstructed images would have negative pixel values and would look completely wrong when we tried to save them!

3.3 reconstruct_error()

Complete the reconstruct_error() function to measuring the reconstruction error.

1. What is the Reconstruction Error?

When we performed PCA in the _do_pca step, we threw away the less important principal components to compress our images from 256 pixels down to a smaller number like 10 or 50. Because we threw away information, the reconstructed images will not be 100% identical to the original images.

The reconstruction error is a mathematical way to measure exactly how much information we lost. We do this by comparing the original image matrix (A) directly against the decompressed/reconstructed image matrix (B). We want to know the total difference between them.

2. The Math: What is the "Square of the Frobenius Norm"?

The assignment asks you to return the square of the Frobenius norm of the matrix $A-B$.

While the name sounds intimidating, the Frobenius norm is simply the matrix equivalent of the standard Euclidean distance (the Pythagorean theorem) used for vectors. Here is how it breaks down in plain English:

1. $A-B$: You subtract the reconstructed matrix from the original matrix. This gives you a new "difference" matrix containing the exact pixel-by-pixel errors.
2. The Frobenius Norm ($|A-B{|}_F$): You square every single individual error value in that difference matrix, sum all of those squared values together, and then take the square root of that total sum.
3. The Square ($|A-B{|}_F^2$): The instructions explicitly ask for the square of the Frobenius norm. This simply means we cancel out the square root from the previous step. We are just left with the total sum of the squared pixel errors across the entire dataset.

3. Implementing the Code

In your assignment's Readme.txt file, your professor provided a hint to use the helper function np.linalg.norm(). By default, when you pass a 2D matrix into this NumPy function, it automatically calculates the Frobenius norm for you.

def reconstruct_error(A, B):
	'''
	To compute the reconstruction error.

	Args:
	A & B: Two matrices needed to be compared with. Should be of same shape.

	Return:
	error: the Frobenius norm's square of the matrix A-B. A scaler number.
	'''
	### YOUR CODE HERE

	# 1. Calculate the difference matrix between the original (A) and reconstructed (B) data
	difference_matrix = A - B

	# 2. Compute the Frobenius norm of this difference matrix
	# np.linalg.norm calculates the Frobenius norm by default for 2D arrays
	frobenius_norm = np.linalg.norm(difference_matrix)

	# 3. Square the Frobenius norm as requested by the instructions
	error = frobenius_norm ** 2

	return error

	### END YOUR CODE

What happens when you run this?

Now that all three functions (_do_pca, get_reduced, and reconstruct_error) are complete, you can run python main.py in your terminal!

The test_pca() function in main.py will automatically pass your original 3000 images (A) into your PCA class, compress them down to p=10, 50, 100, and 200 components, reconstruct them (A_re), and print the reconstruction error using this final function you just wrote.

As you watch the terminal output, you should see that as the number of principal components ($p$) increases, the reconstruction error will rapidly decrease. This is because you are keeping more of the singular values and principal components, allowing the SVD shortcut to build a much more accurate approximation of the original data. The generated images saved in your folder will also look much sharper and closer to the real handwritten digits as $p$ approaches 200!

3.4 Report

Run "main.py" to see the reconstruction results and summarize your observations from the results into a short report. When you run the "main.py" file, a subset (the first two) of the reconstructed images based on $\mathrm{p}=10,50,100,200$ principal components will be automatically saved on the "code" folder. Please attach these images to your report.

PCA Image Reconstruction Report

1. Quantitative Observations (Reconstruction Error) When applying Principal Component Analysis (PCA) to the USPS dataset (which consists of $16\times 16$ pixel images, or 256 total dimensions), the quantitative results clearly demonstrate that the reconstruction error strictly decreases as the number of principal components ($p$) increases.

Based on the algorithm's output, the squared Frobenius norm of the reconstruction errors are:

• $p=10$: 155,351.47
• $p=50$: 41,024.86
• $p=100$: 14,285.81
• $p=200$: 1,371.41

Mathematical Justification: This behavior perfectly aligns with the mathematical theory of PCA. The goal of PCA is to project high-dimensional features to a lower-dimensional space while simultaneously minimizing the reconstruction error. When we reconstruct the data using the first $p$ principal components, the minimum reconstruction error is mathematically equal to the sum of the squared singular values that we discarded ($\sum _{i=p+1}^n{\sigma }_i^2$). Therefore, as $p$ grows larger, we are throwing away fewer singular values (discarding less variance), which naturally forces the reconstruction error to shrink closer to zero.

2. Qualitative Observations (Visual Image Quality) (Note: Attach your generated images for $p=10,50,100,$ and $200$ here)

The mathematical drop in reconstruction error is directly reflected in the visual quality of the reconstructed digit images:

• At $p=10$: The images are highly compressed. Because we are only keeping 10 out of 256 dimensions, the reconstruction only captures the most dominant statistical features of the digits (like overall intensity and basic stroke location). The resulting images will look highly blurry, blocky, and lack fine details.
• At $p=50$ and $p=100$: As we retain more principal components, the digits become much more distinguishable and the edges become sharper. The structural shape of the specific handwritten numbers is mostly restored.
• At $p=200$: Because 200 components is very close to the original 256 dimensions of the data, nearly all of the original variance is retained. The reconstructed images at this stage will appear exceptionally sharp and nearly identical to the original uncompressed data, which corresponds to the very small reconstruction error of $1371.41$.

3. Conclusion The results illustrate the fundamental approximation-generalization tradeoff of PCA dimensionality reduction. By keeping only the first $p$ principal components, we can achieve massive feature reduction (data compression). However, smaller values of $p$ lead to a higher loss of fine detail (higher reconstruction error), whereas higher values of $p$ preserve the image perfectly but result in less compression.