无瘤复发

Question

有没有什么方法可以不使用“在无意义”中进行重复呢？

使用np.add和out关键字对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])

输出：[ 1 1 2 3 5 8 13 21 34 55]

但是，如果dtype更改为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])

输出：[ 1. 1. 2. 1. 3. 1. 4. 1. 5. 1.]

谁能告诉我原因吗？或者其他的递归方法？

Answer 1

下面是一种合理有效的使用枕过滤器的方法：

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)

(请注意，我们不能使用lfilter，因为在信号处理方面，这是一个不稳定的过滤器，使Python解释器崩溃，因此我们必须转换为二阶节形式。)

过滤方法的问题是，我不确定它是否可以推广到解决ODE的实际问题。

另一个解决方案是简单地用Python编写循环并使用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

JIT方法的优点是您可以实现各种花哨的东西，但是如果您坚持使用简单的循环和分支而不调用任何(非JITted) Python函数，那么它的工作效果最好。

时间安排结果：

# 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

方法比纯Python循环快，但比JITted循环慢。我想过滤器函数是相对通用的，可以做一些我们在这里不需要的事情，而JITted循环则编译成一个非常小的循环。

Answer 2

您的add适用于此计算，但必须反复应用，以便非零值传播。我不明白您的计算是如何生成[ 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]
 ...

(编辑-注意，第一个np.add的作用就好像它是完全缓冲的。将**的结果与您的对象和浮动数组进行比较。我在Py3上使用1.13.1版。)

或者强调每一个阶段的价值

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]也做同样的事情。

正如ufunc.at文档中所讨论的，像np.add这样的操作使用缓冲，即使使用out参数也是如此。因此，虽然它们在C级别代码中进行交互，但它们不会累积；也就是说，来自ith步骤的结果不会在i+1步骤中使用。

add.at可用于执行未缓冲的a[i] += b。当索引包含dupicate时，这很方便。但它不允许从更改的a值反馈到b。所以它在这里没有用。

它也是一个add.accumulate (又名np.cumsum)，它可以方便地进行某些迭代定义。但在一般情况下很难适用。

cumprod (乘累加)可以使用qmatrix方法。

对大指数给出错误结果的幂函数

使用np.matrix定义qmatrix，以便*乘表示矩阵乘积(相对于元素)：

In [647]: qmatrix = numpy.matrix([[1, 1], [1, 0]])

使用指向此矩阵的副本(实际上是指针)填充对象矩阵。

In [648]: M = np.empty(10, object)
In [649]: M[:] = [qmatrix for _ in range(10)]

应用cumprod，并从每个矩阵中提取一个元素。

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]

链接答案使用numpy.linalg.matrix_power，它也做重复矩阵积。对于单个电源，matrix_power更快，但对于整个功率范围，cumprod更快。

Answer 3

一个可能不是超级有效但很有趣的解决方案就是像这样滥用scipy.linalg.solve_banded。

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

我们怎么写F_0，F_1，F_2.作为一个向量F和恒等式-F_{i-1} -F_i + F{i+1} =0作为一个很好的带状矩阵，然后求解F.注意这可以适用于其他类似的递归。