function kfda(K, y, regLambda)
% compute the kernel matrix

% index sets of two classes
I_1 = find(y==1);
I_2 = find(y==-1);
m_1 = length(I_1);
m_2 = length(I_2);
m = m_1+m_2;

K_1 = K(:, I_1);
K_2 = K(:, I_2);

% kernel means
M_1 = mean(K_1, 2);
M_2 = mean(K_2, 2);

M = (M_1-M_2) * (M_1-M_2)'; % between class variance in the kernel space

J_1 = (1/sqrt(m_1)) * (eye(m_1)-(1/m_1));
J_2 = (1/sqrt(m_2)) * (eye(m_2)-(1/m_2));

N = K_1 * J_1 * J_1 * K_1 + K_2 * J_2 * J_2 * K_2 + regLambda*eye(m_1+m_2);

% coefficients for the optimal projection direction w
a = [(1/m_1) * ones(m_1,1); -(1/m_2) * ones(m_2, 1)];

J = [J_1, zeros(m_1,m_2); zeros(m_2,m_1), J_2];
Inv = inv(regLambda * eye(m) + J * K * J);
alpha = (1/regLambda) * (eye(m) - J * Inv * J * K) * a;

ratio = (alpha' * M * alpha) / (alpha' * M * alpha);