Nonlinear operators


A recent addition to YALMIP is simplified modeling of nonlinear, often non-differentiable, operators that typically occur in convex programming. Six operators are currently supported: min, max, abs, norm, sumk, sumabsk and detinv. These operators can be used rather intuitively, and YALMIP will automatically try to find out if they are used in a way that enables a convex representation. Although this can simplify the modeling phase significantly in some cases, it is recommended not to use these operators unless you know how to model them by your self, why and when it can be done etc. The text-book [S. Boyd and L. Vandenberghe] should be a suitable introduction for the beginner.

Consider once again the linear regression problem.

a = [1 2 3 4 5 6]';
t = (0:0.2:2*pi)';
x = [sin(t) sin(2*t) sin(3*t) sin(4*t) sin(5*t) sin(6*t)];
y = x*a+(-4+8*rand(length(x),1));
a_hat = sdpvar(6,1);
residuals = y-x*a_hat;

Using abs and max, we can easily solve the L1 and the L problem (Note that the abs operator currently has performance issues and should be avoided. Moreover, explicitly creating absolute values when minimizing the L error is not needed).

solvesdp([],sum(abs(residuals)));
a_L1 = double(a_hat)
solvesdp([],max(abs(residuals)));
a_Linf = double(a_hat)

YALMIP automatically concludes that the objective functions can be modeled using some additional linear inequalities, adds these, and solves the problems. We can simplify the code even more by using the norm operator (this is faster for large-scale problems due to implementation issues in YALMIP). Here we also compute the least-squares solution (note that this norm will generate a second-order cone constraint).

solvesdp([],norm(residuals,1));
a_L1 = double(a_hat)
solvesdp([],norm(residuals,2));
a_L2 = double(a_hat)
solvesdp([],norm(residuals,inf));
a_Linf = double(a_hat)

The following piece of code shows how we easily can solve a regularized problem.

solvesdp([],1e-3*norm(a_hat,2)+norm(residuals,inf));
a_regLinf = double(a_hat)

The norm operator is used exactly as the built-in norm function in MATLAB, both for vectors and matrices. Hence it can be used also to minimize the largest singular value (2-norm in matrix case), or the Frobenious norm of a matrix.

The double command of-course applies also to the nonlinear operators (double(OPERATOR(X)) returns OPERATOR(double(X)).

double(1e-3*norm(a_hat,2)+norm(residuals,inf))
ans =
    3.1175

A construction useful for maximizing determinants of positive definite matrices is the function (det P)1/m where m is the smallest power of 2 larger than or equal to the dimension of P (hence, for matrices with dimension equal to a power of two, the function gives the geometric mean of the eigenvalues of P). This concave function, called geomean2 in YALMIP, is supported as an extended operator.

D = randn(5,5);
P = sdpvar(5,5);
solvesdp(set(P < D*D'),-geomean2(P));

Advanced use

Rather advanced constructions are possible, and YALMIP will try derive an equivalent convex model.

sdpvar x y z
F = set(max(1,x)+max(y^2,z) < 3)+set(max(1,-min(x,y)) < 5)+set(norm([x;y],2) < z);
sol = solvesdp(F,max(x,z)-min(y,z)-z);

If the convexity propagation fails, an error will be issued (error code 14). Note that this does not imply that the model is nonconvex, but only means that the simple sufficient conditions used for checking convexity were violated. Failure is however typically an indication of a bad model, and most often due to an actual nonconvex part in the model. The problem above is convex, but not this problem below, due to the use of the minimizer in the constraint.

sdpvar x y z
F = set(max(1,x)+max(y^2,z) < 3)+set(max(1,min(x,y)) < 5)+set(norm([x;y],2) < z);
sol = solvesdp(F,max(x,z)-min(y,z)-z);
sol.info

ans =
 Convexity propagation failed (YALMIP)

In the same sense, this problem fails due to a nonconvex objective function.

sdpvar x y z
F = set(max(1,x)+max(y^2,z) < 3);
sol = solvesdp(F,-norm([x;y]));
sol.info

ans =
 Convexity propagation failed (YALMIP)

Adding new operators

If you want to add your own operator, all you need to do is to create 1 file. This file should be able to return the numerical value of the operator for a numerical input, and return the epigraph (or hypograph) of the operator when the first input is 'graph'. As an example, the following file implements the nonlinear operator tracenorm. This convex operator returns sum(svd(X)) for matrices X. This value can also be described as the minimizing argument of the optimization problem mint,A,B t subject to set([A X;X' B] > 0) + set(trace(A)+trace(B) < 2*t).

function varargout = tracenorm(varargin)

switch class(varargin{1})    

    case 'double' % What is the numerical value of this argument (needed for displays etc)
        varargout{1} = sum(svd(varargin{1}));

    case 'char'   % YALMIP send 'graph' when it wants the epigraph or hypograph
        if isequal(varargin{1},'graph')
            t = varargin{2}; % Second arg is the extended operator variable
            X = varargin{3}; % Third arg and above are the args user used when defining t.
            A = sdpvar(size(X,1));
            B = sdpvar(size(X,2));
            F = set([A X;X' B] > 0) + set(trace(A)+trace(B) < 2*t);
            varargout{1} = F; % Epigraph model
            varargout{2} = 1; % Convex operator (-1 for concave)            
        else           
        end    

    case 'sdpvar' % Always the same. 
        varargout{1} = yalmip('addextendedvariable',mfilename,varargin{:});    

    otherwise
end

The function sumk.m in YALMIP is implemented using this framework and serve as good examples. The overloaded operator norm.m is also defined using this method, but is a bit more involved, since it supports different norms.

Limitations

The described operators cannot be used in polynomial expressions in the current implementation. The following problem is trivially convex but fails.

sdpvar x y
sol = solvesdp([],norm([x;y])^2);
sol.info

ans =
 Convexity propagation failed (YALMIP)

The following constructions is convex but fails.

sdpvar x y
solvesdp([],norm([x;max(1,x)]));
sol.info

ans =
 Convexity propagation failed (YALMIP)

Another limitation is that the operators not are allowed in cone and semidefinite constraints.

sdpvar x y
solvesdp(set(cone(max(x,y,1),2)),x+y);
sol.info

ans =
 Convexity propagation failed (YALMIP)

In practice, these limitations should not pose a major problem. A better model is possible (and probably recommended) in most cases if these situations occur. Note also the convexity checking algorithm is experimental and very simple,, and will most likely be improved upon in a future release.