# How to efficiently optimize a strictly increasing linear function over convex set

I am trying to minimize the trace of a symmetric matrix $B$ fulfilling $x^T A_i x \leq x^T B x \ \forall x \forall i$ for given (in general) indefinite symmetric matrices $A_i, i = 1, ..., n$.

The straight foward formulation is in d dimensional real space with n matrices would be

d = 2;
n = 3;
As = Normal /@ Symmetrize /@ RandomInteger[{-10, 10}, {n, d, d}];
B = Normal@
SymmetrizedArray[
pos_ :> Subscript[b, pos], {d, d}, {Symmetric[{1, 2}]}];
vars = Variables@B;
f = Tr@B;
x = Array[xc, d];
cond1 = Simplify@And @@
Flatten[Resolve[ForAll[Evaluate@x, 0 <= x.(B - #).x], Reals] & /@
As];
ips = vars /. FindInstance[cond1, vars, Reals, 4];
min1 = NMinimize[{f, cond1}, vars,
Method -> {"Automatic", "InitialPoints" -> ips}]; // AbsoluteTiming
cond1 /. min1[[2]]


False

which sadly usually gives me false results at the end.

I know that the set described by the conditions $x^T A_i x \leq x^T B x \ \forall x \forall i$ is a convex set and the function f = Tr@B is naturally linear (strictly increasing) in B. I want to treat this problem in slightly higher dimensions (d = 7) where the number of conditions is very annoying. Equivalent conditions for $x^T A_i x \leq x^T B x \ \forall x\forall i$ are obtained by consideration of all principal minors of $(B-A_i)$ to be non-negative, see Sylvester's criterion. I tried to treat this problem with those but the numerical solutions still need quite a lot of time (several hours).

My question: do you see a more efficient formulation of this problem or do you know a way to solve this numerically in an time efficient manner?

I am not sure, and I dont know a lot of the subject, but the linearity of the (strictly increasing) objective function and the convexity of the underlying set should be pretty nice properties to work with in efficient algorithms, despite the high number of connstraints describing the convex set created by the positive semidefinite $B-A_i$.

EDIT 1: Using FindMinimum and minors

Based on Henrik's comment, I tried to use FindMinimum and reformulated the conditions with the minors of the respective matrices

d = 3;
n = 3;
As = Normal /@ Symmetrize /@ RandomInteger[{-10, 10}, {n, d, d}];
B = Normal@
SymmetrizedArray[
pos_ :> Subscript[b, pos], {d, d}, {Symmetric[{1, 2}]}];
vars = Variables@B;
f = Tr@B;
cond2 = # >= 0 & /@
Flatten[Table[Minors[B - #, k], {k, Length@B}] & /@ As];
cond2 = And @@ cond2;
cond2 = Simplify@cond2;
min2 = FindMinimum[{f, cond2}, vars,
Method -> "InteriorPoint"]; // AbsoluteTiming
cond2 /. min2[[2]]
x = Array[xc, d];
Bsol = B /. min2[[2]];
And @@ Flatten[
Resolve[ForAll[Evaluate@x, 0 <= x.(Bsol - #).x], Reals] & /@ As]


Sadly, I was wrong, this approach works only for some cases. Sometimes you get a True solution, sometimes you dont :D. I should not be working on this, I am at a birthday, but I can not leave this alone, it's fun :D! Further suggestions?

EDIT 2:

Sorry, my former conditions for the principal minors of the difference matrices $B-A_i$ were wrong, in the former code I took al minors. The present code selects all principal minors. I also tried Henrik's approach (see code below formulation with minors). None of both momentary implementations will give you True solutions for arbitrary cases. But, as you can see in the eigenvalues at the end of the (*Check solutions*) part, they do deliver acceptable results. At least in my computer, the minors formulation (around 13 second) seems slower than the slack variables approach (around 0.2 seconds) for $d=5$ and $n=3$. Although the minimum eigenvalues of $B-A_i$ are sometimes negative for the slack variables approach proposed by Henrik (such that $x^T A_i x \leq x^T B x \ \forall x \forall i$ does NOT hold), some results are at least for me very useful, since these negative eigenvalues are almost negligible for my applications. Thank you Henrik!

(*Input*)
d = 5;
n = 3;
As = Normal /@ Symmetrize /@ RandomReal[{-1, 1}, {n, d, d}];
B = Normal@
SymmetrizedArray[
pos_ :> Subscript[b, pos], {d, d}, {Symmetric[{1, 2}]}];
f = Tr@B;
(*Solution based on principal minors*)
PrincipalMinors[A_] := Diagonal@Minors[A, #] & /@ Range@Length@A;
vars = Variables@B;
MinorsVec = Flatten[PrincipalMinors[B - #] & /@ As];
psd = And @@ Thread[MinorsVec >= 0];
Print["Begin min with Minors"]
minMinors =
FindMinimum[{f, psd}, vars,
Method -> "InteriorPoint"]; // AbsoluteTiming
BMinors = B /. minMinors[[2]];
(*Solution based on slack variables*)
S = UpperTriangularize /@ Array[s, {n, d, d}];
vars = Join[Variables@B, Variables@S];
vars = {#, 1.} & /@ vars;
eqv = DeleteCases[
Flatten[UpperTriangularize /@ (ConstantArray[B, n] -
As - (Transpose[#].# & /@ S))], 0.];
Print["Begin min with slack variables, see Henrik's post"]
minSlack = FindMinimum[{f, eqs}, vars]; // AbsoluteTiming
BSlack = B /. minSlack[[2]];
(*Check solutions*)
Print["Check minimial eigenvalues of difference matrices:"]
ScientificForm@Min@Eigenvalues[BMinors - #] & /@ As
ScientificForm@Min@Eigenvalues[BSlack - #] & /@ As


• Are the $A_i$ positive definite? Commented Jul 28, 2017 at 22:34
• No, sadly I need to deal with even indefinite As Commented Jul 28, 2017 at 22:35
• Well, this is a semidefinite program. See here for more info. I think, an interior point method should work well, here... Commented Jul 28, 2017 at 22:44
• I wonder if the following will work. Diagonalize the A matrices and impose the conditions on B that for all (normalized) eigenvectors xvec[i,j] of each A[i], the inequality x[i,j].(A[i]-B).x[i,j]<=0. That becomes a finite set of constraints. That said, the reference and also the approach by @HenrikSchumacher look good (and got my vote). Commented Jul 29, 2017 at 15:12
• No, it's not as easy as I had hoped. Nor should it be I guess; if semidefinite programming could be set up that way people would be doing it. What I proposed were necessary but not sufficient conditions. Commented Jul 31, 2017 at 18:34

Besides using interior point methods, it is possible to translate this semidefinite problem to an equality constrained problem by introducing "squared" slack variables:

The condition $B-A_i \succeq 0$ can be reformulated as $B-A_i - S_i^\intercal S_i = 0$, where $S_i$ is an upper triangular matrix. Since this is an equation in symmetric matrices, we only have to consider the upper triangular part, yielding $\frac{d(d+1)}{2}$ equations. You will end up in an optimization problem in $(n+1) \frac{d(d+1)}{2}$ variables (i.e. $B$ and the $n$ slack variables $S_1,\dotsc,S_n$) with $n \frac{d(d+1)}{2}$ quadratic constraints. That should be easily attackable by NMinimize or FindMinimum.

In general, squared slack variables have an infamous reputation, so use them with care. For example, methods that primarily search for critical points tend to become lost, as squared slack variables always introduce a plethora of new critical points. See also here for a discussion of squared slack variables for semidefinite programming. (Note that the authors introduce superfluous degrees of freedom by not considering upper triangular matrices, only.)

Moreover, it is quite likely that many methods will have problems since in your case, the objective function is linear so that the constrained Hessian might not be positive definite. Thus, interior point methods (like $\log$-barrier method) might be indeed a better choice.

Edit: Okay, I tried to figure out the details for the slack variables in the meantime. The code often works fine for dimension d=3. This is just the skeleton and may quite likely be improved.

d = 3;
n = 7;
ndofs = d (d + 1)/2;
nslack = n d (d + 1)/2;

A = RandomReal[{-1, 1}, {n, d, d}];
A = A + Transpose /@ A;


Introducing variables x[i] and slack variables z[i].

xx = Table[x[i], {i, 1, ndofs}];
zz = Table[z[i], {i, 1, nslack}];


Producing the the corresponding matrices.

c = 0;
P = With[{pat = DeleteDuplicates[
Flatten[Table[
Table[c++; {{i, j, c}, {j, i, c}}, {j, i, d}], {i, 1, d}], 2]]
},
Normal[SparseArray[pat -> 1., {d, d, ndofs}, 0.]]
];
X = P.xx;
Z = UpperTriangularize /@ Transpose[P.Partition[zz, ndofs]\[Transpose], {2, 3, 1}];


Setting up the optimization problem for FindMinimum:

objective = Tr[X];
constraints =
DeleteCases[
Flatten[
UpperTriangularize /@ (ConstantArray[X, {n}] - A -
Map[Y \[Function] Y\[Transpose].Y, Z])],
0.
];


Just some starting values; not always a good choice.

x0 = 0. xx + 1.;
z0 = 0. zz + 1.;
vars = Join[Transpose[{xx, x0}], Transpose[{zz, z0}]];
sol = FindMinimum[{objective, Thread[constraints == 0.]}, vars]

Norm[constraints /. sol[[2]], \[Infinity]]
Eigenvalues /@ (ConstantArray[X, n] - A /. sol[[2]])


For higher values of d, the default method of FindMinimum does not converge. But for a start, I am quite satisfied.

• I edited my question. Sadly I was wrong, the approach I took works only sometimes. I still have to try yours with the slack variables. Can you provide an example? Commented Jul 29, 2017 at 18:09
• @MauricioLobos: Yes, I tried. Please, have a look. Commented Jul 29, 2017 at 20:07
• Thank you very much! Nice approach! I will ask a friend, if the Gershgorin circle theorem or diagonal dominant conditions help to simplify the problem formulation. Commented Jul 31, 2017 at 13:48
• @MauricioLobos: You're welcome! And let me know about your progress... Commented Jul 31, 2017 at 17:24

This post is rather informative. Henrik's approach is the best one in my opinion so far.

1) It is also possible to formulate the problem based on the Gerschgorin circle theorem. Based on this theorem, I consider the minimum of all real valued Gerschgorin intervals (relevant for real valued symmetric matrices of this problem) of each difference matrix $B-A_i$ and impose that this minimum $m_i$ is nonnegative, $m_i \geq 0 \forall i$. This formulation works with $n$ conditions for the objective function. The problem is that these conditions are very hard to fulfill.

2) Another alternative is to formualte the problem based on diagonal dominance. This approach works with $2dn$ conditions.

I just paste the code which immedeately compares my first approach, the approach of Henrik, the formulation based on Gerschgorin circles and diagonal dominance for real valued symmetric matrices. In the last lines showing the eigenvalues of the difference matrices $B-A_i$ you can see, that the Gerschgorin formulation delivers actual bounds for some cases. This is not generally the case in my computer since Mathematica has problems finding points fulfilling the imposed conditions. The formulation based on diagonal dominance seems problematic for Mathematica. I tried generating initional points with FindInstance but this also takes a loooooot of time.

(*Input*)
d = 5;
n = 3;
As = Normal /@ Symmetrize /@ RandomReal[{-1, 1}, {n, d, d}];
B = Normal@
SymmetrizedArray[
pos_ :> Subscript[b, pos], {d, d}, {Symmetric[{1, 2}]}];
f = Tr@B;
(*Solution based on principal minors*)
PrincipalMinors[A_] := Diagonal@Minors[A, #] & /@ Range@Length@A;
vars = Variables@B;
MinorsVec = Flatten[PrincipalMinors[B - #] & /@ As];
psd = And @@ Thread[MinorsVec >= 0];
Print["Begin min with Minors"]
minMinors =
FindMinimum[{f, psd}, vars,
Method -> "InteriorPoint"]; // AbsoluteTiming
BMinors = B /. minMinors[[2]];
(*Solution based on slack variables*)
S = UpperTriangularize /@ Array[s, {n, d, d}];
vars = Join[Variables@B, Variables@S];
vars = {#, 1.} & /@ vars;
eqv = DeleteCases[
Flatten[UpperTriangularize /@ (ConstantArray[B, n] -
As - (Transpose[#].# & /@ S))], 0.];
Print["Begin min with slack variables, see Henrik's post"]
minSlack = FindMinimum[{f, eqs}, vars]; // AbsoluteTiming
BSlack = B /. minSlack[[2]];
(*Solution based on Gerschgorin circles*)
GerschgorinInt[A_] :=
Table[{A[[i, i]],
Sum[Abs@A[[i, j]], {j, DeleteCases[Range@Length@A, i]}]}, {i,
Length@A}];
GerschgorinIntMin[A_] := #[[1]] - #[[2]] & /@ GerschgorinInt@A;
condGersch = Min@GerschgorinIntMin@(B - #) & /@ As;
vars = Variables@B;
Print["Begin min Gersch"]
minGersch = NMinimize[{f, condGersch}, vars]; // AbsoluteTiming
BGersch = B /. minGersch[[2]];
(*Solution based on diagonal dominance*)
diff = ConstantArray[B, n] - As;
DiagonalDom[A_] :=
Table[{Abs@A[[i, i]],
Sum[Abs@A[[i, j]], {j, DeleteCases[Range@Length@A, i]}]}, {i,
Length@A}];
DiagonalDomCond[A_] := 0 <= #[[2]] <= #[[1]] & /@ DiagonalDom@A;
condDD = Flatten[DiagonalDomCond /@ diff];
Print["Begin min DD"]
minDD = FindMinimum[{f, condDD}, vars]; // AbsoluteTiming
BDD = B /. minDD[[2]];
(*Check solutions*)
Print["Check minimial eigenvalues of difference matrices:"]
ScientificForm@Min@Eigenvalues[BMinors - #] & /@ As
ScientificForm@Min@Eigenvalues[BSlack - #] & /@ As
ScientificForm@Min@Eigenvalues[BGersch - #] & /@ As
ScientificForm@Min@Eigenvalues[BDD - #] & /@ As
`