Trouver efficacement une somme de fonctions éparses

Question

J'ai une situation statistique où je veux trouver la somme de certaines fonctions extrêmement pointues. Autrement dit, étant donné un ensemble de points d'entrée in et les points de sortie out, je veux trouver les numéros sum_in f(in,out), où f est extrêmement fortement pointu.

Pour concrétude, la situation pourrait être quelque chose comme ceci:

import numpy as np 

sample_pts = 10 ** 7 
data_pts = 10 ** 5 

mu = np.random.rand(data_pts) 
x = np.linspace(0, 1, sample_pts) 

def f(mu, x): 
    return np.exp(-1e10 * ((mu - x) ** 2)) 

Ma solution actuelle est juste pour itérer sur les points d'échantillonnage, avec vectorisation sur les points de données:

results = np.zeros(sample_pts) 
for i in range(sample_pts): 
    results[i] = np.sum(f(mu, x[i])) 

Cependant, c'est incroyablement inefficace parce que la grande majorité des nombres calculés ici sont très petits - en fait, beaucoup sont zéro à la précision en virgule flottante! Il devrait y avoir assez de place ici pour obtenir une accélération d'un facteur de mille, au moins. Quel est le moyen le plus rapide pour calculer cette somme en numpy, en supposant qu'une petite erreur numérique (disons un sur un milliard) est acceptable?

Answer 1

Avec cette fonction, vous ne devez itérer

In [833]: mu = np.random.rand(10) 
In [834]: x = np.linspace(0,1,6) 
In [836]: def f(mu, x): 
    ...:  return np.exp(.1* ((mu - x) ** 2)) 
    ...: 
In [837]: f(mu, np.arange(10)) # same shape 
Out[837]: 
array([ 1.00005667e+00, 1.00151080e+00, 1.33174582e+00, 
     2.27563858e+00, 3.14399507e+00, 9.37132782e+00, 
     1.53439871e+01, 6.54667741e+01, 4.85267134e+02, 
     3.11160087e+03]) 
In [838]: f(mu, x[:,None]) # 'outer' broadcasting 
Out[838]: 
array([[ 1.00005667, 1.079973 , 1.00949403, 1.00175693, 1.03860883, 
     1.00729568, 1.06179883, 1.0288728 , 1.00184352, 1.00010102], 
     [ 1.00310927, 1.04691816, 1.00115406, 1.00045585, 1.01741263, 
     1.00048473, 1.03353998, 1.01118532, 1.00041336, 1.00283372], 
     [ 1.01425284, 1.0230266 , 1.00085791, 1.00718177, 1.00465417, 
     1.00170149, 1.01411376, 1.00178422, 1.00700916, 1.01365075], 
     [ 1.03375727, 1.00770978, 1.00859845, 1.02209706, 1.00002398, 
     1.01097526, 1.00304502, 1.00044212, 1.02179017, 1.032814 ], 
     [ 1.06209967, 1.00059511, 1.02456265, 1.04556437, 1.00341039, 
     1.0285303 , 1.00006571, 1.0071267 , 1.0451157 , 1.06079202], 
     [ 1.0999839 , 1.0015108 , 1.04913917, 1.07816138, 1.01489503, 
     1.05479487, 1.00510399, 1.02199931, 1.0775598 , 1.09827913]]) 

et la somme

In [839]: _.sum(axis=1) 
Out[839]: 
array([ 10.22980131, 10.11750708, 10.08823266, 10.1412531 , 
     10.27786262, 10.50142738]) 

j'ai changé la constante; avec le -1e10 tous les termes 0 (avec dans une tolérance).

Je ne connais pas de moyen de compression du calcul autre que l'élimination des valeurs mu et x qui sont connues pour produire de très petits résultats.

Avec le paramètre -1e10, je dois donner un x très proche de mu pour obtenir une somme non nulle:

In [857]: f(mu, mu-.0001).sum() 
Out[857]: 3.7200759760847501e-43