    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017
    018    package org.apache.commons.math3.linear;
    019
    020    import org.apache.commons.math3.exception.DimensionMismatchException;
    021    import org.apache.commons.math3.util.FastMath;
    022
    023    /**
    024     * Calculates the LUP-decomposition of a square matrix.
    025     * <p>The LUP-decomposition of a matrix A consists of three matrices L, U and
    026     * P that satisfy: P&times;A = L&times;U. L is lower triangular (with unit
    027     * diagonal terms), U is upper triangular and P is a permutation matrix. All
    028     * matrices are m&times;m.</p>
    029     * <p>As shown by the presence of the P matrix, this decomposition is
    030     * implemented using partial pivoting.</p>
    031     * <p>This class is based on the class with similar name from the
    032     * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library.</p>
    033     * <ul>
    034     *   <li>a {@link #getP() getP} method has been added,</li>
    035     *   <li>the {@code det} method has been renamed as {@link #getDeterminant()
    036     *   getDeterminant},</li>
    037     *   <li>the {@code getDoublePivot} method has been removed (but the int based
    038     *   {@link #getPivot() getPivot} method has been kept),</li>
    039     *   <li>the {@code solve} and {@code isNonSingular} methods have been replaced
    040     *   by a {@link #getSolver() getSolver} method and the equivalent methods
    041     *   provided by the returned {@link DecompositionSolver}.</li>
    042     * </ul>
    043     *
    044     * @see <a href="http://mathworld.wolfram.com/LUDecomposition.html">MathWorld</a>
    045     * @see <a href="http://en.wikipedia.org/wiki/LU_decomposition">Wikipedia</a>
    046     * @version $Id: LUDecomposition.java 1416643 2012-12-03 19:37:14Z tn $
    047     * @since 2.0 (changed to concrete class in 3.0)
    048     */
    049    public class LUDecomposition {
    050        /** Default bound to determine effective singularity in LU decomposition. */
    051        private static final double DEFAULT_TOO_SMALL = 1e-11;
    052        /** Entries of LU decomposition. */
    053        private final double[][] lu;
    054        /** Pivot permutation associated with LU decomposition. */
    055        private final int[] pivot;
    056        /** Parity of the permutation associated with the LU decomposition. */
    057        private boolean even;
    058        /** Singularity indicator. */
    059        private boolean singular;
    060        /** Cached value of L. */
    061        private RealMatrix cachedL;
    062        /** Cached value of U. */
    063        private RealMatrix cachedU;
    064        /** Cached value of P. */
    065        private RealMatrix cachedP;
    066
    067        /**
    068         * Calculates the LU-decomposition of the given matrix.
    069         * This constructor uses 1e-11 as default value for the singularity
    070         * threshold.
    071         *
    072         * @param matrix Matrix to decompose.
    073         * @throws NonSquareMatrixException if matrix is not square.
    074         */
    075        public LUDecomposition(RealMatrix matrix) {
    076            this(matrix, DEFAULT_TOO_SMALL);
    077        }
    078
    079        /**
    080         * Calculates the LU-decomposition of the given matrix.
    081         * @param matrix The matrix to decompose.
    082         * @param singularityThreshold threshold (based on partial row norm)
    083         * under which a matrix is considered singular
    084         * @throws NonSquareMatrixException if matrix is not square
    085         */
    086        public LUDecomposition(RealMatrix matrix, double singularityThreshold) {
    087            if (!matrix.isSquare()) {
    088                throw new NonSquareMatrixException(matrix.getRowDimension(),
    089                                                   matrix.getColumnDimension());
    090            }
    091
    092            final int m = matrix.getColumnDimension();
    093            lu = matrix.getData();
    094            pivot = new int[m];
    095            cachedL = null;
    096            cachedU = null;
    097            cachedP = null;
    098
    099            // Initialize permutation array and parity
    100            for (int row = 0; row < m; row++) {
    101                pivot[row] = row;
    102            }
    103            even     = true;
    104            singular = false;
    105
    106            // Loop over columns
    107            for (int col = 0; col < m; col++) {
    108
    109                // upper
    110                for (int row = 0; row < col; row++) {
    111                    final double[] luRow = lu[row];
    112                    double sum = luRow[col];
    113                    for (int i = 0; i < row; i++) {
    114                        sum -= luRow[i] * lu[i][col];
    115                    }
    116                    luRow[col] = sum;
    117                }
    118
    119                // lower
    120                int max = col; // permutation row
    121                double largest = Double.NEGATIVE_INFINITY;
    122                for (int row = col; row < m; row++) {
    123                    final double[] luRow = lu[row];
    124                    double sum = luRow[col];
    125                    for (int i = 0; i < col; i++) {
    126                        sum -= luRow[i] * lu[i][col];
    127                    }
    128                    luRow[col] = sum;
    129
    130                    // maintain best permutation choice
    131                    if (FastMath.abs(sum) > largest) {
    132                        largest = FastMath.abs(sum);
    133                        max = row;
    134                    }
    135                }
    136
    137                // Singularity check
    138                if (FastMath.abs(lu[max][col]) < singularityThreshold) {
    139                    singular = true;
    140                    return;
    141                }
    142
    143                // Pivot if necessary
    144                if (max != col) {
    145                    double tmp = 0;
    146                    final double[] luMax = lu[max];
    147                    final double[] luCol = lu[col];
    148                    for (int i = 0; i < m; i++) {
    149                        tmp = luMax[i];
    150                        luMax[i] = luCol[i];
    151                        luCol[i] = tmp;
    152                    }
    153                    int temp = pivot[max];
    154                    pivot[max] = pivot[col];
    155                    pivot[col] = temp;
    156                    even = !even;
    157                }
    158
    159                // Divide the lower elements by the "winning" diagonal elt.
    160                final double luDiag = lu[col][col];
    161                for (int row = col + 1; row < m; row++) {
    162                    lu[row][col] /= luDiag;
    163                }
    164            }
    165        }
    166
    167        /**
    168         * Returns the matrix L of the decomposition.
    169         * <p>L is a lower-triangular matrix</p>
    170         * @return the L matrix (or null if decomposed matrix is singular)
    171         */
    172        public RealMatrix getL() {
    173            if ((cachedL == null) && !singular) {
    174                final int m = pivot.length;
    175                cachedL = MatrixUtils.createRealMatrix(m, m);
    176                for (int i = 0; i < m; ++i) {
    177                    final double[] luI = lu[i];
    178                    for (int j = 0; j < i; ++j) {
    179                        cachedL.setEntry(i, j, luI[j]);
    180                    }
    181                    cachedL.setEntry(i, i, 1.0);
    182                }
    183            }
    184            return cachedL;
    185        }
    186
    187        /**
    188         * Returns the matrix U of the decomposition.
    189         * <p>U is an upper-triangular matrix</p>
    190         * @return the U matrix (or null if decomposed matrix is singular)
    191         */
    192        public RealMatrix getU() {
    193            if ((cachedU == null) && !singular) {
    194                final int m = pivot.length;
    195                cachedU = MatrixUtils.createRealMatrix(m, m);
    196                for (int i = 0; i < m; ++i) {
    197                    final double[] luI = lu[i];
    198                    for (int j = i; j < m; ++j) {
    199                        cachedU.setEntry(i, j, luI[j]);
    200                    }
    201                }
    202            }
    203            return cachedU;
    204        }
    205
    206        /**
    207         * Returns the P rows permutation matrix.
    208         * <p>P is a sparse matrix with exactly one element set to 1.0 in
    209         * each row and each column, all other elements being set to 0.0.</p>
    210         * <p>The positions of the 1 elements are given by the {@link #getPivot()
    211         * pivot permutation vector}.</p>
    212         * @return the P rows permutation matrix (or null if decomposed matrix is singular)
    213         * @see #getPivot()
    214         */
    215        public RealMatrix getP() {
    216            if ((cachedP == null) && !singular) {
    217                final int m = pivot.length;
    218                cachedP = MatrixUtils.createRealMatrix(m, m);
    219                for (int i = 0; i < m; ++i) {
    220                    cachedP.setEntry(i, pivot[i], 1.0);
    221                }
    222            }
    223            return cachedP;
    224        }
    225
    226        /**
    227         * Returns the pivot permutation vector.
    228         * @return the pivot permutation vector
    229         * @see #getP()
    230         */
    231        public int[] getPivot() {
    232            return pivot.clone();
    233        }
    234
    235        /**
    236         * Return the determinant of the matrix
    237         * @return determinant of the matrix
    238         */
    239        public double getDeterminant() {
    240            if (singular) {
    241                return 0;
    242            } else {
    243                final int m = pivot.length;
    244                double determinant = even ? 1 : -1;
    245                for (int i = 0; i < m; i++) {
    246                    determinant *= lu[i][i];
    247                }
    248                return determinant;
    249            }
    250        }
    251
    252        /**
    253         * Get a solver for finding the A &times; X = B solution in exact linear
    254         * sense.
    255         * @return a solver
    256         */
    257        public DecompositionSolver getSolver() {
    258            return new Solver(lu, pivot, singular);
    259        }
    260
    261        /** Specialized solver. */
    262        private static class Solver implements DecompositionSolver {
    263
    264            /** Entries of LU decomposition. */
    265            private final double[][] lu;
    266
    267            /** Pivot permutation associated with LU decomposition. */
    268            private final int[] pivot;
    269
    270            /** Singularity indicator. */
    271            private final boolean singular;
    272
    273            /**
    274             * Build a solver from decomposed matrix.
    275             * @param lu entries of LU decomposition
    276             * @param pivot pivot permutation associated with LU decomposition
    277             * @param singular singularity indicator
    278             */
    279            private Solver(final double[][] lu, final int[] pivot, final boolean singular) {
    280                this.lu       = lu;
    281                this.pivot    = pivot;
    282                this.singular = singular;
    283            }
    284
    285            /** {@inheritDoc} */
    286            public boolean isNonSingular() {
    287                return !singular;
    288            }
    289
    290            /** {@inheritDoc} */
    291            public RealVector solve(RealVector b) {
    292                final int m = pivot.length;
    293                if (b.getDimension() != m) {
    294                    throw new DimensionMismatchException(b.getDimension(), m);
    295                }
    296                if (singular) {
    297                    throw new SingularMatrixException();
    298                }
    299
    300                final double[] bp = new double[m];
    301
    302                // Apply permutations to b
    303                for (int row = 0; row < m; row++) {
    304                    bp[row] = b.getEntry(pivot[row]);
    305                }
    306
    307                // Solve LY = b
    308                for (int col = 0; col < m; col++) {
    309                    final double bpCol = bp[col];
    310                    for (int i = col + 1; i < m; i++) {
    311                        bp[i] -= bpCol * lu[i][col];
    312                    }
    313                }
    314
    315                // Solve UX = Y
    316                for (int col = m - 1; col >= 0; col--) {
    317                    bp[col] /= lu[col][col];
    318                    final double bpCol = bp[col];
    319                    for (int i = 0; i < col; i++) {
    320                        bp[i] -= bpCol * lu[i][col];
    321                    }
    322                }
    323
    324                return new ArrayRealVector(bp, false);
    325            }
    326
    327            /** {@inheritDoc} */
    328            public RealMatrix solve(RealMatrix b) {
    329
    330                final int m = pivot.length;
    331                if (b.getRowDimension() != m) {
    332                    throw new DimensionMismatchException(b.getRowDimension(), m);
    333                }
    334                if (singular) {
    335                    throw new SingularMatrixException();
    336                }
    337
    338                final int nColB = b.getColumnDimension();
    339
    340                // Apply permutations to b
    341                final double[][] bp = new double[m][nColB];
    342                for (int row = 0; row < m; row++) {
    343                    final double[] bpRow = bp[row];
    344                    final int pRow = pivot[row];
    345                    for (int col = 0; col < nColB; col++) {
    346                        bpRow[col] = b.getEntry(pRow, col);
    347                    }
    348                }
    349
    350                // Solve LY = b
    351                for (int col = 0; col < m; col++) {
    352                    final double[] bpCol = bp[col];
    353                    for (int i = col + 1; i < m; i++) {
    354                        final double[] bpI = bp[i];
    355                        final double luICol = lu[i][col];
    356                        for (int j = 0; j < nColB; j++) {
    357                            bpI[j] -= bpCol[j] * luICol;
    358                        }
    359                    }
    360                }
    361
    362                // Solve UX = Y
    363                for (int col = m - 1; col >= 0; col--) {
    364                    final double[] bpCol = bp[col];
    365                    final double luDiag = lu[col][col];
    366                    for (int j = 0; j < nColB; j++) {
    367                        bpCol[j] /= luDiag;
    368                    }
    369                    for (int i = 0; i < col; i++) {
    370                        final double[] bpI = bp[i];
    371                        final double luICol = lu[i][col];
    372                        for (int j = 0; j < nColB; j++) {
    373                            bpI[j] -= bpCol[j] * luICol;
    374                        }
    375                    }
    376                }
    377
    378                return new Array2DRowRealMatrix(bp, false);
    379            }
    380
    381            /** {@inheritDoc} */
    382            public RealMatrix getInverse() {
    383                return solve(MatrixUtils.createRealIdentityMatrix(pivot.length));
    384            }
    385        }
    386    }