上周整理一批实验图像，几百MB的文件堆在硬盘里，磁盘眼看就要见底。换格式、压分辨率多多少少会伤画面细节，我临时决定用 Matlab 跑一把主成分分析（PCA）来做图像压缩和重建。这个思路不算新，但胜在原理清晰、实现门槛低：把一张灰度图的像素矩阵投影到少数几个主要方向上，丢掉贡献微弱的分量，再用这些主成分把图像拼回来。如果你正在学 PCA，或者手头有大量灰度图想要快速压缩，这篇文章可以当一份完整参考。
我尽量把从数学原理到代码落地、从踩坑到调优的整个流程都讲透，不绕弯子。实际动手时你会发现，PCA 图像压缩真正难的不是公式，而是矩阵的尺寸、特征分解的数值问题、以及压缩率和画质之间的取舍。这些坑我都会提到。
1. 先搞清楚：PCA 凭什么能压缩图像
1.1 把图像看成数据矩阵
一张灰度图像，本质就是一个m×n的数值矩阵，每个数字代表对应位置的亮度值。常规图像有个共同特点：相邻像素之间高度相关，行与行、列与列之间往往长得差不多。换句话说，像素矩阵里存在大量"冗余信息"。
PCA 做压缩的思路，是给这个矩阵换一组新的坐标系。原本的坐标轴是"横着第几列、竖着第几行"，主成分分析找出另一组正交基，让图像数据在新基上的投影方差尽可能大。如果绝大多数像素信息能集中在少数几个新坐标轴上，那只需要记录这几个轴上的投影系数，就能近似还原整张图。
我习惯用一个衣柜的类比来解释：一堆衣服平铺在箱子里很占地方，但挂到衣架上之后，衣柜里留出的空间一下大很多。PCA 就是在找那个最合适的"挂衣方向"——你不需要保留每件衣服的完整形状，只需要记录它在衣架上的位置。
1.2 PCA 找的是方差最大的方向
为什么选方差最大的方向？因为方差代表信息量。把数据中心化之后，某个方向上的投影方差越大，说明数据在这个方向上铺得越开，丢失的信息越少。
稍微具体一点：
- 图像矩阵
X的尺寸是m×n，可以看成m个样本点，每个样本点有n个特征； - 先算出每个特征列的均值
meanV，做中心化得到X0 = X - meanV； - 计算协方差矩阵
C = X0' * X0 / (m-1)； - 对
C做特征分解，得到特征值 λ 和对应的特征向量 v； - 特征值越大，说明该特征向量方向上的方差越大，也就是这张图上最重要的信息方向。
把特征值按从大到小排序，取前 K 个对应的特征向量组成变换矩阵Vk。因为自然图像的行列相关性很强，通常前几十个特征值就能贡献 95% 以上的信息量，丢掉后面那些特征值小的方向，视觉上几乎察觉不到差别。
用数学语言表达就是：
Y = X0 * Vk # 投影，得到 K 列主成分系数 X_rec = Y * Vk' + meanV # 重建Y是压缩后的核心数据，尺寸是m×K，远小于原始的m×n。
1.3 整条压缩与重建链路
从读图到重建，完整流程可以拆成八步：
- 读取图像并转成
double类型； - 计算每一列的均值，做中心化；
- 计算协方差矩阵；
- 特征分解，按特征值降序排列；
- 根据累积贡献率选择前 K 个特征向量；
- 投影得到主成分系数；
- 用主成分系数和特征向量重建图像；
- 计算 MSE、PSNR 等指标评估效果。
这八步缺一不可，其中第 2、4、5 步最容易出错，后面我会详细展开。
2. Matlab 实现：核心代码与逐步拆解
2.1 图像读取与数据预处理：uint8 转 double 的坑
Matlab 里读图像最常用imread，但读进来的灰度图是uint8类型，数值范围 0~255。如果直接拿去做协方差计算和矩阵乘法，结果会出现各种匪夷所思的问题：负数被截断、精度丢失、甚至直接报错。所以第一步一定是转成double，并且归一化到 0~1 范围。
clear; close all; clc; % 1. 读取图像并预处理 img = imread('test_gray.png'); % 换成你自己的灰度图 imgD = im2double(img); % uint8 -> double，并归一到[0,1] [m, n] = size(imgD);im2double这个函数会自动把uint8的 0~255 映射到 0~1，不需要自己手动除 255。如果图像是uint16类型，它也能正确处理。我自己刚开始用的时候图省事直接double(img)，结果后续所有的矩阵运算数值都偏大了一个量级，特征值大到离谱，重建图像一片雪花，查了半天才找到原因。
2.2 协方差矩阵与特征分解：为什么别忘了排序
中心化和协方差矩阵的计算如下：
% 2. 中心化：每一列减去该列的均值 meanVec = mean(imgD, 1); % 1×n 的均值向量 X0 = imgD - meanVec; % 去均值后的矩阵 % 3. 协方差矩阵 C = (X0' * X0) / (m - 1); % n×n 协方差矩阵 % 4. 特征分解 [V, D] = eig(C); % V 的每一列是特征向量，D 对角线是特征值 % 5. 按特征值降序重排 [eigVals, idx] = sort(diag(D), 'descend'); V = V(:, idx);这一段我踩过最大的坑是eig的输出顺序。eig返回的特征值并不保证按大小排列，如果你直接取V(:,1:K)，大概率取到的不是最重要的主成分，重建出来图像是一团噪声。所以必须用sort把特征值按降序排列，同时把特征向量的列一起重排。V = V(:, idx)这一步不能省，很多人就是漏了这一行导致结果完全不对。
还有一个小细节：协方差矩阵理论上应该用X0' * X0 / (m-1)，因为样本方差的分母是m-1而不是m。用cov(X0)也可以得到相同结果，不过cov默认会把每列当成一个变量，返回的也是n×n矩阵。两种写法都行，但自写公式更直观，方便调试。
2.3 投影、压缩存储与重建：代码闭环
选完特征向量之后，就可以做投影和重建了：
% 6. 累积贡献率，自动选择 K cumRatio = cumsum(eigVals) / sum(eigVals); K = find(cumRatio >= 0.99, 1, 'first'); if isempty(K); K = n; end % 7. 投影与重建 Vk = V(:, 1:K); Y = X0 * Vk; % m×K 主成分系数 XRec = Y * Vk' + meanVec; % 重建图像 % 8. 误差评估 MSE = mean((imgD(:) - XRec(:)).^2); PSNR = 10 * log10(1 / MSE); fprintf('K=%d, MSE=%.6f, PSNR=%.2f dB\n', K, MSE, PSNR); % 压缩率估算（未量化的 double 存储） originalBytes = m * n; compressedBytes = m * K + K * n + n; ratio = compressedBytes / originalBytes; fprintf('压缩后存储占比: %.2f%%\n', ratio * 100); % 9. 可视化 figure; subplot(1,3,1); imshow(img); title('原图'); subplot(1,3,2); imshow(im2uint8(XRec)); title('PCA重建'); subplot(1,3,3); imshow(im2uint8(abs(imgD - XRec))); title('残差');这里有几个点需要解释。
首先是压缩率。原本存储一张灰度图需要m×n个字节（uint8）；压缩后需要保存三样东西：主成分系数Y（m×K）、特征向量Vk（K×n）、以及均值向量meanVec（1×n）。所以压缩后存储量大约等于K×(m+n)。K 越小压缩率越高，但重建效果越差。
其次是存储格式。上面代码里的Y和Vk都是double，如果直接存成double文件，不仅没有压缩，反而会把原始uint8图像撑大 8 倍。工程上正确的做法是把Y和Vk做量化，比如转成uint16甚至uint8保存，再配合无损压缩，才能真正减小存储空间。这一点在测试代码里不体现，但实际项目中一定要想明白，不然你会遇到"压缩完文件反而变大"的尴尬情况。
3. 压缩率与重建质量的两全之道
3.1 主成分数量 K 怎么定：累积贡献率法
K 的选择是整个流程中最关键的决策。K 太小，重建图像模糊；K 太大，压缩效果变差。最常用的方法是累积贡献率法：
累积贡献率 = 前 K 个特征值之和 / 全部特征值之和一般取 90%、95% 或 99% 作为阈值。比如阈值设成 99%，代码里就是find(cumRatio >= 0.99, 1, 'first')。
我拿一张 256×256 的测试图跑过一次，特征值的分布大概如下：
| K值 | 累积贡献率 | 压缩后存储占比 | PSNR (dB) | 主观效果 |
|---|---|---|---|---|
| 10 | 88.7% | 8.2% | 27.6 | 轮廓清楚，细节模糊 |
| 30 | 95.3% | 23.8% | 33.2 | 主要纹理可辨 |
| 60 | 97.9% | 47.3% | 37.9 | 肉眼基本接近原图 |
| 100 | 99.2% | 78.5% | 42.8 | 几乎看不出差异 |
压缩后存储占比的计算公式是(m*K + K*n + n)/(m*n)，假设 m=n=256。这张表说明一个道理：图像的信息高度集中在前几十个主成分里，K 从 10 增加到 100 时，PSNR 只提升了约 15dB，但存储量翻了将近十倍。实际项目中不需要追求 99% 以上的贡献率，95% 左右往往是性价比最高的甜点区。
注意：上面表格里的数字来自我某一次测试的结果，不同图像曲线差异很大。动手之前先画出特征值下降曲线，看看信息到底集中在几个方向上，再定阈值，远比照搬经验值靠谱。
3.2 大图怎么办：分块 PCA 的工程解法
整图做 PCA 在图像比较小的时候没问题，但一旦图像变成 1024×1024 甚至更大，协方差矩阵就是1024×1024甚至2048×2048，特征分解的时间会指数级上升，内存也扛不住。
工程上最常见的解法是分块 PCA。把图像切成若干个固定大小的小块，比如32×32或64×64，然后对每一块单独做 PCA。分块的好处有两个：
- 块的大小可控，协方差矩阵的尺寸固定为
块宽×块高，计算量不会随整图尺寸爆炸； - 图像的局部相关性更强，小块内像素之间往往更加相似，PCA 压缩效率反而更高。
分块的流程大致是：
- 用
mat2cell或索引把图像切成若干小块； - 将每个小块拉直成一个向量；
- 对所有这些向量做 PCA，或者对每一块单独做 PCA；
- 保存压缩数据时，需要额外保存每一块的均值向量和特征向量。
块大小的选择是个取舍。块越小，局部细节保留得越好，但需要保存更多的均值向量和特征向量，压缩率下降；块越大，压缩率更高，但容易出现局部细节糊掉的情况。我用下来64×64是个不错的起点，既不会让协方差矩阵太大，又保留了足够的局部信息。
3.3 彩色图像：RGB 逐通道处理
灰度图处理完了，自然会想处理彩色图。彩色图读进来是m×n×3的三维矩阵，三个通道分别是 R、G、B。处理方法无非两种：
- 逐通道做 PCA：把 R、G、B 三个通道分别当作三张灰度图，各自做 PCA，最后把重建结果拼回去。实现简单，但三个通道之间的冗余信息没有被利用。
- 转 YCbCr 色彩空间：人类视觉对亮度 Y 分量最敏感，对色度 Cb、Cr 不够敏感。转到 YCbCr 后，可以对 Y 通道保留更多主成分（比如 95% 贡献率），Cb 和 Cr 通道保留更少主成分（比如 85%）。这样可以在同等存储下取得更好的主观画质。
Matlab 里用rgb2ycbcr转换，重建后再用ycbcr2rgb转回去。我实际测试下来，YCbCr 方案比 RGB 逐通道方案在同等压缩率下 PSNR 高 1~2dB，值得一试。
imgRGB = imread('test_color.png'); imgYCbCr = rgb2ycbcr(imgRGB); YCh = imgYCbCr(:,:,1); CbCh = imgYCbCr(:,:,2); CrCh = imgYCbCr(:,:,3); % 对三个通道分别做 PCA，Y 通道用高阈值，Cb/Cr 用低阈值4. 踩坑实录：常见问题与排查技巧
4.1 特征向量符号翻转与排序错位
特征分解有个让人头疼的性质：特征向量的符号不唯一。如果特征向量 v 是解，那么 -v 也同样是解。这意味着你两次运行同一段代码，得到的V矩阵列向量可能符号完全不同。这在数学上没问题，因为重建公式里Y和Vk是成对出现的，如果Vk翻转，Y也会跟着翻转，最终的Y * Vk'不变。
真正危险的操作是单独保存Vk和Y，然后放在两个不同的流程里使用时，偶尔会碰上符号混用。解决方法也简单：重建完成后检查一下与原图的方向一致性，或者干脆不单独处理 Y，重建时直接XRec = Y * Vk' + meanVec，不要手动拆开再组合。
排序错位的问题在前文提过，这里再强调一次：eig返回的特征值不保证降序，如果先[V, D] = eig(C)再直接用V(:, 1:K)，结果完全随机。一定要先排序，再把特征向量跟着一起重排。
4.2 协方差矩阵太大：用 SVD 替代特征分解
当图像尺寸很大时，直接算协方差矩阵再特征分解，内存和时间都很难接受。比如一张4096×4096的灰度图，协方差矩阵是4096×4096，一个 double 矩阵就要占 128MB，再乘上特征分解的临时开销，直接能把普通电脑压垮。
更好的方案是绕开协方差矩阵，直接用奇异值分解 SVD：
[U, S, V] = svd(X0, 'econ'); % 原图中心化矩阵 X0 的尺寸是 m×n % U 是 m×m，S 是 m×n 对角阵，V 是 n×n % 主成分系数 Y = U(:, 1:K) * S(1:K, 1:K)SVD 在数值稳定性上比"先算协方差再特征分解"更好，因为协方差矩阵的计算本身会引入舍入误差。更重要的是，SVD 可以直接对原始数据中心化矩阵操作，不需要显式构造X0' * X0，节省了大量内存。
提示：SVD 和特征分解之间在数学上是等价的。
X0的右奇异向量就是X0'*X0的特征向量，奇异值的平方除以m-1就是特征值。所以如果你已经理解了特征分解的流程，SVD 版本只是换了个更稳的算法，思路完全一样。
4.3 重建图发灰、发暗：均值向量忘加了
这是我见过最多人翻车的错误。PCA 压缩时做了中心化，重建时忘了把均值向量加回去，结果图像整体发灰、发暗，像盖了一层滤镜。代码里XRec = Y * Vk' + meanVec最后那个+ meanVec是必须的，不能丢。
另一种类似问题是归一化方向搞反了。im2double把图像从 0~255 映射到 0~1，如果重建后忘记先转回原始范围，直接用imshow显示，画面会黑到什么都看不见。正确做法是用im2uint8(XRec)转回 uint8 再显示，同时注意im2uint8会自动剪裁到 0~255 范围，所以残差图显示前最好加个abs()，再乘个系数放大，不然细节看不清楚。
4.4 评估指标选择：MSE、PSNR 与主观效果
评估指标方面，最常用的是 MSE 和 PSNR：
MSE = mean((imgD(:) - XRec(:)).^2); PSNR = 10 * log10(1 / MSE); % 因为图像已归一到[0,1]，MAX=1如果图像没有归一化到 0~1，PSNR 公式里的 MAX 要选对。uint8 图像 MAX=255，double 归一化图像 MAX=1。很多人直接拿别人代码里的公式跑，结果 PSNR 算出来偏高或偏低，就是因为这个细节。
不过 PSNR 不是万能指标。它衡量的是像素级的误差，但人眼有时候对局部纹理的失真非常敏感，有时候对整体亮度偏差反而不敏感。我见过一张图，PSNR 达到了 35dB，重建结果主观上看却非常"糊"；反过来也有 PSNR 只有 28dB 但人眼觉得不错的案例。所以评估时一定要主观+客观结合，尤其不能只看 PSNR 选 K。
4.5 其他几个容易踩的小坑
- 图像含透明通道或索引模式：
imread读入 PNG 有时会返回m×n×4（RGBA），或者返回索引图像。处理前一定要先size(img)确认维度，灰度处理只接受m×n，彩色图是m×n×3，遇到四通道要先去掉透明通道。 - 混用
eig和eigs：eigs是专门求前几个最大特征值的稀疏解法器，适合大矩阵，但它的速度和收敛性依赖于矩阵性质，小矩阵直接用eig反而更快。不要盲目迷信eigs，先评估矩阵规模再说。 - 分块边界效应：分块 PCA 重建后，相邻块之间可能出现亮度不连续，俗称"块效应"。缓解方法是块之间重叠几个像素，或者对重建后的块做加权融合。简单的平均融合就能大幅降低边界感。
我在实际项目里最后总结出的通用流程是：小图直接用整图 PCA +eig，大图用分块 +svd，K 的选择统一用累积贡献率 95% + 人工看一眼特征值曲线。先把流程在一张小图上跑通，再切到大图，几乎没翻过车。
最后再分享一个细节，很多人做完压缩就只顾着看重建图和 PSNR，却忘了观察特征值曲线本身。那个曲线其实非常直观：如果前几个特征值就占据了绝对优势，说明这张图信息冗余很大，压缩空间大；如果特征值下降得很平缓，说明图像接近随机噪声，PCA 压缩效果会很差。动手前先看一眼这条曲线，你对这个项目能压到什么程度就有了大概的预期。这一点我觉得比任何自动选 K 的公式都更有用。
