In the plots below the result from the code is plotted. The first plot show loglikelihood, the second shows the x-estimate and the third shows the error compared to x. I used logmvnpdf found here: http://www.mathworks.com/matlabcentral/fileexchange/34064-log-multivariate-normal-distribution-function/content/logmvnpdf.m [![Log-likelihood, x-estimate, x-error][1]][1] %% Get initial guess clc mm = mean(y); x = (A\mm'); mu = A*x; x0 = x; RR = x'*V*x + lambda; sigma = kron(RR, eye(m)); p_old = sum(logmvnpdf(y, mu', sigma)); %% Estimate s = 1e-1; p_old = -1e190; counter = 0; p_save = 0; x_save = x0; while counter < 4*1e3 if mod(counter, 1e3) == 0 s = s/10 end counter = counter + 1; x = x0 + randn(k, 1)*s; mu = A*x; RR = x'*V*x + lambda; sigma = kron(RR, eye(m)); try % fails if sigma is not posdef p_ = sum(logmvnpdf(y, mu', sigma)); catch p_ = -1; counter = max(1, counter - 1); end if p_ > p_old p_old = p_; x0 = x; end p_save(counter) = p_old; x_save(:, counter) = x0; end rr = zeros(counter, 1); for c = 1:counter rr(c, :) = norm(x_save(:, c) - x_true); end [1]: https://i.sstatic.net/Q9tL6.png