1 /*
2 This code was extracted from liblinear-2.2.1 in Feb 2019 and
3 modified for the use with Octave and Matlab
4
5 Copyright (c) 2007-2019 The LIBLINEAR Project.
6 All rights reserved.
7
8 Redistribution and use in source and binary forms, with or without
9 modification, are permitted provided that the following conditions
10 are met:
11
12 1. Redistributions of source code must retain the above copyright
13 notice, this list of conditions and the following disclaimer.
14
15 2. Redistributions in binary form must reproduce the above copyright
16 notice, this list of conditions and the following disclaimer in the
17 documentation and/or other materials provided with the distribution.
18
19 3. Neither name of copyright holders nor the names of its contributors
20 may be used to endorse or promote products derived from this software
21 without specific prior written permission.
22
23
24 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
25 ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
26 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
27 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR
28 CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
29 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
30 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
31 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
32 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
33 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
34 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
35
36 */
37
38 #include <math.h>
39 #include <stdio.h>
40 #include <string.h>
41 #include <stdarg.h>
42 #include "tron.h"
43
44 #ifndef min
min(T x,T y)45 template <class T> static inline T min(T x,T y) { return (x<y)?x:y; }
46 #endif
47
48 #ifndef max
max(T x,T y)49 template <class T> static inline T max(T x,T y) { return (x>y)?x:y; }
50 #endif
51
52 #ifdef __cplusplus
53 extern "C" {
54 #endif
55
56 extern double dnrm2_(int *, double *, int *);
57 extern double ddot_(int *, double *, int *, double *, int *);
58 extern int daxpy_(int *, double *, double *, int *, double *, int *);
59 extern int dscal_(int *, double *, double *, int *);
60
61 #ifdef __cplusplus
62 }
63 #endif
64
default_print(const char * buf)65 static void default_print(const char *buf)
66 {
67 fputs(buf,stdout);
68 fflush(stdout);
69 }
70
uTMv(int n,double * u,double * M,double * v)71 static double uTMv(int n, double *u, double *M, double *v)
72 {
73 const int m = n-4;
74 double res = 0;
75 int i;
76 for (i=0; i<m; i+=5)
77 res += u[i]*M[i]*v[i]+u[i+1]*M[i+1]*v[i+1]+u[i+2]*M[i+2]*v[i+2]+
78 u[i+3]*M[i+3]*v[i+3]+u[i+4]*M[i+4]*v[i+4];
79 for (; i<n; i++)
80 res += u[i]*M[i]*v[i];
81 return res;
82 }
83
info(const char * fmt,...)84 void TRON::info(const char *fmt,...)
85 {
86 char buf[BUFSIZ];
87 va_list ap;
88 va_start(ap,fmt);
89 vsprintf(buf,fmt,ap);
90 va_end(ap);
91 (*tron_print_string)(buf);
92 }
93
TRON(const function * fun_obj,double eps,double eps_cg,int max_iter)94 TRON::TRON(const function *fun_obj, double eps, double eps_cg, int max_iter)
95 {
96 this->fun_obj=const_cast<function *>(fun_obj);
97 this->eps=eps;
98 this->eps_cg=eps_cg;
99 this->max_iter=max_iter;
100 tron_print_string = default_print;
101 }
102
~TRON()103 TRON::~TRON()
104 {
105 }
106
tron(double * w)107 void TRON::tron(double *w)
108 {
109 // Parameters for updating the iterates.
110 double eta0 = 1e-4, eta1 = 0.25, eta2 = 0.75;
111
112 // Parameters for updating the trust region size delta.
113 double sigma1 = 0.25, sigma2 = 0.5, sigma3 = 4;
114
115 int n = fun_obj->get_nr_variable();
116 int i, cg_iter;
117 double delta=0, sMnorm, one=1.0;
118 double alpha, f, fnew, prered, actred, gs;
119 int search = 1, iter = 1, inc = 1;
120 double *s = new double[n];
121 double *r = new double[n];
122 double *g = new double[n];
123
124 const double alpha_pcg = 0.01;
125 double *M = new double[n];
126
127 // calculate gradient norm at w=0 for stopping condition.
128 double *w0 = new double[n];
129 for (i=0; i<n; i++)
130 w0[i] = 0;
131 fun_obj->fun(w0);
132 fun_obj->grad(w0, g);
133 double gnorm0 = dnrm2_(&n, g, &inc);
134 delete [] w0;
135
136 f = fun_obj->fun(w);
137 fun_obj->grad(w, g);
138 double gnorm = dnrm2_(&n, g, &inc);
139
140 if (gnorm <= eps*gnorm0)
141 search = 0;
142
143 fun_obj->get_diag_preconditioner(M);
144 for(i=0; i<n; i++)
145 M[i] = (1-alpha_pcg) + alpha_pcg*M[i];
146 delta = sqrt(uTMv(n, g, M, g));
147
148 double *w_new = new double[n];
149 bool reach_boundary;
150 bool delta_adjusted = false;
151 while (iter <= max_iter && search)
152 {
153 cg_iter = trpcg(delta, g, M, s, r, &reach_boundary);
154
155 memcpy(w_new, w, sizeof(double)*n);
156 daxpy_(&n, &one, s, &inc, w_new, &inc);
157
158 gs = ddot_(&n, g, &inc, s, &inc);
159 prered = -0.5*(gs-ddot_(&n, s, &inc, r, &inc));
160 fnew = fun_obj->fun(w_new);
161
162 // Compute the actual reduction.
163 actred = f - fnew;
164
165 // On the first iteration, adjust the initial step bound.
166 sMnorm = sqrt(uTMv(n, s, M, s));
167 if (iter == 1 && !delta_adjusted)
168 {
169 delta = min(delta, sMnorm);
170 delta_adjusted = true;
171 }
172
173 // Compute prediction alpha*sMnorm of the step.
174 if (fnew - f - gs <= 0)
175 alpha = sigma3;
176 else
177 alpha = max(sigma1, -0.5*(gs/(fnew - f - gs)));
178
179 // Update the trust region bound according to the ratio of actual to predicted reduction.
180 if (actred < eta0*prered)
181 delta = min(alpha*sMnorm, sigma2*delta);
182 else if (actred < eta1*prered)
183 delta = max(sigma1*delta, min(alpha*sMnorm, sigma2*delta));
184 else if (actred < eta2*prered)
185 delta = max(sigma1*delta, min(alpha*sMnorm, sigma3*delta));
186 else
187 {
188 if (reach_boundary)
189 delta = sigma3*delta;
190 else
191 delta = max(delta, min(alpha*sMnorm, sigma3*delta));
192 }
193
194 info("iter %2d act %5.3e pre %5.3e delta %5.3e f %5.3e |g| %5.3e CG %3d\n", iter, actred, prered, delta, f, gnorm, cg_iter);
195
196 if (actred > eta0*prered)
197 {
198 iter++;
199 memcpy(w, w_new, sizeof(double)*n);
200 f = fnew;
201 fun_obj->grad(w, g);
202 fun_obj->get_diag_preconditioner(M);
203 for(i=0; i<n; i++)
204 M[i] = (1-alpha_pcg) + alpha_pcg*M[i];
205
206 gnorm = dnrm2_(&n, g, &inc);
207 if (gnorm <= eps*gnorm0)
208 break;
209 }
210 if (f < -1.0e+32)
211 {
212 info("WARNING: f < -1.0e+32\n");
213 break;
214 }
215 if (prered <= 0)
216 {
217 info("WARNING: prered <= 0\n");
218 break;
219 }
220 if (fabs(actred) <= 1.0e-12*fabs(f) &&
221 fabs(prered) <= 1.0e-12*fabs(f))
222 {
223 info("WARNING: actred and prered too small\n");
224 break;
225 }
226 }
227
228 delete[] g;
229 delete[] r;
230 delete[] w_new;
231 delete[] s;
232 delete[] M;
233 }
234
trpcg(double delta,double * g,double * M,double * s,double * r,bool * reach_boundary)235 int TRON::trpcg(double delta, double *g, double *M, double *s, double *r, bool *reach_boundary)
236 {
237 int i, inc = 1;
238 int n = fun_obj->get_nr_variable();
239 double one = 1;
240 double *d = new double[n];
241 double *Hd = new double[n];
242 double zTr, znewTrnew, alpha, beta, cgtol;
243 double *z = new double[n];
244
245 *reach_boundary = false;
246 for (i=0; i<n; i++)
247 {
248 s[i] = 0;
249 r[i] = -g[i];
250 z[i] = r[i] / M[i];
251 d[i] = z[i];
252 }
253
254 zTr = ddot_(&n, z, &inc, r, &inc);
255 cgtol = eps_cg*sqrt(zTr);
256 int cg_iter = 0;
257
258 while (1)
259 {
260 if (sqrt(zTr) <= cgtol)
261 break;
262 cg_iter++;
263 fun_obj->Hv(d, Hd);
264
265 alpha = zTr/ddot_(&n, d, &inc, Hd, &inc);
266 daxpy_(&n, &alpha, d, &inc, s, &inc);
267
268 double sMnorm = sqrt(uTMv(n, s, M, s));
269 if (sMnorm > delta)
270 {
271 info("cg reaches trust region boundary\n");
272 *reach_boundary = true;
273 alpha = -alpha;
274 daxpy_(&n, &alpha, d, &inc, s, &inc);
275
276 double sTMd = uTMv(n, s, M, d);
277 double sTMs = uTMv(n, s, M, s);
278 double dTMd = uTMv(n, d, M, d);
279 double dsq = delta*delta;
280 double rad = sqrt(sTMd*sTMd + dTMd*(dsq-sTMs));
281 if (sTMd >= 0)
282 alpha = (dsq - sTMs)/(sTMd + rad);
283 else
284 alpha = (rad - sTMd)/dTMd;
285 daxpy_(&n, &alpha, d, &inc, s, &inc);
286 alpha = -alpha;
287 daxpy_(&n, &alpha, Hd, &inc, r, &inc);
288 break;
289 }
290 alpha = -alpha;
291 daxpy_(&n, &alpha, Hd, &inc, r, &inc);
292
293 for (i=0; i<n; i++)
294 z[i] = r[i] / M[i];
295 znewTrnew = ddot_(&n, z, &inc, r, &inc);
296 beta = znewTrnew/zTr;
297 dscal_(&n, &beta, d, &inc);
298 daxpy_(&n, &one, z, &inc, d, &inc);
299 zTr = znewTrnew;
300 }
301
302 delete[] d;
303 delete[] Hd;
304 delete[] z;
305
306 return(cg_iter);
307 }
308
norm_inf(int n,double * x)309 double TRON::norm_inf(int n, double *x)
310 {
311 double dmax = fabs(x[0]);
312 for (int i=1; i<n; i++)
313 if (fabs(x[i]) >= dmax)
314 dmax = fabs(x[i]);
315 return(dmax);
316 }
317
set_print_string(void (* print_string)(const char * buf))318 void TRON::set_print_string(void (*print_string) (const char *buf))
319 {
320 tron_print_string = print_string;
321 }
322