numpy vectorisation à la place de la boucle

je l'équation suivante:numpy vectorisation à la place de la boucle

où v, mu sont | R^3, où Sigma est | R^(3x3) et lorsque le résultat est une valeur scalaire . La mise en œuvre cela en numpy a pas de problème:

result = np.transpose(v - mu) @ Sigma_inv @ (v - mu)

Maintenant, j'ai un tas de v-vecteurs (permet de les appeler V \ en | R^3XN) et je que pour exécuter l'équation ci-dessus d'une manière vectorisé de sorte que, comme un résultat, je reçois un nouveau vecteur Result \ in | R^1xn.

# pseudocode 
Result = np.zeros((n, 1)) 
for i,v in V: 
    Result[i,:] = np.transpose(v - mu) @ Sigma_inv @ (v - mu)

J'ai regardé np.vectorize mais la documentation suggère que son la même chose que boucle sur toutes les entrées que je préférerais ne pas faire. Quelle serait une solution vectorisée élégante?

En tant que nœud latéral: n peut être assez grand et une matrice | R^nxn ne rentrera certainement pas dans ma mémoire!

modifier: le code de travail exemple

import numpy as np 

S = np.array([[1, 2], [3,4]]) 
V = np.array([[10, 11, 12, 13, 14, 15],[20, 21, 22, 23, 24, 25]]) 

Res = np.zeros((V.shape[1], 1)) 
for i in range(V.shape[1]): 
    v = np.transpose(np.atleast_2d(V[:,i])) 
    Res[i,:] = (np.transpose(v) @ S @ v)[0][0] 

print(Res)

Source

2017-09-21 juqa

Il serait plus facile si vous fournir un exemple de travail minimum de ce que vous voulez atteindre, même si vous utilisez des boucles (par exemple, y compris 'importation, on déclarations et définitions de tous les symboles que vous utilisez). – norok2

True: J'ai ajouté un exemple minimal – juqa

Je voudrais regarder dans 'np.einsum', bien que je crois que c'est l'une des fonctions les plus énigmatiques de' numpy', il peut être exactement ce que vous voulez. – norok2

En utilisant une combinaison de matrix-multiplication et np.einsum -

np.einsum('ij,ij->j',V,S.dot(V))

Source

2017-09-21 13:01:44 Divakar

Vous avez probablement besoin d'un remodelage à la fin pour correspondre à la MWE, mais à part ça, c'est excellent! – norok2

Est-ce que ce travail pour vous?

res = np.diag(V.T @ S @ V).reshape(-1, 1)

Il semble fournir le même résultat que vous le souhaitez.

import numpy as np 

S = np.array([[1, 2], [3,4]]) 
V = np.array([[10, 11, 12, 13, 14, 15],[20, 21, 22, 23, 24, 25]]) 

Res = np.zeros((V.shape[1], 1)) 
for i in range(V.shape[1]): 
    v = np.transpose(np.atleast_2d(V[:,i])) 
    Res[i,:] = (np.transpose(v) @ S @ v)[0][0] 

res = np.diag(V.T @ S @ V).reshape(-1, 1) 

print(np.all(np.isclose(Res, res))) 
# output: True

Bien qu'il existe probablement une solution plus efficace en mémoire utilisant np.einsum.

Source

2017-09-21 12:57:32 norok2

merci pour les conseils! Le problème avec votre solution est que "n" dans mon application est assez grand, donc une matrice R^{n x n} ne rentrera certainement pas dans la mémoire. – juqa

La solution de @Divarak est supérieure à cela. – norok2

Voici une solution simple:

import numpy as np 

S = np.array([[1, 2], [3,4]]) 
V = np.array([[10, 11, 12, 13, 14, 15],[20, 21, 22, 23, 24, 25]]) 

Res = np.sum((V.T @ S) * V.T, axis=1)

Source

2017-09-21 13:00:44 jdehesa

Mon problème était que V.T @ S @ V serait trop grand pour rentrer dans la mémoire comme il le sera | R^{n x n} – juqa

@juqa Oui, j'ai changé la réponse pour n'utiliser que n x 2 baies. C'est essentiellement la même chose que la réponse de Divakar, bien que l'einsum puisse être plus rapide que la multiplication et la somme. – jdehesa

merci, j'ai benchmarké les deux solutions pour un vecteur avec 5.000.000 entrées et votre intuition était correcte: einsum a pris ~ 100ms alors que votre solution a pris ~ 150ms (benchmarked sur plusieurs itérations). – juqa

Ce sont multiplications de piles matrice/vecteur. numpy.matmul peut le faire après avoir mis S et V dans la forme correcte:

S = S[np.newaxis, :, :] 
VT = V.T[:, np.newaxis, :] 
V = VT.transpose(0, 2, 1)  

tmp = np.matmul(S, V) 
Res = np.matmul(VT, tmp) 

print(Res) 
#[[[2700]] 
# [[3040]] 
# [[3400]] 
# [[3780]] 
# [[4180]] 
# [[4600]]]

Source

2017-09-21 13:06:55 kazemakase

numpy vectorisation à la place de la boucle

Répondre

Questions connexes