com.joptimizer.algebra.Matrix1NormRescaler.java Source code

Java tutorial

Introduction

Here is the source code for com.joptimizer.algebra.Matrix1NormRescaler.java

Source

/*
 * 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.logging.Log;
import org.apache.commons.logging.LogFactory;

import com.joptimizer.util.ColtUtils;

import cern.colt.function.IntIntDoubleFunction;
import cern.colt.matrix.DoubleFactory1D;
import cern.colt.matrix.DoubleFactory2D;
import cern.colt.matrix.DoubleMatrix1D;
import cern.colt.matrix.DoubleMatrix2D;
import cern.colt.matrix.linalg.Algebra;

/**
 * Calculates the matrix rescaling factors so that the 1-norm of each row and each column of the scaled matrix
 * asymptotically converges to one.
 * 
 * @author alberto trivellato (alberto.trivellato@gmail.com)
 * @see Daniel Ruiz, "A scaling algorithm to equilibrate both rows and columns norms in matrices"
 * @see Philip A. Knight, Daniel Ruiz, Bora Ucar "A Symmetry Preserving Algorithm for Matrix Scaling"
 */
public class Matrix1NormRescaler implements MatrixRescaler {

    private double eps = 1.e-3;
    private Log log = LogFactory.getLog(this.getClass().getName());

    public Matrix1NormRescaler() {
    }

    public Matrix1NormRescaler(double eps) {
        this.eps = eps;
    }

    /**
     * Scaling factors for not singular matrices.
     * @see Daniel Ruiz, "A scaling algorithm to equilibrate both rows and columns norms in matrices"
     * @see Philip A. Knight, Daniel Ruiz, Bora Ucar "A Symmetry Preserving Algorithm for Matrix Scaling"
     */
    public DoubleMatrix1D[] getMatrixScalingFactors(DoubleMatrix2D A) {
        DoubleFactory1D F1 = DoubleFactory1D.dense;
        DoubleFactory2D F2 = DoubleFactory2D.sparse;
        Algebra ALG = Algebra.DEFAULT;
        int r = A.rows();
        int c = A.columns();
        DoubleMatrix1D D1 = F1.make(r, 1);
        DoubleMatrix1D D2 = F1.make(c, 1);
        DoubleMatrix2D AK = A.copy();
        DoubleMatrix1D DR = F1.make(r, 1);
        DoubleMatrix1D DC = F1.make(c, 1);
        DoubleMatrix1D DRInv = F1.make(r);
        DoubleMatrix1D DCInv = F1.make(c);
        log.debug("eps  : " + eps);
        int maxIteration = 50;
        for (int k = 0; k <= maxIteration; k++) {
            double normR = -Double.MAX_VALUE;
            double normC = -Double.MAX_VALUE;
            for (int i = 0; i < r; i++) {
                double dri = ALG.normInfinity(AK.viewRow(i));
                DR.setQuick(i, Math.sqrt(dri));
                DRInv.setQuick(i, 1. / Math.sqrt(dri));
                normR = Math.max(normR, Math.abs(1 - dri));
            }
            for (int j = 0; j < c; j++) {
                double dci = ALG.normInfinity(AK.viewColumn(j));
                DC.setQuick(j, Math.sqrt(dci));
                DCInv.setQuick(j, 1. / Math.sqrt(dci));
                normC = Math.max(normC, Math.abs(1 - dci));
            }

            log.debug("normR: " + normR);
            log.debug("normC: " + normC);
            if (normR < eps && normC < eps) {
                break;
            }

            //D1 = ALG.mult(D1, DRInv);
            for (int i = 0; i < r; i++) {
                double prevD1I = D1.getQuick(i);
                double newD1I = prevD1I * DRInv.getQuick(i);
                D1.setQuick(i, newD1I);
            }
            //D2 = ALG.mult(D2, DCInv);
            for (int j = 0; j < c; j++) {
                double prevD2J = D2.getQuick(j);
                double newD2J = prevD2J * DCInv.getQuick(j);
                D2.setQuick(j, newD2J);
            }
            //log.debug("D1: " + ArrayUtils.toString(D1.toArray()));
            //log.debug("D2: " + ArrayUtils.toString(D2.toArray()));

            if (k == maxIteration) {
                log.warn("max iteration reached");
            }

            //AK = ALG.mult(DRInv, ALG.mult(AK, DCInv));
            AK = ColtUtils.diagonalMatrixMult(DRInv, AK, DCInv);
        }

        return new DoubleMatrix1D[] { D1, D2 };
    }

    /**
     * Scaling factors for symmetric (not singular) matrices.
     * Just the subdiagonal elements of the matrix are required.
     * @see Daniel Ruiz, "A scaling algorithm to equilibrate both rows and columns norms in matrices"
     * @see Philip A. Knight, Daniel Ruiz, Bora Ucar "A Symmetry Preserving Algorithm for Matrix Scaling"
     */
    public DoubleMatrix1D getMatrixScalingFactorsSymm(DoubleMatrix2D A) {
        DoubleFactory1D F1 = DoubleFactory1D.dense;
        DoubleFactory2D F2 = DoubleFactory2D.sparse;
        Algebra ALG = Algebra.DEFAULT;
        int dim = A.columns();
        DoubleMatrix1D D1 = F1.make(dim, 1);
        DoubleMatrix2D AK = A.copy();
        DoubleMatrix2D DR = F2.identity(dim);
        DoubleMatrix1D DRInv = F1.make(dim);
        //log.debug("eps  : " + eps);
        int maxIteration = 50;
        for (int k = 0; k <= maxIteration; k++) {
            double normR = -Double.MAX_VALUE;
            for (int i = 0; i < dim; i++) {
                //double dri = ALG.normInfinity(AK.viewRow(i));
                double dri = this.getRowInfinityNorm(AK, i);
                DR.setQuick(i, i, Math.sqrt(dri));
                DRInv.setQuick(i, 1. / Math.sqrt(dri));
                normR = Math.max(normR, Math.abs(1 - dri));
                if (Double.isNaN(normR)) {
                    throw new IllegalArgumentException("matrix is singular");
                }
            }

            //log.debug("normR: " + normR);
            if (normR < eps) {
                break;
            }

            for (int i = 0; i < dim; i++) {
                double prevD1I = D1.getQuick(i);
                double newD1I = prevD1I * DRInv.getQuick(i);
                D1.setQuick(i, newD1I);
            }
            //logger.debug("D1: " + ArrayUtils.toString(D1.toArray()));

            if (k == maxIteration) {
                log.warn("max iteration reached");
            }

            AK = ColtUtils.diagonalMatrixMult(DRInv, AK, DRInv);
        }

        return D1;
    }

    /**
     * Check if the scaling algorithm returned proper results.
     * Note that AOriginal cannot be only subdiagonal filled, because this check
     * is for both symm and bath notsymm matrices.
     * @param AOriginal the ORIGINAL (before scaling) matrix
     * @param U the return of the scaling algorithm
     * @param V the return of the scaling algorithm
     * @param base
     * @return
     */
    public boolean checkScaling(final DoubleMatrix2D AOriginal, final DoubleMatrix1D U, final DoubleMatrix1D V) {

        int c = AOriginal.columns();
        int r = AOriginal.rows();
        final double[] maxValueHolder = new double[] { -Double.MAX_VALUE };

        IntIntDoubleFunction myFunct = new IntIntDoubleFunction() {
            public double apply(int i, int j, double pij) {
                maxValueHolder[0] = Math.max(maxValueHolder[0], Math.abs(pij));
                return pij;
            }
        };

        DoubleMatrix2D AScaled = ColtUtils.diagonalMatrixMult(U, AOriginal, V);

        //view A row by row
        boolean isOk = true;
        for (int i = 0; isOk && i < r; i++) {
            maxValueHolder[0] = -Double.MAX_VALUE;
            DoubleMatrix2D P = AScaled.viewPart(i, 0, 1, c);
            P.forEachNonZero(myFunct);
            isOk = Math.abs(1. - maxValueHolder[0]) < eps;
        }
        //view A col by col
        for (int j = 0; isOk && j < c; j++) {
            maxValueHolder[0] = -Double.MAX_VALUE;
            DoubleMatrix2D P = AScaled.viewPart(0, j, r, 1);
            P.forEachNonZero(myFunct);
            isOk = Math.abs(1. - maxValueHolder[0]) < eps;
        }
        return isOk;
    }

    /**
     * 
     * @param ASymm symm matrix filled in its subdiagonal elements
     * @param r the index of the row
     * @return
     */
    public double getRowInfinityNorm(final DoubleMatrix2D ASymm, final int r) {

        final double[] maxValueHolder = new double[] { -Double.MAX_VALUE };

        IntIntDoubleFunction myFunct = new IntIntDoubleFunction() {
            public double apply(int i, int j, double pij) {
                //logger.warn("(" + i + "," + j + ")=" + pij);
                maxValueHolder[0] = Math.max(maxValueHolder[0], Math.abs(pij));
                return pij;
            }
        };

        //view A row from starting element to diagonal
        DoubleMatrix2D AR = ASymm.viewPart(r, 0, 1, r + 1);
        AR.forEachNonZero(myFunct);
        //view A col from diagonal to final element
        DoubleMatrix2D AC = ASymm.viewPart(r, r, ASymm.rows() - r, 1);
        AC.forEachNonZero(myFunct);

        return maxValueHolder[0];
    }

}