optimization.java

一个数据挖掘软件ALPHAMINERR的整个过程的JAVA版源代码

	    alam = 1.0;
	
	// Iteration of one newton step, if necessary, backtracking is done
	double initF=fold, // Initial function value
	    hi=alam, lo=alam, newSlope=0, fhi=m_f, flo=m_f;// Variables used for beta condition
	double[] newGrad;  // Gradient on the new variable values
	
	kloop:
	for (k=0;;k++) {
	    if(m_Debug)
		System.err.println("\nLine search iteration: " + k);
	    
	    for (i=0;i<len;i++){
		if(!isFixed[i]){
		    x[i] = xold[i]+alam*direct[i];  // Compute xnew
		    if(!Double.isNaN(nwsBounds[0][i]) && (x[i]<nwsBounds[0][i])){    
			x[i] = nwsBounds[0][i]; //Rounding error	
		    }
		    else if(!Double.isNaN(nwsBounds[1][i]) && (x[i]>nwsBounds[1][i])){		
			x[i] = nwsBounds[1][i]; //Rounding error	
		    }
		}
	    }
	    
	    m_f = objectiveFunction(x);    // Compute fnew
	    if(Double.isNaN(m_f))
		throw new Exception("Objective function value is NaN!");
	
	    while(Double.isInfinite(m_f)){ // Avoid infinity
		if(m_Debug)
		    System.err.println("Too large m_f.  Shrink step by half.");
		alam *= 0.5; // Shrink by half
		if(alam <= m_Epsilon){
		    if(m_Debug)
			System.err.println("Wrong starting points, change them!");
		    return x;
		}
		
		for (i=0;i<len;i++)
		    if(!isFixed[i])
			x[i] = xold[i]+alam*direct[i]; 
		
		m_f = objectiveFunction(x); 
		if(Double.isNaN(m_f))
		    throw new Exception("Objective function value is NaN!");
	
		initF = Double.POSITIVE_INFINITY;
	    }
	    
	    if(m_Debug) {
		System.err.println("obj. function: " + 
				   Utils.doubleToString(m_f, 10, 7));
		System.err.println("threshold: " + 
				   Utils.doubleToString(fold+m_ALF*alam*m_Slope,10,7));
	    }
	    
	    if(m_f<=fold+m_ALF*alam*m_Slope){// Alpha condition: sufficient function decrease
		if(m_Debug)		
		    System.err.println("Sufficient function decrease (alpha condition): "); 
		newGrad = evaluateGradient(x);
		for(newSlope=0.0,i=0; i<len; i++)
		    if(!isFixed[i])
			newSlope += newGrad[i]*direct[i];

		if(newSlope >= m_BETA*m_Slope){ // Beta condition: ensure pos. defnty.	
		    if(m_Debug)		
			System.err.println("Increasing derivatives (beta condition): "); 	

		    if((fixedOne!=-1) && (alam>=alpha)){ // Has bounds and over
			if(direct[fixedOne] > 0){
			    x[fixedOne] = nwsBounds[1][fixedOne]; // Avoid rounding error
			    nwsBounds[1][fixedOne]=Double.NaN; //Add cons. to working set
			}
			else{
			    x[fixedOne] = nwsBounds[0][fixedOne]; // Avoid rounding error
			    nwsBounds[0][fixedOne]=Double.NaN; //Add cons. to working set
			}
			
			if(m_Debug)
			    System.err.println("Fix variable "
					       +fixedOne+" to bound "+ x[fixedOne]+
					       " from value "+ xold[fixedOne]);
			isFixed[fixedOne]=true; // Fix the variable
			wsBdsIndx.addElement(new Integer(fixedOne));
		    }		
		    return x;
		}
		else if(k==0){ // First time: increase alam 
		    // Search for the smallest value not complying with alpha condition
		    double upper = Math.min(alpha,maxalam); 
		    if(m_Debug)
			System.err.println("Alpha condition holds, increase alpha... ");
		    while(!((alam>=upper) || (m_f>fold+m_ALF*alam*m_Slope))){
			lo = alam;
			flo = m_f;
			alam *= 2.0;
			if(alam>=upper)  // Avoid rounding errors
			    alam=upper;

			for (i=0;i<len;i++)
			    if(!isFixed[i])
				x[i] = xold[i]+alam*direct[i];
			m_f = objectiveFunction(x);
			if(Double.isNaN(m_f))
			    throw new Exception("Objective function value is NaN!");
			
			newGrad = evaluateGradient(x);
			for(newSlope=0.0,i=0; i<len; i++)
			    if(!isFixed[i])
				newSlope += newGrad[i]*direct[i];
			
			if(newSlope >= m_BETA*m_Slope){
			    if (m_Debug)		
				System.err.println("Increasing derivatives (beta condition): \n"+
						   "newSlope = "+Utils.doubleToString(newSlope,10,7));
			    
			    if((fixedOne!=-1) && (alam>=alpha)){ // Has bounds and over
				if(direct[fixedOne] > 0){
				    x[fixedOne] = nwsBounds[1][fixedOne]; // Avoid rounding error
				    nwsBounds[1][fixedOne]=Double.NaN; //Add cons. to working set
				}
				else{
				    x[fixedOne] = nwsBounds[0][fixedOne]; // Avoid rounding error
				    nwsBounds[0][fixedOne]=Double.NaN; //Add cons. to working set
				}
				
				if(m_Debug)
				    System.err.println("Fix variable "
						       +fixedOne+" to bound "+ x[fixedOne]+
						       " from value "+ xold[fixedOne]);
				isFixed[fixedOne]=true; // Fix the variable
				wsBdsIndx.addElement(new Integer(fixedOne));
			    }		 				    
			    return x;
			}
		    }
		    hi = alam;
		    fhi = m_f;			
		    break kloop;
		}
		else{
		    if(m_Debug)
			System.err.println("Alpha condition holds.");
		    hi = alam2; lo = alam; flo = m_f;
		    break kloop;
		}		    
	    }        
	    else if (alam < alamin) { // No feasible lambda found       
		if(initF<fold){ 
		    alam = Math.min(1.0,alpha);
		    for (i=0;i<len;i++)
			if(!isFixed[i])
			    x[i] = xold[i]+alam*direct[i]; //Still take Alpha
		    
		    if (m_Debug)
			System.err.println("No feasible lambda: still take"+
					   " alpha="+alam);
		    
		    if((fixedOne!=-1) && (alam>=alpha)){ // Has bounds and over
			if(direct[fixedOne] > 0){
			    x[fixedOne] = nwsBounds[1][fixedOne]; // Avoid rounding error
			    nwsBounds[1][fixedOne]=Double.NaN; //Add cons. to working set
			}
			else{
			    x[fixedOne] = nwsBounds[0][fixedOne]; // Avoid rounding error
			    nwsBounds[0][fixedOne]=Double.NaN; //Add cons. to working set
			}
			
			if(m_Debug)
			    System.err.println("Fix variable "
					       +fixedOne+" to bound "+ x[fixedOne]+
					       " from value "+ xold[fixedOne]);
			isFixed[fixedOne]=true; // Fix the variable
			wsBdsIndx.addElement(new Integer(fixedOne));
		    }		 		    
		}
		else{   // Convergence on delta(x)
		    for(i=0;i<len;i++) 
			x[i]=xold[i];
		    m_f=fold;
		    if (m_Debug)
			System.err.println("Cannot find feasible lambda"); 
		}
		
		return x; 
	    }
	    else { // Backtracking by polynomial interpolation
		if(k==0){ // First time backtrack: quadratic interpolation
		    if(!Double.isInfinite(initF))
			initF = m_f;		    
		    // lambda = -g'(0)/(2*g''(0))
		    tmplam = -0.5*alam*m_Slope/((m_f-fold)/alam-m_Slope);
		}
		else {    // Subsequent backtrack: cubic interpolation 
		    rhs1 = m_f-fold-alam*m_Slope;
		    rhs2 = fhi-fold-alam2*m_Slope;
		    a=(rhs1/(alam*alam)-rhs2/(alam2*alam2))/(alam-alam2);
		    b=(-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2);
		    if (a == 0.0) tmplam = -m_Slope/(2.0*b);
		    else {
			disc=b*b-3.0*a*m_Slope;
			if (disc < 0.0) disc = 0.0;
			double numerator = -b+Math.sqrt(disc);
			if(numerator >= Double.MAX_VALUE){
			    numerator = Double.MAX_VALUE;
			    if (m_Debug)
				System.err.print("-b+sqrt(disc) too large! Set it to MAX_VALUE.");
			}
			tmplam=numerator/(3.0*a);
		    }
		    if (m_Debug)
			System.err.print("Cubic interpolation: \n" + 
					 "a:   " + Utils.doubleToString(a,10,7)+ "\n" +
					 "b:   " + Utils.doubleToString(b,10,7)+ "\n" +
					 "disc:   " + Utils.doubleToString(disc,10,7)+ "\n" +
					 "tmplam:   " + tmplam + "\n" +
					 "alam:   " + Utils.doubleToString(alam,10,7)+ "\n");	
		    if (tmplam>0.5*alam)
			tmplam=0.5*alam;             // lambda <= 0.5*lambda_old
		}
	    }
	    alam2=alam;
	    fhi=m_f;
	    alam=Math.max(tmplam,0.1*alam);          // lambda >= 0.1*lambda_old
	    
	    if(alam>alpha){
		throw new Exception("Sth. wrong in lnsrch:"+
				    "Lambda infeasible!(lambda="+alam+
				    ", alpha="+alpha+", upper="+tmplam+
				    "|"+(-alpha*m_Slope/(2.0*((m_f-fold)/alpha-m_Slope)))+
				    ", m_f="+m_f+", fold="+fold+
				    ", slope="+m_Slope);
	    }	    
	} // Endfor(k=0;;k++)
	
	// Quadratic interpolation between lamda values between lo and hi.
	// If cannot find a value satisfying beta condition, use lo.
	double ldiff = hi-lo, lincr;
	if(m_Debug)
	    System.err.println("Last stage of searching for beta condition (alam between "
			       +Utils.doubleToString(lo,10,7)+" and "
			       +Utils.doubleToString(hi,10,7)+")...\n"+
			       "Quadratic Interpolation(QI):\n"+
			       "Last newSlope = "+Utils.doubleToString(newSlope, 10, 7));
	
	while((newSlope<m_BETA*m_Slope) && (ldiff>=alamin)){
	    lincr = -0.5*newSlope*ldiff*ldiff/(fhi-flo-newSlope*ldiff);

	    if(m_Debug)
		System.err.println("fhi = "+fhi+"\n"+
				   "flo = "+flo+"\n"+
				   "ldiff = "+ldiff+"\n"+
				   "lincr (using QI) = "+lincr+"\n");
	    
	    if(lincr<0.2*ldiff) lincr=0.2*ldiff;
	    alam = lo+lincr;
	    if(alam >= hi){ // We cannot go beyond the bounds, so the best we can try is hi
	    	alam=hi;
		lincr=ldiff;
	    }
	    for (i=0;i<len;i++)
		if(!isFixed[i])
		    x[i] = xold[i]+alam*direct[i];
	    m_f = objectiveFunction(x);
	    if(Double.isNaN(m_f))
		throw new Exception("Objective function value is NaN!");
	
	    if(m_f>fold+m_ALF*alam*m_Slope){ 
		// Alpha condition fails, shrink lambda_upper
		ldiff = lincr;
		fhi = m_f;
	    }	    
	    else{ // Alpha condition holds	    
		newGrad = evaluateGradient(x);
		for(newSlope=0.0,i=0; i<len; i++)
		    if(!isFixed[i])
			newSlope += newGrad[i]*direct[i];
		
		if(newSlope < m_BETA*m_Slope){
		    // Beta condition fails, shrink lambda_lower
		    lo = alam;
		    ldiff -= lincr;
		    flo = m_f;
		}
	    }
	} 
	
	if(newSlope < m_BETA*m_Slope){ // Cannot satisfy beta condition, take lo
	    if(m_Debug)
		System.err.println("Beta condition cannot be satisfied, take alpha condition");
	    alam=lo;
	    for (i=0;i<len;i++)
		if(!isFixed[i])
		    x[i] = xold[i]+alam*direct[i];
	    m_f = flo;
	}
	else if(m_Debug)
	    System.err.println("Both alpha and beta conditions are satisfied. alam="
			       +Utils.doubleToString(alam,10,7));
	
	if((fixedOne!=-1) && (alam>=alpha)){ // Has bounds and over
	    if(direct[fixedOne] > 0){
		x[fixedOne] = nwsBounds[1][fixedOne]; // Avoid rounding error
		nwsBounds[1][fixedOne]=Double.NaN; //Add cons. to working set
	    }
	    else{
		x[fixedOne] = nwsBounds[0][fixedOne]; // Avoid rounding error
		nwsBounds[0][fixedOne]=Double.NaN; //Add cons. to working set
	    }
	    
	    if(m_Debug)
		System.err.println("Fix variable "
				   +fixedOne+" to bound "+ x[fixedOne]+
				   " from value "+ xold[fixedOne]);
	    isFixed[fixedOne]=true; // Fix the variable
	    wsBdsIndx.addElement(new Integer(fixedOne));
	}
	
	return x;
    }
    
    /**
     * Main algorithm.  Descriptions see "Practical Optimization"
     *
     * @param initX initial point of x, assuming no value's on the bound!
     * @param constraints the bound constraints of each variable
     *                    constraints[0] is the lower bounds and 
     *                    constraints[1] is the upper bounds
     * @return the solution of x, null if number of iterations not enough
     * @exception Exception if an error occurs
     */
    public double[] findArgmin(double[] initX, double[][] constraints) 
	throws Exception{
	int l = initX.length;
	
	// Initially all variables are free, all bounds are constraints of
	// non-working-set constraints
	boolean[] isFixed = new boolean[l];
	double[][] nwsBounds = new double[2][l];
	// Record indice of fixed variables, simply for efficiency
	FastVector wsBdsIndx = new FastVector(); 
	// Vectors used to record the variable indices to be freed 	
	FastVector toFree=null, oldToFree=null;	

	// Initial value of obj. function, gradient and inverse of the Hessian
	m_f = objectiveFunction(initX);
	if(Double.isNaN(m_f))
	    throw new Exception("Objective function value is NaN!");
	
	double sum=0;
	double[] grad=evaluateGradient(initX), oldGrad, oldX,
	    deltaGrad=new double[l], deltaX=new double[l],
	    direct = new double[l], x = new double[l];
	Matrix lowerTriangleMatrix = new Matrix(l, l);  // Lower triangle of Cholesky factor 
	double[] diagonalOfCholesky = new double[l];   // Diagonal of Cholesky factor
	
    // Variables used for finding new direction 
	// TWang. Mar 23. 2005.
    Matrix tempLD = new Matrix(l,l); // L*D
    double[] tempDiagonal = new double[l];
    // >> End. TWang
    
	for(int i=0; i<l; i++){
	    lowerTriangleMatrix.setRow(i, new double[l]);
	    lowerTriangleMatrix.setElement(i,i,1.0);
	    diagonalOfCholesky[i] = 1.0;
	    direct[i] = -grad[i];
	    sum += grad[i]*grad[i];
	    x[i] = initX[i];
	    nwsBounds[0][i] = constraints[0][i];
	    nwsBounds[1][i] = constraints[1][i];
	    isFixed[i] = false;
	}	
	double stpmax = m_STPMX*Math.max(Math.sqrt(sum), l);
	
	iterates:
	for(int step=0; step < m_MAXITS; step++){
	    if (m_Debug)
		System.err.println("\nIteration # " + step + ":");	    
	    
	    // Try at most one feasible newton step, i.e. 0<lamda<=alpha
	    oldX = x;
	    oldGrad = grad;
	    
	    // Also update grad
	    if (m_Debug)
		System.err.println("Line search ... ");
	    m_IsZeroStep = false;
	    x=lnsrch(x, grad, direct, stpmax, 
		     isFixed, nwsBounds, wsBdsIndx);
	    if (m_Debug)
		System.err.println("Line search finished.");
	    
	    if(m_IsZeroStep){ // Zero step, simply delete rows/cols of D and L
		for(int f=0; f<wsBdsIndx.size(); f++){
		    int idx=((Integer)wsBdsIndx.elementAt(f)).intValue();
		    lowerTriangleMatrix.setRow(idx, new double[l]);
		    lowerTriangleMatrix.setColumn(idx, new double[l]);
		    diagonalOfCholesky[idx] = 0.0;
		}		
		grad = evaluateGradient(x);
		step--;
	    }
	    else{