python - Recurrence with numpy

Question

Is there any way to make recurrence without using for's in numpy?

Using np.add with out keyword do the trick with 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])

output: [ 1 1 2 3 5 8 13 21 34 55]

However, if dtype is changed to 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])

output: [ 1. 1. 2. 1. 3. 1. 4. 1. 5. 1.]

Could someone tell me why? or any other way to make recursion?

Answer 1

One maybe not super-efficent but fun solution would be to abuse scipy.linalg.solve_banded like so

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])

How this works is that we write F_0, F_1, F_2... as one vector F and the identities -F_{i-1} -F_i + F{i+1} = 0 as a matrix which is nicely banded and then solve for F. Note that this can adapted to other similar recurrences.