dikin.m

Go to the documentation of this file.

%> \file
%>  Algorithm 17.3:  Dikin's method.   Implementation of algorithm 17.3 of \cite Bier15-book
%>
%> @note Tested with \ref run1703.m
%> @note Tested with \ref run1607dikin.m
%>
%> @author Michel Bierlaire
%> @date Sun Mar 22 16:39:13 2015
%> @ingroup Algorithms
%> @ingroup chap17

%> Applies Dikin's  method to solve \f[\min_x c^Tx  \f] subject to \f[Ax = b\f] and \f[ x \geq 0 \f].
%> @param A the constraint matrix
%> @param b the constraint right hand side
%> @param c the cost vector for the objective function
%> @param x0 starting point
%> @param eps algorithm stops if \f$\|d_k\| \leq \varepsilon \f$. 
%> @param beta parameter such that 0 < beta < 1 (default: 0.9)
%> @param maxiter maximum number of iterations (default: 100)
%> @return [solution,iteres,niter] 
%> @return solution: local minimum of the function
%> @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.
%> @return unbounded If 0, the problem is unbounded. If different from 0 the problem is bounded. 
function [solution, iteres, niter, unbounded] =  dikin(A,b,c,x0,eps, beta=0.9, maxiter=100)
  [m,n] = size(A) ;
  if (m != size(b))
    error("The number of rows in A and b do not match") ;
  endif
  if (n != size(c))
    error("The number of columns in A and do not match the size of c") ;
  endif
  xk = x0 ;
  k = 0 ;
  iteres = zeros(1+ maxiter,n+2) ;
  do 
    Hinv = diag(xk)^2 ;
    lambda = (A * Hinv * A') \ (A * Hinv * c) ;
    dk = - Hinv * (c - A' * lambda) ;
    iteres(k+1,:) = [xk' c'*xk  norm(dk) ] ;
    alphamax = realmax ;
    unbounded = 1 ;
    for i=1:n
      if (dk(i) < 0)
    unbounded = 0 ;
    alphai = -xk(i) / dk(i) ;
    if (alphai < alphamax)
      alphamax = alphai ;
    endif
      endif
    endfor
    if (unbounded == 0)
      xk = xk + beta * alphamax * dk ;
      k = k + 1 ;
    endif
  until (norm(dk) <= eps || k >= maxiter || unbounded != 0)
  solution = xk ;
  niter = k+1 ;
endfunction