www.pudn.com > lrr.rar > lrraffine.m, change:2014-12-16,size:2465b

```function [Z,E] = lrraffine(X,A,lambda)
% This routine solves the following nuclear-norm optimization problem.
% Comparing to the algorithm in "lrr.m", here the difference is that
% each column of Z is enforced to sum to one. The goal of this modification
% is to handle the data drawn from affine subspaces. Moreover, this
% modification can improve the performance of LRR even when the
% subspaces are linear.
% --------------------------------
% min |Z|_*+lambda*|E|_2,1
% s.t., X = AZ+E
%       e^tZ = e^t (i.e., each column of Z sums to one)
% inputs:
%        X -- D*N data matrix, D is the data dimension, and N is the number
%             of data vectors.
%        A -- D*M matrix of a dictionary, M is the size of the dictionary
if nargin<3
lambda = 1;
end
tol = 1e-8;
maxIter = 1e6;
[d n] = size(X);
m = size(A,2);
rho = 1.1;
max_mu = 1e30;
mu = 1e-6;
atx = A'*X;
e1 = ones(m,1);
e2 = ones(n,1);
eet = e1*e1';
inv_a = inv(A'*A+eye(m)+eet);
%% Initializing optimization variables
% intialize
J = zeros(m,n);
Z = zeros(m,n);
E = sparse(d,n);

Y1 = zeros(d,n);
Y2 = zeros(m,n);
Y3 = zeros(1,n);
%% Start main loop
iter = 0;
disp(['initial,rank=' num2str(rank(Z))]);
while iter<maxIter
iter = iter + 1;
%update J
temp = Z + Y2/mu;
[U,sigma,V] = svd(temp,'econ');
sigma = diag(sigma);
svp = length(find(sigma>1/mu));
if svp>=1
sigma = sigma(1:svp)-1/mu;
else
svp = 1;
sigma = 0;
end
J = U(:,1:svp)*diag(sigma)*V(:,1:svp)';
%udpate Z
Z = inv_a*(e1*e2'+atx-A'*E+J+(A'*Y1-Y2-e1*Y3)/mu);
%update E
xmaz = X-A*Z;
temp = xmaz+Y1/mu;
E = solve_l1l2(temp,lambda/mu);

leq1 = xmaz-E;
leq2 = Z-J;
stopC = max(max(max(abs(leq1))),max(max(abs(leq2))));
leq3 = e1'*Z-e2';
stopC = max(stopC,max(abs(leq3)));
if iter==1 || mod(iter,50)==0 || stopC<tol
disp(['iter ' num2str(iter) ',mu=' num2str(mu,'%2.1e') ...
',rank=' num2str(rank(Z,1e-3*norm(Z,2))) ',stopALM=' num2str(stopC,'%2.3e')]);
end
if stopC<tol
break;
else
Y1 = Y1 + mu*leq1;
Y2 = Y2 + mu*leq2;
Y3 = Y3 + mu*leq3;
mu = min(max_mu,mu*rho);
end
end

function [E] = solve_l1l2(W,lambda)
n = size(W,2);
E = W;
for i=1:n
E(:,i) = solve_l2(W(:,i),lambda);
end

function [x] = solve_l2(w,lambda)
% min lambda |x|_2 + |x-w|_2^2
nw = norm(w);
if nw>lambda
x = (nw-lambda)*w/nw;
else
x = zeros(length(w),1);
end```