解析N-FINDR：在纯像素假设下对几何线性解混方法的推广

《Analytica Chimica Acta》：Unraveling N-FINDR: A generalization of geometrical linear unmixing approaches under the pure pixel assumption

　　鲁斯塔姆·古利耶夫（Rustam Guliev）| 乌特·诺伊格鲍尔（Ute Neugebauer）| 克劳迪娅·贝莱特斯（Claudia Beleites）
莱布尼茨光子技术研究所（Leibniz Institute of Photonic Technology）是莱布尼茨感染光子研究中心（Leibniz Centre for Photonics in Infection Research, LPI）的成员，地址：德国耶拿市阿尔伯特-爱因斯坦街9号（Albert-Einstein-Str. 9），邮编07745。

**摘要**
**背景**：高光谱解混是高光谱图像分析中的重要任务，其目的是将混合像素分解为纯净的光谱特征。在现有的解混方法中，N-FINDR算法是最广泛使用的方法之一。然而，随着时间的推移，文献中出现了多种实现方式，这些实现方式在参数选择、算法细节和操作假设上往往存在差异。这种缺乏标准化导致了研究之间的混淆和困难。因此，该领域缺乏一个统一且明确的N-FINDR算法表述，无法清晰界定其原理、与相关方法的关系，并确保研究的可重复性。

**结果**：我们提出了一种通用且标准化的N-FINDR算法表述，明确定义了其几何原理、参数选择和操作步骤。以往描述中常用的模糊缩写被替换为具体的参数值，从而实现了透明解释和一致的实施。所提出的框架阐明了N-FINDR与其他纯像素端元提取方法（特别是顶点成分分析（Vertex Component Analysis, VCA）之间的关系，突出了它们之间的共同假设和关键差异。全面的基准测试实验评估了算法变体、初始化策略和参数设置对提取性能的影响。基于这些分析，我们确定了最佳实践配置，这些配置能够在不同数据集上产生稳健且可重复的结果。为了便于采用，我们发布了开源的Python和R语言包来实现推荐的算法变体，使专家和非专家用户都能可靠地使用N-FINDR，而无需进行详细的算法调整。

**意义与创新**：本研究通过提供标准化、透明且可重复的N-FINDR算法表述，解决了长期以来围绕该算法的歧义。通过明确其与其他纯像素方法的关系，并通过系统性的基准测试建立最佳实践，提高了高光谱解混研究的可靠性。开源实现的可用性进一步降低了进入门槛，使得N-FINDR在更广泛的高光谱成像社区中得到更一致的应用。

**1. 综述**
高光谱解混是高光谱数据分析的重要组成部分，它将混合的高光谱数据分解为其组成的纯净光谱特征（即端元）。这一过程包括识别（端元提取）和量化（丰度反演）高光谱数据集中的各个光谱成分。这种分解提供了关于场景中存在的材料或物质、它们的空间分布及其相应丰度的宝贵信息。高光谱解混在许多应用中发挥着关键作用，例如遥感和微生物学，能够识别和表征样本中的不同物种或生化成分[1]、[2]、[3]、[4]。

对于许多类型的线性光谱学（如拉曼强度、红外（IR）吸收和荧光发射），p种化学成分或物种的混合物光谱通常可以近似为这些纯成分光谱的线性组合，权重由相应成分的浓度（丰度）决定。
丰度预期为非负值，并且总和为一（作为成分的摩尔分数）。从几何学角度来看，这一约束意味着所有可能的p种成分混合物都位于一个p?1单纯形（simplex）内部。这一事实促成了基于几何的线性解混方法的发展。自20世纪90年代以来，这些方法根据是否使用所谓的“纯像素假设”分为两个主要分支：
- 在纯像素假设下，假设所有成分的纯光谱都存在于数据集中。其中，迈克尔·温特（Michael Winter）提出的N-FINDR算法[5]被认为是开创性的。其基本思想是提取产生最大单纯形的p个光谱（单纯形体积最大化）：当满足纯像素假设时，纯光谱对应于p?1单纯形的顶点。
- 相反，另一个分支计算的数据中不存在的端元。例如，克雷格（Craig）提出了搜索包含所有光谱的最小单纯形的策略（单纯形体积最小化）。

尽管N-FINDR算法的概念很简单，但随后出现了三个明显的问题：即使是“温特的N-FINDR”这一称呼也存在歧义，因为他的原始论文[5]描述的迭代方法与后续出版物[7]中的伪代码不同。后续对搜索策略的改进导致了多种算法变体的出现，这些变体有许多不同的名称和缩写，如SM N-FINDR、SQ N-FINDR、SC N-FINDR、IN-FINDR、LDU-N-FINDR、LDU-S-N-FINDR、S-N-FINDR、MN-FINDR、MN-FINDR2、RN-FINDR等。不幸的是，有时相同的算法会被重新发明并以不同的名称发表；在其他情况下，相同的术语在算法名称中具有不同的含义[8]、[9]、[10]、[11]。这种歧义和实现不兼容性问题之前已被指出[12]。尽管有人尝试对不同的实现进行分类和比较[8]、[12]、[13]，但问题仍未解决。
关于这些算法及其相互关系的困惑和误解导致了现有实现的重新发明（例如[14]、[15]）以及一些不必要的比较和优化[16]。相反，诸如向量化和平滑丰度反演等重要方面却未得到足够的关注。例如，向量化可以利用针对标准（向量和矩阵）线性代数操作高度优化的库和硬件，从而在使用高级源代码时提供相似或更好的计算时间，在许多情况下甚至不需要进一步并行化。尽管如此，文献中很少提及向量化；而并行化在许多研究中得到了讨论[16]、[17]、[18]、[19]。
所谓顺序算法的出现（其中最著名的是顶点成分分析（Vertex Component Analysis, VCA）[20]和自动目标生成程序（Automatic Target Generation Procedure, ATGP）[21]）进一步增加了这一主题的复杂性。几何解释表明，它们仅在无噪声数据的情况下才与N-FINDR相当（见第4.1.1节）。它们被建议用于N-FINDR的初始化[22]、[23]，然而应用论文似乎将它们用于含噪声的拉曼数据，仿佛它们是潜在的竞争者[1]、[24]。

**2. 算法介绍**
2.1. 术语和缩写
为了提高可读性，我们提供了与本工作相关的术语和缩写：
- **纯像素算法**：指依赖于纯像素假设并基于单纯形几何的搜索策略的算法。它还包括纯变量算法的子组，例如SIMPLISMA（简单易用的交互式自建模混合分析[25]），这些算法通过相同的几何原理识别由单一成分主导的光谱通道。虽然这是一个不同的问题，但通过数据矩阵转置在数学上是相关的（前提是假设成立，例如[26]）。因此，我们将这两种方法视为同一类方法，并统称为纯像素算法。对于那些每次迭代选择一个端元直到达到所需数量的纯像素算法（例如VCA），我们使用“顺序算法”这一术语。
- 其他常用术语包括：
- **视角**：光谱（数据作为光谱集合）、几何（数据作为多维空间中的点）和空间（数据作为像素的空间网格）；
- **性质**：真实量（可能不直接出现在数据中）及其数据衍生的估计值。
- 术语“纯成分光谱”、“顶点”（几何单纯形的角）和“纯像素”（包含单一成分的空间像素）用于指代真实量；对于相应的数据衍生估计值，我们使用“端元光谱”、“端点”和“端元像素”。除非另有说明，顶点和端点定义在搜索空间中，通常是降维后的空间。这种术语可能与相关文献中的惯例不同，但选择它是为了避免可能的歧义。

2.2. 符号说明
从现在开始，将使用以下符号：
- **n**：高光谱图像中的光谱或像素数量。
- **p**：模型中的端元数量，即混合物中的纯成分数量。
- **m**：每个光谱的光谱通道数量（例如波数）。
- **X∈R^n×m**：混合物光谱的数据集（可以预处理）。
- **C∈R≥0^n×p**：成分的丰度（浓度）。
- **S∈R^p×m**：混合物的端元光谱。
- **Z∈R^(n×(p?1)**：降维到p?1的数据集X。
- **E∈R^p×(p?1)**：端点点的矩阵，是Z的子集。
- **E?∈R^p×p**：单纯形体积矩阵，通过添加一列1来扩展矩阵E。
- **det???**：在单纯形体积的上下文中，用于表示矩阵行列式的绝对值。

2.3. 假设
算法基于以下假设[5]：
1. 高光谱数据中的每个光谱都是某些纯成分光谱的线性组合。虽然这对许多类型的光谱数据是一个有效的近似，但这一假设可能因多种原因而被违反。线性光谱学将混合物光谱近似为成分光谱的加性组合；对于许多类型的光谱学（如拉曼或荧光发射），强度或强度的适当变换（如吸收）被近似为与测量体积中相应原子或分子的数量成正比，即假设与摩尔浓度线性相关。质量浓度可能已经表现出显著的非线性。对于具有不同光-物质相互作用的样本（如吸收光谱中的散射损失），也可能出现偏差。光学设置和传感器也可能引入非线性。预处理可以应用于转换原始光谱数据，以使线性近似比原始数据更准确。这可能是使用吸收度A=?logI/I0而不是强度I的变换，或者使用扩展乘法散射校正（Extended Multiplicative Scatter Correction, EMSC）等技术来去除漫反射光谱中的附加和乘法散射效应[31]。
2.4. 线性组合的系数预期为非负值。
3. 对于每个成分，我们假设数据X中至少存在一个纯成分光谱（纯像素假设）。N-FINDR的目标是在众多混合物光谱中定位这些纯成分光谱。这一假设可以在一定程度上放宽[5]。然而，化学相互作用是一个问题。从化学角度来看，可能会产生新的物种（例如表面物种），但仍可以假设线性。然而，除非这种相互作用物种的纯光谱直接存在于数据中，否则纯像素假设就被违反了。这样的数据可能看起来是非线性的，特别是如果缺失的物种仅在特定浓度范围内形成。一般来说，这样的数据不应通过N-FINDR进行分析。
4. 数据中的噪声水平可以忽略不计。算法通常涉及降维步骤。最常见的方法是PCA，它可以减少不相关的噪声。然而，在应用算法之前，其他类型的噪声（如拉曼光谱中的宇宙射线尖峰或化学噪声/杂质）需要单独处理。实际上，算法仍然会收敛，但可能不会得到满意的解决方案。5. 纯组分的数量（p）是已知的。在实践中，组分的准确数量通常是未知的。这是化学计量学中的一个普遍问题，称为化学秩估计。对于解混问题，已经提出了多种算法[32]、[33]。另一种常见的方法是多次运行算法并改变p的值，然后通过光谱解释等方法选择其中一个结果。需要注意的是，通常纯组分的预期数量远小于数据中的光谱数量（n）和原始数据的维度（m）。例如，本研究中的基准数据集的n最多为10^5，m约为200-600，p≤20。2.4. 双线性模型数学上满足上述假设意味着每个光谱x可以表示为以下形式：x=c1s1+c2s2+?+cpsp+?，其中ci≥0，i=1,2,…,p，并且∑i=1 ci=1或者用矩阵形式表示为：(1) X(n×m)=C(n×p)S(p×m)+E(n×m)这种关系是双线性模型用于建模光谱数据的重要性的基础。矩阵C将每种混合物表示为浓度空间中的一个点。换句话说，C中的每一行代表混合物图中的一个点。这个浓度空间是有约束的，因为组分的摩尔分数必须加和为一。分解本身并不是一个唯一的概念。例如，主成分分析（PCA）是一种广泛使用的算法，它采用了类似的分解方法，尽管没有特定的约束。然而，这些算法背后的动机不同。线性解混并不声称提供通过PCA等方法无法获得的全新信息。相反，它应该被视为针对特定用途而定制的不同表示方法。2.5. 单形体的概念在几何学中，单形体是指三角形或四面体概念的推广到任意维度。对于任何p≥1，p-1单形体是p个仿射独立点的凸包。这包括所有形式为∑i=1 ci ei的点，其中ci是非负的并且总和为一。这些点ei称为单形体的顶点。ci是重心坐标，即相对于单形体顶点的坐标。例如，2-单形体对应于一个三角形，即三个顶点、三条边以及所有内部的点。p-1维度称为p-1单形体的内在维度。请注意，数据的维度可以更高。一个简单的例子是在3d空间中的线段（1-单形体）或平面三角形（2-单形体）。当假设得到满足时，所有可能的p组分混合物都位于p纯组分的p维浓度空间内的p-1单形体中。方程（1）的双线性形式意味着混合物光谱X也位于光谱空间中的p-1单形体内，即由光谱数据集的m个光谱通道（例如，波数）所跨越的m维空间。纯组分光谱是相应p-1单形体在光谱空间中的顶点。在这项工作的背景下，另一个重要的概念是单形体体积。以下公式通常用于计算其内在空间中的p-1单形体体积：(2) V simplex = 1/(p?1)!det(E?) = 1/(p?1)!det[1e1,1…e1,p?1] [1e2,1…e2,p?1] ? ? ? ? [1ep,1…ep,p?1]这个公式以及其他计算单形体体积的方法在补充部分S1中有说明。2.6. N-FINDR启发式方法由于纯组分光谱对应于光谱空间中p-1单形体的顶点，因此寻找端元可以简化为识别由数据X生成的最大的单形体所跨越的p个数据点。简而言之，该算法可以描述如下：1. 将数据维度降低到内在维度p-1。这可以通过PCA来实现。假设4（可忽略的噪声）保证了单形体位于前p-1个维度内，因此可以方便地截断数据。经过这种转换后，算法的数学表述更简单，计算要求也更低。然而，PCA并不是唯一可行的方法。最小噪声分数（MNF）[14]和特定波段域[36]也被使用过。还有一种不进行维度缩减的N-FINDR变体[14]，但它并不改善分解效果[14]（另见第4.1.5节）。2. 随机选择p个点作为端元点的初始候选者。3. 试错地用数据矩阵中的其他点替换暂定的端元点。4. 如果单形体体积增加，则保留替换；否则丢弃更改。为此，只需要计算行列式，因为1/(p?1)!因子对于给定的维度是常数，因此对最大化无关紧要。5. 重复步骤3和4，直到没有增加单形体体积的替换，即单形体不再增大体积（收敛到局部最大值）。图1展示了使用N-FINDR算法进行光谱解混的一般工作流程。我们出于将在第4.2节讨论的原因，从传统描述中排除了维度缩减。该图还表明N-FINDR只是数据分析工作流程的一部分。重要的是，N-FINDR并不直接提供丰度估计；相反，它识别出最有可能是纯组分光谱的像素。尽管这个概念很直接，但其实现存在显著的歧义。因此，出现了多种变体，并且有不同的缩写。下载：下载高分辨率图像（592KB）下载：下载全尺寸图像图1. 使用N-FINDR进行光谱解混的一般工作流程。初始的高光谱（HS）图像显示为一个3d数据立方体（a），然后展开为一个2D矩阵（b），其中每一行代表一个光谱。因此，其他光谱数据集也可以用N-FINDR处理。每个光谱可以可视化为m维空间中的一个点（在这种情况下选择了780 cm?1、1000 cm?1和1440 cm?1的波段进行可视化）（c）；图a-c是相同输入数据的不同表示方式。工作流程如下：（1）维度缩减，（2）在降维空间中提取端元（N-FINDR步骤，用黄色突出显示），（3）使用提取的端元进行丰度反演，以及（4）通过为每个端元分配颜色并根据计算出的丰度将它们混合来可视化伪彩色图像。（关于图中颜色参考的解释，请参阅本文的网页版本。）图1还清楚地表明，N-FINDR算法并不涵盖完整的解混工作流程。特别是，N-FINDR并不直接提供丰度估计。突出显示区域之外的步骤是许多端元提取算法的共同步骤。要全面理解这些步骤，需要进行专门的研究，这超出了当前工作的范围。2.7. 理想情况和顺序方法考虑一个所有假设都完美满足的理想情况是有帮助的。为了简化可视化，我们考虑一个有3个纯组分的场景。因此，数据在由m个光谱通道（例如，波数）跨越的m维空间中形成一个完美的三角形。所描述的N-FINDR启发式方法通过迭代增长单形体是有效的，并且明确地产生正确的结果单形体。尽管如此，它可能不是最有效的方法。几乎任何将单形体投影到直线上的情况都会得到一个1-单形体（即线段），其顶点对应于原始单形体的顶点。因此，可以通过将单形体投影到几个随机向量上并选择每个投影的极端点来找到顶点。图2A以图形方式展示了这种简单2D情况。本质上，这个概念构成了顺序方法（如VCA和ATGP）的基础。虽然具体方法在选择投影向量方面有所不同，但基本思想是相同的：顶点对应于投影的极端点。例如，VCA通过将数据投影到与先前找到的端元子空间正交的方向上来迭代识别新的端元点，直到达到所需的端元数量。图2B以图形方式展示了这种方法。ATGP算法在本质上非常相似，但它不需要预先进行维度缩减，并且它选择幅度最大的向量而不是随机向量投影。更详细的描述可以在[22]、[23]、[37]中找到，尽管这些描述并不总是完全一致。请注意，这些方法要求数据的维度至少为p，即比内在单形体维度p-1多一个。图2B还表明，为了找到第三个端元，需要多一个维度，因为VCA不考虑对应于最小值的顶点[20]。下载：下载高分辨率图像（345KB）下载：下载全尺寸图像图2. 在N-FINDR假设完全满足的理想化情况下的模拟PCA投影。每个点代表维度缩减（例如，PCA）后的光谱。红色顶点和连接它们的三角形表示一个单形体。单形体内的黑点对应于混合光谱。黑色箭头代表投影向量。红色箭头是端元向量。这些箭头上的灰色点是数据点的投影。其中，黑色投影点表示投影后选为端元的点。比较了两种投影策略：（A）随机投影：将点投影到随机选择的方向上，并选择两个极端投影（最小值和最大值）作为候选端元。（B）VCA：新的投影向量与先前选定的端元正交，并且每个方向只选择最大投影。（关于图中颜色参考的解释，请参阅本文的网页版本。）几何解释表明，投影向量不需要是正交的。这个概念为PPI算法奠定了基础，该算法识别在多个随机投影下最极端的点[30]。此外，可以通过选择与坐标系匹配的向量来简化投影。换句话说，可以跳过专门的投影，只需选择具有极端坐标的点，无论是在光谱空间中还是例如在PCA进行维度缩减时选择最小值和最大值。需要注意的是，这些方法仅在理想的无噪声情况下有效，即使是轻微的噪声也会显著影响结果。虽然PCA或类似的投影到内在空间的方法可以帮助减少噪声，但这会抵消顺序算法的计算速度优势。我们将在第4.3节讨论与N-FINDR算法相关的更现实的情况。3. 材料与方法3.1. 考虑的算法范围本工作主要关注纯像素算法，特别是各种N-FINDR变体及相关方法，这些方法虽然不使用N-FINDR名称，但基于相同的原则（例如[29]）。顺序方法（VCA、ATGP、OPA等）也与N-FINDR算法一起考虑，即我们不是比较它们之间的差异，而是讨论顺序方法的一般属性。为了说明，我们使用VCA和ATGP作为代表性示例，涵盖了有无维度缩减的情况，以及随机输出和确定性输出的情况。虽然一些算法包括预处理（例如VCA中的维度缩减），但这里的重点在于端元搜索策略。不再进一步讨论预处理，因为它依赖于数据并且通常可以在方法之间转移；例如，VCA中使用的维度缩减也可以在N-FINDR之前应用。在纯像素算法中，像素纯度指数（PPI）[30]没有详细研究，因为它在概念上有所不同，且很少与其他纯像素方法混淆。它根据随机投影中的极端性对像素进行排名，而不是直接识别端元。它有时被用作N-FINDR的预处理步骤；更多细节可以在文献[8]、[22]、[30]中找到。SIMPLISMA[25]也没有在本工作中讨论。最初作为纯变量算法提出，并且也适用于纯组分分析[26]、[27]，它遵循不同的原则，结合了基于行列式的选择和方差加权。可以将SIMPLISMA视为带有先验向量归一化的顺序算法的加权版本，搜索既高度正交又具有高方差的端元。以类似的加权方式扩展N-FINDR是可行的，但这将引出关于解混过程中加权和标准化的更广泛问题，这超出了本工作的范围。其他解混方法，如非线性解混[38]、统计或贝叶斯技术[39]，或者不基于纯像素假设的线性解混方法（例如，迭代约束端元（ICE）[40]、多变量曲线分辨率（MCR-ALS）[41]、神经网络和深度学习[42]等，也超出了本文的范围。尽管我们努力涵盖了广泛的实现和方法，但仍可能存在一些遗漏。尽管如此，所提出的框架能够直接整合额外的方法，并基于特定参数差异进行比较。

3.2. 数据
我们使用两个数据集作为实验数据，用于性能评估和展示特定的算法方面。这些数据集涵盖了光谱解混的两个重要光谱模式和应用领域：生物细胞的拉曼显微光谱成像以及用于地质调查的VIS/NIR反射率遥感。

3.2.1. AVIRIS Cuprite：遥感红外成像
AVIRIS（机载可见光红外成像光谱仪）Cuprite数据集是一个公开可用的遥感高光谱数据集。有关其应用背景、详细描述和全面的光谱及地质分析，请参见[44]。
简而言之，AVIRIS是NASA喷气推进实验室的扫帚式扫描仪，它可以在224个连续的光谱通道中收集光谱，覆盖0.38 μm–2.5 μm的波长范围。这里使用的数据集是在1997年6月由一架飞行在21公里高度的飞机在美国内华达州的Cuprite矿区收集的，生成了一个高光谱VIS/NIR图像，其中每个像素代表了地面上18 m × 18 m区域的光谱测量结果。总共扫描了10.5 km × 37.5 km（614像素 × 2206像素）的区域，分为五个场景（子图像），场景1至4每个覆盖10.5 km × 8.7 km（614像素 × 512像素）的区域，场景5覆盖10.5 km × 2.7 km（614像素 × 158像素）的区域。
这个数据集对于研究包括N-FINDR[5]、[8]、[16]、[45]、[46]、[47]、[48]在内的光谱解混算法非常有价值，因此我们将其包括进来以便与现有文献进行比较。我们假设端元的参考数量为22个，如[22]中所建议的。
尽管有更新的Cuprite区域扫描数据可用，但我们使用的是1997年扫描的场景4的反射率数据，该数据从官方AVIRIS数据分发网页[49]中获取，文件名为f970619t01p02_r02_sc04.a.rfl。这个场景具有地质学意义。我们还采用了与N-FINDR文献中相同的预处理方法：在分析之前移除了光谱带1至6、105至115、150至170和221至224，保留了182个光谱带。由于水的吸收，被移除的光谱带的信噪比较低。之前的N-FINDR算法文献使用未经转换的反射率R[5]、[8]、[16]、[45]、[46]、[47]、[48]。因此，我们坚持这一决定，无论从光谱学角度来看是否更倾向于使用?log(R)的双线性近似。

3.2.2. 巨噬细胞：细胞的拉曼图像扫描
此外，我们还包括了一个数据集，该数据集使用拉曼图像光谱学来研究生物学中的光谱解混，以表征巨噬细胞——一种免疫细胞。该数据集包含65个预处理的巨噬细胞拉曼矩形图像扫描，总计约180,000个光谱。这些巨噬细胞来自三名健康人类的血液样本，并经过处理以获得活化（M1）、交替活化（M2）和非活化（M0）表型。每个拉曼图像扫描的空间步长为0.5 μm，每个光谱的积分时间为1秒。扫描区域不同，以捕捉整个细胞，边长范围从10 μm到100 μm。数据预处理包括手动去除非细胞背景、排除异常光谱、基线校正，以及将光谱区域限制在350 cm-1到1800 cm-1和2650 cm-1到3100 cm-1之间。有关更多详细信息和深入分析，请参见[2]。我们使用端元的参考数量15来匹配原始工作的分析。

3.2.3. 演示数据 - 单个巨噬细胞的拉曼图像
为了演示目的，我们还使用了一个单个活化（M1）巨噬细胞的拉曼图像，该图像是巨噬细胞数据集的一部分。该图像用三个组分进行建模：细胞核、细胞质和水背景。图1中一般光谱解混插图中的PC分数和最终假彩色图像展示了这个数据集。
请注意，数据并没有形成一个完美的单纯形，但这种近似足以覆盖大部分数据，并且提供了高光谱图像的满意可视化。

3.2.4. 模拟拉曼数据集
为了补充实验数据集，我们还使用了一个模拟的拉曼光谱数据集，该数据集旨在模拟在受控条件下的真实光谱混合。这个数据集是基于纯物质的实验测量参考光谱生成的，包括牛血清白蛋白（BSA）、脱氧核糖核酸（DNA）、果糖、D-(+)-葡萄糖、海藻糖、聚苯乙烯和硅（Si）。模拟包含了五种目标组分（BSA、DNA、果糖、葡萄糖和海藻糖）的凸混合体，其中海藻糖的丰度限制在0.3以下，以代表一个不完整的组分。此外，还包含了两个单像素污染物（聚苯乙烯和Si，强度较低），以模拟罕见或人为产生的光谱。添加了真实的噪声源，包括多项式基线、散粒噪声、暗电流和电子噪声，以反映典型的拉曼测量条件。预处理包括基线校正和光谱区域裁剪至400 cm-1到1800 cm-1和2700 cm-1到3100 cm-1。有关模拟模型、参考光谱和预处理的更多详细信息，请参见相应的补充文件。

3.3. 基准测试
默认情况下，基准测试和其他计算是使用R [50]版本4.5.1（2025-06-13）（平台：x86_64-pc-linux-gnu，在Ubuntu 24.04.2 LTS下运行）在Lenovo Thinkpad P14s G5笔记本电脑上进行的，该笔记本电脑配备了Intel(R) Core(TM) Ultra 7 155H处理器（22个逻辑核心，详见规格），64 GB DDR5内存和1TB SSD。基准测试没有使用图形处理单元（GPU）。最终实现也在相同的硬件上使用Python 3.12.0（64位）进行了测试（第4.3.5节）。补充部分S2提供了关于R和Python的基准测试环境的更多详细信息，包括使用的软件包版本、BLAS（基本线性代数子程序）和LAPACK（线性代数包）库，以及硬件规格。所使用的OpenBLAS线性代数库支持多线程，这一点通过在基准测试期间监控CPU负载得到了确认。
基准测试和执行时间的分析是使用bench [51]包进行的：禁用了检查结果（check argument）参数，最小和最大迭代次数（min_iterations和max_iterations）分别设置为20和30，时间单位设置为毫秒，其余参数保持默认值。在总结之前，移除了至少包含一次垃圾收集（自动内存清理）的迭代。因此，当需要计算结果（如指数、体积、替换次数等）时，进行了简单的迭代（使用R中的*apply iterators）。在这些情况下，时间使用system.time进行测量，并在调用垃圾收集器之前立即进行。
随机数生成器（RNG）为每个表达式使用相同的种子进行初始化，以确保结果的可比性，例如，选择了相同的初始端元。
默认情况下，用于基准测试的N-FINDR实现对应于r-hyperspec/unmixR GitHub仓库的版本2.6。特别是，这个版本尽可能地使用了向量化计算。版本1.0（未向量化）被用来通过比较来展示向量化的影响。
除非另有说明，否则在运行N-FINDR之前会对数据应用PCA降维，数据先进行中心化（减去均值）且不进行缩放。尽管这不影响不同策略的比较，但在考虑绝对值或特定优化策略的实际相关性时必须考虑到这一点。PCA是对整个数据集应用一次，然后使用相应的分数矩阵子集作为N-FINDR的输入。PCA是使用svd R包[52]提供的截断SVD分解的PROPACK实现计算的。
最后，当数据集的大小和/或端元数量发生变化时，相应的数据子集可能不再形成一个单纯形。尽管如此，我们仍然认为提供的基准测试结果适合比较和讨论算法的某些方面的计算性能，从而找到最佳的默认实现。

3.4. 代码
所提出的泛化的开源实现已经在R和Python中开发。R实现是从GitHub上的现有光谱解混包unmixR开始重构的，该仓库为r-hyperspec/unmixR。重构前后的版本分别标记为“v1.0”和“v2.6”。算法的Python版本提供了一个新的包pyspc-unmix，可以在r-hyperspec/pyspc-unmix仓库上找到。该实现旨在无缝集成到现有的工作流程中，类似于其他scikit-learn转换（如PCA）的方法。
值得注意的是，这两个包提供的参数集并不相同。unmixR包由于其用于基准测试而包含了更广泛的实现。相比之下，Python版本仅提供了一组经过优化的参数，这些参数证明了最佳性能。
以下R包被用于创建图表：ggplot2 [53]、gridExtra [54]、LearnGeom [55]、scales [56]、plotly [57]、cowplot [58]和ggbeeswarm [59]。其他用于数据操作的R包包括dplyr [60]、tidyr [61]、hyperSpec [62]、pracma [63]、bench [51]和svd [52]。

4. 结果和讨论
4.1. 算法实现的系统检查
在本节中，我们系统地检查了算法的不同方面，以了解实现的多样性、它们之间的关系、有效性，并解决文献中未充分讨论或误解的方面。

4.1.1. 顺序方法及其与N-FINDR的关系
顺序方法（一次选择一个端元的方法）通常被视为N-FINDR的替代方案，因此进行了比较[1]、[24]、[64]、[65]、[66]、[67]。正如引言（第2.7节）中所示，这些方法与理想化的、无噪声的场景密切相关。实际上，这些算法并不声称能产生最优输出（例如[20]）。相比之下，N-FINDR对数据中的噪声更为鲁棒。这主要是由于其迭代细化临时解决方案的能力[12]。因此，可以将顺序算法的解决方案与N-FINDR进行比较，N-FINDR在每次细化一个端元后有一个提前停止的条件。有趣的是，一些顺序方法的作者还提出引入一个外循环来改进解决方案[29]，实际上相当于具有自定义初始化的N-FINDR。
VCA在这种情况下特别相关，因为它在应用工作中经常与N-FINDR进行比较[1]、[24]。然而，直接比较这两种方法可能会产生误导，因为VCA（在其原始形式中）包括两个不同的步骤：基于信噪比（SNR）的条件降维，然后是端元点的顺序提取。因此，当VCA产生的结果优于N-FINDR时，可能不清楚改进是由于降维步骤提供的初始化还是由于顺序提取策略本身。
从理论角度来看，顺序方法在最大化最终单纯形体积方面不期望优于N-FINDR。这一观察结果得到了第4.3.4节中基准测试结果的支持，其中N-FINDR获得的最终单纯形体积可以是VCA识别的初始单纯形的两倍或更多。因此，在实际应用中比较VCA和N-FINDR时应仔细区分初始化的效果和搜索策略本身的效果。
N-FINDR的搜索策略在较大单纯形体积对应于更具代表性的端元的情况下预计表现更好。这一假设在实践中通常成立，尽管根据数据的结构和噪声特性可能会出现例外。另一方面，顺序方法可以被视为为N-FINDR算法提供良好的初始化，这将在下一节中详细讨论。

4.1.2. 初始化
在其原始形式中，算法从选择p个随机数据点作为初始端元开始。许多作者发现这一方面存在问题，并努力寻找更好的初始化方法。在这方面提出了许多修改建议，例如选择相关性最低的像素[68]、过滤数据点（例如，在[8]中使用阈值进行数据球化处理）、使用序列算法[22]、[23]、采用受k-means启发的最大-最小距离方法[22]、[69]以及其他方法[8]、[22]、[70]、[71]、[72]。首先，初始化问题似乎被夸大了。实际上，输出的可重复性并不是一个严重的问题。已经有成熟的方法来处理随机初始化问题，例如设置随机生成器并/或运行多次随机初始化，然后通过单纯形体积选择最佳输出。其次，大多数建议的初始化方法都需要额外的预计算，这可能会完全抵消从更好初始点开始的优势，或者在某些情况下甚至使情况变得更糟[8]。一些研究提出使用序列方法的输出作为N-FINDR的初始化。例如，在[22]中，作者比较了不同的初始化方法，并发现ATGP的性能优于其他初始化策略。后来，[23]也使用ATGP进行初始化，并表明它导致更少的替换。确实，序列方法，特别是ATGP，可以以较低的计算成本有效减少替换次数，因为基本上只使用了投影。然而，需要澄清的是，一旦计算被向量化，序列投影在计算上就与使用最内层循环中的点类似（特别是使用高度体积变化估计器时）。实际上，与高度体积变化的相似性如此之强，以至于它们很容易被混淆（例如，在[73]中的VCA示例）。因此，运行序列方法在计算上与使用最内层循环中的点迭代策略进行单次N-FINDR操作是相似的。

使用单个细胞拉曼图像的演示数据，图3展示了不同的初始化方法，包括第2.7节讨论的投影方法。该图还显示了最终的输出结果。
下载：下载高分辨率图像（615KB）
下载：下载全尺寸图像
图3. 应用于单个巨噬细胞细胞的演示数据的不同初始化方法：VCA - 使用VCA的输出进行初始化，ATGP - 使用ATGP的输出进行初始化，极端坐标和极端坐标（随机）- 顺序或随机选择具有极端坐标的点，随机投影 - 通过生成随机投影向量并选择极端值来选择点，N-FINDR（1）- 使用N-FINDR最外层循环的一次迭代作为初始点。红色圆圈表示初始化方法的输出，而蓝色十字表示最终的N-FINDR输出。（关于此图例中颜色的解释，请参阅本文的网页版本。）

4.1.3. 迭代
初始化之后，临时顶点会不断被替换，直到体积不再增加。然而，有多种遍历数据点的方法。原则上，以下策略是可行的：
•暴力法。即，检查n个点中的所有可能组合p点，并选择体积最大的那个。这种方法在[8]中也被称为SM-N-FINDR（或SiMultaneous N-FINDR）。在实践中这是不可行的，因为它需要(np)=n!p!(n?p)!次单纯形体积的计算。例如，在仅10,000个光谱中搜索5个端元，就已经需要大约8.3×10^17次体积计算。如果使用矩阵行列式进行体积计算，并将其并行化到所有22个核心上，在本研究中使用的计算机上这将需要超过4000年的时间。实际上，N-FINDR应用于包含数百万光谱的数据集，因此使用启发式方法来优化单纯形体积。
•点和端元。另一种方法是使用三个嵌套循环：最外层循环检查收敛性，即在一个完整的最外层迭代中没有任何端元发生变化。接下来的两个循环（即对于每个端元和每个点）分别遍历端元和点。这些循环的顺序可以互换。上述伪代码代表了points-in-innermost-loop（简称points-innermost）策略，而endmembers-innermost策略将使用相反的循环顺序（即对于每个端元在每个点内）。这两种策略成为了N-FINDR不同缩写的来源。在之前的出版物中，points-innermost和endmembers-innermost迭代有不同的称呼（按相应顺序）：SQ-N-FINDR/SC-N-FINDR[8]、MN-FINDR2/MN-FINDR1[11]、LDU-S-N-FINDR/LDU-N-FINDR[10]、Implementation #2/Implementation #1[74]、使用新搜索策略的N-FINDR/原始N-FINDR[75]以及N-FINDR-BR（最佳替换）/N-FINDR[12]。虽然在最理想的无噪声数据情况下最外层循环可能看起来不重要，但在实践中它很重要。对于实验数据，它通常会导致最终单纯形的体积更大（因此，端元提取效果更好）。如果最外层循环的第一次迭代就已经达到了最大体积单纯形，这可以通过在其第二次迭代中没有进一步替换来检测到。最初的论文[5]就指出了这一点的重要性。尽管如此，后续的一些论文考虑了多次或仅一次最外层迭代的N-FINDR，例如SI-N-FINDR（单次迭代）/FI-N-FINDR（全迭代）[12]，或者ISQ/SQ N-FINDR和ISC/SC N-FINDR[8]。Dowler等人[12]详细研究了这一点，并证明使用最外层循环是最一致准确的方法。后来的工作大多展示了完全迭代直到收敛的算法。然而，仍然有一些最近的研究忽略了这一点[16]。
•同时优化点和端元。这种策略将两个最内层循环合并为一个，即考虑所有n×p种可能的单次替换（第j个端元由第i个点替换），并选择体积增加最好的那个。虽然并非所有已知的体积变化计算都适用于这种方法，但有些是可以的[11]、[76]。然而，这种循环策略在文献中并未提及。因此，我们将其作为算法泛化的一部分提供。

补充图S2提供了三种迭代策略的视觉演示：endmembers-innermost、points-innermost以及同时进行这两种迭代。

4.1.4. 体积变化估计器
由于只关心相对体积（即是否增加），可以探索计算体积变化的替代方法。
一类改进源自单纯形体积公式Eq. (2)本身。用另一个数据点替换一个临时端元点意味着只需改变单纯形体积矩阵E?的一行，其余部分保持不变。因此，可以重用一些之前计算过的矩阵行列式部分。基于矩阵/行列式分解的几种改进已被提出：基于下对角线上三角（LDU）矩阵分解的LDU-N-FINDR[10]，以及基于矩阵逆分解、共因子矩阵行列式展开和Woodbury矩阵恒等式的方法I - III[8]、[76]。
然而，这些矩阵分解方法必须谨慎应用。例如，基于LDU的方法依赖于以下公式：det(E?)=det(I0cTA?1IA00d?cTA?1bIA?1b0I)，但没有保证块子矩阵A是可逆的。
其他方法基于算法的几何解释。Ji等人[11]利用了克拉美法则，即E??1zi=det(E?i,1)det(E?)det(E?i,2)det(E?)…det(E?i,n)det(E?)T，其中zi是用于替换的新点，E?i,j是在用第i个数据点替换第j个临时端元点后的单纯形体积矩阵。克拉美法则的良好几何解释可以在[77]、[78]中找到。
另一种几何方法基于另一种单纯形体积公式V=1/p?1Vbase×height，这是众所周知的1/2base×height三角形面积公式的推广。由于单个点的替换不会改变单纯形的底边，因此体积变化只需要计算高度（图4 A）。
虽然这种方法清晰直接，但在文献中经常引起混淆。一些论文从几何角度介绍它[23]、[79]、[80]、[81]、[82]，而其他论文则通过使用Gram–Schmidt正交化得到相同的计算结果[14]、[15]、[83]。实际上，同样的算法以不同的名称为人所知，如GSGA（几何单纯形生长算法）[81]和SVMAX（连续体积最大化）[83]。

下载：下载高分辨率图像（310KB）
下载：下载全尺寸图像
图4. 基于高度的体积变化估计演示。（A）单个细胞拉曼图像（2d情况）。可以通过比较高度来比较体积。所有高度都是通过将所有数据点投影到垂直于底边的直线上来计算的。（B）3d情况的示意图，其中高度是通过将点投影到垂直于底边的平面上来计算的。

4.1.5. 维度降维
虽然行列式单纯形体积公式Eq. (2)要求p?1单纯形处于其内在Rp?1空间中，例如，它不允许计算3维（或更高维）空间中的三角形面积。相比之下，基于高度（或Gram–Schmidt）的方法不需要维度降维，正如Du[14]使用Gram行列式所描述的。图4 B提供了几何解释，通过将点投影到底边的零空间来计算更高维度的高度。
尽管可以在不进行维度降维的情况下计算体积，但恢复的端元光谱通常与维度降维后的计算相比并没有改进，并且计算成本更高。数据中的噪声可能会导致在原始m个光谱通道（例如波数）中产生最大体积的点在投影到降维后的p?1维空间（例如p?1个主成分）时不再具有最大值。此外，端元的数量通常远小于光谱通道的数量（p?m）。因此，不进行维度降维地应用N-FINDR意味着要处理更大的矩阵/向量，从而增加计算成本。最早引入不进行维度降维的N-FINDR的研究之一[14]得出结论：“我们发现这个过程[维度降维]实际上可以提高性能，因为它将重要的数据信息压缩到主要的主成分中，并减少了像素相似性评估中微不足道的光谱变化的敏感性。因此，选择具有维度降维过程的端元提取算法可能更好。”

4.1.6. 优化计算
许多研究将其优化动机归因于行列式计算是一个计算成本高昂的操作。由此衍生出了不同的体积变化估计方法。然而，必须澄清的是，行列式的计算并不像替换次数那么重要。通常，端元的数量p远小于数据点/像素的数量（p?n）。因此，计算一个p×p（例如10 × 10）矩阵的行列式并不会产生显著的计算成本。但这种计算在每次替换时都会执行，而这正是性能差异的关键。由此可以看出，使用最内层循环中的端元进行迭代明显较慢。最内层循环的一次运行只能检查p次体积变化，而使用点和两次迭代分别检查n和pn次体积变化。图4A很好地说明了如何通过单次矩阵-向量乘法检查n次体积变化。进行的基准测试也证明了这一点（见第4.3.1节）。
接下来，最重要的是，向量化（用矩阵操作替换循环）在文献中严重不足。在大多数情况下，最内层循环（点、端元或两者）可以用单次矩阵乘法来替换。例如，使用高度体积变化估计器和点迭代时，所有点都可以通过单次矩阵-向量乘法投影到基的法线上，如图4A所示。这消除了逐个迭代每个点的需要。类似的向量化策略也可以应用于其他方法。
向量化利用了高度优化的库、硬件和现代处理器架构，以高效执行标准矩阵和向量操作。即使对于像C这样的低级编程语言也是如此，向量化操作可以利用优化的BLAS和LAPACK例程。尽管文献中忽略了这一点，但这引发了关于其他优化（如并行计算、随机化或块处理[16]、[19]、[74]）的实际效益的疑问。第4.3.1节中的基准测试表明，向量化可以将性能提高一个数量级——与优化行列式计算所获得的收益相当。
一些研究提出将数据集分成子集，并分别对每个子集应用N-FINDR。然后将这些子集的输出组合起来，作为最终N-FINDR步骤的输入，或者在某些变体中，一个子集的输出作为下一个子集的输入[74]。尽管这种方法通常比优化的线性代数计算效率低，但在数据量过大无法放入内存的情况下，它仍然是一个合理的选择。然而，必须明确的是，即使N-FINDR的假设适用于整个数据集，也不能保证这些假设对任何给定的子集都成立。特别是，子集可能无法形成一个有效的单纯形。作为一种解决方法，可以考虑用凸包搜索（如Quickhull算法[84]）来替代N-FINDR在子集中的应用。另一种常见的优化策略是排除那些不可能是端元点或极不可能成为端元点的点——例如那些位于数据云中心的点[8]、[10]。但是，必须确保这种过滤带来的计算节省能够超过过滤本身的成本。[10]提出了一种有趣的方法，即重用已经用于估计体积变化的计算结果来进行过滤。通常，这样的过滤方法需要额外的参数（例如阈值），这些参数需要用户来设置。阈值的选择可能会使工作流程复杂化并影响结果，因为它可能会排除潜在的端元点。此外，过滤的有效性取决于数据：在某些图像中，大多数点可能聚集在中心附近，而在其他图像中，许多点可能位于边缘（例如，当并非所有纯组分相互混合时）。虽然这是一个值得进一步研究的领域，但由于通过向量化已经实现了改进，因此本工作中没有包含过滤方法。

4.1.7. 丰度反演
正如算法介绍中提到的，N-FINDR本身并不提供丰度信息。丰度反演是一个需要单独讨论的主题，超出了本工作的范围。然而，这里有一些重要的方面值得一提。首先，需要澄清的是，丰度反演并不是一个明确的步骤，它至少需要做出以下决策：
1. 用于丰度反演的数据：原始光谱数据、经过不同预处理的光谱数据、降维数据（例如PCA得分）、重建的数据（例如将PCA得分重建回原始空间），或其他类型的数据。
2. 用于丰度反演的计算模型。文献中提到了几种选项，包括重心坐标（对应于混合物图中的坐标）[73]、无约束的最小二乘法（OLS或UCLS）、总和为一的约束最小二乘法（SCLS，即数值之和必须为一，但允许负值）[85]、非负约束最小二乘法（NCLS或NNLS，即数值必须非负但不必总和为一）[86]，以及同时满足两种约束的完全约束最小二乘法（FCLS）[87]、[88]。
3. 如何处理端元点形成的单纯形之外的点，因为大多数计算模型会产生负值或总和不为一的丰度。例如，可以使用一种约束方法并将数值规范化为总和为一，或者简单地决定相应的光谱不符合给定的端元点集合。

因此，对于应用工作来说，明确所使用的丰度反演方法非常重要；对于软件实现来说，提供一个明确且透明的丰度反演步骤也是必要的。例如，在早期的unmixR R包版本（2.0之前）及其使用的相关工作中，并没有这样做。

4.2. 算法表示的标准化
综上所述，建议使用一种不同的方法来描述算法。与其为每个实现生成一个新的缩写或名称，不如将算法分解为一组参数，这些参数能够明确地定义实现方式。这种方法非常标准，并且已经被用于其他成熟的算法，如k-means。此外，所选择的方法也使得引入新的实现及其与现有实现的关系变得更加容易。输入参数及其默认值的列表基于基准测试：
- X：光谱数据矩阵。可以是降维后的（即n×(p?1)）或原始空间的数据。
- p：要搜索的端元点数量。
- init：使用的初始化策略。可能的选项包括：随机、随机投影，或特定索引的向量（例如，用于使用另一个端元点提取算法的输出，如VCA）。
- iterator：迭代策略（也对应于“向量化操作”。可能的选项包括：最内层循环中的端元点、最内层循环中的点（默认值），以及两者兼有。
- estimator：估计体积变化的方法。可能的选项包括：volume - 直接计算体积而不进行任何优化；height - 利用单纯形体积与高度和底面积的乘积成正比的关系；Cramer - 使用Cramer法则；LDU - 使用LDU矩阵分解；cofactor - 使用余子式展开来计算det(E)。如果X没有降维，则只允许使用volume和height（默认值）；否则，默认使用Cramer值。
- iter_max：最外层迭代的最大次数。默认值为10次。这个值基于经验，通常大约五次最外层迭代就足够了。
- n_init：随机初始化的次数。如果有多次初始化，则N-FINDR将运行n_init次，并选择体积最大的单纯形。默认值为一，这符合对N-FINDR的经典理解。如果提供了特定的端元点列表，则忽略此值。

注意，数据预处理步骤被有意地移出了N-FINDR本身，即没有参数控制要应用哪种降维方法（无、PCA、MNF等）或归一化方法。这提供了一定程度的灵活性，因此无需为每种预处理技术更新N-FINDR的实现。这也使得在N-FINDR之外重新使用预处理结果（例如PCA的载荷和得分）变得更加容易。实际上，这也是合理的，因为预处理（包括降维）通常是非常特定于手头数据的。换句话说，从编程接口的角度来看，降维、归一化和其他预处理步骤并不是N-FINDR的一部分。然而，在展示结果时，为了结果的可重复性，也必须明确提及这些步骤。例如，当将k-means应用于光谱数据集时也是如此。

补充表S1提供了大多数已知N-FINDR实现的建议表示方法。虽然不够全面，但在补充表S2中可以找到类似的分解序列算法（VCA、ATGP等）为参数的表格。

4.3. 基准测试
我们对提出的N-FINDR通用版本进行了基准测试，以比较算法中每个参数的不同选项，并确定默认使用的最佳方法。

4.3.1. 迭代和体积变化估计
首先，我们通过在两个数据集（Cuprite和Macrophages）上测试所有可能的组合来比较各种迭代策略和体积变化估计器。计算性能是根据端元点数量和数据集大小（固定为104个光谱）以及数据集大小和端元点数量（固定为5个）来评估的。所有体积变化估计器都进行了向量化处理。不过，有些组合无法向量化，例如高度估计器与端元点迭代。此外，还包括了直接通过行列式公式计算体积的方法和未向量化的LDU估计器（unmixR R包早期版本的默认值）进行比较。由于这些方法的性能明显低于向量化版本，因此没有对所有情况都进行测试。尽管如此，提供的结果足以展示性能差异。在评估之前，随机种子被固定以保持相同的初始点。

图5展示了根据端元点数量测试的基准结果。正如预期的那样，最内层循环中的端元点迭代器比其他两个迭代器慢得多。其中一个原因是，对端元点的向量化只消除了p次比较，但需要更多的替换操作（见补充图S3）。此外，端元点迭代没有充分利用向量化，因为只使用了一个逻辑处理器核心进行计算。相比之下，最内层循环中的点和两者都利用了所有可用的22个核心。有趣的是，尽管需要较少的替换操作，但两者迭代器的性能仍然较慢。这可能是由于两个主要原因：首先，替换操作的减少不如端元点迭代那么显著；其次，与嵌套循环计算不同，向量化计算需要额外的内存来存储相对体积变化。对于两者迭代器来说，这意味着需要存储一个n×p矩阵，其中每个元素i,j表示用第i个数据点替换第j个端元点时的体积变化。而对于点迭代器，只需要一个长度为n的向量。

图S3还表明，在最终单纯形体积方面没有明显的赢家。

下载：下载高分辨率图像（974KB）
下载：下载全尺寸图像

图5. 在两个数据集（Cuprite（上排）和Macrophages（下排）上比较迭代器的基准结果。性能是根据端元点数量来测试的，数据集大小固定为10000个光谱。每一列代表相应的体积变化估计器的比较结果。每个点对应20到30次重复计算的中位时间。时间轴以对数刻度显示，以便更好地可视化。所有迭代器都使用了向量化实现，除了“points (loop)”，它对应于非向量化的LDU-N-FIDNR实现。为了对比，虚线表示直接使用体积公式时的性能。灰色垂直虚线突出显示了参考端元点数量（来自原始工作的估计值）在相应数据集中的情况。

为了比较体积变化估计器的计算性能，补充图S4按估计器分组展示了相同的结果。此外，补充图S5提供了根据数据集大小的结果。结果强调了向量化的重要性。首先，向量化的LDU实现比其非向量化版本快一个数量级。其次，Cramer法则由于其易于向量化而表现出最佳的整体性能。这种优势在最内层循环中有端元点的情况下最为明显，因为某些实现无法向量化。即使在其他迭代器中，Cramer法则仍然更优越，因为它需要的额外计算较少。例如，使用点迭代器时，高度估计器在切换端元点时必须重新计算法向量，而Cramer方法则重用了之前计算出的逆矩阵。随着数据集大小的增加，这部分计算的重要性会降低。

检查最内层循环中使用点的LDU估计器（图5），可以看出向量化带来的计算改进与优化体积计算带来的改进相当。向量化和非向量化LDU估计器之间的中位数加速比大约为15倍，而非向量化LDU和直接体积计算之间的加速比为23倍。综合来看，总体加速比约为380倍。

4.3.2. 最终单纯形体积和端元点数量
基准测试结果还提供了关于最终单纯形体积与端元点数量之间关系的有趣见解。图6展示了使用随机初始化进行30次N-FINDR运行后获得的最终单纯形体积，这些运行针对不同的数据集、端元点数量和迭代策略组合。它清楚地表明，最终单纯形体积随着端元点数量的增加而增加。这主要是由于使用更多主成分时引入的噪声增加，以及在更高维度空间中操作的影响（例如，每个顶点处5%的偏差会导致体积呈指数级增长，比例为1.05p）。最终结果可能只有多次运行获得的最大体积的25%。这表明，当端元点数量较多时，使用多次初始化是必要的。

值得注意的是，在接近估计的真实端元点数量的情况下，任何指标（替换次数或最终单纯形体积）都没有显著差异（见图S3中的替换次数）。因此，无论是替换次数还是最终单纯形体积都不能用来确定最佳端元点数量。补充模拟通过在具有已知真实值的合成数据上展示不同的端元点数量下的收敛和体积最大化趋势，支持了这些发现。进一步研究表明，较大的单纯形体积并不一定能够产生最准确或最理想的分量谱——尤其是在底层假设被严重违反的情况下。尽管如此，它仍然作为一种有用的启发式方法来提取有意义的分量谱。下载：下载高分辨率图像（939KB）下载：下载全尺寸图像图6. 不同迭代次数和分量数量下的最终单纯形体积。每种条件（由数据集、分量数量和迭代策略定义）在30次N-FINDR运行中进行了评估，每次运行都使用克拉美法则来估计体积变化。迭代策略用颜色编码：最内层的分量（金色）、最外层的点（蓝色）以及两者都有（红色），结果分别显示在不同的列中。体积是相对于每个数据集和分量数量的最大值来表示的（即所有迭代次数合并后的结果）。每个点对应一次N-FINDR运行；黑色点和线表示中值。灰色的垂直虚线突出显示了相应数据集中分量数量的参考值（来自原始工作的估计）。（关于此图例中颜色的解释，请参阅本文的网页版本。）4.3.3. 维度降维为了阐明维度降维的影响，对两个数据集进行了基准测试。图7显示了使用PCA进行维度降维与否的N-FINDR的计算时间和内存分配情况。两种实现都使用了相同的N-FINDR设置：最内层循环中的点和向量化的高度体积变化估计器，并且选择了相同的随机初始像素。PCA的计算时间被包含在总时间中以进行公平比较。此外，还基准测试了不使用N-FINDR的PCA，以显示其对总计算时间的贡献。结果清楚地表明，即使包括PCA的计算时间，使用维度降维的N-FINDR也优于不使用降维的版本。通过考虑内存分配可以更好地理解这一差距。原始光谱空间通常比降维后的PCA数据大得多（p?m），导致所有计算的内存使用量成比例增加。对于选定的单个巨噬细胞示例数据，使用20次随机初始化运行N-FINDR时，使用维度降维的N-FINDR所需的计算时间减少了14倍（包括PCA时间）。尽管图像的大小相对较小（600个光谱）。根据结果，虽然理论上可以在不进行维度降维的情况下运行N-FINDR，但实际上这是不合理的。补充模拟结果也没有显示出在原始光谱空间中操作有任何优势。因此，推荐的用法是使用维度降维，因为它更高效，并且由于噪声减少而产生相似甚至更好的结果。总体而言，这些发现与之前的研究[14]一致。另一个有趣的观察是当N-FINDR应用维度降维时PCA所占的时间比例。基准测试结果显示，PCA计算可能占总计算时间的90%以上。这清楚地表明，一些N-FINDR优化对总分析时间的实际贡献很小。下载：下载高分辨率图像（628KB）下载：下载全尺寸图像图7. 三种方法的中值计算时间（上排）和中值内存分配（下排）的比较：仅使用PCA（红色）、使用PCA的N-FINDR（绿色）以及不使用PCA的N-FINDR（蓝色）。结果显示了两个数据集：Cuprite（左列）和巨噬细胞（右列）。为了清晰起见，两个轴都使用了对数刻度。N-FINDR使用了20次随机初始化，采用了最内层循环中的“点”和向量化的“高度”体积变化估计器。（关于此图例中颜色的解释，请参阅本文的网页版本。）4.3.4. 初始化和顺序方法我们比较了以下初始化策略在速度和输出单纯形体积方面的表现：•VCA - 使用向量成分分析的输出作为N-FINDR的初始分量。•ATGP - 同样，使用自动目标生成程序（ATGP）的输出。•N-FINDR - 为了公平比较，还考虑了最外层N-FINDR循环的第一次迭代（即最大迭代次数1）作为初始化方法。•极端坐标 – 依次选择具有最极端坐标值的像素：首先选择PC1的最小和最大值像素；接下来选择PC2的最小和最大值像素；这个过程继续进行，直到选择了所需数量的唯一像素p。•极端坐标（随机） – 以随机方式选择具有极端坐标值的像素：对于每个p-1个主成分，识别出最小和最大值像素，然后从这个池中随机选择p个唯一像素。•随机投影 – 通过将数据投影到随机生成的向量上来选择极端点，如图2A所示。对于每个随机向量，选择对应于最小和最大投影值的像素。这个过程重复进行，直到获得所需数量的唯一像素p。VCA和ATGP代表了顺序方法的使用。图3展示了这些初始化策略应用于单个巨噬细胞示例数据的情况。每种初始化方法运行了30次，测量了关键指标：初始化后的单纯形体积、最终单纯形体积（即N-FINDR的最终输出）、初始化和后续N-FINDR所花费的时间，以及初始化后进行的替换次数。所应用的N-FINDR实现使用了维度降维，最内层循环中的点，以及向量化的克拉美法则来估计体积变化。图8和补充图S6展示了所选初始化方法在两个基准数据集上的时间性能，数据集中的分量数量不同。可以预见，选择极端坐标是最快的，因为它不需要任何额外的计算。VCA和随机投影的表现相似，并且比最外层N-FINDR循环的第一次迭代更快。图8还表明，在问题的大规模上，初始化时间可能并不重要。尽管初始化时间存在一些相对差异，但总时间差异可以忽略不计。初始化后的计算占用了总时间的99%。当然，所考虑的数据集并不涵盖所有可能的情景。也许可以找到一个初始化对最终结果贡献更大的情况。最明显的情况是数据集接近理想情况，因为在初始化后需要的替换较少（或没有）。下载：下载高分辨率图像（610KB）下载：下载全尺寸图像图8. 比较不同初始化方法的步骤所花费的时间。从左到右的列分别对应：初始化时间、初始化后的N-FINDR时间、总时间（即前两者的总和）。结果展示了两个数据集：Cuprite（上排）和巨噬细胞（下排）。N-FINDR的设置对所有初始化方法都是相同的：最内层循环中的“点”，向量化的*Cramer*体积变化估计器。灰色的垂直虚线突出显示了相应数据集中分量数量的参考值（来自原始工作的估计）。比较初始单纯形体积和最终单纯形体积，图9显示两者之间没有相关性。补充图S7提供了更详细的结果。虽然VCA产生的单纯形体积几乎与仅进行一次分量替换后停止N-FINDR的结果一样好，但所有顺序方法返回的单纯形体积明显小于完全迭代N-FINDR所能达到的体积。这一发现与之前的文献完全一致，这些文献指出不建议早期停止N-FINDR[12]。唯一一种ATGP、VCA和单次迭代N-FINDR立即返回适当解的情况是巨噬细胞数据集，该数据集确实具有非常明显的三角形结构。然而，在这种情况下，完全迭代N-FINDR也将在再迭代一次后检测到收敛。这清楚地显示了N-FINDR相对于顺序方法的优势。下载：下载高分辨率图像（687KB）下载：下载全尺寸图像图9. 比较不同算法的初始单纯形体积和最终单纯形体积。不同算法实现的单纯形体积（中位数：点；四分位数范围：带状线），以及与每个数据集在30次迭代N-FINDR收敛后达到的最大单纯形体积的第一次迭代（“N-FINDR（1）”）进行比较。ATGP和极端坐标是确定性的，因此不会显示随机变化。使用随机化的好处是可以找到更大的单纯形体积，这与ATGP和极端坐标等确定性方法相比具有优势。在随机化方法（即随机投影、VCA和N-FINDR）中，没有发现最终单纯形体积的显著差异。对于更大的最终单纯形，结果表明使用多次随机初始化比使用特定的初始化方法本身更为重要。总体而言，结果表明最佳实践是进行多次随机初始化。随着分量数量的增加，这种方法变得越来越重要，因为最终单纯形体积的变异性也会增加。在随机化方法之间没有观察到明显的优越性。从技术角度来看，正交投影方法（如VCA和ATGP）与N-FINDR的集成更为复杂，因为它们需要不同的数据维度（VCA为n×p，而N-FINDR为n×(p?1)）。尽管如此，所提出的初始化参数可以接受一个明确的初始索引向量，使得VCA（或任何其他纯像素分量提取算法）能够与N-FINDR无缝连接——使用前者的输出作为后者的输入。4.3.5. Python实现最后，我们比较了N-FINDR的Python实现的计算性能。在Cuprite数据集上进行了基准测试，使用了最内层循环中的点、向量化的克拉美体积变化估计器，并进行了30次重复。作为对比，还基准测试了另一种广泛使用的Python实现。该实现基于‘pysptools’包[89]，该包提供了一系列光谱解混算法，包括N-FINDR。‘pysptools’实现使用了最内层循环中的点和体积估计器（即直接行列式计算）。结果显示在图S8中。作为参考，还包括了R实现的相应结果（来自图5）。正如预期的那样，向量化的Python实现明显快于非向量化的Python版本。有趣的是，非向量化的Python实现比其R对应版本更快。总体而言，这些结果与前面章节中的结论和讨论一致。5. 结论在这项研究中，我们全面概述了N-FINDR算法，旨在阐明其基本原理、操作机制及其与其他纯像素分量提取算法（如VCA）的关系。深入分析了N-FINDR的各种实现之间的差异，并进行了说明。因此，我们提出了该算法的通用版本，用明确的值替换了缩写，从而实现了更透明的表示。所提出的通用化有助于清晰理解现有实现及其之间的关系。此外，它简化了新实现的引入。算法的每个组成部分都通过一般分析和基准测试进行了广泛讨论和比较。基于这些比较，我们提供了最佳默认实现。此外，还开发了开源的Python和R包，这些包基于本文中的发现和讨论构建。这些包使非专家能够轻松使用最佳实现来应用该算法，而无需深入研究其复杂细节。同时，算法开发者可以使用这些包将其想法与当前可用的实现进行清晰比较。结果表明，推荐的N-FINDR使用方法是首先进行维度降维，使用最内层循环中的点，并从第一个分量开始重新启动整个过程；使用向量化的克拉美法则方法估计体积变化；并使用不同的随机初始化多次运行算法。初始化方法本身并不重要，只要它是随机的即可。然而，需要注意的是，在实践中，输出结果的质量在很大程度上取决于底层假设的满足程度——尤其是纯度和噪声水平。补充模拟展示了如何仅通过PCA得分图就能对输出形成现实的预期，并评估特定假设的违反如何影响最终结果。关于VCA [20]、ATGP等投影方法，我们认为已经清楚地证明了这些方法在设计上找到的单纯形的体积要么小于N-FINDR，要么最多与其相等。也就是说，在这两种算法家族背后的主要启发式方法中，它们并不比N-FINDR表现得更好。因此，在对这些方法与N-FINDR进行（运行时间）性能比较时应该谨慎对待。由于现在已知有高效的N-FINDR算法，N-FINDR的迭代优化解决方案可以在相对较低的计算成本下获得优势。因此，在应用研究中，除了评估N-FINDR的结果外，很少需要再评估其他投影方法的结果。

在与VCA进行比较时，区分维度缩减步骤的效果和后续搜索策略的效果是很重要的。如果观察到VCA在某些情况下优于N-FINDR，那么研究这一结果背后的原因将非常有价值。特别是，这种情况可能有助于了解单纯形体积与最终端元光谱质量之间的关系。

我们想指出，在向量化之后，优化的重要性其实很小。一旦数据被向量化，维度缩减可能会占用端元提取总时间的90%，而避免维度缩减并不会提高整体性能。鉴于此，我们有理由质疑优化的实际益处。显然，这里引用唐纳德·克努特的一句名言很合适：“毫无疑问，对效率的追求会导致滥用。程序员会浪费大量时间去思考或担心程序中非关键部分的速度……我们应该忘记那些小效率，比如97%的时间：过早的优化是一切问题的根源。然而，我们不应该放弃那3%的关键机会。一个好的程序员不会被这种推理所麻痹，他会明智地仔细检查关键代码；但只有在确定了关键代码之后才会这样做。”[90]

与其仅仅专注于微调算法的特定方面，不如采取一种更有效的方法来进一步发展它，即解决其在应用过程中遇到的实际挑战。这包括通过评估算法假设的满足程度来验证其适用性，并考虑其单纯形的特性，调整算法以适应预处理步骤（如基线校正）的影响，以及处理与噪声相关的问题（如宇宙射线尖峰）。开发估计端元数量的方法是未来研究的一个重要领域。虽然目前已有虚拟维度（VD）等方法，但仍有很大的改进空间。来自化学秩估计领域的进展应该被纳入其中，并可能在这方面带来益处。此外，丰度反演也是一个需要进一步研究的领域。文献中对此主题的讨论相对较少。对不同方法及其影响的深入分析将有助于建立丰度反演的最佳实践，从而填补分析、解释和结果可视化方面的空白。解决这些挑战将提高算法的整体效率和适用性。然而，必须认识到这需要专门且全面的研究，通常需要与非纯像素解混方法的进步相结合。总之，我们乐观地认为这项工作将有助于无缝应用最佳实践，避免算法用户常见的陷阱，并帮助算法开发者明确现有实现与新实现之间的关系。

**作者贡献声明：**
Rustam Guliev：撰写——审阅与编辑、撰写——初稿、可视化、验证、软件开发、调查、形式分析、概念化。
Ute Neugebauer：撰写——审阅与编辑、验证、监督、资源获取、资金筹集。
Claudia Beleites：撰写——审阅与编辑、验证、监督、方法论、概念化。

**关于写作过程中生成式AI和AI辅助技术的声明：**
在准备这项工作时，作者使用了ChatGPT来提高手稿的可读性和语言表达。使用该工具后，作者根据需要对内容进行了审阅和编辑，并对发表文章的内容负全责。