﻿ Euler Math Toolbox - Examples

# Optimization by Affine Scaling

by R. Grothmann

We demonstrate a method be Vanderbei and other authors from 1985-1986 for optimization. There is an implemenation of this method in Euler, which may be useful.

First, we demonstrate inner methods for the example

The solution is [3,0,0].

>shortformat; A=[1,1,1]; b=[3]; c=[1,2,3]; x=[1;1;1];

If we start with this x, we get the target value 6.

>c.x
6

Inner methods proceed from x>0 to x>0 using a direction of descent. To find such a direction, we project c to the kernel of A.

First we project c to the image of A', which is perpendicular to the kernel of A.

>y=fit(A',c')
2

In our case, a Gram-Schmidt method can be used.

>(A.c')/(A.A')
2

Now we find the projection to

>v=c'-A'.y
-1
0
1

We can walk in this direction, until we meet the boundary. Instead we walk only half the way, and stay within the positive quadrant.

>xnew=x-v/2
1.5
1
0.5

The new value is better. Remember, that we minimize!

>c.xnew
5

For the next step, we use the transformation matrix.

This yields an equivalent problem.

The trick is that x is now the vector (1,1,1), which is central in the positive quadrant, and promises a quick descent.

>x=xnew; d=x'; tA=A*d; tc=d*c;

We do the same stuff as above in the equivalent problem.

>y=fit(tA',tc');

We get a direction of descent.

>v=tc'-tA'.y
-0.642857
0.571429
0.785714

We can go as far as the following value.

>lambda=max(1/v')
1.75

But we go only half that far. Note that we are walking away from the (1,1,1) vector.

>txnew=1-lambda*v/2
1.5625
0.5
0.3125

Transform back.

>xnew=txnew*d'
2.34375
0.5
0.15625

The value is indeed smaller.

>c.txnew
3.5

This process is called affine scaling.

We load the file with the function. It is not yet part of the Euler core.

Here are the details of the implementation.

>type affinescaling
function affinescaling (A: real, b: column real, c: vector real,  ..
gamma: number positive, x: none column real, history, infinite ..
)

## Default for gamma : 0.9
## Default for x : none
## Default for history : 0
## Default for infinite : 4.5036e+011

n=cols(A);
if x==none then
u=b-sum(A);
A=A|u;
x=ones(cols(A),1);
c=c|infinite;
endif;
z=c.x;
if history then X=x; endif;
repeat
d=ones(1,cols(A));
d=x';
tA=A*d; tc=c*d;
y=fit(tA',tc');
v=tc'-tA'.y;
while !all(v~=0);
if all(v<=0) then
error("Problem unbounded");
endif;
i=nonzeros(v>0);
lambda=gamma/max(v[i]');
tx=x/d'-lambda*v;
x=d'*tx;
znew=c.x;
if history then X=X|x; endif;
while znew<z;
z=znew;
end
if history then return {x[1:n],X[1:n]};
else return x[1:n];
endif;
endfunction
>{xopt,X}=affinescaling(A,b,c,gamma=0.5,x=[1,1,1]',>history); xopt,
3
0
0

We can plot the points in 3D using Povray. This works only, if Povray is installed and its "bin" directory in the environment variable PATH.

>povstart(distance=9,center=[0,0,1],angle=140°);
>loop 1 to 8; writeln(povpoint(X[:,#]',povlook(red),0.05)); end;
>writeAxes(-0.5,3,-0.5,3,-0.5,3);
>writeln(povplane([1,1,1],3,povlook(blue,0.5),povbox([0,0,0],[4,4,4],look="")));
>povend();

# Another Problem

As a second simple example, we demonstrate

This writes as.

>A=[1,1;4,5]|id(2)
1         1         1         0
4         5         0         1
>b=[1000;4500]
1000
4500
>c=-[5,6,0,0]
[-5,  -6,  0,  0]

The Simplex method has no problems at all.

>simplex(A,b,c,eq=0,>min,>check)
500
500
0
0

The algorithm works too. The result is only an approximation, however. We could start there to find an nearby corner.

>affinescaling(A,b,c)
499.989
500.01
0
0

# Ellipse

Another example is

The angle runs from 0 to 90 degrees. The admissible set is an ellipse.

>phi=linspace(0,pi/2,50)'; A=cos(phi)|2*sin(phi);
>b=ones(rows(A),1);
>c=[1,1];

Using the simplex method.

>xopt=simplex(A,b,c,>max,>check)
0.898138
0.219997

Let us plot the feasible set.

>P=feasibleArea(A,b);
>plot2d(P[1],P[2],a=-0.2,b=1,c=-0.2,d=1,>filled);

Now we fill with slack variables, and use affine scaling.

>tA=A|id(rows(A));
>tc=-c|zeros(1,rows(A));

We choose an inner point by computing all slack variables.

>x=[0.1,0.1]'; x=x_(b-A.x);
>{xopt,X}=affinescaling(tA,b,tc,0.5,x,>history);

We kill the slack variables, and get the same result.

>xopt[1:2]
0.898138
0.219997

But the path the algorithm takes is quite different from the simplex method, which goes around the vertices at the boundary.

Here is the number of steps the algorithm took.

>cols(X)
29

With gamma=0.9, the algorithm takes less steps.

>{xopt,X}=affinescaling(tA,b,tc,0.9,x,>history);
>cols(X)
17

# The Dual Problem

It is better to use the dual problem here, since the matrix is only 2x51.

>lopt=affinescaling(A',c',b');

We search for the two active lambdas.

>{xs,js}=sort(-lopt); j=js[1:cols(A)]
16
15

Then we solve for the primal solution.

>x=A[j]\b[j]
0.898138
0.219997

Examples