Q1：矩阵 \boldsymbol{A} 和 \boldsymbol{X} 的尺寸
高光谱图像数据由一个三维矩阵 \boldsymbol{A} 编码，其维度为：
其中：
- h 是图像的高度（像素数）
- w 是图像的宽度（像素数）
- n_\lambda 是波长数量（光谱通道数）
高光谱数据的结构是三维立方体 f(x, y, \lambda)，其中 (x, y) 是空间坐标，\lambda 是波长（光谱维度）。沿空间维度切片得到单波长图像，沿光谱维度切片得到单像素的光谱曲线。这种数据结构使得我们既可以分析特定波长下的空间分布，也可以分析特定位置的光谱特征。
将三维矩阵展开为二维矩阵 \boldsymbol{X} 时，每个像素的光谱作为一行排列：
其中 n_{pix} = h \times w 是总像素数。
矩阵 \boldsymbol{X} 可以从两个角度理解：
- 从信号处理角度看，每一行是一个像素的光谱信号，共 n_{pix} 个信号，每个信号有 n_\lambda 个采样点
- 从几何角度看，每一行 x_i \in \mathbb{R}^{n_\lambda} 是高维空间中的一个点，n_{pix} 行对应 n_{pix} 个点构成的点云。
在实验代码中，我们对数据文件 Reslice_Stack_NAU2-60X-1.tif 进行加载，获得了矩阵的完整尺寸为: h = 249，w = 696，n_\lambda = 457。
综上，对于本实验数据：
从几何角度看，\boldsymbol{X} 可以理解为457维空间中173304个点构成的点云。
Q2：PCA降维及误差计算
主成分分析（PCA）通过奇异值分解（SVD）实现降维。对中心化后的数据矩阵 \boldsymbol{Z} 进行SVD分解：
其中 \boldsymbol{U} 和 \boldsymbol{V} 满足正交性质
- 矩阵 \boldsymbol{V} 的列向量 v_1, v_2, \ldots 就是主成分方向
- v_1：第一主方向（解释最大方差的方向）
- v_2：第二主方向（与 v_1 正交，解释次大方差）
投影计算：
对于第 i 个数据点 x_i（即 \boldsymbol{Z} 的第 i 行），投影到 v_1 上的坐标为
c_1 就是 x_i 在第一主轴上的坐标。
降维表示
保留前 p 个主成分的降维数据表示为：
其中 \boldsymbol{U}_p、\boldsymbol{\Sigma}_p、\boldsymbol{V}_p 分别是前 p 列/行的截断矩阵。PCA的目标是在高维空间 \mathbb{R}^{457} 中寻找一个低维子空间，使得所有数据点投影到该子空间后，投影误差最小。
确定保留多少个主成分的常用方法是分析奇异值的累积能量，即计算累积和占总和的比例。
重构误差（Frobenius范数）：
相对误差：
这个比值表示丢失的信息比例
逐点误差：
对于点 x_i，设它与投影方向的夹角为 \theta_i：
- 若 \cos^2 \theta_i = 1：点完全在投影轴/平面上，无误差
- 若 \cos^2 \theta_i = 0：点与投影方向正交，完全丢失
图1：累积能量曲线
前10个主成分的累积能量具体数据如下：
| 主成分数 p | 累积能量比例 |
|---|---|
| 1 | 96.63% |
| 2 | 99.37% |
| 3 | 99.74% |
| 4 | 99.84% |
| 5 | 99.88% |
p = 2, 3, 4 时的降维矩阵及误差结果如下：
| p | 降维矩阵 \boldsymbol{Z}_p 尺寸 | 累积能量 | 相对误差 |
|---|---|---|---|
| 2 | 728 \times 457 | 99.37% | 0.63% |
| 3 | 728 \times 457 | 99.74% | 0.26% |
| 4 | 728 \times 457 | 99.84% | 0.16% |
该图展示了PCA分析中前10个主成分的累积能量比例，可见，第一个主成分已经捕获了约96.6%的数据方差，绝大部分信息集中在前几个主成分中，这表明数据可以用少数几个主成分很好地近似原始数据。对于后续的NMF分解中我们选择 p = 3
图3：C_1 vs C_2 vs C_3 三维投影散点图
该图将数据点投影到前三个主成分张成的三维空间中。三维视图清晰地展示了点云沿着一个主方向延伸，右上角的密集区域对应背景像素，而沿对角线向外延伸的稀疏点对应不同强度的荧光信号。三维投影保留了99.7%的数据方差，几乎完整地表征了原始457维数据的结构特征。
图4：\cos^2 矩阵可视化（k=4）
通过计算 \cos^2 矩阵（投影坐标平方后按行归一化），可以分析各主成分对每个观测的贡献程度，辅助判断有效维度。
图中每个像素位置的灰度值表示该像素的光谱向量投影到前4个主成分子空间后保留的能量比例：白色（值接近1）表示该像素被很好地表示，黑色（值接近0）表示表示质量较差。从图中可以看出，大部分区域呈现浅灰色或白色，说明前4个主成分已经能够很好地表示绝大多数像素的光谱信息。
Q3：\boldsymbol{H} 和 \boldsymbol{S} 的尺寸
非负矩阵分解（NMF）将数据矩阵分解为：
其中：
- \boldsymbol{X} \in \mathbb{R}^{n_{pix} \times n_\lambda}：观测数据矩阵
- \boldsymbol{H} \in \mathbb{R}^{n_{pix} \times p}：丰度矩阵（混合系数矩阵）
- \boldsymbol{S} \in \mathbb{R}^{p \times n_\lambda}：源矩阵（端元光谱矩阵）
- p：源的数量
从高光谱解混合的角度理解，每个像素位置的光谱是一个混合信号。一个像素可能同时包含多种荧光标记物，观测到的光谱是各标记物光谱的加权叠加，权重反映各标记物在该像素的浓度。解混合的目标是从这些混合光谱中分离出源信号及其空间丰度图（混合系数）。
结合本实验ROI数据（n_{pix} = 728，n_\lambda = 457）：
当 p = 3 时：
当 p = 4 时：
- \boldsymbol{S} 的每一行是一个纯净的源光谱，长度为457个波长点
- \boldsymbol{H} 的每一行是对应像素中各源的丰度
- \boldsymbol{H} 的每一列重塑为 26 \times 28 图像即为该源的丰度图像
Q4：证明ALS是下降算法
ALS（交替最小二乘法）的目标是最小化目标函数：
- \boldsymbol{H} \geq 0，\boldsymbol{S} \geq 0
算法交替进行两步优化：
**步骤1：**固定 \boldsymbol{S}^{(k)}，优化 \boldsymbol{H}
由于这是关于 \boldsymbol{H} 的凸优化问题（二次目标函数 + 凸约束），存在全局最优解。因为 \boldsymbol{H}^{(k+1)} 是最优解，而 \boldsymbol{H}^{(k)} 是可行解，所以：
**步骤2：**固定 \boldsymbol{H}^{(k+1)}，优化 \boldsymbol{S}
同理，这也是凸优化问题，因此：
综合两步可得：
因此目标函数值在每次完整迭代后单调递减（或保持不变）。又因为 J \geq 0 有下界，根据单调有界定理，序列 \{J^{(k)}\} 必收敛。由此证明ALS是下降算法。
Q5：算法伪代码与停止条件
Q6：推导 \boldsymbol{S} 的乘法更新公式
已知 \boldsymbol{H} 的乘法更新规则为：
其中 \odot 表示逐元素乘法，除法也是逐元素运算。
利用Frobenius范数对转置不变的性质：
对目标函数取转置：
令：
- \tilde{\boldsymbol{X}} = \boldsymbol{X}^T
- \tilde{\boldsymbol{H}} = \boldsymbol{S}^T
- \tilde{\boldsymbol{S}} = \boldsymbol{H}^T
则目标函数变为：
则问题变为：
这与原来关于 \boldsymbol{H} 的优化问题形式完全相同，因此对 \tilde{\boldsymbol{H}} 应用乘法更新规则：
代入变换关系：
取转置得到 \boldsymbol{S} 的乘法更新公式：
其中 \odot 表示逐元素乘法，除法也是逐元素的。
Q7：实现乘法更新规则的NMF算法
本节实现基于乘法更新规则的NMF算法，将数据矩阵分解为 \boldsymbol{X} = \boldsymbol{H}\boldsymbol{S} 的形式，其中 \boldsymbol{H} 和 \boldsymbol{S} 的系数均为非负（避免除零错误）。输入数据为预处理后的ROI数据 \boldsymbol{X}，尺寸为 728 \times 457（n_{pix} \times n_\lambda）。源数量设为 p = 3 ，停止条件参数设为相对误差变化阈值 \varepsilon = 0.0001，最大迭代次数 K_{max} = 1000。最后设置随机种子 rng(42) 保证结果可重复。
优化目标
最小化重建误差：
交替最小化策略
由于同时优化 \boldsymbol{H} 和 \boldsymbol{S} 很困难，采用交替最小化的策略：
-
固定 \boldsymbol{S}，关于 \boldsymbol{H} 最小化
\min_{\boldsymbol{H} \geq 0} \|\boldsymbol{X} - \boldsymbol{H}\boldsymbol{S}^{(k)}\|_F^2 \quad \Rightarrow \quad \boldsymbol{H}^{(k+1)} -
固定 \boldsymbol{H}，关于 \boldsymbol{S} 最小化
\min_{\boldsymbol{S} \geq 0} \|\boldsymbol{X} - \boldsymbol{H}^{(k+1)}\boldsymbol{S}\|_F^2 \quad \Rightarrow \quad \boldsymbol{S}^{(k+1)} -
交替进行，直到误差收敛
算法伪代码
实验结果
图5：NMF乘法更新规则收敛曲线
该图展示了乘法更新规则NMF算法的收敛过程，横轴为迭代次数 k，纵轴为重建误差 \|\boldsymbol{X} - \boldsymbol{H}\boldsymbol{S}\|_F^2。从曲线可以看出，误差在前几次迭代中急剧下降，随后误差趋于收敛值。具体来说，算法在第415次迭代时满足收敛条件，最终重建误差稳定在 8.04 \times 10^8。曲线的单调下降特性验证了Q4结论：乘法更新规则是一个下降算法，目标函数在每次迭代后单调递减。
图6：NMF乘法更新规则结果（源光谱与丰度图像）
该图展示了乘法更新规则分解得到的3个源光谱（左列）及其对应的丰度图像（右列）。源光谱再不同位置有不同的峰值特征，丰度图像显示了每个源在空间上的分布，颜色越红（值越大）表示该源在该像素位置的含量越高，从中可以看出三种成分在ROI区域内呈现不同的空间分布模式
图7：原始光谱与重建光谱对比
然后我们选取了5个代表性像素来对比原始光谱（蓝色实线）与NMF重建光谱（红色虚线)。可以看出，重建光谱与原始光谱曲线完美重合，且原始光谱中的高频噪声在重建光谱中被平滑掉，这是因为NMF将数据分解为3个源的线性组合，起到了去噪效果。不同像素具有不同的光谱形态和强度，但NMF都能够准确重建，说明选择 p=3 个源足以描述数据中的主要光谱变化。相对重建误差仅为0.1249%，验证了分解的准确性。
Q8：实现投影梯度ALS算法及性能比较
NMF问题关于 (\boldsymbol{H}, \boldsymbol{S}) 联合优化是非凸的，但当固定其中一个变量时，关于另一个变量的子问题是凸的。交替最小二乘（ALS）算法利用这一性质，交替进行两步优化：固定 \boldsymbol{S} 优化 \boldsymbol{H}，固定 \boldsymbol{H} 优化 \boldsymbol{S}。每一步都是凸问题，可以高效求解。
已知 \boldsymbol{H} 的更新公式：
其中算子 [z]_+ = \max(z, 0) 表示投影到非负象限。
标量投影
对于标量 z，投影到非负实数的定义为：
矩阵投影
对于矩阵 \boldsymbol{Z}，投影操作是逐元素进行的：
即对矩阵的每个元素分别应用 [\cdot]_+ 操作。
\boldsymbol{H} 的更新公式（已知）
我们先固定 \boldsymbol{S}，然后目标是最小化 \|\boldsymbol{X} - \boldsymbol{H}\boldsymbol{S}\|_F^2。无约束最小二乘解为：
加上非负约束后，投影到非负象限：
\boldsymbol{S} 的更新公式（需要推导）
原问题：固定 \boldsymbol{H}，最小化
无约束最小二乘解：
对目标函数关于 \boldsymbol{S} 求导并令其为零：
整理得正规方程：
求解：
加上非负约束：
综上所述，两个变量的投影梯度ALS更新公式如下：
算法框架
实验结果
图8：NMF投影梯度ALS收敛曲线
该图展示了投影梯度ALS算法的收敛过程，横轴为迭代次数 k，纵轴为重建误差 \|\boldsymbol{X} - \boldsymbol{H}\boldsymbol{S}\|_F^2。与乘法更新规则相比，ALS算法呈现出截然不同的收敛特性：误差在前5次迭代内快速下降然后趋于收敛值。算法仅需61次迭代即满足收敛条件，最终重建误差为 7.84 \times 10^8，优于乘法更新规则的结果。
图9：NMF投影梯度ALS结果（源光谱与丰度图像）
ALS算法具体数值结果如下：
| 指标 | 数值 |
|---|---|
| 相对重建误差 | 0.1217% |
| H最小值 | 0.000000（\geq 0） |
| S最小值 | 0.000000（\geq 0） |
| H尺寸 | 728 \times 3 |
| S尺寸 | 3 \times 457 |
H和S的最小值均为0，满足非负约束条件。
两种算法性能比较
图10：两种算法收敛曲线对比（对数尺度与前50次迭代）
该图通过两个子图对比了乘法更新规则（蓝色实线）和投影梯度ALS（红色虚线）的收敛行为，可见ALS的收敛速度远快于乘法更新规则，ALS仅需61次迭代即可收敛，而乘法更新规则需要415次迭代。
图11：源光谱对比（归一化后）
该图将两种算法得到的源光谱归一化后进行对比，消除了尺度差异的影响。对于源1，乘法更新（蓝色）在450-700nm范围内呈现平坦的宽峰，而ALS（红色虚线）则呈现双峰结构。对于源2和3，两种方法都得到单峰光谱，但峰值位置略有差异。这些差异表明NMF问题存在多个局部最优解，不同算法会收敛到不同的解，但这些解都能以较小的误差重建原始数据。
图12：丰度图像对比（乘法更新 vs ALS）
对比两种算法得到的丰度分布，上排为乘法更新结果，下排为ALS结果。尽管两种方法得到的源光谱有所不同，但丰度图像展示的空间分布都很相似，都成功识别出了ROI区域内的三种荧光成分及其空间位置。
数据汇总如下：
| 指标 | 乘法更新规则 | 投影梯度ALS |
|---|---|---|
| 迭代次数 | 415 | 61 |
| 最终误差 | 8.04 \times 10^8 | 7.84 \times 10^8 |
| 相对重建误差 | 0.1249% | 0.1217% |
从实验结果可以看出，投影梯度ALS算法仅需61次迭代即可收敛，而乘法更新规则需要415次迭代。在最终重建质量方面，ALS的相对重建误差（0.1217%）略优于乘法更新规则（0.1249%），最终误差也更小，这是因为投影梯度ALS每步更新更接近最优解，因此收敛更快且最终结果更优。
由于整体问题非凸，两种算法都是局部最小化算法，最终结果依赖于初始解的选择。这也解释了为什么两种算法收敛到不同的源光谱：NMF问题存在多个局部最优解，不同算法可能收敛到不同的解。