From 31fa9a7faa27fc25f12a186e4b41de0bf823d7d9 Mon Sep 17 00:00:00 2001 From: Mikkel Meyer Andersen Date: Mon, 21 Mar 2011 08:14:17 +0000 Subject: [PATCH] MATH-435 Added "power" method in "RealMatrix" interface and default implementation in "AbstractRealMatrix". git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1083698 13f79535-47bb-0310-9956-ffa450edef68 --- .../math/linear/AbstractRealMatrix.java | 62 +++++++++++++++++++ .../commons/math/linear/RealMatrix.java | 10 +++ .../math/linear/Array2DRowRealMatrixTest.java | 54 ++++++++++++++++ 3 files changed, 126 insertions(+) diff --git a/src/main/java/org/apache/commons/math/linear/AbstractRealMatrix.java b/src/main/java/org/apache/commons/math/linear/AbstractRealMatrix.java index bdc296173..3b4b70ec3 100644 --- a/src/main/java/org/apache/commons/math/linear/AbstractRealMatrix.java +++ b/src/main/java/org/apache/commons/math/linear/AbstractRealMatrix.java @@ -17,6 +17,8 @@ package org.apache.commons.math.linear; +import java.util.ArrayList; + import org.apache.commons.math.exception.NoDataException; import org.apache.commons.math.exception.NotStrictlyPositiveException; import org.apache.commons.math.exception.DimensionMismatchException; @@ -147,6 +149,66 @@ public abstract class AbstractRealMatrix implements RealMatrix { /** {@inheritDoc} */ public RealMatrix preMultiply(final RealMatrix m) { return m.multiply(this); + } + + /** {@inheritDoc} */ + @Override + public RealMatrix power(final int p) { + if (p < 0) { + throw new IllegalArgumentException("p must be >= 0"); + } + + if (!isSquare()) { + throw new NonSquareMatrixException(getRowDimension(), getColumnDimension()); + } + + if (p == 0) { + return MatrixUtils.createRealIdentityMatrix(this.getRowDimension()); + } + + if (p == 1) { + return this.copy(); + } + + final int power = p - 1; + + /* + * Only log_2(p) operations is used by doing as follows: + * 5^214 = 5^128 * 5^64 * 5^16 * 5^4 * 5^2 + * + * In general, the same approach is used for A^p. + */ + + final char[] binaryRepresentation = Integer.toBinaryString(power).toCharArray(); + final ArrayList nonZeroPositions = new ArrayList(); + int maxI = -1; + + for (int i = 0; i < binaryRepresentation.length; ++i) { + if (binaryRepresentation[i] == '1') { + final int pos = binaryRepresentation.length - i - 1; + nonZeroPositions.add(pos); + + // The positions are taken in turn, so maxI is only changed once + if (maxI == -1) { + maxI = pos; + } + } + } + + RealMatrix[] results = new RealMatrix[maxI + 1]; + results[0] = this.copy(); + + for (int i = 1; i <= maxI; ++i) { + results[i] = results[i-1].multiply(results[i-1]); + } + + RealMatrix result = this.copy(); + + for (Integer i : nonZeroPositions) { + result = result.multiply(results[i]); + } + + return result; } /** {@inheritDoc} */ diff --git a/src/main/java/org/apache/commons/math/linear/RealMatrix.java b/src/main/java/org/apache/commons/math/linear/RealMatrix.java index c6bd06fc1..5a4f2ac85 100644 --- a/src/main/java/org/apache/commons/math/linear/RealMatrix.java +++ b/src/main/java/org/apache/commons/math/linear/RealMatrix.java @@ -98,6 +98,16 @@ public interface RealMatrix extends AnyMatrix { * if rowDimension(this) != columnDimension(m) */ RealMatrix preMultiply(RealMatrix m); + + /** + * Returns the result multiplying this with itself p times. + * Depending on the underlying storage, instability for high powers might occur. + * @param p raise this to power p + * @return this^p + * @throws IllegalArgumentException if p < 0 + * NonSquareMatrixException if the matrix is not square + */ + RealMatrix power(final int p); /** * Returns matrix entries as a two-dimensional array. diff --git a/src/test/java/org/apache/commons/math/linear/Array2DRowRealMatrixTest.java b/src/test/java/org/apache/commons/math/linear/Array2DRowRealMatrixTest.java index 4e2394acf..e39a7fa37 100644 --- a/src/test/java/org/apache/commons/math/linear/Array2DRowRealMatrixTest.java +++ b/src/test/java/org/apache/commons/math/linear/Array2DRowRealMatrixTest.java @@ -96,6 +96,7 @@ public final class Array2DRowRealMatrixTest { // tolerances protected double entryTolerance = 10E-16; protected double normTolerance = 10E-14; + protected double powerTolerance = 10E-16; /** test dimensions */ @Test @@ -222,6 +223,59 @@ public final class Array2DRowRealMatrixTest { TestUtils.assertEquals("m3*m4=m5", m3.multiply(m4), m5, entryTolerance); } + public void testPower() { + Array2DRowRealMatrix m = new Array2DRowRealMatrix(testData); + Array2DRowRealMatrix mInv = new Array2DRowRealMatrix(testDataInv); + Array2DRowRealMatrix mPlusInv = new Array2DRowRealMatrix(testDataPlusInv); + Array2DRowRealMatrix identity = new Array2DRowRealMatrix(id); + + TestUtils.assertEquals("m^0", m.power(0), + identity, entryTolerance); + TestUtils.assertEquals("mInv^0", mInv.power(0), + identity, entryTolerance); + TestUtils.assertEquals("mPlusInv^0", mPlusInv.power(0), + identity, entryTolerance); + + TestUtils.assertEquals("m^1", m.power(1), + m, entryTolerance); + TestUtils.assertEquals("mInv^1", mInv.power(1), + mInv, entryTolerance); + TestUtils.assertEquals("mPlusInv^1", mPlusInv.power(1), + mPlusInv, entryTolerance); + + RealMatrix C1 = m.copy(); + RealMatrix C2 = mInv.copy(); + RealMatrix C3 = mPlusInv.copy(); + + for (int i = 2; i <= 10; ++i) { + C1 = C1.multiply(m); + C2 = C2.multiply(mInv); + C3 = C3.multiply(mPlusInv); + + TestUtils.assertEquals("m^" + i, m.power(i), + C1, entryTolerance); + TestUtils.assertEquals("mInv^" + i, mInv.power(i), + C2, entryTolerance); + TestUtils.assertEquals("mPlusInv^" + i, mPlusInv.power(i), + C3, entryTolerance); + } + + try { + Array2DRowRealMatrix mNotSquare = new Array2DRowRealMatrix(testData2T); + mNotSquare.power(2); + Assert.fail("Expecting NonSquareMatrixException"); + } catch (NonSquareMatrixException ex) { + // ignored + } + + try { + m.power(-1); + Assert.fail("Expecting IllegalArgumentException"); + } catch (IllegalArgumentException ex) { + // ignored + } + } + /** test trace */ @Test public void testTrace() {