Affine Euclidean subspace manifold to enforce (strict) inequality constraints

[Ryan Harvey]

Hello,

I have an objective-function that is technically well defined on any x in R^5, where x = (d,r,v). However, it is only physically meaningful when d in S^2, r > 0, and v \in (v_min, v_max) with both v_min and v_max being positive. This is, as far as I'm aware, a feasibility region that works as a smooth manifold as it is an open subset of the product manifold S^2 times R^2+. 

Any suggestions on how one would go about creating an affine submanifold factory which enforces the strict inequality constraints between (v_min,v_max) (provided both are positive and not equal)?

At the moment, I have it working somewhat well in matlab using fmincon's sqp with strict enforcement of feasibility constraints, however I'm interested in exploring using Manopt to solve this kind of problem. It seems like it might lead to a more elegant/interesting solution approach than just throwing strict-feasibility SQP at it and hoping you get convergence. 

Best,

Ryan 

[Nicolas Boumal]

Hello Ryan,

As an aside: strict inequality constraint can be tricky, for the following reason. Imagine minimizing f(x) for x in the open interval (0, 1). Then, either the minimizer of f indeed lies strictly inside (0, 1) -- in which case the constraint is inactive, and so it does not play much of a role -- or f does not attain a minimizer in (0, 1), because its infimum on that interval corresponds to the value f(0) or f(1). In the latter scenario, any optimization algorithm should try to get to x = 0 or x = 1, but it wouldn't be allowed to.

If you are happy to try things out with non-strict inequalities, then here would be my first suggestion: let's work with a product manifold where d is your unit-norm vector in R^3 (that's just a sphere), and q in R^2 has two components, which we use as follows:

r = q(1)^2;

v = (vmin+vmax)/2 + (vmax - vmin)/2 * sin(q(2));

In Manopt, you can create the manifold as follows:

elements.d = spherefactory(3);

elements.q = euclideanfactory(2);

manifold = productmanifold(elements);

problem.M = manifold;

To make the code more transparent, I recommend introducing a change of variable function as follows:

function xx = change_of_variable(x)

    xx.d = x.d;

    q = x.q;

    xx.r = q(1)^2;

    xx.v = (vmin+vmax)/2 + (vmax - vmin)/2 * sin(q(2));

end

problem.cost = @mycost;

function f = mycost(x)

    xx =  change_of_variable(x);

    % ... compute f as a function of xx, which is a structure that contains d, r and v as you normally would

end

From here, you could try to call automatic differentation with:

problem = manoptAD(problem);

This may fail (you need the Deep Learning toolbox, and rather recent version of Matlab, and it can be a bit finnicky).

If so, then you would ideally implement the gradient of your cost function. This must be expressed with respect to d and q (because that is the manifold), so you would compute the gradient of the composition of f with the change of variable (chain rule).

I hope this helps -- don't hesitate to follow up here.

Best,

Nicolas

[Nicolas Boumal]

Of course, you could also try a different change of variable that captures the strict inequalities, simply by mapping R to the strictly positive numbers -- e.g., r = exp(q(1)) -- and mapping R to the open interval (vmin, vmax) -- e.g., with a sigmoid of your choosing).

My prior is that this might not work as well, but in such matters, the only way to know is to try.

If you are going to experiment with various changes of variable, you may find it convenient to abstract away the role of the change of variable and the role of the cost function f, so that f and its derivatives can be implemented once, and the change of variable and its derivatives (necessary for the chain rule) are implemented separately.

In the most recent version of Manopt, the idiomatic way of doing this is via lifts:

https://www.manopt.org/lifts.html

[Ryan Harvey]

Thank you very much for the quick and thorough reply,

I will give your strategy a try and report back.

Best,

Ryan 

[Ryan Harvey]

Hello Nicolas,

The Deep Learning AD had some problems with the way my cost function operates (though, I think with some patient reworking of the inner functions those could be resolved), but it seems to work relatively well with finite differences for the Hessian. The main thing I need to figure out is smart stopping criteria. 

The cost surface is itself fundamentally dependent on a random variable, and the minima can be very ill-conditioned and flat (depending on the realization of the random variable and the fundamental variance covariance matrix of the measurement errors). Do you have any recommendations for implementing stopping criteria that play well with the fact that the cost function is itself fundamentally a function of a random variable? My first-pass at the problem is just implementing the 'mystopfun' custom stopping conditions but with the tolerance a good bit tighter (some back of the napkin math led me to 1e-10 to work well for my case).

Do you have any general rules of thumb for smart design of stopping criteria when the cost surface itself is dependent on a random variable whose covariance matrix/conditioning is not necessarily known ahead of time (except for maybe an approximation)? 

Even with my first-pass at implementing your change-of-variables idea, this is already significantly faster than the fmincon approach, and the guarantee that it has a feasible solution is a plus. So that's already quite good. I guess the next question will be, how are the statistics after I run some trials? Finger's crossed.

Thank you again for all of your help,

Ryan 

[Nicolas Boumal]

Hello Ryan,  good to hear this seems promising.  About handling randomness, and how to integrate that into the stopping criterion, unfortunately I don't have something concrete to offer.