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)
_______________________________________________
relax (http://www.nmr-relax.com)
This is the relax-commits mailing list
relax-commits@xxxxxxx
To unsubscribe from this list, get a password
reminder, or change your subscription options,
visit the list information page at
https://mail.gna.org/listinfo/relax-commits
Re: r24156 - /branches/disp_spin_speed/lib/linear_algebra/matrix_exponential.py