com.joptimizer.algebra.MatrixLogSumRescaler Maven / Gradle / Ivy

Go to download
Show more of this group Show more artifacts with this name
Show all versions of joptimizer Show documentation
JOptimizer is a java library for solving general convex optimization problems, like for example LP, QP, QCQP, SOCP, SDP, GP and many more.
The newest version!
/*
 * Copyright 2011-2014 JOptimizer
 *
 *   Licensed under the Apache License, Version 2.0 (the "License");
 *   you may not use this file except in compliance with the License.
 *   You may obtain a copy of the License at
 *
 *       http://www.apache.org/licenses/LICENSE-2.0
 *
 *   Unless required by applicable law or agreed to in writing, software
 *   distributed under the License is distributed on an "AS IS" BASIS,
 *   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *   See the License for the specific language governing permissions and
 *   limitations under the License.
 */
package com.joptimizer.algebra;

import org.apache.commons.lang3.ArrayUtils;
import org.apache.log4j.Logger;

import cern.colt.function.IntIntDoubleFunction;
import cern.colt.matrix.DoubleMatrix1D;
import cern.colt.matrix.DoubleMatrix2D;
import cern.colt.matrix.impl.DenseDoubleMatrix1D;

/**
 * Calculates the matrix rescaling factors so that the scaled
 * matrix has its entries near to unity in the sense that the sum of the
 * squares of the logarithms of the entries is minimized.
 * 
 * @TODO: add documentation
 * @author alberto trivellato ([email protected])
 * @see Xin Huang, Preprocessing and Postprocessing in Linear Optimization
 * @see A. Chang and J.K. Reid. On the automatic scaling of matrices for Gaussian elimination, 
 * 		Journal of the Institute of Mathematics and Its Applications 10 (1972)
 * @see Gajulapalli, Lasdon "Scaling Sparse Matrices for Optimization Algorithms"
 */
public class MatrixLogSumRescaler implements MatrixRescaler {

	private double base = 10;
	private Logger logger = Logger.getLogger(this.getClass().getName());
	
	public MatrixLogSumRescaler(){
	}
	
	public MatrixLogSumRescaler(double base){
		this.base = base;
	}
	
	/**
	 * Gauss-Seidel scaling for a sparse matrix: 
	 * 
AScaled = U.A.V, with A mxn matrix, U, V diagonal.
	 * Returns the two scaling matrices
	 * 
U[i,i] = base^x[i], i=0,...,m
	 * 
V[i,i] = base^y[i], i=0,...,n
	 * 
	 * @see Gajulapalli, Lasdon "Scaling Sparse Matrices for Optimization Algorithms", algorithms 1 and 2
	 */
	@Override
	public DoubleMatrix1D[] getMatrixScalingFactors(DoubleMatrix2D A) {
		int m = A.rows();
		int n = A.columns();
		final double log10_b = Math.log10(base);

		//Setup for Gauss-Seidel Iterations
		final int[] R = new int[m];
		final int[] C = new int[n];
		final double[] t = new double[1];
		final double[] a = new double[m];
		final double[] b = new double[n];
		final boolean[][] Z = new boolean[m][n];
		A.forEachNonZero(new IntIntDoubleFunction() {
			@Override
			public double apply(int i, int j, double aij) {
				R[i] = R[i] + 1;
				C[j] = C[j] + 1;
				Z[i][j] = true;
				t[0] = -(Math.log10(Math.abs(aij)) / log10_b + 0.5);// log(b, x) = log(k, x) / log(k, b)
				a[i] = a[i] + t[0];
				b[j] = b[j] + t[0];
				return aij;
			}
		});
		
		for(int i=0; i currentColumnIndex) {
				  //sub-diagonal elements
					//log.debug("i:" + i + ", j:" + currentColumnIndex + ": " + pij);
					tHolder[0] = tHolder[0] - 2 * (Math.log10(Math.abs(pij))/log10_b + 0.5) -2*x[i];//log(b, x) = log(k, x) / log(k, b)
					cHolder[0] = cHolder[0] + 2;//- 2*x[i]
				}
				return pij;
			}
		};
		
		//view A column by column
		for (int currentColumnIndex = n - 1; currentColumnIndex >= 0; currentColumnIndex--) {
			//log.debug("currentColumnIndex:" + currentColumnIndex);
			cHolder[0] = 0;//reset
			tHolder[0] = 0;//reset
			currentColumnIndexHolder[0] = currentColumnIndex;
			DoubleMatrix2D P = A.viewPart(0, currentColumnIndex, n, 1);
			P.forEachNonZero(myFunct);
			if(cHolder[0] > 0){
				x[currentColumnIndex] = (int)Math.round(tHolder[0] / cHolder[0]);
			}
		}
		
		//log.debug("x: " + ArrayUtils.toString(x));
		
		DoubleMatrix1D u = new DenseDoubleMatrix1D(n);
		for (int k = 0; k < n; k++) {
			u.setQuick(k, Math.pow(base, x[k]));
		}
		return u;
	}
	
	/**
	 * Check if the scaling algorithm returned proper results.
	 * Note that the scaling algorithm is for minimizing a given objective function of the original matrix elements, and
	 * the check will be done on the value of this objective function. 
	 * @param A the ORIGINAL (before scaling) matrix
	 * @param U the return of the scaling algorithm
	 * @param V the return of the scaling algorithm
	 * @param base
	 * @return
	 */
	@Override
	public boolean checkScaling(final DoubleMatrix2D A, 
			final DoubleMatrix1D U, final DoubleMatrix1D V){
		
		final double log10_2 = Math.log10(base);
		final double[] originalOFValue = {0};
		final double[] scaledOFValue = {0};
		final double[] x = new double[A.rows()];
		final double[] y = new double[A.columns()];
		
		A.forEachNonZero(new IntIntDoubleFunction() {
			@Override
			public double apply(int i, int j, double aij) {
				double v = Math.log10(Math.abs(aij)) / log10_2 + 0.5;
				originalOFValue[0] = originalOFValue[0] + Math.pow(v, 2);
				double xi = Math.log10(U.getQuick(i)) / log10_2;
				double yj = Math.log10(V.getQuick(j)) / log10_2;
				scaledOFValue[0] = scaledOFValue[0] + Math.pow(xi + yj + v, 2);
				x[i] = xi;
				y[j] = yj;
				return aij;
			}
		});
		
		originalOFValue[0] = 0.5 * originalOFValue[0];
		scaledOFValue[0]   = 0.5 * scaledOFValue[0];
		
		logger.debug("x: " + ArrayUtils.toString(x));
		logger.debug("y: " + ArrayUtils.toString(y));
		logger.debug("originalOFValue: " + originalOFValue[0]);
		logger.debug("scaledOFValue  : " + scaledOFValue[0]);
		return !(originalOFValue[0] < scaledOFValue[0]);
	}

}