8.8 KiB
import pyerrors as pe
import numpy as np
import scipy
As an example we look at a symmetric 2x2 matrix which positive semidefinte and has an error on all entries
obs11 = pe.pseudo_Obs(4.1, 0.2, 'e1')
obs22 = pe.pseudo_Obs(1, 0.01, 'e1')
obs12 = pe.pseudo_Obs(-1, 0.1, 'e1')
matrix = np.asarray([[obs11, obs12], [obs12, obs22]])
print(matrix)
We require to use np.asarray
here as it makes sure that we end up with a numpy array of Obs
.
The standard matrix product can be performed with @
print(matrix @ matrix)
Multiplication with unit matrix leaves the matrix unchanged
print(matrix @ np.identity(2))
For large matrices overloading the standard operator @
can become inefficient as pyerrors has to perform a large number of elementary opeations. For these situations pyerrors provides the function linalg.matmul
which optimizes the required automatic differentiation. The function can take an arbitray number of operands.
print(pe.linalg.matmul(matrix, matrix, matrix)) # Equivalent to matrix @ matrix @ matrix but faster for large matrices
Mathematical functions work elementwise
print(np.sinh(matrix))
For a vector of Obs
, we again use np.asarray
to end up with the correct object
vec1 = pe.pseudo_Obs(2, 0.4, 'e1')
vec2 = pe.pseudo_Obs(1, 0.1, 'e1')
vector = np.asarray([vec1, vec2])
for (i), entry in np.ndenumerate(vector):
entry.gamma_method()
print(vector)
The matrix times vector product can then be computed via
product = matrix @ vector
for (i), entry in np.ndenumerate(product):
entry.gamma_method()
print(product)
pyerrors
provides the user with wrappers to the numpy.linalg
functions which work on Obs
valued matrices. We can for example calculate the determinant of the matrix via
det = pe.linalg.det(matrix)
det.gamma_method()
print(det)
The cholesky decomposition can be obtained as follows
cholesky = pe.linalg.cholesky(matrix)
for (i, j), entry in np.ndenumerate(cholesky):
entry.gamma_method()
print(cholesky)
We can now check if the decomposition was succesfull
check = cholesky @ cholesky.T
print(check - matrix)
We can now further compute the inverse of the cholesky decomposed matrix and check that the product with its inverse gives the unit matrix with zero error.
inv = pe.linalg.inv(cholesky)
for (i, j), entry in np.ndenumerate(inv):
entry.gamma_method()
print(inv)
print('Check:')
check_inv = cholesky @ inv
print(check_inv)
Eigenvalues and eigenvectors¶
We can also compute eigenvalues and eigenvectors of symmetric matrices with a special wrapper eigh
e, v = pe.linalg.eigh(matrix)
for (i), entry in np.ndenumerate(e):
entry.gamma_method()
print('Eigenvalues:')
print(e)
for (i, j), entry in np.ndenumerate(v):
entry.gamma_method()
print('Eigenvectors:')
print(v)
We can check that we got the correct result
for i in range(2):
print('Check eigenvector', i + 1)
print(matrix @ v[:, i] - v[:, i] * e[i])