function [ X_pred ] = svd_bias( X )
m = size(X, 1); %users
n = size(X, 2); %items

%define hyperparameters of the model
MAX_ITERS = 100000;
lambda = 0;
f = 1;

[rows, cols, vals] = find(X);
R = [rows, cols, vals];
k = 1;
mu = mean(vals);

bu = k*rand(m, 1);
bi = k*rand(n, 1);
P = k*rand(m, f); % latent vectors for the users
Q = k*rand(n, f); % latent vectors for the items

% bu = [0, 0 ,0]';
% bi = -repmat(mu, 3, 1)-bu;
% disp(bi);
% P = [1 2 3]';
% Q = [1 2 3]';

fprintf('Performing SGD...\n');
for iter=1:MAX_ITERS
    if rem(iter, 1000) == 0
        fprintf('iteration %d\n', iter);
    end
    gamma = 0.05;%1/sqrt(iter);
    %i = i-th observed data point
    i = randi(size(R,1));
    ui = R(i,1);
    ii = R(i,2);
    vi = R(i,3);
    Pui = P(ui,:);
    Qii = Q(ii,:);
    eui = vi - mu - bu(ui) - bi(ii) - dot(Pui,Qii);
    bu(ui) = bu(ui) + gamma*(eui - lambda*bu(ui));
    bi(ii) = bi(ii) + gamma*(eui - lambda*bi(ii));
    Q(ii,:)  = Qii  + gamma*(eui*Pui - lambda*Qii);
    P(ui,:)  = Pui  + gamma*(eui*Qii - lambda*Pui);
end

% disp(bu');
% disp(bi');
% disp(P');
% disp(Q');


X_pred = X;
[zrows, zcols] = find(~X);
I = find(~X);
PQ = P*Q';
disp(repmat(mu, length(zrows), 1) + bu(zrows) + bi(zcols));
X_pred(I) = repmat(mu, length(zrows), 1) + bu(zrows) + bi(zcols) + PQ(I);
%predict all values
%X_pred = repmat(mu, size(PQ)) + repmat(bu,1,n) + repmat(bi',m,1) + PQ;
% disp(mu);
% disp(bu);
% disp(bi);
A = [1 2 3; 2 4 6; 3 6 9];
% disp(X_pred(I)-);
disp(sum((X_pred(I)-A(I)) .^2));

end