You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by er...@apache.org on 2012/09/21 17:41:06 UTC
svn commit: r1388555 - in /commons/proper/math/trunk/src: changes/changes.xml
main/java/org/apache/commons/math3/linear/MatrixUtils.java
test/java/org/apache/commons/math3/linear/MatrixUtilsTest.java
Author: erans
Date: Fri Sep 21 15:41:06 2012
New Revision: 1388555
URL: http://svn.apache.org/viewvc?rev=1388555&view=rev
Log:
MATH-860
Matrix "block inversion".
Modified:
commons/proper/math/trunk/src/changes/changes.xml
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/MatrixUtils.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/MatrixUtilsTest.java
Modified: commons/proper/math/trunk/src/changes/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/changes/changes.xml?rev=1388555&r1=1388554&r2=1388555&view=diff
==============================================================================
--- commons/proper/math/trunk/src/changes/changes.xml (original)
+++ commons/proper/math/trunk/src/changes/changes.xml Fri Sep 21 15:41:06 2012
@@ -52,6 +52,9 @@ If the output is not quite correct, chec
<body>
<release version="3.1" date="TBD" description="
">
+ <action dev="erans" type="add" issue="MATH-860">
+ Added matrix "block inversion" (in "o.a.c.m.linear.MatrixUtils").
+ </action>
<action dev="erans" type="fix" issue="MATH-865">
"CMAESOptimizer": Raise an exception when the variables normalization
would fail due to overflow.
Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/MatrixUtils.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/MatrixUtils.java?rev=1388555&r1=1388554&r2=1388555&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/MatrixUtils.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/MatrixUtils.java Fri Sep 21 15:41:06 2012
@@ -929,4 +929,57 @@ public class MatrixUtils {
}
}
}
+
+ /**
+ * Computes the inverse of the given matrix by splitting it into
+ * 4 sub-matrices.
+ *
+ * @param m Matrix whose inverse must be computed.
+ * @param splitIndex Index that determines the "split" line and
+ * column.
+ * The element corresponding to this index will part of the
+ * upper-left sub-matrix.
+ * @return the inverse of {@code m}.
+ * @throws NonSquareMatrixException if {@code m} is not square.
+ */
+ public static RealMatrix blockInverse(RealMatrix m,
+ int splitIndex) {
+ final int n = m.getRowDimension();
+ if (m.getColumnDimension() != n) {
+ throw new NonSquareMatrixException(m.getRowDimension(),
+ m.getColumnDimension());
+ }
+
+ final int splitIndex1 = splitIndex + 1;
+
+ final RealMatrix a = m.getSubMatrix(0, splitIndex, 0, splitIndex);
+ final RealMatrix b = m.getSubMatrix(0, splitIndex, splitIndex1, n - 1);
+ final RealMatrix c = m.getSubMatrix(splitIndex1, n - 1, 0, splitIndex);
+ final RealMatrix d = m.getSubMatrix(splitIndex1, n - 1, splitIndex1, n - 1);
+
+ final SingularValueDecomposition aDec = new SingularValueDecomposition(a);
+ final RealMatrix aInv = aDec.getSolver().getInverse();
+
+ final SingularValueDecomposition dDec = new SingularValueDecomposition(d);
+ final RealMatrix dInv = dDec.getSolver().getInverse();
+
+ final RealMatrix tmp1 = a.subtract(b.multiply(dInv).multiply(c));
+ final SingularValueDecomposition tmp1Dec = new SingularValueDecomposition(tmp1);
+ final RealMatrix result00 = tmp1Dec.getSolver().getInverse();
+
+ final RealMatrix tmp2 = d.subtract(c.multiply(aInv).multiply(b));
+ final SingularValueDecomposition tmp2Dec = new SingularValueDecomposition(tmp2);
+ final RealMatrix result11 = tmp2Dec.getSolver().getInverse();
+
+ final RealMatrix result01 = aInv.multiply(b).multiply(result11).scalarMultiply(-1);
+ final RealMatrix result10 = dInv.multiply(c).multiply(result00).scalarMultiply(-1);
+
+ final RealMatrix result = new Array2DRowRealMatrix(n, n);
+ result.setSubMatrix(result00.getData(), 0, 0);
+ result.setSubMatrix(result01.getData(), 0, splitIndex1);
+ result.setSubMatrix(result10.getData(), splitIndex1, 0);
+ result.setSubMatrix(result11.getData(), splitIndex1, splitIndex1);
+
+ return result;
+ }
}
Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/MatrixUtilsTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/MatrixUtilsTest.java?rev=1388555&r1=1388554&r2=1388555&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/MatrixUtilsTest.java (original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/MatrixUtilsTest.java Fri Sep 21 15:41:06 2012
@@ -325,5 +325,42 @@ public final class MatrixUtilsTest {
MatrixUtils.solveUpperTriangularSystem(rm, b);
TestUtils.assertEquals( new double[]{-1,3,1} , b.toArray() , 1.0e-12);
}
-}
+ /**
+ * This test should probably be replaced by one that could show
+ * whether this algorithm can sometimes perform better (precision- or
+ * performance-wise) than the direct inversion of the whole matrix.
+ */
+ @Test
+ public void testBlockInverse() {
+ final double[][] data = {
+ { -1, 0, 123, 4 },
+ { -56, 78.9, -0.1, -23.4 },
+ { 5.67, 8, -9, 1011 },
+ { 12, 345, -67.8, 9 },
+ };
+
+ final RealMatrix m = new Array2DRowRealMatrix(data);
+ final int len = data.length;
+ final double tol = 1e-14;
+
+ for (int splitIndex = 0; splitIndex < 3; splitIndex++) {
+ final RealMatrix mInv = MatrixUtils.blockInverse(m, splitIndex);
+ final RealMatrix id = m.multiply(mInv);
+
+ // Check that we recovered the identity matrix.
+ for (int i = 0; i < len; i++) {
+ for (int j = 0; j < len; j++) {
+ final double entry = id.getEntry(i, j);
+ if (i == j) {
+ Assert.assertEquals("[" + i + "][" + j + "]",
+ 1, entry, tol);
+ } else {
+ Assert.assertEquals("[" + i + "][" + j + "]",
+ 0, entry, tol);
+ }
+ }
+ }
+ }
+ }
+}