Bilinear SOS optimization for rank-1 tensor-product model
Simon Wiesheier
Hello everyone,
I consider a rank-1 tensor-product (spline) model
W(x,y) = sum_{i,j} a(i)*b(j)*N_i(x)*N_j(y),
i.e., the spline coefficient matrix C is constrained to satisfy
C = a*b' (rank one).
My goal is to fit the parameters (a,b) to data while enforcing that the Hessian of W(x,y) is positive semidefinite over each spline patch using SOS matrix constraints.
With free spline coefficients C(i,j), all polynomial coefficients entering the SOS constraints are affine in the optimization variables, so the problem is a standard SDP and can be solved with 'sdpt3' or friends.
With the rank-1 parameterization C = a*b', the polynomial coefficients become bilinear in (a,b), and the least-squares data-fitting objective is likewise bilinear in (a,b).
Is there any built-in machinery in YALMIP to handle this type of bilinear SOS problem, beyond implementing a block-alternating optimization manually (optimizing a for bixed b, then optimizing b for fixed a, until convergence)? For example, are there recommended solver interfaces, lifting techniques, or moment/SOS formulations that are practical for this type of model?
Any suggestions would be greatly appreciated.
Thank you and best,
Simon
Johan Löfberg
Supply an example and we can see what can be done'
https://yalmip.github.io/tags/#polynomials
Simon Wiesheier
Hi Johan, Thank you for your reply.
I reduced my application to the attached minimal example. The real problem is the calibration of a function represented by tensor-product B-splines. On each spline patch, I enforce positive semidefiniteness of the Hessian using a localized matrix-SOS certificate,
F = [F, sos(z'*S0*z), sos(z'*S1*z), sos(z'*S2*z)]; R = H - S0 - (1-x^2)*S1 - (1-y^2)*S2; coeff = coefficients(z'*R*z,[x;y;z1;z2]); F = [F, coeff == 0];The attached example replaces the spline basis by a simple polynomial bases on a single patch ((x,y) \in [-1,1]^2), but has the same optimization structure.
The synthetic target is generated from a known rank-one
coefficient matrix with a PSD Hessian on the whole patch.
(i) The first ("full") model optimizes all coefficients of W(x,y)=N_x(x)^T C N_y(y), directly. Therefore, the objective is linear in the decision variables. SDPT3 solves this problem easily and recovers the exact coefficients.
(ii) The second model parameterizes C = ab^T,
which is what I ultimately want (in the real application this corresponds to the product of two univariate spline expansions). This makes both the least-squares objective and the SOS coefficient matching nonlinear in the decision variables. I supplied explicit bounds on all entries of (a) and (b), since this model has an inherent scaling ambiguity.
My
question is: what would you recommend as the appropriate YALMIP
modeling/solver strategy for this rank-one nonlinear SDP? Should I use
BMIBNB or another nonlinear SDP solver, or is there a better built-in
approach in YALMIP for this type of problem?
Ultimately,
the real application involves maybe 10-50 spline patches, each with its
own localized SOS constraints, so the final strategy should also scale
beyond this minimal example. But that should work out I guess.
Many thanks for your time!
Simon Wiesheier
As a follow-up to my previous post, I tried a slightly larger version of the same example (the code is only ~70 lines and inlined below).
Instead of a single local polynomial patch, I simply duplicated the problem. Each patch uses the same local coordinates ((x,y)\in[-1,1]^2), the same localized matrix-SOS certificate, and the target is constructed by repeating the original rank-one factors (repmat), so every local block is identical to the original feasible one-patch example.
The full formulation still solves immediately, but BMIBNB now reports (Upper = Inf), although the duplicated problem is feasible by construction. The lower solver is MOSEK.
Is this problem already too large for BMIBNB? If so, what alternative strategy/solver would you recommend for this type of rank-one nonlinear SDP?
My actual application consists of roughly 10–50 such localized SOS constraints, so I would also appreciate your opinion on whether this overall approach to enforcing positive semidefiniteness via localized matrix-SOS certificates makes sense, or whether you would formulate the problem differently. The motivation to enforce PSD locally was avoiding overly restrictive global coefficient constraints.
Many thanks.
CODE
%% Minimal 25-patch rank-1 Hessian-SOS example
clear; clc;
yalmip('clear');
nPatch = 2; % (nPatch = 1 works perfectly)
nCoef = 3*nPatch;
sdpvar x y z1 z2
z = [z1;z2];
Nx = [1;x;x^2];
Ny = [1;y;y^2];
a = sdpvar(nCoef,1);
b = sdpvar(nCoef,1);
C = a*b.';
% Same PSD rank-1 target as in the one-patch example
aTrue = repmat([1;0.10;0.05],nPatch,1);
bTrue = repmat([1;-0.10;0.05],nPatch,1);
CTrue = aTrue*bTrue.';
objective = sum((C(:)-CTrue(:)).^2);
g1 = 1-x^2;
g2 = 1-y^2;
F = [a(1) == 1, -5 <= a <= 5, -5 <= b <= 5];
for i = 1:nPatch
for j = 1:nPatch
I = 3*i-2:3*i;
J = 3*j-2:3*j;
W = Nx.'*C(I,J)*Ny;
H = [jacobian(jacobian(W,x),x), ...
jacobian(jacobian(W,x),y);
jacobian(jacobian(W,x),y), ...
jacobian(jacobian(W,y),y)];
S0 = sympoly2([x;y],2);
S1 = sympoly2([x;y],0);
S2 = sympoly2([x;y],0);
F = [F, ...
sos(z.'*S0*z), ...
sos(z.'*S1*z), ...
sos(z.'*S2*z)];
R = H-S0-g1*S1-g2*S2;
F = [F, coefficients(z.'*R*z,[x;y;z1;z2]) == 0];
end
end
ops = sdpsettings( ...
'solver','bmibnb', ...
'bmibnb.lowersolver','mosek', ...
'verbose',1, ...
'sos.model',1);
diagnostics = solvesos(F,objective,ops,[a;b]);
disp(yalmiperror(diagnostics.problem))
disp(value(objective))
disp(value(C))
function S = sympoly2(v,d)
[p11,~] = polynomial(v,d);
[p12,~] = polynomial(v,d);
[p22,~] = polynomial(v,d);
S = [p11 p12; p12 p22];
end
Johan Löfberg
The local nonlinear solver in YALMIP is called "Conan the Slayer", short for "Conic and nonlinear semidefinite layer" is an experimental layer mapping a semidefinite model on-the-fly to a model fmincon can work with. It has an algorithmic setting and switching this appears to save the day which is a nice surprise
ops = sdpsettings( ...
'solver','bmibnb', ...
'bmibnb.lowersolver','mosek', ...
'verbose',1, ...
'sos.model',1,'slayer.algorithm','nonconvex');
* Starting YALMIP global branch & bound.
* Upper solver : fmincon
* Lower solver : MOSEK
* LP solver : MOSEK
* -Extracting bounds from model
* -Perfoming root-node bound propagation
* -Calling upper solver (no solution found)
* -Branch-variables : 12
* -More root-node bound-propagation
* -Performing LP-based bound-propagation
* -And some more root-node bound-propagation
* Starting the b&b process
Node Upper Gap(%) Lower Open Time
1 : Inf NaN 3.55066E-09 2 15s
2 : 1.35012E-07 0.00 2.06594E-09 3 17s Solution found by upper solver
* Finished. Cost: 1.3501e-07 (lower bound: 2.0659e-09, relative gap 1.3295e-05%)
* Termination with relative gap satisfied
Also, I would probably avoid a quartic model and change the model to
e = sdpvar(length(C(:)),1);
objective = e'*e;
g1 = 1-x^2;
g2 = 1-y^2;
Simon Wiesheier
Many thanks for your suggestions, they helped.
With the updated formulation (essentially adding 'slayer.algorithm','nonconvex' to ops) I can now solve the duplicated example for nPatch = 3. However, already for nPatch = 5, BMIBNB spends several minutes in the root-node processing and never even starts branching:
* Starting YALMIP global branch & bound.
* Upper solver : fmincon
* Lower solver : MOSEK
* LP solver : MOSEK
* -Extracting bounds from model
* -Perfoming root-node bound propagation
* -Calling upper solver (no solution found)
* -Branch-variables : 30
* -More root-node bound-propagation
* -Performing LP-based bound-propagation
* -And some more root-node bound-propagation
Have we reached the practical limits of BMIBNB here? If so, what approach or solver would you recommend for larger problems of this type?
Many thanks!