Runtime complexity, Recursion, and estimating Feigenbaum's constant for chaos¶
Preamble: Run the cells below to import the necessary Python packages
This notebook created by William Gilpin. Consult the course website for all content and GitHub repository for raw files and runnable online code.
import numpy as np
from IPython.display import Image, clear_output, display
# Import local plotting functions and in-notebook display functions
import matplotlib.pyplot as plt
%matplotlib inline
Finding bifurcations in the logistic map¶
In previous lectures we saw how to use the logistic map to model population growth. We also saw that the logistic map can exhibit chaotic behavior for certain values of the parameter $r$.
In this lecture we will explore the behavior of the logistic map in more detail. In particular, we will look at how the behavior of the logistic map changes as we vary the parameter $r$.
First we will implement classes to represent the logistic map and calculate the bifurcation diagram. Unlike previous lectures, we will implement the Logistic map using an iterable class.
In classes, Python reserves certain method names for special purposes. For example, the __init__
method is called when an object is created. The __str__
method is called when an object is printed. The __call__
method allows an object to be treated like a function.
In a similar vein, objects can be used in for
loops. The __iter__
method is called when an object is used in a for
loop or implicitly in a list comprehension. The __next__
method is called to get the next value in the loop. We will use these methods to implement the logistic map.
class LogisticMap:
"""
Implements the logistic map.
Attributes:
r: The growth rate parameter
x: The current value of x for the logistic map.
"""
def __init__(self, r=3.8, x0=0.65):
self.r = r
self.x = x0
def next(self):
self.x = self.r * self.x * (1 - self.x)
return self.x
def simulate(self, n):
"""Simulate the map for n iterations."""
x = []
for _ in range(n):
x.append(self())
return x
# The __call__ method allows us to call an instance of this class like a function.
# It must return the value that should be returned when the instance is called.
def __call__(self):
return self.next()
# The __str__ method is called when the object is printed. It must return a string.
def __str__(self):
return 'Logistic map with r = {:.2f}'.format(self.r)
# map.__str__()
# map.next()
class BifurcationDiagram:
"""
Find the bifurcation diagram for a map using __iter__ in Python. This stateful
implementation is more efficient than the stateless implementation we used previously
because it does not need to store the entire trajectory in memory.
This implementation will allow us to avoid storing the entire bifurcation diagram in
memory. We can iterate over the map, one value of r at a time, and plot the last
values of x for each value of r after a transient has elapsed
Attributes:
dmap: A discrete map object. Must contain a simulate method that returns a
trajectory of the map.
rmin: The minimum value of r to use.
rmax: The maximum value of r to use.
rsteps: The number of steps to take between rmin and rmax.
n: The number of iterations to use for each value of r.
transient: The number of iterations to discard as transient.
"""
def __init__(self, dmap, rmin, rmax, rsteps, n, transient=100):
self.dmap = dmap
self.rmin = rmin
self.rmax = rmax
self.rsteps = rsteps
self.n = n
self.transient = transient
# The __iter__ method makes this class iterable. It must return an object
# that implements the __next__ method. Iterable objects can be used in for
# loops and list comprehensions.
def __iter__(self):
return self
# The __next__ method is called by the for loop to get the next item in the
# iteration. It must return the next item, or raise StopIteration if there
# are no more items.
def __next__(self):
self.dmap.r = self.rmin
if self.rmin > self.rmax:
raise StopIteration
else:
r = self.rmin
self.rmin += self.rsteps
return self.dmap.simulate(self.n)[self.transient:], r
map = LogisticMap()
# invoke the __print__ method
print(map)
# one step of the map
print(f"Initial condition: {map.x}")
print(f"First iteration: {map()}")
# plot the trajectory
traj = map.simulate(1000)
plt.figure(figsize=(9, 6))
plt.plot(traj)
plt.xlim(0, len(traj))
plt.xlabel('Iteration')
plt.ylabel('x')
Logistic map with r = 3.80 Initial condition: 0.65 First iteration: 0.8644999999999998
Text(0, 0.5, 'x')
diagram = BifurcationDiagram(LogisticMap(), 2.4, 4.0, 0.001, 1000)
plt.figure(figsize=(9, 6))
for traj, r in diagram:
plt.plot([r] * len(traj), traj, ',k', alpha=0.25)
plt.xlabel('r')
plt.ylabel('x')
plt.xlim(2.4, 4.0)
(2.4, 4.0)
When does the logistic map undergo period doubling? We will use a bisecting search to find the values of $r$ at which the period doubles.
# naive search: line scan
def check_period_doubling(prev_traj, current_traj):
"""
Check if the period has doubled from the previous trajectory to the current
trajectory. Uses a basic test of whether the number of unique values in the
current trajectory is greater than the number of unique values in the previous
trajectory.
Args:
prev_traj (array): Previous trajectory
current_traj (array): Current trajectory
threshold (float): A threshold value to determine the similarity between points
Returns:
bool: True if period doubled, False otherwise
"""
prev_unique_vals = len(set(np.round(prev_traj, decimals=4)))
current_unique_vals = len(set(np.round(current_traj, decimals=4)))
return np.abs(current_unique_vals / prev_unique_vals - 2) < 0.3
def line_scan(dmap, rmin, rmax, n_rvals=100, transient=50):
"""
Find the doubling points for a map using a line scan
Scan the values of r between rmin and rmax and find the values of r where the
map doubles in period. This function uses the discrete map's simulate method
to generate a trajectory for each value of r.
Args:
dmap (object): A discrete map object. Must contain a simulate method that
returns a trajectory of the map.
rmin (float): The minimum value of r to use.
rmax (float): The maximum value of r to use.
n_rvals (int): The number of values of r to use between rmin and rmax.
transient (int): The number of iterations to discard as transient.
Returns:
array: An array of values of r where the map doubles in period.
"""
rvals = np.linspace(rmin, rmax, n_rvals)
doubling_rvals = []
prev_traj = None
for r in rvals:
dmap.r = r
traj = np.array(dmap.simulate(1000)[-transient:])
if prev_traj is not None:
if check_period_doubling(prev_traj, traj):
doubling_rvals.append(r)
prev_traj = traj
return np.array(doubling_rvals)
import time
start = time.time()
rvals_c = line_scan(LogisticMap(), 2.4, 3.57, n_rvals=300)
print(rvals_c)
end = time.time()
print(f"Time elapsed: {end - start} seconds")
# # filter out values that change by less than 0.01
# rvals_c = rvals_c[:-1][np.diff(rvals_c) > (3.6 - 2.4) / (1000 - 1)]
# plt.figure()
# gaps = rvals_c[1:] - rvals_c[:-1]
# plt.plot(gaps)
# print(gaps[:-1] / gaps[1:])
# vals = [3.0]
# gap = 0.45252525
# for i in range(5):
# nxt = vals[-1] + gap
# vals.append(nxt)
# gap = gap / 4.667
# rvals_c = np.array(vals)
# gaps = rvals_c[1:] - rvals_c[:-1]
# plt.plot(gaps)
[2.99869565 3.44869565 3.54652174 3.56608696] Time elapsed: 0.1146402359008789 seconds
diagram = BifurcationDiagram(LogisticMap(), 2.4, 4.0, 0.001, 1000)
plt.figure(figsize=(9, 6))
for traj, r in diagram:
plt.plot([r] * len(traj), traj, ',k', alpha=0.25)
plt.plot(rvals_c, [0] * len(rvals_c), 'r.')
plt.xlabel('r')
plt.ylabel('x')
# plt.xlim(2.4, 3.6)
plt.xlim(2.8, 3.6)
(2.8, 3.6)
We can now create a plot of how large the gaps in $r$ are between successive bifurcations.
# filter out values that change by less than 0.01
# rvals_c = rvals_c[:-1][np.diff(rvals_c) > (3.6 - 2.4) / (10000 - 1)]
plt.figure()
gaps = rvals_c[1:] - rvals_c[:-1]
plt.xlabel('Doubling index')
# in latex r_c
plt.ylabel('$r^*_{i+1} - r^*_i$')
plt.plot(gaps)
print(f"Ratios between successive doublings: {gaps[:-1] / gaps[1:]}")
Ratios between successive doublings: [4.6 5. ]
Determining Big-O for an algorithm¶
- It's almost always some combination of a polynomial $N^{\alpha}$ and $\log(N)$. Remember that overwriting values, peeking the front of a queue, etc take constant time or $\mathcal{O}(1)$
Some useful times to know¶
- $\log{(N)}$ usually means you can do something clever (recursion, bisection search, etc)
- $\text{exp}\left(N\right)$ is rare and usually not optimal (describes a process that branches)
- Stirling's approximation: $N! \sim \sqrt{N} e^{-N} N^N$ or $\log{(N!)} \sim N \log{(N)}$
- Sorting a list via mergesort is $\mathcal{O}(N \log(N))$
- Binary search, two-sum are pointers with $\log{(N)}$ time
- Dot product of two length-$N$ vectors is $\mathcal{O}(N)$. $N \times N$ matrix-vector product is $\mathcal{O}(N^2)$. Product of two $N \times N$ matrices is $\mathcal{O}(N^3)$
Binary search¶
- A binary search is a search algorithm that works on a sorted list. It is $\mathcal{O}\left(\log{(N)}\right)$ because it cuts the list in half each time.
def binary_search(a, target):
"""
Find the index of the first element in a *sorted* array that is greater than or equal to target.
We are using a two-pointer approach where one pointer is always greater than the
target and the other is always less than the target. The pointers start at the ends
of the array and move towards each other until they meet. The index of the left
pointer is the index of the first element greater than or equal to the target.
"""
lo = 0
hi = len(a)
while lo < hi:
mid = (lo + hi) // 2
if a[mid] < target:
lo = mid + 1
else:
hi = mid
return lo
a = np.array([1, 2, 3, 4, 7, 8, 9, 10, 15, 43, 99])
print(binary_search(a, 5))
4
import timeit
nvals = np.arange(2, 800)
nvals = np.arange(2, 5000, 3)
all_times = list()
for n in nvals:
all_reps = list()
all_times.append(
timeit.timeit(
"binary_search(np.arange(n), n // 2)", globals=globals(), number=10000
)
)
plt.figure(figsize=(9, 6))
plt.plot(nvals, all_times, 'k')
plt.xlabel('n')
plt.ylabel('Time (s)')
Text(0, 0.5, 'Time (s)')
Improving our calculation of doubling values¶
Now let's return to our logistic map. We will use a binary search to find the values of $r$ at which the period doubles. By using a binary search, we can avoid scanning the large set of intermediate values of $r$ for which the period does not double.
There a are a few nuances to our approach compared to traditional binary search.
We need to perform a three-pointer comparison among the periods of the logistic map at $r_{min}$, $r_{middle}$, and $r_{max}$. This is because the period doubling may not occur exactly at $r$. We thus implement a three-way comparison function.
We need to find several values of $r$ corresponding to critical doubling values, rather than just one value. This requires multiple binary searches over different intervals, as well as a heuristic for modifying the intervals to ensure that we don't repeatedly find the same root. We break down this process as follows:
We search the full $r$ interval.
Upon finding a doubling value, we shrink the interval so that its rightmost value is just below the root we just found by some tolerance $\epsilon$.
We then search the new interval for the next doubling value.
We repeat this process until we have found all the doubling values. This heuristic will fail if the period doubling values are closer than $\epsilon$ apart. In this case, we would need to reduce $\epsilon$ and try again. This same difficulty manifests in methods used to find multiple roots (zero crossings) of continuous functions.
## Three-way comparison function
def check_period_doubling2(left_traj, middle_traj, right_traj):
"""
Check if the period has doubled from the previous trajectory to the current
trajectory. Uses a basic test of whether the number of unique values in the
current trajectory is greater than the number of unique values in the previous
trajectory.
Args:
prev_traj (array): Previous trajectory
current_traj (array): Current trajectory
threshold (float): A threshold value to determine the similarity between points
Returns:
bool: True if period doubled, False otherwise
"""
t_left = len(set(np.round(left_traj, decimals=4)))
t_middle = len(set(np.round(middle_traj, decimals=4)))
t_right = len(set(np.round(right_traj, decimals=4)))
factor_right = t_middle / t_left
factor_left = t_right / t_middle
return factor_right > factor_left
def binary_search_between(dmap, rmin, rmax, transient=50):
"""
Find a doubling between rmin and rmax using a binary search
Scan the values of r between rmin and rmax and find the values of r where the
map doubles in period. This function uses the discrete map's simulate method
to generate a trajectory for each value of r.
"""
# print(f"Searching between {rmin} and {rmax}")
rvals = []
prev_traj = None
while rmax - rmin > 1e-6:
r = (rmax + rmin) / 2
dmap.r = r
traj = np.array(dmap.simulate(1000)[-transient:])
dmap.__init__(r=rmin)
traj1 = np.array(dmap.simulate(1000)[-transient:])
dmap.__init__(r=rmax)
traj2 = np.array(dmap.simulate(1000)[-transient:])
dmap.__init__(r=r)
traj_m = np.array(dmap.simulate(1000)[-transient:])
where_double = check_period_doubling2(traj1, traj_m, traj2)
if where_double:
rmax = r
else:
rmin = r
return r
def binary_search_bifurcation(dmap, rmin, rmax, n_iter, transient=50):
"""
Find doublings using a multipointer binary search
Scan the values of r between rmin and rmax and find the values of r where the
map doubles in period. This function uses the discrete map's simulate method
to generate a trajectory for each value of r.
"""
rvals = []
# narrow down the search between every pair of successive pointers to find a doubling
for i in range(n_iter):
rmax_new = binary_search_between(dmap, rmin, rmax, transient=transient)
if np.abs(rmax - rmax_new) < 1e-3:
rmax -= 1e-4
else:
rvals.append(rmax)
rmax = rmax_new
return np.array(rvals)[::-1]
rvals_c = binary_search_bifurcation(LogisticMap(), 2.8, 3.6, 10, transient=50)
plt.figure()
gaps = rvals_c[1:] - rvals_c[:-1]
plt.xlabel('Doubling index')
# in latex r_c
plt.ylabel('$r^*_{i+1} - r^*_i$')
plt.plot(gaps)
print(f"Ratios between successive doublings: {gaps[:-1] / gaps[1:]}")
Ratios between successive doublings: [4.67187821 4.69514424 0.57676821]
What's going on physically?¶
Many chaotic systems approach chaos through a series of bifurcations, which occur at spacings on the bifurcation diagram equal to Feigenbaum's constant. If $r$ is the parameter of the logistic map, and $r^*_i$ and $r^*_{i+1}$ are the critical parameter values of successive doubling bifurcations, then the ratio of successive bifurcation intervals approaches Feigenbaum's constant:
$$\lim_{i \rightarrow \infty} \frac{r^*_{i+1} - r^*_i}{r^*_{i+2} - r^*_{i+1}} = \delta \approx 4.66920\ldots$$
The gap between successive doublings therefore decreases as a geometric series, with a ratio of $\delta$. A geometric series decreasing by $\delta$ converges to a finite value, which here corresponds to the onset of chaos at $r \approx 3.56995$ in the logistic map.
Recursion and dynamic programming¶
- Problems where the solution to the size $N$ input can be written in terms of the size $N - 1$, $N - 2$, etc inputs
- For example, the Fibonacci sequence: $F_N = F_{N - 1} + F_{N - 2}$, with $F_0 = 0$, $F_1 = 1$
Key ingredients¶
- Base cases: $0! = 1$, $x_0 = 1$ (initial conditions), etc
- Inductive reasoning: write a case for $N$ in terms of $N - 1$, $N - 2$, etc
Recursion vs Dynamic programming¶
Both use a "divide and conquer" approach, but they work in opposite directions
Dynamic programming: start from base cases, and work your way up to $N$
Recursion: start from $N$, then do $N - 1$, etc, until you hit a base case
When working with graphs or trees, we often encounter recursion due to the need to perform depth-first search
When simulating dynamical systems, or processes that don't go backwards, we often run into DP
Dynamic programming is typically $\mathcal{O}(N)$ runtime, $\mathcal{O}(1)$ memory. Recursion is ideally $\mathcal{O}(Log(N))$ runtime, but for the problems we'll see below it's $\mathcal{O}(N)$ runtime, $\mathcal{O}(N)$ memory because we can't skip any steps
Recursion¶
Source: https://betterprogramming.pub/recursive-functions-2b5ce4610c81
We first apply recursion to compute the factorial function, which can be expensive to compute for large $N$.
We implement two approaches to computing this function: recursion and dynamic programming.
Note that the factorial function can be written as follows:
$$N! = N \times (N-1)!$$
$$N_i = i\, N_{i-1}$$
where $i$ starts at $0$ and $0! = 1$ by convention.
def factorial(n):
"""Compute n factorial via recursion"""
# Base case
if n == 0:
return 1
# Recursive case
else:
return n * factorial(n - 1)
### Run a timing test
import timeit
nvals = np.arange(2, 10)
all_times = list()
for n in nvals:
all_times.append(
timeit.timeit("factorial(n)", globals=globals(), number=1000000)
)
plt.figure()
plt.plot(nvals, all_times, 'k')
[<matplotlib.lines.Line2D at 0x2923fafb0>]
def factorial(n):
"""Compute n factorial, for n >= 0, with dynamic programming."""
# Base case
if n in [0, 1]:
return 1
nprev = 1
for i in range(1, n + 1):
nprev *= i
return nprev
def factorial(n):
if n < 0:
return None
prev = 1
for i in range(1, n + 1):
prev = prev * i
return prev
import timeit
nvals = np.arange(2, 10)
all_times = list()
for n in nvals:
all_times.append(
timeit.timeit("factorial(n)", globals=globals(), number=1000000)
)
plt.figure()
plt.plot(nvals, all_times, 'k')
[<matplotlib.lines.Line2D at 0x2946dfbe0>]
Questions¶
Why did we only run timing tests for inputs up to $N = 50$? What complicates our analysis if you try larger values of $N$?
import timeit
nvals = np.arange(2, 800)
all_times = list()
for n in nvals:
all_times.append(
timeit.timeit("factorial(n)", globals=globals(), number=10)
)
plt.figure()
plt.plot(nvals, all_times, 'k')
[<matplotlib.lines.Line2D at 0x2923c3520>]
It turns out that multiplying large floats has unfavorable scaling with N. Gradeschool is $O(N^2)$ where $N$ is the number of digits, and Karatsuba is $O(N^{1.585})$.
We next apply this approach to computing the Fibonacci sequence, which is famously defined by the recurrence relation
$$\text{Fib}(N) = \text{Fib}(N - 1) + \text{Fib}(N - 2)$$
The first few Fibonacci numbers are $1, 1, 2, 3, 5, 8, 13, 21, 34, 55, \ldots$
def fibonacci(n):
"""Compute the nth Fibonacci number via recursion"""
# Base cases
if n == 0:
return 0
elif n == 1:
return 1
# Recursive case
else:
return fibonacci(n - 1) + fibonacci(n - 2)
nvals = np.arange(2, 30)
all_times = list()
for n in nvals:
all_times.append(
timeit.timeit("fibonacci(n)", globals=globals(), number=10)
)
plt.figure(figsize=(9, 6))
plt.semilogy(nvals, all_times)
plt.xlabel("n")
plt.ylabel("Time (s)")
Text(0, 0.5, 'Time (s)')
def fibonacci(n):
"""Compute the nth Fibonacci number via dynamic programming."""
if n == 0:
return 0
elif n == 1:
return 1
n1, n2 = 0, 1
for i in range(2, n + 1):
n1, n2 = n2, n1 + n2
return n2
nvals = np.arange(2, 300)
all_times = list()
for n in nvals:
all_times.append(
timeit.timeit("fibonacci(n)", globals=globals(), number=1)
)
plt.plot(nvals, all_times)
[<matplotlib.lines.Line2D at 0x2920dbb20>]
Questions¶
- Unlike the factorial function, the recursive approach is much slower than the dynamic programming approach, even though they both have the same asymptotic complexity. Why is this?
Recursion versus dynamic programming as discrete-time dynamical systems¶
We can also think of recursion and dynamic programming as dynamical systems. In this case, the input $N$ is the time variable, and the output is the state of the system at time $N$.
For example, the factorial function can be written as the dynamical system
$$x_{N} = x_{N - 1} \cdot N$$
with initial condition $x_0 = 1$.
The Fibonacci sequence can be written as the dynamical system
$$x_{N} = x_{N - 1} + x_{N - 2}$$
with initial conditions $x_0 = 1$ and $x_1 = 1$.
Dynamic programming corresponds to simulating the dynamical system forward in time, starting from the initial condition. Recursion corresponds to simulating the dynamical system backward in time, starting from the final condition.
Depth-first search versus breadth-first search (solving a maze)¶
Breadth-first search (BFS)¶
Source: https://www.geeksforgeeks.org/program-for-nth-fibonacci-number/
Depth-first search (DFS)¶
- Worst case: $\mathcal{O}(N_v + N_e)$, where $N_v$, $N_e$ are the number of vertices and edges, respectively
- The number of edges of a hypercube: $N_e = d N_v / 2$, where $d$ is the dimension
- Implies worst-case $\mathcal{O}(N)$ for dfs in square lattice with N vertices
Source: https://www.geeksforgeeks.org/program-for-nth-fibonacci-number/
Looking ahead to Homework 1: The Abelian sandpile¶
You can implement toppling as either BFS or DFS, or using a purely iterative approach
Physically, what do these different scenarios represent?