A Ruiz / 2014-07-01
powered by 
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
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
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.
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
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
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
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.
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