pr_loqo.cpp

Go to the documentation of this file.
00001 /*
00002  * This program is free software; you can redistribute it and/or modify
00003  * it under the terms of the GNU General Public License as published by
00004  * the Free Software Foundation; either version 3 of the License, or
00005  * (at your option) any later version.
00006  *
00007  * Purpose:     solves quadratic programming problem for pattern recognition
00008  *              for support vectors
00009  *
00010  * Written (W) 1997-1998 Alex J. Smola
00011  * Written (W) 1999-2009 Soeren Sonnenburg
00012  * Written (W) 1999-2008 Gunnar Raetsch
00013  * Copyright (C) 1997-2009 Fraunhofer Institute FIRST and Max-Planck-Society
00014  */
00015 
00016 #include "lib/common.h"
00017 #include "lib/io.h"
00018 #include "lib/lapack.h"
00019 #include "lib/Mathematics.h"
00020 #include "classifier/svm/pr_loqo.h"
00021 
00022 #include <math.h>
00023 #include <time.h>
00024 #include <stdlib.h>
00025 #include <stdio.h>
00026 #ifdef HAVE_LAPACK
00027 extern "C" {
00028 #ifdef HAVE_ATLAS
00029 #include <clapack.h>
00030 #endif
00031 }
00032 #endif
00033 
00034 #define PREDICTOR 1
00035 #define CORRECTOR 2
00036 
00037 /*****************************************************************
00038   replace this by any other function that will exit gracefully
00039   in a larger system
00040   ***************************************************************/
00041 
00042 void nrerror(char error_text[])
00043 {
00044     SG_SDEBUG("terminating optimizer - %s\n", error_text);
00045  // exit(1); 
00046 }
00047 
00048 /*****************************************************************
00049    taken from numerical recipes and modified to accept pointers
00050    moreover numerical recipes code seems to be buggy (at least the
00051    ones on the web)
00052 
00053    cholesky solver and backsubstitution
00054    leaves upper right triangle intact (rows first order)
00055    ***************************************************************/
00056 
00057 #ifdef HAVE_LAPACK
00058 bool choldc(float64_t* a, int32_t n, float64_t* p)
00059 {
00060     if (n<=0)
00061         return false;
00062 
00063     float64_t* a2=new float64_t[n*n];
00064 
00065     for (int32_t i=0; i<n; i++)
00066     {
00067         for (int32_t j=0; j<n; j++)
00068             a2[n*i+j]=a[n*i+j];
00069     }
00070 
00071     /* int for calling external lib */
00072     int result=clapack_dpotrf(CblasRowMajor, CblasUpper, (int) n, a2, (int) n);
00073 
00074     for (int32_t i=0; i<n; i++)
00075         p[i]=a2[(n+1)*i];
00076 
00077     for (int32_t i=0; i<n; i++)
00078     {
00079         for (int32_t j=i+1; j<n; j++)
00080         {
00081             a[n*j+i]=a2[n*i+j];
00082         }
00083     }
00084 
00085     if (result>0)
00086         SG_SDEBUG("Choldc failed, matrix not positive definite\n");
00087 
00088     delete[] a2;
00089     
00090     return result==0;
00091 }
00092 #else
00093 bool choldc(float64_t a[], int32_t n, float64_t p[])
00094 {
00095     void nrerror(char error_text[]);
00096     int32_t i, j, k;
00097     float64_t sum;
00098 
00099     for (i = 0; i < n; i++)
00100     {
00101         for (j = i; j < n; j++)
00102         {
00103             sum=a[n*i + j];
00104 
00105             for (k=i-1; k>=0; k--)
00106                 sum -= a[n*i + k]*a[n*j + k];
00107 
00108             if (i == j)
00109             {
00110                 if (sum <= 0.0)
00111                 {
00112                     SG_SDEBUG("Choldc failed, matrix not positive definite");
00113                     sum = 0.0;
00114                     return false;
00115                 }
00116 
00117                 p[i]=sqrt(sum);
00118 
00119             } 
00120             else
00121                 a[n*j + i] = sum/p[i];
00122         }
00123     }
00124 
00125     return true;
00126 }
00127 #endif
00128 
00129 void cholsb(
00130     float64_t a[], int32_t n, float64_t p[], float64_t b[], float64_t x[])
00131 {
00132   int32_t i, k;
00133   float64_t sum;
00134 
00135   for (i=0; i<n; i++) {
00136     sum=b[i];
00137     for (k=i-1; k>=0; k--) sum -= a[n*i + k]*x[k];
00138     x[i]=sum/p[i];
00139   }
00140 
00141   for (i=n-1; i>=0; i--) {
00142     sum=x[i];
00143     for (k=i+1; k<n; k++) sum -= a[n*k + i]*x[k];
00144     x[i]=sum/p[i];
00145   }
00146 }
00147 
00148 /*****************************************************************
00149   sometimes we only need the forward or backward pass of the
00150   backsubstitution, hence we provide these two routines separately 
00151   ***************************************************************/
00152 
00153 void chol_forward(
00154     float64_t a[], int32_t n, float64_t p[], float64_t b[], float64_t x[])
00155 {
00156   int32_t i, k;
00157   float64_t sum;
00158 
00159   for (i=0; i<n; i++) {
00160     sum=b[i];
00161     for (k=i-1; k>=0; k--) sum -= a[n*i + k]*x[k];
00162     x[i]=sum/p[i];
00163   }
00164 }
00165 
00166 void chol_backward(
00167     float64_t a[], int32_t n, float64_t p[], float64_t b[], float64_t x[])
00168 {
00169   int32_t i, k;
00170   float64_t sum;
00171 
00172   for (i=n-1; i>=0; i--) {
00173     sum=b[i];
00174     for (k=i+1; k<n; k++) sum -= a[n*k + i]*x[k];
00175     x[i]=sum/p[i];
00176   }
00177 }
00178 
00179 /*****************************************************************
00180   solves the system | -H_x A' | |x_x| = |c_x|
00181                     |  A   H_y| |x_y|   |c_y|
00182 
00183   with H_x (and H_y) positive (semidefinite) matrices
00184   and n, m the respective sizes of H_x and H_y
00185 
00186   for variables see pg. 48 of notebook or do the calculations on a
00187   sheet of paper again
00188 
00189   predictor solves the whole thing, corrector assues that H_x didn't
00190   change and relies on the results of the predictor. therefore do
00191   _not_ modify workspace
00192 
00193   if you want to speed tune anything in the code here's the right
00194   place to do so: about 95% of the time is being spent in
00195   here. something like an iterative refinement would be nice,
00196   especially when switching from float64_t to single precision. if you
00197   have a fast parallel cholesky use it instead of the numrec
00198   implementations.
00199 
00200   side effects: changes H_y (but this is just the unit matrix or zero anyway
00201   in our case)
00202   ***************************************************************/
00203 
00204 bool solve_reduced(
00205     int32_t n, int32_t m, float64_t h_x[], float64_t h_y[], float64_t a[],
00206     float64_t x_x[], float64_t x_y[], float64_t c_x[], float64_t c_y[],
00207     float64_t workspace[], int32_t step)
00208 {
00209   int32_t i,j,k;
00210 
00211   float64_t *p_x;
00212   float64_t *p_y;
00213   float64_t *t_a;
00214   float64_t *t_c;
00215   float64_t *t_y;
00216 
00217   p_x = workspace;      /* together n + m + n*m + n + m = n*(m+2)+2*m */
00218   p_y = p_x + n;
00219   t_a = p_y + m;
00220   t_c = t_a + n*m;
00221   t_y = t_c + n;
00222 
00223   if (step == PREDICTOR) {
00224     if (!choldc(h_x, n, p_x))   /* do cholesky decomposition */
00225         return false;
00226 
00227     for (i=0; i<m; i++)         /* forward pass for A' */
00228       chol_forward(h_x, n, p_x, a+i*n, t_a+i*n);
00229                 
00230     for (i=0; i<m; i++)         /* compute (h_y + a h_x^-1A') */
00231       for (j=i; j<m; j++)
00232     for (k=0; k<n; k++) 
00233       h_y[m*i + j] += t_a[n*j + k] * t_a[n*i + k];
00234                 
00235     choldc(h_y, m, p_y);    /* and cholesky decomposition */
00236   }
00237   
00238   chol_forward(h_x, n, p_x, c_x, t_c);
00239                 /* forward pass for c */
00240 
00241   for (i=0; i<m; i++) {     /* and solve for x_y */
00242     t_y[i] = c_y[i];
00243     for (j=0; j<n; j++)
00244       t_y[i] += t_a[i*n + j] * t_c[j];
00245   }
00246 
00247   cholsb(h_y, m, p_y, t_y, x_y);
00248 
00249   for (i=0; i<n; i++) {     /* finally solve for x_x */
00250     t_c[i] = -t_c[i];
00251     for (j=0; j<m; j++)
00252       t_c[i] += t_a[j*n + i] * x_y[j];
00253   }
00254 
00255   chol_backward(h_x, n, p_x, t_c, x_x);
00256   return true;
00257 }
00258 
00259 /*****************************************************************
00260   matrix vector multiplication (symmetric matrix but only one triangle
00261   given). computes m*x = y
00262   no need to tune it as it's only of O(n^2) but cholesky is of
00263   O(n^3). so don't waste your time _here_ although it isn't very
00264   elegant. 
00265   ***************************************************************/
00266 
00267 void matrix_vector(int32_t n, float64_t m[], float64_t x[], float64_t y[])
00268 {
00269   int32_t i, j;
00270 
00271   for (i=0; i<n; i++) {
00272     y[i] = m[(n+1) * i] * x[i];
00273 
00274     for (j=0; j<i; j++)
00275       y[i] += m[i + n*j] * x[j];
00276 
00277     for (j=i+1; j<n; j++) 
00278       y[i] += m[n*i + j] * x[j]; 
00279   }
00280 }
00281 
00282 /*****************************************************************
00283   call only this routine; this is the only one you're interested in
00284   for doing quadratical optimization
00285 
00286   the restart feature exists but it may not be of much use due to the
00287   fact that an initial setting, although close but not very close the
00288   the actual solution will result in very good starting diagnostics
00289   (primal and dual feasibility and small infeasibility gap) but incur
00290   later stalling of the optimizer afterwards as we have to enforce
00291   positivity of the slacks.
00292   ***************************************************************/
00293 
00294 int32_t pr_loqo(
00295     int32_t n, int32_t m, float64_t c[], float64_t h_x[], float64_t a[],
00296     float64_t b[], float64_t l[], float64_t u[], float64_t primal[],
00297     float64_t dual[], int32_t verb, float64_t sigfig_max, int32_t counter_max,
00298     float64_t margin, float64_t bound, int32_t restart)
00299 {
00300   /* the knobs to be tuned ... */
00301   /* float64_t margin = -0.95;     we will go up to 95% of the
00302                    distance between old variables and zero */
00303   /* float64_t bound = 10;         preset value for the start. small
00304                    values give good initial
00305                    feasibility but may result in slow
00306                    convergence afterwards: we're too
00307                    close to zero */
00308   /* to be allocated */
00309   float64_t *workspace;
00310   float64_t *diag_h_x;
00311   float64_t *h_y;
00312   float64_t *c_x;
00313   float64_t *c_y;
00314   float64_t *h_dot_x;
00315   float64_t *rho;
00316   float64_t *nu;
00317   float64_t *tau;
00318   float64_t *sigma;
00319   float64_t *gamma_z;
00320   float64_t *gamma_s;  
00321 
00322   float64_t *hat_nu;
00323   float64_t *hat_tau;
00324 
00325   float64_t *delta_x;
00326   float64_t *delta_y;
00327   float64_t *delta_s;
00328   float64_t *delta_z;
00329   float64_t *delta_g;
00330   float64_t *delta_t;
00331 
00332   float64_t *d;
00333 
00334   /* from the header - pointers into primal and dual */
00335   float64_t *x;
00336   float64_t *y;
00337   float64_t *g;
00338   float64_t *z;
00339   float64_t *s;
00340   float64_t *t;  
00341 
00342   /* auxiliary variables */
00343   float64_t b_plus_1;
00344   float64_t c_plus_1;
00345 
00346   float64_t x_h_x;
00347   float64_t primal_inf;
00348   float64_t dual_inf;
00349 
00350   float64_t sigfig;
00351   float64_t primal_obj, dual_obj;
00352   float64_t mu;
00353   float64_t alfa=-1;
00354   int32_t counter = 0;
00355 
00356   int32_t status = STILL_RUNNING;
00357 
00358   int32_t i,j;
00359 
00360   /* memory allocation */
00361   workspace = (float64_t*) malloc((n*(m+2)+2*m)*sizeof(float64_t));
00362   diag_h_x  = (float64_t*) malloc(n*sizeof(float64_t));
00363   h_y       = (float64_t*) malloc(m*m*sizeof(float64_t));
00364   c_x       = (float64_t*) malloc(n*sizeof(float64_t));
00365   c_y       = (float64_t*) malloc(m*sizeof(float64_t));
00366   h_dot_x   = (float64_t*) malloc(n*sizeof(float64_t));
00367 
00368   rho       = (float64_t*) malloc(m*sizeof(float64_t));
00369   nu        = (float64_t*) malloc(n*sizeof(float64_t));
00370   tau       = (float64_t*) malloc(n*sizeof(float64_t));
00371   sigma     = (float64_t*) malloc(n*sizeof(float64_t));
00372 
00373   gamma_z   = (float64_t*) malloc(n*sizeof(float64_t));
00374   gamma_s   = (float64_t*) malloc(n*sizeof(float64_t));
00375 
00376   hat_nu    = (float64_t*) malloc(n*sizeof(float64_t));
00377   hat_tau   = (float64_t*) malloc(n*sizeof(float64_t));
00378 
00379   delta_x   = (float64_t*) malloc(n*sizeof(float64_t));
00380   delta_y   = (float64_t*) malloc(m*sizeof(float64_t));
00381   delta_s   = (float64_t*) malloc(n*sizeof(float64_t));
00382   delta_z   = (float64_t*) malloc(n*sizeof(float64_t));
00383   delta_g   = (float64_t*) malloc(n*sizeof(float64_t));
00384   delta_t   = (float64_t*) malloc(n*sizeof(float64_t));
00385 
00386   d         = (float64_t*) malloc(n*sizeof(float64_t));
00387 
00388   /* pointers into the external variables */
00389   x = primal;           /* n */
00390   g = x + n;            /* n */
00391   t = g + n;            /* n */
00392 
00393   y = dual;         /* m */
00394   z = y + m;            /* n */
00395   s = z + n;            /* n */
00396 
00397   /* initial settings */
00398   b_plus_1 = 1;
00399   c_plus_1 = 0;
00400   for (i=0; i<n; i++) c_plus_1 += c[i];
00401 
00402   /* get diagonal terms */
00403   for (i=0; i<n; i++) diag_h_x[i] = h_x[(n+1)*i]; 
00404   
00405   /* starting point */
00406   if (restart == 1) {
00407                 /* x, y already preset */
00408     for (i=0; i<n; i++) {   /* compute g, t for primal feasibility */
00409       g[i] = CMath::max(CMath::abs(x[i] - l[i]), bound);
00410       t[i] = CMath::max(CMath::abs(u[i] - x[i]), bound); 
00411     }
00412 
00413     matrix_vector(n, h_x, x, h_dot_x); /* h_dot_x = h_x * x */
00414 
00415     for (i=0; i<n; i++) {   /* sigma is a dummy variable to calculate z, s */
00416       sigma[i] = c[i] + h_dot_x[i];
00417       for (j=0; j<m; j++)
00418     sigma[i] -= a[n*j + i] * y[j];
00419 
00420       if (sigma[i] > 0) {
00421     s[i] = bound;
00422     z[i] = sigma[i] + bound;
00423       }
00424       else {
00425     s[i] = bound - sigma[i];
00426     z[i] = bound;
00427       }
00428     }
00429   }
00430   else {            /* use default start settings */
00431     for (i=0; i<m; i++)
00432       for (j=i; j<m; j++)
00433     h_y[i*m + j] = (i==j) ? 1 : 0;
00434     
00435     for (i=0; i<n; i++) {
00436       c_x[i] = c[i];
00437       h_x[(n+1)*i] += 1;
00438     }
00439     
00440     for (i=0; i<m; i++)
00441       c_y[i] = b[i];
00442     
00443     /* and solve the system [-H_x A'; A H_y] [x, y] = [c_x; c_y] */
00444     solve_reduced(n, m, h_x, h_y, a, x, y, c_x, c_y, workspace,
00445           PREDICTOR);
00446     
00447     /* initialize the other variables */
00448     for (i=0; i<n; i++) {
00449       g[i] = CMath::max(CMath::abs(x[i] - l[i]), bound);
00450       z[i] = CMath::max(CMath::abs(x[i]), bound);
00451       t[i] = CMath::max(CMath::abs(u[i] - x[i]), bound); 
00452       s[i] = CMath::max(CMath::abs(x[i]), bound); 
00453     }
00454   }
00455 
00456   for (i=0, mu=0; i<n; i++)
00457     mu += z[i] * g[i] + s[i] * t[i];
00458   mu = mu / (2*n);
00459 
00460   /* the main loop */
00461   if (verb >= STATUS) {
00462       SG_SDEBUG("counter | pri_inf  | dual_inf  | pri_obj   | dual_obj  | ");
00463       SG_SDEBUG("sigfig | alpha  | nu \n");
00464       SG_SDEBUG("-------------------------------------------------------");
00465       SG_SDEBUG("---------------------------\n");
00466   }
00467   
00468   while (status == STILL_RUNNING) {
00469     /* predictor */
00470     
00471     /* put back original diagonal values */
00472     for (i=0; i<n; i++) 
00473       h_x[(n+1) * i] = diag_h_x[i];
00474 
00475     matrix_vector(n, h_x, x, h_dot_x); /* compute h_dot_x = h_x * x */
00476 
00477     for (i=0; i<m; i++) {
00478       rho[i] = b[i];
00479       for (j=0; j<n; j++)
00480     rho[i] -= a[n*i + j] * x[j];
00481     }
00482     
00483     for (i=0; i<n; i++) {
00484       nu[i] = l[i] - x[i] + g[i];
00485       tau[i] = u[i] - x[i] - t[i];
00486 
00487       sigma[i] = c[i] - z[i] + s[i] + h_dot_x[i];
00488       for (j=0; j<m; j++)
00489     sigma[i] -= a[n*j + i] * y[j];
00490 
00491       gamma_z[i] = - z[i];
00492       gamma_s[i] = - s[i];
00493     }
00494 
00495     /* instrumentation */
00496     x_h_x = 0;
00497     primal_inf = 0;
00498     dual_inf = 0;
00499     
00500     for (i=0; i<n; i++) {
00501       x_h_x += h_dot_x[i] * x[i];
00502       primal_inf += CMath::sq(tau[i]);
00503       primal_inf += CMath::sq(nu[i]);
00504       dual_inf += CMath::sq(sigma[i]);
00505     }
00506     for (i=0; i<m; i++) 
00507       primal_inf += CMath::sq(rho[i]);
00508     primal_inf = sqrt(primal_inf)/b_plus_1;
00509     dual_inf = sqrt(dual_inf)/c_plus_1;
00510 
00511     primal_obj = 0.5 * x_h_x;
00512     dual_obj = -0.5 * x_h_x;
00513     for (i=0; i<n; i++) {
00514       primal_obj += c[i] * x[i];
00515       dual_obj += l[i] * z[i] - u[i] * s[i];
00516     }
00517     for (i=0; i<m; i++)
00518       dual_obj += b[i] * y[i];
00519 
00520     sigfig = log10(CMath::abs(primal_obj) + 1) -
00521              log10(CMath::abs(primal_obj - dual_obj));
00522     sigfig = CMath::max(sigfig, 0.0);
00523            
00524     /* the diagnostics - after we computed our results we will
00525        analyze them */
00526 
00527     if (counter > counter_max) status = ITERATION_LIMIT;
00528     if (sigfig  > sigfig_max)  status = OPTIMAL_SOLUTION;
00529     if (primal_inf > 10e100)   status = PRIMAL_INFEASIBLE;
00530     if (dual_inf > 10e100)     status = DUAL_INFEASIBLE;
00531     if ((primal_inf > 10e100) & (dual_inf > 10e100)) status = PRIMAL_AND_DUAL_INFEASIBLE;
00532     if (CMath::abs(primal_obj) > 10e100) status = PRIMAL_UNBOUNDED;
00533     if (CMath::abs(dual_obj) > 10e100) status = DUAL_UNBOUNDED;
00534 
00535     /* write some nice routine to enforce the time limit if you
00536        _really_ want, however it's quite useless as you can compute
00537        the time from the maximum number of iterations as every
00538        iteration costs one cholesky decomposition plus a couple of 
00539        backsubstitutions */
00540 
00541     /* generate report */
00542     if ((verb >= FLOOD) | ((verb == STATUS) & (status != 0)))
00543      SG_SDEBUG("%7i | %.2e | %.2e | % .2e | % .2e | %6.3f | %.4f | %.2e\n",
00544          counter, primal_inf, dual_inf, primal_obj, dual_obj,
00545          sigfig, alfa, mu);
00546 
00547     counter++;
00548 
00549     if (status == 0) {      /* we may keep on going, otherwise
00550                    it'll cost one loop extra plus a
00551                    messed up main diagonal of h_x */
00552       /* intermediate variables (the ones with hat) */
00553       for (i=0; i<n; i++) {
00554     hat_nu[i] = nu[i] + g[i] * gamma_z[i] / z[i];
00555     hat_tau[i] = tau[i] - t[i] * gamma_s[i] / s[i];
00556     /* diagonal terms */
00557     d[i] = z[i] / g[i] + s[i] / t[i];
00558       }
00559       
00560       /* initialization before the cholesky solver */
00561       for (i=0; i<n; i++) {
00562     h_x[(n+1)*i] = diag_h_x[i] + d[i];
00563     c_x[i] = sigma[i] - z[i] * hat_nu[i] / g[i] - 
00564       s[i] * hat_tau[i] / t[i];
00565       }
00566       for (i=0; i<m; i++) {
00567     c_y[i] = rho[i];
00568     for (j=i; j<m; j++) 
00569       h_y[m*i + j] = 0;
00570       }
00571       
00572       /* and do it */
00573       if (!solve_reduced(n, m, h_x, h_y, a, delta_x, delta_y, c_x, c_y, workspace, PREDICTOR))
00574       {
00575           status=INCONSISTENT;
00576           goto exit_optimizer;
00577       }
00578       
00579       for (i=0; i<n; i++) {
00580     /* backsubstitution */
00581     delta_s[i] = s[i] * (delta_x[i] - hat_tau[i]) / t[i];
00582     delta_z[i] = z[i] * (hat_nu[i] - delta_x[i]) / g[i];
00583     
00584     delta_g[i] = g[i] * (gamma_z[i] - delta_z[i]) / z[i];
00585     delta_t[i] = t[i] * (gamma_s[i] - delta_s[i]) / s[i];
00586     
00587     /* central path (corrector) */
00588     gamma_z[i] = mu / g[i] - z[i] - delta_z[i] * delta_g[i] / g[i];
00589     gamma_s[i] = mu / t[i] - s[i] - delta_s[i] * delta_t[i] / t[i];
00590     
00591     /* (some more intermediate variables) the hat variables */
00592     hat_nu[i] = nu[i] + g[i] * gamma_z[i] / z[i];
00593     hat_tau[i] = tau[i] - t[i] * gamma_s[i] / s[i];
00594     
00595     /* initialization before the cholesky */
00596     c_x[i] = sigma[i] - z[i] * hat_nu[i] / g[i] - s[i] * hat_tau[i] / t[i];
00597       }
00598       
00599       for (i=0; i<m; i++) { /* comput c_y and rho */
00600     c_y[i] = rho[i];
00601     for (j=i; j<m; j++)
00602       h_y[m*i + j] = 0;
00603       }
00604       
00605       /* and do it */
00606       solve_reduced(n, m, h_x, h_y, a, delta_x, delta_y, c_x, c_y, workspace,
00607             CORRECTOR);
00608       
00609       for (i=0; i<n; i++) {
00610     /* backsubstitution */
00611     delta_s[i] = s[i] * (delta_x[i] - hat_tau[i]) / t[i];
00612     delta_z[i] = z[i] * (hat_nu[i] - delta_x[i]) / g[i];
00613     
00614     delta_g[i] = g[i] * (gamma_z[i] - delta_z[i]) / z[i];
00615     delta_t[i] = t[i] * (gamma_s[i] - delta_s[i]) / s[i];
00616       }
00617       
00618       alfa = -1;
00619       for (i=0; i<n; i++) {
00620     alfa = CMath::min(alfa, delta_g[i]/g[i]);
00621     alfa = CMath::min(alfa, delta_t[i]/t[i]);
00622     alfa = CMath::min(alfa, delta_s[i]/s[i]);
00623     alfa = CMath::min(alfa, delta_z[i]/z[i]);
00624       }
00625       alfa = (margin - 1) / alfa;
00626       
00627       /* compute mu */
00628       for (i=0, mu=0; i<n; i++)
00629     mu += z[i] * g[i] + s[i] * t[i];
00630       mu = mu / (2*n);
00631       mu = mu * CMath::sq((alfa - 1) / (alfa + 10));
00632       
00633       for (i=0; i<n; i++) {
00634     x[i] += alfa * delta_x[i];
00635     g[i] += alfa * delta_g[i];
00636     t[i] += alfa * delta_t[i];
00637     z[i] += alfa * delta_z[i];
00638     s[i] += alfa * delta_s[i];
00639       }
00640       
00641       for (i=0; i<m; i++) 
00642     y[i] += alfa * delta_y[i];
00643     }
00644   }
00645 
00646 exit_optimizer: 
00647   if ((status == 1) && (verb >= STATUS)) {
00648       SG_SDEBUG("----------------------------------------------------------------------------------\n");
00649       SG_SDEBUG("optimization converged\n");
00650   }
00651   
00652   /* free memory */
00653   free(workspace);
00654   free(diag_h_x);
00655   free(h_y);
00656   free(c_x);
00657   free(c_y);
00658   free(h_dot_x);
00659   
00660   free(rho);
00661   free(nu);
00662   free(tau);
00663   free(sigma);
00664   free(gamma_z);
00665   free(gamma_s);
00666   
00667   free(hat_nu);
00668   free(hat_tau);
00669     
00670   free(delta_x);
00671   free(delta_y);
00672   free(delta_s);
00673   free(delta_z);
00674   free(delta_g);
00675   free(delta_t);
00676     
00677   free(d);
00678 
00679   /* and return to sender */
00680   return status;
00681 }

SHOGUN Machine Learning Toolbox - Documentation