矩阵指数¶

By Sam Relton

这个notebook简要介绍了矩阵指数，一个在应用数学诸多领域有着广泛应用的有趣的函数。在给出其定义之后，还介绍了一些应用和算法。
注：算法部分介绍了一种可用但相对错误率表现一般的基本算法，此外作者在文末列举了可供使用的软件和进一步阅读的资料，这些内容这里没有收录；微分方程实例后的进一步说明也略去了，可以参阅原文

定义¶

%matplotlib inline
import numpy as np
import scipy as sp
import scipy.linalg as la
import matplotlib.pyplot as plt

矩阵指数是相应的标量函数的直接推广。请记住$\exp(x)$的泰勒级数展开：
$\exp(x) = \sum_{k = 0}^\infty \frac{x^k}{k!}$
我们可以利用矩阵乘法将它应用到方阵：
$\exp(A) = \sum_{k = 0}^\infty \frac{A^k}{k!}$
矩阵指数的2个最有用的属性是：
$\exp((a + b)A) = \exp(aA)\exp(bA)$
且当两个矩阵A和B满足AB = BA时,
$\exp(A + B) = \exp(A)\exp(B)$

应用¶

矩阵指数在包括癌症研究、核燃耗方程、计算机图形学等领域有着广泛而多样的应用。我们将看两个应用实例：线性微分方程的求解和复杂网络中重要节点的发现。

微分方程¶

矩阵指数的一个大亮点是，能够直接给出微分方程的解： $y'(t) = A y(t),\qquad y(0) = y_0$
对于 $y, y_0 \in \mathbb{C}^n$ 且 $A \in \mathbb{C}^{n \times n}$
解为：
$y(t) = \exp(At) y_0$
知道这一点，我们可以直接地计算y(t)的解，而不必使用像欧拉方法那样的时间步进方法。让我们看一个简单的例子：
$y'(t) = \begin{bmatrix} 1 & -20 \\ 3 & 4 \\ \end{bmatrix} y(t), \qquad y(0) = \begin{bmatrix}1\\0\end{bmatrix},\qquad t \in [0, 1]$

import scipy.integrate as sint

A = np.array([[1, -20],
              [3, 4]], dtype=np.double)

def yprime(y, t):
    return sp.dot(A, y)

t = np.linspace(0, 1, 51)
yzero = np.array([1, 0])

# Solve equation using odeint
y = sint.odeint(yprime, yzero, t) # Calls the Fortran library odeint
plt.plot(y[:, 0], y[:, 1]) # Blue line is y(t) for t in [0, 1]

plt.xlim([-4, 9])
plt.ylim([-3, 5])

# Solve for 10 t values using the exponential
tvals = np.linspace(0, 1, 11)

for tval in tvals:
    yval = sp.dot(la.expm(A*tval), yzero)
    plt.plot(yval[0], yval[1], 'rs') # Red squares are evaluated via the matrix exponential

网络中的重要节点¶

矩阵指数可以用于寻找网络中最具影响的节点。所谓网络，包括一组节点，以及一组连接这些节点的边。我们使用邻接矩阵来描述一个网络。
若节点i，j之间有边，则$A_{ij} = 1$，否则$A_{ij} = 0$
我们考虑下面的矩阵：

import networkx as nx

Adj = np.array([[0, 1, 1, 0, 0, 0], 
                [1, 0, 1, 1, 1, 0], 
                [1, 1, 0, 0, 0, 0], 
                [0, 1, 0, 0, 0, 1], 
                [0, 1, 0, 0, 0, 0],
                [0, 0, 0, 1, 0, 0]])

G = nx.from_numpy_matrix(Adj)
nx.draw(G, node_color='y', node_size=1000)

衡量网络中节点重要性的一个常用指标是中心度。用$A^k_{ij}$表示想哦那个节点i到节点j的长度为k的路径条数。我们可以如下计算节点i的中心度：将所有对中心度的贡献值（长度不等的由节点i到自身的路径数）相加

$c(i) = \alpha_1 A_{ii} + \alpha_2 A^2_{ii} + \alpha_3 A^3_{ii} + \cdots$

这里的系数$\alpha_k$待定一般而言，我们假定和长的路径相比，短路径更重要，因此$\alpha_k \ge \alpha_{k+1}$. 已经有了许多计算系数$\alpha_k$的公式，如果取
$\alpha_k = \frac{1}{k!}$ 就有$c(i) = \exp(A)$

看我们之前给出的网络，我们直观上期望节点1就是那个最重要的节点:它的度数(degree)最大，而且与节点3, 4, 0, 2都相连。
现在我们就来计算每个节点的中心度，并对这些节点按照重要性递减的顺序排序：

centralities = np.diag(la.expm(np.array(Adj, dtype=np.double)))
nodeorder = np.argsort(centralities)[::-1]

print np.array([nodeorder, centralities[nodeorder]])

# Note: This is already built into networkx using the following command
# print nx.communicability_centrality_exp(G)

[[ 1.          0.          2.          3.          4.          5.        ]
 [ 4.44723536  2.86427609  2.86427609  2.36018456  1.71615913  1.59432922]]

正如我们所预期的，节点1是最重要的节点.接下来是节点0和2，它们的中心度相同，交换这两个节点，图没有发生变化。

矩阵指数

矩阵指数¶

定义¶

应用¶

微分方程¶

网络中的重要节点¶

results matching ""

No results matching ""