0. Take-home message

通道压缩的过程是把原来的N个物理通道的数据，通过线性组合，变成N_c个虚拟通道的数据，我们希望在这个过程中，信息的损失尽可能小。这个线性组合，其实就是把原来的数据投影到新的基向量上，所以我们的问题变成了寻找最佳的基向量，使得投影之后，信息损失最少（i.e. 数据的方差最大）。而统计学知识告诉我们，协方差矩阵的最大特征向量指向最大方差方向（附录）。所以，新的基向量就是协方差矩阵的特征值比较大的特征向量。怎么求这个特征向量呢？实际上我们不用显式地计算协方差矩阵，只需要对原始数据进行SVD，然后取前N_c个特征向量即可。

1. 引言：我们的目标是什么？

在多通道磁共振成像（如并行成像）中，我们通常会从多个（例如 $N=32$ 个）接收线圈通道中同时采集数据。这会导致数据量非常大，给后续的图像重建和存储带来压力。

通道压缩 (Channel Compression) 的目标非常明确：将原始的 $N$ 个物理通道的数据，通过数学变换，合并成 $N_c$ 个"虚拟通道"（其中 $N_c < N$），同时尽可能多地保留原始信号的能量和信息。

我们的任务是找到一个最佳的变换，将一个维度为 (数据点数 x N) 的矩阵，高效地转换为一个 (数据点数 x N_c) 的矩阵。

2. 核心思想：SVD是实现PCA的强大工具

要找到"最佳"的变换，我们需要一个标准。这个标准就是主成分分析 (Principal Component Analysis, PCA)。

  • PCA的目标：PCA旨在找到一组新的正交基（即主成分或虚拟通道），这些基是原始通道的线性组合。数据在第一个主成分上的方差最大，在第二个主成分上的方差次之，以此类推。通过只保留前 $N_c$ 个最重要的主成分，我们就能在降维的同时最大化保留信息。

  • SVD的角色：那么如何找到这些主成分呢？奇异值分解 (Singular Value Decomposition, SVD) 就是实现PCA最高效、最稳健的数学工具。我们通过对一部分样本数据进行SVD，可以直接得到定义这些主成分的基向量。

简而言之：PCA是我们的"目标"（找到最佳的降维投影），而SVD是我们实现这个目标的"工具"

3. 分步实现：从SVD到通道压缩

我们通过一个四步流程来实现这个目标。

第1步：准备校准数据矩阵 $A$

我们从原始k空间数据中选取一小块区域，通常是k空间中心，这个区域信噪比较高，包含了线圈敏感度的主要信息。我们将这个区域的数据排列成一个矩阵 $A$，其维度为 $M \times N$。

  • $M$ 是校准区域内数据点的总数（例如 kx * ky）。
  • $N$ 是原始物理线圈的通道数。

第2步：对校准矩阵 $A$ 进行SVD

我们对矩阵 $A$ 进行奇异值分解： $A_{M \times N} = U_{M \times M} \Sigma_{M \times N} V^H_{N \times N}$ (在实践中，我们通常使用更高效的"经济型SVD"，分解后 $U$ 是 $M \times N$，$ \Sigma $ 是 $N \times N$ )

  • $U$: 左奇异向量矩阵。
  • $\Sigma$: 对角矩阵，对角线上是按从大到小排列的奇异值 $\sigma_i$。
  • $V^H$: 右奇异向量矩阵的共轭转置。$V$ 的列向量 $v_i$ 是我们最关心的部分。

如附录所示，$V$ 的列向量 $v_i$ 正是数据协方差矩阵 $A^H A$ 的特征向量，它们代表了数据的主成分方向，并且已经根据其重要性（由对应的奇异值 $\sigma_i$ 衡量）排好序。

第3步：构建压缩矩阵 $P$

这是从SVD到通道压缩最关键的桥梁。如果我们想将通道数从 $N$ 压缩到 $N_c$，我们只需选取矩阵 $V$ 的前 $N_c$ 列，来构成我们的压缩/投影矩阵 $P$。

$P = V_{:, 1:N_c}$

这个矩阵 $P$ 的维度是 $N \times N_c$。它的每一列都是一个"虚拟通道"的定义，该定义说明了如何将原始的 $N$ 个物理通道线性组合成一个新的虚拟通道。

第4步：应用压缩矩阵 $P$ (执行投影)

现在，我们可以使用这个压缩矩阵 $P$ 来压缩任意的k空间数据了（无论是校准数据还是整个k空间数据）。

假设我们有完整的k空间数据矩阵 $D$，其维度为 (总数据点数 x N)。压缩后的数据 $D_{comp}$ 通过简单的矩阵乘法得到：

$D_{comp} = D \cdot P$

我们来检查一下维度： (总数据点数 x N) $\times$ (N x N_c) $\rightarrow$ (总数据点数 x N_c)

这个乘法操作的几何意义是：将原始数据 $D$ 中的每一个数据点（一个在 $N$ 维原始通道空间中的向量），投影 (project) 到了由 $P$ 的列向量（即最重要的 $N_c$ 个主成分）所张成的新基空间上。

至此，我们成功地将通道数从 $N$ 减少到了 $N_c$，完成了通道压缩。

4. MATLAB代码示例

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55

% --- 1. Prepare Calibration Data A ---
% M: number of data points in calibration region
% N: number of original channels
M = 2000;
N = 32;
% Target number of compressed channels
N_c = 8;

% Simulate a complex-valued calibration matrix A
% In a real scenario, this comes from the k-space center
A = randn(M, N) + 1j * randn(M, N);
fprintf('Shape of calibration matrix A: %d x %d\n', size(A));

% --- 2. Perform SVD on A ---
% MATLAB's svd function returns U, S, V where V is already the right singular vectors
[U, S, V] = svd(A, 'econ');  % 'econ' for economy SVD

% U shape: (M, N), S shape: (N, N) diagonal matrix, V shape: (N, N)
% Extract singular values as a vector
s = diag(S);

% --- 3. Build the Compression Matrix P ---
% The compression matrix P consists of the first N_c columns of V
P = V(:, 1:N_c);
fprintf('Shape of compression matrix P: %d x %d\n', size(P));

% --- 4. Apply Compression (Projection) ---
% Simulate a full k-space dataset D from the same 32 channels
num_total_points = 256 * 256;
D_original = randn(num_total_points, N) + 1j * randn(num_total_points, N);
fprintf('Shape of original full data D: %d x %d\n', size(D_original));

% Perform compression by right-multiplying with P
% This is the projection step
D_compressed = D_original * P;
fprintf('Shape of compressed data D_comp: %d x %d\n', size(D_compressed));

% The number of channels is now successfully reduced from 32 to 8.

% Optional: Display singular values to understand compression quality
figure;
plot(1:length(s), s, 'bo-', 'LineWidth', 2);
hold on;
plot(1:N_c, s(1:N_c), 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r');
xlabel('Singular Value Index');
ylabel('Singular Value');
title('Singular Values and Compression Selection');
legend('All Singular Values', 'Selected for Compression', 'Location', 'best');
grid on;

% Calculate compression ratio and energy retention
compression_ratio = N_c / N;
energy_retention = sum(s(1:N_c).^2) / sum(s.^2);
fprintf('Compression ratio: %.2f%% (%d/%d channels)\n', compression_ratio*100, N_c, N);
fprintf('Energy retention: %.2f%%\n', energy_retention*100);

5. 结论

通过SVD，我们能够高效地计算出数据的主成分（存储在矩阵 $V$ 的列中）。通过选择最重要的前 $N_c$ 个主成分来构建一个投影矩阵 $P$，我们可以将原始的多通道数据投影到一个维度更低、但信息保留最完备的新空间中，从而实现通道压缩。这个过程的核心就是利用SVD作为工具来执行PCA。

附录：为什么协方差矩阵的最大特征向量指向最大方差方向？

这是一个经典的优化问题，可以通过拉格朗日乘子法求解。

  1. 问题设定

    • 我们有一个中心化（均值为0）的数据矩阵 $X$ ($m$ 个样本, $n$ 个特征)。
    • 我们想找到一个单位投影方向 $w$ ($w^T w = 1$)，使得数据投影到该方向后得到的向量 $z = Xw$ 的方差最大。

  2. 方差表达式 投影后数据的方差 $Var(z)$ 为： $Var(z) = \frac{1}{m-1} z^T z = \frac{1}{m-1} (Xw)^T(Xw) = w^T (\frac{1}{m-1}X^T X) w$

    我们注意到，中间的项 $C_X = \frac{1}{m-1}X^T X$ 正是数据 $X$ 的协方差矩阵。 所以，我们的目标函数是： $Maximize \ f(w) = w^T C_X w$

  3. 拉格朗日乘子法 我们构建拉格朗日函数 $\mathcal{L}$ 来整合约束条件 $w^T w = 1$： $\mathcal{L}(w, \lambda) = w^T C_X w - \lambda(w^T w - 1)$

    对 $w$ 求导并令其为0，以寻找极值点： $\frac{\partial \mathcal{L}}{\partial w} = 2C_X w - 2\lambda w = 0$

    这给出了核心方程： $C_X w = \lambda w$

  4. 最终结论 这个方程的解告诉我们，任何能使方差取得极值的方向 $w$ 都必须是协方差矩阵 $C_X$ 的一个特征向量，而 $\lambda$ 则是对应的特征值。

    那么，这个方向上的方差具体是多少呢？我们将 $C_X w = \lambda w$ 代回方差表达式： $Var(z) = w^T C_X w = w^T (\lambda w) = \lambda (w^T w)$

    因为 $w^T w = 1$，我们得到： $Var(z) = \lambda$

    这意味着，数据在某个特征向量方向上的方差，就等于其对应的特征值。因此，为了使方差最大化，我们必须选择对应最大特征值的那个特征向量。