Solver.h File Reference

#include "IntPoly.h"
#include "Vector.h"

Include dependency graph for Solver.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions
Vector	L2NormMinimizer (Polynomial q, double delta, int *infoOut=NULL, int maxIter=1000, double *lambda1=NULL)
Vector	L2NormMinimizer (Polynomial q, Vector pointXk, double delta, int *infoOut=NULL, int maxIter=1000, double *lambda1=NULL)
Vector	L2NormMinimizer (Polynomial q, Vector pointXk, double delta, int *infoOut, int maxIter, double *lambda1, Vector minusG, Matrix H)
Vector	LAGMAXModified (Polynomial q, double rho, double &VMAX)
Vector	LAGMAXModified (Polynomial q, Vector pointXk, double rho, double &VMAX)
Vector	LAGMAXModified (Vector G, Matrix H, double rho, double &VMAX)
void	CONDOR (double rhoStart, double rhoEnd, int niter, ObjectiveFunction *of, int nnode=0)

Function Documentation

CONDOR()

void CONDOR	(	double	rhoStart,
		double	rhoEnd,
		int	niter,
		ObjectiveFunction *	of,
		int	nnode = 0
	)

Definition at line 75 of file CNLSolver.cpp.

{
    rhoStart=mmax(rhoStart,rhoEnd);
    int dim=of->dim(), info, k, t, nerror;
    double rho=rhoStart, delta=rhoStart, rhoNew,
           lambda1, normD=rhoEnd+1.0, modelStep, reduction, r, valueOF, valueFk, bound, noise;
    Vector d, tmp;
    bool improvement, forceTRStep=true, evalNeeded;
    
    initConstrainedStep(of);

    // pre-create the MultInd indexes to prevent multi-thread problems:
    cacheMultInd.get(dim,1); cacheMultInd.get(dim,2);    

    parallelInit(nnode, dim, of);

    of->initData();
//    points=getFirstPoints(&ValuesF, &nPtsTotal, rhoStart,of);

    InterPolynomial poly(2, rho, Vector::emptyVector, of->data, of);
    if (poly.NewtonBasis==NULL)
    {
        printf("cannot construct lagrange basis.\n");
        exit(255);
    }

    
    
    //Base=poly.vBase;
    k=poly.kbest;
    valueFk=poly.valuesF[k];

    //fprintf(stderr,"init part 1 finished.\n");

/*    InterPolynomial poly(dim,2,of);
    for (;;)
    {
        // find index k of the best (lowest) value of the function 
        k=findK(ValuesF, nPtsTotal, of, points);
        Base=points[k].clone();
    //    Base=of->xStart.clone();
        valueFk=ValuesF[k];
        // translation:
        t=nPtsTotal; while (t--) points[t]-=Base;

        // exchange index 0 and index k (to be sure best point is inside poly):
        tmp=points[k];
        points[k]=points[0];
        points[0]=tmp;
        ValuesF[k]=ValuesF[0];
        ValuesF[0]=valueFk;
        k=0;

        poly=InterPolynomial(2, nPtsTotal, points, ValuesF);
        if (poly.NewtonBasis!=NULL) break;

        // the construction of the first polynomial has failed
        // delete[] points;
        free(ValuesF);

        int kbest=findBest(of->data, of); // return 0 if startpoint is given
        Vector BaseStart=of->data.getLine(kbest,dim);
        double vBaseStart=((double**)of->data)[kbest][dim];
        points=GenerateData(rho, BaseStart, vBaseStart, of, &ValuesF);
        nPtsTotal=n;
    }
*/
    // update M:
    fullUpdateOfM(rho,poly.vBase,of->data,poly);
    
    //fprintf(stderr,"init part 2 finished.\n");
    //fprintf(stderr,"init finished.\n");

    // first of init all variables:
    parallelImprove(&poly, &k, rho, &valueFk, poly.vBase);

    // really start in parallel:
    startParallelThread();

    while (true)
    {
//        fprintf(stderr,"rho=%e; fo=%e; NF=%i\n", rho,valueFk,QP_NF);
        while (true)
        {
            // trust region step
            while (true)
            {
//                poly.print();
                parallelImprove(&poly, &k, rho, &valueFk, poly.vBase);

                niter--;
                if ((niter==0)
                   ||(of->isConstrained&&checkForTermination(poly.NewtonPoints[k], poly.vBase, rhoEnd)))
                {
                    poly.vBase+=poly.NewtonPoints[k];
                    //fprintf(stderr,"rho=%e; fo=%e; NF=%i\n", rho,valueFk,of->getNFE());
                    of->valueBest=valueFk;
                    of->xBest=poly.vBase;
                    // to do : compute H and Lambda

                    Vector vG(dim);
                    Matrix mH(dim,dim);
                    poly.gradientHessian(poly.NewtonPoints[k],vG,mH);
                    of->finalize(vG,mH,FullLambda.clone());
                    return;
                }
   
                // to debug:
                //fprintf(stderr,"Best Value Objective=%e (nfe=%i)\n", valueFk, of->getNFE());
                
                d=ConstrainedL2NormMinimizer(poly,k,delta,&info,1000,&lambda1,poly.vBase,of);

//                if (d.euclidianNorm()>delta)
//                {
//                    printf("Warning d to long: (%e > %e)\n", d.euclidianNorm(), delta);
//                }

                normD=mmin(d.euclidianNorm(), delta);
                d+=poly.NewtonPoints[k];

//              next line is equivalent to reduction=valueFk-poly(d); 
//              BUT is more precise (no rounding error)
                reduction=-poly.shiftedEval(d,valueFk);

                //if (normD<0.5*rho) { evalNeeded=true; break; }
                if ((normD<0.5*rho)&&(!forceTRStep)) { evalNeeded=true; break; }

                //  IF THE MODEL REDUCTION IS SMALL, THEN WE DO NOT SAMPLE FUNCTION
                //  AT THE NEW POINT. WE THEN WILL TRY TO IMPROVE THE MODEL.

                noise=0.5*mmax(of->noiseAbsolute*(1+of->noiseRelative), condorAbs(valueFk)*of->noiseRelative);
                if ((reduction<noise)&&(!forceTRStep)) { evalNeeded=true; break; }
                forceTRStep=false; evalNeeded=false;

                if (quickHack) (*quickHack)(poly,k);
                tmp=poly.vBase+d; nerror=0; valueOF=of->eval(tmp, &nerror); 
                of->saveValue(tmp,valueOF,nerror);
                if (nerror)
                {
                    // evaluation failed
                    delta*=0.5;
                    if (normD>=2*rho) continue;
                    break;
                }
                if (!of->isFeasible(tmp, &r))
                {
                    printf("violation: %e\n",r);
                }

                // update of delta:
                r=(valueFk-valueOF)/reduction;
                if (r<=0.1) delta=0.5*normD;
                else if (r<0.7) delta=mmax(0.5*delta, normD);
                     else delta=mmax(rho+ normD, mmax(1.25*normD, delta));
            // powell's heuristics:
                if (delta<1.5*rho) delta=rho;

                if (valueOF<valueFk) 
                {
                    t=poly.findAGoodPointToReplace(-1, rho, d,&modelStep);
                    k=t; valueFk=valueOF; 
                    improvement=true;
//                    fprintf(stderr,"Value Objective=%e\n", valueOF);
                } else
                {
                    t=poly.findAGoodPointToReplace(k, rho, d,&modelStep);
                    improvement=false;
//                    fprintf(stderr,".");
                };

                if (t<0) { poly.updateM(d, valueOF); break; }

                // If we are along constraints, it's more important to update
                // the polynomial with points which increase its quality.
                // Thus, we will skip this update to use only points coming
                // from checkIfValidityIsInBound

                if ((!of->isConstrained)||(improvement)||(reduction>0.0)||(normD<rho)) poly.replace(t, d, valueOF); 

                if (improvement) continue;
//                if (modelStep>4*rho*rho) continue;
                if (modelStep>2*rho) continue;
                if (normD>=2*rho) continue;
                break;
            }
            // model improvement step
            forceTRStep=true;

//            fprintf(stderr,"improvement step\n");
            bound=0.0;
            if (normD<0.5*rho) 
            {
                bound=0.5*sqr(rho)*lambda1;
                if (poly.nUpdateOfM<10) bound=0.0;
            }

            parallelImprove(&poly, &k, rho, &valueFk, poly.vBase);

            // !! change d (if needed):
            t=poly.checkIfValidityIsInBound(d, k, bound, rho );
            if (t>=0)
            {
                if (quickHack) (*quickHack)(poly,k);
                tmp=poly.vBase+d; nerror=0; valueOF=of->eval(tmp, &nerror); 
                if (nerror) 
                {
                    Vector GXk(dim);
                    Matrix H(dim,dim);
                    poly.NewtonBasis[t].gradientHessian(poly.NewtonPoints[k],GXk,H);
                    double rhot=rho,vmax;

                    while (nerror)
                    {
                        rhot*=.5;
                        d=LAGMAXModified(GXk,H,rhot,vmax);
                        d+=poly.NewtonPoints[k];
                        tmp=poly.vBase+d; nerror=0; valueOF=of->eval(tmp, &nerror); 
                        of->saveValue(tmp,valueOF,nerror);
                    }
                }
                poly.replace(t, d, valueOF); 
                if ((valueOF<valueFk)&&
                    (of->isFeasible(tmp))) { k=t; valueFk=valueOF; };
                continue;
            }

            // the model is perfect for this value of rho:
            // OR
            // we have crossed a non_linear constraint which prevent us to advance
            if ((normD<=rho)||(reduction<0.0)) break;
        }

        
        
        // change rho because no improvement can now be made:
        if (rho<=rhoEnd) break;

        //fprintf(stderr,"rho=%e; fo=%e; NF=%i\n", rho,valueFk,of->getNFE());

        if (rho<16*rhoEnd) rhoNew=rhoEnd;
        else if (rho<250*rhoEnd) rhoNew=sqrt(rho*rhoEnd);
             else rhoNew=0.1*rho;
        delta=mmax(0.5*rho,rhoNew);
        rho=rhoNew;

        
        // update of the polynomial: translation of x[k].
        // replace BASE by BASE+x[k]
        poly.translate(poly.NewtonPoints[k]);
    }
    parallelFinish();
    
    Vector vG(dim);
    Matrix mH(dim,dim);

    
    
    if (evalNeeded)
    {
        tmp=poly.vBase+d; nerror=0; valueOF=of->eval(tmp,&nerror); 
        of->saveValue(tmp,valueOF,nerror);
        if ((nerror)||(valueOF<valueFk))
        {
            poly.vBase+=poly.NewtonPoints[k];
            poly.gradientHessian(poly.NewtonPoints[k],vG,mH);
        }
        else
        {
            valueFk=valueOF; poly.vBase=tmp;
            poly.gradientHessian(d,vG,mH);
        }
    } else 
    {
        poly.vBase+=poly.NewtonPoints[k];
        poly.gradientHessian(poly.NewtonPoints[k],vG,mH);
    }


//    delete[] points; :not necessary: done in destructor of poly which is called automatically:
    //fprintf(stderr,"rho=%e; fo=%e; NF=%i\n", rho,valueFk,of->getNFE());
    
    of->valueBest=valueFk;
    of->xBest=poly.vBase;
    of->finalize(vG,mH,FullLambda.clone());
    
}

L2NormMinimizer() [1/3]

Vector L2NormMinimizer	(	Polynomial	q,
		double	delta,
		int *	infoOut = NULL,
		int	maxIter = 1000,
		double *	lambda1 = NULL
	)

Definition at line 307 of file UTRSSolver.cpp.

{
    return L2NormMinimizer(q, Vector::emptyVector, delta, infoOut, maxIter, lambda1);
}

L2NormMinimizer() [2/3]

Vector L2NormMinimizer	(	Polynomial	q,
		Vector	pointXk,
		double	delta,
		int *	infoOut = NULL,
		int	maxIter = 1000,
		double *	lambda1 = NULL
	)

Definition at line 297 of file UTRSSolver.cpp.

{
    int n=q.dim();
    Matrix H(n,n);
    Vector vG(n);
    q.gradientHessian(pointXk,vG,H);
    return L2NormMinimizer(q,pointXk,delta,infoOut,maxIter,lambda1,vG,H);
}

L2NormMinimizer() [3/3]

Vector L2NormMinimizer	(	Polynomial	q,
		Vector	pointXk,
		double	delta,
		int *	infoOut,
		int	maxIter,
		double *	lambda1,
		Vector	minusG,
		Matrix	H
	)

Definition at line 147 of file UTRSSolver.cpp.

{
// lambda1>0.0 if interior convergence

    const double theta=0.01;
//    const double kappaEasy=0.1, kappaHard=0.2;
    const double kappaEasy=0.01, kappaHard=0.02;

    double normG,lambda,lambdaCorrection, lambdaPlus, lambdaL, lambdaU, 
           uHu, alpha, normS;
    int info=0, n=minusG.sz();
    Matrix HLambda(n,n);
    MatrixTriangle L(n);
    Vector s(n), omega(n), u(n), sFinal;
    bool gIsNull, choleskyFactorAlreadyComputed=false;

//    printf("\nG= "); minusG.print();
//    printf("\nH=\n"); H.print();

    gIsNull=minusG.isNull();
    normG=minusG.euclidianNorm();
    lambda=normG/delta;
    minusG.multiply(-1.0);

    lambdaL=initLambdaL(normG,delta,H);
    lambdaU=initLambdaU(normG,delta,H);

    //    Special case: parl = paru.
    lambdaU= mmax(lambdaU,(1+kappaEasy)*lambdaL);

    lambda=mmax(lambda, lambdaL);
    lambda=mmin(lambda, lambdaU);

    while (maxIter--)
    {
        if (!choleskyFactorAlreadyComputed)
        {
            if (!H.cholesky(L, lambda, &lambdaCorrection))
            {
       //         lambdaL=mmax(mmax(lambdaL,lambda),lambdaCorrection);
                lambdaL=mmax(lambdaL,lambda+lambdaCorrection);
                lambda=mmax(sqrt(lambdaL*lambdaU), lambdaL+theta*(lambdaU-lambdaL));
                continue;
            }
        } else choleskyFactorAlreadyComputed=false;


        // cholesky factorization successful : solve Hlambda * s = -G
        s.copyFrom(minusG);
        L.solveInPlace(s);
        L.solveTransposInPlace(s);
        normS=s.euclidianNorm();

        // check for termination
#ifndef POWEL_TERMINATION
        if (condorAbs(normS-delta)<kappaEasy*delta) 
        { 
            s.multiply(delta/normS);
            info=1; 
            break; 
        }
#else
//      powell check !!!
        HLambda.copyFrom(H);
        HLambda.addUnityInPlace(lambda);
        double sHs=s.scalarProduct(HLambda.multiply(s));
        if (sqr(delta/normS-1)<kappaEasy*(1+lambda*delta*delta/sHs)) 
        {
            s.multiply(delta/normS);
            info=1; 
            break;
        }
#endif

        if (normS<delta)
        {
            // check for termination
            // interior convergence; maybe break;
            if (lambda==0) { info=1; break; }
            lambdaU=mmin(lambdaU,lambda); 
        } else lambdaL=mmax(lambdaL,lambda);

//        if (lambdaU-lambdaL<kappaEasy*(2-kappaEasy)*lambdaL) { info=3; break; };

        omega.copyFrom(s);
        L.solveInPlace(omega);
        lambdaPlus=lambda+(normS-delta)/delta*sqr(normS)/sqr(omega.euclidianNorm());
        lambdaPlus=mmax(lambdaPlus, lambdaL);
        lambdaPlus=mmin(lambdaPlus, lambdaU);

        if (normS<delta)
        {
            L.LINPACK(u);
#ifndef POWEL_TERMINATION
            HLambda.copyFrom(H);
            HLambda.addUnityInPlace(lambda);
#endif
            uHu=u.scalarProduct(HLambda.multiply(u));
            lambdaL=mmax(lambdaL,lambda-uHu);

            alpha=findAlpha(s,u,delta,q,pointXk,sFinal,minusG,H);
            // check for termination
#ifndef POWEL_TERMINATION
            if (sqr(alpha)*uHu<
                kappaHard*(s.scalarProduct(HLambda.multiply(s))   ))//  +lambda*sqr(delta)))
#else
            if (sqr(alpha)*uHu+sHs<
                kappaHard*(sHs+lambda*sqr(delta)))
#endif
            {
                s=sFinal; info=2; break;
            }
        }
        if ((normS>delta)&&(!gIsNull)) { lambda=lambdaPlus; continue; };

        if (H.cholesky(L, lambdaPlus, &lambdaCorrection)) 
        { 
            lambda=lambdaPlus; 
            choleskyFactorAlreadyComputed=true; 
            continue; 
        }

        lambdaL=mmax(lambdaL,lambdaPlus);
        // check lambdaL for interior convergence
//      if (lambdaL==0) return s;
        lambda=mmax(sqrt(lambdaL*lambdaU), lambdaL+theta*(lambdaU-lambdaL));
    }

    if (infoOut) *infoOut=info;
    if (lambda1)
    {
        if (lambda==0.0)
        {
            // calculate the value of the lowest eigenvalue of H
            // to check
            lambdaL=0; lambdaU=initLambdaU2(H);
            while (lambdaL<0.99*lambdaU)
            {
                lambda=0.5*(lambdaL+lambdaU);
                if (H.cholesky(L,-lambda)) lambdaL=lambda;
//                if (H.cholesky(L,-lambda,&lambdaCorrection)) lambdaL=lambda+lambdaCorrection;
                else lambdaU=lambda;
            }
            *lambda1=lambdaL;
        } else *lambda1=0.0;
    }
    return s;
}

LAGMAXModified() [1/3]

Vector LAGMAXModified	(	Polynomial	q,
		double	rho,
		double &	VMAX
	)

Definition at line 217 of file MSSolver.cpp.

{
    return LAGMAXModified(q,Vector::emptyVector,rho,VMAX);
}

LAGMAXModified() [2/3]

Vector LAGMAXModified	(	Polynomial	q,
		Vector	pointXk,
		double	rho,
		double &	VMAX
	)

Definition at line 208 of file MSSolver.cpp.

{
    int n=q.dim();
    Matrix H(n,n);
    Vector G(n), D(n);
    q.gradientHessian(pointXk,G,H);
    return LAGMAXModified(G,H,rho,VMAX);
}

LAGMAXModified() [3/3]

Vector LAGMAXModified	(	Vector	G,
		Matrix	H,
		double	rho,
		double &	VMAX
	)

Definition at line 40 of file MSSolver.cpp.

{
    //not tested in depth but running already quite good

//    SUBROUTINE LAGMAX (N,G,H,RHO,D,V,VMAX)
//    IMPLICIT REAL*8 (A-H,O-Z)
//    DIMENSION G(*),H(N,*),D(*),V(*)
//
//    N is the number of variables of a quadratic objective function, Q say.
//    G is the gradient of Q at the origin.
//    H is the symmetric Hessian matrix of Q. Only the upper triangular and
//      diagonal parts need be set.
//    RHO is the trust region radius, and has to be positive.
//    D will be set to the calculated vector of variables.
//    The array V will be used for working space.
//    VMAX will be set to |Q(0)-Q(D)|.
//
//     Calculating the D that maximizes |Q(0)-Q(D)| subject to ||D|| .EQ. RHO
//    requires of order N**3 operations, but sometimes it is adequate if
//    |Q(0)-Q(D)| is within about 0.9 of its greatest possible value. This
//    subroutine provides such a solution in only of order N**2 operations,
//    where the claim of accuracy has been tested by numerical experiments. 
/*
    int n=G.sz();
    Vector D(n), V(n);
    lagmax_(&n, (double *)G, *((double**)H), &rho, 
                    (double*)D, (double*)V, &VMAX);
    return D;
*/
    int i,n=G.sz();
    Vector D(n);

    Vector V=H.getMaxColumn();
    D=H.multiply(V);
    double vv=V.square(),
           dd=D.square(),
           vd=V.scalarProduct(D), 
           dhd=D.scalarProduct(H.multiply(D)),
           *d=D, *v=V, *g=G;

//
//     Set D to a vector in the subspace spanned by V and HV that maximizes
//     |(D,HD)|/(D,D), except that we set D=HV if V and HV are nearly parallel.
//
    if (sqr(vd)<0.9999*vv*dd)
    {
        double a=dhd*vd-dd*dd,
               b=.5*(dhd*vv-dd*vd),
               c=dd*vv-vd*vd,
               tmp1=-b/a;
        if (b*b>a*c)
        {
            double tmp2=sqrt(b*b-a*c)/a, dd1, dd2, dhd1, dhd2;
            Vector D1=D.clone();
            D1.multiply(tmp1+tmp2);
            D1+=V;
            dd1=D1.square();
            dhd1=D1.scalarProduct(H.multiply(D1));

            Vector D2=D.clone();
            D2.multiply(tmp1-tmp2);
            D2+=V;
            dd2=D2.square();
            dhd2=D2.scalarProduct(H.multiply(D2));

            if (condorAbs(dhd1/dd1)>condorAbs(dhd2/dd2)) { D=D1; dd=dd1; dhd=dhd1; } 
            else { D=D2; dd=dd2; dhd=dhd2; }
            d=(double*)D;
        }
    };


//
//     We now turn our attention to the subspace spanned by G and D. A multiple
//     of the current D is returned if that choice seems to be adequate.
//
    double gg=G.square(),
           normG=sqrt(gg),
           gd=G.scalarProduct(D), 
           temp=gd/gg,
           scale=sign(rho/sqrt(dd), gd*dhd);

    i=n; while (i--) v[i]=d[i]-temp*g[i];
    vv=V.square();

    if ((normG*dd)<(0.5-2*rho*condorAbs(dhd))||(vv/dd<1e-4))
    {
        D.multiply(scale);
        VMAX=condorAbs(scale*(gd+0.5*scale*dhd));
        return D;
    }

//
//     G and V are now orthogonal in the subspace spanned by G and D. Hence
//     we generate an orthonormal basis of this subspace such that (D,HV) is
//     negligible or zero, where D and V will be the basis vectors.
//
    H.multiply(D,G); //  D=HG;
    double ghg=G.scalarProduct(D),
           vhg=V.scalarProduct(D),
           vhv=V.scalarProduct(H.multiply(V));
    double theta, cosTheta, sinTheta;

    if (condorAbs(vhg)<0.01*mmax(condorAbs(vhv),condorAbs(ghg)))
    {
        cosTheta=1.0;
        sinTheta=0.0;
    } else
    {
        theta=0.5*atan(0.5*vhg/(vhv-ghg));
        cosTheta=cos(theta);
        sinTheta=sin(theta);
    }
    i=n;
    while(i--)
    {
        d[i]= cosTheta*g[i]+ sinTheta*v[i];
        v[i]=-sinTheta*g[i]+ cosTheta*v[i];
    };

//
//     The final D is a multiple of the current D, V, D+V or D-V. We make the
//     choice from these possibilities that is optimal.
//

    double norm=rho/D.euclidianNorm();
    D.multiply(norm);
    dhd=(ghg*sqr(cosTheta)+vhv*sqr(sinTheta))*sqr(norm);

    norm=rho/V.euclidianNorm();
    V.multiply(norm);
    vhv=(ghg*sqr(sinTheta)+vhv*sqr(cosTheta)*sqr(norm));

    double halfRootTwo=sqrt(0.5),   // =sqrt(2)/2=cos(PI/4)
           t1=normG*cosTheta*rho,   // t1=condorAbs(D.scalarProduct(G));
           t2=normG*sinTheta*rho,   // t2=condorAbs(V.scalarProduct(G));
           at1=condorAbs(t1),
           at2=condorAbs(t2),
           t3=0.25*(dhd+vhv),
           q1=condorAbs(at1+0.5*dhd),
           q2=condorAbs(at2+0.5*vhv),
           q3=condorAbs(halfRootTwo*(at1+at2)+t3),
           q4=condorAbs(halfRootTwo*(at1-at2)+t3);
    if ((q4>q3)&&(q4>q2)&&(q4>q1))
    {
        double st1=sign(t1*t3), st2=sign(t2*t3);
        i=n; while (i--) d[i]=halfRootTwo*(st1*d[i]-st2*v[i]);
        VMAX=q4;
        return D;
    }
    if ((q3>q2)&&(q3>q1))
    {
        double st1=sign(t1*t3), st2=sign(t2*t3);
        i=n; while (i--) d[i]=halfRootTwo*(st1*d[i]+st2*v[i]);
        VMAX=q3;
        return D;
    }
    if (q2>q1)
    {
        if (t2*vhv<0) V.multiply(-1);
        VMAX=q2;
        return V;
    }
    if (t1*dhd<0) D.multiply(-1);
    VMAX=q1;
    return D;
}