conjugateGradient.m

Go to the documentation of this file.

%> \file
%> Algorithm 9.2: Conjugate gradient algorithm for quadratic problems. Implementation of algorithm 9.2 of \cite Bier15-book
%> 
%> @author <a href="http://people.epfl.ch/michel.bierlaire">Michel Bierlaire</a>
%> @date Sat Apr  5 23:32:26 2014
%> @ingroup Algorithms
%> @ingroup chap09

%> @note Tested with \ref run0908cg.m
%> @note Called by \ref newtonLocalQuadratic

%> Applies the conjugate gradient method to solve  \f[\min_x \frac{1}{2} x^T Q x + b^T x\f] where \f$Q \in \mathbb{R}^n\times\mathbb{R}^n \f$ and \f$b \in \mathbb{R}^n\f$.
%> @param Q matrix of size \f$n \times n \f$.
%> @param b vector of size \f$n\f$.
%> @param x0 starting point
%> @param printlevel if different from 0, informations are printed at each iteration (Default: 0)
%> @return [D, solution] 
%> @return D: matrix gathering all directions generated by the algorithms
%> @return  solution: optimal solution
function [D, solution] = conjugateGradient(Q,b,x0,printlevel=0)
  n = size(x0,1) ;
  xk = x0 ;
  gk = Q * xk + b ;
  dk = -gk ;
  D = dk ;
  betak = 0 ;
  for k=1:n
    xprint = xk ;
    gprint = gk ;
    dprint = dk ;
    if (printlevel != 0)
      printf("%3d\t%+10.5e\t%+10.5e\t%+10.5e",k,xprint(1),gprint(1),dprint(1)) ;
    endif
    denom = dk' * Q * dk ;
    if (denom <= 0)
      error("The matrix must de positive definite") ;
    endif
    alphak = - (dk' * gk) / denom;
    xk = xk + alphak * dk ;
    gkp1 = Q * xk + b ;
    if (printlevel != 0)
      if (betak == 0)
    printf("\t%+10.5e\n",alphak) ;
      else
    printf("\t%+10.5e\t%+10.5e\n",alphak,betak) ;
      endif
      for i=2:n
    printf("\t%+10.5e\t%+10.5e\t%+10.5e\t\t\n",xprint(i),gprint(i),dprint(i)) ;
      endfor
    endif
    betak = (gkp1' * gkp1) / (gk' * gk) ;
    dk = -gkp1 + betak * dk ;
    D = [D  dk] ;
    gk = gkp1 ;
%    printf("%d %f\n",k,norm(gk)) ;
  endfor  
  if (printlevel != 0)
    printf("%3d\t%+10.5e\t%+10.5e\t\t\t\n",n+1,xk(i),gk(i)) ;
    for i=2:n
      printf("\t%+10.5e\t%+10.5e\t\t\t\n",xk(i),gk(i)) ;
    endfor  
  endif
  solution = xk ;
endfunction