Search code examples
pythonscipymathematical-optimizationnumerical-methodsrecurrence

How do I fit recurrence relations using scipy.optimize.fmin_l_bfgs_b?


I recently watched a video by Veritasium where he went through a fascinating discussion of the logistic map. It got me thinking about recurrence relations, and how parametric recurrence relations might be fit to data.

I would like to parametrize $\theta$ in $\hat y_{n+1} = \theta \hat y_{n}(1-\hat y_{n})$ on a sequence $y_k \in [0, 1]$ where $k=n+1$ by minimizing the mean square error using L-BFGS-B algorithm available in Scipy. This example would be instructive for me to generalize to fitting other recurrence relations to real-world data. How do I implement the objective function where the predicted values are the output of a recurrence relation in such a way that it can be passed to the func argument of fmin_l_bfgs_b?


Solution

  • If I understand your question correctly, you want to find the value of $\theta$ minimizing $\sum_{n=0}^{N-1} (\hat y_{n+1} - \theta \hat y_{n} (1 - \hat y_{n}))^2$, given a sequence ${\hat y_k}_{k=0}^N$. If so, assuming your data is ys and your initial guess x0, you could do that through

    def f(l): 
        t = l[0] 
        return ((ys[1:] - (t * ys[:-1] * (1 - ys[:-1])))**2).sum()
    
    fmin_l_bfgs_b(f, x0=(x0,), approx_grad=True)
    

    For example, if we create some data for which theta is approximately 3:

    In [43]: import numpy as np 
        ...: ys = [0.3] 
        ...: theta = 3 
        ...: for _ in range(100): 
        ...:     ys.append((np.random.uniform(-0.02, 0.02) + theta)*ys[-1] * (1 - ys[-1])) 
        ...: ys = np.array(ys) 
        ...:                                                                                        
    
    In [44]: def f(l): 
        ...:     t = l[0] 
        ...:     return ((ys[1:] - (t * ys[:-1] * (1 - ys[:-1])))**2).sum() 
        ...: fmin_l_bfgs_b(f, x0=(0.5,), approx_grad=True)                                          
    Out[44]: 
    (array([2.99949454]),
     0.0006258145273212467,
     {'grad': array([-5.70908338e-07]),
      'task': b'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL',
      'funcalls': 6,
      'nit': 2,
      'warnflag': 0})
    

    Here, of course, you could also supply the gradient; I was just a bit lazy.

    However, if that is indeed what you are trying to do, chances are you'll want something tailored to least squares problems (such as Levenberg--Marquardt); in SciPy, such methods are available in scipy.optimize.least_squares. With those, your problem boils down to the following:

    def F(t): 
        return ys[1:] - (t * ys[:-1] * (1 - ys[:-1]))
    
    least_squares(F, x0=x0)
    

    With the data from above:

    In [53]: def F(t): 
        ...:     return ys[1:] - (t * ys[:-1] * (1 - ys[:-1])) 
        ...:                                                                                        
    
    In [54]: least_squares(F, x0=0.5)                                                               
    Out[54]: 
     active_mask: array([0.])
            cost: 0.00031290726365087974
             fun: ...
            grad: array([-2.43698605e-09])
             jac: ...
         message: '`gtol` termination condition is satisfied.'
            nfev: 4
            njev: 4
      optimality: 2.4369860459044074e-09
          status: 1
         success: True
               x: array([2.9994946])