l20范数最小化求解系数方程_贪婪组稀疏方法(Greedy group sparsity)

本文章部分参考Fast group sparse classification

前言

之前的一篇文章介绍了l0范数最小化求解方法中的正交追踪匹配法，其基本思想是不求整体的最优解，而是快速得到某个局部最优解。求解l0-minimization的稀疏解在人脸识别中有很广泛的应用，最著名的文章之一便是John Wright和Ma Yi发表在PAMI上的 Robust Face Recognition via Sparse Representation，但是利用l0-minimization的求解方法在面对不同类之间高度相关的情况下会产生错误，最稀疏的解不一定是分类问题最优解，于是Ehsan Elhamifar提出了用结构化稀疏解来分类的算法Robust classification using structured sparse representation。具体的稀疏表达分类算法在以后的博客中会写到，这篇文章关注结构化稀疏解的算法——Block Orthogonal Matching Pursuit (BOMP)，Group Orthogonal Matching Pursuit (GOMP). 和其他凸优化方法相比贪婪组稀疏方法在速度上有很大优势。

这篇博客的matlab实现已上传到github.

BOMP

输入 观测数据向量y∈Rm和字典矩阵A∈Rm×n.
输出 稀疏系数向量 x∈Rn.
Step 1 初始化 令标签集Ω0=∅,初始残差向量r0=y，令k=1.
Step 2 辨识 求矩阵A中与残差向量rk−1最强相关的列，并将该列所在的组的所有列加入标签集Ωk
jk∈argmaxj|<rk−1,aj>|,lk={j|j∈group(jk)},Ωk=Ωk−1∪lk
Step 3 估计 最小化问题 minx∥y−AΩkx∥的解由
xk=(AHΩkAΩk)−1AHΩky给出，其中
AΩk=[aω1,...,aωk],ω1,...,ωk∈Ωk
Step 4 更新残差
rk=y−AΩkxk
Step 5 令k=k+1，并重复Step 2 至Step 4。若某个停止判据满足，则停止迭代。
Step 6 输出系数向量
x(i)=xk(i), if i∈Ωk
x(i)=0, otherwise

注: AH为共轭转置

下面是三种常用的停止判据
1. 运行到某个固定的迭代步数后停止。
2. 残差能量小于某个预先给定值ε。
∥rk∥2≤ε
3. 当字典矩阵A的任何一列都没有残差向量rk的明显能量时
∥AHrk∥∞≤ε

GOMP

GOMP算法和BOMP算法类似，唯一的不同是在Step2，GOMP不求矩阵A中与残差向量rk−1最强相关的列，而是求残差向量rk−1与整个组中的所有原子的相关度平均值，平均值最高的组整体加入标签集Ω. 详见github。

BOMP的matlab实现

function [x, Out] = BOMP(A,y,groups,varargin)% This is a l21 minimization solver for block orthogonal matching pursuit:%   minimiza ||x||_{2,1}%   subject to y = Ax,% -------------------------------------------------------------------------% Author: Beier ZHU%         Tsinghua University% -------------------------------------------------------------------------% % ==========================Required inputs================================% A -- an m x n matrix% y -- an m x 1 vector% group -- an n-entry vector whose i-th entry is the group number of x_i% ==========================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)) ~= 3        error('Wrong number of required parameters');    end    %----------------------------------------------------------------------    % Set parameters to their defaults    %----------------------------------------------------------------------    opts.tol = 1e-6;    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_A] = size(A);    n = length(groups);    [m_y, n_y] = size(y);    if n ~= n_A        error(['group length must equal to the # of columns']);    end    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-6;    while true        [~, max_index] = max(abs(A'*r));        group_id = groups(max_index);        group_index = find(groups==group_id);        Omega = [Omega group_index'];        A_Omega = [A_Omega A(:,group_index')];        x_k = A_Omega\y;        r = y - A_Omega*x_k;        iter = iter + 1;        if (iter > opts.maxIter) || (norm(r) <= opts.tol) || (norm(A'*r, inf) <= varepsilon)                break;        end    end    x(Omega) = x_k;    Out.iter = iter;    Out.residual = norm(r);    Out.energy = norm(A'*r, inf);    % Out.tmp = Omega;end