Optimization: principles and algorithms, by Michel Bierlaire
truncatedConjugateGradient.m
Go to the documentation of this file.
1 %> \file
2 %> Algorithm 12.3: Truncated conjugate gradients method to find an approximation of the trust region subproblem. Implementation of algorithm 12.3 of \cite Bier15-book
3 %>
4 %> @author Michel Bierlaire
5 %> @date Sat Mar 21 15:49:09 2015
6 %> @ingroup Algorithms
7 %> @ingroup chap12
8 
9 %> Given a current iterate \f$\widehat{x}\f$, consider the trust region subproblem \f[ \min_d d^T \nabla f(\widehat{x}) + \frac{1}{2} d^T \nabla^2 f(\widehat{x}) d \f] subject to \f[\|d\|_2 \leq \Delta. \f] Apply the truncated conjugate gradients method to find an approximate solution.
10 
11 %> @note Tested with \ref run1203tcg.m
12 %> @note Called by \ref newtonTrustRegion
13 %> @note Called by \ref symmetricRankOne
14 
15 %> @param g the gradient \f$ \nabla f(\widehat{x})\f$
16 %> @param H the hessian \f$ \nabla^2 f(\widehat{x})\f$
17 %> @param delta the radius \f$ \Delta\f$ of the trust region
18 %> @return [dstar,type]
19 %> @return dstar: the approximate solution
20 %> @return type: the type of outcome of the method, based on the following code:
21 %> - type 1: convergence
22 %> - type 2: out of the trust region
23 %> - type 3: negative curvature detected
24 function [step,type] = truncatedConjugateGradient(g,H,delta)
25  tol = 1.0e-6 ;
26  n = size(g,1) ;
27  xk = zeros(n,1) ;
28  gk = g ;
29  dk = -gk ;
30  for k=1:n
31  if (dk' * H * dk <= 0)
32  type = 3 ;
33  a = dk' * dk ;
34  b = 2 * xk' * dk ;
35  c = xk' * xk - delta * delta ;
36  rho = b * b - 4 * a * c ;
37  step = xk + ((-b + sqrt(rho)) / (2*a)) * dk ;
38  return ;
39  endif
40  alphak = - (dk' * gk) / (dk' * H * dk) ;
41  xkp1 = xk + alphak * dk ;
42  if (norm(xkp1) > delta)
43  type = 2 ;
44  a = dk' * dk ;
45  b = 2 * xk' * dk ;
46  c = xk' * xk - delta * delta ;
47  rho = b * b - 4 * a * c ;
48  step = xk + ((-b + sqrt(rho)) / (2*a)) * dk ;
49  return ;
50  endif
51  xk = xkp1 ;
52  gkp1 = H * xk + g ;
53  betak = (gkp1' * gkp1) / (gk' * gk) ;
54  dk = -gkp1 + betak * dk ;
55  gk = gkp1 ;
56  if (norm(gkp1) <= tol)
57  type = 1 ;
58  step = xk ;
59  return ;
60  endif
61  endfor
62  type = 1 ;
63  step = xk ;
64 endfunction
function symmetricRankOne(in obj, in x0, in delta0, in eps, in tr)
Applies SR1 algorithm with trust region to solve where . The parameters of the method are taken from...
function truncatedConjugateGradient(in g, in H, in delta)
Given a current iterate , consider the trust region subproblem subject to Apply the truncated conju...
function newtonTrustRegion(in obj, in x0, in delta0, in eps, in tr)
Applies Newton&#39;s algorithm with trust region to solve where . The parameters of the method are taken...
Copyright 2015-2016 Michel Bierlaire