Numerical roundoff error in bilinear alternation leads to infeasibility

Question

We are trying an alteration optimization strategy to solve Lyapunov problems.

We break down our decision variables into two sets, Set 1 and Set 2.

We were perplexed how it was possible that, after getting a solution to Set 1, and plugging in those solved variables into the optimization over Set 2, the transferred variables would not be feasible.

The constraints that fail are those due to the SOS coefficient matching equality constraints.

Here, we print in each row the constraint that failed, and the value of our Initial Guess. We can see that the Initial Guess is off only a very small amount compared to the constraints.

LinearEqualityConstraint
(2 * Symmetric(97,40) + 2 * Symmetric(96,41)) == 9.50028
 [9.50027496]
LinearEqualityConstraint
(2 * Symmetric(97,47) + 2 * Symmetric(96,48)) == 234.465
 [234.4647013]
LinearEqualityConstraint
(2 * Symmetric(97,54) + 2 * Symmetric(96,55)) == -234.463
 [-234.46336504]
LinearEqualityConstraint
(2 * Symmetric(97,61) + 2 * Symmetric(96,62)) == 12.7962
 [12.79618825]
LinearEqualityConstraint
(2 * Symmetric(97,68) + 2 * Symmetric(96,69)) == -12.7964
 [-12.79637068]
LinearEqualityConstraint
(2 * Symmetric(97,75) + 2 * Symmetric(96,76)) == -51.4061
 [-51.40605828]
LinearEqualityConstraint
(2 * Symmetric(97,81) + 2 * Symmetric(96,82)) == 51.406
 [51.40604213]
LinearEqualityConstraint
(2 * Symmetric(97,86) + 2 * Symmetric(96,87)) == 192.794
 [192.79430158]
LinearEqualityConstraint
(2 * Symmetric(97,90) + 2 * Symmetric(96,91)) == -141.924
 [-141.92366183]
LinearEqualityConstraint
(2 * Symmetric(97,93) + 2 * Symmetric(96,94)) == -37.6674
 [-37.66740401]
InitialGuess V_sos and

Our guess for what's happening is:

1. When you extract the solution from one optimization using result.GetSolution(var), you lose some precision.

2. Or, when you set the previous solution using prog.SetInitialGuess(np_array) you lose some precision.

What's the solution here? Should we just keep feeding the solution back in even though it says infeasible?

Answer 1

This is a partial cookbook when I debug SOS problem, especially when working with Lyapunov problems:

---

Choose the right monomial basis

The main idea is to remove the 0-th order monomial 1 from the monomial basis of the sos polynomial. Here is a quick explanation:

The mathematical problem is

Find ?(x)
?Vdot ? ?(x) * (? ? V) is sos
 ?(x) is sos

Namely you want to prove that V?? ? Vdot ? 0

So first I would suggest to re-write your dynamics to make sure that 0 is the goal state (you can always shift your state).

Second you can see that since x=0 is the equilibrium point, then both V(0) = 0 and Vdot(0) = 0 (Because x=0 is the global minimal of V(x), hence ?V/?x=0 at x=0, indicating Vdot(0) = 0), now your sos polynomial p(x) = ?Vdot ? ?(x) * (? ? V) must satisfy p(0) = -?(0) * ?. But ?(x) >= 0 and ? > 0, so we know ?(0) = 0.

Lemma If a sos polynomial s(x) satisfies s(0) = 0, then its monomial basis cannot contain the 0-th order monomial (namely 1).

Proof Remember that s(x) is a sos polynomial, namely

s(x) = m(x)?Qm(x)

where m(x) contains the monomial basis, and Q is a psd matrix. Now let's decompose the monomial basis m(x) into two parts, the 0-th order monomial 1 and the remaining monomials mbar(x). For example, if m(x) = [x1, x2, 1], then mbar(x) = [x1, x2]. We also decompose the psd matrix Q accordingly

s(x) = [mbar(x)]? [Q11 Q10] [mbar(x)]
       [      1]  [Q10 Q00] [      1]

Since s(0) = Q00 = 0, we also know that Q10 = 0, so now we can use a smaller psd matrix Q11 rather than Q. Equivalently we write s(x) = mbar(x)? * Q11 * mbar(x), where mbar(x) is the monomial basis that doesn't contain the 0-th order monomial, QED.

So why removing the 0-th order monomial from the monomial basis and use the smaller psd matrix Q11 is a good idea when your sos polynomial s(x) satisfies s(0) = 0? The reason is that if the monomial basis contains 1, then your psd matrix Q has to be on the boundary of the psd cone, namely your SDP problem doesn't have a strict interior. This could leads to violation of Slater's condition, which also breaks the strong duality. One example is that if your s(x) = x², by including the 0'th order monomial, it is written as

x² = [x] [1 0] [x]
     [1] [0 0] [1]

And you see that the Gram matrix [[1 0], [0, 0]] is on the boundary of the psd cone (with one eigen value equal to 0). But if you remove 1 from the monomial basis, then its Gram matrix is just Q11=1, strictly in the interior of the psd cone.

In Drake, after removing 1 from the monomial basis, you can create your sos polynomial ?(x) as

lambda_poly, lambda_gram = prog.NewSosPolynomial(monomial_basis)

whee monomial_basis doesn't contain the 0-th order monomial.

---

Backoff during bilinear alternation

This is a typical problem in bilinear alternation. The issue is that when you solve a conic optimization problem with an objective function, the optimal solution always occurs at the boundary of the cone, namely it is very close to being infeasible. Then when you fix some variables to this solution at the cone boundary, in the next iteration the problem is very likely infeasible due to numerical roundoff error.

A typical solution is that after solving the optimization problem on variable Set 1 with an objective, now "backoff" a little bit by solving a feasibility problem on variable Set 1. This new solution is often strictly feasible (namely it is inside the strict interior of the cone), now pass this strictly feasible solution Set 1 to the next iteration and search for Set 2.

More concretely, suppose at one iteration you solve the following optimization problem

min c'*x
s.t constraint_on_x

and denote the optimal cost as p. Now solve a new feasibility problem

find x
s.t c'*x <= p + epsilon
    constraint_on_x

where epsilon can be a small positive number. This new solution will be used in the next iteration to search for a different set of variables.

You can check if your solution is on the boundary of the positive semidefinite cone by checking the Eigen value of your psd matrix. Here is the pseudo-code

for binding : prog.positive_semidefinite_constraints():
    psd_sol = result.GetSolution(binding.variables())
    psd_sol.reshape((binding.evaluator().matrix_rows(), binding.evaluator().matrix_rows()))
    print(f"minimal eigenvalue {np.linalg.eig(psd_sol)[0].min()}")

You should see that before doing this "backoff" some of the minimal eigen value is almost 0. After "backoff" the minimal eigen value gets larger.