矢量迭代不动点搜索

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

一般的方法是使用布尔索引来计算那些还没有达到平衡的。你知道吗

我调整了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的版本可能会引发Runtimeerror，当到达fixpoint之前的迭代次数很大时（>；1000）。其他两种解决方案没有这个缺点。你知道吗

然而，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 3

首先为了通用性，我修改了标准以适应定点定义：当|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()