自编码器 Autoencoder 数据降维实验

问题

考虑把高维数据压缩到低维空间时，经典的做法是收集数据快照并使用 SVD 来获取最优的线性子空间，也就是主成分分析（PCA）的做法。这种方法需要假设实际数据就是落在一个低维线性子空间中，对于非线性的情况并不能很好地处理。

机器学习中的自编码器 Autoencoder 通过神经网络直接学习非线性映射，可以实现非线性降阶。

下面构造一个具有两个真实自由度、但观测维度为 100 的数据集：

$$ X\in \mathbb{R} ^{N\times d},\qquad N=2000,\quad d=100. $$

首先在 $[-3,3]^2$ 上均匀采样二维潜在变量 $z^*$。随后通过非线性映射生成 $\mathbb{R}^{100}$ 中的样本，映射中的矩阵 $A$、$B$ 的元素服从标准高斯分布，外面套一层三角函数并线性组合，最后加上小幅高斯噪声：

$$ X=\sin\left(z^* A^T\right) +0.5\cos\left(z^* B^T\right)+\varepsilon. $$

因此数据的真实自由度是 2，但它们通常不位于某个二维线性平面上，而是位于一个非线性流形上。

下面通过实验分别展示 PCA 和 Autoencoder 的降维和重构过程

$$ X \in \mathbb{R}^{100} \longrightarrow z \in \mathbb{R}^2 \longrightarrow \hat{X} \in \mathbb{R}^{100}, $$

并比较二维潜在空间和相对重构误差。相对重构误差的计算如下

$$ e_{\mathrm{rel}}=\frac{||X - \hat{X} ||_F}{||X ||_F}. $$

方法介绍

PCA 在给定维数 $r$ 时基于 SVD 分解寻找最优线性子空间。 设数据均值为 $\bar{x}$，前 $r$ 个主成分组成列正交矩阵 $W\in \mathbb{R}^{d\times r}$，压缩和解码过程分别为

$$ z=W^T(x-\bar{x}),\qquad \hat{x}=\bar{x}+Wz. $$

其中 $z\in \mathbb{R} ^r$ 是低维表示。矩阵 $W$ 使训练数据在线性编码和解码下的重构误差最小：

$$ \min_{W^TW=I}\frac{1}{N}\sum_{i=1}^{N} ||x_i-\bar{x}-WW^T(x_i-\bar{x})||_2^2. $$

这个问题是有显式解的，可以直接通过 SVD 分解得到结果。

Autoencoder 将 PCA 的线性编码器和解码器替换成神经网络 $E$ 和 $D$，压缩和解码过程分别为

$$ z=E(x),\qquad \hat{x}=D(z). $$

两个网络的参数通过最小化训练数据的重构误差得到：

$$ L=\frac{1}{N}\sum_{i=1}^{N} ||x_i-D\left(E(x_i)\right)||_2^2. $$

Autoencoder 基于神经网络提供的是非线性映射，而 PCA 的编码和解码都是线性映射。

代码演示

实验依赖 NumPy、Matplotlib、scikit-learn 和 PyTorch，可以在项目虚拟环境中安装：

uv pip install numpy matplotlib scikit-learn torch

完整代码如下：

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
import torch
import torch.nn as nn
from torch.utils.data import TensorDataset, DataLoader

import matplotlib.pyplot as plt

plt.rcParams.update(
    {
        "savefig.dpi": 300,
        "figure.autolayout": True,
        "text.usetex": True,
    }
)
# ============================================================
# 1. 构造非线性低维流形数据
# ============================================================

np.random.seed(0)

N = 2000  # snapshot 数量
d = 100  # 原始维度
r = 2  # 降维后的维度

z_true = np.random.uniform(-3, 3, (N, r))
A = np.random.randn(d, r)
B = np.random.randn(d, r)

X = np.sin(z_true @ A.T) + 0.5 * np.cos(z_true @ B.T) + 0.02 * np.random.randn(N, d)

print("data shape:", X.shape)


# ============================================================
# 2. PCA 降维和重构
# ============================================================

pca = PCA(n_components=r)
Z_pca = pca.fit_transform(X)
X_pca = pca.inverse_transform(Z_pca)

pca_error = np.linalg.norm(X - X_pca) / np.linalg.norm(X)
print("PCA relative error:", pca_error)


# ============================================================
# 3. Autoencoder 降维和重构
# ============================================================


class AutoEncoder(nn.Module):
    def __init__(self):
        super().__init__()
        self.encoder = nn.Sequential(
            nn.Linear(d, 64),
            nn.Tanh(),
            nn.Linear(64, 16),
            nn.Tanh(),
            nn.Linear(16, r),
        )
        self.decoder = nn.Sequential(
            nn.Linear(r, 16),
            nn.Tanh(),
            nn.Linear(16, 64),
            nn.Tanh(),
            nn.Linear(64, d),
        )

    def forward(self, x):
        z = self.encoder(x)
        return self.decoder(z)


device = "cuda" if torch.cuda.is_available() else "cpu"
model = AutoEncoder().to(device)

dataset = TensorDataset(torch.tensor(X, dtype=torch.float32))
loader = DataLoader(dataset, batch_size=128, shuffle=True)

optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)
loss_history = []

for epoch in range(300):
    total_loss = 0

    for (x,) in loader:
        x = x.to(device)
        xhat = model(x)
        loss = torch.mean((xhat - x) ** 2)

        optimizer.zero_grad()
        loss.backward()
        optimizer.step()

        total_loss += loss.item()

    loss_history.append(total_loss / len(loader))

    if epoch % 50 == 0:
        print(f"epoch {epoch}, loss {loss_history[-1]:.6e}")

with torch.no_grad():
    Xt = torch.tensor(X, dtype=torch.float32).to(device)
    X_ae = model(Xt).cpu().numpy()
    Z_ae = model.encoder(Xt).cpu().numpy()

ae_error = np.linalg.norm(X - X_ae) / np.linalg.norm(X)
print("AE relative error :", ae_error)


# ============================================================
# 4. 绘图
# ============================================================

plt.figure(figsize=(8, 4))
plt.plot(loss_history)
plt.yscale("log")
plt.xlabel("epoch")
plt.ylabel("MSE")
plt.title("Autoencoder training loss")
plt.grid()
plt.savefig("training-loss.png")

plt.figure(figsize=(10, 8))

plt.subplot(2, 2, 1)
plt.scatter(Z_pca[:, 0], Z_pca[:, 1], c=z_true[:, 0], s=5)
plt.title("PCA latent space ($z^*_1$)")
plt.xlabel("$z_1$")
plt.ylabel("$z_2$")

plt.subplot(2, 2, 2)
plt.scatter(Z_pca[:, 0], Z_pca[:, 1], c=z_true[:, 1], s=5)
plt.title("PCA latent space ($z^*_2$)")
plt.xlabel("$z_1$")
plt.ylabel("$z_2$")

plt.subplot(2, 2, 3)
plt.scatter(Z_ae[:, 0], Z_ae[:, 1], c=z_true[:, 0], s=5)
plt.title("Autoencoder latent space ($z^*_1$)")
plt.xlabel("$z_1$")
plt.ylabel("$z_2$")

plt.subplot(2, 2, 4)
plt.scatter(Z_ae[:, 0], Z_ae[:, 1], c=z_true[:, 1], s=5)
plt.title("Autoencoder latent space ($z^*_2$)")
plt.xlabel("$z_1$")
plt.ylabel("$z_2$")

plt.savefig("latent-space.png")

idx = np.random.choice(N, 6, replace=False)

plt.figure(figsize=(10, 8))

for i, j in enumerate(idx):
    plt.subplot(3, 2, i + 1)
    plt.plot(X[j], label="original")
    plt.plot(X_pca[j], "--", label="PCA")
    plt.plot(X_ae[j], ":", label="AE")
    plt.legend()

plt.savefig("reconstruction.png")

plt.figure(figsize=(10, 8))

for i, j in enumerate(idx):
    plt.subplot(3, 2, i + 1)
    plt.plot(X_pca[j] - X[j], label="PCA error")
    plt.plot(X_ae[j] - X[j], label="AE error")
    plt.legend()

plt.savefig("reconstruction-error.png")

实验结果

PCA 的重构误差：

PCA relative error: 0.777887741358002

Autoencoder 的重构误差：

AE relative error : 0.11560801506149158

Autoencoder 的训练损失持续下降，说明网络逐渐学到了数据的低维结构。

下面是潜在空间的可视化，其中的颜色依次表示真实变量 $z^*$ 的两个分量。

随机选取六个样本比较重构结果。PCA 只能用二维平面逼近原始曲线，因此丢失了较多细节；Autoencoder 的重构曲线与原始曲线更接近。

重构误差如下图。