Hi.

I am stuck at a similar problem in cvx to solve for Laplacain matrix L in

maximize {trace(S*X) - 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(S*X) - 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_matrix*X) - lambda*norm(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 ?