Sean: $$A = \left( \begin{array}{ccc} 1 & 0 & -1 \\ 2 & 1 & 3 \\ 0 & 1 & -1 \end{array} \right) \qquad B = \left( \begin{array}{ccc} 0 & -2 & 1 \\ 3 & -1 & 0 \\ 2 & 0 & -1 \end{array} \right) \qquad v = (-1, 4, -2)$$ Computar en Python la siguiente operación: $A^tA^2v^t+ BAv^t + <v, A_{1 \cdot}>v^t$ siendo $A_{1 \cdot}$ la primera fila de $A$
Respuesta: $(40, 17, 57)$
Posible implementación: A.T @ (A @ A) @ v + B @ A @ v + (v @ A[0,:])*v
Lema: sea $A\in \mathbb{R}^{n\times n}$ simétrica definida positiva, existe una única $L\in \mathbb{R}^{n\times n}$ triangular inferior tal que $A = LL^t$ y $l_{ii}>0 \; \forall 1\leq i \leq n$
Veamos qué ocurre en el caso de una matriz $A\in \mathbb{R}^{3\times 3}$ simétrica definida positiva: $$ A = \left( \begin{array}{ccc} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{array} \right) $$ Queremos escribirla como el producto $LL^t$ donde: $$ L = \left( \begin{array}{ccc} l_{11} & 0 & 0 \\ l_{21} & l_{22} & 0 \\ l_{31} & l_{32} & l_{33} \end{array} \right) $$
Luego: $$ \left( \begin{array}{ccc} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{array} \right) = \left( \begin{array}{ccc} l_{11} & 0 & 0 \\ l_{21} & l_{22} & 0 \\ l_{31} & l_{32} & l_{33} \end{array} \right) \left( \begin{array}{ccc} l_{11} & l_{21} & l_{31} \\ 0 & l_{22} & l_{32} \\ 0 & 0 & l_{33} \end{array} \right) = \left( \begin{array}{ccc} l_{11}^2 & l_{21}l_{11} & l_{31}l_{11} \\ l_{21}l_{11} & l_{21}^2+l_{22}^2 & l_{31}l_{21}+l_{32}l_{22} \\ l_{31}l_{11} & l_{31}l_{21}+l_{32}l_{22} & l_{31}^2 + l_{32}^2 + l_{33}^2 \end{array} \right) $$ Entonces: $$L = \left( \begin{array}{ccc} \sqrt{a_{11}} & 0 & 0 \\ \dfrac{a_{21}}{l_{11}} & \sqrt{a_{22}-l_{21}^2} & 0 \\ \dfrac{a_{31}}{l_{11}} & \dfrac{a_{32}-l_{31}l_{21}}{l_{22}} & \sqrt{a_{33} - l_{31}^2-l_{32}^2} \end{array} \right) $$
Si tenemos un sistema de ecuaciones lineales representado por $Ax=b$ donde $A \in\mathbb{R}^{n\times n}$ es triangular inferior y $a_{ii}\neq0$, entonces el sistema ya está triangulado y los valores de $x$ son fácilmente calculables.
Si tenemos un sistema de ecuaciones lineales representado por $Ax=b$ donde $A \in\mathbb{R}^{n\times n}$ es triangular superior y $a_{ii}\neq0$, entonces el sistema ya está triangulado y los valores de $x$ son fácilmente calculables.
Se puede iterar de manera decreciente utilizando:
for i in reversed(range(n))
Supongamos que tenemos el sistema de ecuaciones lineales $Ax=b$ donde $A$ es simétrica definida positiva. Por descomposición de Cholesky, existe $L$ triangular inferior tal que $A = LL^t$ y la diagonal de $L$ no tiene elementos nulos. Entonces: $$Ax=b \Leftrightarrow LL^tx=b$$ Sea $y=L^tx$, se tiene que, para resolver el problema original, basta resolver: $$\begin{cases} Ly = b \\ L^tx = y \end{cases} $$
Ejercicio: implementar una función que resuelva sistemas de ecuaciones lineales con matriz asociada simétrica definida positiva
Lema: sea $A\in\mathbb{R}^{n\times n}$, luego se puede escribir como el producto de una matriz ortonormal $Q$ (es decir, $Q$ tal que $Q^t = Q^{-1}$) y una matriz triangular superior $R$. Si $A$ es inversible, la descomposición es única y la diagonal de $R$ es postiva.
Sea $A=[a_1|a_2|\dots|a_n]$ cuyas columnas son linealmente independientes, les aplicaremos el método de Gram-Schmidt. Sea la proyección $$proy_ua = \dfrac{<u,a>}{<u,u>}u$$ Hallamos la base ortonormal $\{w_1,\dots,w_n\}$ de la siguiente manera:
$\begin{array}{cclrcl} u_1 &=& a_1 & w_1 & = &\dfrac{u_1}{||u_1||}\\ u_2 &=& a_2 - proy_{u_1}a_2 & w_2 & = &\dfrac{u_2}{||u_2||}\\ u_3 &=& a_3 - proy_{u_1}a_3 - proy_{u_2}a_3 & w_3 & = &\dfrac{u_3}{||u_3||}\\ \vdots & & \\ u_k & = & a_k - \displaystyle \sum_{j=1}^{k-1}proy_{u_j}a_k & w_k & = &\dfrac{u_k}{||u_k||} \end{array}$
Se tiene entonces que: $$Q = [w_1|w_2|\dots|w_n] \qquad \displaystyle R={\begin{pmatrix}\langle w_{1},{a} _{1}\rangle &\langle w_{1},{a} _{2}\rangle &\langle w _{1},{a} _{3}\rangle &\ldots \\0&\langle w_{2},{a} _{2}\rangle &\langle w_{2},{a} _{3}\rangle &\ldots \\0&0&\langle w_{3},{a} _{3}\rangle &\ldots \\\vdots &\vdots &\vdots &\ddots \end{pmatrix}}$$
Sea $A\in\mathbb{R}^{n\times n}$, entonces podemos aplicar el siguiente procedimiento para hallar, de existir, la solución de $Ax=b$
$$Ax=b \Leftrightarrow QRx=b \Leftrightarrow Rx=Q^tb$$Sea $\{(x_1,y_1),\dots,(x_n,y_n)\}$ un conjunto de datos en el plano, el objetivo es hallar $P\in\mathcal{P}_k$ que mejor aproxime a los datos en el sentido de minimizar la siguiente expresión: $$\displaystyle\sum_{i=1}^n (y_i-P(x_i))^2$$
Consideremos $$A = \left( \begin{array}{ccccc} 1 & x_1 & x_1^2 & \dots & x_1^k \\ 1 & x_2 & x_2^2 & \dots & x_2^k \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & x_n^2 & \dots & x_n^k \end{array} \right) $$
Lo ideal sería encontrar $\alpha=(\alpha_1,\dots,\alpha_{k+1})^t$ tal que suceda lo siguiente: $$\left( \begin{array}{ccccc} 1 & x_1 & x_1^2 & \dots & x_1^k \\ 1 & x_2 & x_2^2 & \dots & x_2^k \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & x_n^2 & \dots & x_n^k \end{array} \right) \left( \begin{array}{c} \alpha_0 \\ \alpha_2 \\ \vdots \\ \alpha_{k} \end{array} \right) = \left( \begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_n \end{array} \right) $$ Pero eso rara vez sucede.
$\{1,x,\dots,x^k\}$ es una base de $\mathcal{P}_k$. Nos gustaría encontrar $P=\displaystyle\sum_{j=0}^{k}\alpha_j x^j$ que minimice
$$E(\alpha) = \displaystyle\sum_{i=1}^n \left( y_i-\displaystyle\sum_{j=0}^{k}\alpha_j x_i^j \right)^2$$Entonces buscamos $\alpha$ tal que: $$\dfrac{\partial E}{\partial \alpha_m} = \displaystyle\sum_{i=1}^n 2\left( y_i- \displaystyle\sum_{j=1}^k \alpha_j x_i^j \right)\cdot x_i^m = 0 \qquad \forall 0\leq m \leq k$$ $$ \Leftrightarrow \displaystyle\sum_{i=1}^n\displaystyle\sum_{j=1}^k \alpha_j x_i^jx_i^m = \displaystyle\sum_{i=1}^n y_ix_i^m \qquad \forall 0\leq m \leq k $$ Esto es justamente pedir que $\alpha$ sea la solución de $A^tA\alpha = A^ty$ con $y=(y_1,\dots,y_n)^t$
Notar que si $A\in\mathbb{R}^{n\times m}$, entonces $A^tA\in\mathbb{R}^{m\times m}$ y es definida positiva. Por lo visto anteriormente, dados los datos $\{(x_1,y_1),\dots,(x_n,y_n)\}$, para hallar los coeficientes del polinomio de grado $k$ que mejor los aproxima en sentido de cuadrados mínimos basta resolver:
$$A^tA\alpha = A^ty$$Siendo:
$$A = \left( \begin{array}{ccccc} 1 & x_1 & x_1^2 & \dots & x_1^k \\ 1 & x_2 & x_2^2 & \dots & x_2^k \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & x_n^2 & \dots & x_n^k \end{array} \right) \qquad y = \left( \begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_n \end{array} \right) \qquad \alpha = \left( \begin{array}{c} \alpha_0 \\ \alpha_2 \\ \vdots \\ \alpha_{k} \end{array} \right) $$Implementar dos funciones: ambas deben tomar como input el grado del polinomio y los datos y devolver una lista con los coeficientes del polinomio que mejor los aproxima. Una de las funciones debe resolver el sistema usando descomposición de Cholesky y la otra aplicando descomposición QR.
Vamos a comparar gráficamente el resultado de ambos métodos. Para eso, utilizaremos el archivo datos.txt.
Importamos numpy y matplotlib.pyplot. También importamos la función polyval que permite evaluar polinomios a partir de sus coeficientes de la siguiente manera: polyval(x,c) donde x puede ser un número o un array y c es la lista de coeficientes del polinomio ordenados de la siguiente manera: c[i] es el coeficiente de $x^i$.
import numpy as np
from numpy.polynomial.polynomial import polyval
import matplotlib.pyplot as plt
datos_X = np.loadtxt('datos.txt', usecols=0) # Como los datos de X están en la primera columna, usamos usecols=0
datos_Y = np.loadtxt('datos.txt', usecols=1) # Como los datos de X están en la segunda columna, usamos usecols=1
Ahora podemos aplicar las funciones implementadas para hallar los coeficientes del los polinomios. La idea es graficar los datos (plt.plot(X,Y)) y los dos polinomios obtenidos, todo en el mismo gráfico.
Necesitamos la librería image de matplotlib.
import matplotlib.image as mimg
import matplotlib.pyplot as plt
archivo = 'casita.png'
A = mimg.imread(archivo) # Lee el archivo de la imagen y lo convierte en matriz
plt.imshow(A, cmap='gray') # Estos dos comandos permiten mostrar la imagen en pantalla
plt.show()
mimg.imsave('casita_2.png', A, cmap='gray') # Guarda la imagen representada por A en un archivo