Optimization: principles and algorithms, by Michel Bierlaire
longSteps.m
Go to the documentation of this file.
1 %> \file
2 %> Algorithm 18.4: long steps interior point algorithm. Implementation of algorithm 18.4 of \cite Bier15-book
3 %>
4 %> @note Tested with \ref run1807longSteps.m
5 %>
6 %> @author Michel Bierlaire
7 %> @date Mon Mar 23 14:52:36 2015
8 %> @ingroup Algorithms
9 %> @ingroup chap18
10 
11 %> Applies the long step interior point algorithm to solve \f[\min_x c^Tx \f] subject to \f[Ax = b\f] and \f[ x \geq 0 \f].
12 %> @param A the constraint matrix
13 %> @param b the constraint right hand side
14 %> @param c the cost vector for the objective function
15 %> @param x0 starting primal point (nx1)
16 %> @param lambda0 starting dual point for equality constraints (mx1)
17 %> @param mu0 starting dual point for inequality constraints (nx1)
18 %> @param eps algorithm stops if \f$\|d_k\| \leq \varepsilon \f$.
19 %> @param gamma parameter such that gamma > 0 (default: 1.0e-3)
20 %> @param sigma parameter such that 0 < sigma < 1 (default: 0.1)
21 %> @param maxiter maximum number of iterations (default: 100)
22 %> @return [x,lambda,mu]
23 %> @return x: primal solution
24 %> @return lambda: dual solution for equality constraints
25 %> @return mu: dual solution for inequality constraints
26 function [x,lambda,mu] = longSteps(A,b,c,x0,lambda0,mu0,eps,gamma=1.0e-3,sigma=0.1,maxiter=100)
27 
28  [m,n] = size(A) ;
29  if (m != size(b))
30  error("The number of rows in A and b do not match") ;
31  endif
32  if (n != size(c))
33  error("The number of columns in A and do not match the size of c") ;
34  endif
35  if (n != size(x0))
36  error("The number of columns in A and do not match the size of x0") ;
37  endif
38  if (m != size(lambda0))
39  error("The number of rows in A and do not match the size of lambda0") ;
40  endif
41  if (n != size(mu0))
42  error("The number of columns in A and do not match the size of mu0") ;
43  endif
44 
45  primalcons = norm(A*x0 - b) ;
46  if (primalcons > eps)
47  error("Initial point is not primal feasible") ;
48  endif
49  dualcons = norm(A' * lambda0 + mu0 - c) ;
50  if (dualcons > eps)
51  error("Initial point is not dual feasible") ;
52  endif
53 
54  if (min(x0) <= 0)
55  error("x0 must be strictly positive")
56  endif
57  if (min(mu0) <= 0)
58  error("mu0 must be strictly positive")
59  endif
60 
61  x = x0 ;
62  lambda = lambda0 ;
63  mu = mu0 ;
64  e = ones(n,1) ;
65  k=0 ;
66  do
67  nu = (mu' * x ) / n ;
68  S = diag(mu) ;
69  X = diag(x) ;
70  voisinage = max(x .* mu) / nu;
71  if (voisinage < gamma)
72  error("Iterate not in the vicinity of the central path: %f",voisinage)
73  endif
74  if (k==0)
75  printf("%d\t%e\t\t%e\n",k,voisinage,nu) ;
76  else
77  printf("%d\t%e\t%e\t%e\t%e\n",k,voisinage,norm(alpha*d),nu,alpha);
78  endif
79  Q = [ A zeros(m,m) zeros(m,n) ; zeros(n,n) A' eye(n,n) ; S zeros(n,m) X] ;
80  lefthand = [ zeros(m+n,1) ; -X*S*e + nu * sigma * e] ;
81  d = Q \ lefthand ;
82  alpha = 1.0 ;
83  ok = 0 ;
84  while (ok != 1)
85  # Candidate
86  cx = x + alpha * d(1:n) ;
87  clambda = lambda + alpha * d(n+1:n+m) ;
88  cmu = mu + alpha * d(n+m+1:n+m+n) ;
89  if (min(cx) <= 0 || min(cmu) <= 0)
90  ok = 0 ;
91  alpha = alpha / 2.0 ;
92  else
93  cnu = (cmu' * cx ) / n ;
94  if (cnu <= eps)
95  # Solution found
96  x = cx ;
97  lambda = clambda ;
98  mu = cmu ;
99  return ;
100  endif
101  cvoisinage = max(cx .* cmu) / cnu;
102  if (cvoisinage >= gamma)
103  x = cx ;
104  lambda = clambda ;
105  mu = cmu ;
106  ok = 1 ;
107  else
108  ok = 0 ;
109  alpha = alpha / 2.0 ;
110  endif
111  endif
112  endwhile
113  k=k+1 ;
114  until (nu <= eps || k >= maxiter)
115 endfunction
116 
function longSteps(in A, in b, in c, in x0, in lambda0, in mu0, in eps, in gamma)
Applies the long step interior point algorithm to solve subject to and .
Copyright 2015-2016 Michel Bierlaire