考虑线性系统
寻找解决方案的形式为:
如果数据向量在的列空间中,则上式具有精确解:
如果具有零空间分量,我们将无法满足差分方程,因此我们放宽了要求,并要求尽可能小。
这导致最小二乘问题的定义:
一般最小二乘法解是
有多种解决方案。 例如,正规方程式:
提供解决方案:
正规方程的构造,右边的向量,显然在的列空间中。 实际上,我们给出了组合列向量的方案:
现在,我们可以使用高斯消除和分解之类的工具来解决该问题。
在数值问题中,正规方程可能是一个较差的选择,而和奇异值分解(SVD)等方法较为好。
用SVD下的基础投影,投影在范围空间:
算子将个向量投影到范围空间(列空间,图像)上。
SVD奇异值分解为
这意味着伪逆是
正规方程的解是
MATLAB的代码:
function [s, residual] = OLS(A, y, m, err)
% Orthogonal Least Squares [1] for Sparse Signal Reconstruction
% Input
% A = N X d dimensional measurment matrix
% y = N dimensional observation vector
% m = sparsity of the underlying signal
% Output
% s = estimated sparse signal
% r = residual
% [1] T. Blumensath, M. E. Davies; "On the Difference Between Orthogonal
% Matching Pursuit and Orthogonal Least Squares", manuscript 2007
if nargin < 4
err = 1e-5;
end
s = zeros(size(A,2),1);
r(:,1) = y; L = []; Psi = [];
normA=(sum(A.^2,1)).^0.5;
NI = 1:size(A,2);
for i = 2:m+1
DR = A'*r(:,i-1);
[v, I] = max(abs(DR(NI))./normA(NI)');
INI = NI(I);
L = [L' INI']';
NI(I)=[];
Psi = A(:,L);
x = Psi\y;
yApprox = Psi*x;
r(:,i) = y - yApprox;
if norm(r(:,end)) < err
break;
end
end
s(L) = x;
residual = r(:,end);
例子:
% Demo
% Generate Sparse Signal
signal = GenSig(150, 10);
% Create Measurement Matrix
A=matrix_normalizer(randn(100,150));
% Create Measured Signal
y = A*signal;
figure
% reconstruct Signal by OLS
[s1, residual] = OLS(A, y, 10);
stem(signal), hold on, stem(s1,'r+')
结果: