Contents

% Distributed L1 regularized logistic regression
% (compared against l1_logreg package)

Generate problem data

rand('seed', 0);
randn('seed', 0);

n = 200;
m = 200;
N = 100;

w = sprandn(n, 1, 100/n);       % N(0,1), 10% sparse
v = randn(1);                  % random intercept

X0 = sprandn(m*N, n, 10/n);           % data / observations
btrue = sign(X0*w + v);

% noise is function of problem size use 0.1 for large problem
b0 = sign(X0*w + v + sqrt(0.1)*randn(m*N, 1)); % labels with noise

% packs all observations in to an m*N x n matrix
A0 = spdiags(b0, 0, m*N, m*N) * X0;

ratio = sum(b0 == 1)/(m*N);
mu = 0.1*1/(m*N) * norm((1-ratio)*sum(A0(b0==1,:),1) + ratio*sum(A0(b0==-1,:),1), 'inf');

x_true = [v; w];

Solve problem

[x history] = distr_l1_logreg(A0, b0, mu, N, 1.0, 1.0);
iter	    # bfgs	    r norm	   eps pri	    s norm	  eps dual	 objective
  1	      2034	   57.0681	    0.6804	   19.7440	    0.5849	  11824.14
  2	      1833	   26.5753	    0.5584	   28.1814	    0.8198	   9973.21
  3	      1773	   16.2131	    0.7186	   21.4621	    0.9342	   9219.04
  4	      1696	   11.9694	    0.8572	   15.0765	    1.0017	   8914.00
  5	      1661	    9.2421	    0.9576	   10.5568	    1.0475	   8782.13
  6	      1560	    7.2694	    1.0297	    7.5041	    1.0813	   8720.73
  7	      1513	    5.8004	    1.0819	    5.4344	    1.1076	   8690.23
  8	      1486	    4.6913	    1.1200	    3.9979	    1.1288	   8674.40
  9	      1424	    3.8447	    1.1482	    2.9738	    1.1462	   8665.95
 10	      1397	    3.1897	    1.1691	    2.2269	    1.1608	   8661.38
 11	      1359	    2.6766	    1.1846	    1.6725	    1.1730	   8658.87
 12	      1305	    2.2700	    1.1963	    1.2562	    1.1835	   8657.51
 13	      1255	    1.9447	    1.2050	    0.9416	    1.1924	   8656.76
 14	      1178	    1.6828	    1.2114	    0.7029	    1.2001	   8656.36
 15	      1100	    1.4707	    1.2162	    0.5215	    1.2068	   8656.14
 16	      1030	    1.2984	    1.2198	    0.3843	    1.2127	   8656.02
 17	       968	    1.1578	    1.2223	    0.2810	    1.2179	   8655.96
Elapsed time is 6.819467 seconds.

Reporting

K = length(history.objval);

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)');