Étant donné une séquence de nombres a [i] pour i = 0 à N-1, j'essaye de calculer la somme suivante:
sum = 0 rolling_product = 1 for i from 0 to (G-1): rolling_product *= a[i] sum += rolling_product for i from G to (N-1): rolling_product /= a[i-G] rolling_product *= a[i] sum += rolling_product
Je voudrais faire la taille G d'un groupe multiplié (dans l'exemple ci-dessus 3) un paramètre variable. Ensuite, le résultat peut être naïvement obtenu en utilisant un simple algorithme O (N * G), qui peut être écrit en pseudocode comme ceci:
sum = 0
for i from 0 to (N-G-1):
group_contribution = 1
for j from 0 to (G-1):
group_contribution *= a[i+j]
sum += group_contribution
Cependant, pour les grands G c'est évident que cet algorithme est terriblement inefficace, surtout en supposant que les nombres de la séquence a [i] ne sont pas connus à l'avance et doivent être calculés de manière coûteuse à l'exécution.
Pour cette raison, j'ai envisagé d'utiliser l'algorithme de complexité O (N + G) suivant qui recycle les valeurs de la séquence a [i] en calculant un produit roulant:
a[0] * a[1] * a[2] + a[1] * a[2] * a[3] + a[2] * a[3] * a[4] + ... a[N-4] * a[N-3] * a[N-2] + a[N-3] * a[N-2] * a[N-1]
Je suis cependant préoccupé par le numérique stabilité de la division dans la représentation en virgule flottante standard.
Je serais intéressé de savoir s'il existe un moyen stable et plus rapide de calculer cette somme. Cela me semble être une tâche numérique de base, mais actuellement je ne sais pas comment le faire efficacement.
Merci pour vos idées!
4 Réponses :
Votre produit roulant est une bonne idée, mais il a des problèmes de stabilité comme vous le dites. Je réglerais cela de cette façon:
a [i] , calculez une somme glissante de log (abs (a [i])) , à l'exclusion des zéros. Ensuite, lorsque vous avez besoin du produit, c'est (num_zeros> 0? 0.0: exp (log_sum)) * signe . Cela réglera l'instabilité majeure. G , écrasez votre llog_sum glissante par ce nombre et remettez-le à zéro. Cela supprimera toutes les erreurs d'arrondi qui s'accumuleraient à long terme. Il me semble que le produit roulant avec division sera un algorithme plus stable que le log et l'exponentiation en termes de précision.
En guise de préambule, vous pouvez envisager d'exécuter des cas de test sur les deux algorithmes et comparer les résultats (comme une erreur relative, par exemple).
Ensuite, si vous avez la mémoire et le temps supplémentaires, voici une approche stable en temps O ( N log 2 G ) et Mémoire. Elle est similaire à l'approche en temps constant, espace linéaireithmique au problème de requête minimale de plage .
Produits de précalcul des plages de puissance de deux
Soit B [ i ] [ j ] le produit de 2 j éléments de a commençant à la position i , donc
B [ i ] [ j ] = a [ i ] × a [ i + 1] × ... × a [ i + 2 j - 1]
Nous nous intéressons aux valeurs N log 2 G en B , à savoir celles pour 0 ≤ j ≤ log 2 G . Nous pouvons calculer chacune de ces valeurs dans O (1), puisque
B [ i ] [ j ] = B [ i ] [ j - 1] × B [ i + 2 j - 1 ] [ j - 1]
Calcul de la somme
Pour calculer un terme dans la somme, nous décomposons G en blocs de taille puissance de deux. Par exemple, si G = 13, alors le premier terme est
a [0] × ... × a [12] = ( a [0] × ... × a [7]) × ( a [8] × ... × a [11]) × a [ 12] = B [0] [3] × B [8] [2] × B [12] [0]
Chacun des termes O ( N ) peut être calculé en temps O (log 2 G ), donc la complexité totale de la recherche la somme est O ( N log 2 G ).
Créez une nouvelle séquence b où b [0] = a [0] * a [1] * a [2] * ... a [G-1] etc.
Maintenant, vous avez un problème plus simple où calculer la somme des valeurs b, vous pouvez maintenir un total et chaque fois que vous ajoutez une valeur, vous soustrayez b [0] et ajoutez la nouvelle valeur et faites-les glisser tous vers le bas (en utilisant circulaire tampon donc rien ne bouge). [Code type de moyenne mobile de fenêtre glissante typique]
Gardez un cache des dernières valeurs G a [] et le calcul de la nouvelle valeur à ajouter à la fin n'est que des opérations O (G) et vous ne calculez un [i] qu'une seule fois.
Oui, si vous calculez soigneusement les produits partiels inversés, vous n'avez pas besoin de diviser.
def window_products(seq, g):
lst = list(seq)
reverse_products = lst[:]
for i in range(len(lst) - 2, -1, -1):
if i % g != len(lst) % g:
reverse_products[i] *= reverse_products[i + 1]
product = 1
for i in range(len(lst) - g + 1):
yield reverse_products[i] * product
if i % g == len(lst) % g:
product = 1
else:
product *= lst[i + g]
print(list(window_products(range(10), 1)))
print(list(window_products(range(10), 2)))
print(list(window_products(range(10), 3)))