diff --git a/others/VBLRMat/Readme.txt b/others/VBLRMat/Readme.txt
new file mode 100644
index 0000000..9b1378a
--- /dev/null
+++ b/others/VBLRMat/Readme.txt
@@ -0,0 +1,26 @@
+These MATLAB codes implement the variational Bayesian low-rank
+matrix completion and robust PCA methods described in the paper
+
+S. D. Babacan, M. Luessi, R. Molina, and A. K. Katsaggelos,
+"Sparse Bayesian Methods for Low-Rank Matrix Estimation,"
+IEEE Transactions on Signal Processing, 2012.
+
+VBMC: Variational Bayesian low-rank matrix completion main code
+VBRPCA: Variational Bayesian robust PCA main code
+
+Example runs of these codes are shown in the scripts RunVBMC and RunVBRPCA. These scripts explain the general use of these codes, additional details can be found in the main codes.
+
+Please contact S. Derin Babacan at dbabacan [at] gmail.com for bug reports/ comments/ suggestions/ questions.
+
+Copyright (c): S. Derin Babacan, 2011
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+
diff --git a/others/VBLRMat/RunVBMC.m b/others/VBLRMat/RunVBMC.m
new file mode 100644
index 0000000..366d08f
--- /dev/null
+++ b/others/VBLRMat/RunVBMC.m
@@ -0,0 +1,74 @@
+%% An example run of Variational Matrix Low-Rank Matrix Completion
+clear
+clc
+
+%% Create a low-rank matrix
+% Dimensions
+m = 300;
+n = 300;
+r = 10; % rank
+
+A = randn(m,r);
+B = randn(r,n);
+X_true = A*B; % low rank
+
+% random sparse sampling
+df = r*(m+n-r); % Degrees of freedom
+p = 0.2; % percentage measurements
+pmn = p*m*n;
+osdf = pmn / df; % Oversampling degrees of freedom
+
+Omega = randperm(m*n); Omega = Omega(1:pmn);
+P = zeros(m,n); P(Omega) = 1;
+
+% Observations
+Y = P.*X_true;
+
+% Add noise?
+sigma = 0; % no noise
+sigma = 1e-3; % some noise
+% sigma = sqrt(1e-3); % more noise
+
+Y = Y + sigma*randn(size(Y));
+
+SNR = 10*log10( sum(sum((P.*X_true).^2)) / sum(sum((Y-P.*X_true).^2)));
+
+%% Run VBMC
+% all options are *optional*, everything will be set automatically
+% you can modify these options to get better performance
+options.verbose = 1;
+options.MAXITER = 100;
+options.DIMRED = 1; % Reduce dimensionality during iterations?
+% you can also set the threshold to reduce dimensions
+% options.DIMRED_THR = 1e3;
+%Estimate noise variance? (beta is inverse noise variance)
+options.UPDATE_BETA = 1;
+% options.UPDATE_BETA_START = 1;% iteration number to start estimating noise variance
+
+% If the noise inv. variance is not to be estimated, set
+% options.UPDATE_BETA = 0; % and set beta using
+% options.beta = 1e3;
+% Manually tuning this parameter can give significantly better results.
+
+% options.initial_rank = 'auto'; % This sets to the maximum possible rank
+options.initial_rank = 50; % or we can set a value.
+
+% Provide original matrix for simulations
+options.X_true = X_true;
+
+
+
+fprintf('Running VB matrix completion...\n');
+tic
+[X_hat, A_hat, B_hat] = VBMC(P, Y, options);
+t_total = toc;
+
+% Show results
+err = norm( X_hat - X_true, 'fro' ) / norm( X_true, 'fro');
+r_hat = rank(X_hat);
+fprintf('\n\n-------------------------------------------------------------------\n')
+fprintf('Dimensions (%d, %d), true rank = %d, measurements = %g %%, SNR = %g\n',m, n, r, p*100, SNR);
+fprintf('Rel. error = %g, estimated rank = %d, time = %g\n', err, r_hat, t_total);
+fprintf('-------------------------------------------------------------------\n')
+
+
diff --git a/others/VBLRMat/RunVBRPCA.m b/others/VBLRMat/RunVBRPCA.m
new file mode 100644
index 0000000..358db0f
--- /dev/null
+++ b/others/VBLRMat/RunVBRPCA.m
@@ -0,0 +1,81 @@
+%% An example run of VBRPCA
+clear
+clc
+
+%% Create a low-rank + sparse matrix
+% Dimensions
+m = 500;
+n = 500;
+r = 25; % rank of the low-rank component
+
+A = randn(m,r);
+B = randn(r,n);
+X_true = A*B; % low rank component
+
+SparseRatio = 0.05;
+E_true = zeros(size(X_true));
+Eomega = randsample(m*n, round(m*n*SparseRatio));
+E_true(Eomega) = -10+20*rand(length(Eomega),1); % sparse component
+
+% Observation
+Y = X_true + E_true;
+
+% Add noise?
+sigma = 0; % no noise
+% sigma = 1e-3; % some noise
+% sigma = sqrt(1e-3); % more noise
+
+Y = Y + sigma*randn(size(Y));
+
+SNR = 10*log10( sum(sum((X_true+E_true).^2)) / sum(sum((Y-X_true-E_true).^2)));
+
+
+%% Run VRBPCA
+% all options are *optional*, everything will be set automatically
+% you can modify these options to get better performance
+options.verbose = 1;
+options.initial_rank = 'auto'; % This sets to the maximum possible rank
+% options.initial_rank = 300; % or we can use a value.
+options.X_true = X_true;
+options.E_true = E_true;
+options.inf_flag = 2; % inference flag for the sparse component
+% 1 for standard VB, 2 for MacKay. MacKay generally converges faster.
+
+options.MAXITER = 200;
+%Estimate noise variance? (beta is inverse noise variance)
+options.UPDATE_BETA = 1;
+% If the noise inv. variance is not to be estimated, set
+% options.UPDATE_BETA = 0; % and set beta using
+% options.beta = 1e3;
+
+% Select the optimization mode:
+% 'VB': fully Bayesian inference (default)
+% 'VB_app': fully Bayesian with covariance approximation
+% 'MAP': maximum a posteriori (covariance is set to 0)
+options.mode = 'VB';
+% For large scale problems, set to 'VB_app'.
+% options.mode = 'VB_app';
+
+
+fprintf('Running VBRPCA...\n')
+tic
+[X_hat, A_hat, B_hat, E_hat] = VBRPCA(Y,options);
+t_total = toc;
+
+% Results
+Xerr = norm( X_hat - X_true, 'fro' ) / norm( X_true, 'fro');
+Eerr = norm( E_hat - E_true, 'fro' ) / norm( E_true, 'fro');
+
+r_hat = rank(X_hat);
+% the model does not allow for "exact" sparsity, but only "soft" sparsity.
+% That is, most coefficients are almost zero (very small values), but not
+% exactly zero. If the sparse coefficients are needed, one can use
+% thresholding such as
+nonzero_E = length(find(abs(E_hat)>1e-7));
+
+% Show results
+fprintf('\n\n-------------------------------------------------------------------\n')
+fprintf('Dimensions (%d, %d), true rank = %d, nonzero E elements = %d, SNR = %g\n',m, n, r, length(Eomega), SNR);
+fprintf('low rank comp. error = %g, sparse comp. error %g, est. rank = %d,\nnonzero E elements = %d, running time = %g\n', Xerr, Eerr, r_hat, nonzero_E, t_total);
+fprintf('-------------------------------------------------------------------\n')
+
diff --git a/others/VBLRMat/VBMC.m b/others/VBLRMat/VBMC.m
new file mode 100644
index 0000000..b415bcd
--- /dev/null
+++ b/others/VBLRMat/VBMC.m
@@ -0,0 +1,395 @@
+function [X, A, B, t, varargout] = VBMC(P, Y, options)
+% Variational Bayesian Low Rank Matrix Completion
+% Author (a.k.a. person to blame): S. Derin Babacan
+%
+% Last updated: January 31, 2012
+% Usage:
+% [X, A, B, t, varargout] = VSBL_MAT(P, Y, options)
+% Solves for X = AB' in
+% Y = P(X) + N = P(AB') + N where X, A, B are low rank,
+% N is dense Gaussian noise, and P is a binary sampling matrix.
+%
+% If you're using this code, please acknowledge
+% Reference:
+% S. D. Babacan, M. Luessi, R. Molina, and A. K. Katsaggelos,
+% "Sparse Bayesian Methods for Low-Rank Matrix Estimation,"
+% IEEE Transactions on Signal Processing, 2012.
+%
+% This code is not optimized for speed. It implements the full variational
+% Bayesian inference described in the paper, without any manipulation of
+% the matrices. These are described in the paper above but not implemented
+% in this version of the code. These manipulations might lead to significant
+% decrease in running times depending on the size of matrix Y.
+%
+% -----------------------------------------------------------------------
+% INPUTS
+% P : binary (0-1) sampling matrix same size as Y
+% Y : input matrix
+% options : All options are *optional*. The default values are set
+% automatically.
+%
+% verbose : output the progress? (0/1) default: 0.
+% init : Initialization method.
+% 'rand': initialize A and B with random matrices
+% 'ml' : Apply SVD to Y and initialize A and B using its factors.
+% (default)
+% a_gamma0
+% b_gamma0 : hyperparameter values for gamma. default: 1e-6
+%
+%
+% a_beta0
+% b_beta0 : hyperparameter values for beta. default: 0
+%
+% inf_flag : Flag for inference of components of E. default: 1
+% 1: Standard
+% 2: fixed point MacKay
+% MAXITER : max number of iterations. default: 100
+% UPDATE_BETA : update noise variance? (0/1). default: 1
+% UPDATE_BETA_START : iteration number to start updating the noise variance
+% DIM_RED : prune irrelevant dimensions during iterations? (0/1). default: 1
+% DIMRED_THR : threshold to remove columns from A and B. Used only if
+% DIM_RED is 1. default: 1e4
+% initial_rank : The initial rank of A and B to start. Default is full
+% rank, calculated from the SVD of Y. A smaller value is
+% helpful if the data dimensions is large.
+% X_true,E_true : Original X and E matrices for simulations.
+% If supplied and verbose=1,
+% the errors are reported at each iteration
+% -----------------------------------------------------------------------
+% OUTPUTS
+% X : low rank component
+% A, B : low rank factors of X = AB'
+% optional outputs (see the paper for the full definitions):
+% varargout{1} = gamma --- hyperparameters for A and B
+% varargout{2} = beta --- noise inverse variance
+% varargout{3} = it --- number of iterations
+% varargout{4} = t --- running time
+% varargout{5} = Sigma_A --- covariance matrix of A
+% varargout{6} = Sigma_B --- covariance matrix of B
+% -----------------------------------------------------------------------
+%
+% Copyright (c): S. Derin Babacan, 2011
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+
+%% Check input variables,
+if ~exist('options','var')
+ options = populate_vars([]);
+else
+ options = populate_vars(options);
+end
+
+verbose = options.verbose;
+a_gamma0 = options.a_gamma0;
+b_gamma0 = options.b_gamma0;
+a_beta0 = options.a_beta0;
+b_beta0 = options.b_beta0;
+MAXITER = options.MAXITER;
+inf_flag = options.inf_flag;
+UPDATE_BETA = options.UPDATE_BETA;
+DIMRED = options.DIMRED;
+
+if DIMRED,
+ if isfield(options, 'DIMRED_THR')
+ DIMRED_THR = options.DIMRED_THR;
+ else
+ DIMRED_THR = 1e4;
+ end
+end
+
+UPDATE_BETA_START = options.UPDATE_BETA_START;
+
+
+% For synthetic simulations
+if isfield(options, 'X_true')
+ X_true = options.X_true;
+end
+
+
+
+%% Initialization
+[m n] = size(Y);
+mn = m*n;
+pmn = length(find(P == 1));
+p = pmn/mn;
+
+Y2sum = sum(Y(:).^2);
+scale2 = Y2sum / (mn);
+scale = sqrt(scale2);
+
+% Initialize A and B
+switch options.init,
+
+ % Maximum likelihood
+ case 'ml'
+ [U, S, V] = svd(Y, 'econ');
+
+ if strcmp(options.initial_rank, 'auto')
+ r = min([m,n]);
+ else
+ r = options.initial_rank;
+ end
+
+ A = U(:,1:r)*(S(1:r,1:r)).^(0.5);
+ B = (S(1:r,1:r)).^(0.5)*V(:,1:r)'; B = B';
+
+ Sigma_A = repmat( scale*eye(r,r), [1 1 m] );
+ Sigma_B = repmat( scale*eye(r,r), [1 1 n] );
+
+ gammas = 1./diag(S.^2);
+ gammas = gammas(1:r);
+ gammas = ones(r,1)*scale;
+
+ ssscale2 = 1;
+ beta = ssscale2./scale2;
+
+ if UPDATE_BETA == 0 & isfield(options, 'beta')
+ beta = options.beta;
+ end
+
+ % Random initialization
+ case 'rand'
+
+ if strcmp(options.initial_rank, 'auto')
+ r = min([m,n]);
+ else
+ r = options.initial_rank;
+ end
+
+ A = randn(m,r) * sqrt(scale);
+ B = randn(n,r) * sqrt(scale);
+ gammas = scale*ones(r,1);
+
+ Sigma_A = repmat( scale*eye(r,r), [1 1 m] );
+ Sigma_B = repmat( scale*eye(r,r), [1 1 n] );
+
+ if UPDATE_BETA == 0 & isfield(options, 'beta')
+ beta = options.beta;
+ else
+ beta = 1./scale2;
+ end
+
+
+end
+
+X = A*B';
+
+
+%% Iterations
+tic
+for it=1:MAXITER,
+
+ old_X = X;
+
+ Aw = diag(gammas);
+
+ %% A step
+ diagsA = 0;
+ for i=1:m, %iterate over rows
+
+ observed = find(P(i,:));
+ Bi = B(observed,:);
+ Sigma_A(:,:,i) = (beta*Bi'*Bi + beta*sum(Sigma_B(:,:,observed),3) + Aw)^(-1);
+ A(i,:) = beta*Y(i,observed)*Bi*Sigma_A(:,:,i);
+ diagsA = diagsA + diag( Sigma_A(:,:,i) );
+
+ end
+
+
+ %% B step
+ diagsB = 0;
+ for j=1:n, %Iterate over cols
+
+ observed = find(P(:,j));
+ Aj = A(observed,:);
+ Sigma_B(:,:,j) = (beta*Aj'*Aj + beta*sum(Sigma_A(:,:,observed),3) + Aw)^(-1);
+ B(j,:) = beta*Y(observed,j)'*Aj*Sigma_B(:,:,j);
+ diagsB = diagsB + diag( Sigma_B(:,:,j) );
+
+ end
+
+ %% update X
+ X = A*B';
+
+ %% estimate gammas
+ % Choose inference method (fixed point Mackay or standard)
+ if inf_flag == 1, %Standard
+ gammas = (m + n + a_gamma0)./( diag(B'*B) + diag(sum(Sigma_B,3)) + diag(A'*A)+ diag(sum(Sigma_A,3))+ b_gamma0);
+ elseif inf_flag == 2, %Mackay
+ gammas = (m + n + a_gamma0 - gammas.*( diag(sum(Sigma_A,3)) + diag(sum(Sigma_B,3)) ))./( diag(A'*A) + diag(B'*B)+b_gamma0);
+ end
+
+
+ %% estimate beta
+ if UPDATE_BETA & it>=UPDATE_BETA_START,
+
+
+ err = sum(sum( abs(Y - P.*X).^2 ) );
+
+ for l = 1: m
+ observed = find(P(l, :));
+ err = err + trace(B(observed,:)'*B(observed,:)*Sigma_A(:, :, l)) ...
+ + trace(A(l, :)'*A(l, :) *sum(Sigma_B(:, :, observed), 3)) ...
+ + trace( Sigma_A(:, :, l)*sum(Sigma_B(:, :, observed), 3));
+
+ end
+
+ beta = (pmn + a_beta0)/(err+b_beta0);
+ end
+
+ %% Prune irrelevant dimensions?
+ if DIMRED,
+ MAX_GAMMA = min(gammas) * DIMRED_THR;
+
+ if sum(find(gammas > MAX_GAMMA)),
+
+ indices = find(gammas <= MAX_GAMMA);
+
+ A = A(:,indices);
+ B = B(:,indices);
+ gammas = gammas(indices);
+
+ Sigma_A = Sigma_A(indices,indices,:);
+ Sigma_B = Sigma_B(indices,indices,:);
+
+ [m r] = size(A);
+ [n r] = size(B);
+ end
+
+ end
+
+ %% Check convergence and display progress
+ conv = sum(sum( ( old_X - X ).^2 ) ) / mn;
+
+ if verbose,
+
+ if exist('X_true','var')
+ % For synthetic simulations
+ err = sum(sum( ( X - X_true ).^2 ) ) / sum(sum( ( X_true ).^2 ) );
+ rms = sqrt( sum(sum( ( X - X_true ).^2 ) )/mn );
+ relrec = norm(X_true-X,'fro')/norm(X_true,'fro');
+
+ fprintf('it %d: Error = %g, beta = %g, conv = %g, r = %d\n', it, relrec, beta, conv, r);
+
+ else
+
+ fprintf('it %d: beta = %g, conv = %g, r = %d\n', it, beta, conv, r);
+
+ end
+ end
+
+ % Check for convergence (do at least 5 iterations)
+ if it > 5 & conv < options.thr
+ break;
+ end
+
+
+end
+t = toc;
+
+varargout{1} = gammas;
+varargout{2} = beta;
+varargout{3} = it;
+varargout{4} = t;
+varargout{5} = Sigma_A;
+varargout{6} = Sigma_B;
+
+
+%% Function to populate options.
+function [options] = populate_vars(options)
+
+
+if isempty(options)
+ options.init = 'ml';
+ options.verbose = 1;
+ options.a_gamma0 = 1e-6;
+ options.b_gamma0 = 1e-6;
+ options.a_beta0 = 0;
+ options.b_beta0 = 0;
+ options.MAXITER = 100;
+ options.DIMRED = 1;
+ options.UPDATE_BETA = 1;
+ options.mode = 'VB';
+ options.thr = 1e-7;
+ options.initial_rank = 'auto';
+ options.inf_flag = 1;
+ options.UPDATE_BETA_START = 1;
+else
+
+ if ~isfield(options, 'init')
+ options.init = 'ml';
+ end
+
+ if ~strcmp(options.init,'ml') | ~strcmp(options.init,'rand')
+ options.init = 'ml';
+ end
+
+ if ~isfield(options, 'verbose')
+ options.verbose = 1;
+ end
+
+ if ~isfield(options, 'a_gamma0')
+ options.a_gamma0 = 1e-6;
+ end
+
+ if ~isfield(options, 'b_gamma0')
+ options.b_gamma0 = 1e-6;
+ end
+
+ if ~isfield(options, 'a_beta0')
+ options.a_beta0 = 0;
+ end
+
+ if ~isfield(options, 'b_beta0')
+ options.b_beta0 = 0;
+ end
+
+ if ~isfield(options, 'MAXITER')
+ options.MAXITER = 100;
+ end
+
+ if ~isfield(options, 'DIMRED')
+ options.DIMRED = 1;
+ end
+
+ if ~isfield(options, 'inf_flag')
+ options.inf_flag = 1;
+ end
+
+ if options.inf_flag ~= 1 & options.inf_flag ~= 2,
+ options.inf_flag = 1; % Standard inference
+ end
+
+ if ~isfield(options, 'UPDATE_BETA')
+ options.UPDATE_BETA = 1;
+ end
+
+ if ~isfield(options, 'UPDATE_BETA_START')
+ options.UPDATE_BETA_START = 1;
+ end
+
+
+ if ~isfield(options, 'mode')
+ options.mode = 'VB';
+ end
+
+ if ~isfield(options, 'thr')
+ options.thr = 1e-7;
+ end
+
+ if ~isfield(options, 'initial_rank')
+ options.initial_rank = 'auto';
+ end
+
+end
+
+
+
diff --git a/others/VBLRMat/VBRPCA.m b/others/VBLRMat/VBRPCA.m
new file mode 100644
index 0000000..e7755af
--- /dev/null
+++ b/others/VBLRMat/VBRPCA.m
@@ -0,0 +1,392 @@
+function [X, A, B, E, varargout] = VBRPCA(Y, options)
+% Variational Bayesian Robust Principal Component Analysis
+% Author (a.k.a. person to blame): S. Derin Babacan
+% Last updated: January 31, 2012
+% Usage:
+% [X, A, B, E, varargout] = VBRPCA(Y, options)
+% Solves for A,B,E in
+% Y = X + E + N
+% = AB' + E + N
+% where X, A, B are low rank and E is sparse, N is
+% dense Gaussian noise
+%
+% If you're using this code, please acknowledge
+% S. D. Babacan, M. Luessi, R. Molina, and A. K. Katsaggelos,
+% "Sparse Bayesian Methods for Low-Rank Matrix Estimation,"
+% IEEE Transactions on Signal Processing, 2012.
+%
+% This code is not optimized for speed.
+% -----------------------------------------------------------------------
+% INPUTS
+% Y : input matrix
+% options : All options are *optional*. The default values are set
+% automatically.
+%
+% verbose : output the progress? (0/1) default: 0.
+% init : Initialization method.
+% 'rand': initialize A and B with random matrices
+% 'ml' : Apply SVD to Y and initialize A and B using its factors.
+% (default)
+% a_gamma0
+% b_gamma0 : hyperparameter values for gamma. default: 1e-6
+%
+% a_alpha0
+% b_alpha0 : hyperparameter values for alpha. default: 0
+%
+% a_beta0
+% b_beta0 : hyperparameter values for beta. default: 0
+%
+% inf_flag : Flag for inference of components of E. default: 1
+% 1: Standard
+% 2: fixed point MacKay
+% MAXITER : max number of iterations. default: 100
+% UPDATE_BETA : update noise variance? (0/1). default: 1
+% DIM_RED : prune irrelevant dimensions during iterations? (0/1). default: 1
+% DIMRED_THR : threshold to remove columns from A and B. Used only if
+% DIM_RED is 1. default: 1e4
+% initial_rank : The initial rank of A and B to start. Default is full
+% rank, calculated from the SVD of Y. A smaller value is
+% helpful if the data dimensions is large.
+% X_true,E_true : Original X and E matrices for simulations.
+% If supplied and verbose=1,
+% the errors are reported at each iteration
+% -----------------------------------------------------------------------
+% OUTPUTS
+% X : low rank component of Y
+% A, B : low rank factors in X = AB'
+% E : sparse component of Y
+% optional outputs (see the paper for the full definitions):
+% varargout{1} = gamma --- hyperparameters for A and B
+% varargout{2} = alpha --- hyperparameters for E
+% varargout{3} = beta --- noise inverse variance
+% varargout{4} = it --- number of iterations
+% varargout{5} = t --- running time
+% varargout{6} = Sigma_A --- covariance matrix of A
+% varargout{7} = Sigma_B --- covariance matrix of B
+% -----------------------------------------------------------------------
+%
+%
+% Copyright (c): S. Derin Babacan, 2011
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+
+
+
+%% Check input variables
+if ~exist('options','var')
+ options = populate_vars([]);
+else
+ options = populate_vars(options);
+end
+
+verbose = options.verbose;
+a_gamma0 = options.a_gamma0;
+b_gamma0 = options.b_gamma0;
+a_beta0 = options.a_beta0;
+b_beta0 = options.b_beta0;
+a_alpha0 = options.a_alpha0;
+b_alpha0 = options.b_alpha0;
+MAXITER = options.MAXITER;
+inf_flag = options.inf_flag;
+UPDATE_BETA = options.UPDATE_BETA;
+DIMRED = options.DIMRED;
+
+if DIMRED,
+ if isfield(options, 'DIMRED_THR')
+ DIMRED_THR = options.DIMRED_THR;
+ else
+ DIMRED_THR = 1e4;
+ end
+end
+
+% For synthetic simulations
+if isfield(options, 'X_true') & isfield(options, 'E_true')
+ X_true = options.X_true;
+ E_true = options.E_true;
+end
+
+
+%% Initialization
+[m n] = size(Y);
+mn = m*n;
+
+Y2sum = sum(Y(:).^2);
+scale2 = Y2sum / (mn);
+scale = sqrt(scale2);
+
+% Initialize A, B and E
+switch options.init,
+ % Maximum likelihood
+ case 'ml'
+
+ [U, S, V] = svd(Y, 'econ');
+
+ if strcmp(options.initial_rank, 'auto')
+ r = min([m,n]);
+ else
+ r = options.initial_rank;
+ end
+
+ A = U(:,1:r)*(S(1:r,1:r)).^(0.5);
+ B = (S(1:r,1:r)).^(0.5)*V(:,1:r)'; B = B';
+
+ Sigma_A = scale*eye(r,r);
+ Sigma_B = scale*eye(r,r);
+ gammas = scale*ones(r,1);
+
+ beta = 1./scale2;
+
+ if UPDATE_BETA == 0 & isfield(options, 'beta')
+ beta = options.beta;
+ end
+
+ Sigma_E = scale*ones(r,r);
+ alphas = ones(m,n)*scale;
+
+ % Random initialization
+ case 'rand'
+ r = min([m,n]);
+ A = randn(m,r) * sqrt(scale);
+ B = randn(n,r) * sqrt(scale);
+ gammas = scale*ones(r,1);
+ Sigma_A = eye(r) * scale;
+ Sigma_B = eye(r) * scale;
+ E = randn(m, n) * sqrt(scale);
+ beta = 1./scale2;
+ Sigma_E = scale*ones(r,r);
+ alphas = ones(m,n)*scale;
+
+end
+
+X = A*B';
+E = Y-X;
+
+%% Iterations
+tic
+for it = 1:MAXITER,
+ old_X = X;
+ old_E = E;
+
+ %% update X
+ W = diag(gammas);
+ %% A step
+
+ if strcmp(options.mode, 'VB')
+ Sigma_A = (beta*B'*B + W + beta*n*Sigma_B)^(-1);
+ A = beta*(Y-E)*B*Sigma_A;
+ elseif strcmp(options.mode, 'VB_app')
+ Sigma_A = (beta*B'*B + W + beta*n*Sigma_B);
+ A = (beta*(Y-E)*B ) / Sigma_A;
+ Sigma_A = diag(1./diag(Sigma_A));
+ elseif strcmp(options.mode, 'MAP')
+ Sigma_A = (beta*B'*B + W + beta*n*Sigma_B);
+ A = (beta*(Y-E)*B ) / Sigma_A;
+ Sigma_A = zeros(size( Sigma_A ) );
+ end
+
+ %% B step
+ if strcmp(options.mode, 'VB')
+ Sigma_B = (beta*A'*A + W + beta*m*Sigma_A)^(-1);
+ B = beta*(Y-E)'*A*Sigma_B;
+ elseif strcmp(options.mode, 'VB_app')
+ Sigma_B = (beta*A'*A + W + beta*m*Sigma_A);
+ B = (beta*(Y-E)'*A)/Sigma_B;
+ Sigma_B = diag(1./diag(Sigma_B));
+ elseif strcmp(options.mode, 'MAP')
+ Sigma_B = (beta*A'*A + W + beta*m*Sigma_A);
+ B = (beta*(Y-E)'*A)/Sigma_B;
+ Sigma_B = zeros(size( Sigma_B ) );
+ end
+
+ X = A*B';
+
+ %% E-step
+ Sigma_E = 1./(alphas+beta);
+ E = beta*(Y-X).*Sigma_E;
+
+
+ %% Estimate alphas
+ %----------------------------------------------------------------------------------------
+ % Choose inference method (fixed point Mackay or standard)
+ if inf_flag == 1,
+ % Standard
+ alphas = 1./(E.^2 + Sigma_E);
+ elseif inf_flag == 2,
+ % Mackay
+ alphas = (1-alphas.*Sigma_E + a_alpha0)./ (E.^2 + eps + b_alpha0);
+ end
+
+
+ %----------------------------------------------------------------------------------------
+ %% estimate gammas
+ gammas = (m + n + a_gamma0)./( diag(B'*B) + diag(A'*A)+ m*diag(Sigma_A) + n*diag(Sigma_B)+ b_gamma0);
+
+ % MacKay update for gamma? (not tested extensively)
+% gammas = ( m+n - m*diag(Sigma_A) - n*diag(Sigma_B) + a_gamma0 )./ (diag(B'*B) + diag(A'*A)+ b_gamma0);
+
+ %% estimate beta?
+ if UPDATE_BETA,
+ err = sum(sum( abs(Y - X - E).^2 ) ) + n*trace(A'*A*Sigma_B) + m*trace(B'*B*Sigma_A) + m*n*trace(Sigma_A*Sigma_B) + sum(sum((Sigma_E)));
+ beta = (m*n + a_beta0)/(err+b_beta0);
+ end
+
+ %% Prune irrelevant dimensions?
+ if DIMRED,
+ MAX_GAMMA = min(gammas) * DIMRED_THR;
+
+ if sum(find(gammas > MAX_GAMMA)),
+ indices = find(gammas <= MAX_GAMMA);
+
+ A = A(:,indices);
+ B = B(:,indices);
+ gammas = gammas(indices);
+
+ Sigma_A = Sigma_A(indices,indices);
+ Sigma_B = Sigma_B(indices,indices);
+
+
+ [m r] = size(A);
+ [r n] = size(B'); % r is the current rank estimate
+ end
+
+ end
+
+
+ %% Check convergence and display progress
+ Xconv = norm(old_X-X,'fro')/norm(old_X,'fro');
+
+ if verbose,
+
+ Econv = norm(old_E - E,'fro')/ norm(old_E,'fro');
+
+ if exist('X_true','var') % For synthetic simulations
+ Xerr = norm( X - X_true , 'fro') / norm( X_true , 'fro');
+ Eerr = norm( E - E_true , 'fro') / norm( E_true , 'fro');
+ fprintf('it %d: X Error = %g, Xconv = %g, E Error = %g, Econv = %g, beta = %g, r = %d\n', it, Xerr, Xconv, Eerr, Econv, beta, r);
+
+ else
+ fprintf('it %d: Xconv = %g, Econv = %g, beta = %g, r = %d\n', it, Xconv, Econv, beta, r);
+
+ end
+ end
+
+ % Check for convergence (do at least 5 iterations)
+ if it > 5 & Xconv < options.thr
+ break;
+ end
+
+end
+t = toc;
+
+varargout{1} = gammas;
+varargout{2} = alphas;
+varargout{3} = beta;
+varargout{4} = it;
+varargout{5} = t;
+varargout{6} = Sigma_A;
+varargout{7} = Sigma_B;
+
+
+
+%% Function to populate options.
+function [options] = populate_vars(options)
+
+if isempty(options),
+
+ options.init = 'ml';
+ options.verbose = 1;
+ options.a_gamma0 = 1e-6;
+ options.b_gamma0 = 1e-6;
+ options.a_alpha0 = 0;
+ options.b_alpha0 = 0;
+ options.a_beta0 = 0;
+ options.b_beta0 = 0;
+ options.MAXITER = 100;
+ options.DIMRED = 1;
+ options.inf_flag = 2;
+ options.UPDATE_BETA = 1;
+ options.mode = 'VB';
+ options.thr = 1e-5;
+ options.initial_rank = 'auto';
+
+else
+
+ if ~isfield(options, 'init')
+ options.init = 'ml';
+ end
+
+ if ~strcmp(options.init,'ml') | ~strcmp(options.init,'rand')
+ options.init = 'ml';
+ end
+
+ if ~isfield(options, 'verbose')
+ options.verbose = 0;
+ end
+
+ if ~isfield(options, 'a_gamma0')
+ options.a_gamma0 = 1e-6;
+ end
+
+ if ~isfield(options, 'b_gamma0')
+ options.b_gamma0 = 1e-6;
+ end
+
+ if ~isfield(options, 'a_alpha0')
+ options.a_alpha0 = 0;
+ end
+
+ if ~isfield(options, 'b_alpha0')
+ options.b_alpha0 = 0;
+ end
+
+ if ~isfield(options, 'a_beta0')
+ options.a_beta0 = 0;
+ end
+
+ if ~isfield(options, 'b_beta0')
+ options.b_beta0 = 0;
+ end
+
+ if ~isfield(options, 'MAXITER')
+ options.MAXITER = 100;
+ end
+
+ if ~isfield(options, 'DIMRED')
+ options.DIMRED = 1;
+ end
+
+ if ~isfield(options, 'inf_flag')
+ options.inf_flag = 2;
+ end
+
+ if options.inf_flag ~= 1 | options.inf_flag ~= 2,
+ options.inf_flag = 2;
+ end
+
+ if ~isfield(options, 'UPDATE_BETA')
+ options.UPDATE_BETA = 1;
+ end
+
+ if ~isfield(options, 'mode')
+ options.mode = 'VB';
+ end
+
+ if ~isfield(options, 'thr')
+ options.thr = 1e-5;
+ end
+
+ if ~isfield(options, 'initial_rank')
+ options.initial_rank = 'auto';
+ end
+
+end
+
+