FieldLUDecomposition.java
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You 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 org.apache.commons.math4.legacy.field.linalg;
import org.apache.commons.numbers.field.Field;
import org.apache.commons.math4.legacy.linear.SingularMatrixException;
/**
* Calculates the LUP-decomposition of a square matrix.
*
* <p>The LUP-decomposition of a matrix A consists of three matrices
* L, U and P that satisfy: PA = LU, L is lower triangular, and U is
* upper triangular and P is a permutation matrix. All matrices are
* m×m.</p>
*
* <p>Since {@link Field field} elements do not provide an ordering
* operator, the permutation matrix is computed here only in order to
* avoid a zero pivot element, no attempt is done to get the largest
* pivot element.</p>
*
* @param <T> Type of the field elements.
*
* @see <a href="http://mathworld.wolfram.com/LUDecomposition.html">MathWorld</a>
* @see <a href="http://en.wikipedia.org/wiki/LU_decomposition">Wikipedia</a>
*
* @since 4.0
*/
public final class FieldLUDecomposition<T> {
/** Field to which the elements belong. */
private final Field<T> field;
/** Entries of LU decomposition. */
private final FieldDenseMatrix<T> mLU;
/** Pivot permutation associated with LU decomposition. */
private final int[] pivot;
/** Singularity indicator. */
private final boolean isSingular;
/** Parity of the permutation associated with the LU decomposition. */
private final boolean isEven;
/**
* Calculates the LU-decomposition of the given {@code matrix}.
*
* @param matrix Matrix to decompose.
* @throws IllegalArgumentException if the matrix is not square.
*/
private FieldLUDecomposition(FieldDenseMatrix<T> matrix) {
matrix.checkMultiply(matrix);
field = matrix.getField();
final int m = matrix.getRowDimension();
pivot = new int[m];
// Initialize permutation array and parity.
for (int row = 0; row < m; row++) {
pivot[row] = row;
}
mLU = matrix.copy();
boolean even = true;
boolean singular = false;
// Loop over columns.
for (int col = 0; col < m; col++) {
T sum = field.zero();
// Upper.
for (int row = 0; row < col; row++) {
sum = mLU.get(row, col);
for (int i = 0; i < row; i++) {
sum = field.subtract(sum,
field.multiply(mLU.get(row, i),
mLU.get(i, col)));
}
mLU.set(row, col, sum);
}
// Lower.
int nonZero = col; // Permutation row.
for (int row = col; row < m; row++) {
sum = mLU.get(row, col);
for (int i = 0; i < col; i++) {
sum = field.subtract(sum,
field.multiply(mLU.get(row, i),
mLU.get(i, col)));
}
mLU.set(row, col, sum);
if (mLU.get(nonZero, col).equals(field.zero())) {
// try to select a better permutation choice
++nonZero;
}
}
// Singularity check.
if (nonZero >= m) {
singular = true;
} else {
// Pivot if necessary.
if (nonZero != col) {
T tmp = field.zero();
for (int i = 0; i < m; i++) {
tmp = mLU.get(nonZero, i);
mLU.set(nonZero, i, mLU.get(col, i));
mLU.set(col, i, tmp);
}
int temp = pivot[nonZero];
pivot[nonZero] = pivot[col];
pivot[col] = temp;
even = !even;
}
// Divide the lower elements by the "winning" diagonal element.
final T luDiag = mLU.get(col, col);
for (int row = col + 1; row < m; row++) {
mLU.set(row, col, field.divide(mLU.get(row, col),
luDiag));
}
}
}
isSingular = singular;
isEven = even;
}
/**
* Factory method.
*
* @param <T> Type of the field elements.
* @param m Matrix to decompose.
* @return a new instance.
*/
public static <T> FieldLUDecomposition<T> of(FieldDenseMatrix<T> m) {
return new FieldLUDecomposition<>(m);
}
/**
* @return {@code true} if the matrix is singular.
*/
public boolean isSingular() {
return isSingular;
}
/**
* Builds the "L" matrix of the decomposition.
*
* @return the lower triangular matrix.
* @throws SingularMatrixException if the matrix is singular.
*/
public FieldDenseMatrix<T> getL() {
if (isSingular) {
throw new SingularMatrixException();
}
final int m = pivot.length;
final FieldDenseMatrix<T> mL = FieldDenseMatrix.zero(field, m, m);
for (int i = 0; i < m; i++) {
for (int j = 0; j < i; j++) {
mL.set(i, j, mLU.get(i, j));
}
mL.set(i, i, field.one());
}
return mL;
}
/**
* Builds the "U" matrix of the decomposition.
*
* @return the upper triangular matrix.
* @throws SingularMatrixException if the matrix is singular.
*/
public FieldDenseMatrix<T> getU() {
if (isSingular) {
throw new SingularMatrixException();
}
final int m = pivot.length;
final FieldDenseMatrix<T> mU = FieldDenseMatrix.zero(field, m, m);
for (int i = 0; i < m; i++) {
for (int j = i; j < m; j++) {
mU.set(i, j, mLU.get(i, j));
}
}
return mU;
}
/**
* Builds the "P" matrix.
*
* <p>P is a matrix with exactly one element set to {@link Field#one() one} in
* each row and each column, all other elements being set to {@link Field#zero() zero}.
* The positions of the "one" elements are given by the {@link #getPivot()
* pivot permutation vector}.</p>
* @return the "P" rows permutation matrix.
* @throws SingularMatrixException if the matrix is singular.
*
* @see #getPivot()
*/
public FieldDenseMatrix<T> getP() {
if (isSingular) {
throw new SingularMatrixException();
}
final int m = pivot.length;
final FieldDenseMatrix<T> mP = FieldDenseMatrix.zero(field, m, m);
for (int i = 0; i < m; i++) {
mP.set(i, pivot[i], field.one());
}
return mP;
}
/**
* Gets the pivot permutation vector.
*
* @return the pivot permutation vector.
*
* @see #getP()
*/
public int[] getPivot() {
return pivot.clone();
}
/**
* Return the determinant of the matrix.
* @return determinant of the matrix
*/
public T getDeterminant() {
if (isSingular) {
return field.zero();
} else {
final int m = pivot.length;
T determinant = isEven ?
field.one() :
field.negate(field.one());
for (int i = 0; i < m; i++) {
determinant = field.multiply(determinant,
mLU.get(i, i));
}
return determinant;
}
}
/**
* Creates a solver for finding the solution {@code X} of the linear
* system of equations {@code A X = B}.
*
* @return a solver.
* @throws SingularMatrixException if the matrix is singular.
*/
public FieldDecompositionSolver<T> getSolver() {
if (isSingular) {
throw new SingularMatrixException();
}
return new Solver<>(mLU, pivot);
}
/**
* Specialized solver.
*
* @param <T> Type of the field elements.
*/
private static final class Solver<T> implements FieldDecompositionSolver<T> {
/** Field to which the elements belong. */
private final Field<T> field;
/** LU decomposition. */
private final FieldDenseMatrix<T> mLU;
/** Pivot permutation associated with LU decomposition. */
private final int[] pivot;
/**
* Builds a solver from a LU-decomposed matrix.
*
* @param mLU LU matrix.
* @param pivot Pivot permutation associated with the decomposition.
*/
private Solver(final FieldDenseMatrix<T> mLU,
final int[] pivot) {
field = mLU.getField();
this.mLU = mLU.copy();
this.pivot = pivot.clone();
}
/** {@inheritDoc} */
@Override
public FieldDenseMatrix<T> solve(final FieldDenseMatrix<T> b) {
mLU.checkMultiply(b);
final FieldDenseMatrix<T> bp = b.copy();
final int nColB = b.getColumnDimension();
final int m = pivot.length;
// Apply permutations.
for (int row = 0; row < m; row++) {
final int pRow = pivot[row];
for (int col = 0; col < nColB; col++) {
bp.set(row, col,
b.get(row, col));
}
}
// Solve LY = b
for (int col = 0; col < m; col++) {
for (int i = col + 1; i < m; i++) {
for (int j = 0; j < nColB; j++) {
bp.set(i, j,
field.subtract(bp.get(i, j),
field.multiply(bp.get(col, j),
mLU.get(i, col))));
}
}
}
// Solve UX = Y
for (int col = m - 1; col >= 0; col--) {
for (int j = 0; j < nColB; j++) {
bp.set(col, j,
field.divide(bp.get(col, j),
mLU.get(col, col)));
}
for (int i = 0; i < col; i++) {
for (int j = 0; j < nColB; j++) {
bp.set(i, j,
field.subtract(bp.get(i, j),
field.multiply(bp.get(col, j),
mLU.get(i, col))));
}
}
}
return bp;
}
/** {@inheritDoc} */
@Override
public FieldDenseMatrix<T> getInverse() {
return solve(FieldDenseMatrix.identity(field, pivot.length));
}
}
}