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
Fabian Joswig 01ef97ffdf
docs: example 5 validation updated.
2023-07-20 12:20:06 +01:00

8 KiB
Raw Blame History

None <html> <head> </head>
In [1]:
import pyerrors as pe
import numpy as np

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)

[[4.10(20) -1.00(10)]
 [-1.00(10) 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)

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

Multiplication with unit matrix leaves the matrix unchanged

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

[[4.1 -1.0]
 [-1.0 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

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

Mathematical functions work elementwise

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

[[30.16185746098009 -1.1752011936438014]
 [-1.1752011936438014 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])
pe.gm(vector)
print(vector)

[2.00(40) 1.00(10)]

The matrix times vector product can then be computed via

In [8]:
product = matrix @ vector
pe.gm(product)
print(product)

[7.2(1.7) -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)
pe.gm(cholesky)
print(cholesky)

[[2.025(49) 0.0]
 [-0.494(51) 0.870(29)]]

We can now check if the decomposition was succesfull

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

[[-8.881784197001252e-16 0.0]
 [0.0 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)
pe.gm(inv)
print(inv)
print('Check:')
check_inv = cholesky @ inv
print(check_inv)

[[0.494(12) 0.0]
 [0.280(40) 1.150(39)]]
Check:
[[1.0 0.0]
 [0.0 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)
pe.gm(e)
print('Eigenvalues:')
print(e)
pe.gm(v)
print('Eigenvectors:')
print(v)

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

We can check that we got the correct result

In [14]:
print(matrix @ v == e * v)

[[ True  True]
 [ True  True]]
In [ ]:
</html>