function P = mddm(X,Y,d)
%MDDM Multilabel Dimensionality reduction via Dependence Maximization
%   X: DxN feature matrix
%   Y: QxN label matrix
%   d: the dimensionality to be reduced to

MAX_ROWS = 1000; 

D = size(X,1);
N = size(X,2);

L = Y'*Y;
H = eye(N)-1/N; % <= 4000

tic; S = H*L*H; toc; % <= 3500
tic; A = X*S; toc; % <= 1300
tic; A = prodmat(A,X'); toc;

tic;
[P,~,~] = svds(A,d);
toc;

end

