Skip to content

Instantly share code, notes, and snippets.

@pelagisk

pelagisk/Matrisfaktorisering.ipynb

Created November 15, 2017 10:25
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save pelagisk/1b0ef0d77ab8c0ec98a4fc2079a45557 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save pelagisk/1b0ef0d77ab8c0ec98a4fc2079a45557 to your computer and use it in GitHub Desktop.
Download ZIP
Notebook for teaching about matrix factorisation, mostly in Swedish!
Display the source blob
Display the rendered blob
Raw
Matrisfaktorisering.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Matrisfaktorisering\n",
"\n",
"Matrisfaktorisering är ett samlingsnamn på metoder som \"delar upp\" matriser i matrisprodukter. Vi kommer diskutera två fall, LU-faktorisering och QR-faktorisering. Det finns många fler: polär dekomposition, singulärvärdes-dekomposition (SVD), Cholesky-faktorisering och tensorfaktoriseringar m.m....\n",
"\n",
"### Varför är det viktigt? Och varför behöver vi bry oss?\n",
"\n",
"- Det är ett problem som är lätt att formulera med för tråkigt/långsamt för människor att göra för hand. \n",
"- Nästan alla problem i naturvetenskap har att göra med matriser! Det är vanligt att behöva jobba med matriser av dimension 100.000 x 100.000. Utan matrisfaktoriseringar skulle det vara omöjligt att lösa såna problem.\n",
"- Matrisfaktoriseringar används överallt! Datorgrafik. När Netflix presenterar de mest populära filmerna för dig så är det matrisfaktorisering. [Det blir också allt viktigare för kvantfysik.](https://en.wikipedia.org/wiki/Matrix_product_state)\n",
"- Det är ett bra, lagom svårt problem att öva på! "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Linjära ekvationssystem\n",
"\n",
"Ett *system av linjära ekvationer*, t.ex. \n",
"\n",
"$2w + u = -4$\n",
"\n",
"$3w + 4u = 3$\n",
"\n",
"kan alltid skrivas om som en ekvation för en okänd vektor $\\mathbf{x}=(w,u)^T$:\n",
"\n",
"$\\left( {\\begin{array}{cc}\n",
" 2 & 1 \\\\\n",
" 3 & 4 \\\\\n",
" \\end{array} } \\right) \\left( {\\begin{array}{c}\n",
" w \\\\\n",
" u \\\\\n",
" \\end{array} } \\right) = \\left( {\\begin{array}{c}\n",
" -4 \\\\\n",
" 3 \\\\\n",
" \\end{array} } \\right) \\Leftrightarrow \\mathbf{A} \\mathbf{x} = \\mathbf{b}$\n",
"\n",
"Ni har lärt er att lösa linjära ekvationssystem med Gauss-eliminering. I praktiken är det opraktiskt att utföra för stora matriser, men det är lätt att låta en dator utföra beräkningen. Följande kod är baserad på exempel 6.1 i Newmans bok (sida 218):"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"from copy import copy\n",
"import numpy as np\n",
"\n",
"\n",
"def gaussian_elimination(A, b):\n",
" \"\"\"Solves the matrix equation A x = b for x\"\"\"\n",
" N, _ = A.shape\n",
" U = copy(A)\n",
" c = copy(b)\n",
" for i in range(N):\n",
" \n",
" # divide by diagonal element\n",
" div = U[i, i]\n",
" U[i, :] /= div\n",
" c[i] /= div\n",
"\n",
" # subtract from lower rows\n",
" for j in range(i+1, N):\n",
" mult = U[j, i]\n",
" U[j, :] -= mult * U[i, :]\n",
" c[j] -= mult * c[i]\n",
"\n",
" # backsubstitution\n",
" x = np.empty(N, float)\n",
" for i in range(N-1, -1, -1):\n",
" x[i] = c[i]\n",
" for j in range(i+1, N):\n",
" x[i] -= U[i, j] * x[j]\n",
"\n",
" return x"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Det finns redan en funktion för att lösa linjära ekvationssystem i NumPy, som helt enkelt heter `solve`. För att testa vår kod jämför vi med `numpy.linalg.solve`:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A = \n",
"[[ 2. 1.]\n",
" [ 3. 4.]]\n",
"b = \n",
"[-4. 3.]\n",
"---- Using Newman's gaussian elimination:\n",
"x = [-3.8 3.6]\n",
"A*x = [-4. 3.]\n",
"max |b-A*x| = 1.7763568394e-15\n",
"---- Using NumPy's solve:\n",
"x = [-3.8 3.6]\n",
"A*x = [-4. 3.]\n",
"max |b-A*x| = 8.881784197e-16\n"
]
}
],
"source": [
"A = np.array([[2., 1.], \n",
" [3., 4.]])\n",
"b = np.array([-4., 3.])\n",
"\n",
"\n",
"print(\"A = \", A, sep=\"\\n\")\n",
"print(\"b = \", b, sep=\"\\n\")\n",
"\n",
"print(\"---- Using Newman's gaussian elimination:\")\n",
"\n",
"x = gaussian_elimination(A, b)\n",
"print(\"x = \", x)\n",
"print(\"A*x = \", np.dot(A, x))\n",
"print(\"max |b-A*x| = \", np.max(np.abs(b - np.dot(A, x))))\n",
"\n",
"print(\"---- Using NumPy's solve:\")\n",
"\n",
"x = np.linalg.solve(A, b)\n",
"\n",
"print(\"x = \", x)\n",
"print(\"A*x = \", np.dot(A, x))\n",
"print(\"max |b-A*x| = \", np.max(np.abs(b - np.dot(A, x))))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# LU-faktorisering\n",
"\n",
"En matrisfaktorisering $\\mathbf{A} = \\mathbf{L} \\mathbf{U}$ som har att göra med Gauss-eliminering. $\\mathbf{L}$ är en matris som är undertriangulär (bara nollor över diagonalen) och $\\mathbf{U}$ är övertriangulär (bara nollor under diagonalen). \n",
"\n",
"Antag att vi redan gjort en sån faktorisering: då är det lätt att lösa ett linjärt ekvationssystem som:\n",
"\n",
"$\\mathbf{A} \\mathbf{x} = \\mathbf{v} \\Rightarrow \\mathbf{L} \\mathbf{U} \\mathbf{x} = \\mathbf{v}$\n",
"\n",
"för alla tänkbara vektorer $\\mathbf{v}$! Börja med att räkna ut $\\mathbf{U} \\mathbf{x} = \\mathbf{y}$. Resultatet är en ekvation $\\mathbf{L} \\mathbf{y} = \\mathbf{v}$ där $\\mathbf{L}$ är i nedåttriangulär form. Det är enkelt att lösa en sån ekvation genom att börja med den första ekvationen (som är en linjär ekvation i **en** variabel eftersom $\\mathbf{L}$ är nedåttriangulär), substituera i den andra och så vidare precis som i delen \"backsubstitution\" i funktionen `gaussian_elimination` gör! \n",
"\n",
"Koden i `gaussian_elimination` kan återanvändas för att göra en LU-faktorisering med några få ändringar."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A = \n",
"[[ 0.628 0.376 0.612 0.534 0.186 0.328 0.542 0.251 0.784 0.326]\n",
" [ 0.716 0.672 0.974 0.15 0.642 0.371 0.608 0.452 0.648 0.194]\n",
" [ 0.003 0.968 0.353 0.029 0.473 0.785 0.795 0.064 0.976 0.715]\n",
" [ 0.068 0.125 0.989 0.53 0.721 0.153 0.107 0.187 0.568 0.316]\n",
" [ 0.528 0.882 0.895 0.72 0.215 0.402 0.053 0.846 0.925 0.839]\n",
" [ 0.539 0.661 0.785 0.25 0.813 0.593 0.746 0.035 0.863 0.051]\n",
" [ 0.947 0.924 0.382 0.192 0.194 0.951 0.558 0.2 0.51 0.012]\n",
" [ 0.565 0.489 0.988 0.866 0.39 0.86 0.249 0.292 0.469 0.779]\n",
" [ 0.026 0.327 0.898 0.953 0.77 0.329 0.645 0.929 0.123 0.273]\n",
" [ 0.694 0.627 0.25 0.987 0.665 0.998 0.924 0.349 0.086 0.807]]\n",
"---- Using my own LU algorithm:\n",
"L = \n",
"[[ 0.628 0. 0. 0. 0. 0. 0. 0. 0. 0. ]\n",
" [ 0.716 0.244 0. 0. 0. 0. 0. 0. 0. 0. ]\n",
" [ 0.003 0.966 -0.747 0. 0. 0. 0. 0. 0. 0. ]\n",
" [ 0.068 0.084 0.827 2.672 0. 0. 0. 0. 0. 0. ]\n",
" [ 0.528 0.566 -0.262 0.69 -0.299 0. 0. 0. 0. 0. ]\n",
" [ 0.539 0.338 -0.124 0.122 0.298 -0.264 0. 0. 0. 0. ]\n",
" [ 0.947 0.357 -0.946 -2.276 0.153 0.099 -1.24 0. 0. 0. ]\n",
" [ 0.565 0.15 0.266 1.326 -0.079 0.46 -1.577 0.056 0. 0. ]\n",
" [ 0.026 0.311 0.519 2.796 0.208 -0.456 0.945 0.968 2.836 0. ]\n",
" [ 0.694 0.211 -0.666 -0.85 0.928 -1.009 0.12 0.978 2.79 1.32 ]]\n",
"U = \n",
"[[ 1. 0.599 0.975 0.85 0.296 0.523 0.862 0.4 1.248 0.518]\n",
" [ 0. 1. 1.135 -1.881 1.766 -0.015 -0.037 0.678 -1.009 -0.725]\n",
" [-0. -0. 1. -2.468 1.651 -1.069 -1.109 0.793 -2.608 -1.894]\n",
" [ 0. 0. 0. 1. -0.305 0.375 0.363 -0.207 1.02 0.714]\n",
" [-0. -0. -0. -0. 1. 1.351 3.081 -2.012 1.834 0.041]\n",
" [-0. -0. -0. -0. -0. 1. 3.048 -1.186 1.748 1.198]\n",
" [-0. -0. -0. -0. -0. -0. 1. -0.227 0.734 0.412]\n",
" [ 0. 0. 0. 0. 0. 0. 0. 1. -4.349 4.533]\n",
" [ 0. 0. 0. 0. 0. 0. 0. 0. 1. -1.68 ]\n",
" [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. ]]\n",
"L*U = \n",
"[[ 0.628 0.376 0.612 0.534 0.186 0.328 0.542 0.251 0.784 0.326]\n",
" [ 0.716 0.672 0.974 0.15 0.642 0.371 0.608 0.452 0.648 0.194]\n",
" [ 0.003 0.968 0.353 0.029 0.473 0.785 0.795 0.064 0.976 0.715]\n",
" [ 0.068 0.125 0.989 0.53 0.721 0.153 0.107 0.187 0.568 0.316]\n",
" [ 0.528 0.882 0.895 0.72 0.215 0.402 0.053 0.846 0.925 0.839]\n",
" [ 0.539 0.661 0.785 0.25 0.813 0.593 0.746 0.035 0.863 0.051]\n",
" [ 0.947 0.924 0.382 0.192 0.194 0.951 0.558 0.2 0.51 0.012]\n",
" [ 0.565 0.489 0.988 0.866 0.39 0.86 0.249 0.292 0.469 0.779]\n",
" [ 0.026 0.327 0.898 0.953 0.77 0.329 0.645 0.929 0.123 0.273]\n",
" [ 0.694 0.627 0.25 0.987 0.665 0.998 0.924 0.349 0.086 0.807]]\n",
"max |A-L*U| = 6.66133814775e-16\n",
"---- Using SciPy LU\n",
"P = \n",
"[[ 0. 0. 0. 0. 0. 0. 1. 0. 0. 0.]\n",
" [ 0. 0. 0. 0. 0. 0. 0. 0. 1. 0.]\n",
" [ 0. 1. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
" [ 0. 0. 1. 0. 0. 0. 0. 0. 0. 0.]\n",
" [ 0. 0. 0. 0. 1. 0. 0. 0. 0. 0.]\n",
" [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 1.]\n",
" [ 1. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
" [ 0. 0. 0. 0. 0. 1. 0. 0. 0. 0.]\n",
" [ 0. 0. 0. 0. 0. 0. 0. 1. 0. 0.]\n",
" [ 0. 0. 0. 1. 0. 0. 0. 0. 0. 0.]]\n",
"L = \n",
"[[ 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]\n",
" [ 0.003 1. 0. 0. 0. 0. 0. 0. 0. 0. ]\n",
" [ 0.072 0.06 1. 0. 0. 0. 0. 0. 0. 0. ]\n",
" [ 0.733 -0.052 -0.012 1. 0. 0. 0. 0. 0. 0. ]\n",
" [ 0.557 0.381 0.583 0.353 1. 0. 0. 0. 0. 0. ]\n",
" [ 0.597 -0.065 0.832 0.381 0.708 1. 0. 0. 0. 0. ]\n",
" [ 0.663 -0.246 0.473 0.199 0.388 0.042 1. 0. 0. 0. ]\n",
" [ 0.028 0.312 0.827 0.601 0.416 0.101 0.606 1. 0. 0. ]\n",
" [ 0.756 -0.027 0.739 -0.438 -0.376 -0.717 0.661 0.504 1. 0. ]\n",
" [ 0.569 0.14 0.551 -0.173 -0.538 -0.556 0.278 -0.072 0.515 1. ]]\n",
"U = \n",
"[[ 0.947 0.924 0.382 0.192 0.194 0.951 0.558 0.2 0.51 0.012]\n",
" [ 0. 0.965 0.351 0.028 0.473 0.782 0.793 0.063 0.975 0.715]\n",
" [ 0. 0. 0.94 0.515 0.678 0.037 0.019 0.168 0.472 0.272]\n",
" [ 0. 0. 0. 0.854 0.556 0.342 0.557 0.208 -0.232 0.84 ]\n",
" [ 0. 0. 0. 0. -0.665 -0.568 -0.768 0.539 0.076 0.105]\n",
" [ 0. 0. 0. 0. 0. 0.584 0.283 -0.424 -0.131 0.198]\n",
" [ 0. 0. 0. 0. 0. 0. 0.533 -0.178 0.484 0.149]\n",
" [ 0. 0. 0. 0. 0. 0. 0. 0.566 -0.758 -0.833]\n",
" [ 0. 0. 0. 0. 0. 0. 0. 0. -0.164 0.876]\n",
" [ 0. 0. 0. 0. 0. 0. 0. 0. 0. -0.446]]\n",
"L*U = \n",
"[[ 0.628 0.376 0.612 0.534 0.186 0.328 0.542 0.251 0.784 0.326]\n",
" [ 0.716 0.672 0.974 0.15 0.642 0.371 0.608 0.452 0.648 0.194]\n",
" [ 0.003 0.968 0.353 0.029 0.473 0.785 0.795 0.064 0.976 0.715]\n",
" [ 0.068 0.125 0.989 0.53 0.721 0.153 0.107 0.187 0.568 0.316]\n",
" [ 0.528 0.882 0.895 0.72 0.215 0.402 0.053 0.846 0.925 0.839]\n",
" [ 0.539 0.661 0.785 0.25 0.813 0.593 0.746 0.035 0.863 0.051]\n",
" [ 0.947 0.924 0.382 0.192 0.194 0.951 0.558 0.2 0.51 0.012]\n",
" [ 0.565 0.489 0.988 0.866 0.39 0.86 0.249 0.292 0.469 0.779]\n",
" [ 0.026 0.327 0.898 0.953 0.77 0.329 0.645 0.929 0.123 0.273]\n",
" [ 0.694 0.627 0.25 0.987 0.665 0.998 0.924 0.349 0.086 0.807]]\n",
"max |A-L*U| = 2.22044604925e-16\n"
]
}
],
"source": [
"import numpy as np\n",
"from scipy.linalg import lu as scipy_lu\n",
"from luqr import axels_lu as my_lu\n",
"\n",
"\n",
"np.set_printoptions(precision=3) # only show 3 decimal places in the matrices\n",
"\n",
"# start from a random matrix\n",
"A = np.random.rand(10, 10)\n",
"print(\"A = \", A, sep=\"\\n\")\n",
"\n",
"print(\"---- Using my own LU algorithm:\")\n",
"\n",
"L, U = my_lu(A)\n",
"print(\"L = \", L, \"U = \", U, \"L*U = \", L.dot(U), sep=\"\\n\") \n",
"print(\"max |A-L*U| =\", np.max(np.abs(A-L.dot(U)))) \n",
"\n",
"print(\"---- Using SciPy LU\")\n",
"\n",
"P, L, U = scipy_lu(A)\n",
"print(\"P = \", P, \"L = \", L, \"U = \", U, \"L*U = \", P.dot(L).dot(U), sep=\"\\n\") \n",
"print(\"max |A-L*U| =\", np.max(np.abs(A-P.dot(L).dot(U))))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Parentes: pivoting\n",
"\n",
"Säg att vi ville göra en LU-faktorisering av\n",
"\n",
"$\\left( {\\begin{array}{cc}\n",
" 0 & 1 \\\\\n",
" 3 & 4 \\\\\n",
" \\end{array} } \\right)$\n",
" \n",
"Med vår algoritm skulle vi direkt få problem eftersom vi delar element med diagonalelementen. En enkel lösning på problemet är att byta plats på raderna och hitta LU-faktoriseringen av den, och komma ihåg att byta tillbaka raderna när vi använder faktoriseringen senare.\n",
"\n",
"SciPys funktion `scipy.linalg.lu` gör precis det: den returnerar en tupel `(P, L, U)`. Det första värdet är en **permutationsmatris**, en matris av ettor och nollor som är till för att byta plats på raderna i produkten `L.dot(U)`. En bra hemuppgift är att programmera en egen algoritm för pivoting enligt sida 221 i Newmans bok."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# QR-faktorisering\n",
"\n",
"$\\mathbf{A} = \\mathbf{Q} \\mathbf{R}$ där $\\mathbf{Q}$ är en *ortogonal* matris (det betyder att $\\mathbf{Q} \\mathbf{Q}^T = \\mathbf{Q}^T \\mathbf{Q} = \\mathbb{1}$) och $\\mathbf{R}$ är övertriangulär. $\\mathbf{A}$ måste vara en kvadratisk matris! Den här faktoriseringen används för att lösa egenvärdesproblem och för att anpassa kurvor till mätdata med mera.\n",
"\n",
"QR-algoritmen är så här:\n",
"\n",
"1. För en matris $\\mathbf{A}$, låt oss skriva $\\mathbf{A}_{:,i} = \\mathbf{a}_i$ för den i:te kolumnen i matrisen.\n",
"\n",
"2. Beräkna vektorerna\n",
"\n",
" $\\mathbf{u}_i = \\mathbf{a}_i - \\sum_{j=0}^{i-1} (\\mathbf{q}_i \\cdot \\mathbf{a}_i) \\mathbf{q_j}$ \n",
"\n",
" och \n",
"\n",
" $\\mathbf{q}_i = \\frac{\\mathbf{u}_i}{\\left|{\\mathbf{u_i}}\\right|}$\n",
"\n",
"3. Matrisen $\\mathbf{Q}$ har $\\mathbf{q}_i$ som sina kolumner.\n",
"\n",
"4. Skapa en matris $\\mathbf{R}$ som ser ut så här:\n",
"\n",
" $\\mathbf{R} = \\left( {\\begin{array}{cccc}\n",
" \\left| \\mathbf{u}_0 \\right| & \\mathbf{q}_0 \\cdot \\mathbf{a}_1 & \\mathbf{q}_0 \\cdot \\mathbf{a}_2 & \\dots \\\\\n",
" 0 & \\left| \\mathbf{u}_1 \\right| & \\mathbf{q}_1 \\cdot \\mathbf{a}_2 & \\dots \\\\\n",
" 0 & 0 & \\left| \\mathbf{u}_2 \\right| & \\dots \\\\\n",
" \\dots & \\dots & \\dots & \\dots \\\\\n",
" \\end{array} } \\right)$\n",
" \n",
" eller formulerat på ett annat sätt:\n",
"\n",
" $\\mathbf{R}_{ij} = \\begin{cases}\\mathbf{q}_0 \\cdot \\mathbf{a}_1 \\qquad (i > j) \\\\ \\left| \\mathbf{u}_i \\right| \\qquad\\quad (i = j) \\\\ 0 \\qquad\\qquad (i < j) \\end{cases}$ \n",
"\n",
"5. Returnera $(\\mathbf{Q}, \\mathbf{R})$"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A = \n",
"[[ 0.349 0.378 0.495 0.175 0.786 0.258 0.763 0.775 0.752 0.751]\n",
" [ 0.954 0.682 0.518 0.044 0.759 0.143 0.954 0.603 0.669 0.325]\n",
" [ 0.295 0.905 0.902 0.654 0.096 0.063 0.457 0.013 0.35 0.148]\n",
" [ 0.606 0.112 0.05 0.465 0.154 0.453 0.039 0.065 0.272 0.464]\n",
" [ 0.228 0.109 0.433 0.646 0.834 0.048 0.604 0.188 0.047 0.614]\n",
" [ 0.31 0.154 0.521 0.295 0.152 0.667 0.22 0.144 0.095 0.914]\n",
" [ 0.78 0.359 0.856 0.104 0.639 0.791 0.238 0.397 0.514 0.264]\n",
" [ 0.955 0.228 0.292 0.306 0.136 0.585 0.874 0.132 0.874 0.16 ]\n",
" [ 0.022 0.427 0.906 0.626 0.032 0.585 0.43 0.901 0.422 0.984]\n",
" [ 0.559 0.271 0.36 0.543 0.95 0.307 0.604 0.42 0.193 0.082]]\n",
"---- Using my own QR algorithm:\n",
"Q = \n",
"[[ 0.187 0.2 0.039 -0.134 0.442 0.424 0.372 -0.117 0.573 0.223]\n",
" [ 0.512 0.18 -0.419 -0.324 0.196 -0.113 0.11 0.308 -0.451 0.251]\n",
" [ 0.158 0.777 -0.205 0.189 -0.23 -0.182 -0.245 -0.343 0.161 -0.064]\n",
" [ 0.326 -0.219 -0.151 0.474 -0.231 0.353 -0.41 0.289 0.192 0.363]\n",
" [ 0.123 -0.013 0.356 0.449 0.427 -0.608 0.037 0.022 0.022 0.319]\n",
" [ 0.166 -0.012 0.393 -0.025 -0.18 0.293 0.053 -0.566 -0.466 0.398]\n",
" [ 0.419 -0.059 0.472 -0.5 0.003 -0.118 -0.501 0.057 0.224 -0.175]\n",
" [ 0.513 -0.292 -0.05 0.104 -0.437 -0.24 0.521 -0.15 0.2 -0.236]\n",
" [ 0.012 0.432 0.502 0.137 -0.249 0.188 0.295 0.573 -0.134 -0.112]\n",
" [ 0.3 -0.028 0.008 0.367 0.438 0.295 -0.092 -0.126 -0.281 -0.626]]\n",
"R = \n",
"[[ 1.862 0.993 1.284 0.81 1.352 1.127 1.562 0.895 1.376 0.938]\n",
" [ 0. 0.962 1.115 0.596 0.232 0.044 0.559 0.57 0.372 0.563]\n",
" [ 0. 0. 0.815 0.482 0.344 0.787 0.121 0.525 0.102 0.981]\n",
" [ 0. 0. 0. 0.854 0.158 0.008 0.212 -0.093 -0.141 0.344]\n",
" [ 0. 0. 0. 0. 1.119 -0.341 0.406 0.401 -0.077 0.074]\n",
" [ 0. 0. 0. 0. 0. 0.375 -0.136 0.424 0.14 0.453]\n",
" [ 0. 0. 0. 0. 0. 0. 0.702 0.435 0.467 0.393]\n",
" [ 0. 0. 0. 0. 0. 0. 0. 0.498 0.139 0.137]\n",
" [ 0. 0. 0. 0. 0. 0. 0. 0. 0.373 -0.079]\n",
" [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.723]]\n",
"Q*R = \n",
"[[ 0.349 0.378 0.495 0.175 0.786 0.258 0.763 0.775 0.752 0.751]\n",
" [ 0.954 0.682 0.518 0.044 0.759 0.143 0.954 0.603 0.669 0.325]\n",
" [ 0.295 0.905 0.902 0.654 0.096 0.063 0.457 0.013 0.35 0.148]\n",
" [ 0.606 0.112 0.05 0.465 0.154 0.453 0.039 0.065 0.272 0.464]\n",
" [ 0.228 0.109 0.433 0.646 0.834 0.048 0.604 0.188 0.047 0.614]\n",
" [ 0.31 0.154 0.521 0.295 0.152 0.667 0.22 0.144 0.095 0.914]\n",
" [ 0.78 0.359 0.856 0.104 0.639 0.791 0.238 0.397 0.514 0.264]\n",
" [ 0.955 0.228 0.292 0.306 0.136 0.585 0.874 0.132 0.874 0.16 ]\n",
" [ 0.022 0.427 0.906 0.626 0.032 0.585 0.43 0.901 0.422 0.984]\n",
" [ 0.559 0.271 0.36 0.543 0.95 0.307 0.604 0.42 0.193 0.082]]\n",
"max |A-Q*R| = 2.22044604925e-16\n",
"max |I-Q*Q^T| 3.64604248608e-15\n",
"---- Using SciPy QR\n",
"Q = \n",
"[[-0.187 -0.2 -0.039 0.134 -0.442 -0.424 0.372 -0.117 0.573 -0.223]\n",
" [-0.512 -0.18 0.419 0.324 -0.196 0.113 0.11 0.308 -0.451 -0.251]\n",
" [-0.158 -0.777 0.205 -0.189 0.23 0.182 -0.245 -0.343 0.161 0.064]\n",
" [-0.326 0.219 0.151 -0.474 0.231 -0.353 -0.41 0.289 0.192 -0.363]\n",
" [-0.123 0.013 -0.356 -0.449 -0.427 0.608 0.037 0.022 0.022 -0.319]\n",
" [-0.166 0.012 -0.393 0.025 0.18 -0.293 0.053 -0.566 -0.466 -0.398]\n",
" [-0.419 0.059 -0.472 0.5 -0.003 0.118 -0.501 0.057 0.224 0.175]\n",
" [-0.513 0.292 0.05 -0.104 0.437 0.24 0.521 -0.15 0.2 0.236]\n",
" [-0.012 -0.432 -0.502 -0.137 0.249 -0.188 0.295 0.573 -0.134 0.112]\n",
" [-0.3 0.028 -0.008 -0.367 -0.438 -0.295 -0.092 -0.126 -0.281 0.626]]\n",
"R = \n",
"[[-1.862 -0.993 -1.284 -0.81 -1.352 -1.127 -1.562 -0.895 -1.376 -0.938]\n",
" [ 0. -0.962 -1.115 -0.596 -0.232 -0.044 -0.559 -0.57 -0.372 -0.563]\n",
" [ 0. 0. -0.815 -0.482 -0.344 -0.787 -0.121 -0.525 -0.102 -0.981]\n",
" [ 0. 0. 0. -0.854 -0.158 -0.008 -0.212 0.093 0.141 -0.344]\n",
" [ 0. 0. 0. 0. -1.119 0.341 -0.406 -0.401 0.077 -0.074]\n",
" [ 0. 0. 0. 0. 0. -0.375 0.136 -0.424 -0.14 -0.453]\n",
" [ 0. 0. 0. 0. 0. 0. 0.702 0.435 0.467 0.393]\n",
" [ 0. 0. 0. 0. 0. 0. 0. 0.498 0.139 0.137]\n",
" [ 0. 0. 0. 0. 0. 0. 0. 0. 0.373 -0.079]\n",
" [ 0. 0. 0. 0. 0. 0. 0. 0. 0. -0.723]]\n",
"Q*R = \n",
"[[ 0.349 0.378 0.495 0.175 0.786 0.258 0.763 0.775 0.752 0.751]\n",
" [ 0.954 0.682 0.518 0.044 0.759 0.143 0.954 0.603 0.669 0.325]\n",
" [ 0.295 0.905 0.902 0.654 0.096 0.063 0.457 0.013 0.35 0.148]\n",
" [ 0.606 0.112 0.05 0.465 0.154 0.453 0.039 0.065 0.272 0.464]\n",
" [ 0.228 0.109 0.433 0.646 0.834 0.048 0.604 0.188 0.047 0.614]\n",
" [ 0.31 0.154 0.521 0.295 0.152 0.667 0.22 0.144 0.095 0.914]\n",
" [ 0.78 0.359 0.856 0.104 0.639 0.791 0.238 0.397 0.514 0.264]\n",
" [ 0.955 0.228 0.292 0.306 0.136 0.585 0.874 0.132 0.874 0.16 ]\n",
" [ 0.022 0.427 0.906 0.626 0.032 0.585 0.43 0.901 0.422 0.984]\n",
" [ 0.559 0.271 0.36 0.543 0.95 0.307 0.604 0.42 0.193 0.082]]\n",
"max |A-Q*R| = 5.55111512313e-16\n",
"max |I-Q*Q^T| 7.77156117238e-16\n"
]
}
],
"source": [
"import numpy as np\n",
"from scipy.linalg import qr as scipy_qr\n",
"from luqr import axels_qr as my_qr\n",
"\n",
"np.set_printoptions(precision=3) # only show 3 decimal places in the matrices\n",
"\n",
"# start from a random square matrix\n",
"A = np.random.rand(10, 10)\n",
"print(\"A = \", A, sep=\"\\n\")\n",
"I = np.identity(A.shape[0])\n",
"\n",
"print(\"---- Using my own QR algorithm:\")\n",
"\n",
"Q, R = my_qr(A)\n",
"print(\"Q = \", Q, \"R = \", R, \"Q*R = \", Q.dot(R), sep=\"\\n\") \n",
"print(\"max |A-Q*R| =\", np.max(np.abs(A-Q.dot(R))))\n",
"print(\"max |I-Q*Q^T|\", np.max(np.abs(I-Q.dot(Q.transpose()))))\n",
"\n",
"print(\"---- Using SciPy QR\")\n",
"\n",
"Q, R = scipy_qr(A)\n",
"print(\"Q = \", Q, \"R = \", R, \"Q*R = \", Q.dot(R), sep=\"\\n\") \n",
"print(\"max |A-Q*R| =\", np.max(np.abs(A-Q.dot(R))))\n",
"print(\"max |I-Q*Q^T|\", np.max(np.abs(I-Q.dot(Q.transpose()))))"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"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.6.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment