本帖最后由 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 上传

运行结果