# [Numpy-svn] r4968 - in trunk/numpy: core linalg linalg/tests

numpy-svn@scip... numpy-svn@scip...
Sat Apr 5 21:34:35 CDT 2008

```Author: stefan
Date: 2008-04-05 21:34:19 -0500 (Sat, 05 Apr 2008)
New Revision: 4968

Modified:
trunk/numpy/core/defmatrix.py
trunk/numpy/linalg/info.py
trunk/numpy/linalg/linalg.py
trunk/numpy/linalg/tests/test_linalg.py
Log:
Factor out matrix_multiply from defmatrix.  Based on a patch by
Anne Archibald.

Modified: trunk/numpy/core/defmatrix.py
===================================================================
--- trunk/numpy/core/defmatrix.py	2008-04-06 00:40:20 UTC (rev 4967)
+++ trunk/numpy/core/defmatrix.py	2008-04-06 02:34:19 UTC (rev 4968)
@@ -2,7 +2,8 @@

import sys
import numeric as N
-from numeric import concatenate, isscalar, binary_repr
+from numeric import concatenate, isscalar, binary_repr, identity
+from numpy.lib.utils import issubdtype

# make translation table
_table = [None]*256
@@ -46,8 +47,83 @@
"""
return matrix(data, dtype=dtype, copy=False)

+def matrix_power(M,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 type as M is returned. If n<0, the inverse is computed
+    and raised to the exponent.

+    Parameters
+    ----------
+    M : array-like
+        Must be a square array (that is, of dimension two and with
+        equal sizes).
+    n : integer
+        The exponent can be any integer or long integer, positive
+        negative or zero.
+
+    Returns
+    -------
+    M to the power n
+        The return value is a an array the same shape and size as M;
+        if the exponent was positive or zero then the type of the
+        elements is the same as those of M. If the exponent was negative
+        the elements are floating-point.
+
+    Raises
+    ------
+    LinAlgException
+        If the matrix is not numerically invertible, an exception is raised.
+
+    --------
+    The matrix() class provides an equivalent function as the exponentiation
+    operator.
+
+    Examples
+    --------
+    >>> matrix_power(array([[0,1],[-1,0]]),10)
+    array([[-1,  0],
+           [ 0, -1]])
+    """
+    if len(M.shape) != 2 or M.shape[0] != M.shape[1]:
+        raise ValueError("input must be a square array")
+    if not issubdtype(type(n),int):
+        raise TypeError("exponent must be an integer")
+
+    from numpy.linalg import inv
+
+    if n==0:
+        M = M.copy()
+        M[:] = identity(M.shape[0])
+        return M
+    elif n<0:
+        M = inv(M)
+        n *= -1
+
+    result = M
+    if n <= 3:
+        for _ in range(n-1):
+            result=N.dot(result,M)
+        return result
+
+    # binary decomposition to reduce the number of Matrix
+    # multiplications for n > 3.
+    beta = binary_repr(n)
+    Z,q,t = M,0,len(beta)
+    while beta[t-q-1] == '0':
+        Z = N.dot(Z,Z)
+        q += 1
+    result = Z
+    for k in range(q+1,t):
+        Z = N.dot(Z,Z)
+        if beta[t-k-1] == '1':
+            result = N.dot(result,Z)
+    return result
+
+
class matrix(N.ndarray):
"""mat = matrix(data, dtype=None, copy=True)

@@ -195,39 +271,7 @@
return self

def __pow__(self, other):
-        shape = self.shape
-        if len(shape) != 2 or shape[0] != shape[1]:
-            raise TypeError, "matrix is not square"
-        if type(other) in (type(1), type(1L)):
-            if other==0:
-                return matrix(N.identity(shape[0]))
-            if other<0:
-                x = self.I
-                other=-other
-            else:
-                x=self
-            result = x
-            if other <= 3:
-                while(other>1):
-                    result=result*x
-                    other=other-1
-                return result
-            # binary decomposition to reduce the number of Matrix
-            #  Multiplies for other > 3.
-            beta = binary_repr(other)
-            t = len(beta)
-            Z,q = x.copy(),0
-            while beta[t-q-1] == '0':
-                Z *= Z
-                q += 1
-            result = Z.copy()
-            for k in range(q+1,t):
-                Z *= Z
-                if beta[t-k-1] == '1':
-                    result *= Z
-            return result
-        else:
-            raise TypeError, "exponent must be an integer"
+        return matrix_power(self, other)

def __rpow__(self, other):
return NotImplemented

Modified: trunk/numpy/linalg/info.py
===================================================================
--- trunk/numpy/linalg/info.py	2008-04-06 00:40:20 UTC (rev 4967)
+++ trunk/numpy/linalg/info.py	2008-04-06 02:34:19 UTC (rev 4968)
@@ -10,6 +10,7 @@
- lstsq           Solve linear least-squares problem
- pinv            Pseudo-inverse (Moore-Penrose) calculated using a singular
value decomposition
+- matrix_power    Integer power of a square matrix

Eigenvalues and decompositions:

Modified: trunk/numpy/linalg/linalg.py
===================================================================
--- trunk/numpy/linalg/linalg.py	2008-04-06 00:40:20 UTC (rev 4967)
+++ trunk/numpy/linalg/linalg.py	2008-04-06 02:34:19 UTC (rev 4968)
@@ -9,7 +9,7 @@
zgetrf, dpotrf, zpotrf, dgeqrf, zgeqrf, zungqr, dorgqr.
"""

-__all__ = ['solve', 'tensorsolve', 'tensorinv',
+__all__ = ['matrix_power', 'solve', 'tensorsolve', 'tensorinv',
'inv', 'cholesky',
'eigvals',
'eigvalsh', 'pinv',
@@ -26,6 +26,7 @@
isfinite, size
from numpy.lib import triu
from numpy.linalg import lapack_lite
+from numpy.core.defmatrix import matrix_power

fortran_int = intc

@@ -134,6 +135,7 @@
if size(a) == 0:
raise LinAlgError("Arrays cannot be empty")

+
# Linear equations

def tensorsolve(a, b, axes=None):
@@ -326,6 +328,7 @@
a, wrap = _makearray(a)
return wrap(solve(a, identity(a.shape[0], dtype=a.dtype)))

+
# Cholesky decomposition

def cholesky(a):
@@ -1053,6 +1056,7 @@
sign = add.reduce(pivots != arange(1, n+1)) % 2
return (1.-2.*sign)*multiply.reduce(diagonal(a), axis=-1)

+
# Linear Least Squares

def lstsq(a, b, rcond=-1):

Modified: trunk/numpy/linalg/tests/test_linalg.py
===================================================================
--- trunk/numpy/linalg/tests/test_linalg.py	2008-04-06 00:40:20 UTC (rev 4967)
+++ trunk/numpy/linalg/tests/test_linalg.py	2008-04-06 02:34:19 UTC (rev 4968)
@@ -6,6 +6,7 @@
from numpy import array, single, double, csingle, cdouble, dot, identity, \
multiply, atleast_2d
from numpy import linalg
+from linalg import matrix_power
restore_path()

old_assert_almost_equal = assert_almost_equal
@@ -46,38 +47,38 @@
except linalg.LinAlgError, e:
pass

-class test_solve(LinalgTestCase):
+class TestSolve(LinalgTestCase):
def do(self, a, b):
x = linalg.solve(a, b)
assert_almost_equal(b, dot(a, x))

-class test_inv(LinalgTestCase):
+class TestInv(LinalgTestCase):
def do(self, a, b):
a_inv = linalg.inv(a)
assert_almost_equal(dot(a, a_inv), identity(a.shape[0]))

-class test_eigvals(LinalgTestCase):
+class TestEigvals(LinalgTestCase):
def do(self, a, b):
ev = linalg.eigvals(a)
evalues, evectors = linalg.eig(a)
assert_almost_equal(ev, evalues)

-class test_eig(LinalgTestCase):
+class TestEig(LinalgTestCase):
def do(self, a, b):
evalues, evectors = linalg.eig(a)
assert_almost_equal(dot(a, evectors), evectors*evalues)

-class test_svd(LinalgTestCase):
+class TestSVD(LinalgTestCase):
def do(self, a, b):
u, s, vt = linalg.svd(a, 0)
assert_almost_equal(a, dot(u*s, vt))

-class test_pinv(LinalgTestCase):
+class TestPinv(LinalgTestCase):
def do(self, a, b):
a_ginv = linalg.pinv(a)
assert_almost_equal(dot(a, a_ginv), identity(a.shape[0]))

-class test_det(LinalgTestCase):
+class TestDet(LinalgTestCase):
def do(self, a, b):
d = linalg.det(a)
if a.dtype.type in (single, double):
@@ -87,7 +88,7 @@
assert_almost_equal(d, multiply.reduce(ev))

-class test_lstsq(LinalgTestCase):
+class TestLstsq(LinalgTestCase):
def do(self, a, b):
u, s, vt = linalg.svd(a, 0)
x, residuals, rank, sv = linalg.lstsq(a, b)
@@ -95,5 +96,58 @@
assert_equal(rank, a.shape[0])
assert_almost_equal(sv, s)

+class TestMatrixPower(ParametricTestCase):
+    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]])
+
+    large = identity(10)
+    t = large[1,:].copy()
+    large[1,:] = large[0,:]
+    large[0,:] = t
+
+    def test_large_power(self):
+        assert_equal(matrix_power(self.R90,2L**100+2**10+2**5+1),self.R90)
+    def test_large_power_trailing_zero(self):
+        assert_equal(matrix_power(self.R90,2L**100+2**10+2**5),identity(2))
+
+    def testip_zero(self):
+        def tz(M):
+            mz = matrix_power(M,0)
+            assert_equal(mz, identity(M.shape[0]))
+            assert_equal(mz.dtype, M.dtype)
+        for M in [self.Arb22, self.arbfloat, self.large]:
+            yield tz, M
+
+    def testip_one(self):
+        def tz(M):
+            mz = matrix_power(M,1)
+            assert_equal(mz, M)
+            assert_equal(mz.dtype, M.dtype)
+        for M in [self.Arb22, self.arbfloat, self.large]:
+            yield tz, M
+
+    def testip_two(self):
+        def tz(M):
+            mz = matrix_power(M,2)
+            assert_equal(mz, dot(M,M))
+            assert_equal(mz.dtype, M.dtype)
+        for M in [self.Arb22, self.arbfloat, self.large]:
+            yield tz, M
+
+    def testip_invert(self):
+        def tz(M):
+            mz = matrix_power(M,-1)
+            assert_almost_equal(identity(M.shape[0]), dot(mz,M))
+        for M in [self.R90, self.Arb22, self.arbfloat, self.large]:
+            yield tz, M
+
+    def test_invert_noninvertible(self):
+        import numpy.linalg
+        self.assertRaises(numpy.linalg.linalg.LinAlgError,
+                lambda: matrix_power(self.noninv,-1))
+
+
if __name__ == '__main__':
NumpyTest().run()

```