This is an automated email from the ASF dual-hosted git repository. erans pushed a commit to branch master in repository https://gitbox.apache.org/repos/asf/commons-math.git
commit ec77f54bef7e60de920307bd75ce44a2fa104670 Author: Gilles Sadowski <[email protected]> AuthorDate: Sun Oct 27 15:10:11 2019 +0100 MATH-1500: LU decomposition for a matrix whose entries are elements of a field. --- .../field/linalg/FieldDecompositionSolver.java | 59 ++++ .../math4/field/linalg/FieldLUDecomposition.java | 361 +++++++++++++++++++++ .../field/linalg/FieldLUDecompositionTest.java | 267 +++++++++++++++ 3 files changed, 687 insertions(+) diff --git a/src/main/java/org/apache/commons/math4/field/linalg/FieldDecompositionSolver.java b/src/main/java/org/apache/commons/math4/field/linalg/FieldDecompositionSolver.java new file mode 100644 index 0000000..0136062 --- /dev/null +++ b/src/main/java/org/apache/commons/math4/field/linalg/FieldDecompositionSolver.java @@ -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(); +} diff --git a/src/main/java/org/apache/commons/math4/field/linalg/FieldLUDecomposition.java b/src/main/java/org/apache/commons/math4/field/linalg/FieldLUDecomposition.java new file mode 100644 index 0000000..ab53e29 --- /dev/null +++ b/src/main/java/org/apache/commons/math4/field/linalg/FieldLUDecomposition.java @@ -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×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)); + } + } +} diff --git a/src/test/java/org/apache/commons/math4/field/linalg/FieldLUDecompositionTest.java b/src/test/java/org/apache/commons/math4/field/linalg/FieldLUDecompositionTest.java new file mode 100644 index 0000000..b756643 --- /dev/null +++ b/src/test/java/org/apache/commons/math4/field/linalg/FieldLUDecompositionTest.java @@ -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]); + } + } +}
