convex optimization

convex optimization

This page shows a few initial experiments with the basic algorithms described in the excellent Convex Optimization book by Boyd and Vandenberghe.

In the following we use the following standard parameters:

optimize = interiorBarrierII α β μ0 μ
  where
    α = 0.01
    β = 0.8
    μ0 = 10
    μ = 10

phaseI = interiorBarrierI α β μ0 μ
  where
    α = 0.01
    β = 0.8
    μ0 = 10
    μ = 10

optimization without constraints

Minimize the quadratic function

$$(x-2)^2 + (y-2)^2 + (z-2)^2$$

g0 = mkQuadAt (vector [2,2,2]) 0 (vector [0,0,0]) (diagl [2,2,2])

Any starting point is feasible and we get the exact solution in just one step:

sh 6 $ optimize 1E-4 g0 [] (matrix 3 []) (vector[]) (vector[0.5,0,0])
1x3
0.500000  0.000000  0.000000
1x3
2  2  2

equality constraints

Minimize

$$(x-2)^2 + (y-2)^2 + (z-2)^2$$

subject to

$$x+y+z=1$$

sh 6 $ optimize 1E-4 g0 [] (matrix 3 [1,1,1]) (vector[1]) (vector[0,0,0])
1x3
0  0  0
1x3
0.333333  0.333333  0.333333

The objective function is exactly quadratic, so the solution is obtained again in a single step (even though in this case the starting point is not feasible).

We can add another equality constraint, for example:

$$y=0.5$$

optimize 1E-4 g0 [] (matrix 3 [1,1,1,0,1,0]) (vector[1,0.5]) (vector [0,0,0])
1x3
0  0  0
1x3
0.250000  0.500000  0.250000

The feasible set is a line, and by symmetry we see that the solution is correct.

inequality constraints

Minimize

$$(x-2)^2 + (y-2)^2 + (z-2)^2$$

subject to

$$x + y + z \leq 1 $$

This constraint has the following representation:

g1 = mkLin (vector [1,1,1]) (-1)
sh 6 $ optimize 1E-4 g0 [g1] (matrix 3 []) (vector[]) (vector [0.5,0,0])
1x4
0.500000  0.000000  0.000000  -0.500000
7x4
0.600663  0.167772  0.167772  -0.063792
0.388893  0.302315  0.302315  -0.006478
0.329482  0.329482  0.329482  -0.011553
0.327120  0.327120  0.327120  -0.018641
0.324787  0.324787  0.324787  -0.025639
0.323587  0.323587  0.323587  -0.029238
0.323396  0.323396  0.323396  -0.029811
3x4
0.332458  0.332458  0.332458  -0.002625
0.332349  0.332349  0.332349  -0.002952
0.332334  0.332334  0.332334  -0.002997
5x4
0.333294  0.333294  0.333294  -0.000119
0.333270  0.333270  0.333270  -0.000191
0.333247  0.333247  0.333247  -0.000260
0.333235  0.333235  0.333235  -0.000295
0.333233  0.333233  0.333233  -0.000300
5x4
0.333330  0.333330  0.333330  -0.000010
0.333328  0.333328  0.333328  -0.000017
0.333325  0.333325  0.333325  -0.000024
0.333324  0.333324  0.333324  -0.000029
0.333323  0.333323  0.333323  -0.000030

The first three columns are the approximations to the solution, and the remaining ones are the values of the constraints, which approach zero if they are active. Again by symmetry we see that the solution is correct.

In the barrier method we choose a target relative tolerance $\epsilon$ (here we use $\epsilon = 10^{-4}$). We also show the inner loops for the central path with increasing values of the barrier parameter.

As another example, to minimize the above objective function inside a half-sphere, we use the following constraints:

$$x^2 + y^2 + z^2 \leq 1^2 $$

$$x > 0$$

g2 = mkQuadAt (vector [0,0,0]) (-1) (vector [0,0,0]) (diagl[2,2,2])

g3 = mkLin (vector [-1,0,0]) 0
optimize 1E-4 g0 [g2,g3] (matrix 3 []) (vector []) (vector [0.5,0,0])
1x5
0.500000  0.000000  0.000000  -0.750000  -0.500000
6x5
0.726099  0.370086  0.370086  -0.198853  -0.726099
0.699322  0.487516  0.487516  -0.035605  -0.699322
0.596246  0.559690  0.559690  -0.017985  -0.596246
0.588870  0.559310  0.559310  -0.027577  -0.588870
0.583679  0.558495  0.558495  -0.035486  -0.583679
0.582063  0.558022  0.558022  -0.038426  -0.582063
4x5
0.587699  0.569825  0.569825  -0.005208  -0.587699
0.576059  0.576489  0.576489  -0.003477  -0.576059
0.577665  0.575477  0.575477  -0.003956  -0.577665
0.577839  0.575354  0.575354  -0.004037  -0.577839
4x5
0.578511  0.576690  0.576690  -0.000181  -0.578511
0.577927  0.576940  0.576940  -0.000281  -0.577927
0.577529  0.577102  0.577102  -0.000367  -0.577529
0.577409  0.577147  0.577147  -0.000402  -0.577409
4x5
0.577473  0.577281  0.577281  -0.000018  -0.577473
0.577412  0.577308  0.577308  -0.000028  -0.577412
0.577369  0.577325  0.577325  -0.000037  -0.577369
0.577356  0.577330  0.577330  -0.000040  -0.577356

In this case only the quadratic constraint becomes active.

(The initial point must be feasible (we are now testing only phase II), otherwise the process will not finish.)

We can also minimize the same function inside the opposite half-sphere:

$$x^2 + y^2 + z^2 \leq 1^2 $$

$$x < 0$$

g4 = mkLin (vector [1,0,0]) 0
optimize 1E-4 g0 [g2,g4] (matrix 3 []) (vector []) (vector [-0.5,0,0])
1x5
-0.500000  0.000000  0.000000  -0.750000  -0.500000
6x5
-0.045344  0.462607  0.462607  -0.569933  -0.045344
-0.037309  0.674906  0.674906  -0.087613  -0.037309
-0.017136  0.696855  0.696855  -0.028492  -0.017136
-0.022103  0.692108  0.692108  -0.041484  -0.022103
-0.023974  0.688938  0.688938  -0.050155  -0.023974
-0.024151  0.688138  0.688138  -0.052350  -0.024151
4x5
-0.003729  0.704086  0.704086  -0.008512  -0.003729
-0.001878  0.705790  0.705790  -0.003717  -0.001878
-0.002340  0.705371  0.705371  -0.004897  -0.002340
-0.002482  0.705196  0.705196  -0.005391  -0.002482
4x5
-0.000128  0.706995  0.706995  -0.000316  -0.000128
-0.000190  0.706948  0.706948  -0.000450  -0.000190
-0.000236  0.706920  0.706920  -0.000529  -0.000236
-0.000249  0.706914  0.706914  -0.000546  -0.000249
5x5
-0.000010  0.707100  0.707100  -0.000020  -0.000010
-0.000016  0.707095  0.707095  -0.000032  -0.000016
-0.000021  0.707091  0.707091  -0.000046  -0.000021
-0.000024  0.707088  0.707088  -0.000053  -0.000024
-0.000025  0.707087  0.707087  -0.000055  -0.000025

In this case both constraints become active at the solution.

Finally, we can add an equality constraint $z=0$ to work in a half-disc:

optimize 1E-4 g0 [g2,g4] (matrix 3 [0,0,1]) (vector [0]) (vector [-0.5,0,-0.5])
1x5
-0.500000  0.000000  -0.500000  -0.500000  -0.500000
10x5
-0.008480  0.546133  -0.336160  -0.588663  -0.008480
-0.010764  0.972421  -0.198469  -0.014892  -0.010764
-0.014590  0.993978  -0.071449  -0.006689  -0.014590
-0.020164  0.995888  -0.000000  -0.007800  -0.020164
-0.023367  0.992231   0.000000  -0.014932  -0.023367
-0.024289  0.985881   0.000000  -0.027448  -0.024289
-0.024380  0.976066  -0.000000  -0.046702  -0.024380
-0.024380  0.964286  -0.000000  -0.069558  -0.024380
-0.024377  0.955709  -0.000000  -0.086026  -0.024377
-0.024377  0.953226  -0.000000  -0.090766  -0.024377
6x5
-0.003658  0.984073  -0.000000  -0.031587  -0.003658
-0.002958  0.997994  -0.000000  -0.003999  -0.002958
-0.002409  0.996800   0.000000  -0.006383  -0.002409
-0.002491  0.995661  -0.000000  -0.008653  -0.002491
-0.002494  0.995113  -0.000000  -0.009744  -0.002494
-0.002494  0.995035  -0.000000  -0.009899  -0.002494
4x5
-0.000117  0.999598  -0.000000  -0.000804  -0.000117
-0.000179  0.999519   0.000000  -0.000961  -0.000179
-0.000230  0.999501   0.000000  -0.000998  -0.000230
-0.000248  0.999500   0.000000  -0.000999  -0.000248
5x5
-0.000010  0.999981   0.000000  -0.000039  -0.000010
-0.000017  0.999969  -0.000000  -0.000062  -0.000017
-0.000022  0.999957   0.000000  -0.000086  -0.000022
-0.000025  0.999951   0.000000  -0.000098  -0.000025
-0.000025  0.999950   0.000000  -0.000100  -0.000025

general nonlinear functions

Minimize

$$(x-2)^2 + (y-2)^2$$

subject to

$$ y < \cos(x)$$

This function is not convex, so we add more constraints to get a convex feasible region:

$$ x < \frac{\pi}{2} $$

$$ x > -\frac{\pi}{2} $$

f0 = mkQuadAt (vector [2,2]) 0 (vector [0,0]) (diagl [1,1])
nlc = FunAp{..}
  where
    evF (elems -> [x,y]) = y - cos x
    evG (elems -> [x,y]) = vector [sin x, 1]
    evH (elems -> [x,y]) = matrix 2 [ cos x, 0
                                    ,   0  , 0 ]

nlc1 = mkLin (vector [1,0]) (-pi/2)

nlc2 = mkLin (vector [-1,0]) (-pi/2)
optimize 1E-4 f0 [nlc, nlc1, nlc2] (matrix 2 []) (vector []) (vector [0,-1.5])
1x5
0.000000  -1.500000  -2.500000  -1.570796  -1.570796
4x5
0.913424  0.243622  -0.367417  -0.657373  -2.484220
0.963981  0.505977  -0.064277  -0.606816  -2.534777
0.851086  0.592588  -0.066579  -0.719710  -2.421882
0.847392  0.591234  -0.070706  -0.723405  -2.418188
4x5
0.883771  0.629749  -0.004491  -0.687025  -2.454567
0.898622  0.616567  -0.006122  -0.672174  -2.469418
0.901385  0.613479  -0.007045  -0.669411  -2.472182
0.901505  0.613223  -0.007207  -0.669291  -2.472302
4x5
0.905236  0.617058  -0.000441  -0.665560  -2.476033
0.907140  0.615390  -0.000612  -0.663657  -2.477936
0.907506  0.615008  -0.000705  -0.663290  -2.478302
0.907523  0.614977  -0.000722  -0.663273  -2.478320
5x5
0.907907  0.615370  -0.000027  -0.662889  -2.478703
0.908055  0.615237  -0.000044  -0.662742  -2.478851
0.908121  0.615167  -0.000061  -0.662675  -2.478917
0.908132  0.615149  -0.000070  -0.662664  -2.478928
0.908132  0.615147  -0.000072  -0.662664  -2.478929




 Convex Optimization 



























 

 










































































































 -2 
 -1 
 0 
 1 
 2 
 3 

 x 




 -2 
 -1 
 0 
 1 
 2 
 3 







 y 






































 

 








































 constraints 
 objective function 
 central path 





Phase I

The iterative process can be stopped as soon as all constraints are verified (otherwise the optimization tries to find the center of the feasible region, moving to a very far location if it is unbounded)

optimIOk 1E-2 f0 [nlc, nlc1, nlc2] (matrix 2 []) (vector []) (vector [1,2])
2x3
1.000   2.000  10.000
0.706  -6.639   1.305

We add another constraint and to show how the center of the region is found:

phaseIb 1E-2 f0 [nlc, nlc1, nlc2, mkLin (vector [0,-1]) (-1)] (matrix 2 []) (vector []) (vector [1,2])
1x8
1.000  2.000  10.000  -8.540  -10.571  -12.571  -13.000  -10.000
4x8
-0.011  -0.760  1.734  -3.494  -3.315  -3.294  -1.974  -1.734
-0.001   0.065  0.304  -1.239  -1.877  -1.874  -1.369  -0.304
-0.000   0.003  0.094  -1.091  -1.665  -1.665  -1.097  -0.094
-0.000   0.000  0.126  -1.126  -1.697  -1.697  -1.126  -0.126
2x8
-0.000  0.000  0.009  -1.009  -1.580  -1.580  -1.009  -0.009
-0.000  0.000  0.010  -1.010  -1.581  -1.581  -1.010  -0.010




 Convex Optimization 



























 

 









































































 -2 
 -1 
 0 
 1 
 2 
 3 

 x 




 -2 
 -1 
 0 
 1 
 2 
 3 







 y 






































 

 





























 constraints 
 central path 





There is a problem in phase I with the nonconvex constraint: if the starting point is located in a concave region we may get a negative hessian, and hence a negative Newton decrease which stops the optimization. Using the $l_1$-style phase I method we find the solution but the path is somewhat erratic.





 Convex Optimization 



























 

 














































































 -2 
 -1 
 0 
 1 
 2 
 3 

 x 




 -2 
 -1 
 0 
 1 
 2 
 3 







 y 






































 

 





























 constraints 
 central path 





With convex constraints phase I is straightforward:

constr = [mkLin (vector [cos t, sin t]) (-1) | t <- ts]
  where
    ts = init $ toList (linspace 20 (0, 2*pi))
optimIOk 1E-2 f0 constr (matrix 2 []) (vector []) (vector [2,2])
3x3
 2.000   2.000  10.000
-1.252  -1.252   1.484
-0.033  -0.033   0.483




 Convex Optimization 



























 

 










































































































































































































































 -2 
 -1 
 0 
 1 
 2 
 3 

 x 




 -2 
 -1 
 0 
 1 
 2 
 3 







 y 






































 

 


















 central path 






A Ruiz / 2014-07-01
powered by hmatrix