From bb8a14e8153064873953abf20acd5f3f2f33c7f1 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Fri, 4 Apr 2014 14:38:17 +0200 Subject: [PATCH] LA: positive integer matrix power #205 --- .../LinearAlgebra/Matrix.Arithmetic.cs | 113 ++++++++++++++++++ .../Double/DenseMatrixTests.cs | 23 ++++ 2 files changed, 136 insertions(+) diff --git a/src/Numerics/LinearAlgebra/Matrix.Arithmetic.cs b/src/Numerics/LinearAlgebra/Matrix.Arithmetic.cs index 5ef604b2..1eeacb99 100644 --- a/src/Numerics/LinearAlgebra/Matrix.Arithmetic.cs +++ b/src/Numerics/LinearAlgebra/Matrix.Arithmetic.cs @@ -1014,6 +1014,119 @@ namespace MathNet.Numerics.LinearAlgebra return result; } + private static Matrix IntPower(int exponent, Matrix x, Matrix y, Matrix work) + { + // We try to be smart about not allocating more matrices than needed + // and to minimize the number of multiplications (not optimal on either though) + + // TODO: For large or non-integer exponents we could diagonalize the matrix with + // a similarity transform (eigenvalue decomposition) + + // return y*x + if (exponent == 1) + { + // return x + if (y == null) + { + return x; + } + + if (work == null) work = y.Multiply(x); else y.Multiply(x, work); + return work; + } + + // return y*x^2 + if (exponent == 2) + { + if (work == null) work = x.Multiply(x); else x.Multiply(x, work); + + // return x^2 + if (y == null) + { + return work; + } + + y.Multiply(work, x); + return x; + } + + // recursive n <-- n/2, y <-- y, x <-- x^2 + if (exponent.IsEven()) + { + // we store the new x in work, keep the y as is and reuse the old x as new work matrix. + if (work == null) work = x.Multiply(x); else x.Multiply(x, work); + return IntPower(exponent/2, work, y, x); + } + + // recursive n <-- (n-1)/2, y <-- x, x <-- x^2 + if (y == null) + { + // we store the new x in work, directly use the old x as y. no work matrix. + if (work == null) work = x.Multiply(x); else x.Multiply(x, work); + return IntPower((exponent - 1)/2, work, x, null); + } + + // recursive n <-- (n-1)/2, y <-- y*x, x <-- x^2 + // we store the new y in work, the new x in y, and reuse the old x as work + if (work == null) work = y.Multiply(x); else y.Multiply(x, work); + x.Multiply(x, y); + return IntPower((exponent - 1)/2, y, work, x); + } + + /// + /// Raises this square matrix to a positive integer exponent and places the results into the result matrix. + /// + /// The positive integer exponent to raise the matrix to. + /// The result of the power. + public void Power(int exponent, Matrix result) + { + if (RowCount != ColumnCount || result.RowCount != RowCount || result.ColumnCount != ColumnCount) + { + throw DimensionsDontMatch(this, result); + } + if (exponent < 0) + { + throw new ArgumentException(Resources.ArgumentNotNegative); + } + if (exponent == 0) + { + Build.DiagonalIdentity(RowCount, ColumnCount).CopyTo(result); + return; + } + if (exponent == 1) + { + CopyTo(result); + return; + } + if (exponent == 2) + { + Multiply(this, result); + return; + } + + var res = IntPower(exponent, Clone(), null, result); + if (!ReferenceEquals(res, result)) + { + res.CopyTo(result); + } + } + + /// + /// Multiplies this square matrix with another matrix and returns the result. + /// + /// The positive integer exponent to raise the matrix to. + public Matrix Power(int exponent) + { + if (RowCount != ColumnCount) throw new ArgumentException(Resources.ArgumentMatrixSquare); + if (exponent < 0) throw new ArgumentException(Resources.ArgumentNotNegative); + + if (exponent == 0) return Build.DiagonalIdentity(RowCount, ColumnCount); + if (exponent == 1) return this; + if (exponent == 2) return Multiply(this); + + return IntPower(exponent, Clone(), null, null); + } + /// /// Negate each element of this matrix. /// diff --git a/src/UnitTests/LinearAlgebraTests/Double/DenseMatrixTests.cs b/src/UnitTests/LinearAlgebraTests/Double/DenseMatrixTests.cs index 02289517..fdb452ca 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/DenseMatrixTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/DenseMatrixTests.cs @@ -173,5 +173,28 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double { Assert.That(() => Matrix.Build.DenseIdentity(order), Throws.TypeOf()); } + + [Test] + public void MatrixPower() + { + var d = Matrix.Build.Random(3, 3, 1); + var d2 = d*d; + var d3 = d2*d; + var d4 = d3*d; + var d5 = d4*d; + var d6 = d5*d; + var d7 = d6*d; + var d8 = d7*d; + + AssertHelpers.AlmostEqual(Matrix.Build.DiagonalIdentity(3), d.Power(0), 10); + AssertHelpers.AlmostEqual(d, d.Power(1), 10); + AssertHelpers.AlmostEqual(d2, d.Power(2), 10); + AssertHelpers.AlmostEqual(d3, d.Power(3), 10); + AssertHelpers.AlmostEqual(d4, d.Power(4), 10); + AssertHelpers.AlmostEqual(d5, d.Power(5), 10); + AssertHelpers.AlmostEqual(d6, d.Power(6), 10); + AssertHelpers.AlmostEqual(d7, d.Power(7), 10); + AssertHelpers.AlmostEqual(d8, d.Power(8), 10); + } } }