Example usage for org.apache.commons.math3.exception.util LocalizedFormats UNABLE_TO_PERFORM_QR_DECOMPOSITION_ON_JACOBIAN

Introduction

In this page you can find the example usage for org.apache.commons.math3.exception.util LocalizedFormats UNABLE_TO_PERFORM_QR_DECOMPOSITION_ON_JACOBIAN.

Prototype

LocalizedFormats UNABLE_TO_PERFORM_QR_DECOMPOSITION_ON_JACOBIAN

To view the source code for org.apache.commons.math3.exception.util LocalizedFormats UNABLE_TO_PERFORM_QR_DECOMPOSITION_ON_JACOBIAN.

Click Source Link

Usage

From source file:be.ugent.maf.cellmissy.analysis.doseresponse.impl.DoseResponseLMOptimizer.java

/**
 * Copy of super source code. Decompose a matrix A as A.P = Q.R using
 * Householder transforms.// w w w. j  ava 2 s  . c  om
 * <p>
 * As suggested in the P. Lascaux and R. Theodor book
 * <i>Analyse num&eacute;rique matricielle appliqu&eacute;e &agrave; l'art
 * de l'ing&eacute;nieur</i> (Masson, 1986), instead of representing the
 * Householder transforms with u<sub>k</sub> unit vectors such that:
 * <pre>
 * H<sub>k</sub> = I - 2u<sub>k</sub>.u<sub>k</sub><sup>t</sup>
 * </pre> we use <sub>k</sub> non-unit vectors such that:
 * <pre>
 * H<sub>k</sub> = I - beta<sub>k</sub>v<sub>k</sub>.v<sub>k</sub><sup>t</sup>
 * </pre> where v<sub>k</sub> = a<sub>k</sub> - alpha<sub>k</sub>
 * e<sub>k</sub>. The beta<sub>k</sub> coefficients are provided upon exit
 * as recomputing them from the v<sub>k</sub> vectors would be costly.</p>
 * <p>
 * This decomposition handles rank deficient cases since the tranformations
 * are performed in non-increasing columns norms order thanks to columns
 * pivoting. The diagonal elements of the R matrix are therefore also in
 * non-increasing absolute values order.</p>
 *
 * @param jacobian Weighted Jacobian matrix at the current point.
 * @param solvedCols Number of solved point.
 * @return data used in other methods of this class.
 * @throws ConvergenceException if the decomposition cannot be performed.
 */
private InternalData qrDecomposition(RealMatrix jacobian, int solvedCols) throws ConvergenceException {
    // Code in this class assumes that the weighted Jacobian is -(W^(1/2) J),
    // hence the multiplication by -1.
    final double[][] weightedJacobian = jacobian.scalarMultiply(-1).getData();

    final int nR = weightedJacobian.length;
    final int nC = weightedJacobian[0].length;

    final int[] permutation = new int[nC];
    final double[] diagR = new double[nC];
    final double[] jacNorm = new double[nC];
    final double[] beta = new double[nC];

    // initializations
    for (int k = 0; k < nC; ++k) {
        permutation[k] = k;
        double norm2 = 0;
        for (int i = 0; i < nR; ++i) {
            double akk = weightedJacobian[i][k];
            norm2 += akk * akk;
        }
        jacNorm[k] = FastMath.sqrt(norm2);
    }

    // transform the matrix column after column
    for (int k = 0; k < nC; ++k) {

        // select the column with the greatest norm on active components
        int nextColumn = -1;
        double ak2 = Double.NEGATIVE_INFINITY;
        for (int i = k; i < nC; ++i) {
            double norm2 = 0;
            for (int j = k; j < nR; ++j) {
                double aki = weightedJacobian[j][permutation[i]];
                norm2 += aki * aki;
            }
            if (Double.isInfinite(norm2) || Double.isNaN(norm2)) {
                throw new ConvergenceException(LocalizedFormats.UNABLE_TO_PERFORM_QR_DECOMPOSITION_ON_JACOBIAN,
                        nR, nC);
            }
            if (norm2 > ak2) {
                nextColumn = i;
                ak2 = norm2;
            }
        }
        if (ak2 <= getRankingThreshold()) {
            return new InternalData(weightedJacobian, permutation, k, diagR, jacNorm, beta);
        }
        int pk = permutation[nextColumn];
        permutation[nextColumn] = permutation[k];
        permutation[k] = pk;

        // choose alpha such that Hk.u = alpha ek
        double akk = weightedJacobian[k][pk];
        double alpha = (akk > 0) ? -FastMath.sqrt(ak2) : FastMath.sqrt(ak2);
        double betak = 1.0 / (ak2 - akk * alpha);
        beta[pk] = betak;

        // transform the current column
        diagR[pk] = alpha;
        weightedJacobian[k][pk] -= alpha;

        // transform the remaining columns
        for (int dk = nC - 1 - k; dk > 0; --dk) {
            double gamma = 0;
            for (int j = k; j < nR; ++j) {
                gamma += weightedJacobian[j][pk] * weightedJacobian[j][permutation[k + dk]];
            }
            gamma *= betak;
            for (int j = k; j < nR; ++j) {
                weightedJacobian[j][permutation[k + dk]] -= gamma * weightedJacobian[j][pk];
            }
        }
    }

    return new InternalData(weightedJacobian, permutation, solvedCols, diagR, jacNorm, beta);
}

---