# Lyapunov + Symmetries = Infeasible, Lyapunov - Symmetries = Feasible

Hi again @blegat @abulak

I am using SOS programming to find a Lyapunov function, obeying some natural symmetry. The time derivative of V also obeys said symmetry. The feasibility of the Lyapunov function is dependent on the largeness of a parameter specified at the beginning of the problem which determines the dynamics of the differential equation. I will outline a sort of pseudocode here.

Truncated system:

``````@polyvar a[1:6]
Vbasis = (some vector of symmetry-invariant polynomial terms in a)
V_coeffs = @variable(model, [1:length(Vbasis)])
V = sum(V_coeffs.*Vbasis)
dVdt = f.*differentiate(V,a)
model = SOSModel(Mosek.Optimizer)
con_Vsos = @constraint(model, V - epsilon*dot(a,a) >= 0, symmetry=pattern)
minusdVdtsos = -dVdt - epsilon*dot(a,a)
con_dv = @constraint(model, minusdVdtsos >= 0, symmetry=pattern)
optimize!(model)
``````

This works beautifully, the symmetry is imposed and outperforms the setting where we enforce sparsity=SignSymmetry() instead. Now, I add another term to the Lyapunov function V for which I only have estimates on dV/dt.

``````@polyvar a[1:6], q
Vbasis = (some vector of symmetry-invariant polynomial terms in a, q^2)
V_coeffs = @variable(model, [1:length(Vbasis)])
V = V_coeffs.*Vbasis

sbasis = (some vector of polynomial terms in a, q^2)
s_coeffs = @variable(model, [1:length(sbasis)])
s = sum(s_coeffs.*sbasis
bound = constant*s*q^2 - (positive term)

Mis =  differentiate(V,a) - 2*dV/d(q^2)*a

model = SOSModel(Mosek.Optimizer)

# These constraints enforce the fact that |Mis| <= s; we do this for the length of Mis (6) but I only put in one constraint for simplicity
@constraint(model, s + Mis>=0)
@constraint(model, s - Mis >=0)

dVdt = f.*differentiate(V,a) + bound
con_Vsos = @constraint(model, V - epsilon*dot([a;q],[a;q]) >= 0, symmetry=pattern)
minusdVdtsos = -dVdt - epsilon*dot([a;q],[a;q])
con_dv = @constraint(model, minusdVdtsos >= 0, symmetry=pattern)
optimize!(model)
``````

If I do not enforce the natural symmetry on V, the program runs beautifully, and indeed the threshold parameter I am using to determine feasibility is equivalent between this Julia program and the old Matlab program I am running. Moreover, if I enforce instead sparsity=SignSymmetry() on any one of the constraints (turning them on and off like circuit breakers), it also preserves the feasibility result below the threshold parameter. However, the minute I enforce the natural symmetry, Mosek returns an infeasibility result unless my parameter is way below the threshold parameter. There are similar symmetries that can also be imposed on the bounds, but those also behave badly.

Here are multiple checks that have been made so far:

1. The natural symmetry with q is the same as the natural symmetry without q, we do not enforce any symmetries on q so the pattern generation is identical except for the q term. I have verified that
- the homomorphism is indeed a homomorphism
- both V and minusdVdtsos are invariant under the defined homomorphism
2. I have ensured that any coefficients on any variables for the program are all above a certain threshold tolerance so there is no extra noise being sent to the Mosek solver. (There was actually a non-trivial observation here that caused this to be not as straightforward as I expected, and I may post about it in a different thread to ask for an error request for other users at some point, hopefully it’s a useful observation.)
3. The size of the bases for the constraints in both the unsymmetrized (monomial basis) and symmetrized (nonmonomial basis) problems are of the same length, so the problem has many fewer constraints and variables that the original non-symmetrized problem. (Easy to be seen in the Mosek output as well.)
4. I have minimized the number of free variables for the SOS program; for example, if V(a, q) = c_1(a1^2 + a2^2) + c2q^2, the free variables are also written this way.

My best intuition is there is some way in which JuMP is formulating the problem to hand to Mosek that causes Mosek to behave much worse numerically than it should. Supporting this intuition, I will sometimes get warnings from Mosek that the matrix ‘A’ is sparse. Indeed, this error also sometimes ocurs when I enforce the restriction sparsity=SignSymmetry(), but it does not strongly hamper the feasibility results like my handwritten symmetry. I was concerned about the basis size for the original problem vs the symmetrized problem, if one was longer than the other, but they were the same size so I was not worried about the fact that the block symmetry does not enforce Newton polytope reduction. Any time you have to peek at this would be greatly appreciated. Thanks.

2 Likes

Here is a minimalish reproduction of this problem, or at least a very similar one:

``````using SumOfSquares
using DynamicPolynomials
using LinearAlgebra
using MosekTools

sparsity = Sparsity.SignSymmetry()
@polyvar x[1:3] z[1:3]

# Some ODE system with symmetries
f = -[x+x*x, x, x]
Df = differentiate(f,x)

# Build bases for our two unknown polynomials respecting the symmetries
basisV = [x; x*x; x*x; x*x; x*x]
basisrho = [1; basisV]

model = SOSModel(Mosek.Optimizer)

B = @variable(model)

# Lyapunov function to find
Vc = @variable(model,[1:length(basisV)])
V = dot(Vc,basisV)
dVdx = differentiate(V,x)

# Second unknown polynomial
rhoc = @variable(model,[1:length(basisrho)])
rho = dot(rhoc,basisrho)

@constraint(model, B - dot(z,Df*z) - dot(f, dVdx) - rho*(1-dot(z,z)) >= 0, sparsity=sparsity)
@objective(model, Min, B)

optimize!(model)
display(solution_summary(model))

# Should be 1
display(value(B))
``````

The problem is reported as infeasible, but changing to Sparsity.NoPattern(), the problem is easily solved. This issue seems to exist regardless of whether I use Mosek or some other solvers I have tried.

Is there a bug? Or am I doing something wrong? This same problem recreated with YALMIP works fine with Mosek, successfully exploiting the symmetries.

2 Likes

If look at model.moi_backend.optimizer.model.task I get

``````    Minimize
#obj: + #x1
Subject to
#c1: + #x1 - #x7 + #MX1 ⋅ #X̄4 = -0.0
#c2: + #MX2 ⋅ #X̄4 = 0.0
#c3: + #x2 - #x8 + #MX3 ⋅ #X̄4 = -0.0
#c4: + #x7 + #MX4 ⋅ #X̄4 = -0.0
#c5: = -2.0
#c6: + #x7 + #MX5 ⋅ #X̄3 = -0.0
#c7: + #MX6 ⋅ #X̄3 = 0.0
#c8: + #x7 + #MX7 ⋅ #X̄3 = -1.0
#c9: + #x5 - #x12 + #MX8 ⋅ #X̄2 = -0.0
#c10: + #x2 + 2.0 #x4 + 2.0 #x6 - #x11 + #MX9 ⋅ #X̄2 = -0.0
#c11: + #x5 - #x10 + #MX10 ⋅ #X̄2 = -0.0
#c12: + #MX11 ⋅ #X̄4 = 0.0
#c13: + 2.0 #x3 - #x9 + #MX12 ⋅ #X̄4 = -0.0
#c14: = -1.0
#c15: + #MX13 ⋅ #X̄2 = 0.0
#c16: = -1.0
#c17: + #MX14 ⋅ #X̄2 = 0.0
#c18: + #MX15 ⋅ #X̄2 = 0.0
#c19: + #x8 + #MX16 ⋅ #X̄4 = -0.0
#c20: + #x8 + #MX17 ⋅ #X̄3 = -0.0
#c21: + #MX18 ⋅ #X̄3 = 0.0
#c22: + #x8 + #MX19 ⋅ #X̄3 = -0.0
#c23: + 2.0 #x3 = -0.0
#c24: + #MX20 ⋅ #X̄4 = 0.0
#c25: + #x12 + #MX21 ⋅ #X̄2 = -0.0
#c26: + #x12 + #MX22 ⋅ #X̄1 = -0.0
#c27: + #MX23 ⋅ #X̄1 = 0.0
#c28: + #x12 + #MX24 ⋅ #X̄1 = -0.0
#c29: + #x11 + #MX25 ⋅ #X̄2 = -0.0
#c30: + #x11 + #MX26 ⋅ #X̄1 = -0.0
#c31: + #MX27 ⋅ #X̄1 = 0.0
#c32: + #x11 + #MX28 ⋅ #X̄1 = -0.0
#c33: + #x10 + #MX29 ⋅ #X̄2 = -0.0
#c34: + #x10 + #MX30 ⋅ #X̄1 = -0.0
#c35: + #MX31 ⋅ #X̄1 = 0.0
#c36: + #x10 + #MX32 ⋅ #X̄1 = -0.0
#c37: + #x9 + #MX33 ⋅ #X̄4 = -0.0
#c38: + #x9 + #MX34 ⋅ #X̄3 = -0.0
#c39: + #MX35 ⋅ #X̄3 = 0.0
#c40: + #x9 + #MX36 ⋅ #X̄3 = -0.0
``````

c5, c14 and c16 are particularly suspicious.

If anyone can help with this it would be greatly appreciated

1 Like

I cannot offer much help because I am not familiar with SOS optimization. This is pretty much just running different options from the documentation to see what happens.

I did make the example run with `sparsity = Sparsity.SignSymmetry()` by changing the basis:

``````# basisV = [x; x*x; x*x; x*x; x*x] ####Default, yields infeasible
basisV = [x,x*x,x*x, x*x, x*x,x*z,x*z] ####Works
``````

I also wanted to see the decomposition for sign sparsity and it turns out that changing the basis led to no decomposition (new basis+SignSparsity()):

It seems that with the new basis, the “decomposition” is the same as when NoPattern() is used (default basis+NoPattern()):

And this is what happens for the default basis and SignSparsity():

I don’t think we should be getting zero matrices in the decomposition, right? So maybe either sign sparsity cannot be enforced here, or something is off in the decomposition?

Other types of decomposition work OK for the default basis, for instance:

``````sparsity = Sparsity.Monomial(ChordalCompletion())
``````

gives:

The code for the sake of completeness:

Summary
``````using SumOfSquares
using DynamicPolynomials
using LinearAlgebra
using COSMO ####Same for CSDP solver

# sparsity = Sparsity.NoPattern()
sparsity = Sparsity.SignSymmetry()
#sparsity = Sparsity.Monomial(ChordalCompletion())

@polyvar x[1:3] z[1:3]

# Some ODE system with symmetries
f = -[x+x*x, x, x]
Df = differentiate(f,x)

# Build bases for our two unknown polynomials respecting the symmetries
# basisV = [x; x*x; x*x; x*x; x*x]
basisV = [x,x*x,x*x, x*x, x*x,x*z,x*z]

basisrho = [1; basisV]

model = SOSModel(COSMO.Optimizer)

B = @variable(model)

# Lyapunov function to find
Vc = @variable(model,[1:length(basisV)])
V = dot(Vc,basisV)
dVdx = differentiate(V,x)

# Second unknown polynomial
rhoc = @variable(model,[1:length(basisrho)])
rho = dot(rhoc,basisrho)

cst = @constraint(model, B - dot(z,Df*z) - dot(f, dVdx) - rho*(1-dot(z,z)) >= 0.0, sparsity=sparsity)
@objective(model, Min, B)

optimize!(model)
display(solution_summary(model))

# Should be 1
display(value(B))
``````
1 Like

I have ensured that any coefficients on any variables for the program are all above a certain threshold tolerance so there is no extra noise being sent to the Mosek solver. (There was actually a non-trivial observation here that caused this to be not as straightforward as I expected, and I may post about it in a different thread to ask for an error request for other users at some point, hopefully it’s a useful observation.)

Yes please do ! I think we should have a presolve layer doing these kind of things. I currently considering reviving MathOptPresolve for that.

Thanks for the MRE, it seems we’re close, I’ll take a look

2 Likes

This one should be fixed by Fix and document XORSpace by blegat · Pull Request #314 · jump-dev/SumOfSquares.jl · GitHub, thanks for reporting it! EDIT: it should be fixed in SumOfSquares v0.7.1

1 Like

Ok I will post that in a separate thread!

The new SumOfSquares version works wonderfully for @jezzaparker’s minimal example. Unfortunately, his problem is unrelated to the OG post at this juncture, as the resulting type of the way I compiled my bases is the same as that recommended by @blob.

I have compared all relevant outputs for the problem that is handed over to Mosek for my particular system, and unfortunately have not found anything seemingly suspicious (though this could be my untrained eye). The last thing I can think at this point is, V obeys a certain symmetry, and dVdt obeys that same symmetry, but if you compute this explicitly you see that

``````V - SymbolicWedderburn.action(act(), GroupElement(1,G), V)
Output: 0

negdVdt - SymbolicWedderburn.action(act(), GroupElement(1,G), negdVdt)
Output: # polynomial with coefficient values of 0 or of order ~10^-16, so machine 0
``````

This is not a problem in my truncated version, but is it possible this could be problematic in the full version? Especially since the bounds are added to negdVdt and obey their own set of symmetries? (When added together the bounds still allow negdVdt to obey the original symmetry per the above calculation as well.). I’m not sure this is “the problem”, as I can choose to impose symmetry only on V (neglecting all symmetries on every other constraint) and it still returns that the problem is infeasible (even though it is feasible without imposing said symmetry).

In order to dig into this a little more, I believe the only way to make progress via forum would be through another minimal example, but unfortunately this would have to be for a completely different system from the problem I am considering. Some of the minimal examples I have written do not reproduce the error, but it is possible they are too minimal. It would probably be most effective to either converse with someone familiar with this via email or a video chat where we can look over the code privately. In the meantime, if anyone is aware of any literature imposing multiple different types of symmetries for different constraints in Julia, that would be greatly appreciated.

Hi @ilovemath, any update on this ? Did you check the group correctness with conformance tests (see rewrite the definitions of groups in symmetry.jl by kalmarek · Pull Request #32 · blegat/CondensedMatterSOS.jl · GitHub). I just found a bug that I fixed in Fix simultaneous block diagonalization for Complex by blegat · Pull Request #320 · jump-dev/SumOfSquares.jl · GitHub, is this helping ?

2 Likes

Hi @blegat ,

I have updated JuMP, SumOfSquares, Mosek, and MosekTools, and re-ran the code (after restarting Julia of course), and unfortunately the issue persists. And yes, I have double checked my group structures and generators, they are good, and match up with automatically generated invariant polynomials. The symmetry enforcement works for case where I have no SOS inequality constraints. But, darn, it appears I didn’t state something very important as clearly as I should have. The problem’s feasibility is dependent on a parameter, R. This parameter R tunes the coefficient values for the f in `dVdt = f.*differentiate(V,a) + (bounds)`. For the full system (or even the system with sign symmetries implemented), there is a critical R for which Mosek claims there is a feasible polynomial. If we enforce the stronger symmetry for this R value, Mosek claims the problem is infeasible (while sometimes giving errors that there are nearly zero elements in sparse matrices handed to Mosek), which we know is not true because it is feasible without that symmetry enforcement. If I choose R significantly below the critical R value, feasibility holds. As I said in my original post,

My best intuition is there is some way in which JuMP is formulating the problem to hand to Mosek that causes Mosek to behave much worse numerically than it should. Supporting this intuition, I will sometimes get warnings from Mosek that the matrix ‘A’ is sparse (I meant to say here, that it has near-zero entries in sparse matrices). Indeed, this error also sometimes occurs when I enforce the restriction sparsity=SignSymmetry(), but it does not strongly hamper the feasibility results like my handwritten symmetry.

This near-zero matrix error that JuMP sends me occurs for the enforcement of sign symmetry, but this is not the case for YALMIP’s sign symmetry enforcement in Matlab. I’m not sure why using a different basis for the gram matrices would cause this much grief with how the problem is being sent to Mosek.

I had not responded yet because I have been spending a long time (amidst moving to a new job and international travel) trying to construct a sufficiently minimal case (because I don’t want to post all ~400 lines of code here, that’s egregious), but in the most accessible minimal example I have, the symmetry works for the critical R value. So thank you for reaching out, the last thing I was going to finish doing before reposting was to figure out how the symmetry enforcement + bridges changes the problem sent to Mosek. We have decided to email you the smallest and clearest possible minimal example to illustrate what is happening, as it is too large to post here (over 400 lines). The minimal example we are sending does not have this issue of near-zero matrix entries, but we can send you the script for that version as well if necessary. We will send more details in the email.

2 Likes