iconEuler Examples

Affine Scaling

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

Affine Scaling

Affine Scaling

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
Nowe we find the projection to

Affine Scaling

>v=c'-A'.y
           -1 
            0 
            1 
We can walk in this direction, untol 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.

Affine Scaling

This yields an equivalent problem.

Affine Scaling

Affine Scaling

The trick is that x is now the vector (1,1,1), which is central in the
positive quadrant, and promisses 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.
>load affinescaling;
Here are the details of the implemenation.
>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.
>load povray;
>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])));
>povend();

Affine Scaling

Another Problem

As a second simple example, we demonstrate

Affine Scaling

Affine Scaling

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

Affine Scaling

Affine Scaling

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);
>plot2d(xopt[1],xopt[2],>points,>add);
>insimg;

Affine Scaling

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.
>plot2d(X[1],X[2],>add); insimg;

Affine Scaling

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)
15
>plot2d(X[1],X[2],color=blue,>add); insimg;

Affine Scaling

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 Homepage