r24156 - /branches/disp_spin_speed/lib/linear_algebra/matrix_exponential.py -- June 19, 2014 - 17:42

Author: tlinnet
Date: Thu Jun 19 17:41:59 2014
New Revision: 24156

URL: http://svn.gna.org/viewcvs/relax?rev=24156&view=rev
Log:
Added function to compute the matrix exponential for higher dimensional data 
of shape [NE][NS][NM][NO][ND][7][7].

This is done by using numpy.einsum, to make the dot product of the last two 
axis.

Task #7807 (https://gna.org/task/index.php?7807): Speed-up of dispersion 
models for Clustered analysis.

Modified:
    branches/disp_spin_speed/lib/linear_algebra/matrix_exponential.py

Modified: branches/disp_spin_speed/lib/linear_algebra/matrix_exponential.py
URL: 
http://svn.gna.org/viewcvs/relax/branches/disp_spin_speed/lib/linear_algebra/matrix_exponential.py?rev=24156&r1=24155&r2=24156&view=diff
==============================================================================
--- branches/disp_spin_speed/lib/linear_algebra/matrix_exponential.py   
(original)
+++ branches/disp_spin_speed/lib/linear_algebra/matrix_exponential.py   Thu 
Jun 19 17:41:59 2014
@@ -23,7 +23,7 @@
 """Module for the calculation of the matrix exponential."""
 
 # Python module imports.
-from numpy import array, diag, dot, exp
+from numpy import array, any, diag, dot, einsum, eye, exp, iscomplex, 
newaxis, multiply, tile
 from numpy.linalg import eig, inv
 
 # relax module imports.
@@ -55,3 +55,56 @@
     # Return only the real part.
     else:
         return array(eA.real)
+
+
+def matrix_exponential_rankN(A):
+    """Calculate the exact matrix exponential using the eigenvalue 
decomposition approach, for higher dimensional data.
+
+    @param A:   The square matrix to calculate the matrix exponential of.
+    @type A:    numpy float array of rank [NE][NS][NM][NO][ND][7][7]
+    @return:    The matrix exponential.  This will have the same 
dimensionality as the A matrix.
+    @rtype:     numpy float array of rank [NE][NS][NM][NO][ND][7][7]
+    """
+
+    NE, NS, NM, NO, ND, Row, Col = A.shape
+
+    # Is the original matrix real?
+    complex_flag = any(iscomplex(A))
+
+    # The eigenvalue decomposition.
+    W, V = eig(A)
+
+    # W: The eigenvalues, each repeated according to its multiplicity.
+    # The eigenvalues are not necessarily ordered.
+    # The resulting array will be always be of complex type. Shape 
[NE][NS][NM][NO][ND][7]
+    # V: The normalized (unit 'length') eigenvectors, such that the column 
v[:,i]
+    # is the eigenvector corresponding to the eigenvalue w[i]. Shape 
[NE][NS][NM][NO][ND][7][7]
+
+    # Calculate the exponential of all elements in the input array. Shape 
[NE][NS][NM][NO][ND][7]
+    # Add one axis, to allow for broadcasting multiplication.
+    W_exp = exp(W).reshape(NE, NS, NM, NO, ND, Row, 1)
+
+    # Make a eye matrix, with Shape [NE][NS][NM][NO][ND][7][7]
+    eye_mat = tile(eye(7)[newaxis, newaxis, newaxis, newaxis, newaxis, ...], 
(NE, NS, NM, NO, ND, 1, 1) )
+
+    # Transform it to a diagonal matrix, with elements from vector down the 
diagonal.
+    W_exp_diag = multiply(W_exp, eye_mat )
+
+    # Make dot products for higher dimension.
+    # "...", the Ellipsis notation, is designed to mean to insert as many 
full slices (:)
+    # to extend the multi-dimensional slice to all dimensions.
+    dot_V_W = einsum('...ij,...jk', V, W_exp_diag)
+
+    # Compute the (multiplicative) inverse of a matrix.
+    inv_V = inv(V)
+
+    # Calculate the exact exponential.
+    eA = einsum('...ij,...jk', dot_V_W, inv_V)
+
+    # Return the complex matrix.
+    if complex_flag:
+        return array(eA)
+
+    # Return only the real part.
+    else:
+        return array(eA.real)