LUDecomposition.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      https://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */

  18. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */

  23. package org.hipparchus.linear;

  25. import org.hipparchus.exception.LocalizedCoreFormats;
  26. import org.hipparchus.exception.MathIllegalArgumentException;
  27. import org.hipparchus.util.FastMath;

  29. /**
  30.  * Calculates the LUP-decomposition of a square matrix.
  31.  * <p>The LUP-decomposition of a matrix A consists of three matrices L, U and
  32.  * P that satisfy: P&times;A = L&times;U. L is lower triangular (with unit
  33.  * diagonal terms), U is upper triangular and P is a permutation matrix. All
  34.  * matrices are m&times;m.</p>
  35.  * <p>As shown by the presence of the P matrix, this decomposition is
  36.  * implemented using partial pivoting.</p>
  37.  * <p>This class is based on the class with similar name from the
  38.  * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library.</p>
  39.  * <ul>
  40.  *   <li>a {@link #getP() getP} method has been added,</li>
  41.  *   <li>the {@code det} method has been renamed as {@link #getDeterminant()
  42.  *   getDeterminant},</li>
  43.  *   <li>the {@code getDoublePivot} method has been removed (but the int based
  44.  *   {@link #getPivot() getPivot} method has been kept),</li>
  45.  *   <li>the {@code solve} and {@code isNonSingular} methods have been replaced
  46.  *   by a {@link #getSolver() getSolver} method and the equivalent methods
  47.  *   provided by the returned {@link DecompositionSolver}.</li>
  48.  * </ul>
  49.  *
  50.  * @see <a href="http://mathworld.wolfram.com/LUDecomposition.html">MathWorld</a>
  51.  * @see <a href="http://en.wikipedia.org/wiki/LU_decomposition">Wikipedia</a>
  52.  */
  53. public class LUDecomposition {
  54.     /** Default bound to determine effective singularity in LU decomposition. */
  55.     private static final double DEFAULT_TOO_SMALL = 1e-11;
  56.     /** Entries of LU decomposition. */
  57.     private final double[][] lu;
  58.     /** Pivot permutation associated with LU decomposition. */
  59.     private final int[] pivot;
  60.     /** Parity of the permutation associated with the LU decomposition. */
  61.     private boolean even;
  62.     /** Singularity indicator. */
  63.     private boolean singular;
  64.     /** Cached value of L. */
  65.     private RealMatrix cachedL;
  66.     /** Cached value of U. */
  67.     private RealMatrix cachedU;
  68.     /** Cached value of P. */
  69.     private RealMatrix cachedP;

  71.     /**
  72.      * Calculates the LU-decomposition of the given matrix.
  73.      * This constructor uses 1e-11 as default value for the singularity
  74.      * threshold.
  75.      *
  76.      * @param matrix Matrix to decompose.
  77.      * @throws MathIllegalArgumentException if matrix is not square.
  78.      */
  79.     public LUDecomposition(RealMatrix matrix) {
  80.         this(matrix, DEFAULT_TOO_SMALL);
  81.     }

  83.     /**
  84.      * Calculates the LU-decomposition of the given matrix.
  85.      * @param matrix The matrix to decompose.
  86.      * @param singularityThreshold threshold (based on partial row norm)
  87.      * under which a matrix is considered singular
  88.      * @throws MathIllegalArgumentException if matrix is not square
  89.      */
  90.     public LUDecomposition(RealMatrix matrix, double singularityThreshold) {
  91.         if (!matrix.isSquare()) {
  92.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NON_SQUARE_MATRIX,
  93.                                                    matrix.getRowDimension(), matrix.getColumnDimension());
  94.         }

  96.         final int m = matrix.getColumnDimension();
  97.         lu = matrix.getData();
  98.         pivot = new int[m];
  99.         cachedL = null;
  100.         cachedU = null;
  101.         cachedP = null;

  103.         // Initialize permutation array and parity
  104.         for (int row = 0; row < m; row++) {
  105.             pivot[row] = row;
  106.         }
  107.         even     = true;
  108.         singular = false;

  110.         // Loop over columns
  111.         for (int col = 0; col < m; col++) {

  113.             // upper
  114.             for (int row = 0; row < col; row++) {
  115.                 final double[] luRow = lu[row];
  116.                 double sum = luRow[col];
  117.                 for (int i = 0; i < row; i++) {
  118.                     sum -= luRow[i] * lu[i][col];
  119.                 }
  120.                 luRow[col] = sum;
  121.             }

  123.             // lower
  124.             int max = col; // permutation row
  125.             double largest = Double.NEGATIVE_INFINITY;
  126.             for (int row = col; row < m; row++) {
  127.                 final double[] luRow = lu[row];
  128.                 double sum = luRow[col];
  129.                 for (int i = 0; i < col; i++) {
  130.                     sum -= luRow[i] * lu[i][col];
  131.                 }
  132.                 luRow[col] = sum;

  134.                 // maintain best permutation choice
  135.                 if (FastMath.abs(sum) > largest) {
  136.                     largest = FastMath.abs(sum);
  137.                     max = row;
  138.                 }
  139.             }

  141.             // Singularity check
  142.             if (FastMath.abs(lu[max][col]) < singularityThreshold) {
  143.                 singular = true;
  144.                 return;
  145.             }

  147.             // Pivot if necessary
  148.             if (max != col) {
  149.                 final double[] luMax = lu[max];
  150.                 final double[] luCol = lu[col];
  151.                 for (int i = 0; i < m; i++) {
  152.                     final double tmp = luMax[i];
  153.                     luMax[i] = luCol[i];
  154.                     luCol[i] = tmp;
  155.                 }
  156.                 int temp = pivot[max];
  157.                 pivot[max] = pivot[col];
  158.                 pivot[col] = temp;
  159.                 even = !even;
  160.             }

  162.             // Divide the lower elements by the "winning" diagonal elt.
  163.             final double luDiag = lu[col][col];
  164.             for (int row = col + 1; row < m; row++) {
  165.                 lu[row][col] /= luDiag;
  166.             }
  167.         }
  168.     }

  170.     /**
  171.      * Returns the matrix L of the decomposition.
  172.      * <p>L is a lower-triangular matrix</p>
  173.      * @return the L matrix (or null if decomposed matrix is singular)
  174.      */
  175.     public RealMatrix getL() {
  176.         if ((cachedL == null) && !singular) {
  177.             final int m = pivot.length;
  178.             cachedL = MatrixUtils.createRealMatrix(m, m);
  179.             for (int i = 0; i < m; ++i) {
  180.                 final double[] luI = lu[i];
  181.                 for (int j = 0; j < i; ++j) {
  182.                     cachedL.setEntry(i, j, luI[j]);
  183.                 }
  184.                 cachedL.setEntry(i, i, 1.0);
  185.             }
  186.         }
  187.         return cachedL;
  188.     }

  190.     /**
  191.      * Returns the matrix U of the decomposition.
  192.      * <p>U is an upper-triangular matrix</p>
  193.      * @return the U matrix (or null if decomposed matrix is singular)
  194.      */
  195.     public RealMatrix getU() {
  196.         if ((cachedU == null) && !singular) {
  197.             final int m = pivot.length;
  198.             cachedU = MatrixUtils.createRealMatrix(m, m);
  199.             for (int i = 0; i < m; ++i) {
  200.                 final double[] luI = lu[i];
  201.                 for (int j = i; j < m; ++j) {
  202.                     cachedU.setEntry(i, j, luI[j]);
  203.                 }
  204.             }
  205.         }
  206.         return cachedU;
  207.     }

  209.     /**
  210.      * Returns the P rows permutation matrix.
  211.      * <p>P is a sparse matrix with exactly one element set to 1.0 in
  212.      * each row and each column, all other elements being set to 0.0.</p>
  213.      * <p>The positions of the 1 elements are given by the {@link #getPivot()
  214.      * pivot permutation vector}.</p>
  215.      * @return the P rows permutation matrix (or null if decomposed matrix is singular)
  216.      * @see #getPivot()
  217.      */
  218.     public RealMatrix getP() {
  219.         if ((cachedP == null) && !singular) {
  220.             final int m = pivot.length;
  221.             cachedP = MatrixUtils.createRealMatrix(m, m);
  222.             for (int i = 0; i < m; ++i) {
  223.                 cachedP.setEntry(i, pivot[i], 1.0);
  224.             }
  225.         }
  226.         return cachedP;
  227.     }

  229.     /**
  230.      * Returns the pivot permutation vector.
  231.      * @return the pivot permutation vector
  232.      * @see #getP()
  233.      */
  234.     public int[] getPivot() {
  235.         return pivot.clone();
  236.     }

  238.     /**
  239.      * Return the determinant of the matrix
  240.      * @return determinant of the matrix
  241.      */
  242.     public double getDeterminant() {
  243.         if (singular) {
  244.             return 0;
  245.         } else {
  246.             final int m = pivot.length;
  247.             double determinant = even ? 1 : -1;
  248.             for (int i = 0; i < m; i++) {
  249.                 determinant *= lu[i][i];
  250.             }
  251.             return determinant;
  252.         }
  253.     }

  255.     /**
  256.      * Get a solver for finding the A &times; X = B solution in exact linear
  257.      * sense.
  258.      * @return a solver
  259.      */
  260.     public DecompositionSolver getSolver() {
  261.         return new Solver();
  262.     }

  264.     /** Specialized solver. */
  265.     private class Solver implements DecompositionSolver {

  267.         /** {@inheritDoc} */
  268.         @Override
  269.         public boolean isNonSingular() {
  270.             return !singular;
  271.         }

  273.         /** {@inheritDoc} */
  274.         @Override
  275.         public RealVector solve(RealVector b) {
  276.             final int m = pivot.length;
  277.             if (b.getDimension() != m) {
  278.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  279.                                                        b.getDimension(), m);
  280.             }
  281.             if (singular) {
  282.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.SINGULAR_MATRIX);
  283.             }

  285.             final double[] bp = new double[m];

  287.             // Apply permutations to b
  288.             for (int row = 0; row < m; row++) {
  289.                 bp[row] = b.getEntry(pivot[row]);
  290.             }

  292.             // Solve LY = b
  293.             for (int col = 0; col < m; col++) {
  294.                 final double bpCol = bp[col];
  295.                 for (int i = col + 1; i < m; i++) {
  296.                     bp[i] -= bpCol * lu[i][col];
  297.                 }
  298.             }

  300.             // Solve UX = Y
  301.             for (int col = m - 1; col >= 0; col--) {
  302.                 bp[col] /= lu[col][col];
  303.                 final double bpCol = bp[col];
  304.                 for (int i = 0; i < col; i++) {
  305.                     bp[i] -= bpCol * lu[i][col];
  306.                 }
  307.             }

  309.             return new ArrayRealVector(bp, false);
  310.         }

  312.         /** {@inheritDoc} */
  313.         @Override
  314.         public RealMatrix solve(RealMatrix b) {

  316.             final int m = pivot.length;
  317.             if (b.getRowDimension() != m) {
  318.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  319.                                                        b.getRowDimension(), m);
  320.             }
  321.             if (singular) {
  322.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.SINGULAR_MATRIX);
  323.             }

  325.             final int nColB = b.getColumnDimension();

  327.             // Apply permutations to b
  328.             final double[][] bp = new double[m][nColB];
  329.             for (int row = 0; row < m; row++) {
  330.                 final double[] bpRow = bp[row];
  331.                 final int pRow = pivot[row];
  332.                 for (int col = 0; col < nColB; col++) {
  333.                     bpRow[col] = b.getEntry(pRow, col);
  334.                 }
  335.             }

  337.             // Solve LY = b
  338.             for (int col = 0; col < m; col++) {
  339.                 final double[] bpCol = bp[col];
  340.                 for (int i = col + 1; i < m; i++) {
  341.                     final double[] bpI = bp[i];
  342.                     final double luICol = lu[i][col];
  343.                     for (int j = 0; j < nColB; j++) {
  344.                         bpI[j] -= bpCol[j] * luICol;
  345.                     }
  346.                 }
  347.             }

  349.             // Solve UX = Y
  350.             for (int col = m - 1; col >= 0; col--) {
  351.                 final double[] bpCol = bp[col];
  352.                 final double luDiag = lu[col][col];
  353.                 for (int j = 0; j < nColB; j++) {
  354.                     bpCol[j] /= luDiag;
  355.                 }
  356.                 for (int i = 0; i < col; i++) {
  357.                     final double[] bpI = bp[i];
  358.                     final double luICol = lu[i][col];
  359.                     for (int j = 0; j < nColB; j++) {
  360.                         bpI[j] -= bpCol[j] * luICol;
  361.                     }
  362.                 }
  363.             }

  365.             return new Array2DRowRealMatrix(bp, false);
  366.         }

  368.         /**
  369.          * Get the inverse of the decomposed matrix.
  370.          *
  371.          * @return the inverse matrix.
  372.          * @throws MathIllegalArgumentException if the decomposed matrix is singular.
  373.          */
  374.         @Override
  375.         public RealMatrix getInverse() {
  376.             return solve(MatrixUtils.createRealIdentityMatrix(pivot.length));
  377.         }

  379.         /** {@inheritDoc} */
  380.         @Override
  381.         public int getRowDimension() {
  382.             return lu.length;
  383.         }

  385.         /** {@inheritDoc} */
  386.         @Override
  387.         public int getColumnDimension() {
  388.             return lu[0].length;
  389.         }

  391.     }

  393. }
Created with JaCoCo 0.8.12.202403310830