¿Hay alguna manera de hacer recurrencia sin usar for’s in numpy?
Usando np.add out palabra clave, haz el truco con dtype="int"
import numpy as np N = 100 fib = np.zeros(N, dtype=np.int) fib[:2] = 1. np.add(fib[:-2], fib[1:-1], out=fib[2:]) print(fib[:10])
salida: [ 1 1 2 3 5 8 13 21 34 55]
Sin embargo, si dtype se cambia a np.float
import numpy as np N = 100 fib = np.zeros(N, dtype=np.float) fib[:2] = 1. np.add(fib[:-2], fib[1:-1], out=fib[2:]) print(fib[:10])
salida: [ 1. 1. 2. 1. 3. 1. 4. 1. 5. 1.]
¿Podría alguien decirme por qué? o alguna otra forma de hacer recursion?
Aquí hay un método razonablemente eficiente usando un filtro scipy:
import numpy as np from scipy import signal def scipy_fib(n): x = np.zeros(n, dtype=np.uint8) x[0] = 1 sos = signal.tf2sos([1], [1, -1, -1]) return signal.sosfilt(sos, x)
(Tenga en cuenta que no podemos usar lfilter porque, en términos de procesamiento de la señal, este es un filtro inestable que hace que el intérprete de Python se bloquee, por lo que tenemos que convertirlo al formato de secciones de segundo orden).
El problema con el enfoque de filtrado es que no estoy seguro de si se puede generalizar al problema real de resolver una EDO.
Otra solución es simplemente escribir el bucle en Python y acelerarlo con un comstackdor JIT :
import numba as nb @nb.jit(nopython=True) def numba_fib(n): y = np.empty(n) y[:2] = 1 for i in range(2, n): y[i] = y[i-1] + y[i-2] return y
La ventaja del enfoque JIT es que puede implementar todo tipo de cosas sofisticadas, pero funciona mejor si se adhiere a bucles y twigs simples sin tener que llamar a ninguna función de Python (sin JIT).
Resultados de tiempo:
# Baseline without JIT: %timeit numba_fib(10000) # 100 loops, best of 3: 5.47 ms per loop %timeit scipy_fib(10000) # 1000 loops, best of 3: 719 µs per loop %timeit numba_fib(10000) # 10000 loops, best of 3: 33.8 µs per loop # For fun, this is how Paul Panzer's solve_banded approach compares: %timeit banded_fib(10000) # 1000 loops, best of 3: 1.33 ms per loop
El enfoque del filtro scipy es más rápido que el bucle de Python puro pero más lento que el bucle JITted. Supongo que la función de filtro es relativamente genérica y hace cosas que no necesitamos aquí, mientras que el bucle JITted se comstack en un bucle muy pequeño.
scipy.linalg.solve_banded solución no muy eficiente pero divertida sería abusar de scipy.linalg.solve_banded como tal
import numpy as np from scipy import linalg N = 50 a = np.zeros((N,)) + np.array([[1, -1, -1]]).T b = np.zeros((N,)) b[0] = 1 linalg.solve_banded((2, 0), a, b) # array([ 1.00000000e+00, 1.00000000e+00, 2.00000000e+00, # 3.00000000e+00, 5.00000000e+00, 8.00000000e+00, # 1.30000000e+01, 2.10000000e+01, 3.40000000e+01, # 5.50000000e+01, 8.90000000e+01, 1.44000000e+02, # 2.33000000e+02, 3.77000000e+02, 6.10000000e+02, # 9.87000000e+02, 1.59700000e+03, 2.58400000e+03, # 4.18100000e+03, 6.76500000e+03, 1.09460000e+04, # 1.77110000e+04, 2.86570000e+04, 4.63680000e+04, # 7.50250000e+04, 1.21393000e+05, 1.96418000e+05, # 3.17811000e+05, 5.14229000e+05, 8.32040000e+05, # 1.34626900e+06, 2.17830900e+06, 3.52457800e+06, # 5.70288700e+06, 9.22746500e+06, 1.49303520e+07, # 2.41578170e+07, 3.90881690e+07, 6.32459860e+07, # 1.02334155e+08, 1.65580141e+08, 2.67914296e+08, # 4.33494437e+08, 7.01408733e+08, 1.13490317e+09, # 1.83631190e+09, 2.97121507e+09, 4.80752698e+09, # 7.77874205e+09, 1.25862690e+10])
Cómo funciona esto es que escribimos F_0, F_1, F_2 … como un vector F y las identidades -F_ {i-1} -F_i + F {i + 1} = 0 como una matriz que está bien agrupada y luego resolvemos para F. Tenga en cuenta que esto se puede adaptar a otras recurrencias similares.
Dado que las recurrencias lineales tienen soluciones analíticas ( aquí para Fibonacci), otro método rápido de escisión es:
def fib_scipy2(N): r5=math.sqrt(5) phi=(1+r5)/2 a= (-phi*ones(N)).cumprod() return (np.abs(a)-1/a)/r5
Carreras :
In [413]: fib_scipy2(8) Out[413]: array([ 1., 1., 2., 3., 5., 8., 13., 21.]) In [414]: %timeit fib(10**4) 103 µs ± 888 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
Puedes ajustarlo a cualquier ecuación lineal.
Su add funciona para este cálculo, pero debe aplicarse repetidamente, de modo que los valores distintos de cero se propaguen. No veo cómo se generó su cálculo [ 1. 1. 2. 1. 3. 1. 4. 1. 5. 1.] .
In [619]: fib=np.zeros(10,int);fib[:2]=1 In [620]: for _ in range(10): ...: np.add(fib[:-2], fib[1:-1], out=fib[2:]) ...: print(fib) ...: [1 1 2 1 0 0 0 0 0 0] # ** [1 1 2 3 3 1 0 0 0 0] [1 1 2 3 5 6 4 1 0 0] [ 1 1 2 3 5 8 11 10 5 1] [ 1 1 2 3 5 8 13 19 21 15] ...
(Editar: np.add cuenta que el primer np.add actúa como si np.add completamente almacenado en búfer. Compare el resultado en ** con sus dos matrices de objetos y flotantes. Estoy usando la versión 1.13.1 en Py3).
O para resaltar los buenos valores en cada etapa.
In [623]: fib=np.zeros(20,int);fib[:2]=1 In [624]: for i in range(10): ...: np.add(fib[:-2], fib[1:-1], out=fib[2:]) ...: print(fib[:(i+2)]) ...: [1 1] [1 1 2] [1 1 2 3] [1 1 2 3 5] [1 1 2 3 5 8] [ 1 1 2 3 5 8 13] [ 1 1 2 3 5 8 13 21] [ 1 1 2 3 5 8 13 21 34] [ 1 1 2 3 5 8 13 21 34 55] [ 1 1 2 3 5 8 13 21 34 55 89]
fib[2:] = fib[:-2]+fib[1:-1] hace lo mismo.
Como se explica en la documentación de ufunc.at , las operaciones como np.add utilizan el almacenamiento en búfer, incluso con el parámetro out . Así que, si bien ellos intervienen en el código de nivel C , no se acumulan; es decir, los resultados del ith paso no se utilizan en el paso i+1 .
add.at se puede utilizar para realizar a[i] += b búfer a[i] += b . Eso es útil cuando los índices contienen duplicados. pero no permite la retroalimentación de los valores de a a b . Así que no es útil aquí.
También es un add.accumulate (también np.cumsum como np.cumsum ) que puede ser útil para ciertas definiciones iterativas. Pero es difícil de aplicar en casos generales.
cumprod (multiplicar acumular) puede trabajar con el enfoque qmatrix
La función matrix_power de Numpy da resultados incorrectos para grandes exponentes
Use np.matrix para definir qmatrix , de modo que * multiplicar signifique el producto de la matriz (a diferencia de elementwise):
In [647]: qmatrix = numpy.matrix([[1, 1], [1, 0]])
Rellene una matriz de objetos con copias (punteros en realidad) a esta matriz.
In [648]: M = np.empty(10, object) In [649]: M[:] = [qmatrix for _ in range(10)]
Aplique cumprod y extraiga un elemento de cada matriz.
In [650]: m1=np.cumprod(M) In [651]: m1 Out[651]: array([matrix([[1, 1], [1, 0]]), matrix([[2, 1], [1, 1]]), matrix([[3, 2], [2, 1]]), matrix([[5, 3], [3, 2]]), matrix([[8, 5], [5, 3]]), matrix([[13, 8], [ 8, 5]]), matrix([[21, 13], [13, 8]]), matrix([[34, 21], [21, 13]]), matrix([[55, 34], [34, 21]]), matrix([[89, 55], [55, 34]])], dtype=object) In [652]: [x[0,1] for x in m1] Out[652]: [1, 1, 2, 3, 5, 8, 13, 21, 34, 55]
La respuesta vinculada utiliza numpy.linalg.matrix_power que también hace productos matriciales repetidos. Para una sola potencia, matrix_power es más rápido, pero para una amplia gama de poderes, el cumprod es más rápido.