Hi.
I am stuck at a similar problem in cvx to solve for Laplacain matrix L in
maximize {trace(SX) - log det X + lambda||X||_1}
where S = covariance matrix
X = precision matrix
||X||_1 = 1-norm of X defined as: norm(L,1) -(1/variance)*N)
N = no. of points in the data matrix
The problem is to be solved using ADMM algorithm in cvx.
I used the following code:
[X, history] = covsel(Covariance_matrix, lambda, 1, 1);
K = length(history.objval);
X_admm = X;
h = figure;
plot(1:K, history.objval, ‘k’, ‘MarkerSize’, 10, ‘LineWidth’, 2);
ylabel(‘f(x^k) + g(z^k)’); xlabel(‘iter (k)’);
g = figure;
subplot(2,1,1);
semilogy(1:K, max(1e-8, history.r_norm), ‘k’, …
1:K, history.eps_pri, ‘k–’, ‘LineWidth’, 2);
ylabel(’||r||_2’);
subplot(2,1,2);
semilogy(1:K, max(1e-8, history.s_norm), ‘k’, …
1:K, history.eps_dual, ‘k–’, ‘LineWidth’, 2);
ylabel(’||s||_2’); xlabel(‘iter (k)’);
function [Z, history] = covsel(Covariance_matrix, lambda, rho, alpha)
% covsel Sparse inverse covariance selection via ADMM
%
% [X, history] = covsel(D, lambda, rho, alpha)
%
% Solves the following problem via ADMM:
%
% minimize trace(SX) - log det X + lambda||X||_1
%
% with variable X, where S is the empirical covariance of the data
% matrix D (training observations by features).
%
% The solution is returned in the matrix X.
%
% history is a structure that contains the objective value, the primal and
% dual residual norms, and the tolerances for the primal and dual residual
% norms at each iteration.
%
% rho is the augmented Lagrangian parameter.
%
% alpha is the over-relaxation parameter (typical values for alpha are
% between 1.0 and 1.8).
%
t_start = tic;
% Global constants and defaults
QUIET = 0;
MAX_ITER = 1000;
ABSTOL = 1e-4;
RELTOL = 1e-2;
S = Covariance_matrix;
N = size(S,1);
X = Precision_matrix;
% ADMM solver
X = zeros(N);
Z = zeros(N);
U = zeros(N);
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,L] = eig(rho*(Z - U) - S);
es = diag(L);
xi = (es + sqrt(es.^2 + 4*rho))./(2*rho);
X = Q*diag(xi)*Q';
% z-update with relaxation
Z_old = Z;
X_hat = alpha*X + (1 - alpha)*Z_old;
Z = shrinkage(X_hat + U, lambda/rho);
U = U + (X_hat - Z);
% diagnostics, reporting, termination checks
history.objval(k) = objective( X, S, Z, lambda);
history.r_norm(k) = norm(X - Z, 'fro');
history.s_norm(k) = norm(-rho*(Z - Z_old),'fro');
history.eps_pri(k) = sqrt(N*N)*ABSTOL + RELTOL*max(norm(X,'fro'), norm(Z,'fro'));
history.eps_dual(k)= sqrt(N*N)*ABSTOL + RELTOL*norm(rho*U,'fro');
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
end
function obj = objective(X, Covariance_matrix, Z, lambda)
obj = log(det(X)) - trace(Covariance_matrixX) - lambdanorm(Z(:), 1);
end
function y = shrinkage(a, kappa)
y = max(0, a-kappa) - max(0, -a-kappa);
end
Use this from :https://web.stanford.edu/~boyd/papers/admm/covsel/covsel.html
But there is an error in : history.objval(k) = objective( X, S, Z, lambda)
and in [X, history] = covsel(Covariance_matrix, lambda, 1, 1) ,
stated below:
“Unable to perform assignment because the indices on the left side are not compatible with the size of the right side”.
Any suggestions as to where the mistake in the code is ?