l1范数最小化求解系数方程:正交匹配追踪法（orthogonal matching pursuit）

正交匹配追踪算法

Step 1 初始化 令标签集Ω0=$\Omega_0=\varnothing$,初始残差向量r0=y$r_0=y$，令k=1$k=1$.
Step 2 辨识 求矩阵A$A$中与残差向量rk1$r_{k-1}$最强相关的列
jk$j_k\in$argmaxj|<rk1,ϕj>|,Ωk=Ωk1jk$_j ||, \Omega_k=\Omega_{k-1}\cup {j_k}$
Step 3 估计 最小化问题 minxyAΩkx$_x\|y-A_{\Omega_k}x\|$的解由
xk=(AHΩkAΩk)1AHΩky$x_k=(A_{\Omega_k}^H A_{\Omega_k})^{-1}A_{\Omega_k}^Hy$给出，其中
AΩk=[aω1,...,aωk],ω1,...,ωkΩk$A_{\Omega_k}=[a_{\omega_1},...,a_{\omega_k}],\omega_1,...,\omega_k \in \Omega_k$
Step 4 更新残差
rk=yAΩkxk$r_k=y-A_{\Omega_k}x_k$
Step 5k=k+1$k=k+1$，并重复Step 2Step 4。若某个停止判据满足，则停止迭代。
Step 6 输出系数向量
x(i)=xk(i),$x(i)=x_k(i),$ if iΩk$i\in\Omega_k$
x(i)=0,$x(i)=0,$ otherwise

1. 运行到某个固定的迭代步数后停止。
2. 残差能量小于某个预先给定值ε$\varepsilon$
rk2ε$\| r_k \|_2 \leq \varepsilon$
3. 当字典矩阵A$A$的任何一列都没有残差向量rk$r_k$的明显能量时
AHrkε$\| A^H r_k \|_\infty \leq \varepsilon$

Matlab 实现

function [x, Out] = OMP(A,y,varargin)
% This is a l1 minimization solver for orthogonal matching pursuit:
%   minimize     ||x||_1
%   subject to   y = Ax,
% -----------------------------------------------------------------
% Author: Beier ZHU
%         Tsinghua University
% -----------------------------------------------------------------
%
% =========================Required inputs=========================
% A -- an m x n matrix
% y -- an m x 1 vector
% =========================Optional inputs=========================
% 'maxIter' -- maximum number of iterations
% 'StopTolerance' -- stopping tolerance
% ===========================Outputs===============================
% x -- last iterate (hopefully an approximate solution)
% Out.iter -- # of iterations taken

% Test for number of required parametres
if (nargin-length(varargin)) ~= 2
error('Wrong number of required parameters');
end
%--------------------------------------------------------------
% Set parameters to their defaults
%--------------------------------------------------------------
opts.tol = 1e-3;
opts.maxIter = 1e3;
%--------------------------------------------------------------
% Set parameters to user specified values
%--------------------------------------------------------------
if (rem(length(varargin),2)==1)
error('Options should be given in pairs');
else
for i=1:2:(length(varargin)-1)
switch upper(varargin{i})
case 'STOPTOLERANCE'
opts.tol = varargin{i+1};
case 'MAXITER'
opts.maxit = varargin{i+1};
otherwise
error(['Unrecognized option: ''' varargin{i} '''']);
end
end
[m_A, n] = size(A);
[m_y, n_y] = size(y);
if n_y ~= 1
error(['y must be a column vector']);
end
if m_A ~= m_y
error(['A and y must have same # of row'])
end
% --------------------------------------------------------------
x = zeros(n,1);
Omega = [];
A_Omega = [];
r = y;
iter = 1;
varepsilon = 1e-3; % step 1

while true
[~, max_index] = max(abs(A'*r));
Omega = [Omega max_index];
A_Omega = [A_Omega A(:,max_index)]; % step 2
x_k = A_Omega\y; % step 3
r = y - A_Omega*x_k; % step 3
iter = iter + 1; % step 4
if (iter > opts.maxIter) || (norm(r) <= opts.tol) || (norm(A'*r, inf) <= varepsilon)
break; % 终止条件
end
end
x(Omega) = x_k; %step 5
Out.iter = iter;
end
'''

[0]: 矩阵分析与应用 第二版 张贤达