tamaño - suma y multiplicacion de matrices en python

NumPy Broadcasting: cálculo de la suma de las diferencias al cuadrado entre dos matrices

Question

tengo el siguiente código. Está tomando una eternidad en Python. Debe haber una manera de traducir este cálculo en una transmisión ...

def euclidean_square(a,b):
    squares = np.zeros((a.shape[0],b.shape[0]))
    for i in range(squares.shape[0]):
        for j in range(squares.shape[1]):
            diff = a[i,:] - b[j,:]
            sqr = diff**2.0
            squares[i,j] = np.sum(sqr)
    return squares

Answer 1

Otra solución además de usar cdist es la siguiente

difference_squared = np.zeros((a.shape[0], b.shape[0]))
for dimension_iterator in range(a.shape[1]):
    difference_squared = difference_squared + np.subtract.outer(a[:, dimension_iterator], b[:, dimension_iterator])**2.

Answer 2

Puede usar np.einsum después de calcular las diferencias de forma broadcasted way , de esta manera:

ab = a[:,None,:] - b
out = np.einsum(''ijk,ijk->ij'',ab,ab)

O use scipy''s cdist con su argumento métrico opcional establecido como ''sqeuclidean'' para darnos las distancias cuadradas euclidianas según sea necesario para nuestro problema, así:

from scipy.spatial.distance import cdist
out = cdist(a,b,''sqeuclidean'')

Answer 3

Recolecté los diferentes métodos propuestos aquí, y en other dos questions , y medí la velocidad de los diferentes métodos:

import numpy as np
import scipy.spatial
import sklearn.metrics

def dist_direct(x, y):
    d = np.expand_dims(x, -2) - y
    return np.sum(np.square(d), axis=-1)

def dist_einsum(x, y):
    d = np.expand_dims(x, -2) - y
    return np.einsum(''ijk,ijk->ij'', d, d)

def dist_scipy(x, y):
    return scipy.spatial.distance.cdist(x, y, "sqeuclidean")

def dist_sklearn(x, y):
    return sklearn.metrics.pairwise.pairwise_distances(x, y, "sqeuclidean")

def dist_layers(x, y):
    res = np.zeros((x.shape[0], y.shape[0]))
    for i in range(x.shape[1]):
        res += np.subtract.outer(x[:, i], y[:, i])**2
    return res

# inspired by the excellent https://github.com/droyed/eucl_dist
def dist_ext1(x, y):
    nx, p = x.shape
    x_ext = np.empty((nx, 3*p))
    x_ext[:, :p] = 1
    x_ext[:, p:2*p] = x
    x_ext[:, 2*p:] = np.square(x)

    ny = y.shape[0]
    y_ext = np.empty((3*p, ny))
    y_ext[:p] = np.square(y).T
    y_ext[p:2*p] = -2*y.T
    y_ext[2*p:] = 1

    return x_ext.dot(y_ext)

# https://.com/a/47877630/648741
def dist_ext2(x, y):
    return np.einsum(''ij,ij->i'', x, x)[:,None] + np.einsum(''ij,ij->i'', y, y) - 2 * x.dot(y.T)

Utilizo timeit para comparar la velocidad de los diferentes métodos. Para la comparación, uso vectores de longitud 10, con 100 vectores en el primer grupo y 1000 vectores en el segundo grupo.

import timeit

p = 10
x = np.random.standard_normal((100, p))
y = np.random.standard_normal((1000, p))

for method in dir():
    if not method.startswith("dist_"):
        continue
    t = timeit.timeit(f"{method}(x, y)", number=1000, globals=globals())
    print(f"{method:12} {t:5.2f}ms")

En mi computadora portátil, los resultados son los siguientes:

dist_direct   5.07ms
dist_einsum   3.43ms
dist_ext1     0.20ms  <-- fastest
dist_ext2     0.35ms
dist_layers   2.82ms
dist_scipy    0.60ms
dist_sklearn  0.67ms

Si bien los dos métodos dist_ext1 y dist_ext2 , ambos basados ​​en la idea de escribir (xy)**2 como x**2 - 2*x*y + y**2 , son muy rápidos, hay una desventaja: cuando la distancia entre x e y es muy pequeño, debido a un error de cancelación, el resultado numérico a veces puede ser (muy levemente) negativo.