One approach for exploiting sparsity in polynomial optimization is based on correlative sparsity and proposed by Waki, Kim, Kojima and Muramatsu. We motivate this approach by considering the following problem,
We first setup and solve the standard second-order moment relaxation
In [1]:
using Polyopt
In [2]:
x1,x2,x3 = variables(["x1","x2","x3"]);
In [3]:
f = x1-x2+x3 - (x1^3+x2^3+x3^3);
In [4]:
g = [-1-x1*x2, 1-x1*x2,
-1-x2*x3, 1-x2*x3,
1-x1^2,
1-x2^2,
1-x3^2 ];
In [5]:
prob = momentprob(2, f, g);
In [6]:
X, Z, t, y, solsta = solve_mosek(prob);
We inspect the solution
In [7]:
[ prob.basis y ]
Out[7]:
which is feasible
In [8]:
xo = Polyopt.vectorize([x1,x2,x3],4)*y
Out[8]:
In [9]:
[ Polyopt.evalpoly(gi, xo) for gi=g ]
Out[9]:
so we have found an optimal solution.
The correlative sparsity pattern $C$ is a symmetric matrix with the number variables as dimension. An element $C_{ij}$ is nonzero if either
For the simple example, let us recall
In [10]:
f
Out[10]:
In [11]:
g
Out[11]:
Let us compute the correlative sparsity pattern
In [12]:
C = Polyopt.correlative_sparsity(f, g);
In [13]:
full(C)
Out[13]:
This is a chordal matrix with 2 cliques,
In [14]:
K = Array{Int,1}[ [1,2], [2,3] ]
Out[14]:
In the standard moment relaxation we consider a monomial basis vector
In [15]:
u = monomials(2, [x1,x2,x3])
Out[15]:
and the second-order moment matrix $M_2(y)\succeq 0$ corresponds to
In [16]:
u*u'
Out[16]:
In the chordal relaxation we consider basis vectors formed over the cliques
In [17]:
uk = [ monomials(2, [x1,x2,x3][ki]) for ki=K ]
Out[17]:
and the moment constraint $M_2(y)\succeq 0$ is replaced by two smaller constraints $M_2(y, K_i)\succeq 0$ corresponding to
In [18]:
uk[1]*uk[1]'
Out[18]:
In [19]:
uk[2]*uk[2]'
Out[19]:
Similarly we consider the localization matrices for $g_i(x)\geq 0$. For the standard second-order moment relexation we have
In [20]:
v = monomials(1, [x1,x2,x3])
Out[20]:
In [21]:
g[1]
Out[21]:
and the localization matrix for $g_1(x) \geq 0$ is $M_1(g_1 y)\succeq 0$ corresponding to
In [22]:
g[1]*v*v'
Out[22]:
In the chordal relaxation we have
In [23]:
vk = [ monomials(1, [x1,x2,x3][ki]) for ki=K ]
Out[23]:
and we replace $M_1(g_1 y)\succeq 0$ by $M_1(g_1 y, K_1)\succeq 0$ corresponding to
In [24]:
g[1]*vk[1]*vk[1]'
Out[24]:
Similarly we replace $M_1(g_2y)$ by
In [25]:
g[2]*vk[1]*vk[1]'
Out[25]:
and so on. The function $g_6(x)=1-x_2^2$ belongs to both clique 1 and 2, so we replace $M_1(g_6 y)\succeq 0$ by two smaller inequalities
In [26]:
g[6]*vk[1]*vk[1]'
Out[26]:
and
In [27]:
g[6]*vk[2]*vk[2]'
Out[27]:
We form the chordal moment relaxation as
In [28]:
probc = Polyopt.momentprob_chordal(2, K, f, g);
and we solve it as
In [29]:
Xc, Zc, tc, yc, solstac = solve_mosek(probc);
In [30]:
full([ probc.basis yc ])
Out[30]:
In this case we also recover the optimal value at the second order relaxation. In general this is not always the case, except when all polynomials are quadratic, but under mild assumptions the hierarchy of chordal moment relaxations converges to the optimal solution (just as the regular moment relaxations).
The sums-of-squares certificate is also evaluted over the cliques, in this case we have
In [31]:
r = dot(uk[1], Xc[1]*uk[1]) + dot(uk[2], Xc[2]*uk[2]) +
g[1]*dot(vk[1], Xc[3]*vk[1]) +
g[2]*dot(vk[1], Xc[4]*vk[1]) +
g[3]*dot(vk[2], Xc[5]*vk[2]) +
g[4]*dot(vk[2], Xc[6]*vk[2]) +
g[5]*dot(vk[1], Xc[7]*vk[1]) +
g[6]*(dot(vk[1], Xc[8]*vk[1]) + dot(vk[2], Xc[9]*vk[2])) +
g[7]*dot(vk[2], Xc[10]*vk[2]);
In [32]:
truncate( f - tc - r, 1e-6 )
Out[32]:
When the correlative sparsity pattern is non-chordal it is difficult to find the cliques (it is NP-hard in general). In that case we embed the correlative sparsity pattern in a larger chordal graph for which we can easily detect the cliques. The chordal embedding procedure essentially corresponds to a symbolic Cholesky factorization of the correlative sparsity pattern, and the amount of generated fill-in depends on the quality of a symmetric matrix reordering (performed by CHOLMOD).
To demonstrate this procedure we consider a problem from http://gamsworld.org/global/globallib/ex5_4_2.htm, where we have rescaled the problem such that $x_i\leq 1$ to make the relaxations easier to solve.
In [33]:
x1,x2,x3,x4,x5,x6,x7,x8 = variables(["x1","x2","x3","x4","x5","x6","x7","x8"]);
In [34]:
f = x1 + x2 + x3;
In [35]:
g = [ 400/1000 - (x4 + x6),
300/1000 - (-x4 + x5 + x7),
100/1000 - (-x5 + x8),
.083333 - (0.01*x1 - 10*x1*x6 + 0.83333*x4),
-(x2*x4 - x2*x7 - 0.125*x4 + 0.125*x5),
-0.125 - (x3*x5 - x3*x8 - 0.25*x5),
x1-100/10000, 1-x1,
x2-1000/10000, 1-x2,
x3-1000/10000, 1-x3,
x4-10/1000, 1-x4,
x5-10/1000, 1-x5,
x6-10/1000, 1-x6,
x7-10/1000, 1-x7,
x8-10/1000, 1-x8 ];
For reference, let us first solve the third-order standard moment relaxation.
In [36]:
prob = momentprob(3, f, g);
In [37]:
X, Z, t, y, solsta = solve_mosek(prob);
and the extracted solution
In [38]:
[ prob.basis y ]
Out[38]:
is feasible
In [39]:
xo = Polyopt.vectorize([x1,x2,x3,x4,x5,x6,x7,x8],6)*y
Out[39]:
In [40]:
[ Polyopt.evalpoly(gi, xo) for gi=g ]
Out[40]:
Now, let us form the chordal moment relaxation using a chordal embedding,
In [41]:
probc = momentprob_chordalembedding(3, f, g);
In [42]:
Xc, Zc, tc, yc, solstac = solve_mosek(probc);
In [43]:
full([ prob.basis y yc ])
Out[43]:
In [ ]: