lattice/pyerrors
1
0
Fork 0
mirror of https://github.com/fjosw/pyerrors.git synced 2025-03-15 14:50:25 +01:00
Code Issues Projects Releases Packages Wiki Activity
pyerrors/examples/05_matrix_operations.ipynb

8.8 KiB
Raw Permalink Blame History

None <html> <head> </head>
In [1]:
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

In [2]:
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)

[[Obs[4.10(20)] Obs[-1.00(10)]]
 [Obs[-1.00(10)] Obs[1.000(10)]]]

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 @

In [3]:
print(matrix @ matrix)

[[Obs[17.81] Obs[-5.1]]
 [Obs[-5.1] Obs[2.0]]]

Multiplication with unit matrix leaves the matrix unchanged

In [4]:
print(matrix @ np.identity(2))

[[Obs[4.1] Obs[-1.0]]
 [Obs[-1.0] Obs[1.0]]]

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.

In [5]:
print(pe.linalg.matmul(matrix, matrix, matrix))  # Equivalent to matrix @ matrix @ matrix but faster for large matrices

[[Obs[78.12099999999998] Obs[-22.909999999999997]]
 [Obs[-22.909999999999997] Obs[7.1]]]

Mathematical functions work elementwise

In [6]:
print(np.sinh(matrix))

[[Obs[30.161857460980094] Obs[-1.1752011936438014]]
 [Obs[-1.1752011936438014] Obs[1.1752011936438014]]]

For a vector of Obs, we again use np.asarray to end up with the correct object

In [7]:
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)

[Obs[2.00(40)] Obs[1.00(10)]]

The matrix times vector product can then be computed via

In [8]:
product = matrix @ vector
for (i), entry in np.ndenumerate(product):
    entry.gamma_method()
print(product)

[Obs[7.2(1.7)] Obs[-1.00(46)]]

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

In [9]:
det = pe.linalg.det(matrix)
det.gamma_method()
print(det)

3.10(28)

The cholesky decomposition can be obtained as follows

In [10]:
cholesky = pe.linalg.cholesky(matrix)
for (i, j), entry in np.ndenumerate(cholesky):
    entry.gamma_method()
print(cholesky)

[[Obs[2.025(49)] Obs[0.0]]
 [Obs[-0.494(50)] Obs[0.870(29)]]]

We can now check if the decomposition was succesfull

In [11]:
check = cholesky @ cholesky.T
print(check - matrix)

[[Obs[-8.881784197001252e-16] Obs[0.0]]
 [Obs[0.0] Obs[0.0]]]

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.

In [12]:
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)

[[Obs[0.494(12)] Obs[0.0]]
 [Obs[0.280(40)] Obs[1.150(39)]]]
Check:
[[Obs[1.0] Obs[0.0]]
 [Obs[0.0] Obs[1.0]]]

Eigenvalues and eigenvectors

We can also compute eigenvalues and eigenvectors of symmetric matrices with a special wrapper eigh

In [13]:
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)

Eigenvalues:
[Obs[0.705(57)] Obs[4.39(19)]]
Eigenvectors:
[[Obs[-0.283(26)] Obs[-0.9592(75)]]
 [Obs[-0.9592(75)] Obs[0.283(26)]]]

We can check that we got the correct result

In [14]:
for i in range(2):
    print('Check eigenvector', i + 1)
    print(matrix @ v[:, i] - v[:, i] * e[i])

Check eigenvector 1
[Obs[-5.551115123125783e-17] Obs[0.0]]
Check eigenvector 2
[Obs[0.0] Obs[-2.220446049250313e-16]]
In [ ]:
</html>