OLS 正交最小二乘算法:稀疏信号重建

考虑线性系统

\mathbf{A} x = b

寻找解决方案的形式为:

\mathbf{A} x - b = 0.

如果数据向量b\mathbf{A}的列空间中,则上式具有精确解:

x = \mathbf{A}^{-1}b.

如果b具有零空间分量,我们将无法满足差分方程,因此我们放宽了要求,并要求\mathbf{A}x - b尽可能小。

这导致最小二乘问题的定义:

x_{LS} = \left\{ x\in\mathbb{C}^{n} \colon \arg\min\lVert \mathbf{A} x_{LS} - b \rVert_{2}^{2} \right\}

一般最小二乘法解是

x_{LS} = \mathbf{A}^{-1} b + \left( \mathbf{I}_{n} - \mathbf{A}^{-1}\mathbf{A} \right)y, \quad y \in\mathbb{C}^{n}.

有多种解决方案。 例如,正规方程式:

\mathbf{A}^{T} \mathbf{A} x = \mathbf{A}^{T} b

提供解决方案:

x_{LS} = \left( \mathbf{A}^{T}\mathbf{A} \right)^{-1} \mathbf{A}^{T} b

正规方程的构造,右边的向量\mathbf{A}^{T}b,显然在\mathcal{A}^{T}的列空间中。 实际上,我们给出了组合列向量的方案:

b_{1} \left[ \mathbf{A}^{T} \right]_{1} + b_{2} \left[ \mathbf{A}^{T} \right]_{2} + \dots b_{m} \left[ \mathbf{A}^{T} \right]_{m}.

现在,我们可以使用高斯消除和\mathbf{L}\mathbf{U}分解之类的工具来解决该问题。

在数值问题中,正规方程可能是一个较差的选择,而\mathbf{Q}\mathbf{R}和奇异值分解(SVD)等方法较为好。

用SVD下的基础投影,投影在\mathbf{A}范围空间:

\mathbf{P}_{\color{blue}{\mathcal{R}\left( \mathbf{A} \right)}} = \mathbf{A}\mathbf{A}^{-1}

算子将m个向量投影到范围空间(列空间,图像)上。

SVD奇异值分解为

\mathbf{A} = \mathbf{U}\, \Sigma \, \mathbf{V}^{T} = \underbrace{\left[ \begin{array}{cc} \color{blue}{\mathbf{U}_{\mathcal{R}}} & \color{red}{\mathbf{U}_{\mathcal{N}}} \end{array} \right]}_{m\times m} \underbrace{\left[ \begin{array}{c} \mathbf{S}_{\rho \times \rho} \\\mathbf{0} \end{array} \right]}_{m\times n} \underbrace{\left[ \begin{array}{c} \color{blue}{\mathbf{V}_{\mathcal{R}}} \end{array} \right]^{T}}_{n\times n} .

这意味着伪逆是

\mathbf{A}^{-1} = \mathbf{V}\, \Sigma^{-1} \mathbf{U}^{T} = \underbrace{\left[ \begin{array}{c} \color{blue}{\mathbf{V}_{\mathcal{R}}} \end{array} \right]}_{n\times n} \underbrace{\left[ \begin{array}{cc} \mathbf{S}_{\rho \times \rho}^{-1} & \mathbf{0} \end{array} \right]}_{n\times m} \underbrace{\left[ \begin{array}{cc} \color{blue}{\mathbf{U}_{\mathcal{R}}} & \color{red}{\mathbf{U}_{\mathcal{N}}} \end{array} \right]^{T}}_{m\times m}

正规方程的解是

\begin{align*} \left( \mathbf{A}^{T}\mathbf{A} \right)^{-1} \mathbf{A}^{T} &= \left( \left( \mathbf{V}\Sigma^{\mathrm{T}} \mathbf{U}^{T} \right) \left( \mathbf{U}\Sigma\mathbf{V}^{T} \right) \right)^{-1} \mathbf{V}\Sigma^{\mathrm{T}} \mathbf{U}^{T} \\ &=\left( \mathbf{V}\, \mathbf{S}^{2} \mathbf{V}^{T} \right)^{-1} \mathbf{V}\, \Sigma^{\mathrm{T}} \mathbf{U}^{T} \\ &=\left( \mathbf{V}\, \mathbf{S}^{-2} \mathbf{V}^{T} \right) \mathbf{V}\, \Sigma^{\mathrm{T}} \mathbf{U}^{T} \\ &=\mathbf{V}\Sigma^{-1} \mathbf{U}^{T} \\ &=\mathbf{A}^{-1} \end{align*}

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+')

结果:

  • 0
    点赞
  • 9
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值