1 #include <math.h>
2 #include <stdio.h>
3 #include <string.h>
4 #include <stdarg.h>
5 #include "newton.h"
6
7 #ifndef min
min(T x,T y)8 template <class T> static inline T min(T x,T y) { return (x<y)?x:y; }
9 #endif
10
11 #ifndef max
max(T x,T y)12 template <class T> static inline T max(T x,T y) { return (x>y)?x:y; }
13 #endif
14
15 #ifdef __cplusplus
16 extern "C" {
17 #endif
18
19 extern double dnrm2_(int *, double *, int *);
20 extern double ddot_(int *, double *, int *, double *, int *);
21 extern int daxpy_(int *, double *, double *, int *, double *, int *);
22 extern int dscal_(int *, double *, double *, int *);
23
24 #ifdef __cplusplus
25 }
26 #endif
27
default_print(const char * buf)28 static void default_print(const char *buf)
29 {
30 fputs(buf,stdout);
31 fflush(stdout);
32 }
33
34 // On entry *f must be the function value of w
35 // On exit w is updated and *f is the new function value
linesearch_and_update(double * w,double * s,double * f,double * g,double alpha)36 double function::linesearch_and_update(double *w, double *s, double *f, double *g, double alpha)
37 {
38 double gTs = 0;
39 double eta = 0.01;
40 int n = get_nr_variable();
41 int max_num_linesearch = 20;
42 double *w_new = new double[n];
43 double fold = *f;
44
45 for (int i=0;i<n;i++)
46 gTs += s[i] * g[i];
47
48 int num_linesearch = 0;
49 for(num_linesearch=0; num_linesearch < max_num_linesearch; num_linesearch++)
50 {
51 for (int i=0;i<n;i++)
52 w_new[i] = w[i] + alpha*s[i];
53 *f = fun(w_new);
54 if (*f - fold <= eta * alpha * gTs)
55 break;
56 else
57 alpha *= 0.5;
58 }
59
60 if (num_linesearch >= max_num_linesearch)
61 {
62 *f = fold;
63 return 0;
64 }
65 else
66 memcpy(w, w_new, sizeof(double)*n);
67
68 delete [] w_new;
69 return alpha;
70 }
71
info(const char * fmt,...)72 void NEWTON::info(const char *fmt,...)
73 {
74 char buf[BUFSIZ];
75 va_list ap;
76 va_start(ap,fmt);
77 vsprintf(buf,fmt,ap);
78 va_end(ap);
79 (*newton_print_string)(buf);
80 }
81
NEWTON(const function * fun_obj,double eps,double eps_cg,int max_iter)82 NEWTON::NEWTON(const function *fun_obj, double eps, double eps_cg, int max_iter)
83 {
84 this->fun_obj=const_cast<function *>(fun_obj);
85 this->eps=eps;
86 this->eps_cg=eps_cg;
87 this->max_iter=max_iter;
88 newton_print_string = default_print;
89 }
90
~NEWTON()91 NEWTON::~NEWTON()
92 {
93 }
94
newton(double * w)95 void NEWTON::newton(double *w)
96 {
97 int n = fun_obj->get_nr_variable();
98 int i, cg_iter;
99 double step_size;
100 double f, fold, actred;
101 double init_step_size = 1;
102 int search = 1, iter = 1, inc = 1;
103 double *s = new double[n];
104 double *r = new double[n];
105 double *g = new double[n];
106
107 const double alpha_pcg = 0.01;
108 double *M = new double[n];
109
110 // calculate gradient norm at w=0 for stopping condition.
111 double *w0 = new double[n];
112 for (i=0; i<n; i++)
113 w0[i] = 0;
114 fun_obj->fun(w0);
115 fun_obj->grad(w0, g);
116 double gnorm0 = dnrm2_(&n, g, &inc);
117 delete [] w0;
118
119 f = fun_obj->fun(w);
120 fun_obj->grad(w, g);
121 double gnorm = dnrm2_(&n, g, &inc);
122 info("init f %5.3e |g| %5.3e\n", f, gnorm);
123
124 if (gnorm <= eps*gnorm0)
125 search = 0;
126
127 while (iter <= max_iter && search)
128 {
129 fun_obj->get_diag_preconditioner(M);
130 for(i=0; i<n; i++)
131 M[i] = (1-alpha_pcg) + alpha_pcg*M[i];
132 cg_iter = pcg(g, M, s, r);
133
134 fold = f;
135 step_size = fun_obj->linesearch_and_update(w, s, &f, g, init_step_size);
136
137 if (step_size == 0)
138 {
139 info("WARNING: line search fails\n");
140 break;
141 }
142
143 fun_obj->grad(w, g);
144 gnorm = dnrm2_(&n, g, &inc);
145
146 info("iter %2d f %5.3e |g| %5.3e CG %3d step_size %4.2e \n", iter, f, gnorm, cg_iter, step_size);
147
148 if (gnorm <= eps*gnorm0)
149 break;
150 if (f < -1.0e+32)
151 {
152 info("WARNING: f < -1.0e+32\n");
153 break;
154 }
155 actred = fold - f;
156 if (fabs(actred) <= 1.0e-12*fabs(f))
157 {
158 info("WARNING: actred too small\n");
159 break;
160 }
161
162 iter++;
163 }
164
165 if(iter >= max_iter)
166 info("\nWARNING: reaching max number of Newton iterations\n");
167
168 delete[] g;
169 delete[] r;
170 delete[] s;
171 delete[] M;
172 }
173
pcg(double * g,double * M,double * s,double * r)174 int NEWTON::pcg(double *g, double *M, double *s, double *r)
175 {
176 int i, inc = 1;
177 int n = fun_obj->get_nr_variable();
178 double one = 1;
179 double *d = new double[n];
180 double *Hd = new double[n];
181 double zTr, znewTrnew, alpha, beta, cgtol, dHd;
182 double *z = new double[n];
183 double Q = 0, newQ, Qdiff;
184
185 for (i=0; i<n; i++)
186 {
187 s[i] = 0;
188 r[i] = -g[i];
189 z[i] = r[i] / M[i];
190 d[i] = z[i];
191 }
192
193 zTr = ddot_(&n, z, &inc, r, &inc);
194 double gMinv_norm = sqrt(zTr);
195 cgtol = min(eps_cg, sqrt(gMinv_norm));
196 int cg_iter = 0;
197 int max_cg_iter = max(n, 5);
198
199 while (cg_iter < max_cg_iter)
200 {
201 cg_iter++;
202
203 fun_obj->Hv(d, Hd);
204 dHd = ddot_(&n, d, &inc, Hd, &inc);
205 // avoid 0/0 in getting alpha
206 if (dHd <= 1.0e-16)
207 break;
208
209 alpha = zTr/dHd;
210 daxpy_(&n, &alpha, d, &inc, s, &inc);
211 alpha = -alpha;
212 daxpy_(&n, &alpha, Hd, &inc, r, &inc);
213
214 // Using quadratic approximation as CG stopping criterion
215 newQ = -0.5*(ddot_(&n, s, &inc, r, &inc) - ddot_(&n, s, &inc, g, &inc));
216 Qdiff = newQ - Q;
217 if (newQ <= 0 && Qdiff <= 0)
218 {
219 if (cg_iter * Qdiff >= cgtol * newQ)
220 break;
221 }
222 else
223 {
224 info("WARNING: quadratic approximation > 0 or increasing in CG\n");
225 break;
226 }
227 Q = newQ;
228
229 for (i=0; i<n; i++)
230 z[i] = r[i] / M[i];
231 znewTrnew = ddot_(&n, z, &inc, r, &inc);
232 beta = znewTrnew/zTr;
233 dscal_(&n, &beta, d, &inc);
234 daxpy_(&n, &one, z, &inc, d, &inc);
235 zTr = znewTrnew;
236 }
237
238 if (cg_iter == max_cg_iter)
239 info("WARNING: reaching maximal number of CG steps\n");
240
241 delete[] d;
242 delete[] Hd;
243 delete[] z;
244
245 return cg_iter;
246 }
247
set_print_string(void (* print_string)(const char * buf))248 void NEWTON::set_print_string(void (*print_string) (const char *buf))
249 {
250 newton_print_string = print_string;
251 }
252