Maximization of a nasty Gaussian likelihood

Question

I asked this question in math.SE before. One answer so far, and we were unable to reach a conclusion. It is more related to statistics, so I wanted to post this here.

I have a Gaussian likelihood function, $$p(y|x) = \mathcal{N}(y; Ax, (x^\top V x + \lambda) \otimes I)$$ where $A,V,\lambda$ is known, and $\otimes$ is the Kronecker product. (the notation indicates that covariance is a scalar times identity matrix -- scalar is: $x^\top V x + \lambda$). Note that $A$ is a rectangular matrix, say $m\times k$ with $m>k$, hence $V$ is a $k\times k$ matrix. $\lambda > 0$. $x$ is $k\times 1$ and y is $m\times 1$.

I would like to maximise this with respect to $x$, in other words, solve the following problem, $$x^* = \arg \max_x p(y|x) = \arg \max_x \mathcal{N}(y;Ax,(x^\top V x + \lambda) \otimes I)$$ I tried to take derivative of the log-likelihood and set it to zero, however I was unable to leave out $x$ and obtain an exact solution.

I wonder if there is an exact solution, and if not: what the best numerical scheme (one suggests) is to overcome this problem.

Any help is greatly appreciated. Thanks!

PS: Pseudoinverse is not the solution, according to 2D numerical simulations! And another empirical observation from 2D simulations: As $\lambda \to \infty$ (for very large values), pseudoinverse solution becomes more and more accurate, so this hints about structure of the solution a bit.

Answer 1

This is basically a regression problem, only with a constraint on the variance. If we didn't have $\sigma^2 = x^TVx + \lambda$, the solution would be $(A^TA)^{-1}A^TY$, but that solution almost undoubtedly won't satisfy $s^2=\hat{x}^TV\hat{x} + \lambda$, where $s^2$ is the mean square error from the regression.

There must be some way to iterate towards a solution that starts with OLS and moves towards something that satisfies, or is close to, the variance condition. Would a Lagrange multiplier work?

Minimize $$ \sum_{i=1}^N (y_i - A_ix)^2 + \alpha(s^2 - x^TVx + \lambda) $$

where $\alpha$ is the Lagrange multiplier. Taking the derivative of $s^2$ would be painful, but you might hold that constant at each iteration, updating it as you update $x$.

Answer 2

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