【机器学习】LatLRR

LatLRR

参考论文：Latent Low-Rank Representation for Subspace Segmentation and Feature Extraction
作者：Guangcan Liu， Shuicheng Yan
Electrical and Computer Engineering, National University of Singapore

LRR

 LRR的主要思想是获取给定字典的样本的低秩表示，LRR问题可以描述为：$$mi{n}_z||Z|{|}_{\ast },s.t.X_o=AZ$$其中$X_o$表示观测到的数据矩阵（每一列是一个观测样本），$A$是字典。对于子空间分割，通常用观测到的数据本身作为字典，所以问题变为：$$mi{n}_z||Z|{|}_{\ast },s.t.X_o=X_{o}Z①$$然而，数据常被多种噪声干扰，为增强鲁棒性，上式可改为：$$mi{n}_z||Z|{|}_{\ast }+\lambda ||E|{|}_l,s.t.X_o=X_{o}Z+E$$但是，LRR由于直接把观测到的数据本身作为字典，所以会有一些问题：

• 首先，为了能够表示底层子空间，字典（$A=X_O$）必须包含从子空间中采样的足够的数据样本，否则，$Z=I$可能是①的唯一可行解，因此LRR失败。所以，如果数据采样不够，就不能用$X_O$作为字典来表示子空间；
• 其次，为了实现鲁棒的分割，LRR需要字典$A$中有足够的无噪声数据，当$A=X_O$时，这个假设可能是无效的，这会降低LRR的鲁棒性。

也可参见自本人的另一篇归纳文章： 【机器学习】LPP\NPE\SPP\CRP\RPCA\LRR\LRPP\LRPE

LatLRR

所以，为了解决上述两个问题，考虑以下模型（LatLRR）：$$mi{n}_z||Z|{|}_{\ast }+\lambda ||E|{|}_l,s.t.X_o=[X_o,X_H]Z+E$$其中$X_O$是观测到的数据，$X_H$是隐藏的数据，

令$L_{H|O}^{\ast }=U\Sigma V_H^T{V}_H{\Sigma }^{-1}{U}^T$，则：

假设$X_O和X_H$取样于同一个低秩子空间集合，并且子空间集合的秩为$r$，然后，

所以，$Z_{O|H}^{\ast }和L_{H|O}^{\ast }$都是低秩的，我们可以通过下面式子来恢复$Z_{O|H}^{\ast }$

所以：

$L$是投影矩阵：

注意：$L\in R^{d\ast d},X\in R^{d\ast n}$,所以$Y$和$X$是同样维数的。

算法如下：

代码：

%% Latent Low-Rank Representation
function [Z,L,E] = latent_lrr(X,lambda)
% Latent Low-Rank Representation for Subspace Segmentation and Feature Extraction
% Guangcan Liu, Shuicheng Yan. ICCV 2011.
% Problem:
%   min_Z,L,E ||Z||_* + ||L||_* +¡¡lambda||E||_1,
%           s.t. X = XZ + LX +E.
% Solning problem by Inexact ALM

A = X;
tol = 1e-6;
rho = 1.1;
max_mu = 10;
mu = 1e-6;
maxIter = 1e6;
[d n] = size(X);
m = size(A,2);
atx = X'*X;
inv_a = inv(A'*A+eye(m));
inv_b = inv(A*A'+eye(d));
%% Initializing optimization variables
J = zeros(m,n);
Z = zeros(m,n);
L = zeros(d,d);
S = zeros(d,d);

% E = sparse(d,n);
E = zeros(d,n);

Y1 = zeros(d,n);
Y2 = zeros(m,n);
Y3 = zeros(d,d);

%% Start main loop
iter = 0;
disp('initial');

while iter<maxIter
    iter = iter + 1;
%     disp(['iter====>>>>' num2str(iter)]);
    %updating J by the Singular Value Thresholding(SVT) operator
    temp_J = Z + Y2/mu;
    [U_J,sigma_J,V_J] = svd(temp_J,'econ');
    sigma_J = diag(sigma_J);
    svp_J = length(find(sigma_J>1/mu));
    if svp_J>=1
        sigma_J = sigma_J(1:svp_J)-1/mu;
    else
        svp_J = 1;
        sigma_J = 0;
    end
    J = U_J(:,1:svp_J)*diag(sigma_J)*V_J(:,1:svp_J)';
    
    %updating S by the Singular Value Thresholding(SVT) operator
    temp_S = L + Y3/mu;
    [U_S,sigma_S,V_S] = svd(temp_S,'econ');
    sigma_S = diag(sigma_S);
    svp_S = length(find(sigma_S>1/mu));
    if svp_S>=1
        sigma_S = sigma_S(1:svp_S)-1/mu;
    else
        svp_S = 1;
        sigma_S = 0;
    end
    S = U_S(:,1:svp_S)*diag(sigma_S)*V_S(:,1:svp_S)';
    
    %udpate Z
    Z = inv_a*(atx-X'*L*X-X'*E+J+(X'*Y1-Y2)/mu);
    
    %udpate L
    L = ((X-X*Z-E)*X'+S+(Y1*X'-Y3)/mu)*inv_b;
    
    %update E
    xmaz = X-X*Z-L*X;
    temp = xmaz+Y1/mu;
    E = max(0,temp - lambda/mu)+min(0,temp + lambda/mu);
    
    leq1 = xmaz-E;
    leq2 = Z-J;
    leq3 = L-S;
    max_l1 = max(max(abs(leq1)));
    max_l2 = max(max(abs(leq2)));
    max_l3 = max(max(abs(leq3)));
    
    stopC1 = max(max_l1, max_l2);
    stopC = max(stopC1, max_l3);
    if stopC<tol
        disp('LRR done.');
        break;
    else
        Y1 = Y1 + mu*leq1;
        Y2 = Y2 + mu*leq2;
        Y3 = Y3 + mu*leq3;
        mu = min(max_mu,mu*rho);
    end
    
end
end