Lecture 6: PCA and Dimensionality Reduction

Code

import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots

# Set random seed for reproducibility
np.random.seed(42)

# Generate 2D MVN data with correlation 0.7
n_points = 500
mean = [0, 0]
cov = [[1, 0.7], [0.7, 1]]
X_original = np.random.multivariate_normal(mean, cov, n_points)

# Define transformation functions
def get_transformation(name):
    if name == 'Identity':
        return np.eye(2), np.array([0, 0])
    elif name == 'Rotate 45°':
        c, s = np.cos(np.pi/4), np.sin(np.pi/4)
        return np.array([[c, -s], [s, c]]), np.array([0, 0])
    elif name == 'Rotate 90°':
        c, s = np.cos(np.pi/2), np.sin(np.pi/2)
        return np.array([[c, -s], [s, c]]), np.array([0, 0])
    elif name == 'Rotate 180°':
        c, s = np.cos(np.pi), np.sin(np.pi)
        return np.array([[c, -s], [s, c]]), np.array([0, 0])
    elif name == 'Rotate 270°':
        c, s = np.cos(3*np.pi/2), np.sin(3*np.pi/2)
        return np.array([[c, -s], [s, c]]), np.array([0, 0])
    elif name == 'Stretch X by 2':
        return np.array([[2, 0], [0, 1]]), np.array([0, 0])
    elif name == 'Rotate 45° + Stretch X':
        c, s = np.cos(np.pi/4), np.sin(np.pi/4)
        R = np.array([[c, -s], [s, c]])
        S = np.array([[2, 0], [0, 1]])
        return S @ R, np.array([0, 0])
    elif name == 'Shift by (1,1)':
        return np.eye(2), np.array([1, 1])
    elif name == 'Decorrelate (PCA)':
        eigenvalues, eigenvectors = np.linalg.eig(cov)
        return eigenvectors.T, np.array([0, 0])

# Create figure with subplots
fig = make_subplots(
    rows=1, cols=2,
    subplot_titles=('Original Data', 'Transformed Data'),
    horizontal_spacing=0.15
)

# Transformation options
transformations = ['Identity', 'Rotate 45°', 'Rotate 90°', 'Rotate 180°', 
                   'Rotate 270°', 'Stretch X by 2', 'Rotate 45° + Stretch X', 
                   'Shift by (1,1)', 'Decorrelate (PCA)']

# Add traces for each transformation (initially invisible)
for i, trans_name in enumerate(transformations):
    A, b = get_transformation(trans_name)
    X_transformed = (A @ X_original.T).T + b
    
    # Original data (always visible)
    fig.add_trace(
        go.Scatter(
            x=X_original[:, 0],
            y=X_original[:, 1],
            mode='markers',
            marker=dict(size=4, color='blue', opacity=0.6),
            name='Original',
            showlegend=False,
            visible=(i == 0)  # Only first one visible
        ),
        row=1, col=1
    )
    
    # Transformed data
    fig.add_trace(
        go.Scatter(
            x=X_transformed[:, 0],
            y=X_transformed[:, 1],
            mode='markers',
            marker=dict(size=4, color='red', opacity=0.6),
            name='Transformed',
            showlegend=False,
            visible=(i == 0)  # Only first one visible
        ),
        row=1, col=2
    )

# Get initial transformation for annotation
A_init, b_init = get_transformation('Identity')

# Create dropdown menu
buttons = []
for i, trans_name in enumerate(transformations):
    A, b = get_transformation(trans_name)
    
    # Format matrix A and vector b as text - horizontal layout
    matrix_text = (f"A = [{A[0,0]:6.3f}  {A[0,1]:6.3f}]     "
                   f"b = [{b[0]:6.3f}]<br>"
                   f"    [{A[1,0]:6.3f}  {A[1,1]:6.3f}]         "
                   f"    [{b[1]:6.3f}]")
    
    # Create visibility array: show traces for this transformation
    visible = [False] * (len(transformations) * 2)
    visible[i * 2] = True      # Original data for this transformation
    visible[i * 2 + 1] = True  # Transformed data for this transformation
    
    buttons.append(
        dict(
            label=trans_name,
            method='update',
            args=[
                {'visible': visible},
                {
                    'title.text': f'Affine Transformation: {trans_name}',
                    'annotations': [
                        dict(
                            text=matrix_text,
                            x=0.5,
                            y=1.05,  # Moved above the plot
                            xref='paper',
                            yref='paper',
                            showarrow=False,
                            font=dict(size=13, family='monospace'),
                            align='center',
                            xanchor='center',
                            yanchor='bottom'
                        )
                    ]
                }
            ]
        )
    )

# Initial annotation text - horizontal layout
matrix_text_init = (f"A = [{A_init[0,0]:6.3f}  {A_init[0,1]:6.3f}]     "
                    f"b = [{b_init[0]:6.3f}]<br>"
                    f"    [{A_init[1,0]:6.3f}  {A_init[1,1]:6.3f}]         "
                    f"    [{b_init[1]:6.3f}]")

# Add dropdown menu
fig.update_layout(
    updatemenus=[
        dict(
            buttons=buttons,
            direction='down',
            showactive=True,
            x=0.5,
            xanchor='center',
            y=0.98,  # Moved down slightly to make room for annotation
            yanchor='top'
        )
    ],
    title_text='Affine Transformation: Identity',
    height=550,
    showlegend=False,
    margin=dict(t=120),  # Increased top margin for annotation
    annotations=[
        dict(
            text=matrix_text_init,
            x=0.5,
            y=1.05,  # Above the plot
            xref='paper',
            yref='paper',
            showarrow=False,
            font=dict(size=13, family='monospace'),
            align='center',
            xanchor='center',
            yanchor='bottom'
        )
    ]
)

# Set axes to be equal and centered
axis_range = [-4, 4]
fig.update_xaxes(range=axis_range, constrain='domain', row=1, col=1)
fig.update_yaxes(range=axis_range, scaleanchor="x", scaleratio=1, row=1, col=1)
fig.update_xaxes(range=axis_range, constrain='domain', row=1, col=2)
fig.update_yaxes(range=axis_range, scaleanchor="x2", scaleratio=1, row=1, col=2)

# Add grid
fig.update_xaxes(showgrid=True, gridwidth=1, gridcolor='lightgray')
fig.update_yaxes(showgrid=True, gridwidth=1, gridcolor='lightgray')

fig.show()

MNIST Digits Dataset

Let’s start with the MNIST handwritted digits. Out goal is to explore how PCA can be used to reduce the dimension of image data, and how this affects classification accuracy.

Let’s start with two baseline models: multinomial logistic regression and XGBoost on the full 784-dimensional data.

Code

import torch
import torchvision
import numpy as np
import xgboost as xgb
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score

# 1. Load MNIST dataset
print("Loading MNIST dataset...")
transform = torchvision.transforms.ToTensor()
mnist_train = torchvision.datasets.MNIST(root='./data', train=True, download=True, transform=transform)
mnist_test = torchvision.datasets.MNIST(root='./data', train=False, download=True, transform=transform)

# Convert to numpy and flatten (28x28 -> 784)
X_train_mnist = mnist_train.data.numpy().reshape(-1, 784).astype(np.float32) / 255.0
y_train_mnist = mnist_train.targets.numpy()
X_test_mnist = mnist_test.data.numpy().reshape(-1, 784).astype(np.float32) / 255.0
y_test_mnist = mnist_test.targets.numpy()

# 2. Center the data
print("Centering MNIST data...")
mean_mnist = np.mean(X_train_mnist, axis=0)
X_train_mnist = X_train_mnist - mean_mnist
X_test_mnist = X_test_mnist - mean_mnist

print(f"Training set shape: {X_train_mnist.shape}")
print(f"Test set shape: {X_test_mnist.shape}\n")

# 3. Multinomial Logistic Regression
print("Training Multinomial Logistic Regression on MNIST...")
log_reg_mnist = LogisticRegression(
    multi_class='multinomial',
    solver='lbfgs',
    max_iter=1000,
    random_state=42
)
log_reg_mnist.fit(X_train_mnist, y_train_mnist)

train_acc_lr = log_reg_mnist.score(X_train_mnist, y_train_mnist)
test_acc_lr = log_reg_mnist.score(X_test_mnist, y_test_mnist)

print(f"Logistic Regression Training Accuracy: {train_acc_lr:.4f}")
print(f"Logistic Regression Test Accuracy:     {test_acc_lr:.4f}\n")

# 4. XGBoost
print("Training XGBoost on MNIST (this may take a minute)...")
xgb_clf_mnist = xgb.XGBClassifier(
    n_estimators=100,
    max_depth=6,
    learning_rate=0.1,
    random_state=42,
    eval_metric='mlogloss',
    tree_method='hist'  # Faster traiing
)
xgb_clf_mnist.fit(X_train_mnist, y_train_mnist)

train_acc_xgb = xgb_clf_mnist.score(X_train_mnist, y_train_mnist)
test_acc_xgb = xgb_clf_mnist.score(X_test_mnist, y_test_mnist)

print(f"XGBoost Training Accuracy: {train_acc_xgb:.4f}")
print(f"XGBoost Test Accuracy:     {test_acc_xgb:.4f}")

Loading MNIST dataset...
Centering MNIST data...
Training set shape: (60000, 784)
Test set shape: (10000, 784)

Training Multinomial Logistic Regression on MNIST...
Logistic Regression Training Accuracy: 0.9393
Logistic Regression Test Accuracy:     0.9256

Training XGBoost on MNIST (this may take a minute)...
XGBoost Training Accuracy: 0.9960
XGBoost Test Accuracy:     0.9702

Now, figure out the reconstruction error with PCA:

Code

from sklearn.decomposition import PCA

import matplotlib.pyplot as plt

# Fit PCA on the full MNIST training data to get all components
print("Fitting PCA with all components on MNIST to compute reconstruction error...")
pca_mnist_full = PCA(random_state=42)
pca_mnist_full.fit(X_train_mnist)

# Get all eigenvalues (explained variance)
eigenvalues_mnist = pca_mnist_full.explained_variance_

# Define range of k values to test
k_values_mnist = list(range(1, 785))

# Compute reconstruction error for each k
reconstruction_errors_mnist = []

for k in k_values_mnist:
    # Average of eigenvalues from k onwards (the ones we didn't keep)
    if k < len(eigenvalues_mnist):
        leftover_eigenvalues_mnist = eigenvalues_mnist[k:]
        reconstruction_error_mnist = np.mean(leftover_eigenvalues_mnist)
    else:
        reconstruction_error_mnist = 0.0
    
    reconstruction_errors_mnist.append(reconstruction_error_mnist)

print(f"Computed reconstruction errors for k = 1 to {len(k_values_mnist)}")

# Plot reconstruction error vs k
plt.figure(figsize=(12, 6))
plt.plot(k_values_mnist, reconstruction_errors_mnist, 'b-', linewidth=2)
plt.xlabel('Number of PCA Components (k)', fontsize=12)
plt.ylabel('Reconstruction Error (Avg. Leftover Eigenvalues)', fontsize=12)
plt.title('MNIST: PCA Reconstruction Error vs Number of Components', fontsize=14, fontweight='bold')
plt.grid(True, alpha=0.3)
plt.xlim([1, 784])
plt.tight_layout()
plt.show()

Fitting PCA with all components on MNIST to compute reconstruction error...
Computed reconstruction errors for k = 1 to 784

Let’s see if we can improve our accuracy for MNL using PCA to reduce dimension. We’ll try dimensions from 20 to 500 in steps of 20, plus the full 784 dimensions. We’ll then assess the accuracy on the test set.

Code

print("Running PCA + Multinomial Logistic Regression for different dimensions on MNIST...")
print("=" * 80)

# Define the range of PCA components: 20 to 500 in steps of 20
k_values_acc = list(range(20, 501, 20))
# Add full dimension (784)
k_values_acc.append(784)

# Store accuracies
train_accuracies_mnist_pca = []
test_accuracies_mnist_pca = []

for k in k_values_acc:
    if k == 784:
        # No PCA needed for full dimension, use original data
        X_train_reduced = X_train_mnist
        X_test_reduced = X_test_mnist
    else:
        # Fit PCA
        pca = PCA(n_components=k, random_state=42)
        X_train_reduced = pca.fit_transform(X_train_mnist)
        X_test_reduced = pca.transform(X_test_mnist)
    
    # Train Logistic Regression
    # Increasing max_iter slightly to ensure convergence for all k
    lr = LogisticRegression(
        multi_class='multinomial',
        solver='lbfgs',
        max_iter=1000,
        random_state=42
    )
    lr.fit(X_train_reduced, y_train_mnist)
    
    # Calculate accuracy
    train_acc = lr.score(X_train_reduced, y_train_mnist)
    test_acc = lr.score(X_test_reduced, y_test_mnist)
    
    train_accuracies_mnist_pca.append(train_acc)
    test_accuracies_mnist_pca.append(test_acc)
    
    print(f"k = {k:3d} | Train Acc: {train_acc:.4f} | Test Acc: {test_acc:.4f}")

print("=" * 80)

# Plot accuracy vs k
plt.figure(figsize=(12, 6))
plt.plot(k_values_acc, train_accuracies_mnist_pca, 'b-o', label='Training Accuracy', linewidth=2, markersize=5)
plt.plot(k_values_acc, test_accuracies_mnist_pca, 'r-o', label='Test Accuracy', linewidth=2, markersize=5)
plt.xlabel('Number of PCA Components (k)', fontsize=12)
plt.ylabel('Accuracy', fontsize=12)
plt.title('MNIST: Logistic Regression Accuracy vs PCA Dimensions', fontsize=14, fontweight='bold')
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Find best test accuracy
best_idx_mnist = np.argmax(test_accuracies_mnist_pca)
best_k_mnist = k_values_acc[best_idx_mnist]
best_test_acc_mnist = test_accuracies_mnist_pca[best_idx_mnist]

print(f"Best test accuracy: {best_test_acc_mnist:.4f} at k = {best_k_mnist}")

Running PCA + Multinomial Logistic Regression for different dimensions on MNIST...
================================================================================
k =  20 | Train Acc: 0.8753 | Test Acc: 0.8810
k =  40 | Train Acc: 0.9030 | Test Acc: 0.9071
k =  60 | Train Acc: 0.9138 | Test Acc: 0.9165
k =  80 | Train Acc: 0.9187 | Test Acc: 0.9187
k = 100 | Train Acc: 0.9227 | Test Acc: 0.9213
k = 120 | Train Acc: 0.9239 | Test Acc: 0.9224
k = 140 | Train Acc: 0.9261 | Test Acc: 0.9242
k = 160 | Train Acc: 0.9270 | Test Acc: 0.9248
k = 180 | Train Acc: 0.9297 | Test Acc: 0.9254
k = 200 | Train Acc: 0.9304 | Test Acc: 0.9266
k = 220 | Train Acc: 0.9311 | Test Acc: 0.9251
k = 240 | Train Acc: 0.9320 | Test Acc: 0.9258
k = 260 | Train Acc: 0.9321 | Test Acc: 0.9257
k = 280 | Train Acc: 0.9330 | Test Acc: 0.9256
k = 300 | Train Acc: 0.9334 | Test Acc: 0.9252
k = 320 | Train Acc: 0.9339 | Test Acc: 0.9249
k = 340 | Train Acc: 0.9350 | Test Acc: 0.9241
k = 360 | Train Acc: 0.9355 | Test Acc: 0.9247
k = 380 | Train Acc: 0.9359 | Test Acc: 0.9249
k = 400 | Train Acc: 0.9363 | Test Acc: 0.9256
k = 420 | Train Acc: 0.9371 | Test Acc: 0.9255
k = 440 | Train Acc: 0.9378 | Test Acc: 0.9256
k = 460 | Train Acc: 0.9381 | Test Acc: 0.9260
k = 480 | Train Acc: 0.9383 | Test Acc: 0.9265
k = 500 | Train Acc: 0.9385 | Test Acc: 0.9255
k = 784 | Train Acc: 0.9393 | Test Acc: 0.9256
================================================================================
Best test accuracy: 0.9266 at k = 200

Back to CIFAR-10

Let’s start by loading in the CIFAR-10 data again and using XGBoost for the ten class classification problem. Judge the accuracy on a held out data set of size 10k.

Code

import numpy as np
import matplotlib.pyplot as plt
import torch
import torchvision
import torchvision.transforms as transforms

# Load CIFAR-10 train and test sets
print("Loading CIFAR-10 dataset...")
transform = transforms.ToTensor()
trainset = torchvision.datasets.CIFAR10(root='./data', train=True,
                                        download=True, transform=transform)
testset = torchvision.datasets.CIFAR10(root='./data', train=False,
                                       download=True, transform=transform)

# CIFAR-10 class names
class_names = ['airplane', 'automobile', 'bird', 'cat', 'deer', 
               'dog', 'frog', 'horse', 'ship', 'truck']

print(f"Training set size: {len(trainset)} images")
print(f"Test set size: {len(testset)} images")
print(f"Image shape: 32×32×3 = 3,072 dimensions")
print(f"Number of classes: {len(class_names)}")

# Find one example from each class
fig, axes = plt.subplots(2, 5, figsize=(6, 2.5))  # Much smaller figure
axes = axes.flatten()

# Get one example per class
examples_found = {i: None for i in range(10)}
for idx in range(len(trainset)):
    img, label = trainset[idx]
    if examples_found[label] is None:
        examples_found[label] = img
    if all(v is not None for v in examples_found.values()):
        break

# Display
for class_idx, (ax, class_name) in enumerate(zip(axes, class_names)):
    img = examples_found[class_idx]
    # Convert from CxHxW to HxWxC for matplotlib
    img_np = img.permute(1, 2, 0).numpy()
    
    ax.imshow(img_np)
    ax.set_title(class_name, fontsize=8, fontweight='bold')  # Smaller font
    ax.axis('off')

plt.suptitle('CIFAR-10: One Example from Each Class (p = 3,072 dimensions)', 
             fontsize=10, fontweight='bold', y=1.00)  # Smaller title
plt.tight_layout()
plt.show()

print(f"\nDimensionality: 3,072 features (32×32×3 RGB pixels)")
print(f"Training samples: {len(trainset)} images")
print(f"Test samples: {len(testset)} images")

Loading CIFAR-10 dataset...
Training set size: 50000 images
Test set size: 10000 images
Image shape: 32×32×3 = 3,072 dimensions
Number of classes: 10

Dimensionality: 3,072 features (32×32×3 RGB pixels)
Training samples: 50000 images
Test samples: 10000 images

Fit the XGBoost model with reasonable hyperparameter values.

Code

import xgboost as xgb
from sklearn.metrics import confusion_matrix
import seaborn as sns

# Flatten CIFAR-10 images for logistic regression
print("Preparing CIFAR-10 data for Logistic Regression...")

# Convert datasets to numpy arrays and flatten
X_train_cifar = []
y_train_cifar = []
for img, label in trainset:
    X_train_cifar.append(img.numpy().flatten())
    y_train_cifar.append(label)

X_test_cifar = []
y_test_cifar = []
for img, label in testset:
    X_test_cifar.append(img.numpy().flatten())
    y_test_cifar.append(label)

X_train_cifar = np.array(X_train_cifar)
y_train_cifar = np.array(y_train_cifar)
X_test_cifar = np.array(X_test_cifar)
y_test_cifar = np.array(y_test_cifar)

# Center the data
print("Centering the data...")
mean_cifar = np.mean(X_train_cifar, axis=0)
X_train_cifar = X_train_cifar - mean_cifar
X_test_cifar = X_test_cifar - mean_cifar

Preparing CIFAR-10 data for Logistic Regression...
Centering the data...

Code

# Fit XGBoost Classifier on CIFAR-10
print("Fitting XGBoost Classifier on CIFAR-10...")
print("This may take a few minutes...\n")

xgb_clf_cifar = xgb.XGBClassifier(
    n_estimators=100,
    max_depth=6,
    learning_rate=0.1,
    random_state=42,
    eval_metric='mlogloss',
    tree_method='hist'  # Faster for large datasets
)

xgb_clf_cifar.fit(X_train_cifar, y_train_cifar)

# Compute accuracies
train_acc_cifar_xgb = xgb_clf_cifar.score(X_train_cifar, y_train_cifar)
test_acc_cifar_xgb = xgb_clf_cifar.score(X_test_cifar, y_test_cifar)

print(f"XGBoost on CIFAR-10")
print(f"Training Accuracy: {train_acc_cifar_xgb:.4f}")
print(f"Test Accuracy: {test_acc_cifar_xgb:.4f}\n")

# Generate predictions for confusion matrix
y_pred_cifar_xgb = xgb_clf_cifar.predict(X_test_cifar)

# Compute confusion matrix
cm_xgb = confusion_matrix(y_test_cifar, y_pred_cifar_xgb)

# Plot confusion matrix
plt.figure(figsize=(10, 8))
sns.heatmap(cm_xgb, annot=True, fmt='d', cmap='Blues', 
            xticklabels=class_names, yticklabels=class_names)
plt.title(f'CIFAR-10 XGBoost Confusion Matrix\nTest Accuracy: {test_acc_cifar_xgb:.4f}')
plt.ylabel('True Label')
plt.xlabel('Predicted Label')
plt.tight_layout()
plt.show()

# Print per-class accuracy
print("\nPer-class accuracy:")
for i, class_name in enumerate(class_names):
    class_accuracy = cm_xgb[i, i] / cm_xgb[i, :].sum()
    print(f"{class_name:12s}: {class_accuracy:.4f}")

Fitting XGBoost Classifier on CIFAR-10...
This may take a few minutes...

XGBoost on CIFAR-10
Training Accuracy: 0.8401
Test Accuracy: 0.5252


Per-class accuracy:
airplane    : 0.5880
automobile  : 0.6040
bird        : 0.3860
cat         : 0.3490
deer        : 0.4260
dog         : 0.4430
frog        : 0.6370
horse       : 0.5350
ship        : 0.6690
truck       : 0.6150

To start, let’s check the reconstruction error of PCA for CIFAR-10 as we vary the number of components used.

Code

from sklearn.decomposition import PCA

# Fit PCA on the full training data to get all components
print("Fitting PCA with all components to compute reconstruction error...")
pca_full = PCA(random_state=42)
pca_full.fit(X_train_cifar)

# Get all eigenvalues (explained variance)
eigenvalues = pca_full.explained_variance_

# Define range of k values to test
k_values = list(range(1, 3073))

# Compute reconstruction error for each k
# Reconstruction error = average of leftover eigenvalues (those not used)
reconstruction_errors = []

for k in k_values:
    # Average of eigenvalues from k onwards (the ones we didn't keep)
    if k < len(eigenvalues):
        leftover_eigenvalues = eigenvalues[k:]
        reconstruction_error = np.mean(leftover_eigenvalues)
    else:
        reconstruction_error = 0.0
    
    reconstruction_errors.append(reconstruction_error)

print(f"Computed reconstruction errors for k = 1 to {len(k_values)}")

# Plot reconstruction error vs k
plt.figure(figsize=(12, 6))
plt.plot(k_values, reconstruction_errors, 'b-', linewidth=2)
plt.xlabel('Number of PCA Components (k)', fontsize=12)
plt.ylabel('Reconstruction Error (Avg. Leftover Eigenvalues)', fontsize=12)
plt.title('CIFAR-10: PCA Reconstruction Error vs Number of Components', fontsize=14, fontweight='bold')
plt.grid(True, alpha=0.3)
plt.xlim([1, 3072])
plt.tight_layout()
plt.show()

Fitting PCA with all components to compute reconstruction error...
Computed reconstruction errors for k = 1 to 3072

Looking at this plot, we can see that the sweet spot seems to be somewhere between 300 and 600.

Now, let’s do PCA on the data to reduce the dimension from 3072 to something much smaller. We’re going to try (300,320,…,600). For each of these values, we’ll project the data down to that dimension using PCA, fit a multinomial logistic regression on the training data, and then judge the accuracy on the same validation set.

Code

from sklearn.decomposition import PCA
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
# Fit multinomial logistic regression on original features (full dimensionality)
print("Fitting Multinomial Logistic Regression on full features (3072 dimensions)...")

log_reg_full = LogisticRegression(
    multi_class='multinomial',
    solver='lbfgs',
    max_iter=1000,
    random_state=42
)

# Fit on the centered original data
log_reg_full.fit(X_train_cifar, y_train_cifar)

# Compute accuracies
train_acc_full = log_reg_full.score(X_train_cifar, y_train_cifar)
test_acc_full = log_reg_full.score(X_test_cifar, y_test_cifar)

print(f"\nLogistic Regression on Full CIFAR-10 (3072 features)")
print(f"Training Accuracy: {train_acc_full:.4f}")
print(f"Test Accuracy:     {test_acc_full:.4f}")

Fitting Multinomial Logistic Regression on full features (3072 dimensions)...

Logistic Regression on Full CIFAR-10 (3072 features)
Training Accuracy: 0.4991
Test Accuracy:     0.3856

/home/kmcalist/.local/lib/python3.10/site-packages/sklearn/linear_model/_logistic.py:460: ConvergenceWarning:

lbfgs failed to converge (status=1):
STOP: TOTAL NO. of ITERATIONS REACHED LIMIT.

Increase the number of iterations (max_iter) or scale the data as shown in:
    https://scikit-learn.org/stable/modules/preprocessing.html
Please also refer to the documentation for alternative solver options:
    https://scikit-learn.org/stable/modules/linear_model.html#logistic-regression

Code

# Use the full training set (already centered) and test set for PCA
# No need to split into train/validation since we're using the existing test set

# Define the range of PCA components to try
n_components_list = list(range(300, 601, 20))

# Store accuracies
train_accuracies = []
test_accuracies = []

print("Running PCA + Multinomial Logistic Regression for different dimensions...")
print("=" * 80)

for n_comp in n_components_list:
    # Apply PCA
    pca = PCA(n_components=n_comp, random_state=42)
    X_train_pca = pca.fit_transform(X_train_cifar)
    X_test_pca = pca.transform(X_test_cifar)
    
    # Fit multinomial logistic regression
    log_reg = LogisticRegression(
        multi_class='multinomial',
        solver='lbfgs',
        max_iter=1000,
        random_state=42
    )
    log_reg.fit(X_train_pca, y_train_cifar)
    
    # Compute accuracies
    train_acc = log_reg.score(X_train_pca, y_train_cifar)
    test_acc = log_reg.score(X_test_pca, y_test_cifar)
    
    train_accuracies.append(train_acc)
    test_accuracies.append(test_acc)
    
    if n_comp % 100 == 0 or n_comp == 20:
        print(f"n_components = {n_comp:4d} | Train: {train_acc:.4f} | Test: {test_acc:.4f}")

print("=" * 80)
print("\nCompleted!\n")

Running PCA + Multinomial Logistic Regression for different dimensions...
================================================================================
n_components =  300 | Train: 0.4343 | Test: 0.4071
n_components =  400 | Train: 0.4417 | Test: 0.4091
n_components =  500 | Train: 0.4459 | Test: 0.4044
n_components =  600 | Train: 0.4504 | Test: 0.4028
================================================================================

Completed!

Code

# Plot the results
plt.figure(figsize=(12, 6))

plt.plot(n_components_list, train_accuracies, 'b-', label='Training Accuracy', linewidth=2)
plt.plot(n_components_list, test_accuracies, 'r-', label='Test Accuracy', linewidth=2)

plt.xlabel('Number of PCA Components', fontsize=12)
plt.ylabel('Accuracy', fontsize=12)
plt.title('CIFAR-10: Multinomial Logistic Regression Accuracy vs PCA Dimensions', fontsize=14, fontweight='bold')
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)
plt.xlim([300, 600])
plt.tight_layout()
plt.show()

# Find best test accuracy
best_idx = np.argmax(test_accuracies)
best_n_comp = n_components_list[best_idx]
best_test_acc = test_accuracies[best_idx]

print(f"Best test accuracy: {best_test_acc:.4f} at {best_n_comp} components")
print(f"Corresponding test accuracy: {best_test_acc:.4f}")

Best test accuracy: 0.4121 at 340 components
Corresponding test accuracy: 0.4121

Code

print("Running PCA + XGBoost for different dimensions...")
print("=" * 80)

# Define the range of PCA components to try for XGBoost (50, 100, ..., 500)
n_components_xgb_list = list(range(50, 501, 50))

# Store accuracies
train_accuracies_xgb_pca = []
test_accuracies_xgb_pca = []

for n_comp in n_components_xgb_list:
    # Apply PCA
    pca_xgb = PCA(n_components=n_comp, random_state=42)
    X_train_pca = pca_xgb.fit_transform(X_train_cifar)
    X_test_pca = pca_xgb.transform(X_test_cifar)
    
    # Fit XGBoost Classifier
    # Using same parameters as before for consistency
    xgb_clf_pca = xgb.XGBClassifier(
        n_estimators=100,
        max_depth=6,
        learning_rate=0.1,
        random_state=42,
        eval_metric='mlogloss',
        tree_method='hist' 
    )
    
    xgb_clf_pca.fit(X_train_pca, y_train_cifar)
    
    # Compute accuracies
    train_acc = xgb_clf_pca.score(X_train_pca, y_train_cifar)
    test_acc = xgb_clf_pca.score(X_test_pca, y_test_cifar)
    
    train_accuracies_xgb_pca.append(train_acc)
    test_accuracies_xgb_pca.append(test_acc)
    
    print(f"n_components = {n_comp:4d} | Train: {train_acc:.4f} | Test: {test_acc:.4f}")

print("=" * 80)
print("\nCompleted!\n")

# Plot the results
plt.figure(figsize=(12, 6))

plt.plot(n_components_xgb_list, train_accuracies_xgb_pca, 'b-', label='Training Accuracy', linewidth=2)
plt.plot(n_components_xgb_list, test_accuracies_xgb_pca, 'r-', label='Test Accuracy', linewidth=2)

plt.xlabel('Number of PCA Components', fontsize=12)
plt.ylabel('Accuracy', fontsize=12)
plt.title('CIFAR-10: XGBoost Accuracy vs PCA Dimensions', fontsize=14, fontweight='bold')
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)
plt.xticks(n_components_xgb_list)
plt.tight_layout()
plt.show()

# Find best test accuracy
best_idx_xgb = np.argmax(test_accuracies_xgb_pca)
best_n_comp_xgb = n_components_xgb_list[best_idx_xgb]
best_test_acc_xgb = test_accuracies_xgb_pca[best_idx_xgb]

print(f"Best XGBoost test accuracy: {best_test_acc_xgb:.4f} at {best_n_comp_xgb} components")

Running PCA + XGBoost for different dimensions...
================================================================================
n_components =   50 | Train: 0.6903 | Test: 0.4694
n_components =  100 | Train: 0.7354 | Test: 0.4836
n_components =  150 | Train: 0.7528 | Test: 0.4857
n_components =  200 | Train: 0.7635 | Test: 0.4861
n_components =  250 | Train: 0.7709 | Test: 0.4866
n_components =  300 | Train: 0.7751 | Test: 0.4856
n_components =  350 | Train: 0.7822 | Test: 0.4870
n_components =  400 | Train: 0.7842 | Test: 0.4883
n_components =  450 | Train: 0.7876 | Test: 0.4883
n_components =  500 | Train: 0.7913 | Test: 0.4850
================================================================================

Completed!

Best XGBoost test accuracy: 0.4883 at 400 components

It’s a nonlinear manifold embedded in a higher dimensional space! PCA can only get us so far - dimensionality reduction can only get us so far if we aren’t capturing the true underlying structure of the data.

Instead, let’s do it end to end with a neural network.

Code

import torch
from torch.utils.data import DataLoader, TensorDataset
import pytorch_lightning as pl
from pytorch_lightning.callbacks import ModelCheckpoint

import torch.nn as nn
import torch.nn.functional as F

# Define the neural network model
class CIFAR10Classifier(pl.LightningModule):
    def __init__(self, input_dim=3072, hidden_dim=1024, num_classes=10, learning_rate=0.01, dropout=0.1):
        super().__init__()
        self.save_hyperparameters()
        
        # 4 hidden layers with batch normalization
        self.fc1 = nn.Linear(input_dim, 1024)
        self.bn1 = nn.BatchNorm1d(1024)
        
        self.fc2 = nn.Linear(1024, 512)
        self.bn2 = nn.BatchNorm1d(512)
        
        self.fc3 = nn.Linear(512, 256)
        self.bn3 = nn.BatchNorm1d(256)
        
        #self.fc4 = nn.Linear(hidden_dim, hidden_dim)
        #self.bn4 = nn.BatchNorm1d(hidden_dim)
        
        self.fc5 = nn.Linear(256, num_classes)
        
        self.dropout = nn.Dropout(dropout)
        self.learning_rate = learning_rate
        
    def forward(self, x):
        x = self.fc1(x)
        x = self.bn1(x)
        x = F.relu(x)
        x = self.dropout(x)
        
        x = self.fc2(x)
        x = self.bn2(x)
        x = F.relu(x)
        x = self.dropout(x)
        
        x = self.fc3(x)
        x = self.bn3(x)
        x = F.relu(x)
        x = self.dropout(x)
        
        #x = self.fc4(x)
        #x = self.bn4(x)
        #x = F.relu(x)
        x = self.dropout(x)
        
        x = self.fc5(x)  # Logits for cross-entropy
        return x
    
    def training_step(self, batch, batch_idx):
        x, y = batch
        logits = self(x)
        loss = F.cross_entropy(logits, y)
        acc = (logits.argmax(dim=1) == y).float().mean()
        
        self.log('train_loss', loss, prog_bar=True)
        self.log('train_acc', acc, prog_bar=True)
        return loss
    
    def validation_step(self, batch, batch_idx):
        x, y = batch
        logits = self(x)
        loss = F.cross_entropy(logits, y)
        acc = (logits.argmax(dim=1) == y).float().mean()
        
        self.log('val_loss', loss, prog_bar=True)
        self.log('val_acc', acc, prog_bar=True)
        return loss
    
    def test_step(self, batch, batch_idx):
        x, y = batch
        logits = self(x)
        loss = F.cross_entropy(logits, y)
        acc = (logits.argmax(dim=1) == y).float().mean()
        
        self.log('test_loss', loss)
        self.log('test_acc', acc)
        return loss
    
    def configure_optimizers(self):
        optimizer = torch.optim.AdamW(
            self.parameters(), 
            lr=self.learning_rate,
            weight_decay=0.01  # Add weight decay (L2 regularization)
        )
        return optimizer

# Prepare data
print("Preparing data for PyTorch Lightning...")

# Convert to PyTorch tensors
X_train_tensor = torch.FloatTensor(X_train_cifar)
y_train_tensor = torch.LongTensor(y_train_cifar)
X_test_tensor = torch.FloatTensor(X_test_cifar)
y_test_tensor = torch.LongTensor(y_test_cifar)

# Create datasets
train_dataset = TensorDataset(X_train_tensor, y_train_tensor)
test_dataset = TensorDataset(X_test_tensor, y_test_tensor)

# Create data loaders
train_loader = DataLoader(train_dataset, batch_size=256, shuffle=True, num_workers=4)
test_loader = DataLoader(test_dataset, batch_size=256, shuffle=False, num_workers=4)

# Initialize model
model = CIFAR10Classifier(input_dim=3072, hidden_dim=128, num_classes=10, learning_rate=0.001)

# Set up trainer
trainer = pl.Trainer(
    max_epochs=50,
    accelerator='gpu' if torch.cuda.is_available() else 'cpu',
    devices=1,
    callbacks=[
        ModelCheckpoint(
            monitor='val_acc',
            mode='max',
            filename='cifar10-{epoch:02d}-{val_acc:.4f}'
        )
    ],
    log_every_n_steps=50
)

# Train the model
print("\nTraining neural network on CIFAR-10...")
print(f"Using device: {'GPU' if torch.cuda.is_available() else 'CPU'}")
trainer.fit(model, train_loader, test_loader)

# Test the model
print("\nTesting the model...")
test_results = trainer.test(model, test_loader)

print(f"\nFinal Test Accuracy: {test_results[0]['test_acc']:.4f}")
print(f"Final Test Loss: {test_results[0]['test_loss']:.4f}")

GPU available: True (cuda), used: True
TPU available: False, using: 0 TPU cores
HPU available: False, using: 0 HPUs
You are using a CUDA device ('NVIDIA GeForce RTX 5090') that has Tensor Cores. To properly utilize them, you should set `torch.set_float32_matmul_precision('medium' | 'high')` which will trade-off precision for performance. For more details, read https://pytorch.org/docs/stable/generated/torch.set_float32_matmul_precision.html#torch.set_float32_matmul_precision
2026-01-29 12:57:34.818428: I tensorflow/core/util/port.cc:113] oneDNN custom operations are on. You may see slightly different numerical results due to floating-point round-off errors from different computation orders. To turn them off, set the environment variable `TF_ENABLE_ONEDNN_OPTS=0`.
2026-01-29 12:57:34.832501: E external/local_xla/xla/stream_executor/cuda/cuda_dnn.cc:9261] Unable to register cuDNN factory: Attempting to register factory for plugin cuDNN when one has already been registered
2026-01-29 12:57:34.832520: E external/local_xla/xla/stream_executor/cuda/cuda_fft.cc:607] Unable to register cuFFT factory: Attempting to register factory for plugin cuFFT when one has already been registered
2026-01-29 12:57:34.833045: E external/local_xla/xla/stream_executor/cuda/cuda_blas.cc:1515] Unable to register cuBLAS factory: Attempting to register factory for plugin cuBLAS when one has already been registered
2026-01-29 12:57:34.835982: I tensorflow/core/platform/cpu_feature_guard.cc:182] This TensorFlow binary is optimized to use available CPU instructions in performance-critical operations.
To enable the following instructions: AVX2 AVX512F AVX512_VNNI AVX512_BF16 AVX_VNNI FMA, in other operations, rebuild TensorFlow with the appropriate compiler flags.
2026-01-29 12:57:35.321123: W tensorflow/compiler/tf2tensorrt/utils/py_utils.cc:38] TF-TRT Warning: Could not find TensorRT
LOCAL_RANK: 0 - CUDA_VISIBLE_DEVICES: [0]

  | Name    | Type        | Params | Mode 
------------------------------------------------
0 | fc1     | Linear      | 3.1 M  | train
1 | bn1     | BatchNorm1d | 2.0 K  | train
2 | fc2     | Linear      | 524 K  | train
3 | bn2     | BatchNorm1d | 1.0 K  | train
4 | fc3     | Linear      | 131 K  | train
5 | bn3     | BatchNorm1d | 512    | train
6 | fc5     | Linear      | 2.6 K  | train
7 | dropout | Dropout     | 0      | train
------------------------------------------------
3.8 M     Trainable params
0         Non-trainable params
3.8 M     Total params
15.236    Total estimated model params size (MB)
/home/kmcalist/.local/lib/python3.10/site-packages/pytorch_lightning/utilities/_pytree.py:21: `isinstance(treespec, LeafSpec)` is deprecated, use `isinstance(treespec, TreeSpec) and treespec.is_leaf()` instead.
`Trainer.fit` stopped: `max_epochs=50` reached.
LOCAL_RANK: 0 - CUDA_VISIBLE_DEVICES: [0]

Preparing data for PyTorch Lightning...

Training neural network on CIFAR-10...
Using device: GPU

Testing the model...

Final Test Accuracy: 0.5739
Final Test Loss: 2.3761

┏━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
┃        Test metric        ┃       DataLoader 0        ┃
┡━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━┩
│         test_acc          │    0.5738999843597412     │
│         test_loss         │    2.3761208057403564     │
└───────────────────────────┴───────────────────────────┘

Code

import os
from pytorch_lightning.callbacks import ModelCheckpoint
import glob

# Load the best model checkpoint

# Find the best checkpoint file
checkpoint_files = glob.glob('lightning_logs/version_0/checkpoints/*.ckpt')

if checkpoint_files:
    # Get the most recent checkpoint (or you can sort by validation accuracy)
    best_checkpoint = max(checkpoint_files, key=lambda x: os.path.getctime(x))
    
    print(f"Loading checkpoint: {best_checkpoint}")
    
    # Load the model from checkpoint
    loaded_model = CIFAR10Classifier.load_from_checkpoint(best_checkpoint)
    
    # Evaluate on test set
    trainer_eval = pl.Trainer(
        accelerator='gpu' if torch.cuda.is_available() else 'cpu',
        devices=1
    )
    
    results = trainer_eval.test(loaded_model, test_loader)
    
    print(f"\nValidation (Test) Accuracy: {results[0]['test_acc']:.4f}")
    print(f"Validation (Test) Loss: {results[0]['test_loss']:.4f}")
else:
    print("No checkpoint files found. Using the last trained model.")
    print(f"\nValidation (Test) Accuracy: {test_results[0]['test_acc']:.4f}")
    print(f"Validation (Test) Loss: {test_results[0]['test_loss']:.4f}")

Loading checkpoint: lightning_logs/version_0/checkpoints/cifar10-epoch=46-val_acc=0.5788.ckpt

Validation (Test) Accuracy: 0.5788
Validation (Test) Loss: 2.3064

GPU available: True (cuda), used: True
TPU available: False, using: 0 TPU cores
HPU available: False, using: 0 HPUs
LOCAL_RANK: 0 - CUDA_VISIBLE_DEVICES: [0]

┏━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
┃        Test metric        ┃       DataLoader 0        ┃
┡━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━┩
│         test_acc          │    0.5788000226020813     │
│         test_loss         │    2.3063502311706543     │
└───────────────────────────┴───────────────────────────┘