Week 1: Python intro
Week 2: Matrices, vectors, norms, ranks
Week 3: Linear systems, eigenvectors, eigenvalues
Week 4: Singular value decomposition + test + homework seminar
Week 5: Sparse & structured matrices
Week 6: Iterative methods, preconditioners, matrix functions
Week 7: Test week
Week 8: App Period
Today we will talk about matrix functions.
Book to read: Nick Higham, Functions of matrices
In [3]:
import numpy as np
import scipy.linalg as sla
A = 2 * np.eye(2)
A = np.array([[4, 4], [4, 4]])
sla.funm(A, lambda x: np.sqrt(x) )
Out[3]:
It is very easy to define a matrix polynomial as
$$
P(A) = \sum_{k=0}^n c_k A^k.
$$
Side-note: Hamilton-Cayley theorem states that $F(A) = 0$ where $F(\lambda) = \det(A - \lambda I)$.
We can define a function of the matrix by Taylor series:
$$
f(A) = \sum_{k=0}^{\infty} c_k A^k.
$$
The convergence is understood as the convergence in some matrix norm.
Example of such series is the Neumann series
$$
(I - F)^{-1} = \sum_{k=0}^{\infty} F^k,
$$
which is well defined for $\rho(F) < 1$.
The most well-known matrix function is matrix exponential. In the scalar case,
$$
e^x = 1 + x + \frac{x^2}{2} + \frac{x^3}{6} + \ldots = \sum_{k=0}^{\infty} \frac{x^k}{k!},
$$
and it directly translates to the matrix case:
$$
e^A = \sum_{k=0}^{\infty} \frac{A^k}{k!},
$$
the series that always converges (why?).
A lot of practical problems are reduced to a system of linear ODEs of the form
$$
\frac{dy}{dt} = Ay, \quad y(0) = y_0.
$$
The formal solution is given by $y(t) = e^{At} y_0$, so if we know
$e^{At}$ (or can compute matrix-by-vector product fast) there is a big gain over the
time-stepping schemes.
Let us look how matrices are diagonalizable:
In [21]:
import numpy as np
eps = 0.5
p = 4
a = np.eye(p)
for i in xrange(p-1):
a[i, i+1] = 1
a[p-1, 2] = eps
a = np.array(a)
val, vec = np.linalg.eig(a)
#print a
print np.linalg.norm(a - vec.dot(np.diag(val)).dot(np.linalg.inv(vec)))
#print 'S * D * S^{-1}:'
#print vec.dot(np.diag(val)).dot(np.linalg.inv(vec))
print a
Now we can compute a function
In [23]:
import numpy as np
eps = 1e-10
p = 5
a = np.eye(p)
for i in xrange(p-1):
a[i, i+1] = 1
a[p-1, 0] = eps
a = np.array(a)
val, vec = np.linalg.eig(a)
print np.linalg.norm(a - vec.dot(np.diag(val)).dot(np.linalg.inv(vec)))
fun = lambda x: np.exp(x)
#Using diagonalization
fun_diag = vec.dot(np.diag(fun(val))).dot(np.linalg.inv(vec))
#Using Schur
import scipy
fun_m = scipy.linalg.expm(a)
np.linalg.norm(fun_m - fun_diag)
Out[23]:
Therefore, $F(A)=U F(T) U^*$
Matrix exponential is well approximated by rational function:
$$\exp(x) \approx \frac{p(x)}{q(x)},$$where $p(x)$ and $q(x)$ are polynomials
and computation of a rational function of a matrix is reduced to matrix-matrix products and matrix inversions.
The rational form is also very useful when only a product of a matrix exponential by vector is needed, since
evaluation reduces to matrix-by-vector products and linear systems solvers
A short demo on Pade approximation via sympy package
In [31]:
import sympy.mpmath
%matplotlib inline
from sympy.mpmath import pade, taylor, polyval
import matplotlib.pyplot as plt
x = np.linspace(0, 1, 128)
a = taylor(sympy.mpmath.exp, 0, 100)
p, q = pade(a, 4, 2)
plt.plot(x, polyval(p[::-1], x)/polyval(q[::-1], x) - np.exp(x))
Out[31]:
One of the recent application is the computation. This example was taken from here. Take the following network.
In [32]:
import networkx as nx
import numpy as np
%matplotlib inline
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, with_labels=True )
One measure of the importance of each node is its centrality.
We can count the number of paths of different lengths from $i$ to $i$,
and add them up to get centrality of the node
$$c(i) = \alpha_1 A_{ii} + \alpha_2 A^2_{ii} + \alpha_3 A^3_{ii} + \cdots,$$
where the coefficients $\alpha_k$ remain to be chosen.
With $\alpha_k = \frac{1}{k!}$ we get
$$
c = \mathrm{diag}(\exp(A))
$$
In [33]:
centralities = np.diag(sla.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)
In [29]:
from IPython.core.display import HTML
def css_styling():
styles = open("./styles/custom.css", "r").read()
return HTML(styles)
css_styling()
Out[29]:
In [ ]: