function for computing powers of a matrix

[LevGivon]

It would be nice to add the following function to numpy to allow for the raising of matricies to arbitrary exponents (and perhaps modify matrix.__pow__() to use this approach when one attempts to raise a matrix to a noninteger exponent):

from numpy import diag,dot,shape,eye
from numpy.linalg import eig,inv

def mpower(x,y):                                                                
    '''Compute x raised to the power y when x is a square matrix and y          
    is a scalar.'''                                                             
                                                                                
    s = shape(x)                                                                
    if len(s) != 2 or s[0] != s[1]:                                
        raise ValueError('matrix must be square')                               
    if y == 0:                                                                  
        return eye(s[0]) 
    [e,v] = eig(x)                                                              
    d = diag(e)                                                                 
    return dot(dot(v,d**y),inv(v)) 

[nils]

Your approach is only valid for non-defective dense matrices.

[LevGivon]

I suppose one could modify the proposed function to check whether the matrix is defective within an appropriate tolerance.

Is the power of a defective matrix mathematically defined?

[charris]

Is the power of a defective matrix mathematically defined?

Arbitrary powers? No, try to find the square root of [[0, 1],[0,0]].

[LevGivon]

I suppose one could modify the proposed function to check whether the matrix is defective within an appropriate tolerance.

For the sake of comparison, Matlab's matrix exponentation mechanism apparently uses decomposition without checking for defectiveness (and encounters numerical issues when the matrix is defective). Octave returns nans and infinities when non-integer exponents are specified even for non-defective matricies.

[stefan]

The approach used in numpy.matlib is a binary decomposition. If the SVD can't be used for some reason, it would be useful to have this exposed to numpy in general.

The equivalent code looks something like:

def mpower(A,n):
    B = A
    p = int(np.floor(np.log(n)/np.log(2))) # 2^? is n
    for i in range(p):
        B = np.dot(B,B)

    for i in range(n-2**p):
        B = np.dot(B,A)

    return B

See defmatrix.py:215.

[charris]

Powers can be defined for normal matrices, i.e., matrices for which A*transpose(A) = transpose(A)*A, but normality itself is subject to numerical fuzziness in computational practice. And for normal matrices you still need to deal with square roots of negative reals and the fact that there can be a multiplicity of roots, and so on and so forth. I think the user needs to know what sort of matrices they are dealing with and what they intend before they use a function like this.

This functionality could be part of a package, perhaps in a package for symmetric positive indefinite matrices; those matrices are like the non-negative reals. But I don't think numpy should try to implement such a specialized function.

However, raising matrices to (small) non-negative integer powers is possible and might cover most needs. A special integer power function might do the trick.

[peridot]

patch to expose defmatrix's patch as a function

[peridot]

The matrix [[1.1].[0.1]] has a square root ([[1.0.5],[0,1]]). It is unfortunately rather difficult to work with this sort of thing numerically. I suggest non-integer powers be left to scipy's expm and friends. Integer powers are easily dealt with.

[stefan]

Fixed in r4968.