分数阶累加的Python实现

分数阶累加是分数阶差分的逆运算，它不仅可用于分数阶差分方程的分析，也可以用于建立分数阶灰色模型。然而许多初学者在动手实现分数阶灰色模型时经常发现非常困难，究其原因其实是对定义公式的分析不够，对相应程序语言的特性不熟悉。

本文将从分数阶累加的定义出发，深入分析其计算过程，结合Python语言的特性，详细讲解其实现过程。

1、分数阶累加的定义

对任意原始序列 x(0)(1),x(0)(2),...,x(0)(n)x^{(0)}(1),x^{(0)}(2),...,x^{(0)}(n)x(0)(1),x(0)(2),...,x(0)(n)，其分数阶累加定义为：
x(r)(k)=∑i=1k(k−i+r−1k−i)x(0)(i)(1)x^{(r)}(k)=\sum_{i=1}^{k}\left(\begin{array}{c} k-i+r-1 \\ k-i \end{array}\right) x^{(0)}(i) \tag{1} x(r)(k)=i=1∑k(k−i+r−1k−i)x(0)(i)(1)

其中：
(k−i+r−1k−i)=(r+k−i−1)(r+k−i−2)⋯(r+1)r(k−i)!(2)\left(\begin{array}{c} k-i+r-1 \\ k-i \end{array}\right)=\frac{(r+k-i-1)(r+k-i-2) \cdots(r+1) r}{(k-i) !} \tag{2} (k−i+r−1k−i)=(k−i)!(r+k−i−1)(r+k−i−2)⋯(r+1)r(2)

通俗地讲，分数阶累加实际上就是一种加权求和，而在具体实现时关键是求出其加权系数 (2)，该系数称为广义二项式系数或广义牛顿二项式系数。通常初学者在看到该定义的第一反应即是用循环来实现，这种方式直接简单，但效率较低。尤其在使用类似Python这类脚本语言时其速度非常之慢。

2、广义牛顿二项式系数的计算

要实现分数阶累加的快速算法，关键的问题就在于如何快速算出其加权系数。对于一些不常用科学计算类程序的人而言，通常不太容易想到简单的解决方法。事实上，Python中提供了一个函数binom正好可以解决该问题。

binom函数是Scipy库中的一个函数。在Scipy的官方文档中对该函数并没有过多的介绍，那么我们先来简单写几个测试看看它的用法：

from scipy.special import binom# integers
print(binom(10,3))# decimals
print(binom(10,0.3))# array_like data
print(binom([5,8,10],[0.1,2,3.1]))

运行结果：

120.0
2.246507496036618
[  1.24554362  28.         129.20102313]

从上述代码中不难看出，binom函数支持的输入是十分简单的，第一个变量是二项式系数中的n，第二个变量是k，即是在n中选出k个的组合数。同时它也支持小数和数组，这为我们进一步简化计算过程提供了极大的便利。

3、分数阶累加的矩阵形式及向量计算

本文既然是要讨论快速算法，那么自然就会尽量避免循环。虽然在目前看来，完全可以直接使用binom函数对其进行计算，但这样仍然要使用2层循环，并且计算n(n+1)2\frac{n(n+1)}{2}2n(n+1)次。试想假设我们有10个点要计算，那么就得调用5050次，这个量级对于Python的循环是很低效的。

实际上Python中的numpy库是具有十分完善的矩阵计算功能，其计算效率非常高。我们要想办法利用这一特性，那么很自然地就需要引入它的矩阵定义：

[x(r)(1),x(r)(2),⋯,x(r)(n)]=[x(0)(1),x(0)(2),⋯,x(0)(n)][1(r1)⋯(r+n−3n−2)(r+n−2n−1)01⋯(r+n−4n−3)(r+n−3n−2)⋮⋮⋮⋮⋮00⋯1(r1)00⋯01](3)\left[ \begin{matrix} x^{(r)}(1), x^{(r)}(2), \cdots, x^{(r)}(n) \end{matrix} \right] =\left[ \begin{matrix} x^{(0)}(1), x^{(0)}(2),\cdots , x^{(0)}(n)\end{matrix} \right] \left[\begin{matrix} 1 & \left(\begin{matrix} r \\ 1 \end{matrix}\right) & \cdots & \left(\begin{matrix} r+n-3 \\ n-2 \end{matrix}\right) & \left(\begin{matrix} r+n-2 \\ n-1 \end{matrix}\right) \\ 0 & 1 & \cdots & \left(\begin{matrix} r+n-4 \\ n-3 \end{matrix}\right) & \left(\begin{matrix} r+n-3 \\ n-2 \end{matrix}\right) \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & \cdots & 1 & \left(\begin{matrix} r \\ 1 \end{matrix}\right)\\ 0 & 0 & \cdots & 0 &1 \end{matrix}\right] \tag{3} [x(r)(1),x(r)(2),⋯,x(r)(n)]=[x(0)(1),x(0)(2),⋯,x(0)(n)]⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡10⋮00(r1)1⋮00⋯⋯⋮⋯⋯(r+n−3n−2)(r+n−4n−3)⋮10(r+n−2n−1)(r+n−3n−2)⋮(r1)1⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤(3)

这里有一个特别要注意的点：处在同一右对角线上的所有元素都相等。进一步我们可以发现一个更直接有利于我们编程的事实：从右向左的每一列都相当于其右列去掉最后一个元素，再向左移动一位。那么我们容易得出一个非常简单且有用的结论：我们在编程实现时只需要计算最后一列。此时的时间复杂度已经从O(n2)O(n^2)O(n2)降到了线性时间复杂度O(n)O(n)O(n)，同时，其内存消耗也由原来的n2n^2n2单位变成了nnn个单位。

要计算最后一列系数，实际上只需要考虑定义式 (3) 中k=nk=nk=n的情形，利用Python的向量计算格式，可以很容易实现这一目的：

import numpy as np
from scipy.special import binomr = 0.1
k = 10
i = np.linspace(1,k,k)fnum = k - i + r -1
ki = k - icoefs = binom(fnum,ki)print(coefs)

运行结果：

[0.01447564 0.01608405 0.01812287 0.02079674 0.02446675 0.02983750.0385     0.055      0.1        1.        ]

4、分数阶累加的实现

结合上述分析，这里我们得到的系数coefs实际上是矩阵式 (3) 中的第一行。再次观察矩阵式注意到在计算累加时是采用原序列的行左乘矩阵的列的方式进行。从编程的角度考虑则不难发现要计算第kkk个分数阶累加值，即是将前kkk个原序列点中的值组成的向量对应乘以系数的倒数前kkk个值再求和。这一步实现的时候就需要用到一层循环了。

x = np.random.rand(10)
xr = np.zeros(10)
for k in range(10):xr[k] = sum(x[:k+1]*coefs[-(k+1):])print(x)
print(xr)

运行结果：

[0.42343781 0.53237537 0.75292602 0.32316064 0.6931979  0.496077060.74748506 0.3630722  0.29325613 0.95218414]
[0.42343781 0.57471915 0.82945264 0.44403624 0.80005567 0.638403230.89195739 0.5386026  0.45048115 1.09707706]

5、分数阶累加的完整代码

最后为方便使用，我们将该方法封装为完整的函数

from scipy.special import binom
import numpy as np# compute the r-order accumulation
def fago(x,r):n = len(x)i = np.linspace(1,n,n)coefs = binom(n - i + r -1,n - i)xr = np.zeros(10)for k in range(10):xr[k] = sum(x[:k+1]*coefs[-(k+1):])return xr

简单测试一下这个函数：

import time t1 = time.time()
print(fago(x,0.1))
t2 = time.time()print('Runtime:',t2-t1,'s')

运行结果：

array([0.42343781, 0.57471915, 0.82945264, 0.44403624, 0.80005567,0.63840323, 0.89195739, 0.5386026 , 0.45048115, 1.09707706])
Runtime: 0.0010225772857666016 s

可以看到和刚才的结果一模一样，运行时间只有0.001秒。

更多关于分数阶累加的应用可查看相应论文：

A_novel_fractional_grey_system_model_and_its_application
A novel fractional time delayed grey model with Grey Wolf Optimizer and its applications in forecasting the natural gas and coal consumption in Chongqing China
The novel fractional discrete multivariate grey system model and its applications