diff options
Diffstat (limited to 'numpy/linalg')
| -rw-r--r-- | numpy/linalg/linalg.py | 117 | ||||
| -rw-r--r-- | numpy/linalg/tests/test_linalg.py | 33 |
2 files changed, 131 insertions, 19 deletions
diff --git a/numpy/linalg/linalg.py b/numpy/linalg/linalg.py index 5ee230f92..5757b1827 100644 --- a/numpy/linalg/linalg.py +++ b/numpy/linalg/linalg.py @@ -16,20 +16,20 @@ __all__ = ['matrix_power', 'solve', 'tensorsolve', 'tensorinv', 'inv', 'svd', 'eig', 'eigh', 'lstsq', 'norm', 'qr', 'cond', 'matrix_rank', 'LinAlgError', 'multi_dot'] +import operator import warnings from numpy.core import ( array, asarray, zeros, empty, empty_like, intc, single, double, - csingle, cdouble, inexact, complexfloating, newaxis, ravel, all, Inf, dot, - add, multiply, sqrt, maximum, fastCopyAndTranspose, sum, isfinite, size, - finfo, errstate, geterrobj, longdouble, moveaxis, amin, amax, product, abs, - broadcast, atleast_2d, intp, asanyarray, object_, ones, matmul, - swapaxes, divide, count_nonzero, ndarray, isnan + csingle, cdouble, inexact, complexfloating, newaxis, all, Inf, dot, + add, multiply, sqrt, fastCopyAndTranspose, sum, isfinite, + finfo, errstate, geterrobj, moveaxis, amin, amax, product, abs, + atleast_2d, intp, asanyarray, object_, matmul, + swapaxes, divide, count_nonzero, isnan ) from numpy.core.multiarray import normalize_axis_index -from numpy.lib import triu, asfarray +from numpy.lib.twodim_base import triu, eye from numpy.linalg import lapack_lite, _umath_linalg -from numpy.matrixlib.defmatrix import matrix_power # For Python2/3 compatibility _N = b'N' @@ -532,6 +532,109 @@ def inv(a): return wrap(ainv.astype(result_t, copy=False)) +def matrix_power(a, n): + """ + Raise a square matrix to the (integer) power `n`. + + For positive integers `n`, the power is computed by repeated matrix + squarings and matrix multiplications. If ``n == 0``, the identity matrix + of the same shape as M is returned. If ``n < 0``, the inverse + is computed and then raised to the ``abs(n)``. + + Parameters + ---------- + a : (..., M, M) array_like + Matrix to be "powered." + n : int + The exponent can be any integer or long integer, positive, + negative, or zero. + + Returns + ------- + a**n : (..., M, M) ndarray or matrix object + The return value is the same shape and type as `M`; + if the exponent is positive or zero then the type of the + elements is the same as those of `M`. If the exponent is + negative the elements are floating-point. + + Raises + ------ + LinAlgError + For matrices that are not square or that (for negative powers) cannot + be inverted numerically. + + Examples + -------- + >>> from numpy.linalg import matrix_power + >>> i = np.array([[0, 1], [-1, 0]]) # matrix equiv. of the imaginary unit + >>> matrix_power(i, 3) # should = -i + array([[ 0, -1], + [ 1, 0]]) + >>> matrix_power(i, 0) + array([[1, 0], + [0, 1]]) + >>> matrix_power(i, -3) # should = 1/(-i) = i, but w/ f.p. elements + array([[ 0., 1.], + [-1., 0.]]) + + Somewhat more sophisticated example + + >>> q = np.zeros((4, 4)) + >>> q[0:2, 0:2] = -i + >>> q[2:4, 2:4] = i + >>> q # one of the three quaternion units not equal to 1 + array([[ 0., -1., 0., 0.], + [ 1., 0., 0., 0.], + [ 0., 0., 0., 1.], + [ 0., 0., -1., 0.]]) + >>> matrix_power(q, 2) # = -np.eye(4) + array([[-1., 0., 0., 0.], + [ 0., -1., 0., 0.], + [ 0., 0., -1., 0.], + [ 0., 0., 0., -1.]]) + + """ + a = asanyarray(a) + _assertRankAtLeast2(a) + _assertNdSquareness(a) + + try: + n = operator.index(n) + except TypeError: + raise TypeError("exponent must be an integer") + + if n == 0: + a = empty_like(a) + a[...] = eye(a.shape[-2], dtype=a.dtype) + return a + + elif n < 0: + a = inv(a) + n = abs(n) + + # short-cuts. + if n == 1: + return a + + elif n == 2: + return matmul(a, a) + + elif n == 3: + return matmul(matmul(a, a), a) + + # Use binary decomposition to reduce the number of matrix multiplications. + # Here, we iterate over the bits of n, from LSB to MSB, raise `a` to + # increasing powers of 2, and multiply into the result as needed. + z = result = None + while n > 0: + z = a if z is None else matmul(z, z) + n, bit = divmod(n, 2) + if bit: + result = z if result is None else matmul(result, z) + + return result + + # Cholesky decomposition def cholesky(a): diff --git a/numpy/linalg/tests/test_linalg.py b/numpy/linalg/tests/test_linalg.py index b3dd2e4ae..5ed1ff1c0 100644 --- a/numpy/linalg/tests/test_linalg.py +++ b/numpy/linalg/tests/test_linalg.py @@ -10,8 +10,8 @@ import traceback import pytest import numpy as np -from numpy import array, single, double, csingle, cdouble, dot, identity -from numpy import multiply, atleast_2d, inf, asarray +from numpy import array, single, double, csingle, cdouble, dot, identity, matmul +from numpy import multiply, atleast_2d, inf, asarray, matrix from numpy import linalg from numpy.linalg import matrix_power, norm, matrix_rank, multi_dot, LinAlgError from numpy.linalg.linalg import _multi_dot_matrix_chain_order @@ -445,8 +445,7 @@ def identity_like_generalized(a): a = asarray(a) if a.ndim >= 3: r = np.empty(a.shape, dtype=a.dtype) - for c in itertools.product(*map(range, a.shape[:-2])): - r[c] = identity(a.shape[-2]) + r[...] = identity(a.shape[-2]) return r else: return identity(a.shape[0]) @@ -927,16 +926,21 @@ class TestMatrixPower(object): R90 = array([[0, 1], [-1, 0]]) Arb22 = array([[4, -7], [-2, 10]]) noninv = array([[1, 0], [0, 0]]) - arbfloat = array([[0.1, 3.2], [1.2, 0.7]]) + arbfloat = array([[[0.1, 3.2], [1.2, 0.7]], + [[0.2, 6.4], [2.4, 1.4]]]) large = identity(10) t = large[1, :].copy() - large[1, :] = large[0,:] + large[1, :] = large[0, :] large[0, :] = t def test_large_power(self): assert_equal( matrix_power(self.R90, 2 ** 100 + 2 ** 10 + 2 ** 5 + 1), self.R90) + assert_equal( + matrix_power(self.R90, 2 ** 100 + 2 ** 10 + 1), self.R90) + assert_equal( + matrix_power(self.R90, 2 ** 100 + 2 + 1), -self.R90) def test_large_power_trailing_zero(self): assert_equal( @@ -945,7 +949,7 @@ class TestMatrixPower(object): def testip_zero(self): def tz(M): mz = matrix_power(M, 0) - assert_equal(mz, identity(M.shape[0])) + assert_equal(mz, identity_like_generalized(M)) assert_equal(mz.dtype, M.dtype) for M in [self.Arb22, self.arbfloat, self.large]: tz(M) @@ -961,7 +965,7 @@ class TestMatrixPower(object): def testip_two(self): def tz(M): mz = matrix_power(M, 2) - assert_equal(mz, dot(M, M)) + assert_equal(mz, matmul(M, M)) assert_equal(mz.dtype, M.dtype) for M in [self.Arb22, self.arbfloat, self.large]: tz(M) @@ -969,14 +973,19 @@ class TestMatrixPower(object): def testip_invert(self): def tz(M): mz = matrix_power(M, -1) - assert_almost_equal(identity(M.shape[0]), dot(mz, M)) + assert_almost_equal(matmul(mz, M), identity_like_generalized(M)) for M in [self.R90, self.Arb22, self.arbfloat, self.large]: tz(M) def test_invert_noninvertible(self): - import numpy.linalg - assert_raises(numpy.linalg.linalg.LinAlgError, - lambda: matrix_power(self.noninv, -1)) + assert_raises(LinAlgError, matrix_power, self.noninv, -1) + + def test_invalid(self): + assert_raises(TypeError, matrix_power, self.R90, 1.5) + assert_raises(TypeError, matrix_power, self.R90, [1]) + assert_raises(LinAlgError, matrix_power, np.array([1]), 1) + assert_raises(LinAlgError, matrix_power, np.array([[1], [2]]), 1) + assert_raises(LinAlgError, matrix_power, np.ones((4, 3, 2)), 1) class TestBoolPower(object): |