opérations de matrice NumPy

Question

Je veux calculer les valeurs suivantes pour tous i et j:

M_ki = Sum[A_ij - A_ik - A_kj + A_kk, 1 <= j <= n] 

Comment puis-je le faire en utilisant Numpy (Python) sans une boucle explicite?

Merci!

Answer 1

Voici une stratégie générale pour résoudre ce genre de problème.

D'abord, écrire un petit script, avec la boucle écrit explicitement dans deux fonctions différentes, et un test à la fin en vous assurant que les deux fonctions sont exactement les mêmes:

import numpy as np 
from numpy import newaxis 

def explicit(a): 
    n = a.shape[0] 
    m = np.zeros_like(a) 
    for k in range(n): 
     for i in range(n): 
      for j in range(n): 
       m[k,i] += a[i,j] - a[i,k] - a[k,j] + a[k,k] 
    return m 

def implicit(a): 
    n = a.shape[0] 
    m = np.zeros_like(a) 
    for k in range(n): 
     for i in range(n): 
      for j in range(n): 
       m[k,i] += a[i,j] - a[i,k] - a[k,j] + a[k,k] 
    return m 

a = np.random.randn(10,10) 
assert np.allclose(explicit(a), implicit(a), atol=1e-10, rtol=0.) 

Ensuite, vectoriser la fonction étape par étape en éditant implicit, l'exécution du script à chaque étape pour vous assurer qu'ils continuent de rester le même:

étape 1

def implicit(a): 
    n = a.shape[0] 
    m = np.zeros_like(a) 
    for k in range(n): 
     for i in range(n): 
      m[k,i] = (a[i,:] - a[k,:]).sum() - n*a[i,k] + n*a[k,k] 
    return m 

Étape 2

def implicit(a): 
    n = a.shape[0] 
    m = np.zeros_like(a) 
    m = - n*a.T + n*np.diag(a)[:,newaxis] 
    for k in range(n): 
     for i in range(n): 
      m[k,i] += (a[i,:] - a[k,:]).sum() 
    return m 

Étape 3

def implicit(a): 
    n = a.shape[0] 
    m = np.zeros_like(a) 
    m = - n*a.T + n*np.diag(a)[:,newaxis] 
    m += (a.T[newaxis,...] - a[...,newaxis]).sum(1) 
    return m 

Et le tour est joué '! Pas de boucles dans le dernier. Pour vectoriser ce genre d'équations, broadcasting est le chemin à parcourir!

Attention: assurez-vous que explicit est l'équation que vous vouliez vectoriser. Je ne savais pas si les termes qui ne dépendent pas de j devraient aussi être additionnés.