Beyond Convex for Global Optimization

Paritosh Mokhasi,  
Kernel Developer, Algorithms R&D

(RowDefault[, (Spacer1), , (Spacer1), ])  #WolframTechConf

ABSTRACT

In the field of optimization, convex optimization holds special status because of its property that the minimum is always a global minimum and there are highly efficient solvers available to solve convex problems. However, not all optimization problems can be formulated as purely convex problems. This talk will focus on generalizations of the convex optimization algorithms that allows us to solve a much broader category of problems. We will discuss the construction of generalized convex functions and provide details on how such problems are solved within the Mathematica optimization framework.

Convex Optimization

A convex optimization problem is one of the form
			
				minimize    (Subscript[f, 0])(x)
				subject to    (Subscript[f, i])(x)≤0, i=1,2,...,m,
				
	where the functions (Subscript[f, 0])(x),...,(Subscript[f, m])(x): (ℝ^n)->ℝ are convex, i.e., satisfy
	
				(Subscript[f, i])(α x+β y)≤ α (Subscript[f, i])(x)+β (Subscript[f, i])(y)
				
	for all x,y∈(ℝ^n) and all α,β∈ℝ with α+β=1,α≥0,β≥0.

Convex Optimization

Minimize 	         f(x),
		     subject to	        (Subscript[a, 0]).x+(Subscript[b, 0])=0,
(Subscript[a, j]).x+(Subscript[b, j])(Subscript[\[VectorGreaterEqual], (Subscript[κ, j])])0, j=1... k

Beyond Convex Optimization: Quasiconvex

A quasiconvex optimization problem is one of the form
			
				minimize    (Subscript[f, 0])(x)
				subject to    (Subscript[f, i])(x)≤0, i=1,2,...,m,
				
	where (Subscript[f, 1])(x),...,(Subscript[f, m])(x): (ℝ^n)->ℝ are convex, and (Subscript[f, 0])(x) is quasiconvex: 
	
				(Subscript[f, 0])(α x+β y)≤ max.{(Subscript[f, 0])(x),(Subscript[f, 0])(y)}
				
	for all x,y∈(ℝ^n) and all α,β∈ℝ with α+β=1,α≥0,β≥0.

Quasiconvex functions: Properties.

The function is quasiconvex if its domain and all its sublevel sets are convex:    (Subscript[S, α])={x∈dom. f|f(x)≤α}, α∈ℝ

The function is nondecreasing.

The function is nonincreasing.

There is a point c∈dom.f such that for t≤c, f is nonincreasing and for t≥c, f is nondecreasing.

Solving Quasiconvex Problems

Consider the following example:

refSol=NMinimize[{-x y, x≥0,y≥0,Norm[{x-1,y-2}]≤1},{x,y}]

NMinimize[{-x y, x≥0,y≥0,Norm[{x-1,y-2}]≤1},{x,y},Method->#]&/@{"QuasiconvexBisection","QuasiconvexSubgradient"}

ContourPlot[-x y,{x,0,2},{y,1,3},Contours->15,ImageSize->400,ContourLabels->True,RegionFunction->Function[{x,y},Norm[{x-1,y-2}]≤1]]

Solving Quasiconvex Problem: Bisection Method

This quasiconvex problem can be solved by seeking a feasible point for a specified sublevel set.

pfun = ParametricConvexOptimization[0,{-x y≤α, x≥0,y≥0,Norm[{x-1,y-2}]≤1},
{x,y},α,{"PrimalMinimizerRules","SolutionStatus"}];

Quiet[{pfun[-2],pfun[-5]}]

{lb,ub}={-5.,-2.};
Reap[Do[
midSol=Quiet[pfun[Sow[midpt = (lb+ub)/2]]];
If[midSol[[2,1]]==="Infeasible", lb=midpt, ub=midpt];
,{iter,10}]][[2,1]]

Solving Quasiconvex Problems: Subgradient Method.

Subgradient methods are very simple: Get to the optimal solution using a simple iteration of the form: (x^((k+1)))=Π((x^((k)))-(Subscript[α, k])g((x^((k))))), where g(x) is called the subgradient and Π is the projection of point on the convex set.

refSol=NMinimize[{-x y, x≥0,y≥0,Norm[{x-1,y-2}]≤1},{x,y}]

pfun=ParametricConvexOptimization[Norm[p-q],{Inactive[Norm][Inactive[Plus][p,-{1,2}]]≤1,
p\[VectorGreaterEqual]0},p∈Vectors[2,Reals],q∈Vectors[2,Reals], "PrimalMinimizerVector"];

f[{x_,y_}]:=-x*y; df[{x_,y_}]:= {-y,-x};
xx={1,1.1};fbest=Infinity;
res=Reap[Do[
Sow[xx = pfun[xx]];
{fx,gx}={f[xx],df[xx]};
fbest = Min[fbest,fx];
xx -=(((fx-fbest)+ (1/(1+iter)))/gx^2)*gx;
,{iter,30}]][[2,1]];

ContourPlot[-x y,{x,0,2},{y,1,3},Contours->15,ImageSize->400,ContourLabels->True,RegionFunction->Function[{x,y},Norm[{x-1,y-2}]≤1],Epilog->{Arrow[Partition[res,2,1]],{PointSize[0.03],Red,Point[{x,y}/.refSol[[2]]]}}]

Beyond Convex Optimization: Difference of convex.

A difference of convex optimization problem is one of the form
			
				minimize    (Subscript[f, 0])(x)-(Subscript[g, 0])(x)
				subject to    (Subscript[f, i])(x)-(Subscript[g, i])(x)≤0, i=1,2,...,m,
				
	where the functions (Subscript[f, 0])(x),...,(Subscript[f, m])(x) and (Subscript[g, 0])(x),...,(Subscript[g, m])(x): (ℝ^n)->ℝ are convex
	
				(Subscript[f, i])(α x+β y)≤ α (Subscript[f, i])(x)+β (Subscript[f, i])(y),
(Subscript[g, i])(α x+β y)≤ α (Subscript[g, i])(x)+β (Subscript[g, i])(y)
				
	for all x,y∈(ℝ^n) and all α,β∈ℝ with α+β=1,α≥0,β≥0.

region=1<=Norm[{x,y}]<=3&&1/2<=Norm[{x-1/2,y-1/2}]<=2&&1/2<=Norm[{x+1/2,y-1/2}]<=2;
RegionPlot[region,{x,-2,2},{y,-2,3}, ImageSize->450]

Difference of convex: Some examples.

The objective has a difference of convex form

res=NMinimize[{x^2-y^2,Norm[{x,y}]<=2},{x,y}]
Show[Plot3D[x^2-y^2,{x,-2,2},{y,-2,2},],Graphics3D[{PointSize[0.04],Point[{x,y,x^2-y^2}/.res[[2]]]}], ImageSize->400]

One of the constraints has a difference of convex form

res=NMinimize[{x^2+y,1<=Norm[{x,y}], Norm[{x,y}]<=2},{x,y}]
Show[Plot3D[x^2+y,{x,-2,2},{y,-2,2},],Graphics3D[{PointSize[0.04],Point[{x,y,x^2+y}/.res[[2]]]}]]

Difference of convex: Solution Strategy

If a function h(x)=f(x)-g(x) and f,g are convex, then h(x) can be approximated as h(x)=f(x)+g((Subscript[x, 0]))+∇g.(x-(Subscript[x, 0])) making the function convex. Consider:

res=NMinimize[{x^2+y,1<=Norm[{x,y}], Norm[{x,y}]<=2},{x,y}]

ContourPlot[x^2+y,{x,-2,2},{y,-2,2},]

pfun = ParametricConvexOptimization[x^2+y + w*alp,{Norm[{x,y}]<=2,
-Norm[{x0,y0}]-({x0,y0}/Sqrt[x0^2+y0^2]).({x,y}-{x0,y0})<=-1+alp, alp>=0},{x,y,alp},{x0,y0,w}];

pfun[0.5,1.1,1]

Difference of Convex: Path Planning

Find a path through circular obstacles such that the distance between the start and end points is minimized:

p = {{2,2},{7,6},{3,5},{4,8},{8,9},{6,2.5}};
r={1,1.3,1.5,1,0.7,1.5};
{start,end}={{0,0},{10,10}};
domain=Graphics[{MapThread[Circle[#1,#2]&,{p,r}],
{PointSize[0.02],Point[{start,end}]}},Frame->True, ImageSize->450]

Difference of Convex: Path Planning

The path is discretized into n different points with distance of l/n between points, where l is the trajectory length being minimized:

n=30;
distanceConstraints=Table[Inactive[Norm][x[i]-x[i-1]]<=l/n,{i,n}];

The points cannot be inside the circular objects:

objectConstraints=Table[(InactiveHead[Norm])[-p[[j]]+x[i]]>=r[[j]],{i,1,n},{j,1,6}];

The start and end points are known:

positionConstraints={x[0]==start,x[n]==end};

Collect the variables:

vars=Append[Table[x[i]∈Vectors[2,Reals],{i,0,n}],l∈Reals];

Difference of Convex: Path Planning

Minimize the length l subject to the constraints:

Short[res=NMinimize[{l,distanceConstraints,objectConstraints,positionConstraints},vars],5]

Visualize the result:

pts=vars[[1;;-2,1]]/.res[[2]];
Show[domain, Graphics[{Red,PointSize[0.02],Point[pts]}], ImageSize->450]

Difference of Convex: Circle Packing

Find the smallest square that can contain n circles of given radius (Subscript[r, i]) for i=1,...,n and do not overlap. Specify the number of circles and the radius of each circle:

n=20;
r = BlockRandom[RandomReal[{0.5,1.5},n],RandomSeeding->123];

If (Subscript[c, i]) is the center of circle i, then the objective is to minimize (max_(i = 1,..,n))(||(Subscript[c, i])(Subscript[||, ∞])+ (Subscript[r, i])). The objective can be transformed so as to minimize s and -s<=||(Subscript[c, i])(Subscript[||, ∞])+(Subscript[r, i])<=s, i=1,...,n:

objectiveConstraint = Table[-s<=Norm[c[i],Infinity]+r[[i]]<=s,{i,n}];

The circles must not overlap:

nonOverlapConstraint = Table[Norm[c[i]-c[j]]>=r[[i]]+r[[j]],{i,1,n},{j,i+1,n}];

Collect the variables:

vars=Append[Table[Element[c[i],Vectors[2,Reals]],{i,n}],Element[s,Reals]];

Difference of Convex: Circle Packing

Minimize the objective s:

Short[res=NMinimize[{s,objectiveConstraint,nonOverlapConstraint},vars],5]

The circles are contained in the square [-s,s] × [-s,s] :

Graphics[{Table[Circle[c[i],r[[i]]],{i,n}],{FaceForm[None],EdgeForm[Black],Rectangle[{-s,-s},{s,s}]}}]/.res[[2]]

Difference of Convex: Data Fitting

Minimize ||a.x-b|| subject to the constraints (Subscript[x, i])∈{-1,1}, i=1,2,…,n for a given matrix a and vector b:

SeedRandom[123];
n=20; m=2n;
{a,b,(Subscript[x, ref])}=Module[{a,b,xref},];

The constraint (Subscript[x, i])∈{-1,1} can be modeled as (Subscript[x, i]^2)=1 or  1<=((Subscript[x, i])^2)<=1. Find the solution:

res=NMinimize[{Inactive[Norm][Inactive[Plus][a.x,-b]],1\[VectorLessEqual]x^2\[VectorLessEqual]1},x∈Vectors[n,Reals]]

Thank you for your attention!