SPA是一种经典的特征变量选择算法，广泛应用于光谱分析（如近红外、红外光谱）中。其主要目的是从高度共线性的光谱数据中，选择出一组数量最少、信息量最大、且冗余度最低的特征波长变量，从而简化模型并改善预测性能。
一、 算法核心思想
SPA通过一系列的投影操作，从一个初始波长开始，迭代地寻找与之前已选波长向量共线性最小的新波长点。
- 初始化： 选择一个初始波长（通常是光谱信号强度最大的点，或者随机选择/手动指定）。
- 投影操作： 将剩余的未选波长点向量投影到已选波长向量所张成的子空间的正交补空间上。
- 选择最大投影点： 在投影后的向量中，选择投影向量模长最大的那个波长点。模长大意味着该点携带的信息与已选点集的信息冗余度最小。
- 迭代： 将新选的点加入已选点集，重复步骤2和3，直到选取了预定数量
K的波长点。
简单理解：就好比你在已有几个观点的前提下，SPA帮你找一个与现有观点最不重复、最能提供新信息的新观点。
二、 讲解
函数定义：selected_variables = spa(X, K, initial_var)
- 输入：
X: 输入数据矩阵，大小为n samples x p variables（例如，n个样本的光谱，每个光谱有p个波长点）。K: 想要选择的变量（波长）数量。initial_var: （可选）指定的初始变量索引。如果未提供或为0，则自动选择信号范数最大的点作为起点。
- 输出：
selected_variables: 一个包含最终选择的K个变量索引的向量。
function selected_variables = spa(X, K, initial_var)% SPA - 连续投影算法% 输入：% X - 数据矩阵 (n x p)% K - 需要选择的变量数% initial_var - (可选) 初始变量位置，默认为0（自动选择）% 输出：% selected_variables - 被选中的变量索引向量[n, p] = size(X); % 获取矩阵大小selected_variables = zeros(1, K); % 初始化输出向量% 步骤 1: 数据标准化（中心化），消除基线影响X_c = X - mean(X, 1);% 步骤 2: 选择初始变量if nargin < 3 || initial_var == 0% 计算每个波长点信号的欧几里得范数（强度）norms = sqrt(sum(X_c.^2, 1));[~, initial_index] = max(norms); % 找到范数最大的那个点的索引selected_variables(1) = initial_index;else% 使用用户指定的初始变量selected_variables(1) = initial_var;end% 初始化工作向量和投影矩阵% 从所有变量索引中移除已选的第一个变量candidate_set = 1:p;candidate_set(selected_variables(1)) = [];% 初始化一个投影向量，开始时为初始变量对应的光谱列x = X_c(:, selected_variables(1));iter = 2; % 从第二个变量开始迭代% 步骤 3: 迭代选择剩余的 K-1 个变量while iter <= K% 3.1 从候选集中提取出所有候选变量的光谱数据X_j = X_c(:, candidate_set);% 3.2 核心操作：正交投影% 计算投影算子 P = I - x*(x'x)^-1*x'% 这个算子可以将任何向量投影到x的正交补空间上P = eye(n) - (x * x') / (x' * x);% 将所有的候选变量向量进行投影PX = P * X_j;% 3.3 计算投影后向量的范数（模长）norms_PX = sqrt(sum(PX.^2, 1));% 3.4 选择投影范数最大的那个候选变量的索引[~, max_index] = max(norms_PX);% 3.5 更新选中变量列表和候选集% 注意：max_index 是 candidate_set 中的索引，不是原始索引selected_index = candidate_set(max_index);selected_variables(iter) = selected_index;% 从候选集中移除已选变量candidate_set(max_index) = [];% 3.6 更新投影向量为最新选中的变量向量，用于下一次迭代x = X_c(:, selected_index);iter = iter + 1;end
end
三、如何使用该函数
以下代码演示如何在一个模拟的光谱数据集上使用SPA函数。
% 示例：使用SPA选择特征波长
clear; clc;% 1. 生成模拟光谱数据
nSamples = 100;
nWavelengths = 500;
X = randn(nSamples, nWavelengths); % 随机噪声作为基线% 模拟几个高斯峰作为有效信号
peak_positions = [100, 200, 300, 450]; % 峰的位置
peak_intensity = [5, 8, 3, 6]; % 峰的强度
for i = 1:length(peak_positions)pos = peak_positions(i);intensity = peak_intensity(i);% 在每个峰位置附近添加信号X = X + intensity * exp(-((1:nWavelengths) - pos).^2 / (2*10^2));
end% 2. 应用SPA算法选择10个特征波长
K = 10;
selected_idx = spa(X, K, 0); % 自动选择初始点% 3. 可视化结果
figure;
% 3.1 绘制平均光谱和被SPA选中的波长点
plot(mean(X, 1), 'b-', 'LineWidth', 1.5);
hold on;
plot(selected_idx, mean(X(:, selected_idx), 1), 'ro', 'MarkerFaceColor', 'r', 'MarkerSize', 8);
title('SPA特征波长选择结果');
xlabel('波长点索引');
ylabel('吸光度/强度');
legend('平均光谱', ['SPA选中的', num2str(K), '个特征波长'], 'Location', 'best');
grid on;
hold off;% 4. 输出选中的波长索引
disp('SPA选中的特征波长索引为：');
disp(selected_idx);
参考代码 连续投影算法 www.youwenfan.com/contentcnh/55071.html
代码解释：
- 生成数据： 创建了一个带有随机噪声和4个模拟高斯峰的光谱数据集。
- 调用SPA： 使用上面编写的
spa函数从500个波长点中选择10个最重要的。 - 可视化： 绘制平均光谱图，并在图上用红色圆圈标记出SPA算法选中的波长点。你会发现，这些点很可能就位于我们事先设置的4个高斯峰及其相关区域附近。
- 输出结果： 在命令窗口打印出选中波长的索引号。