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
018package org.apache.commons.math3.linear;
019
020import org.apache.commons.math3.exception.DimensionMismatchException;
021import 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 1566017 2014-02-08 14:13:34Z tn $
047 * @since 2.0 (changed to concrete class in 3.0)
048 */
049public 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        /**
382         * Get the inverse of the decomposed matrix.
383         *
384         * @return the inverse matrix.
385         * @throws SingularMatrixException if the decomposed matrix is singular.
386         */
387        public RealMatrix getInverse() {
388            return solve(MatrixUtils.createRealIdentityMatrix(pivot.length));
389        }
390    }
391}