# Optimization with matrix determinant as constraint

I'm solving a constrained optimization for matrix $\mathbf{A}$ with dimension 6x6, where one of the constraints is $\mathrm{det}(\mathbf{A})>0$. I use the NLopt package to solve my problem and provide value&&gradient routines of the objective and constraints. As described in the matrix cookbook, the gradient of matrix determinant is computed as $\frac{\partial \mathrm{det}(\mathbf{A})}{\partial \mathbf{A}} = \mathrm{det}(\mathbf{A})(\mathbf{A}^{-1})^T$ and involves matrix inverse. During the optimization iterations, one intermediate solution might violates the constraint and leads to singular matrix $\mathbf{A}$. The singular $\mathbf{A}$ will halt the optimization because evaluation of constraint gradient throws a "not invertible" exception.

One possible solution is to use optimization algorithms that do not need the gradient of constraints. But I wonder if there's any workaround if I use gradient-based algorithms. As far as I know, methods using gradient information would converge faster.

=============== Edit =================

The objective function is of the form:

$$\min_{\mathbf{A}}\quad\sum_i\frac{1}{2}(\mathbf{E}_i^T\mathbf{A}\mathbf{E}_i-t_i)^2$$

where $\mathbf{E}_i$ are known 6x1 vectors and $t_i$ are knwon scalars. Other constraints enforced on $\mathbf{A}$ are: (1) $\mathbf{A}$ is symmetric matrix (2) entries of $\mathbf{A}$ is non-negative

• Your formula ${\partial \mathrm{det}(\mathbf{A})}/{\partial \mathbf{A}} = \mathrm{det}(\mathbf{A})(\mathbf{A}^{-1})^T$ is numerically unstable when $\mathrm{det}(\mathbf{A})=0$: you are trying to express a finite quantity as $0\cdot\infty$, and this is not possible by naively computing the product of a small term times a very big (ill conditioned) one. Apart from the very pertinent answer by @ArnoldNeumaier, other computing strategies could be devised, if you provide some more background. How big is $\mathbf{A}$, are you trying a global or a local optimisation? – Stefano M Jun 4 '15 at 20:57
• Indeed the gradient is numerically unstable and that's why I'm seraching for a solution that avoids naively computing the unstable gradient. $\mathbf{A}$ is 6$\times$6 in our case, and the objective is quadratic with respect to $\mathbf{A}$. – Fei Zhu Jun 5 '15 at 2:24
• The constraint $\det A \gt 0$ is very mild. It would be helpful to know what objective function is being optimized and all the constraints so that suggestions can be more targeted. – hardmath Jun 5 '15 at 16:49
• How many $\mathbf{E_i}$ vectors are involved? – hardmath Jun 28 '15 at 4:22
• Generally no more than 20. – Fei Zhu Jun 30 '15 at 5:23

Replace $A$ by a factored form from which the determinant can be computed stably and cheaply. E.g., add new triangular variables $L$ and $R$ and the constraint $A=PLR$ for some fixed permutation matrix $P$; often the identity permutation $P=I$ is enough. This means that you treat the nontrivial entries of $L$ and $R$ as additional variables. Then determinants and their gradients are trivially to compute. If $A$ is known to be symmetric, you may use instead $A=LDL^T$.

Alternatively you may eliminate $A$ by substituting $PLR$ for it everywhere it occurs. Thus in place of minimizing $f(A,\det A)$ subject to $F(A)=0$, say, you would minimize $f(LR,\det R)$ subject to $F(LR)=0$, (assuming for simplicity that $L$ is unit lower triangular and no permutation was applied), and use $\det R=\prod R_{ii}$.

In case the solution has no such factorization, you need to change $P$ after a preliminary round of minimization with $P=I$.

• Embarrassed to admit that I didn't fully understand your answer. My background is computer engineering and I only scratch the surface of numerical computing. If I'm interpreting it correctly, did you mean solving for $\mathbf{L}$ and $\mathbf{R}$ instead of solving for $\mathbf{A}$? The determinants of the factored matrices may be trivially to compute, while would the gradients still involve matrix inverse? Is this solution a commonplace strategy? If so, could you please provide some references (e.g. papers) so that I can figure it out by reading the papers. – Fei Zhu Jun 5 '15 at 2:40

The constraint $\det(A) > 0$ is very unstable and absolutely horrible. Would you believe that on practical, non-contrived 6 by 6 symmetric matrices, changing one element of the matrix by 1e-15 can make the determinant computed in double precision swing from 1e20 to -1e20? Good luck in trying to have a meaningful determination of whether $\det(A)$ is > or < 0 at any point. The gradient is the least of your worries. The constraint itself just will not work unless you compute in ultra-high precision, or perhaps you are only dealing with a 2 or 3 dimensional problem. Prof Neumaier was giving you a more numerically stable workaround. But I can say from my experience, imposing a Cholesky factorization type constraint does not seem to work well, and makes the problem non-convex, even if the original problem is convex; you may succeed, but you may wind up adding a bunch of spurious local optima (and even saddlepoints?) to your problem, and that is not a good thing.

Is the matrix $A$ (constrained to be) symmetric, and do you really wish to impose that $\det(A) > 0$ and that all leading principal submatrices have $\det > 0$, i.e., that $A$ be positive definite? If so, instead of using the constraint $\det(A) > 0$ (which apart from its terrible numerical properties, is not sufficient to ensure positive definiteness of $A$, unless it is imposed (as perhaps you do?) on all leading principal submatrices), use a semidefinite constraint, $A > 0$ in Linear Matrix Inequality (LMI) sense. I.e., a constraint that A is positive definite. However, practical optimization software either does not accept strict inequality constraints, or if it does, treats them as non-strict inequality constraints. If positive semidefinite is o.k., your constraint would be $A \ge 0$ (in LMI sense). If you really need A to be strictly positive definite (det(A) strictly positive), then you need to determine a minimum acceptable eigenvalue, min_eig, and then impose the constraint A - min_eig * (Identity matrix of correct size) >= 0 (in LMI sense). You do need to make some allowance for solver tolerance used in constraint feasibility assessment (i.e., solver will allow constraint to be violated by a little, for example, by 1e-6 or 1e-8).

At this point, it would help to know what your objective function and the rest of the constraints are. If the objective function and the rest of the constraints are convex (presuming a minimization problem), you could use CVX or YALMIP (both are free) to solve your problem. You said in a comment (belongs in the question, so you should edit) that the objective is quadratic, but is it convex? If non-convex, you may still be able to use YALMIP, perhaps in conjunction with having YALMIP call the free solver PENLAB. If you use CVX or YALMIP, they take care of any differentiation for you.

Edit in response to edit of Question: As shown in the edited version of your question, $A$ is constrained to be symmetric nonnegative and have $\det(A) > 0$, but $A$ is not constrained to be positive definite. It is in fact possible for an nonnegative symmetric matrix of dimension >= 3 to have positive determinant yet not be positive definite. Do you wish to allow that situation? If so, my answer above the edit above does not apply. If in fact you require $A$ be positive definite, then my answer above the edit does apply, and you should further edit your question to reflect the constraint that $A$ be positive definite.

• It would be great if you provided an example in your answer, or showed a summary of the analysis that leads to an estimate of the condition number of the determinant. – David Ketcheson Jun 28 '15 at 5:55
• I don't have the data any more. But there were several matrix entries which were optimization variables, and the optimization drove the matrix to near singularity, so very high condition number. That resulted in the huge oscillations in computed determinant. which made the use of det(matrix) >= 0 constraint useless. This was just an experiment on my part to see whether it was viable -I knew it might not be.There are many optimization problems with a det or psd constraint for which the optimum will be right at the boundary, so this situation occurs quite frequently, and apparently with the OP. – Mark L. Stone Jun 28 '15 at 16:10