mirror of
https://github.com/fjosw/pyerrors.git
synced 2025-03-15 14:50:25 +01:00
404 lines
8 KiB
Text
404 lines
8 KiB
Text
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 1,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"import pyerrors as pe\n",
|
|
"import numpy as np"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"As an example we look at a symmetric 2x2 matrix which positive semidefinte and has an error on all entries"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 2,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"[[4.10(20) -1.00(10)]\n",
|
|
" [-1.00(10) 1.000(10)]]\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"obs11 = pe.pseudo_Obs(4.1, 0.2, 'e1')\n",
|
|
"obs22 = pe.pseudo_Obs(1, 0.01, 'e1')\n",
|
|
"obs12 = pe.pseudo_Obs(-1, 0.1, 'e1')\n",
|
|
"matrix = np.asarray([[obs11, obs12], [obs12, obs22]])\n",
|
|
"print(matrix)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"We require to use `np.asarray` here as it makes sure that we end up with a numpy array of `Obs`."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"The standard matrix product can be performed with `@`"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 3,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"[[17.81 -5.1]\n",
|
|
" [-5.1 2.0]]\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"print(matrix @ matrix)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Multiplication with unit matrix leaves the matrix unchanged"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 4,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"[[4.1 -1.0]\n",
|
|
" [-1.0 1.0]]\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"print(matrix @ np.identity(2))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"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."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 5,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"[[78.12099999999998 -22.909999999999997]\n",
|
|
" [-22.909999999999997 7.1]]\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"print(pe.linalg.matmul(matrix, matrix, matrix)) # Equivalent to matrix @ matrix @ matrix but faster for large matrices"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Mathematical functions work elementwise"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 6,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"[[30.16185746098009 -1.1752011936438014]\n",
|
|
" [-1.1752011936438014 1.1752011936438014]]\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"print(np.sinh(matrix))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"For a vector of `Obs`, we again use `np.asarray` to end up with the correct object"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 7,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"[2.00(40) 1.00(10)]\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"vec1 = pe.pseudo_Obs(2, 0.4, 'e1')\n",
|
|
"vec2 = pe.pseudo_Obs(1, 0.1, 'e1')\n",
|
|
"vector = np.asarray([vec1, vec2])\n",
|
|
"pe.gm(vector)\n",
|
|
"print(vector)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"The matrix times vector product can then be computed via"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 8,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"[7.2(1.7) -1.00(46)]\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"product = matrix @ vector\n",
|
|
"pe.gm(product)\n",
|
|
"print(product)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"`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"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 9,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"3.10(28)\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"det = pe.linalg.det(matrix)\n",
|
|
"det.gamma_method()\n",
|
|
"print(det)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"The cholesky decomposition can be obtained as follows"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 10,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"[[2.025(49) 0.0]\n",
|
|
" [-0.494(51) 0.870(29)]]\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"cholesky = pe.linalg.cholesky(matrix)\n",
|
|
"pe.gm(cholesky)\n",
|
|
"print(cholesky)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"We can now check if the decomposition was succesfull"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 11,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"[[-8.881784197001252e-16 0.0]\n",
|
|
" [0.0 0.0]]\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"check = cholesky @ cholesky.T\n",
|
|
"print(check - matrix)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"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."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 12,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"[[0.494(12) 0.0]\n",
|
|
" [0.280(40) 1.150(39)]]\n",
|
|
"Check:\n",
|
|
"[[1.0 0.0]\n",
|
|
" [0.0 1.0]]\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"inv = pe.linalg.inv(cholesky)\n",
|
|
"pe.gm(inv)\n",
|
|
"print(inv)\n",
|
|
"print('Check:')\n",
|
|
"check_inv = cholesky @ inv\n",
|
|
"print(check_inv)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Eigenvalues and eigenvectors\n",
|
|
"We can also compute eigenvalues and eigenvectors of symmetric matrices with a special wrapper `eigh`"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 13,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Eigenvalues:\n",
|
|
"[0.705(57) 4.39(19)]\n",
|
|
"Eigenvectors:\n",
|
|
"[[-0.283(26) -0.9592(76)]\n",
|
|
" [-0.9592(76) 0.283(26)]]\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"e, v = pe.linalg.eigh(matrix)\n",
|
|
"pe.gm(e)\n",
|
|
"print('Eigenvalues:')\n",
|
|
"print(e)\n",
|
|
"pe.gm(v)\n",
|
|
"print('Eigenvectors:')\n",
|
|
"print(v)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"We can check that we got the correct result"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 14,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"[[ True True]\n",
|
|
" [ True True]]\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"print(matrix @ v == e * v)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": []
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"display_name": "Python 3 (ipykernel)",
|
|
"language": "python",
|
|
"name": "python3"
|
|
},
|
|
"language_info": {
|
|
"codemirror_mode": {
|
|
"name": "ipython",
|
|
"version": 3
|
|
},
|
|
"file_extension": ".py",
|
|
"mimetype": "text/x-python",
|
|
"name": "python",
|
|
"nbconvert_exporter": "python",
|
|
"pygments_lexer": "ipython3",
|
|
"version": "3.11.0"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 4
|
|
}
|