Vectoriced iterative fixpoint search

Question

I have a function that will always converge to a fixpoint, e.g. f(x)= (x-a)/2+a. I have a function that will find this fixpoint through repetive invoking of the function:

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

Which works fine, now I want to do this for a vectoriced version;

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

However this is quite inefficient, if most of the instances only need about 10 iterations and one needs 1000. What is a fast method to remove `x that already have been found?

The code can use numpy or scipy.

There's an example of f in the first sentence. `f(x) = (x-a)/2+a` which will have fixpoint `a` — Jürg W. Spaak, Aug 18 '17 at 14:26
Well then there is no problem, since they all converge the same? — Jonas Adler, Aug 18 '17 at 14:29
@JonasAdler No, they don't, e.g. `x = 2**np.arange(100)`. The entries will converge with `np.arange(100)` iterations. — Jürg W. Spaak, Aug 18 '17 at 14:30
Is `f(x)= (x-a)/2+a` the actual function that you are working with or is it just a sample here? — Divakar, Aug 18 '17 at 14:37
The largest value in `x` will always determine the number of iterations whether you remove elements that meet the criteria or not. — wwii, Aug 18 '17 at 15:17

Jonas Adler · Answer 1 · 2017-08-18T14:42:00.187

1

One way to solve this would be to use recursion:

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

With this implementation, we only call f on the points which need to be updated.

Note that by using recursion we avoid having to do the check x > 0.1 all the time, with each call working on smaller and smaller arrays.

%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

edited Aug 18 '17 at 14:42

answered Aug 18 '17 at 14:34

Jonas Adler

10,365
5
46
73

That's a great solution (especially as it's super small). However this will only work when the number of iterations needed is small (less than about 1000), because afterwards the maximal recursion depth has been reached. @Jonas Furthermore this [post](https://stackoverflow.com/questions/3323001/what-is-the-maximum-recursion-depth-in-python-and-how-to-increase-it) says recursion is not very efficient (however somebody in the comments says it is...) – Jürg W. Spaak Aug 18 '17 at 18:00

score 1 · Answer 2 · answered Aug 19 '17 at 09:10

First for generality,I change the criteria to fit with the fix-point definition : we stop when |x-f(x)|<=epsilon.

You can mix boolean indexing and integer indexing to keep each time the active points. Here a way to do that :

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

An example with a lot of fix points :

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

Thanks for your solution (which is a lot easier to understand. However is about half as fast as the one proposed by me, have a look at the review. — Jürg W. Spaak, Aug 21 '17 at 09:04

Jürg W. Spaak · Accepted Answer · 2017-08-21T09:09:33.830

The general method is to use boolean indexing to compute only the ones, that did not yet reach equilibrium.

I adapted the algorithm given by Jonas Adler to avoid maximal recursion depth:

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

Edit: Here I review the 3 solutions proposed. I will call the different fucntions according to their proposer, find_fix_Jonas, find_fix_Jurg, find_fix_BM. I changed the fixpoint condition in all functions according to BM (see updated fucntion of Jonas at the end).

Speed:

%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

According to readability I think the version of Jonas is the easiest one to understand, so should be chosen when speed does not matter very much.

Jonas's version however might raise a Runtimeerror, when the number of iterations until fixpoint is reached is large (>1000). The other two solutions do not have this drawback.

The verion of B.M. however might be easier to understand than the version proposed by me.

#

Version of Jonas used:

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

score -2 · Answer 4 · answered Aug 18 '17 at 15:08

-2

...remove `x that already have been found?

Create a new array using boolean indexing with your condition.

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

answered Aug 18 '17 at 15:08

wwii

23,232
7
37
77

yes, I'm quite sure this does not help with my problem. Your comment to the question is correct, but I'm of course working with a function, that can't be solved analytically, so there's no priori knowledge to how often I'll have to use the iteration. And even if there were, I'd still like to get rid of the ones that already reached the fix point – Jürg W. Spaak Aug 18 '17 at 17:54

Vectoriced iterative fixpoint search

4 Answers4