trainLFM.m

function [X,Y] = trainLFM(P, f, alpha, lambda)
    %trainLFM Trains Latent Factor Model for CF with implicit feedback
    %   Takes a playcount matrix P and its parameters alpha and lambda.
    %   and uses the C_{u,s} = 1 + \alpha P_{u,s} as confidence model

    [num_users, num_songs] = size(P);

    % Initialize decomposisiton X and Y with values from 0 to 1
    X = rand(num_users, f);
    Y = rand(num_songs, f);

    % TODO add convergence check.

    % Update the user rows with the least squares expression
    % x_{u} = (Y^{T} C^{u} Y - \lambda I)^{-1} X^{T} C^{u} t_{u}

    % Use this equality to precompute Y^{T} Y for all users
    % Y^{T} C^{u} Y = Y^{T} Y + Y^{T} (C^{u} - I) Y
    A = Y'*Y;

    % Update all users
    parfor u = 1:num_users
        AU = A + (Y' * diag(alpha*P(u,:)) * Y) - lambda * eye(f);
        tu = (P(u,:) > 0)';
        b = Y' * (diag(alpha*P(u,:)) + diag(tu)) * tu;
        X(u,:) = AU\b;
    end

    % Update the song rows with the least squares expression
    % y_{s} = (X^{T} C^{s} X - \lambda I)^{-1} X^{T} C^{s} X^{T} t_{s}
    
    % Use this equality to precompute Y^{T} Y for all users
    % X^{T} C^{u} X = X^{T} X + X^{T} (C^{s} - I) X
    A = X'*X;
    parfor s = 1:num_songs
        AS = A + (X' * diag(alpha * P(:,s)) * X) - lambda * eye(f);
        ts = (P(:,s) > 0);
        b = X' * (diag(alpha * P(:,s)) + diag(ts)) * ts;
        Y(s,:) = AS\b;
    end
end