Optimization: principles and algorithms, by Michel Bierlaire
dikin.m
Go to the documentation of this file.
1 %> \file
2 %> Algorithm 17.3: Dikin's method. Implementation of algorithm 17.3 of \cite Bier15-book
3 %>
4 %> @note Tested with \ref run1703.m
5 %> @note Tested with \ref run1607dikin.m
6 %>
7 %> @author Michel Bierlaire
8 %> @date Sun Mar 22 16:39:13 2015
9 %> @ingroup Algorithms
10 %> @ingroup chap17
11 
12 %> Applies Dikin's method to solve \f[\min_x c^Tx \f] subject to \f[Ax = b\f] and \f[ x \geq 0 \f].
13 %> @param A the constraint matrix
14 %> @param b the constraint right hand side
15 %> @param c the cost vector for the objective function
16 %> @param x0 starting point
17 %> @param eps algorithm stops if \f$\|d_k\| \leq \varepsilon \f$.
18 %> @param beta parameter such that 0 < beta < 1 (default: 0.9)
19 %> @param maxiter maximum number of iterations (default: 100)
20 %> @return [solution,iteres,niter]
21 %> @return solution: local minimum of the function
22 %> @return iteres: sequence of iterates generated by the algorithm. It contains n+2 columns. Columns 1:n contains the value of the current iterate. Column n+1 contains the value of the objective function. Column n+2 contains the value of the norm of the gradient. It contains maxiter rows, but only the first niter rows are meaningful.
23 %> @return unbounded If 0, the problem is unbounded. If different from 0 the problem is bounded.
24 function [solution, iteres, niter, unbounded] = dikin(A,b,c,x0,eps, beta=0.9, maxiter=100)
25  [m,n] = size(A) ;
26  if (m != size(b))
27  error("The number of rows in A and b do not match") ;
28  endif
29  if (n != size(c))
30  error("The number of columns in A and do not match the size of c") ;
31  endif
32  xk = x0 ;
33  k = 0 ;
34  iteres = zeros(1+ maxiter,n+2) ;
35  do
36  Hinv = diag(xk)^2 ;
37  lambda = (A * Hinv * A') \ (A * Hinv * c) ;
38  dk = - Hinv * (c - A' * lambda) ;
39  iteres(k+1,:) = [xk' c'*xk norm(dk) ] ;
40  alphamax = realmax ;
41  unbounded = 1 ;
42  for i=1:n
43  if (dk(i) < 0)
44  unbounded = 0 ;
45  alphai = -xk(i) / dk(i) ;
46  if (alphai < alphamax)
47  alphamax = alphai ;
48  endif
49  endif
50  endfor
51  if (unbounded == 0)
52  xk = xk + beta * alphamax * dk ;
53  k = k + 1 ;
54  endif
55  until (norm(dk) <= eps || k >= maxiter || unbounded != 0)
56  solution = xk ;
57  niter = k+1 ;
58 endfunction
function dikin(in A, in b, in c, in x0, in eps, in beta)
Applies Dikin&#39;s method to solve subject to and .
Copyright 2015-2016 Michel Bierlaire