Constrained Optimization for Multivariate Likelihoods

Why Constraints Matter in Likelihood Optimization

A numerical optimizer can produce a matrix that looks close to a covariance estimate but is not a valid covariance matrix. For a nonsingular multivariate Gaussian model, the covariance matrix must be symmetric positive-definite: every nonzero vector \(v\) must satisfy \(v^\mathsf{T}\Sigma v>0\). Positive-definiteness is a condition for the model to be valid, not a condition for every matrix operation to exist. A symmetric indefinite matrix can be invertible and can have positive determinant; a singular matrix has no inverse, and the real logarithm of its determinant is undefined when the determinant is nonpositive.

The practical task is therefore not just to optimize a likelihood, but to do so over a feasible set and in a way that remains numerically stable. This tutorial develops the main choices for covariance constraints, explains line-search and trust-region steps, and shows how to diagnose whether an algorithm has actually converged. The objective will usually be the negative log-likelihood, so optimization means minimizing it.

Cholesky Coordinates for Covariance Matrices

One direct way to enforce positive-definiteness is to optimize over a lower-triangular matrix \(L\) and set \(\Sigma=LL^\mathsf{T}\). If \(L\) is nonsingular, then \(v^\mathsf{T}LL^\mathsf{T}v=\|L^\mathsf{T}v\|^2>0\) for all nonzero \(v\), so the resulting covariance is positive-definite. Taking each diagonal entry of \(L\) to be an exponential of an unconstrained parameter ensures that it is positive.

The bijection follows from the Cholesky factorization. For completeness, its existence can be seen by induction on dimension. In dimension one, a positive scalar \(a\) factors as \((\sqrt a)^2\). For a larger positive-definite matrix, partition it as \(\begin{pmatrix}a&r^\mathsf{T}\\r&B\end{pmatrix}\), where \(a>0\). Its Schur complement \(B-rr^\mathsf{T}/a\) is positive-definite: for nonzero \(z\), apply positive-definiteness of the original matrix to \((-r^\mathsf{T}z/a,z^\mathsf{T})^\mathsf{T}\). By induction, that Schur complement has a Cholesky factor \(T\). Then \(L=\begin{pmatrix}\sqrt a&0\\r/\sqrt a&T\end{pmatrix}\) satisfies \(LL^\mathsf{T}=\Sigma\). This construction also determines each entry successively, giving uniqueness. Replacing each positive diagonal element by its logarithm gives unrestricted real coordinates. \(\square\)

Cholesky coordinates enforce positive-definiteness, but they do not automatically enforce every possible structural constraint. For example, a bound on a variance or correlation becomes a constraint on the entries of \(L\) and may need to be handled separately. Nor does setting selected off-diagonal entries of \(L\) to zero generally impose the corresponding covariance entries to be zero; the product \(LL^\mathsf{T}\) determines the covariance structure. Choose a parameterization that matches the actual model constraints.

A Gaussian Likelihood Example and Its Constrained Version

For independent Gaussian observations with known mean and positive-definite sample scatter \(S\), minimizing the negative log-likelihood over \(\Sigma\), apart from a positive constant factor, amounts to minimizing \(f(\Sigma)=\log\det\Sigma+\operatorname{tr}(\Sigma^{-1}S)\). Here \(S\) is the average outer product of the observations centered at the known mean. The next result identifies the unconstrained optimum and provides a useful benchmark for checking an optimizer.

Proof. Let \(A=S^{-1/2}\Sigma S^{-1/2}\), which is positive-definite, and let its eigenvalues be \(\lambda_1,\ldots,\lambda_p>0\). The determinant identity gives \(\log\det\Sigma=\log\det S+\log\det A\). The cyclic trace identity gives \(\operatorname{tr}(\Sigma^{-1}S)=\operatorname{tr}(A^{-1})\). Therefore \(f(\Sigma)-f(S)=\sum_{j=1}^{p}(\log\lambda_j+1/\lambda_j-1)\). For \(h(\lambda)=\log\lambda+1/\lambda-1\), its derivative is \(h'(\lambda)=(\lambda-1)/\lambda^2\). Thus \(h\) decreases up to \(1\), increases after \(1\), and has unique minimum \(h(1)=0\). The difference is nonnegative and is zero only when every eigenvalue of \(A\) is one, which for a symmetric positive-definite \(A\) means \(A=I\), or \(\Sigma=S\). \(\square\)

Constraints can change the optimum, even when the unconstrained estimate is known. A common example is a lower bound on eigenvalues to keep a covariance away from singularity. The constraint must be part of the statistical or numerical specification; it is not appropriate to add one silently just to make an optimizer behave.

Line Searches and Sufficient Decrease

Given a feasible parameter vector \(x\), an optimizer may choose a direction \(p\) and test points \(x+\alpha p\). For minimization, \(p\) is a descent direction when \(\nabla f(x)^\mathsf{T}p<0\). The Armijo condition accepts a trial step when \(f(x+\alpha p)\leq f(x)+c\alpha\nabla f(x)^\mathsf{T}p\), where \(0<c<1\). It requires a specified amount of decrease relative to the trial length and directional derivative.

A standard backtracking scheme starts from a chosen \(\alpha_0>0\), then repeatedly multiplies the trial length by \(\rho\in(0,1)\) until the Armijo condition holds. The cap on tested lengths in this procedure comes from the choice of \(\alpha_0\), not from Armijo itself. If constraints are handled in the original coordinates, every trial also needs to remain feasible; with a Cholesky parameterization, positive-definiteness is preserved for every finite parameter vector.

Proof. The Lipschitz-gradient condition implies the descent lemma \(f(x+\alpha p)\leq f(x)+\alpha\nabla f(x)^\mathsf{T}p+\frac{M}{2}\alpha^2\|p\|^2\). Write \(d=-\nabla f(x)^\mathsf{T}p>0\). The Armijo condition follows whenever \(\frac{M}{2}\alpha^2\|p\|^2\leq(1-c)\alpha d\). For \(M>0\), every positive \(\alpha\leq 2(1-c)d/(M\|p\|^2)\) satisfies this inequality. Backtracking produces trial lengths \(\alpha_0\rho^k\), which eventually fall below any fixed positive threshold. If \(M=0\), the descent lemma already gives the Armijo inequality for every positive trial length. Thus a trial is eventually accepted, provided its segment remains in the stated neighborhood. \(\square\)

Trust Regions and Convergence Diagnostics

A line search chooses a direction and then scales it. A trust-region method instead limits the length of a proposed step and checks whether the objective behaves as its local model predicted. For example, for a quadratic model \(m(p)=f(x)+g^\mathsf{T}p+\frac12p^\mathsf{T}Bp\), define the ratio \(\rho_{\mathrm{TR}}=[f(x)-f(x+p)]/[m(0)-m(p)]\), when the predicted decrease \(m(0)-m(p)\) is positive. A small or negative ratio means the model predicted improvement that did not occur, so the step is rejected or the trust region is reduced; a ratio near one indicates a good prediction and may justify enlarging the region. Steps must still satisfy any constraints not already guaranteed by the parameterization.

A small change in objective alone is weak evidence of convergence: an algorithm can stagnate with a nonzero gradient, or make tiny steps because of poor scaling. A useful diagnostic set for constrained likelihood optimization includes:

• The objective and its change across accepted iterations, with a check that minimization has not increased it beyond the method's stated tolerance.
• The gradient norm in the coordinates actually optimized, together with a constraint-aware stationarity check when constraints are active.
• Feasibility residuals, such as violations of bounds or equality constraints; for a covariance, also check symmetry and the smallest eigenvalue or a reliable Cholesky factorization.
• Step size, line-search rejections or trust-region ratios, and warnings about ill-conditioning or overflow in transformed parameters.
• Whether the stopping rule is relative to objective and parameter scales, and whether a different reasonable initialization reaches a compatible solution when the problem may have multiple local optima.

For a smooth problem with constraints, the gradient need not be zero at a constrained optimum. Instead, stationarity is expressed using feasible directions or, under appropriate regularity conditions, Karush–Kuhn–Tucker conditions. For example, if constraints are written \(g_i(x)\leq0\) and \(h_j(x)=0\), the KKT conditions include primal feasibility, nonnegative inequality multipliers, complementary slackness, and stationarity of the Lagrangian. A diagnostic that checks only the unconstrained gradient can therefore misclassify a valid boundary optimum.

Common Pitfalls

• Assuming an indefinite matrix makes every likelihood calculation undefined. An indefinite matrix may be invertible and may have positive determinant. It is still not a valid Gaussian covariance. Check the model constraint, not just whether an inverse happens to exist.
• Assuming Armijo puts an upper bound on the accepted step. It imposes a relative decrease requirement, not a universal maximum length. For \(f(x)=-x\), \(p=1\), and \(0<c<1\), the condition \(f(x+\alpha p)\leq f(x)+c\alpha f'(x)p\) reduces to \(-\alpha\leq-c\alpha\), which holds for every \(\alpha>0\). A particular backtracking algorithm's initial trial length \(\alpha_0\) supplies a cap on its tested steps.
• Treating a small step as proof of optimality. A small step can result from scaling, a restrictive trust region, or repeated line-search contraction. Check stationarity and feasibility as well as step size.
• Adding a covariance floor without recognizing its effect. A floor changes the feasible set and can change the estimate, as the diagonal example demonstrates. It should reflect the model or a stated regularization choice.