Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

MATH-1500: LU decomposition for a matrix whose entries are elements of a field.

This commit is contained in:
Gilles Sadowski 2019-10-27 15:10:11 +01:00
parent 407f1e780e
commit ec77f54bef
3 changed files with 687 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

59
src/main/java/org/apache/commons/math4/field/linalg/FieldDecompositionSolver.java Normal file
View File

@ -0,0 +1,59 @@
/*
* 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.field.linalg;
import java.util.List;
/**
* Interface handling decomposition algorithms that can solve {@code A X = B}.
*
* <p>Decomposition algorithms decompose an A matrix has a product of several specific
* matrices from which they can solve the above system of equations in a least-squares
* sense: Find X such that {@code ||A X - B||} is minimal.</p>
*
* <p>Some solvers like {@link FieldLUDecomposition} can only find the solution for
* square matrices and when the solution is an exact linear solution, i.e. when
* {@code ||A X - B||} is exactly 0.
* Other solvers can also find solutions with non-square matrix {@code A} and with
* non-zero minimal norm.
* If an exact linear solution exists it is also the minimal norm solution.</p>
*
* @param <T> Type of the field elements.
*
* @since 4.0
*/
public interface FieldDecompositionSolver<T> {
/**
* Solves the linear equation {@code A X = B}.
*
* <p>Matrix {@code A} is implicit: It is provided by the underlying
* decomposition algorithm.</p>
*
* @param b Right-hand side of the equation.
* @return the matrix {@code X} that minimizes {@code ||A X - B||}.
* @throws IllegalArgumentException if the dimensions do not match.
*/
FieldDenseMatrix<T> solve(final FieldDenseMatrix<T> b);
/**
* Computes the inverse of a decomposed (square) matrix.
*
* @return the inverse matrix.
*/
FieldDenseMatrix<T> getInverse();
}

361
src/main/java/org/apache/commons/math4/field/linalg/FieldLUDecomposition.java Normal file
View File

@ -0,0 +1,361 @@
/*
* 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.field.linalg;
import java.util.List;
import java.util.ArrayList;
import org.apache.commons.numbers.field.Field;
import org.apache.commons.math4.linear.SingularMatrixException;
import org.apache.commons.math4.exception.DimensionMismatchException;
/**
* 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&times;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 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 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 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));
}
}
}

267
src/test/java/org/apache/commons/math4/field/linalg/FieldLUDecompositionTest.java Normal file
View File

@ -0,0 +1,267 @@
/*
* 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.field.linalg;
import org.apache.commons.numbers.fraction.Fraction;
import org.apache.commons.numbers.field.FractionField;
import org.junit.Test;
import org.junit.Assert;
import org.apache.commons.math4.linear.NonSquareMatrixException;
import org.apache.commons.math4.linear.SingularMatrixException;
public class FieldLUDecompositionTest {
private final Fraction[][] testData = {
{ Fraction.of(1), Fraction.of(2), Fraction.of(3)},
{ Fraction.of(2), Fraction.of(5), Fraction.of(3)},
{ Fraction.of(1), Fraction.of(0), Fraction.of(8)}
};
private final Fraction[][] testDataMinus = {
{ Fraction.of(-1), Fraction.of(-2), Fraction.of(-3)},
{ Fraction.of(-2), Fraction.of(-5), Fraction.of(-3)},
{ Fraction.of(-1), Fraction.of(0), Fraction.of(-8)}
};
private final Fraction[][] luData = {
{ Fraction.of(2), Fraction.of(3), Fraction.of(3) },
{ Fraction.of(2), Fraction.of(3), Fraction.of(7) },
{ Fraction.of(6), Fraction.of(6), Fraction.of(8) }
};
// singular matrices
private final Fraction[][] singular = {
{ Fraction.of(2), Fraction.of(3) },
{ Fraction.of(2), Fraction.of(3) }
};
private final Fraction[][] bigSingular = {
{ Fraction.of(1), Fraction.of(2), Fraction.of(3), Fraction.of(4) },
{ Fraction.of(2), Fraction.of(5), Fraction.of(3), Fraction.of(4) },
{ Fraction.of(7), Fraction.of(3), Fraction.of(256), Fraction.of(1930) },
{ Fraction.of(3), Fraction.of(7), Fraction.of(6), Fraction.of(8) }
}; // 4th row = 1st + 2nd
/**
* @param data Matrix.
* @return a {@link FieldDenseMatrix} instance.
*/
private FieldDenseMatrix<Fraction> create(Fraction[][] data) {
final FieldDenseMatrix<Fraction> m = FieldDenseMatrix.create(FractionField.get(),
data.length,
data[0].length);
for (int i = 0; i < data.length; i++) {
for (int j = 0; j < data.length; j++) {
m.set(i, j, data[i][j]);
}
}
return m;
}
/** test dimensions */
@Test
public void testDimensions() {
FieldDenseMatrix<Fraction> matrix = create(testData);
FieldLUDecomposition<Fraction> LU = FieldLUDecomposition.of(matrix);
Assert.assertEquals(testData.length, LU.getL().getRowDimension());
Assert.assertEquals(testData.length, LU.getU().getRowDimension());
Assert.assertEquals(testData.length, LU.getP().getRowDimension());
}
/** test PA = LU */
@Test
public void testPAEqualLU() {
FieldDenseMatrix<Fraction> matrix = create(testData);
FieldLUDecomposition<Fraction> lu = FieldLUDecomposition.of(matrix);
FieldDenseMatrix<Fraction> l = lu.getL();
FieldDenseMatrix<Fraction> u = lu.getU();
FieldDenseMatrix<Fraction> p = lu.getP();
Assert.assertEquals(p.multiply(matrix), l.multiply(u));
matrix = create(testDataMinus);
lu = FieldLUDecomposition.of(matrix);
l = lu.getL();
u = lu.getU();
p = lu.getP();
Assert.assertEquals(p.multiply(matrix), l.multiply(u));
matrix = FieldDenseMatrix.identity(FractionField.get(), 17);
lu = FieldLUDecomposition.of(matrix);
l = lu.getL();
u = lu.getU();
p = lu.getP();
Assert.assertEquals(p.multiply(matrix), l.multiply(u));
}
/** test that L is lower triangular with unit diagonal */
@Test
public void testLLowerTriangular() {
FieldDenseMatrix<Fraction> matrix = create(testData);
FieldDenseMatrix<Fraction> l = FieldLUDecomposition.of(matrix).getL();
for (int i = 0; i < l.getRowDimension(); i++) {
Assert.assertEquals(Fraction.ONE, l.get(i, i));
for (int j = i + 1; j < l.getColumnDimension(); j++) {
Assert.assertEquals(Fraction.ZERO, l.get(i, j));
}
}
}
/** test that U is upper triangular */
@Test
public void testUUpperTriangular() {
FieldDenseMatrix<Fraction> matrix = create(testData);
FieldDenseMatrix<Fraction> u = FieldLUDecomposition.of(matrix).getU();
for (int i = 0; i < u.getRowDimension(); i++) {
for (int j = 0; j < i; j++) {
Assert.assertEquals(Fraction.ZERO, u.get(i, j));
}
}
}
/** test that P is a permutation matrix */
@Test
public void testPPermutation() {
FieldDenseMatrix<Fraction> matrix = create(testData);
FieldDenseMatrix<Fraction> p = FieldLUDecomposition.of(matrix).getP();
FieldDenseMatrix<Fraction> ppT = p.multiply(p.transpose());
FieldDenseMatrix<Fraction> id = FieldDenseMatrix.identity(FractionField.get(),
p.getRowDimension());
Assert.assertEquals(id, ppT);
for (int i = 0; i < p.getRowDimension(); i++) {
int zeroCount = 0;
int oneCount = 0;
int otherCount = 0;
for (int j = 0; j < p.getColumnDimension(); j++) {
final Fraction e = p.get(i, j);
if (e.equals(Fraction.ZERO)) {
++zeroCount;
} else if (e.equals(Fraction.ONE)) {
++oneCount;
} else {
++otherCount;
}
}
Assert.assertEquals(p.getRowDimension() - 1, zeroCount);
Assert.assertEquals(1, oneCount);
Assert.assertEquals(0, otherCount);
}
for (int j = 0; j < p.getRowDimension(); j++) {
int zeroCount = 0;
int oneCount = 0;
int otherCount = 0;
for (int i = 0; i < p.getColumnDimension(); i++) {
final Fraction e = p.get(i, j);
if (e.equals(Fraction.ZERO)) {
++zeroCount;
} else if (e.equals(Fraction.ONE)) {
++oneCount;
} else {
++otherCount;
}
}
Assert.assertEquals(p.getRowDimension() - 1, zeroCount);
Assert.assertEquals(1, oneCount);
Assert.assertEquals(0, otherCount);
}
}
@Test
public void testIsSingular1() {
FieldLUDecomposition<Fraction> lu = FieldLUDecomposition.of(create(testData));
Assert.assertFalse(lu.isSingular());
lu.getSolver();
}
@Test(expected=SingularMatrixException.class)
public void testIsSingular2() {
FieldLUDecomposition<Fraction> lu = FieldLUDecomposition.of(create(singular));
Assert.assertTrue(lu.isSingular());
lu.getSolver();
}
@Test(expected=SingularMatrixException.class)
public void testIsSingular3() {
FieldLUDecomposition<Fraction> lu = FieldLUDecomposition.of(create(bigSingular));
Assert.assertTrue(lu.isSingular());
lu.getSolver();
}
@Test
public void testMatricesValues1() {
FieldLUDecomposition<Fraction> lu = FieldLUDecomposition.of(create(testData));
FieldDenseMatrix<Fraction> lRef = create(new Fraction[][] {
{ Fraction.of(1), Fraction.of(0), Fraction.of(0) },
{ Fraction.of(2), Fraction.of(1), Fraction.of(0) },
{ Fraction.of(1), Fraction.of(-2), Fraction.of(1) }
});
FieldDenseMatrix<Fraction> uRef = create(new Fraction[][] {
{ Fraction.of(1), Fraction.of(2), Fraction.of(3) },
{ Fraction.of(0), Fraction.of(1), Fraction.of(-3) },
{ Fraction.of(0), Fraction.of(0), Fraction.of(-1) }
});
FieldDenseMatrix<Fraction> pRef = create(new Fraction[][] {
{ Fraction.of(1), Fraction.of(0), Fraction.of(0) },
{ Fraction.of(0), Fraction.of(1), Fraction.of(0) },
{ Fraction.of(0), Fraction.of(0), Fraction.of(1) }
});
int[] pivotRef = { 0, 1, 2 };
// check values against known references
FieldDenseMatrix<Fraction> l = lu.getL();
Assert.assertEquals(lRef, l);
FieldDenseMatrix<Fraction> u = lu.getU();
Assert.assertEquals(uRef, u);
FieldDenseMatrix<Fraction> p = lu.getP();
Assert.assertEquals(pRef, p);
int[] pivot = lu.getPivot();
for (int i = 0; i < pivotRef.length; ++i) {
Assert.assertEquals(pivotRef[i], pivot[i]);
}
}
@Test
public void testMatricesValues2() {
FieldLUDecomposition<Fraction> lu = FieldLUDecomposition.of(create(luData));
FieldDenseMatrix<Fraction> lRef = create(new Fraction[][] {
{ Fraction.of(1), Fraction.of(0), Fraction.of(0) },
{ Fraction.of(3), Fraction.of(1), Fraction.of(0) },
{ Fraction.of(1), Fraction.of(0), Fraction.of(1) }
});
FieldDenseMatrix<Fraction> uRef = create(new Fraction[][] {
{ Fraction.of(2), Fraction.of(3), Fraction.of(3) },
{ Fraction.of(0), Fraction.of(-3), Fraction.of(-1) },
{ Fraction.of(0), Fraction.of(0), Fraction.of(4) }
});
FieldDenseMatrix<Fraction> pRef = create(new Fraction[][] {
{ Fraction.of(1), Fraction.of(0), Fraction.of(0) },
{ Fraction.of(0), Fraction.of(0), Fraction.of(1) },
{ Fraction.of(0), Fraction.of(1), Fraction.of(0) }
});
int[] pivotRef = { 0, 2, 1 };
// check values against known references
FieldDenseMatrix<Fraction> l = lu.getL();
Assert.assertEquals(lRef, l);
FieldDenseMatrix<Fraction> u = lu.getU();
Assert.assertEquals(uRef, u);
FieldDenseMatrix<Fraction> p = lu.getP();
Assert.assertEquals(pRef, p);
int[] pivot = lu.getPivot();
for (int i = 0; i < pivotRef.length; ++i) {
Assert.assertEquals(pivotRef[i], pivot[i]);
}
}
}