UTRSSolver.cpp File Reference

#include <stdio.h>
#include <memory.h>
#include "Matrix.h"
#include "tools.h"
#include "Poly.h"

Include dependency graph for UTRSSolver.cpp:

Go to the source code of this file.

Functions
double	findAlpha (Vector s, Vector u, double delta, Polynomial &q, Vector pointXk, Vector &output, Vector minusG, Matrix H)
double	initLambdaL (double normG, double delta, Matrix H)
double	initLambdaU (double normG, double delta, Matrix H)
double	initLambdaU2 (Matrix H)
Vector	L2NormMinimizer (Polynomial q, Vector pointXk, double delta, int *infoOut, int maxIter, double *lambda1, Vector minusG, Matrix H)
Vector	L2NormMinimizer (Polynomial q, Vector pointXk, double delta, int *infoOut, int maxIter, double *lambda1)
Vector	L2NormMinimizer (Polynomial q, double delta, int *infoOut, int maxIter, double *lambda1)

Function Documentation

findAlpha()

double findAlpha	(	Vector	s,
		Vector	u,
		double	delta,
		Polynomial &	q,
		Vector	pointXk,
		Vector &	output,
		Vector	minusG,
		Matrix	H
	)

Definition at line 35 of file UTRSSolver.cpp.

{
    static Vector v1, v2, tmp;
    double a=0,b=0,c=-sqr(delta), *sp=s, *up=u;
    int n=s.sz();
    while (n--)
    {
        a+=sqr(*up); b+=*up * *sp; c+=sqr(*sp);
        sp++; up++;
    }
    double tmp1=-b/a, tmp2= sqrt(b*b-a*c)/a, q1, q2;

    n=s.sz();
    v1.setSize(n); v2.setSize(n); tmp.setSize(n);
    if (!(q==Polynomial::emptyPolynomial))
    {
        v1.copyFrom(u);
        v1.multiply(tmp1+tmp2);
        v1+=s;
        tmp.copyFrom(v1);
        tmp+=pointXk;
        q1=q(tmp);

    // !!! don't do this:
    //    output=v1;
    //    return tmp1+tmp2;

        v2.copyFrom(u);
        v2.multiply(tmp1-tmp2);
        v2+=s;
        tmp.copyFrom(v2);
        tmp+=pointXk;
        q2=q(tmp);

    } else
    {
        v1.copyFrom(u);
        v1.multiply(tmp1+tmp2);
        v1+=s;
        H.multiply(tmp,v1);
        q1=-minusG.scalarProduct(v1)+0.5*v1.scalarProduct(tmp);

        v2.copyFrom(u);
        v2.multiply(tmp1-tmp2);
        v2+=s;
        H.multiply(tmp,v2);
        q2=-minusG.scalarProduct(v2)+0.5*v2.scalarProduct(tmp);
    }
    if (q1>q2) { output=v1; return tmp1+tmp2; }
    output=v2; return tmp1-tmp2;
}

initLambdaL()

double initLambdaL	(	double	normG,
		double	delta,
		Matrix	H
	)

Definition at line 90 of file UTRSSolver.cpp.

{
    int n=H.nLine(),i,j;
    double **h=H, sum,l,a=INF;

    for (i=0; i<n; i++) a=mmin(a,h[i][i]);
    l=mmax(0.0,-a);

    a=0;
    for (i=0; i<n; i++) 
    {
        sum=h[i][i];
        for (j=0; j<n; j++) if (j!=i) sum+=condorAbs(h[i][j]);
        a=mmax(a,sum);
    }
    a=mmin(a,H.frobeniusNorm());
    a=mmin(a,H.LnftyNorm());

    l=mmax(l,normG/delta-a);
    return l;
}

initLambdaU()

double initLambdaU	(	double	normG,
		double	delta,
		Matrix	H
	)

Definition at line 112 of file UTRSSolver.cpp.

{
    int n=H.nLine(),i,j;
    double **h=H, sum,l,a=-INF;

    for (i=0; i<n; i++) 
    {
        sum=-h[i][i];
        for (j=0; j<n; j++) if (j!=i) sum+=condorAbs(h[i][j]);
        a=mmax(a,sum);
    }
    a=mmin(a,H.frobeniusNorm());
    a=mmin(a,H.LnftyNorm());

    l=mmax(0.0,normG/delta+a);
    return l;
}

initLambdaU2()

double initLambdaU2

(

Matrix

H

)

Definition at line 130 of file UTRSSolver.cpp.

{
    int n=H.nLine(),i,j;
    double **h=H, sum,a=-INF;

    for (i=0; i<n; i++) 
    {
        sum=h[i][i];
        for (j=0; j<n; j++) if (j!=i) sum+=condorAbs(h[i][j]);
        a=mmax(a,sum);
    }
    a=mmin(a,H.frobeniusNorm());
    return mmin(a,H.LnftyNorm());
}

L2NormMinimizer() [1/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;
}

L2NormMinimizer() [2/3]

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

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,
		double	delta,
		int *	infoOut,
		int	maxIter,
		double *	lambda1
	)

Definition at line 307 of file UTRSSolver.cpp.

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