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

醉雨轩^_^

于 2020-04-16 19:05:27 发布

阅读量1.6k

点赞数

分类专栏：算法信号处理

本文链接：https://blog.csdn.net/miao0967020148/article/details/105564781

版权

信号处理同时被 2 个专栏收录

25 篇文章 0 订阅

订阅专栏

算法

12 篇文章 3 订阅

订阅专栏

考虑线性系统

$\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
点赞
踩
8

收藏

觉得还不错? 一键收藏
0
评论
OLS 正交最小二乘算法：稀疏信号重建

考虑线性系统寻找解决方案的形式为：如果数据向量在的列空间中，则上式具有精确解：如果具有零空间分量，我们将无法满足差分方程，因此我们放宽了要求，并要求尽可能小。这导致最小二乘问题的定义：一般最小二乘法解是有多种解决方案。例如，提到的正规方程式：提供解决方案:正规方程的构造，右边的向量，显然在的列空间中。实际上，我们给出了组合列向量的方案...
复制链接

扫一扫