blackbox

blackbox can be used to apply (almost) arbitrary operators on sdpvar objects.

Syntax

y = blackbox(x,f)
y = blackbox(x,f,'derivative',df,'convex','increasing','definiteness','positive','bounds',@mybounder)

Examples

General nonlinear functions are supported through the callback nonlinear operator framework. Most functions (exponentials, logarithms, trigonometrics, etc.) are available from start, but the user can use blackbox as a last resort to define other (elementwise) functions, without going through the hassle of creating operators from scratch. It can also be a trick to circumvent symblic handling of very large complex operations which involves many functions but a small number of decision variables.

Basic use

The following example shows how we can use a nonlinear solver to find a local optimizer to a small trigonometric problem. All involved operators here are already supported natively, so the use of blackbox is only for illustration. The important difference is that YALMIP never sees cosine, sine or quadratics in the first approach. It will only be given iterates of \(x\) from the solver, and then call the black-box function and return the function value to the solver. In the standard second model, YALMIP will explicitly work with the computational tree of the expression.

sdpvar x
y = blackbox(x,@(x)(sin(cos(x^2))))
optimize([1 <= x <= 2],y)

optimize([1 <= x <= 2],sin(cos(x^2)))

If you have the function implemented in a file, you can use a function handle to the file.

y = blackbox(x,@myfile)

Vectorized elementwise functions are supported.

x = sdpvar(2,3);
y = blackbox(x,@(x)(sin(cos(x.^2))))
optimize([1 <= x <= 2],sum(sum(y)))

In the following example, we have a small number of decision variables, but the introduction of the exponential operator on a large linear function will introduce a large amount of symbolic expressions and auxilliary variables which will cause unnecessary performance loss.

A = randn(1000,2);b = randn(1000,1);
y = exp(A*[3;3] + b);

% Nonlinear least squares, slow version where YALMIP (due to its design) introduces
% intermediate variables  \\( z \\) to replace the term in the exponential and 
% append  and a constraint  \\( z==Ax+b \\)
x = sdpvar(2,1);
e = exp(A*x + b) - y;
optimize([],e'*e);

% Define operator which computes fit, YALMIP and solver only sees \\( x \\).
J = @(x)(sum((exp(A*x + b)-y).^2));
f = blackbox(x,J);
optimize([],f)

Supplying derivatives

As an example of a more complicated scenario, let us maximize \(\int_0^a (x^2-x^3)dx\) with respect to the integration limit \(a \geq 0\). To implement the integral, we use the MATLAB command integral which takes a function handle and the two limits, which in our case means one of the limit arguments is a decision variable.

f = @(a)(integral(@(x) x.^2-x.^3,0,a))
sdpvar a
Area = blackbox(a,f);
optimize(a>=0,-Area)

In our examples so far, the solver will use numerical differentiation, but if we have a gradient available, we can supply it. Calculus to the rescue!

f = @(a)(integral(@(x) x.^2-x.^3,0,a))
sdpvar a
df = @(a)(a^2-a^3);
Area = blackbox(a,f,'derivative',df);
optimize(a>=0,-Area)

Advanced use in global optimization

A black-box operator can be equipped with properties that will enable efficient use in the global solver BMIBNB. Consider the following nonconvex problem that will be cast as a MILP by YALMIP due to the nonconvex use of absolute values.

n = 15;
c = randn(n,1);
A = randn(3*n,n);
b = rand(3*n,1);
x = sdpvar(n,1);
Domain = [-1 <= x <= 1];
optimize([Domain, A*x <= b, c'*abs(x) == 0],c'*x)      

If we want to solve this using a nonlinear programming solver instead, we create an alternative computation (with a generous definition of derivative for this non-smooth operator)

f = @(x)(abs(x));
df = @(x)sign(x);
my_abs = blackbox(x,f,'derivative',df);
optimize([Domain, A*x <= b, c'*my_abs == 0],c'*x)      

This will now be solved using a general nonlinear solver, but as the problem is nonconvex, no guarantees on global optimality can be given. To combat this, we can try to use BMIBNB.

To search for the global optima, BMIBNB performs a plethora of strategies based on bound propagation and envelope approximations of nonlinear operators. This is applied also on black-box operators, but for this to work they have to have properties broadcasted that BMIBNB can use. At the moment, the only such properties are the function values and derivatives, which is too little (it will force BMIBNB to use sampling, which means all guarantees lost, and performance suffer).

The most basic property is the range of the function in an interval. We can define this in a function and create a handle to it

function [L,U] = my_bounds(xL,xU)
if xL<=0 && xU >= 0
   L = 0;
   U = max(abs(xL),abs(xU));
else
   L = min(abs(xL),abs(xU));
   U = max(abs(xL),abs(xU));
end

Alternatively use an anonymous function directly (anonymous functions do not support multiple outputs so instead a vector is returned). Some acrobatics leads us to a one-liner

B = @(L,U)([~(L<0 & U>0)*min(abs(L),abs(U)) max(abs(L),abs(U))]);

With this operator, BMIBNB will be able to process nodes in a safe manner, but performance will be poor as the bounds alone give very little information. The most important information is an outer approximation of the convex hull of the operator on an interval. Luckily, we do not have to define this manually, but can simply append information that the function is convex, and BMIBNB will be able to derive this envelope automatically using derivative information.

f = @(x)(abs(x));
df = @(x)sign(x);
B = @(L,U)([~(L<0 & U>0)*min(abs(L),abs(U)) max(abs(L),abs(U))]);

my_abs = blackbox(x,@(x)(abs(x)),'derivative',df,'bounds',B,'convex');

ops = sdpsettings('solver','bmibnb');
optimize([Domain, A*x <= b, 
          c'*my_abs == 0],c'*x,ops)

To help with selection of branching points, positions of stationary points can be described. This can help BMIBNB as a split in the origin will lead to two nodes where the operator is linear in both, thus having a trivially tight envelope approximation.

my_abs = blackbox(x,@(x)(abs(x)),'derivative',df,'bounds',B,'convex','stationary',0);

ops = sdpsettings('solver','bmibnb');
optimize([Domain, A*x <= b, 
          sum(my_abs) >= n/2],c'*x,ops)

With a known stationary point and a convex function, the range of the function on an interval can easily be deduced. BMIBNB understands this and the bound generator is no longer required. Hence, we can just leave it at

my_abs = blackbox(x,@(x)(abs(x)),'derivative',df,'convex','stationary',0);

ops = sdpsettings('solver','bmibnb');
optimize([Domain, A*x <= b, 
          sum(my_abs) >= n/2],c'*x,ops)

A final note: The work above is just for illustration. The MILP model with a strong MILP solver will be the best approach for this particular case.

Comments

The command also supports \(R^n \rightarrow R \) functions, but this is not recommended when using BMIBNB as no convex hull generation is performed.