Provides classes for calculating with krylov spaces and matrix functions
Project description
matrix-functions
The krylov space is a vector space generated by a square matrix and the Cayley-Hamilton theorem which states that every square matrix is a root of its own characteristic polynomial. From this it follows that $A^n$ is a linear combination of $\unicode{x1D7D9}=A^0,\ A,\ A^2,\ \dots,\ A^{n-1}$ and therefore every polynomial in A is such a linear combination:
$$ A^m = \sum_{k=0}^{n-1} \alpha_{mk} A^k $$
It turns out that $\alpha_{mk}$ only depends on the eigenvalues $\lambda_1,\ \dots,\ \lambda_n$. Because of this every matrix function can be expressed in the krylov space if the function is analytical in the eigenvalues:
$$ A^m = \sum_{k=0}^{n-1} A^k \sum_{j=1}^{k+1} \Lambda_{k+1-j} \sum_{l=1}^r \sum_{p=0}^{\min(\mu_l-1,m-j)} \bar\beta_{lp} \lambda_l^{m-j-p} \frac{(m-j)!}{(m-j-p)!} $$
$$ f(A) = \sum_{k=0}^{n-1} A^k \sum_{j=1}^{k+1} \Lambda_{k+1-j} \sum_{l=1}^r \sum_{p=0}^{\mu_l-1} \bar\beta_{lp} \sum_{q=0}^p \binom pq (-1)^{p-q}\frac{(j-1+p-q)!}{(j-1)!} \lambda_l^{-j-p+q} f^{(q)}(\lambda_l) $$
Example
# One can show that sin(x) = 2cos(a)sin(x-a) - sin(x-2a). So for `x=0` and `a=.2` we can define a
# matrix `M` such that [sin(x+na), sin(x+(n-1)a)] = matrix_power(M, n) @ [sin(x), sin(x-a)]:
import numpy as np
from matrixfuncs import *
import matplotlib.pyplot as plt
a = .2
M = np.array([[2*np.cos(a), -1], [1,0]])
for n in [0, 1, 10, 100, 1000, 10000]:
sinNa = [1,0] @ np.linalg.matrix_power(M, n) @ np.sin([0,-a])
error = (sinNa - np.sin(n*a))**2
print('n:', error)
# n: 0.0
# n: 0.0
# n: 4.930380657631324e-30
# n: 1.199857436841159e-27
# n: 9.247078067402441e-26
# n: 3.397767010759534e-24
# In order to generalize the discrete n to a continous variable we have to replace
# matrix_power(M, n) with something like exp(n*log(M)). For calculating the function of a matrix
# we can use MFunc.
krylovM = KrylovSpace(M)
mfn = MFunc.fromKrylov(krylovM)
mfn0 = [1,0] @ mfn @ np.sin([0,-a])
# Compare the numeric evaluation of sin to numpy's implementation in 1000 points between 0 and 2pi:
ts = np.linspace(0, 2*np.pi, 1000)
sinFromMFunc = mfn0(lambda evs: np.exp(np.outer(np.log(evs), ts/a))).coeffs
sints = np.sin(ts)
error = np.sum((sinFromMFunc - sints)**2)
print(error)
# 1.1051029979680703e-26
plt.plot(ts, sinFromMFunc)
plt.savefig('sinFromMFunc.png', dpi=200)
Project details
Release history Release notifications | RSS feed
Download files
Download the file for your platform. If you're not sure which to choose, learn more about installing packages.
Source Distribution
Built Distribution
Hashes for matrixfuncs-0.1.0-py3-none-any.whl
Algorithm | Hash digest | |
---|---|---|
SHA256 | 834524375fd5962a2720fb902347242348105a75d4eb9ea942e5b3a3b9c64902 |
|
MD5 | aa920bb26ee4746c054b3062c95b4b69 |
|
BLAKE2b-256 | 2b8f25bac7d50eb1a5aa0b41e670d929b19e951f97f36f9eadc7b3d5e6e369f2 |