# Global optimizers handling minimization of an expression arising from the likelihood of a multivariate normal

I am interested in converting the following optimisation problem into a form that an exponential cone and/or SDP solver such as MOSEK can handle. This is a multivariate version of the question I previously asked here.

We are given the following:

• $$A\in\mathbb{R}^{M\times N}$$ (with $$N\ge M$$), a full rank matrix.
• $$b\in\mathbb{R}^M$$.
• $$C\in\mathbb{R}^{N\times K}$$ (with $$N\ge K$$), a full rank matrix.
• $$d\in\mathbb{R}^M$$.

(If it makes the problem any easier, you may assume that all elements of $$C$$ are either $$0$$ or $$1$$, with precisely one $$1$$ per row.)

$$z\in\mathbb{R}^K$$ is our optimisation variable.

Let $$u:=Cz+d$$ and let $$D\in\mathbb{R}^{N\times N}$$ be the diagonal matrix with diagonal $$\exp(u_1),\dots,\exp(u_N)$$.

We want to choose $$z$$ to minimize: $$\log{\det{\left(ADA^\top\right)}}+b^\top{\left(ADA^\top\right)}^{-1}b.$$

If $$N=M$$, then this problem is easy, as the objective simplifies to: $$\log{{\left(\det{A}\right)}^2}+\sum_{n=1}^N{u_n}+b^\top A^{-\top} D^{-1} A^{-1} b,$$ which is a standard exponential cone problem. YALMIP + MOSEK seems to handle it fine:

z=sdpvar(3,1); C=randn(10,3); d=randn(10,1); u=C*z+d; A=randn(4,10); b=randn(4,1); InvAb2=(A\b).^2;
optimize([],sum(u)+sum(InvAb2.*exp(-u)),sdpsettings('solver','mosek','verbose',1))
...
info: 'Successfully solved (MOSEK-CONE)'


My question is whether we can convert the objective into a form MOSEK (or similar) can handle when $$N>M$$.

The first term can be handled using the formula for the determinant here. This implies (using notation from here, but with sets in bold): $$\det{\left(ADA^\top\right)}=\sum_{\textbf{S} \in \binom{[N]}{M}} {\left[\det(A_{[M],\textbf{S}})^2\exp{\sum_{s \in \textbf{S}}{u_s}}\right]},$$ so $$\log{\det{\left(ADA^\top\right)}}$$ can be written as a $$\operatorname{logsumexp}$$ expression.

Replacing the second term in the original objective with a dummy exponential, this seems to work fine:

z=sdpvar(3,1); C=randn(10,3); d=randn(10,1); u=C*z+d; A=randn(4,10); b=randn(4,1); S=nchoosek(1:10,4); logdetADATerms=sdpvar(size(S,1),1);
for s = 1 : size(S,1); logdetADATerms(s)=log(det(A(:,S(s,:)))^2)+sum(u(S(s,:))); end;
...
info: 'Successfully solved (MOSEK-CONE)'


But of course this dummy exponential is not the original second term ($$b^\top{\left(ADA^\top\right)}^{-1}b$$).

Idea 1:

One idea to tackle the second term I had was to introduce a new auxiliary optimization variable $$f\in\mathbb{R}^N$$, with the constraint $$Af=b$$. Then the objective is equivalent to: $$\log{\det{\left(ADA^\top\right)}}+\sum_{n=1}^N{f_n^2 \exp{(-u_n)}}.$$ Perhaps unsurprisingly given the product term, this does not work though.

z=sdpvar(3,1); f=sdpvar(10,1); C=randn(10,3); d=randn(10,1); u=C*z+d; A=randn(4,10); b=randn(4,1); S=nchoosek(1:10,4); logdetADATerms=sdpvar(size(S,1),1);
for s = 1 : size(S,1); logdetADATerms(s)=log(det(A(:,S(s,:)))^2)+sum(u(S(s,:))); end;
Warning: Solver not applicable (mosek does not support multi-term monomial equalities)


Note though that with a (wrong!) modified objective without the exp(-u) it will solve:

z=sdpvar(3,1); f=sdpvar(10,1); C=randn(10,3); d=randn(10,1); u=C*z+d; A=randn(4,10); b=randn(4,1); S=nchoosek(1:10,4); logdetADATerms=sdpvar(size(S,1),1);
for s = 1 : size(S,1); logdetADATerms(s)=log(det(A(:,S(s,:)))^2)+sum(u(S(s,:))); end;
...
info: 'Successfully solved (MOSEK-CONE)'


Idea 2:

Another idea I had to handle the second term was to combine the expression for the inverse of a matrix in terms of its minors, with the Cauchy-Binet formula. This produce a linear combination of exponentials, divided by the determinant we already found. See the following code:

z=sdpvar(3,1); C=randn(10,3); d=randn(10,1); u=C*z+d; A=randn(4,10); b=randn(4,1); S=nchoosek(1:10,4); logdetADATerms=sdpvar(size(S,1),1);
for s = 1 : size(S,1); logdetADATerms(s)=log(det(A(:,S(s,:)))^2)+sum(u(S(s,:))); end;
T=nchoosek(1:10,3); logdetbbTerms=sdpvar(size(T,1),1); for t = 1 : size(T,1); bbTerm=0; for m1 = 1 : 4; L=setdiff(1:4,m1); for m2 = 1 : 4; R=setdiff(1:4,m2); bbTerm=bbTerm+b(m1)*b(m2)*(-1)^(m1+m2)*det(A(L,T(t,:)))*det(A(R,T(t,:))); end; end; assert( bbTerm > 0 ); logdetbbTerms(t)=log(bbTerm)+sum(u(T(t,:))); end;
...
info: 'Convexity requirements not met (MOSEK-CONE)'


It is easy to produce a convex relaxation of this problem though by placing the logsumexp into constraints. This gives:

z=sdpvar(3,1); C=randn(10,3); d=randn(10,1); u=C*z+d; A=randn(4,10); b=randn(4,1); S=nchoosek(1:10,4); logdetADATerms=sdpvar(size(S,1),1);
for s = 1 : size(S,1); logdetADATerms(s)=log(det(A(:,S(s,:)))^2)+sum(u(S(s,:))); end;
T=nchoosek(1:10,3); logdetbbTerms=sdpvar(size(T,1),1); for t = 1 : size(T,1); bbTerm=0; for m1 = 1 : 4; L=setdiff(1:4,m1); for m2 = 1 : 4; R=setdiff(1:4,m2); bbTerm=bbTerm+b(m1)*b(m2)*(-1)^(m1+m2)*det(A(L,T(t,:)))*det(A(R,T(t,:))); end; end; assert( bbTerm > 0 ); logdetbbTerms(t)=log(bbTerm)+sum(u(T(t,:))); end;
value( logbb-logsumexp(logdetbbTerms) )
...
info: 'Successfully solved (MOSEK-CONE)'
problem: 0
ans =
0.56829
ans =
7.9628e-09


So MOSEK can handle this convex relaxation, but it is not equivalent to the original problem in general, as the bound is not tight. (With these random coefficients, it appears to be tight over half the time at least.)

However, even if this was always tight, this would not really be a viable approach. This tiny problem ends up with 669 constraints, 331 cones and 1328 scalar variables. Plausible problems have $$M\approx 400$$ and $$N\approx 1600$$. This results in around $${10}^{389}$$ variables!

So perhaps my question ought to be:

Can we can convert the objective into a form MOSEK (or similar) can handle when $$N>M$$ with only polynomial in $$N,M,K$$ variables?

In the special case in which $$C=I$$ and $$d=0$$, the first order condition states that the diagonal of the matrix $$A^\top {\left(ADA^\top\right)}^{-1} {\left( ADA^\top - b b^\top \right)} {\left(ADA^\top\right)}^{-1} A$$ is zero.