% Lp regularised multiple kernel Fisher discriminant analysis via 
% semi-infinite programming (SIP). Multiclass version.
% The QCLP problem in the beta step is solved using the Mosek toolbox.
% Fei Yan, f.yan@surrey.ac.uk, University of Surrey, UK, October 2009.

% Inputs
% "p": norm for regularising the kernel weights. must be greater than 1
% "G": the cell containing centered kernel matrices
% "label": labels of the training samples
% "lambda": regularization parameter
% "epsilon": tolerance parameter
% "max_iter": maximal number of iterations

% Outputs
% "alpha": the learnt coefficients of training samples
% "beta": the learnt weights of kernels
% "theta": value of the objective function in SIP
% "violation": normalised constraint violation

% References: 
% [1] Multiclass Discriminant Kernel Learning via Convex Programming. JMLR 2008
% [2] Lp Norm Multiple Kernel Fisher Discriminant Analysis for Object and Image Categorisationss. 
%     CVPR 2010


function [alpha beta theta violation] = lp_mk_fda_multiclass(p, G, label, lambda, epsilon, max_iter)


m = length(label);
n = length(G);

unique_label = unique(label); % labels in "unique_label" are sorted in ascending order
c = length(unique_label);
for k = 1:c
	mm(k) = length(find(label==unique_label(k)));
end

h = zeros(m, c);
for i = 1:m
	for k = 1:c
		if(label(i)==unique_label(k))
			h(i, k) = sqrt(m/mm(k)) - sqrt(mm(k)/m);
		else
			h(i, k) = -sqrt(mm(k)/m);
		end
	end
end

S = 1;
theta = -inf;
beta = ones(n, 1);
beta = beta/((sum(beta.^p))^(1/p));
violation = zeros(max_iter, 1);

prob.c = [zeros(n, 1); -1];

prob.a = sparse(eye(n+1)); % prob.a must be sparse
prob.a(n+1, n+1) = 0;
prob.blc = [zeros(n, 1); -inf];	% the last one has to be -inf, not 0. becuase 	
				% quadratic constrait cannot have finite lower 
				% and upper bound at the same time
prob.buc = ones(n+1, 1);

			
t = 0;
while(t < max_iter)
	
	t = t+1;

	% alpha step

	G_Sum = zeros(m, m);
	for j = 1:n
		G_Sum = G_Sum + beta(j, 1) * G{j};
	end
	b = 1/2 * eye(m) + 1/(2*lambda) * G_Sum;
	alpha = b \ h; % alpha is m x c
	
	u = zeros(n, 1);
	for j = 1:n
		for k = 1:c
			u(j, 1) = u(j, 1) + 1/(4*lambda) * (alpha(:, k)'*G{j}*alpha(:, k));
		end
	end
	v = 0;
	for k = 1:c
		v = v + 1/4*alpha(:, k)'*alpha(:, k) - alpha(:, k)'*h(:, k);
	end
	S = beta' * u + v;

	violation(t, 1) = abs(1-S/theta);
	if(violation(t, 1) <= epsilon)
		break;
	end


	% beta step
	
	prob.a = [prob.a; [u' -1]];
	prob.blc = [prob.blc; -v];
	prob.buc = [prob.buc; +inf];
		
	% modify the {n+1}^{th} constraint, which is the norm constraint
	
	prob.qcsubk = [(n+1)*ones(1, n)];
	prob.qcsubi = [1:n];
	prob.qcsubj = [1:n];
	prob.qcval = beta'.^(p-2) * p * (p-1);
	
	prob.a(n+1, :) = [- beta'.^(p-1) * p * (p-2) 0];
	prob.blc(n+1, 1) = -inf;
	prob.buc(n+1, 1) = - p * (p-3) / 2;
	
			
	[r, res] = mosekopt('minimize echo(0)', prob);
	beta = res.sol.itr.xx(1:n);
	theta = res.sol.itr.xx(n+1);
end
