view release on metacpan or  search on metacpan

src/liblinear/newton.cpp  view on Meta::CPAN

	const double alpha_pcg = 0.01;
	double *M = new double[n];

	// calculate gradient norm at w=0 for stopping condition.
	double *w0 = new double[n];
	for (i=0; i<n; i++)
		w0[i] = 0;
	fun_obj->fun(w0);
	fun_obj->grad(w0, g);
	double gnorm0 = dnrm2_(&n, g, &inc);
	delete [] w0;

	f = fun_obj->fun(w);
	fun_obj->grad(w, g);
	double gnorm = dnrm2_(&n, g, &inc);
	info("init f %5.3e |g| %5.3e\n", f, gnorm);

	if (gnorm <= eps*gnorm0)
		search = 0;

	while (iter <= max_iter && search)
	{
		fun_obj->get_diag_preconditioner(M);
		for(i=0; i<n; i++)
			M[i] = (1-alpha_pcg) + alpha_pcg*M[i];
		cg_iter = pcg(g, M, s, r);

		fold = f;
		step_size = fun_obj->linesearch_and_update(w, s, &f, g, init_step_size);

		if (step_size == 0)
		{
			info("WARNING: line search fails\n");
			break;
		}

		fun_obj->grad(w, g);
		gnorm = dnrm2_(&n, g, &inc);

		info("iter %2d f %5.3e |g| %5.3e CG %3d step_size %4.2e \n", iter, f, gnorm, cg_iter, step_size);
		
		if (gnorm <= eps*gnorm0)
			break;
		if (f < -1.0e+32)
		{
			info("WARNING: f < -1.0e+32\n");
			break;
		}
		actred = fold - f;
		if (fabs(actred) <= 1.0e-12*fabs(f))
		{
			info("WARNING: actred too small\n");
			break;
		}

		iter++;
	}

	if(iter >= max_iter)
		info("\nWARNING: reaching max number of Newton iterations\n");

	delete[] g;
	delete[] r;
	delete[] s;
	delete[] M;
}

int NEWTON::pcg(double *g, double *M, double *s, double *r)
{
	int i, inc = 1;
	int n = fun_obj->get_nr_variable();
	double one = 1;
	double *d = new double[n];
	double *Hd = new double[n];
	double zTr, znewTrnew, alpha, beta, cgtol, dHd;
	double *z = new double[n];
	double Q = 0, newQ, Qdiff;

	for (i=0; i<n; i++)
	{
		s[i] = 0;
		r[i] = -g[i];
		z[i] = r[i] / M[i];
		d[i] = z[i];
	}

	zTr = ddot_(&n, z, &inc, r, &inc);
	double gMinv_norm = sqrt(zTr);
	cgtol = min(eps_cg, sqrt(gMinv_norm));
	int cg_iter = 0;
	int max_cg_iter = max(n, 5);

	while (cg_iter < max_cg_iter)
	{
		cg_iter++;

		fun_obj->Hv(d, Hd);
		dHd = ddot_(&n, d, &inc, Hd, &inc);
		// avoid 0/0 in getting alpha
		if (dHd <= 1.0e-16)
			break;
		
		alpha = zTr/dHd;
		daxpy_(&n, &alpha, d, &inc, s, &inc);
		alpha = -alpha;
		daxpy_(&n, &alpha, Hd, &inc, r, &inc);

		// Using quadratic approximation as CG stopping criterion
		newQ = -0.5*(ddot_(&n, s, &inc, r, &inc) - ddot_(&n, s, &inc, g, &inc));
		Qdiff = newQ - Q;
		if (newQ <= 0 && Qdiff <= 0)
		{
			if (cg_iter * Qdiff >= cgtol * newQ)
				break;
		}
		else
		{
			info("WARNING: quadratic approximation > 0 or increasing in CG\n");
			break;
		}