Finding the self-consistent solution to an equation

Question

At the bottom of this question are a set of functions transcribed from a published neural-network model. When I call R, I get the following error:

RuntimeError: maximum recursion depth exceeded while calling a Python object

Note that within each call to R, a recursive call to R is made for every other neuron in the network. This is what causes the recursion depth to be exceeded. Each return value for R depends on all the others (with the network involving N = 512 total values.) Does anyone have any idea what method should be used to compute the self-consistent solution for R? Note that R itself is a smooth function. I've tried treating this as a vector root-solving problem -- but in this case the 512 dimensions are not independent. With so many degrees of freedom, the roots are never found (using the scipy.optimize functions). Does Python have any tools that can help with this? Maybe it would be more natural to solve R using something like Mathematica? I don't know how this is normally done.

"""Recurrent model with strong excitatory recurrence."""


import numpy as np


l = 3.14


def R(x_i):

    """Steady-state firing rate of neuron at location x_i.

    Parameters
    ----------
    x_i : number
      Location of this neuron.

    Returns
    -------
    rate : float
      Firing rate.

    """

    N = 512
    T = 1

    x = np.linspace(-2, 2, N)
    sum_term = 0
    for x_j in x:
        sum_term += J(x_i - x_j) * R(x_j)

    rate = I_S(x_i) + I_A(x_i) + 1.0 / N * sum_term - T

    if rate < 0:
        return 0

    return rate


def I_S(x):
    """Sensory input.

    Parameters
    ----------
    x : number
      Location of this neuron.

    Returns
    -------
    float
      Sensory input to neuron at x.

    """
    S_0 = 0.46
    S_1 = 0.66
    x_S = 0
    sigma_S = 1.31
    return S_0 + S_1 * np.exp(-0.5 * (x - x_S) ** 2 / sigma_S ** 2)


def I_A(x):
    """Attentional additive bias.

    Parameters
    ----------
    x : number
      Location of this neuron.

    Returns
    -------
    number
      Additive bias for neuron at x.

    """
    x_A = 0
    A_1 = 0.089
    sigma_A = 0.35
    A_0 = 0
    sigma_A_prime = 0.87
    if np.abs(x - x_A) < l:
        return (A_1 * np.exp(-0.5 * (x - x_A) ** 2 / sigma_A ** 2) +
                A_0 * np.exp(-0.5 * (x - x_A) ** 2 / sigma_A_prime ** 2))
    return 0


def J(dx):
    """Connection strength.

    Parameters
    ----------
    dx : number
      Neuron i's distance from neuron j.

    Returns
    -------
    number
      Connection strength.

    """
    J_0 = -2.5
    J_1 = 8.5
    sigma_J = 1.31
    if np.abs(dx) < l:
        return J_0 + J_1 * np.exp(-0.5 * dx ** 2 / sigma_J ** 2)
    return 0


if __name__ == '__main__':

    pass

You could rewrite the algorithm so it is iterative instead of recursive. — spectras, Jul 05 '16 at 00:40
I don't know your specific algorithm, but it has been proven that all recursive algorithm can be expressed iteratively. See this question for details: http://stackoverflow.com/questions/931762/can-every-recursion-be-converted-into-iteration — spectras, Jul 05 '16 at 00:57
@spectras isn't that a bit like saying anything can in theory be computed using a turing machine? it doesn't exactly suggest an implementable solution . . . — abcd, Jul 05 '16 at 01:01
What's your question then? Point you at an alternative algorithm? That's off topic here. Implement it for you? I'm afraid that's too broad a question for SO. — spectras, Jul 05 '16 at 01:07

score 1 · Accepted Answer · answered Jul 05 '16 at 01:11

1

This recursion never ends since there is no termination condition before recursive call, adjusting maximum recursion depth does not help

def R(x_i): 
   ...
   for x_j in x:
       sum_term += J(x_i - x_j) * R(x_j)

Perhaps you should be doing something like

# some suitable initial guess
state = guess

while True: # or a fixed number of iterations
   next_state = compute_next_state(state)

   if some_condition_check(state, next_state):
       # return answer
       return state 

   if some_other_check(state, next_state):
       # something wrong, terminate
      raise ...

answered Jul 05 '16 at 01:11

J.J. Hakala

6,136
6
27
61

thanks, yeah, this is on the right track for sure. i'd imagine there must be formal proposals for optimal ways to achieve this sort of thing. do you know what buzz words i might Google to find algorithms that might help me? maybe there are optimization functions in scipy that would be appropriate here? – abcd Jul 05 '16 at 01:16
i think standard root-solving methods would work here. start w/ an array of 512 values (guesses). instead of making recursive calls w/in `R`, use these values. then take as my error the difference between `R`'s return values and the guess values. continue guessing from there, finding values that minimize the error. – abcd Jul 05 '16 at 01:26
i'll code this up. if you don't mind my fleshing out your answer, i'll do that. otherwise, i'll post my own. – abcd Jul 05 '16 at 01:28
I think OP's algorithm does terminate because the `for` loops terminate. But it does something like `512**512` operations which is kinda steep for 21st century technologies. – tdelaney Jul 05 '16 at 01:28
@tdelaney even if i did all those operations, though, would the code as written produce a sensible result? i don't think it would. – abcd Jul 05 '16 at 01:32

score 0 · Answer 2 · answered Jul 05 '16 at 00:58

0

Change the maximum recursion depth using sys.setrecursionlimit

import sys
sys.setrecursionlimit(10000)

def rec(i):
    if i > 1000:
        print 'i is over 1000!'
        return
   rec(i + 1)

rec(0)

More info: https://docs.python.org/3/library/sys.html#sys.setrecursionlimit`

answered Jul 05 '16 at 00:58

Tomohiro Koana

782
1
7
17

this strikes me as a naive answer, but i could be wrong about that. it seems to me that, because every value of `R` depends on every other value, there is no way for the solution to be computed the way i've written the functions, no matter how deep i allow recursion to go. is this suggestion coming from a place of understanding? – abcd Jul 05 '16 at 00:59
just confirmed this proposal does not solve the problem. – abcd Jul 05 '16 at 01:04
Some time of stopping condition such as `if i > 1000` in the above sample code is needed in every recursive function, which is lacking in `R` – Tomohiro Koana Jul 05 '16 at 01:05
Every call to `R` will lead to 512 call to `R`. It's stuck in an infinite loop. – Tomohiro Koana Jul 05 '16 at 01:09
ah, yes, true. that is the problem. thanks. still, i see no solution. – abcd Jul 05 '16 at 01:11

Finding the self-consistent solution to an equation

2 Answers2