## Contents

function [xave, history] = linear_svm(A, lambda, p, rho, alpha)

% linear_svm   Solve linear support vector machine (SVM) via ADMM
%
% [x, history] = linear_svm(A, lambda, p, rho, alpha)
%
% Solves the following problem via ADMM:
%
%   minimize   (1/2)||w||_2^2 + \lambda sum h_j(w, b)
%
% where A is a matrix given by [-y_j*x_j -y_j], lambda is a
% regularization parameter, and p is a partition of the observations in to
% different subsystems.
%
% The function h_j(w, b) is a hinge loss on the variables w and b.
% It corresponds to h_j(w,b) = (Ax + 1)_+, where x = (w,b).
%
% This function implements a *distributed* SVM that runs its updates
% serially.
%
% The solution is returned in the vector x = (w,b).
%
% 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;


## Data preprocessing

[m, n] = size(A);
N = max(p);
% group samples together
for i = 1:N,
tmp{i} = A(p==i,:);
end
A = tmp;


x = zeros(n,N);
z = zeros(n,N);
u = zeros(n,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
for i = 1:N,
cvx_begin quiet
variable x_var(n)
minimize ( sum(pos(A{i}*x_var + 1)) + rho/2*sum_square(x_var - z(:,i) + u(:,i)) )
cvx_end
x(:,i) = x_var;
end
xave = mean(x,2);

% z-update with relaxation
zold = z;
x_hat = alpha*x +(1-alpha)*zold;
z = N*rho/(1/lambda + N*rho)*mean( x_hat + u, 2 );
z = z*ones(1,N);

% u-update
u = u + (x_hat - z);

% diagnostics, reporting, termination checks
history.objval(k)  = objective(A, lambda, p, x, z);

history.r_norm(k)  = norm(x - z);
history.s_norm(k)  = norm(-rho*(z - zold));

history.eps_pri(k) = sqrt(n)*ABSTOL + RELTOL*max(norm(x), norm(-z));
history.eps_dual(k)= sqrt(n)*ABSTOL + RELTOL*norm(rho*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

end

function obj = objective(A, lambda, p, x, z)
obj = hinge_loss(A,x) + 1/(2*lambda)*sum_square(z(:,1));
end

function val = hinge_loss(A,x)
val = 0;
for i = 1:length(A)
val = val + sum(pos(A{i}*x(:,i) + 1));
end
end