python – Suma numéricamente estable y rápida de los últimos K elementos en secuencia

Pregunta:

Supongamos que tengo una secuencia larga, posiblemente infinita $ x: = [x_1, x_2, …] $ , y quiero usarla para calcular otra secuencia $ y: = [y_1, y_2, …] $ donde cada elemento es la suma de los últimos K elementos de la secuencia de entrada. es decir

$ y_i = \ sum_ {j = max (1, iK)} ^ i x_j $

La forma ingenua e ineficiente de hacer (en Python) esto sería:

def sum_of_last_k(x: Sequence[float], k: int):
    buffer = [0. ] * k  # Initialize a buffer of k zeros
    for i, xi in enumerate(x):
        buffer[i % k] = xi  # Where % is modulo
        yield sum(buffer)

… Lo que tendría $ \ mathcal O (K) $ eficiencia por iteración. Sin embargo, quiero hacerlo en línea y de manera eficiente. Podríamos hacer esto en $ \ mathcal O (1) $ manteniendo una suma acumulada:

def sum_of_last_k(x: Sequence[float], k: int):
   buffer = [0. ] * (k-1)  # Initialize a buffer of k-1 zeros
   running_sum = 0
   for i, xi in enumerate(x):
       ix_wrapped = i % (k-1)
       old_value = buffer[ix_wrapped]
       buffer[ix_wrapped] = xi
       running_sum = running_sum + xi - old_value
       yield running_sum

Esto tiene una eficiencia de $ \ mathcal O (1) $ , pero todavía hay un problema: debido a errores de precisión de punto flotante, puede haber un error de redondeo en running_sum que se acumula con el tiempo.

Ahora sé que podría resolver esto running_sum calcular running_sum desde cero cada M iteraciones, donde M es un número grande, y hacer que el tiempo de ejecución promedio sea $ \ mathcal O ((M-1 + K) / M) $ que puede ser muy cerca de 1 … Pero quiero que el peor tiempo de ejecución sea $ \ mathcal O (1) $ , no $ \ mathcal O (K) $ .

Entonces, ¿hay alguna manera de calcular esto en $ \ mathcal O (1) $ tiempo por iteración sin dejar de ser numéricamente estable?

Respuesta:

¡Pregunta muy interesante!

Estrategia adaptativa inspirada en LAPACK

Esto me recuerda un error que se encontró en una rutina LAPACK (QR revelador de rango) relacionada con normas de 'desactualización': esencialmente, se le da la norma de un vector v , y desea calcular en cada iteración en O (1 ) la norma del mismo vector después de cortar su entrada inicial: v[1:], v[2:], ... (en notación Python). Hay transformaciones ortogonales entre estas sumas que dificultan el uso de diferentes estrategias. Consulte http://www.netlib.org/lapack/lawnspdf/lawn176.pdf para obtener más detalles.

Por lo que tengo entendido, la solución que adoptaron los investigadores fue rastrear la acumulación de error de punto flotante explícitamente y volver a calcular la suma cuando es demasiado grande.

Las versiones simples de esta udea no son demasiado complicadas de implementar: por ejemplo, realice un seguimiento del número más grande que apareció en su cálculo y vuelva a calcular la suma desde cero si el resultado actual es significativamente menor que él. Esta estrategia sigue siendo O (K) en el peor de los casos, pero al menos paga O (K) solo cuando realmente la necesita por razones de estabilidad y O (1) en caso contrario.

O (\ sqrt {K}) estrategia de tiempo

De todos modos, para ir determinísticamente por debajo de O (K), una posible idea es la siguiente:

• suponga que $ K $ es un cuadrado perfecto, por simplicidad.
• divide los elementos cíclicamente en $ \ sqrt {K} $ cubos: el cubo $ j $ contiene todos los elementos x[i] con i % sqrt(K) == j . Para cada cubo, memorice su contribución a la suma acumulada en una variable temporal.
• en cada paso, cuando cambia la ventana, todos los depósitos no cambian excepto uno de ellos, en el que se reemplaza un elemento: vuelva a calcular la suma de este depósito desde cero; luego calcule la suma de sumas. Ambas operaciones requieren $ O (\ sqrt {K}) $ tiempo.

Tenga en cuenta que aún necesita espacio O (K) con esta estrategia si solo tiene acceso en línea a la secuencia, pero en cada iteración solo accede a $ O (\ sqrt {K}) $ .

Posiblemente, una versión recursiva de esta estrategia con tamaños de cubos elegidos adecuadamente podría llegar a $ O (\ log K) $ por $ K $ suficientemente grandes.