任务：手动实现PCA算法并应用于手写数字数据集
一、简介
手写数字数据集（MNIST dataset）是深度学习领域中一个非常经典的数据集，包含了70,000张28×28像素的手写数字图像，每个图像有784个特征（每个像素点的灰度值），并且分为10个类别（0-9）。PCA可以用于降维，减少计算复杂度，同时还可以用于数据可视化。
二、步骤
1. 数据加载与预处理
· 使用sklearn.datasets加载MNIST数据集。
· 将数据集分为特征矩阵X和目标标签y。
· 对特征矩阵X进行归一化处理（将像素值缩放到[0, 1]区间）。
2. 手动实现PCA算法
· 编写一个函数manual_pca(X, n_components)，手动实现PCA算法，将数据降维到指定的维度n_components。
· PCA算法步骤：
1. 对数据进行中心化（减去均值）。
2. 计算协方差矩阵。
3. 对协方差矩阵进行特征分解，获取特征值和特征向量。
4. 选择最大的n_components个特征值对应的特征向量，组成降维矩阵。
5. 使用降维矩阵对原始数据进行降维。
3. 降维与可视化
· 使用手动实现的PCA函数将MNIST数据集降维到2维。
· 使用matplotlib库绘制降维后的数据点，并根据目标标签y对不同类别的数据点使用不同的颜色进行区分。
· 添加合适的标题、坐标轴标签和图例，使图形清晰易懂。
4. 性能对比
· 使用sklearn.decomposition.PCA类对MNIST数据集进行降维到2维。
· 比较手动实现的PCA和sklearn的PCA在降维结果上的差异（可以通过可视化对比）。
· 分析两种方法的运行时间和内存占用情况。
5. 解释方差
· 打印手动实现的PCA模型的解释方差比率（explained_variance_ratio_），分析前两个主成分所占的方差比例。
三、代码分析
3.1 数据加载函数
def load_mnist(path, kind='train'):
"""从idx文件加载MNIST数据"""
raw_path = os.path.join(path, "MNIST", "raw")
if not os.path.exists(raw_path):
raise FileNotFoundError(f"未找到raw目录：{raw_path}")
# 解析标签文件（.idx1格式）
labels_path = os.path.join(raw_path, f"{kind}-labels-idx1-ubyte")
with open(labels_path, 'rb') as f:
magic, n = struct.unpack('>II', f.read(8)) # 前4字节为魔数，后4字节为样本数
labels = np.fromfile(f, dtype=np.uint8) # 读取标签数据
# 解析图像文件（.idx3格式）
images_path = os.path.join(raw_path, f"{kind}-images-idx3-ubyte")
with open(images_path, 'rb') as f:
magic, num, rows, cols = struct.unpack('>IIII', f.read(16))
images = np.fromfile(f, dtype=np.uint8).reshape(num, rows*cols) # 展平为784维特征
return images, labels
3.2 数据预处理
def load_and_preprocess_data(data_root):
"""加载数据并归一化"""
print(f"正在从 {data_root} 加载MNIST数据...")
X_train, y_train = load_mnist(data_root, kind='train')
X_test, y_test = load_mnist(data_root, kind='t10k')
print(f"训练集大小: {X_train.shape[0]}样本，特征维度: {X_train.shape[1]}")
print(f"测试集大小: {X_test.shape[0]}样本")
# 归一化到[0, 1]区间
scaler = MinMaxScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
return X_train_scaled, X_test_scaled, y_train, y_test
3.3 手动实现 PCA（核心步骤）
def manual_pca(X, n_components):
"""手动实现PCA算法"""
# 1. 数据中心化
X_mean = np.mean(X, axis=0)
X_centered = X - X_mean
# 2. 计算协方差矩阵（样本数多时有优化空间）
cov_matrix = np.cov(X_centered, rowvar=False) # rowvar=False表示样本为行，特征为列
# 3. 特征值分解（对称矩阵用eigh效率更高）
eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)
# 4. 按特征值降序排序
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]
# 5. 选择主成分并降维
top_eigenvectors = eigenvectors[:, :n_components] # 前k个特征向量
X_reduced = np.dot(X_centered, top_eigenvectors) # 投影到低维空间
# 计算解释方差比率
total_variance = np.sum(eigenvalues)
explained_variance_ratio = eigenvalues[:n_components] / total_variance
return X_reduced, explained_variance_ratio
3.4 主函数逻辑
def main():
# 1. 数据加载与预处理
X_train, X_test, y_train, y_test = load_and_preprocess_data(data_root)
# 2. 选取部分样本加速实验（MNIST样本量较大）
sample_size = 5000
X_sample = X_train[:sample_size]
y_sample = y_train[:sample_size]
print(f"\n使用 {sample_size} 个样本进行分析")
# 3. 手动PCA实验
manual_time, manual_mem, X_manual, ev_manual = run_pca(manual_pca, X_sample, 2)
# 4. sklearn PCA实验
sklearn_time, sklearn_mem, X_sklearn, ev_sklearn = run_pca(sklearn_pca, X_sample, 2)
# 5. 结果可视化与对比
plot_comparison(X_manual, X_sklearn, y_sample)
plot_explained_variance(ev_manual, ev_sklearn)
3.5 性能统计函数
def run_pca(pca_func, X, n_components):
"""运行PCA并统计时间和内存"""
start_time = time.time()
mem_before = psutil.Process(os.getpid()).memory_info().rss / 1024**2 # MB
X_reduced, ev_ratio = pca_func(X, n_components)
mem_after = psutil.Process(os.getpid()).memory_info().rss / 1024**2
run_time = time.time() - start_time
mem_usage = mem_after - mem_before
print(f"运行时间: {run_time:.4f}秒 | 内存占用: {mem_usage:.2f} MB")
print(f"解释方差比率: {ev_ratio.round(4)}")
print(f"前两主成分方差和: {np.sum(ev_ratio):.4f}\n")
return run_time, mem_usage, X_reduced, ev_ratio
四、完整代码
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import MinMaxScaler
from sklearn.decomposition import PCA
import time
import psutil
import os
from struct import unpack
# 设置中文显示
plt.rcParams["font.family"] = ["Microsoft YaHei", "SimHei"]
plt.rcParams["axes.unicode_minus"] = False # 解决负号显示问题
def load_mnist(path, kind='train'):
"""
从指定路径的raw子目录加载MNIST数据（idx格式）
参数:
path: 数据集根目录（包含raw子目录）
kind: 'train'或't10k'，分别表示训练集和测试集
返回:
images: 图像数据，形状为(n_samples, 784)
labels: 标签数据，形状为(n_samples,)
"""
raw_path = os.path.join(path, "MNIST", "raw") # 指向raw子目录
if not os.path.exists(raw_path):
raise FileNotFoundError(f"未找到raw目录：{raw_path}")
labels_path = os.path.join(raw_path, f"{kind}-labels-idx1-ubyte")
images_path = os.path.join(raw_path, f"{kind}-images-idx3-ubyte")
# 读取标签文件
with open(labels_path, 'rb') as lbpath:
magic, n = unpack('>II', lbpath.read(8)) # 解析文件头
labels = np.fromfile(lbpath, dtype=np.uint8) # 读取标签数据
# 读取图像文件
with open(images_path, 'rb') as imgpath:
magic, num, rows, cols = unpack('>IIII', imgpath.read(16)) # 解析文件头
images = np.fromfile(imgpath, dtype=np.uint8).reshape(len(labels), 784) # 展平为784维特征
return images, labels
def load_and_preprocess_data(data_root):
"""加载并预处理MNIST数据集"""
print(f"正在从 {data_root} 加载MNIST数据...")
X_train, y_train = load_mnist(data_root, kind='train')
X_test, y_test = load_mnist(data_root, kind='t10k')
print(f"训练集样本数: {X_train.shape[0]}, 特征维度: {X_train.shape[1]}")
print(f"测试集样本数: {X_test.shape[0]}")
# 数据归一化（缩放到[0, 1]）
scaler = MinMaxScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
return X_train_scaled, X_test_scaled, y_train, y_test
def manual_pca(X, n_components):
"""手动实现PCA算法（同前）"""
X_mean = np.mean(X, axis=0)
X_centered = X - X_mean
cov_matrix = np.cov(X_centered, rowvar=False)
eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]
top_eigenvectors = eigenvectors[:, :n_components]
X_reduced = np.dot(X_centered, top_eigenvectors)
total_variance = np.sum(eigenvalues)
explained_variance_ratio = eigenvalues[:n_components] / total_variance
return X_reduced, explained_variance_ratio
def plot_pca_results(X_reduced, y, title, ax=None):
"""可视化降维结果（同前）"""
if ax is None:
fig, ax = plt.subplots(figsize=(10, 8))
y_int = y.astype(int)
colors = plt.cm.tab10.colors
for digit in range(10):
mask = (y_int == digit)
ax.scatter(X_reduced[mask, 0], X_reduced[mask, 1],
label=str(digit), alpha=0.6, s=10, c=[colors[digit]])
ax.set_title(title, fontsize=14)
ax.set_xlabel("第一主成分", fontsize=12)
ax.set_ylabel("第二主成分", fontsize=12)
ax.legend(title="数字标签", loc='best', frameon=True, markerscale=2)
ax.grid(True, linestyle='--', alpha=0.7)
return ax
def get_memory_usage():
"""获取内存占用（同前）"""
process = psutil.Process(os.getpid())
return process.memory_info().rss / 1024 / 1024 # MB
def main():
# 你的数据集根目录（包含MNIST/raw子目录）
data_root = r""
# 加载数据并预处理
X_train, X_test, y_train, y_test = load_and_preprocess_data(data_root)
# 取部分样本加速计算（MNIST训练集共60000样本，取5000测试）
sample_size = 5000
X_sample = X_train[:sample_size]
y_sample = y_train[:sample_size]
print(f"\n使用 {sample_size} 个训练样本进行分析")
# ---------------------- 手动PCA ----------------------
print("\n---------- 手动实现PCA ----------")
start_time = time.time()
mem_before = get_memory_usage()
X_reduced_manual, explained_var_manual = manual_pca(X_sample, n_components=2)
manual_time = time.time() - start_time
manual_mem = get_memory_usage() - mem_before
print(f"运行时间: {manual_time:.4f} 秒 | 内存占用: {manual_mem:.2f} MB")
print(f"解释方差比率: {explained_var_manual.round(4)}")
print(f"前两主成分方差和: {np.sum(explained_var_manual):.4f}")
# ---------------------- sklearn PCA ----------------------
print("\n---------- sklearn实现PCA ----------")
start_time = time.time()
mem_before = get_memory_usage()
pca_sklearn = PCA(n_components=2)
X_reduced_sklearn = pca_sklearn.fit_transform(X_sample)
sklearn_time = time.time() - start_time
sklearn_mem = get_memory_usage() - mem_before
print(f"运行时间: {sklearn_time:.4f} 秒 | 内存占用: {sklearn_mem:.2f} MB")
print(f"解释方差比率: {pca_sklearn.explained_variance_ratio_.round(4)}")
print(f"前两主成分方差和: {np.sum(pca_sklearn.explained_variance_ratio_):.4f}")
# ---------------------- 可视化对比 ----------------------
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 8))
plot_pca_results(X_reduced_manual, y_sample, "手动PCA降维结果", ax1)
plot_pca_results(X_reduced_sklearn, y_sample, "sklearn PCA降维结果", ax2)
plt.suptitle("PCA降维效果对比", fontsize=16, y=1.02)
plt.tight_layout()
plt.savefig("pca_comparison.png", dpi=300)
plt.show()
# ---------------------- 解释方差对比图 ----------------------
plt.figure(figsize=(10, 6))
bars = plt.bar(
["手动PC1", "手动PC2", "sklearn PC1", "sklearn PC2"],
[explained_var_manual[0], explained_var_manual[1],
pca_sklearn.explained_variance_ratio_[0], pca_sklearn.explained_variance_ratio_[1]],
color=["#1f77b4", "#aec7e8", "#ff7f0e", "#ffbb78"],
width=0.4
)
for bar in bars:
height = bar.get_height()
plt.text(bar.get_x() + bar.get_width()/2, height, f"{height:.4f}", ha="center", va="bottom")
plt.title("主成分解释方差比率", fontsize=14)
plt.ylabel("解释方差比例", fontsize=12)
plt.grid(axis="y", linestyle="--", alpha=0.7)
plt.tight_layout()
plt.savefig("explained_variance.png", dpi=300)
plt.show()
if __name__ == "__main__":
main()
