{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Factorización de matrices\n",
"#### https://meet.noysi.com/metodosnumericos1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Veamos algunas funciones de Sage para resolver sistemas lineales."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/latex": [
"$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrr}\n",
"1.0 & 2.0 & 1.0 \\\\\n",
"3.0 & 2.0 & 1.0 \\\\\n",
"2.0 & 4.0 & 4.0\n",
"\\end{array}\\right)$$"
],
"text/plain": [
"[1.0 2.0 1.0]\n",
"[3.0 2.0 1.0]\n",
"[2.0 4.0 4.0]"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Definimos una matrix de números en doble precisión\n",
"A = matrix(RDF,[[1,2,1],[3,2,1],[2,4,4]])\n",
"show(A)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Podemos calcular una factorización de Schur"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/latex": [
"$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrr}\n",
"0.30080722495605444 & 0.5333732643197461 & -0.790587107359547 \\\\\n",
"0.39331621941084755 & -0.8245806870099741 & -0.40665592600934425 \\\\\n",
"0.8688022588383368 & 0.18862569156680933 & 0.45782418406889225\n",
"\\end{array}\\right) \\left(\\begin{array}{rrr}\n",
"6.5033076346549805 & -1.507281212000403 & -2.6133261715360243 \\\\\n",
"0.0 & -0.8882360088434105 & 1.377907714295652 \\\\\n",
"0.0 & 0.0 & 1.384928374188435\n",
"\\end{array}\\right)$$"
],
"text/plain": [
"[ 0.30080722495605444 0.5333732643197461 -0.790587107359547]\n",
"[ 0.39331621941084755 -0.8245806870099741 -0.40665592600934425]\n",
"[ 0.8688022588383368 0.18862569156680933 0.45782418406889225] [ 6.5033076346549805 -1.507281212000403 -2.6133261715360243]\n",
"[ 0.0 -0.8882360088434105 1.377907714295652]\n",
"[ 0.0 0.0 1.384928374188435]"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"Q,T = A.schur()\n",
"show(Q,T)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Comprobamos que es la factorización buscada"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/latex": [
"$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrr}\n",
"1.0000000000000024 & 2.0000000000000013 & 1.0000000000000013 \\\\\n",
"3.0000000000000018 & 2.000000000000004 & 1.0000000000000033 \\\\\n",
"2.0000000000000044 & 4.000000000000004 & 4.0000000000000036\n",
"\\end{array}\\right) \\left(\\begin{array}{rrr}\n",
"0.9999999999999999 & 3.695068889668288 \\times 10^{-16} & 2.0753696586837751 \\times 10^{-16} \\\\\n",
"3.695068889668288 \\times 10^{-16} & 1.0000000000000002 & 2.8058042907128213 \\times 10^{-16} \\\\\n",
"2.0753696586837751 \\times 10^{-16} & 2.8058042907128213 \\times 10^{-16} & 1.0000000000000004\n",
"\\end{array}\\right)$$"
],
"text/plain": [
"[1.0000000000000024 2.0000000000000013 1.0000000000000013]\n",
"[3.0000000000000018 2.000000000000004 1.0000000000000033]\n",
"[2.0000000000000044 4.000000000000004 4.0000000000000036] [ 0.9999999999999999 3.695068889668288e-16 2.0753696586837751e-16]\n",
"[ 3.695068889668288e-16 1.0000000000000002 2.8058042907128213e-16]\n",
"[2.0753696586837751e-16 2.8058042907128213e-16 1.0000000000000004]"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"show(Q*T*Q.conjugate_transpose(),Q*Q.conjugate_transpose())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"También una factorización LU (de Doolittle). \n",
"**¡Ojo!** Sage no devuelve en la forma $PA=LU$, sino $A=PLU$."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/latex": [
"$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrr}\n",
"0.0 & 0.0 & 1.0 \\\\\n",
"1.0 & 0.0 & 0.0 \\\\\n",
"0.0 & 1.0 & 0.0\n",
"\\end{array}\\right) \\left(\\begin{array}{rrr}\n",
"1.0 & 0.0 & 0.0 \\\\\n",
"0.6666666666666666 & 1.0 & 0.0 \\\\\n",
"0.3333333333333333 & 0.5 & 1.0\n",
"\\end{array}\\right) \\left(\\begin{array}{rrr}\n",
"3.0 & 2.0 & 1.0 \\\\\n",
"0.0 & 2.666666666666667 & 3.3333333333333335 \\\\\n",
"0.0 & 0.0 & -1.0\n",
"\\end{array}\\right) \\left(\\begin{array}{rrr}\n",
"1.0 & 2.0 & 1.0 \\\\\n",
"3.0 & 2.0 & 1.0 \\\\\n",
"2.0 & 4.0 & 4.0\n",
"\\end{array}\\right)$$"
],
"text/plain": [
"[0.0 0.0 1.0]\n",
"[1.0 0.0 0.0]\n",
"[0.0 1.0 0.0] [ 1.0 0.0 0.0]\n",
"[0.6666666666666666 1.0 0.0]\n",
"[0.3333333333333333 0.5 1.0] [ 3.0 2.0 1.0]\n",
"[ 0.0 2.666666666666667 3.3333333333333335]\n",
"[ 0.0 0.0 -1.0] [1.0 2.0 1.0]\n",
"[3.0 2.0 1.0]\n",
"[2.0 4.0 4.0]"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"P,L,U = A.LU()\n",
"show(P,L,U,P*L*U)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(\n",
"[1 0 0] [1 0 0] [ 1 2 1]\n",
"[0 1 0] [3 1 0] [ 0 -4 -2]\n",
"[0 0 1], [2 0 1], [ 0 0 2]\n",
")"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Definimos una matrix de números enteros\n",
"B = matrix([[1,2,1],[3,2,1],[2,4,4]])\n",
"B.LU(pivot='nonzero')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"También podemos resolver directamente el sistema $Ax = b$ usando la función *solve_right* (usa eliminación gaussiana)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.5, 0.0, 0.5)"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"b = vector(RDF,[1,2,3])\n",
"A.solve_right(b)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Otra manera de denotarlo es con el operador \"\\\\\"."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.5, 0.0, 0.5)"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A\\b"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(1.0, 2.0, 3.0)"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Comprobamos\n",
"A*(A\\b)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Para medir el tiempo que tarda en realizar una operación, usaremos la función *time* de la librería *time*. Para ello, primero cargamos la librería, después guardamos el tiempo inicial, ponemos nuestro código y guardamos el tiempo final. "
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.0013327598571777344"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import time\n",
"st = time.time() # Tiempo inicial \n",
"\n",
"# Código\n",
"factorial(10000)\n",
"\n",
"et = time.time() # Tiempo final\n",
"et - st # Diferencia"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Vamos a calcular cuánto tiempo tarda en resolver un sistema de 100x100"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"# Generamos una matriz aleatoria de 100x100\n",
"A = matrix.random(RDF,100,100)\n",
"# Generamos un vector aleatorio de 100 elementos\n",
"b = vector([random() for k in range(100)])"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.6425516605377197"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import time\n",
"st = time.time() # Tiempo inicial \n",
"\n",
"# Resolvemos el sistema, pero no lo mostramos (para no calcular ese tiempo)\n",
"for _ in range(1000):\n",
" xs = A\\b\n",
"\n",
"et = time.time() # Tiempo final\n",
"et - st # Tiempo que ha tardado"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
" Ejercicio 1. \n",
"Consideremos la matriz:\n",
"$$\n",
"A=\\left(\\begin{array}{rrr}\n",
"1 & 2 & 1 \\\\\n",
"3 & 2 & 1 \\\\\n",
"-2 & -4 & 4\n",
"\\end{array}\\right)$$\n",
" \n",
"**a)** Obtener la factorización de Schur. ¿Qué ocurre? Calcular los autovalores de la matriz $A$ y discutir porqué no factoriza en los reales como una matriz unitaria por una matriz triangular. Obtener la factorización en los números complejos en doble precisión (CDF)\n",
"\n",
"**b)** Resolver el sistem $Ax=b$ con $b=(1,2,3)$ usando la factorización anterior.\n",
"\n",
"**c)** Consideremos ahora la aplicación lineal dada por\n",
"$$A=\\left(\\begin{array}{rrr}\n",
"4 & 2 & -2 \\\\\n",
"2 & 2 & -4 \\\\\n",
"-2 & -4 & 11 \\\\\n",
"\\end{array}\\right)$$\n",
"Calcular su factorización de Schur. Deducir que la matriz es definida positiva. Representar las columnas de la matriz $U$ y también las de $A U$. Comprobar que en ambos casos son vectores ortogonales.\n",
"\n",
"d) Obtener la factorización de Schur de cualquiera de las matrices anteriores \"paso a paso\", siguiendo la demostración vista en clase. Para calcular un autovector se puede usar la función eigenvectors_right y para extender a una base ortonormal, se puede utilizar la función gram_schmidt de Sage. \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/latex": [
"$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrr}\n",
"1.0 & 2.0 & 1.0 \\\\\n",
"3.0 & 2.0 & 1.0 \\\\\n",
"-2.0 & -4.0 & 4.0\n",
"\\end{array}\\right)$$"
],
"text/plain": [
"[ 1.0 2.0 1.0]\n",
"[ 3.0 2.0 1.0]\n",
"[-2.0 -4.0 4.0]"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"A = matrix(RDF, [ [ 1, 2, 1],\n",
" [ 3, 2, 1],\n",
" [-2,-4, 4]])\n",
"show(A)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(\n",
"[ -0.7318717503574348 -0.16272623270313655 0.6617279759984384]\n",
"[ 0.6448378884837489 0.14860231667063856 0.7497342522893958]\n",
"[ 0.22033574064269745 -0.9754165903194197 0.0038257973117043517],\n",
"\n",
"[ -1.0632187495592365 -1.0556113258123003 -0.28402152320817975]\n",
"[ 0.0 4.031609374779613 4.3672968471306906]\n",
"[ 0.0 -1.446911310750419 4.031609374779613]\n",
")"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A.schur()"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[-1.0632187495592365,\n",
" 4.031609374779613 + 2.513780261979563*I,\n",
" 4.031609374779613 - 2.513780261979563*I]"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A.eigenvalues()"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"A = matrix(CDF, [ [ 1, 2, 1],\n",
" [ 3, 2, 1],\n",
" [-2,-4, 4]])"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"Q,T = A.schur()"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5.528169885528619e-16"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(Q*Q.conjugate_transpose() - matrix.identity(3)).norm(2)"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [],
"source": [
"b = vector([1,2,3])\n",
"x0 = Q.conjugate_transpose()*b\n",
"x1 = T\\x0\n",
"x2 = Q*x1"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.0785460984167022e-15"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"( A*x2 - b ).norm(2)"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/latex": [
"$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrr}\n",
"4.0 & 2.0 & -2.0 \\\\\n",
"2.0 & 2.0 & -4.0 \\\\\n",
"-2.0 & -4.0 & 11.0\n",
"\\end{array}\\right)$$"
],
"text/plain": [
"[ 4.0 2.0 -2.0]\n",
"[ 2.0 2.0 -4.0]\n",
"[-2.0 -4.0 11.0]"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"A = matrix(RDF, [ [ 4, 2, -2],\n",
" [ 2, 2, -4],\n",
" [-2,-4, 11]])\n",
"show(A)"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/latex": [
"$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrr}\n",
"-0.27149455714159304 & -0.9065005265405895 & -0.32333805965911955 \\\\\n",
"-0.3649468832088046 & -0.21390284460681533 & 0.906123250725469 \\\\\n",
"0.8905641346240427 & -0.36400874776679004 & 0.2727510837202741\n",
"\\end{array}\\right) \\left(\\begin{array}{rrr}\n",
"13.24888536294344 & 1.4918621893400541 \\times 10^{-15} & -9.114968950728652 \\times 10^{-16} \\\\\n",
"0.0 & 3.668823351415331 & 1.1650132811216793 \\times 10^{-15} \\\\\n",
"0.0 & 0.0 & 0.08229128564122203\n",
"\\end{array}\\right)$$"
],
"text/plain": [
"[-0.27149455714159304 -0.9065005265405895 -0.32333805965911955]\n",
"[ -0.3649468832088046 -0.21390284460681533 0.906123250725469]\n",
"[ 0.8905641346240427 -0.36400874776679004 0.2727510837202741] [ 13.24888536294344 1.4918621893400541e-15 -9.114968950728652e-16]\n",
"[ 0.0 3.668823351415331 1.1650132811216793e-15]\n",
"[ 0.0 0.0 0.08229128564122203]"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"U,D = A.schur()\n",
"show(U,D)"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.0000000000000007"
]
},
"execution_count": 44,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"U.column(0)*U.column(0)"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4.4640274541081684e-14"
]
},
"execution_count": 48,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(A*U).column(0)*(A*U).column(1)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
" Ejercicio 2. \n",
" \n",
"**a)** Obtener la factorización LU con pivote de la matriz \n",
"$$\n",
"A=\\left(\\begin{array}{rrr}\n",
"1 & 2 & 1 \\\\\n",
"3 & 2 & 1 \\\\\n",
"3 & -4 & 4\n",
"\\end{array}\\right)$$\n",
"\n",
"**b)** Utilizar dicha factorización para obtener la inversa de $A$.\n",
"\n",
"c) Obtener la factorización LU con pivote de $A^tA$. A partir de dicha factorización, obtener la factorización de la forma $L D L^t$ y a partir de ella la factorización de Choleski. \n",
"\n",
"d) Representar en una gráfica el tiempo de ejecución de la resolución de los sistemas lineales generados aleatoriamente para matrices de tamaño hasta 1000x1000. Debería ser aproximadamente la gráfica de una función cúbica. \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [],
"source": [
"A = matrix([[1,2,1],[3,2,1],[3,-4,4]])"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(\n",
"[0 0 1] [ 1 0 0] [ 3 2 1]\n",
"[1 0 0] [ 1 1 0] [ 0 -6 3]\n",
"[0 1 0], [ 1/3 -2/9 1], [ 0 0 4/3]\n",
")"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P,L,U = A.LU()\n",
"P,L,U"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Resolver $A X = I$, $P L U X = I$.\n",
"\n",
"Podemos hacer\n",
"\n",
"$X = (x1,x2,x3)$ y resolver $P L U x1 = (1,0,0)$, ..."
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"((1, 0, 0), (0, 1, 0), (0, 0, 1))"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"e1,e2,e3 = matrix.identity(3)\n",
"e1,e2,e3"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(-1/2, 3/8, 3/4)"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"z1 = P.transpose()*e1 # Transponiendo\n",
"y1 = L\\z1 # Sustitución progresiva\n",
"x1 = U\\y1 # Sustitución regresiva\n",
"x1"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(1/2, -1/24, -5/12)"
]
},
"execution_count": 42,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"z2 = P.transpose()*e2 # Transponiendo\n",
"y2 = L\\z2 # Sustitución progresiva\n",
"x2 = U\\y2 # Sustitución regresiva\n",
"x2"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0, -1/12, 1/6)"
]
},
"execution_count": 43,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"z3 = P.transpose()*e3 # Transponiendo\n",
"y3 = L\\z3 # Sustitución progresiva\n",
"x3 = U\\y3 # Sustitución regresiva\n",
"x3"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[ -1/2 1/2 0]\n",
"[ 3/8 -1/24 -1/12]\n",
"[ 3/4 -5/12 1/6]"
]
},
"execution_count": 44,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Ai = matrix([x1,x2,x3]).transpose()\n",
"Ai"
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/latex": [
"$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrr}\n",
"1 & 0 & 0 \\\\\n",
"0 & 1 & 0 \\\\\n",
"0 & 0 & 1\n",
"\\end{array}\\right)$$"
],
"text/plain": [
"[1 0 0]\n",
"[0 1 0]\n",
"[0 0 1]"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"show(Ai*A)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "SageMath 9.5",
"language": "sage",
"name": "sagemath"
},
"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.2"
}
},
"nbformat": 4,
"nbformat_minor": 4
}