0 引言

张量动力系统已广泛应用于Volterra系统识别^[1]、张量乘积(TP)模型变换^[2]、人体动作识别^[3]和塑性模型^[4]等各个领域.由于在张量常微分方程的求解过程中, 必须用到张量指数函数, 其计算问题已经成为一个重要的研究领域.

考虑如下常微分方程^[4]:

(1)

其中A和Y₀是给定的张量, 一般是非对称的, 则关于系统(1)的张量指数函数exp(·)有如下唯一解:

(2)

对于一般的张量A, 其指数函数通常表示为幂级数展开形式, 有

(3)

目前计算式(2)中张量指数函数的方法通常是将无穷级数(3)前若干项截断求和^[5], 但是该方法的计算精度和效率受限于其舍入方式和终止准则.截断法取前n_max项, 得到近似解, 即

(4)

选取的项数n_max受限于所需精度

(5)

显然, 截断法的精度与张量指数函数选取的项数n_max有关, 保留的项数越多, 精度越高, 但是需要进行的张量乘积次数也越多; 保留的项数越少, 计算量越少, 但是精度也越低.因此, 在一定程度上, 这将限制该算法的实用性.

为研究张量指数函数的计算问题, 本文定义了张量的一种广义逆, 并以此为基础构造了张量广义逆Padé逼近(GITPA)的一种可以递推计算的连分式算法, 其特点在于:在计算过程中, 不必用到张量的乘积, 也不需要计算张量的逆, 另外, 该方法对奇异张量也是适用的.目前, 关于张量的逆和广义逆的计算问题，还没有找到一种比较可行的计算方法.作为GITPA方法的一个重要应用, 后文将给出计算张量指数函数的数值实验, 以说明该连分式算法的有效性.

1 张量的一种广义逆 1.1 张量简介

张量是一种多维数组, 其中向量是一阶张量, 矩阵是二阶张量, 特别地, 一个N阶张量拥有N个下标, 是N个拥有独立坐标系的向量空间的外积(张量积).一个N阶n₁×n₂×...×n_N维张量可以表示为

文献[6]提出了张量的切片方法.如对一个三阶张量, 可以固定其中任意一个下标, 从而得到该张量的一种表示形式.设A=(a_i1i2i3)∈R^2×2×3, 固定其第3个下标, 则张量可以表示为

下面定义两个三阶张量的t-积, 可以通过递归自然地推广到高阶张量的t-积.

定义1(块循环矩阵)^[7]   设A∈R^l×p×n, 那么A的块循环矩阵定义为

定义一个展开算子unfold(·)^[7], 用以下方式将一个p×m×n的张量展开成一个pn×m的矩阵:

fold(·)^[7]是其逆算子, 它会将一个pn×m的矩阵转化成p×m×n的张量.因此有

定义2(张量t-积)^[7-8]   设A是一个l×p×n的张量, B是一个p×m×n的张量, 则t-积A*B是一个l×m×n的张量, 定义如下:

例1   设A, B∈R^2×2×3, 现固定其第3个下标, 分别产生

由定义2得到

类似地，可以定义更高阶张量t-积.

定义3(高阶张量t-积)^[9]   设张量A∈R^{n₁×n₂×...×n_p}, B∈R^{n₂×l×...×n_p}, 则高阶张量t-积A*B是一个n₁×l×...×n_p的p阶张量, 定义为

设A∈R^{n₁×n₂×...×n_p}, 张量的范数等于它所有元素平方和的平方根^[6], 即

(6)

这与矩阵的Frobenius-模是类似的.两个大小相同的张量A, B∈R^{n₁×n₂×...×n_p}的内积等于它们对应元素的乘积的和^[6], 即

(7)

于是有(A, A)=||A||².

1.2 一种张量广义逆

参考复数、向量和矩阵的倒数, 容易得到以下结论:若复数c∈C, 则cc^*=|c|², 1/c=c^-1=c^*/|c|²; 若向量v∈Cⁿ, 则vv^*=|v|², 1/v=v^-1=v^*/|v|^2[10]; 若矩阵A=a_ij, B=b_ij∈C^s×t, 则, 1/A= A^-1=A^*/||A||^2[11-12].因此, 下面张量的广义逆可以看作是上述向量的倒数和矩阵的倒数在张量上的推广.

定义4   设A, B∈ C^{n₁×n₂×...×n_p}, 内积为

若A·A^*=||A||², 则张量A的广义逆定义为

(8)

由定义3, 这里张量A的广义逆A_r^-1实质上是指张量A的倒数.由上面定义易证以下引理.

引理1   设A, B∈ C^{n₁×n₂×...×n_p}, A, B≠0且c∈R, c≠0, 则下列关系成立:

1)

2)

2 广义逆张量Padé逼近的连分式递推算法 2.1 广义逆张量Padé逼近的定义

设f(z)是给定的张量幂级数, 系数为张量, 即

(9)

其中(c_i^{(i₁i₂...i_p)})∈R^{n₁×n₂×...×n_p}包含了所有p阶n₁×n₂×...× n_p维的张量.

定义5   C^{n₁×n₂×...×n_p}[z]是一个p阶张量多项式的集合, 维数分别为n₁×n₂×...×n_p.一个张量多项式

是m阶的, 表示为∂A(z)=m, 若满足

且

定义6   阶数为[n/2k]型广义逆张量Padé逼近(GITPA)是张量的一个有理函数, 有

(10)

其中: P(z)为一个张量多项式, Q(z)为一个数量多项式, 满足以下条件:

1)

2)

3)

P(z)=(p_(i1i2...ip)(z))∈C^{n₁×n₂×...×n_p}, 整除条件2)表明分母数量多项式Q(z)能够整除张量分子多项式P(z)的范数的平方, 而P(z)的范数||P(z)||²由定义(6)给出, *表示取复共轭.

2.2 连分式算法格式

设张量幂级数为式(9), 利用张量广义逆(8), 递推地定义Thiele型张量连分式如下:

其中:张量B_i∈C^{n₁×n₂×... ×n_p}, z ∈R, H(z)的n项截断分式定义为

(11)

利用张量广义逆的定义, 下面给出一个用于计算机编程的计算截断连分式(9)的递推算法.

算法1   张量Thiele型连分式算法.

输入:给定张量幂级数系数C_i, 逼近阶n.

Step 1:计算

Step 2:计算

令i=2.

Step 3:计算

Step 4:若i < n, 执行Step 3.若i=n, 计算

Step 5:将计算得到的B_i(i=0, ..., n)代入逼近公式(11), 并输出R_n(z).

在算法1中, D表示张量函数对于z的导数.

下面证明恒等定理.该定理表明, 根据GITPA的定义5, 由算法1得到的连分式截断分式R_n(z)是一个阶数为[n/2k]型的GITPA, 其分母多项式的阶数分为两种情况, n=2k或n=2k+1.

定理1(恒等定理)   设利用算法1计算截断连分式R_n(z)的过程中没有出现分母为0的情况, 则满足

(12)

证明   通过归纳法证明.当n=0, k=0时, 有

当n=1, k=0时, 有

其中: ∂P₁=1, ∂Q₁=0, Q₁|||P₁||².可得

(13)

当n=2, k=1时, 有

其中.由式(13)得到.

令, 则有∂g₁=2.由此可知

其中Q₂(z)=g₁.可知∂P₂(z)=2, ∂Q₂(z)=2.又由

可知Q₂|||P₂||²成立.由此可得

(14)

设当n-1=2k或n-1=2k+1(k=0, 1, ...)时, 下式成立：

(15)

当n=2k(k=0, 1, ...)时, 可以写作

(16)

(17)

由式(15)可知

(18)

设式(18)中, 易知∂g_n-1=2k.将式(17)代入(16), 可得

(19)

其中Q_n(z)=g_n-1.由式(18)得到

(20)

因为

从而有

(21)

由式(19) ~ (21)可知, 当n=2k(k=0, 1, ...)时, 有

(22)

同理可以证明, 当n=2k+1(k=0, 1, ...)时, 有

(23)

综上, 定理1成立.

2.3 算法算例

例2   设张量A∈R^2×2×2, 张量幂级数展开式为

(24)

计算其[2/2]型GITPA.

根据算法1, 计算结果如下:

其中

由定义6, R₂(z)=P₂(z)/Q₂(z)是式(24)的[2/2]型GITPA, 满足以下条件:

1)

2) Q₂(z)|||P₂(z)||².其中

3) Q₂(z)exp(Az)-P₂(t)=O(t³).

3 数值实验

设将要计算的张量指数函数为

(25)

目前, 计算张量指数函数(25)通常使用的方法是级数截断法^[5], 为了便于数值比较和验算, 下面给出它的算法程序.

算法2^[5]   级数截断法.

输入:给定张量A、自变量z和误差限的值.

Step 1:初始化n=0和exp(A):=I.

Step 2:n:=n+1.

Step 3:计算和Aⁿ.

Step 4:将上一步计算结果加入, 有

Step 5:判断是否停止, 若, 则停止, 否则转至Step 2.

输出:exp(At).

由连分式算法1计算得到的逼近值记作A值, 由截断法(即算法2)取前30项计算得到的近似值为理论值, 记作E值.定义残差RES为张量逼近值和理论值差的模的平方, 定义如下：

其中模算子||·||由式(6)定义.

例3   设给定张量指数函数exp(Az), A中元素为, , , , 其余均为0.

现在计算该张量指数函数的[4/4]型GITRA, 其计算结果在表 1中记作A值.首先, 利用定义2, 即张量t-积将exp(Az)展开成幂级数, 得到

分别计算该张量下标为(1, 2, 1), (1, 2, 2), (2, 2, 1)和(2, 2, 2)在点0.1, 0.2, 0.3, 0.4, 0.5的A值和E值, 计算结果见表 1.

从表 1中张量指数函数的连分式算法得到的逼近值(即A值)与截断法得到的理论值(即E值)的残差来看, 本文给出的张量连分式算法是有效的.

例4   设给定张量指数函数exp(Az), 其中A为如下2×2×2×2四阶张量：

现在计算该张量指数函数的[4/4]型GITRA.首先利用定义3, 即高阶张量t-积将exp(Az)展开成幂级数, 再利用本文给出的算法1计算张量指数函数, 将其计算结果与由截断法(即算法2)取前30项计算得到值进行比较, 得到指数函数每个元素在点0.1, 0.3, 0.5的A值和E值, 计算结果见表 2.在表 2中, 张量A的元素表示为a_i1i2i3i4.

从表 2中本文算法得到的逼近值与截断法得到的理论值的残差来看, 本文给出的张量连分式算法对高阶张量也是有效的.

4 结论

本文提出的一种张量广义逆是一个实用的计算方法, 它是向量广义逆^[10]和矩阵广义逆^[11-15]在张量上的推广.在该张量广义逆的基础上, 得到了计算张量指数函数的张量Thiele型连分式算法.从计算张量指数函数的数值实验来看, 所给出的算法是有效的.由张量连分式算法Step 1 ~ Step 4利用张量广义逆分别计算B_i(i=1, 2, ..., n)可以看出, 当张量指数函数(25)中张量A维数较大时, 连分式算法的计算效果会更加明显.以后的研究工作考虑将所提出的张量连分式方法用于控制论中的模型简化问题.