function K = correlation(Y)

fast = 1;
if fast
    m     = size(Y, 1);
    Yi    = isnan(Y);
    Y(Yi) = 0;
    known = m - sum(Yi, 1);

    % mean
    mu = sum(Y, 1)./max(known, 1);

    % standard deviation
    H = (Y - repmat(mu, m, 1)).^2;
    H(Yi) = 0;
    sigma = sqrt(sum(H, 1)./max(known - 1, 1));
    sigma(sigma < eps) = eps;

    % z-normalization
    Y = (Y - repmat(mu, m, 1))./repmat(sigma, m, 1);
    Y(Yi) = 0;

    K = (Y'*Y)./max((1 - Yi)'*(1 - Yi) - 1, 1);
else
    [m,n] = size(Y);

    for i = 1:n
        % known ratings
        k = find(isnan(Y(:,i)) == 0); %find(~isnan(Y(:,i)))    
        Y(k,i) = (Y(k,i) - mean(Y(k,i)))/max(std(Y(k,i)), eps);    

        % (normalized) 0-ratings
        z = find(abs(Y(:,i)) < eps);
        Y(z,i) = eps;

        % unknown ratings
        u = find(isnan(Y(:,i)) == 1);
        Y(u,i) = 0;
    end

    K = zeros(n, n);
    for i = 1:n
        Z = Y .* repmat(Y(:,i), 1, n);
        K(i,:) = sum(Z)./max(sum(Z ~= 0) - 1, eps);
    end
end
