矢量迭代固定点搜索

Question

我有一个函数总是会收敛到一个固定点，例如f(x)= (x-a)/2+a。我有一个函数可以通过重复调用该函数来找到这个固定点：

def find_fix_point(f,x):
    while f(x)>0.1:
        x = f(x)
    return x

效果很好，现在我想为矢量版本执行此操作；

def find_fix_point(f,x):
    while (f(x)>0.1).any():
        x = f(x)
    return x

但是，如果大多数实例只需要大约 10 次迭代而一个需要 1000 次迭代，那么这是非常低效的。有什么快速的方法可以删除已经找到的 `x？

代码可以使用numpy或scipy。

Answer 1

解决这个问题的一种方法是使用递归：

def find_fix_point_recursive(f, x):
    ind = x > 0.1
    if ind.any():
        x[ind] = find_fix_point_recursive(f, f(x[ind]))
    return x

通过这个实现，我们只在需要更新的点上调用f。

请注意，通过使用递归，我们避免了一直检查x > 0.1，每次调用都处理越来越小的数组。

%timeit x = np.zeros(10000); x[0] = 10000; find_fix_point(f, x)
1000 loops, best of 3: 1.04 ms per loop

%timeit x = np.zeros(10000); x[0] = 10000; find_fix_point_recursive(f, x)
10000 loops, best of 3: 141 µs per loop

Answer 2

首先，为了一般性，我更改了标准以适应固定点定义：我们在|x-f(x)|<=epsilon 时停止。

你可以混合boolean indexing and integer indexing来保持每次的活跃点。这是一种方法：

def find_fix_point(f,x,epsilon):
    ind=np.mgrid[:len(x)] # initial indices.
    while ind.size>0:
        xind=x[ind]  # integer indexing
        yind=f(xind)
        x[ind]=yind
        ind=ind[abs(yind-xind)>epsilon]  # boolean indexing

一个有很多固定点的例子：

from matplotlib.pyplot import plot,show
x0=np.linspace(0,1,1000) 
x = x0.copy() 
def f(x): return x*np.sin(1/x)     
find_fix_point(f,x,1e-5)        
plot(x0,x,'.');show()    

Answer 3

一般的方法是使用布尔索引来仅计算尚未达到平衡的那些。

我调整了 Jonas Adler 给出的算法以避免最大递归深度：

def find_fix_point_vector(f,x):
    x = x.copy()
    x_fix = np.empty(x.shape)
    unfixed = np.full(x.shape, True, dtype = bool)
    while unfixed.any():
        x_new = f(x) #iteration
        x_fix[unfixed] = x_new # copy the values
        cond = np.abs(x_new-x)>1
        unfixed[unfixed] = cond #find out which ones are fixed
        x = x_new[cond] # update the x values that still need to be computed
    return x_fix

编辑： 在这里，我回顾了提出的 3 个解决方案。我将根据他们的提议者find_fix_Jonas, find_fix_Jurg, find_fix_BM 调用不同的功能。我根据 BM 更改了所有函数中的定点条件（参见最后 Jonas 的更新函数）。

速度：

%timeit find_fix_BM(f, np.linspace(0,100000,10000),1)
100 loops, best of 3: 2.31 ms per loop

%timeit find_fix_Jonas(f, np.linspace(0,100000,10000))
1000 loops, best of 3: 1.52 ms per loop

%timeit find_fix_Jurg(f, np.linspace(0,100000,10000))
1000 loops, best of 3: 1.28 ms per loop

根据可读性，我认为 Jonas 的版本是最容易理解的，所以应该在速度不是很重要的情况下选择。

但是，当到达固定点之前的迭代次数很大 (>1000) 时，Jonas 的版本可能会引发 Runtimeerror。其他两种解决方案没有这个缺点。

B.M.的版本但是可能比我提出的版本更容易理解。

#

使用的 Jonas 版本：

def find_fix_Jonas(f, x):
    fx = f(x)
    ind = np.abs(fx-x)>1
    if ind.any():
        fx[ind] = find_fix_Jonas(f, fx[ind])
    return fx

Answer 4

...删除已经找到的`x？

根据您的条件使用布尔索引创建一个新数组。

>>> a = np.array([3,1,6,3,9])
>>> a != 3
array([False,  True,  True, False,  True], dtype=bool)
>>> b = a[a != 3]
>>> b
array([1, 6, 9])
>>>