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

%% 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

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

%% 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 rand*0 < (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

%% 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
        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

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

%% 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