本帖最后由 zhangxueqing 于 2017-3-26 16:36 编辑
谢谢!但是如果不放循环里,那么图形画不出来?
t_start= tic;
%Global constants and defaults
QUIET = 0;
MAX_ITER = 1000;
ABSTOL = 1e-5;
RELTOL = 1e-3;
%Data preprocessing
[m, n] = size(A);
fprintf('(1) GADMM1 for lasso problem \n');
% save a matrix-vector multiply
Atb = A'*b;
%ADMM solver
x = zeros(n,1);
z = zeros(n,1);
u = zeros(n,1);
I=eye(n);
% cache the factorization
[L U] = factor(A, rho1);
if ~QUIET
fprintf('%3s\t%10s\t%10s\t%10s\t%10s\t%10s\n', 'iter', ...
'r norm', 'eps pri', 's norm', 'eps dual', 'objective');
end
for k = 1:MAX_ITER
% x-update
q = Atb + rho1*(z + u); % temporary value
if( m >= n ) % if skinny
x = U \ (L \ q);
else % if fat
x = q/rho1 - (A'*(U \ ( L \ (A*q) )))/rho1^2;
end
% z-update with relaxation
zold = z;
z = shrinkage(alpha1*x+(1-alpha1)*zold-u, tau/rho1);
% u-update
u = u - (alpha1*x+(1-alpha1)*zold - z);
% diagnostics, reporting, termination checks
history.objval(k) = objective(A, b, tau, x, z);
history.r_norm(k) = norm(x - z);
history.s_norm(k) = norm(-rho1*(z - zold));
history.eps_pri(k) = sqrt(n)*ABSTOL + RELTOL*max(norm(x), norm(-z));
history.eps_dual(k)= sqrt(n)*ABSTOL + RELTOL*norm(rho1*u);
if ~QUIET
fprintf('%3d\t%10.4f\t%10.4f\t%10.4f\t%10.4f\t%10.2f\n', k, ...
history.r_norm(k), history.eps_pri(k), ...
history.s_norm(k), history.eps_dual(k), history.objval(k));
end
if (history.r_norm(k) < history.eps_pri(k) && ...
history.s_norm(k) < history.eps_dual(k))
break;
end
end
if ~QUIET
toc(t_start);
end
h1=k;
y1=1/2*norm((A*x - b).^2) + tau*norm(z,1) ; %每步目标函数值
hold on
subplot(2,2,1)
H1=plot(h1,y1,'--.r','linewidth',2)
xlabel('Iteration')
ylabel('The objective function value (alpha=0.5)')
set(H1,'linestyle','-','color','b');
end
function p = objective(A, b, tau, x, z)
p = ( 1/2*norm((A*x - b).^2) + tau*norm(z,1) );
end
function z = shrinkage(x, kappa)
z = max( 0, x - kappa ) - max( 0, -x - kappa );
end
function [L U] = factor(A, rho1)
[m, n] = size(A);
if ( m >= n ) % if skinny
L = chol( A'*A + rho1*speye(n), 'lower' );
else % if fat
L = chol( speye(m) + 1/rho1*(A*A'), 'lower' );
end
% force matlab to recognize the upper / lower triangular structure
L = sparse(L);
U = sparse(L');
end
2017-3-26 16:36 上传
运行结果