source: orange/source/orange/liblinear/tron.cpp @ 8978:f0ba3673b0a7

Revision 8978:f0ba3673b0a7, 5.1 KB checked in by ales_erjavec <ales.erjavec@…>, 3 years ago (diff)

Refactored the code for LibSVM and LIBLINEAR. The sources are unmodified in liblinear and libsvm directories. Can be excluded from the build and instead linked to a system libraries instead.
The most significant change is in the way custom kernels are handled for SVMLearner - it now uses 'PRECOMPUTED' kernel functionality from LibSVM (as a bonus the results seem to be better then before (was there a bug in the previous code?), also faster)

Line 
1#include <math.h>
2#include <stdio.h>
3#include <string.h>
4#include <stdarg.h>
5#include "tron.h"
6
7#ifndef min
8template <class T> static inline T min(T x,T y) { return (x<y)?x:y; }
9#endif
10
11#ifndef max
12template <class T> static inline T max(T x,T y) { return (x>y)?x:y; }
13#endif
14
15#ifdef __cplusplus
16extern "C" {
17#endif
18
19extern double dnrm2_(int *, double *, int *);
20extern double ddot_(int *, double *, int *, double *, int *);
21extern int daxpy_(int *, double *, double *, int *, double *, int *);
22extern int dscal_(int *, double *, double *, int *);
23
24#ifdef __cplusplus
25}
26#endif
27
28static void default_print(const char *buf)
29{
30    fputs(buf,stdout);
31    fflush(stdout);
32}
33
34void TRON::info(const char *fmt,...)
35{
36    char buf[BUFSIZ];
37    va_list ap;
38    va_start(ap,fmt);
39    vsprintf(buf,fmt,ap);
40    va_end(ap);
41    (*tron_print_string)(buf);
42}
43
44TRON::TRON(const function *fun_obj, double eps, int max_iter)
45{
46    this->fun_obj=const_cast<function *>(fun_obj);
47    this->eps=eps;
48    this->max_iter=max_iter;
49    tron_print_string = default_print;
50}
51
52TRON::~TRON()
53{
54}
55
56void TRON::tron(double *w)
57{
58    // Parameters for updating the iterates.
59    double eta0 = 1e-4, eta1 = 0.25, eta2 = 0.75;
60
61    // Parameters for updating the trust region size delta.
62    double sigma1 = 0.25, sigma2 = 0.5, sigma3 = 4;
63
64    int n = fun_obj->get_nr_variable();
65    int i, cg_iter;
66    double delta, snorm, one=1.0;
67    double alpha, f, fnew, prered, actred, gs;
68    int search = 1, iter = 1, inc = 1;
69    double *s = new double[n];
70    double *r = new double[n];
71    double *w_new = new double[n];
72    double *g = new double[n];
73
74    for (i=0; i<n; i++)
75        w[i] = 0;
76
77        f = fun_obj->fun(w);
78    fun_obj->grad(w, g);
79    delta = dnrm2_(&n, g, &inc);
80    double gnorm1 = delta;
81    double gnorm = gnorm1;
82
83    if (gnorm <= eps*gnorm1)
84        search = 0;
85
86    iter = 1;
87
88    while (iter <= max_iter && search)
89    {
90        cg_iter = trcg(delta, g, s, r);
91
92        memcpy(w_new, w, sizeof(double)*n);
93        daxpy_(&n, &one, s, &inc, w_new, &inc);
94
95        gs = ddot_(&n, g, &inc, s, &inc);
96        prered = -0.5*(gs-ddot_(&n, s, &inc, r, &inc));
97                fnew = fun_obj->fun(w_new);
98
99        // Compute the actual reduction.
100            actred = f - fnew;
101
102        // On the first iteration, adjust the initial step bound.
103        snorm = dnrm2_(&n, s, &inc);
104        if (iter == 1)
105            delta = min(delta, snorm);
106
107        // Compute prediction alpha*snorm of the step.
108        if (fnew - f - gs <= 0)
109            alpha = sigma3;
110        else
111            alpha = max(sigma1, -0.5*(gs/(fnew - f - gs)));
112
113        // Update the trust region bound according to the ratio of actual to predicted reduction.
114        if (actred < eta0*prered)
115            delta = min(max(alpha, sigma1)*snorm, sigma2*delta);
116        else if (actred < eta1*prered)
117            delta = max(sigma1*delta, min(alpha*snorm, sigma2*delta));
118        else if (actred < eta2*prered)
119            delta = max(sigma1*delta, min(alpha*snorm, sigma3*delta));
120        else
121            delta = max(delta, min(alpha*snorm, sigma3*delta));
122
123        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);
124
125        if (actred > eta0*prered)
126        {
127            iter++;
128            memcpy(w, w_new, sizeof(double)*n);
129            f = fnew;
130                fun_obj->grad(w, g);
131
132            gnorm = dnrm2_(&n, g, &inc);
133            if (gnorm <= eps*gnorm1)
134                break;
135        }
136        if (f < -1.0e+32)
137        {
138            info("warning: f < -1.0e+32\n");
139            break;
140        }
141        if (fabs(actred) <= 0 && prered <= 0)
142        {
143            info("warning: actred and prered <= 0\n");
144            break;
145        }
146        if (fabs(actred) <= 1.0e-12*fabs(f) &&
147            fabs(prered) <= 1.0e-12*fabs(f))
148        {
149            info("warning: actred and prered too small\n");
150            break;
151        }
152    }
153
154    delete[] g;
155    delete[] r;
156    delete[] w_new;
157    delete[] s;
158}
159
160int TRON::trcg(double delta, double *g, double *s, double *r)
161{
162    int i, inc = 1;
163    int n = fun_obj->get_nr_variable();
164    double one = 1;
165    double *d = new double[n];
166    double *Hd = new double[n];
167    double rTr, rnewTrnew, alpha, beta, cgtol;
168
169    for (i=0; i<n; i++)
170    {
171        s[i] = 0;
172        r[i] = -g[i];
173        d[i] = r[i];
174    }
175    cgtol = 0.1*dnrm2_(&n, g, &inc);
176
177    int cg_iter = 0;
178    rTr = ddot_(&n, r, &inc, r, &inc);
179    while (1)
180    {
181        if (dnrm2_(&n, r, &inc) <= cgtol)
182            break;
183        cg_iter++;
184        fun_obj->Hv(d, Hd);
185
186        alpha = rTr/ddot_(&n, d, &inc, Hd, &inc);
187        daxpy_(&n, &alpha, d, &inc, s, &inc);
188        if (dnrm2_(&n, s, &inc) > delta)
189        {
190            info("cg reaches trust region boundary\n");
191            alpha = -alpha;
192            daxpy_(&n, &alpha, d, &inc, s, &inc);
193
194            double std = ddot_(&n, s, &inc, d, &inc);
195            double sts = ddot_(&n, s, &inc, s, &inc);
196            double dtd = ddot_(&n, d, &inc, d, &inc);
197            double dsq = delta*delta;
198            double rad = sqrt(std*std + dtd*(dsq-sts));
199            if (std >= 0)
200                alpha = (dsq - sts)/(std + rad);
201            else
202                alpha = (rad - std)/dtd;
203            daxpy_(&n, &alpha, d, &inc, s, &inc);
204            alpha = -alpha;
205            daxpy_(&n, &alpha, Hd, &inc, r, &inc);
206            break;
207        }
208        alpha = -alpha;
209        daxpy_(&n, &alpha, Hd, &inc, r, &inc);
210        rnewTrnew = ddot_(&n, r, &inc, r, &inc);
211        beta = rnewTrnew/rTr;
212        dscal_(&n, &beta, d, &inc);
213        daxpy_(&n, &one, r, &inc, d, &inc);
214        rTr = rnewTrnew;
215    }
216
217    delete[] d;
218    delete[] Hd;
219
220    return(cg_iter);
221}
222
223double TRON::norm_inf(int n, double *x)
224{
225    double dmax = fabs(x[0]);
226    for (int i=1; i<n; i++)
227        if (fabs(x[i]) >= dmax)
228            dmax = fabs(x[i]);
229    return(dmax);
230}
231
232void TRON::set_print_string(void (*print_string) (const char *buf))
233{
234    tron_print_string = print_string;
235}
Note: See TracBrowser for help on using the repository browser.