1 /****************************************************************************/
2 /* This file is part of FreeFEM. */
3 /* */
4 /* FreeFEM is free software: you can redistribute it and/or modify */
5 /* it under the terms of the GNU Lesser General Public License as */
6 /* published by the Free Software Foundation, either version 3 of */
7 /* the License, or (at your option) any later version. */
8 /* */
9 /* FreeFEM is distributed in the hope that it will be useful, */
10 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
11 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
12 /* GNU Lesser General Public License for more details. */
13 /* */
14 /* You should have received a copy of the GNU Lesser General Public License */
15 /* along with FreeFEM. If not, see <http://www.gnu.org/licenses/>. */
16 /****************************************************************************/
17 // SUMMARY : Conjugate gradient class
18 // LICENSE : LGPLv3
19 // ORG : LJLL Universite Pierre et Marie Curie, Paris, FRANCE
20 // AUTHORS : ...
21 // E-MAIL : ...
22
23 // Definition of the conjugate gradient class
24 // Non-linear conjugate gradient algorithm
25 // author: Wenceslau Gouveia
26 // modified: H. Lydia Deng, 02/23/94, 03/14/94
27
28 // .NAME ConjugateGradient class
29 // .LIBRARY Base
30 // .HEADER Optimization Algorithms
31 // .INCLUDE defs.hpp
32 // .FILE CG.hpp
33
34 // .SECTION Description
35 // .B ConjugateGradient()
36 // The conjugate gradient procedure implemented in this object is an extension of
37 // the conjugate gradient used in linear system solving to handle non quadratic
38 // objective functions. Such extensions amount basically to the inclusion of a
39 // line search step and a modification in the computation of the conjugate directions.
40 // Details can be found in the classic (and tough!) book Practical Optimization by
41 // Powell.
42 //
43 // .SECTION Description
44 // Public Operations
45 // Constructors:
46 // ConjugateGradient(LineSearch* ls, int iter, double tol)
47 // Here:
48 // ls: Defines the line search used. At the present version you should use the
49 // CubicLineSearch procedure.
50 // iter: Maximum number of iterations
51 // tol: Minimum accepted module of the gradient at optimum solution
52 // Methods:
53 // Model<>ConjugateGradient::optimizer(Model<double>&)
54 // Here:
55 // model0: Initial model for the conjugate gradient procedure.
56 //
57 // The optimum model is returned by the function.
58 //
59 // .SECTION Caveats
60 // This procedure requires the derivatives of the objective function. At the present version the
61 // COOOL library cannot handle analytic derivatives provided by the user. Here the derivatives
62 // are computed numerically by a centered finite differencing scheme in a automatic fashion.
63 // The hability to support user-provided derivatives will hopefully be
64 // implemented in the near future.
65
66 #ifndef CONJUGATE_GRADIENT_HH
67 #define CONJUGATE_GRADIENT_HH
68
69 #include "Optima.hpp"
70 #include "defs.hpp"
71 #include "Param.hpp"
72
73 // Nonlinear Conjugate Gradient
74 /*
75 This conjugate gradient procedure implemented in this
76 object is an extension of the conjugate gradient
77 used in linear system solving to handle non quadratic
78 objective functions. Such extensions amount basically
79 to the inclusion of a line search step and a
80 modification in the computation of the conjugate directions.
81 Details can be found in the classic (and tough!) book
82 Practical Optimization by Powell.
83 */
84
85 template< class LS >
86 class ConjugateGradient : public Optima< LS > {
87 typedef typename LS::Real Real;
88 typedef typename LS::Param Param;
89 typedef typename LS::Vect Vect;
90 typedef typename LS::VMat VMat;
91 typedef LS LineSearch;
92
93 public:
94 // a constructor
95 ConjugateGradient(LS* ls, // Pointer to the line-search object
96 int iter, // Maximum number of iterations
97 Real tol, // Minimum accepted gradient at optimum solution
98 int verb = 0 // verbose or quiet
99 );
100 /* At the present version, you must use the CubicLineSearch procedure */
101 // a constructor
~ConjugateGradient()102 ~ConjugateGradient( ) { ; }
103
104 // Conjugate gradient search starting from m0, returns an optimum Model
105 Param optimizer(Param& m0);
106 };
107
108 template< class LS >
ConjugateGradient(LS * p,int it,Real eps,int verb)109 ConjugateGradient< LS >::ConjugateGradient(LS* p, int it, Real eps, int verb) : Optima< LS >(verb) {
110 ls = p;
111 iterMax = it;
112 tol = eps;
113 iterNum = 0;
114 }
115
116 // Needs :
117 // A scalar product on Vect: (u,v)
118 // A Real Vect product: operator*
119 // A Real Vect product: operator*=
120 // A Vect difference: operator-
121 template< class LS >
optimizer(Param & model0)122 ConjugateGradient< LS >::Param ConjugateGradient< LS >::optimizer(Param& model0) {
123 // Reset the residue history for every new optimizer
124 iterNum = 0;
125 isSuccess = 0;
126 if (residue != NULL) {
127 delete residue;
128 residue = new list< Real >;
129 }
130
131 int n = model0.size( );
132 Param model1(model0); // new model
133 Vect search(n); // search direction
134 Vect g0(n); // old gradient vector
135 Vect g1(n); // new gradient vector
136 double beta; // beta parameter
137 double lambda = .025; // line search parameter
138 double descent = 0.; // descent direction
139
140 // Beginning iterations
141 g0 = *(ls->gradient(model0));
142
143 // Check the gradient, in case the initial model is the optimal
144 Real err = (Real)sqrt((g0, g0));
145 if (isVerbose) cout << "Initial residue : " << err << endl;
146 appendResidue(err); // residual
147 if (err < tol) {
148 if (isVerbose) cout << "Initial guess was great!" << endl;
149 isSuccess = 1;
150 return model0;
151 }
152
153 // Considering first iteration
154 search = -1. * g0;
155 descent = (search, g0);
156
157 // We use CubicLineSearch
158 // cerr << "Line Search" << endl;
159 // cerr << "model0 " << model0;
160 // cerr << "search " << (Param)search; // we cast only for display
161 // cerr << "descent " << descent << endl;
162 // cerr << "lambda " << lambda << endl;
163 // cerr << endl;
164
165 model1 = ls->search(model0, search, descent, lambda);
166 g1 = *(ls->gradient(model1)); // Gradient at new model
167 err = (Real)sqrt((g1, g1));
168 if (isVerbose)
169 cout << "Iteration (0) : "
170 << "current value of the objective function: " << ls->currentValue( )
171 << "\t current residue: " << err << endl;
172 appendResidue(err); // residual
173
174 iterNum = 0;
175 Real temp;
176 do {
177 iterNum++;
178
179 temp = 1. / ((g0, g0));
180 beta = ((g1 - g0), g1);
181 beta *= temp; // computation Polak & Ribiere
182
183 search = beta * search - g1; // search direction
184
185 descent = (search, g1); // descent
186 if (descent > 0.) {
187 if (isVerbose) cout << "Reset search direction to gradient vector!" << endl;
188 search = -1. * g1;
189 descent = (search, g1);
190 }
191
192 model0 = model1;
193 g0 = g1; // save the old model and gradient before new search
194
195 // We use CubicLineSearch
196 // cerr << "Line Search" << endl;
197 // cerr << "model0 " << model0;
198 // cerr << "search " << (Param)search;
199 // cerr << "descent " << descent << endl;
200 // cerr << "lambda " << lambda << endl;
201 // cerr << endl;
202
203 model1 = ls->search(model0, search, descent, lambda); // line search
204 g1 = *(ls->gradient(model1));
205
206 err = (Real)sqrt((g1, g1));
207 if (isVerbose)
208 cout << "Iteration (" << iterNum << ") : "
209 << "current value of the nrj : " << ls->currentValue( ) << "\t current residue : " << err
210 << endl;
211 appendResidue(err); // residual
212
213 } while (finalResidue( ) > tol && iterNum < iterMax); // stopping criterion
214
215 if (finalResidue( ) <= tol) isSuccess = 1;
216
217 return (model1); // hopefully answer
218 }
219
220 #endif
221