Contents
function h = matrix_decomp
Problem instance settings
m = 20;
n = 50;
use_cvx = 1;
Problem data
s = RandStream.create('mt19937ar','seed',5489);
RandStream.setDefaultStream(s);
N = 3;
r = 4;
L = randn(m,r) * randn(r,n);
S = sprandn(m,n,0.05);
S(S ~= 0) = 20*binornd(1,0.5,nnz(S),1)-10;
V = 0.01*randn(m,n);
A = S + L + V;
g2_max = norm(A(:),inf);
g3_max = norm(A);
g2 = 0.15*g2_max;
g3 = 0.15*g3_max;
CVX
if use_cvx
tic;
cvx_begin
cvx_precision low
variables X_1(m,n) X_2(m,n) X_3(m,n)
minimize(0.5*square_pos(norm(X_1,'fro')) + g2*norm(X_2(:),1) + g3*norm_nuc(X_3))
subject to
X_1 + X_2 + X_3 == A;
cvx_end
h.cvx_toc = toc;
h.p_cvx = cvx_optval;
h.X1_cvx = X_1;
h.X2_cvx = X_2;
h.X3_cvx = X_3;
X_2(abs(X_2) < 1e-4) = 0;
rhat = sum(svd(X_3) > 1e-4);
fprintf('CVX (vs true):\n');
fprintf('|V| = %.2f; |X_1| = %.2f\n', norm(V, 'fro'), norm(X_1,'fro'));
fprintf('nnz(S) = %d; nnz(X_2) = %d\n', nnz(S), nnz(X_2));
fprintf('rank(L) = %d; rank(X_3) = %d\n', rank(L), rhat);
end
Calling sedumi: 5491 variables, 1002 equality constraints
------------------------------------------------------------
SeDuMi 1.21 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
eqs m = 1002, order n = 2077, dim = 7908, blocks = 1004
nnz(A) = 3005 + 0, nnz(ADA) = 1002004, nnz(L) = 501503
it : b*y gap delta rate t/tP* t/tD* feas cg cg prec
0 : 1.48E+04 0.000
1 : 7.70E+02 1.31E+04 0.000 0.8795 0.9000 0.9000 4.73 1 1 7.6E+00
2 : 1.47E+03 1.12E+04 0.000 0.8559 0.9000 0.9000 3.65 1 1 5.5E+00
3 : 1.78E+03 8.88E+03 0.000 0.7945 0.9000 0.9000 2.86 1 1 3.8E+00
4 : 2.01E+03 6.86E+03 0.000 0.7724 0.9000 0.9000 2.26 1 1 2.6E+00
5 : 1.99E+03 5.88E+03 0.000 0.8575 0.9000 0.9000 1.84 1 1 2.2E+00
6 : 2.24E+03 4.05E+03 0.000 0.6891 0.9000 0.9000 1.65 1 1 1.4E+00
7 : 2.18E+03 2.88E+03 0.000 0.7119 0.9000 0.9000 1.41 1 1 1.0E+00
8 : 2.20E+03 1.65E+03 0.000 0.5726 0.9000 0.9000 1.27 1 1 6.0E-01
9 : 2.18E+03 1.01E+03 0.000 0.6124 0.9000 0.9000 1.13 1 1 3.9E-01
10 : 2.14E+03 3.02E+02 0.000 0.2985 0.9000 0.9000 1.08 1 1 1.3E-01
11 : 2.13E+03 7.40E+01 0.000 0.2452 0.9000 0.9000 1.02 1 1 3.3E-02
12 : 2.13E+03 1.76E+01 0.000 0.2372 0.9036 0.9000 1.01 1 1 8.1E-03
13 : 2.12E+03 3.51E+00 0.000 0.1998 0.9086 0.9000 1.00 1 1 1.8E-03
14 : 2.12E+03 6.93E-01 0.000 0.1977 0.9084 0.9000 1.00 1 1 4.2E-04
15 : 2.12E+03 1.60E-01 0.000 0.2314 0.9049 0.9000 1.00 1 1 1.1E-04
16 : 2.12E+03 8.96E-03 0.000 0.0558 0.9000 0.0000 1.00 1 1 4.8E-05
17 : 2.12E+03 5.97E-09 0.000 0.0000 0.9134 0.9000 1.00 1 1 1.3E-05
18 : 2.12E+03 1.25E-09 0.000 0.2092 0.9000 0.9000 1.00 1 1 2.7E-06
19 : 2.12E+03 3.02E-10 0.000 0.2419 0.9000 0.9000 1.00 1 1 6.4E-07
iter seconds digits c*x b*y
19 7.6 Inf 2.1233920957e+03 2.1233920997e+03
|Ax-b| = 1.6e-06, [Ay-c]_+ = 1.1E-06, |x|= 7.0e+02, |y|= 3.5e+02
Detailed timing (sec)
Pre IPM Post
3.000E-02 7.610E+00 2.000E-02
Max-norms: ||b||=1.453096e+01, ||c|| = 4.948966e+00,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 153.012.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +2123.39
CVX (vs true):
|V| = 0.31; |X_1| = 26.24
nnz(S) = 49; nnz(X_2) = 53
rank(L) = 4; rank(X_3) = 4
ADMM
MAX_ITER = 100;
ABSTOL = 1e-4;
RELTOL = 1e-2;
tic;
lambda = 1;
rho = 1/lambda;
X_1 = zeros(m,n);
X_2 = zeros(m,n);
X_3 = zeros(m,n);
z = zeros(m,N*n);
U = zeros(m,n);
fprintf('\n%3s\t%10s\t%10s\t%10s\t%10s\t%10s\n', 'iter', ...
'r norm', 'eps pri', 's norm', 'eps dual', 'objective');
for k = 1:MAX_ITER
B = avg(X_1, X_2, X_3) - A./N + U;
X_1 = (1/(1+lambda))*(X_1 - B);
X_2 = prox_l1(X_2 - B, lambda*g2);
X_3 = prox_matrix(X_3 - B, lambda*g3, @prox_l1);
x = [X_1 X_2 X_3];
zold = z;
z = x + repmat(-avg(X_1, X_2, X_3) + A./N, 1, N);
U = B;
h.objval(k) = objective(X_1, g2, X_2, g3, X_3);
h.r_norm(k) = norm(x - z,'fro');
h.s_norm(k) = norm(-rho*(z - zold),'fro');
h.eps_pri(k) = sqrt(m*n*N)*ABSTOL + RELTOL*max(norm(x,'fro'), norm(-z,'fro'));
h.eps_dual(k) = sqrt(m*n*N)*ABSTOL + RELTOL*sqrt(N)*norm(rho*U,'fro');
if k == 1 || mod(k,10) == 0
fprintf('%4d\t%10.4f\t%10.4f\t%10.4f\t%10.4f\t%10.2f\n', k, ...
h.r_norm(k), h.eps_pri(k), h.s_norm(k), h.eps_dual(k), h.objval(k));
end
if h.r_norm(k) < h.eps_pri(k) && h.s_norm(k) < h.eps_dual(k)
break;
end
end
h.admm_toc = toc;
h.admm_iter = k;
h.X1_admm = X_1;
h.X2_admm = X_2;
h.X3_admm = X_3;
fprintf('\nADMM (vs true):\n');
fprintf('|V| = %.2f; |X_1| = %.2f\n', norm(V, 'fro'), norm(X_1,'fro'));
fprintf('nnz(S) = %d; nnz(X_2) = %d\n', nnz(S), nnz(X_2));
fprintf('rank(L) = %d; rank(X_3) = %d\n', rank(L), rank(X_3));
if use_cvx
fprintf('\nADMM vs CVX solutions (in Frobenius norm):\n');
fprintf('X_1: %.2e; X_2: %.2e; X_3: %.2e\n', ...
norm(h.X1_cvx - X_1,'fro'), norm(h.X2_cvx - X_2,'fro'), norm(h.X3_cvx - X_3,'fro'));
end
iter r norm eps pri s norm eps dual objective
1 41.7766 0.6214 61.5891 0.6077 584.62
10 9.7749 0.7995 7.9580 0.6145 2538.44
20 3.9687 0.8679 2.3887 0.4883 2598.00
30 1.2071 0.8445 0.6616 0.4595 2477.69
ADMM (vs true):
|V| = 0.31; |X_1| = 26.18
nnz(S) = 49; nnz(X_2) = 52
rank(L) = 4; rank(X_3) = 4
ADMM vs CVX solutions (in Frobenius norm):
X_1: 3.59e-01; X_2: 6.14e-01; X_3: 5.30e-01
end
function x = avg(varargin)
N = length(varargin);
x = 0;
for k = 1:N
x = x + varargin{k};
end
x = x/N;
end
function p = objective(X_1, g_2, X_2, g_3, X_3)
p = norm(X_1,'fro').^2 + g_2*norm(X_2(:),1) + g_3*norm(svd(X_3),1);
end