function [lambda_new, psi_new] = mstep(beta, s, lambda_old, psi_old, eye_k)

bs = beta*s;
lambda = bs'*inv(eye_k - beta*lambda_old + bs * beta');
psi = diag(diag(s - lambda*bs));

lambda_new = lambda;
psi_new = psi;