Le tableau numpy 3D A contient une série (dans cet exemple, je choisis 3) de tableau numpy 2D D de forme 2 x 2 code >. La matrice D est la suivante: D = np.array([[1,2],[3,4]])
idx = np.arange(3)
A = np.zeros((3,2,2))
A[idx,:,:] = D # This gives A = [[[1,2],[3,4]],[[1,2],[3,4]],\
# [[1,2],[3,4]]]
# In mathematical notation: A = {D, D, D}
n = np.array([0,1,2])
Result = np.power(A,n) # ------> Not the desired output.
A est initialisé et affecté comme ci-dessous: idx = np.arange(3)
A = np.zeros((3,2,2))
A[idx,:,:] = D # This gives A = [[[1,2],[3,4]],[[1,2],[3,4]],\
# [[1,2],[3,4]]]
# In mathematical notation: A = {D, D, D}
4 Réponses :
D'accord, j'ai passé beaucoup de temps sur ce problème mais je n'arrivais pas à trouver une solution vectorisée comme vous le souhaiteriez. Je voudrais donc plutôt proposer d'abord une solution de base, puis peut-être une optimisation si vous avez besoin de trouver des puissances continues.
La fonction que vous recherchez s'appelle numpy.linalg.matrix_power
def all_powers_matrix(D, *rangeargs): #end exclusive
''' Expects 2D matrix.
Use as all_powers_matrix(D, end) or
all_powers_matrix(D, start, end)
'''
if len(rangeargs) == 1:
start = 0
end = rangeargs[0]
elif len(rangeargs) == 2:
start = rangeargs[0]
end = rangeargs[1]
else:
print("incorrect args")
return None
result = np.zeros((end - start, *D.shape))
result[0] = np.linalg.matrix_power(A[0], start)
for i in range(start + 1, end):
result[i] = np.linalg.multi_dot([result[i - 1], D])
return result
return result
result = all_powers_matrix(D, 3)
#Output:
array([[[ 1., 0.],
[ 0., 1.]],
[[ 1., 2.],
[ 3., 4.]],
[[ 7., 10.],
[15., 22.]]])
Cependant, si vous remarquez, vous finissez par calculer plusieurs puissances pour la même matrice de base. Nous pourrions plutôt utiliser les résultats intermédiaires et partir de là, en utilisant numpy.linalg.multi_dot
def all_powers_arr_of_matrix(A):
result = np.zeros(A.shape)
result[0] = np.linalg.matrix_power(A[0], 0)
for i in range(1, A.shape[0]):
result[i] = np.linalg.multi_dot([result[i - 1], A[i]])
return result
result = all_powers_arr_of_matrix(A)
#Output:
array([[[ 1., 0.],
[ 0., 1.]],
[[ 1., 2.],
[ 3., 4.]],
[[ 7., 10.],
[15., 22.]]])
De plus, nous pouvons éviter de créer la matrice A entièrement, en économisant un certain temps.
import numpy as np
D = np.matrix([[1,2],[3,4]])
idx = np.arange(3)
A = np.zeros((3,2,2))
A[idx,:,:] = D # This gives A = [[[1,2],[3,4]],[[1,2],[3,4]],\
# [[1,2],[3,4]]]
# In mathematical notation: A = {D, D, D}
np.zeros(A.shape)
n = np.array([0,1,2])
result = [np.linalg.matrix_power(D, i) for i in n]
np.array(result)
#Output:
array([[[ 1, 0],
[ 0, 1]],
[[ 1, 2],
[ 3, 4]],
[[ 7, 10],
[15, 22]]])
Notez que vous devrez ajouter la gestion des erreurs si vous décidez d'utiliser ces fonctions telles quelles.
Merci d'avoir pris votre temps pour résoudre ce problème. Je le mets pour timeit test et je vais vous expliquer le résultat.
Un produit cumulatif existe dans numpy, mais pas pour les matrices. Par conséquent, vous devez créer votre propre fonction «matcumprod». Vous pouvez utiliser np.dot pour cela, mais np.matmul (ou @ ) est spécialisé pour la multiplication matricielle.
Puisque vous indiquez que vos puissances vont toujours de 0 à some_power, je suggère la fonction suivante :
def matcumprod(D, upto):
Res = np.empty((upto, *D.shape), dtype=A.dtype)
Res[0, :, :] = np.eye(D.shape[0])
Res[1, :, :] = D.copy()
for i in range(1,upto):
Res[i, :, :] = Res[i-1,:,:] @ D
return Res
Au fait, une boucle surpasse souvent une fonction numpy intégrée si cette dernière utilise beaucoup de mémoire, alors ne vous inquiétez pas si vos pouvoirs restez dans les limites ...
Je n'ai pas de solution complète, mais je voulais mentionner certaines choses qui sont un peu trop longues pour les commentaires.
Vous pouvez d'abord examiner la l'exponentiation de la chaîne d'addition si vous calculez de grandes puissances de grandes matrices. Il s'agit essentiellement de demander combien de multiplications matricielles sont nécessaires pour calculer A ^ k pour un k donné. Par exemple A ^ 5 = A (A ^ 2) ^ 2 donc vous n'avez besoin que de trois multiplications matricielles: A ^ 2 et (A ^ 2) ^ 2 et A (A ^ 2) ^ 2 . C'est peut-être le moyen le plus simple de gagner en efficacité, mais vous devrez probablement toujours utiliser des boucles explicites.
Votre question est également liée au problème de calcul de Ax, A ^ 2x, ..., A ^ kx pour un A et un x
Pour calculer la puissance de la matrice D , une façon pourrait être de trouver les valeurs propres et les vecteurs propres droits de celle-ci avec np.linalg.eig puis augmentez la puissance de la matrice diagonale comme c'est plus simple, puis après quelques manipulations, vous pouvez utiliser deux np.einsum pour calculer A
#get eigvalues and eigvectors
eigval, eigvect = np.linalg.eig(D)
# to check how it works, you can do:
print (np.dot(eigvect*eigval,np.linalg.inv(eigvect)))
#[[1. 2.]
# [3. 4.]]
# so you get back on D
#use power as ufunc of outer with n on the eigenvalues to get all the one you want
arrp = np.power.outer( eigval, n).T
#apply_along_axis to create the diagonal matrix along the last axis
diagp = np.apply_along_axis( np.diag, axis=-1, arr=arrp)
#finally use two np.einsum to calculate with the subscript to get what you want
A = np.einsum('lij,jk -> lik',
np.einsum('ij,kjl -> kil',eigvect,diagp), np.linalg.inv(eigvect)).round()
print (A)
print (A.shape)
#[[[ 1. 0.]
# [-0. 1.]]
#
# [[ 1. 2.]
# [ 3. 4.]]
#
# [[ 7. 10.]
# [15. 22.]]]
#
#(3, 2, 2)
Pouvez-vous expliquer le résultat attendu? is This
D0 = [[1,0], [0,1]], D1 = [[1,2], [3,4]], D2 = [[7,10], [15,22 ]]la sortie attendue? Si tel est le cas, cela ne semble pas être une fonction de puissance.Je veux D0, D1, D2 dans A avec ce résultat. Avant la mise sous tension: A = [[[1,2], [3,4]], [[1,2], [3,4]], [[1,2], [3,4]]] . Après application de la puissance, il devrait être A = [[[1,0], [0,1]], [[1,2], [3,4]], [[7,10], [15,22] ]].
pourquoi le fait d'élever
[[1,2], [3,4]]à zéro le transforme-t-il en[[1,0], [0,1]]. Ne devrait-il pas être[[1,1], [1,1]]? Quelle est la logique de l'opération?Non, cela devrait se révéler être une matrice d'identité. Dans le cas d'un nombre scalaire, a (b ^ 0) = a car b ^ 0 = 1. C'est la propriété d'identité. De même pour la matrice, toute matrice carrée à la puissance 0 est une matrice d'identité, c'est-à-dire dans ce cas [[1,0], [0,1]] (disons que c'est B). Donc, si vous multipliez maintenant une matrice A par B (AB), vous obtenez A, qui peut également être vérifié manuellement.
Ah je vois, merci d'avoir expliqué