l1範數最小化求解係數方程:正交匹配追蹤法（orthogonal matching pursuit）

本文部分參考[0]

目錄

貪婪演算法

正交匹配追蹤法是訊號處理中擬合稀疏模型的一種貪婪分佈最小二乘法(greedy stepwise least squares)法。

貪婪演算法(greedy algorithm)的基本思想是：不求整體最優解，而是試圖儘快找到某種意義上的區域性最優解。貪婪法雖然不能對所有問題得到整體最優解，但對範圍相當廣泛的許多問題能產生整體最優解或者整體最優解的近似解。和使用凸優化演算法的l1範數相比，貪婪演算法在速度上具有很大優勢。典型的貪婪演算法有匹配追蹤法和正交匹配追蹤法：

匹配追蹤法

匹配追蹤法(matching pursuit， MP) 是由法蘭西大牛Mallat於1993年提出。其基本思想是，不針對某個代價函式進行最小化，而是考慮迭代地構造一個稀疏解x: 只使用字典矩陣A的少數列向量的線性組合對觀測向量x實現稀疏逼近Ax=y，其中字典矩陣A被選擇的列向量所組成的集合是以逐列的方式建立的。在每一步迭代，字典矩陣中通當前殘差向量r=Ax−y中最相似的列向量被選擇作為作用集的新一列。如果殘差隨著迭代的進行遞減，則可以保證演算法的收斂。

正交匹配追蹤

匹配追蹤只能保證殘差向量與每一步迭代所選的字典矩陣列向量正交，但與以前選擇的列向量一般不正交。正交匹配追蹤(orthogonal matching pursuit, OMP)則能保證每步迭代後殘差向量與以前選擇的所有列向量正交，以保證迭代的最優性，從而減少了迭代次數，效能也更穩健。

正交匹配追蹤演算法

輸入 觀測資料向量y∈Rm和字典矩陣A∈Rm×n.
輸出 稀疏係數向量 x∈Rn.
Step 1 初始化 令標籤集Ω0=∅,初始殘差向量r0=y，令k=1.
Step 2 辨識 求矩陣A中與殘差向量rk−1最強相關的列
jk∈argmaxj|<rk−1,ϕj>|,Ωk=Ωk−1∪jk
Step 3 估計 最小化問題 minx∥y−AΩkx∥的解由
xk=(AHΩkAΩk)−1AHΩky給出，其中
AΩk=[aω1,...,aωk],ω1,...,ωk∈Ωk
Step 4 更新殘差
r

注: AH為共軛轉置

下面是三種常用的停止判據
1. 執行到某個固定的迭代步數後停止。
2. 殘差能量小於某個預先給定值ε。
∥rk∥2≤ε
3. 當字典矩陣A的任何一列都沒有殘差向量rk的明顯能量時
∥AHrk∥∞≤ε

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]: 矩陣分析與應用 第二版 張賢達