Optimization with matrix determinant as constraint - Computational Science Stack Exchange most recent 30 from scicomp.stackexchange.com 2022-01-27T23:05:05Z https://scicomp.stackexchange.com/feeds/question/19840 https://creativecommons.org/licenses/by-sa/4.0/rdf https://scicomp.stackexchange.com/q/19840 4 Optimization with matrix determinant as constraint Fei Zhu https://scicomp.stackexchange.com/users/1976 2015-06-04T08:34:11Z 2018-03-19T21:03:41Z <p>I'm solving a constrained optimization for matrix $\mathbf{A}$ with dimension 6x6, where one of the constraints is $\mathrm{det}(\mathbf{A})&gt;0$. I use the <a href="http://ab-initio.mit.edu/wiki/index.php/NLopt" rel="nofollow">NLopt</a> package to solve my problem and provide value&amp;&amp;gradient routines of the objective and constraints. As described in <a href="http://www.eattardo.com/wp-content/uploads/2013/10/MatrixCookbook2012.pdf" rel="nofollow">the matrix cookbook</a>, 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.</p> <p>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.</p> <p>=============== Edit =================</p> <p>The objective function is of the form:</p> <p>$$\min_{\mathbf{A}}\quad\sum_i\frac{1}{2}(\mathbf{E}_i^T\mathbf{A}\mathbf{E}_i-t_i)^2$$ </p> <p>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</p> https://scicomp.stackexchange.com/questions/19840/-/19846#19846 4 Answer by Arnold Neumaier for Optimization with matrix determinant as constraint Arnold Neumaier https://scicomp.stackexchange.com/users/1117 2015-06-04T15:46:40Z 2015-06-05T09:47:31Z <p>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$.</p> <p>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}$.</p> <p>In case the solution has no such factorization, you need to change $P$ after a preliminary round of minimization with $P=I$.</p> https://scicomp.stackexchange.com/questions/19840/-/20008#20008 3 Answer by Mark L. Stone for Optimization with matrix determinant as constraint Mark L. Stone https://scicomp.stackexchange.com/users/16756 2015-06-26T14:10:21Z 2018-03-19T21:03:41Z <p>The constraint $\det(A) &gt; 0$ is <strong>very unstable</strong> and <strong>absolutely horrible</strong>. 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 &lt; 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.</p> <p>Is the matrix $A$ (constrained to be) symmetric, and do you really wish to impose that $\det(A) &gt; 0$ and that all leading principal submatrices have $\det &gt; 0$, i.e., that $A$ be positive definite? If so, instead of using the constraint $\det(A) &gt; 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 &gt; 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, <code>min_eig</code>, and then impose the constraint <code>A - min_eig * (Identity matrix of correct size) &gt;= 0</code> (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).</p> <p>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.</p> <p><strong>Edit in response to edit of Question</strong>: As shown in the edited version of your question, $A$ is constrained to be symmetric nonnegative and have $\det(A) &gt; 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. </p>