A = [
1 2
1 3
]
M = [
4 -1 1
-1 4.25 2.75
1 2.75 3.5
];
[P, D] = eig(M)
Pinv = inv(P);
P*D*Pinv
No todas las matrices son diagonalizables¶
Algunas matrices, como $$ A = \left(\begin{array}{cc} 0 & 1\\ 0 & 0 \end{array}\right) $$ no son diagonalizables. En el caso de la matriz $A$, el polinomio característico es $p_A(\lambda)=\lambda^2$, y el único autovalor es $\lambda=0$. Si $A$ fuera diagonalizable: $$ A = P\cdot D \cdot P^{-1}, $$ entonces tendríamos que $p_A(\lambda)=p_D(\lambda)$, porque: $$ p_A(\lambda) = \det(A - \lambda Id) = \det(P^{-1})\cdot \det(A - \lambda Id) \cdot \det(P)= \det(P^{-1}\cdot (A - \lambda Id)\cdot P) = \det(P^{-1}\cdot A\cdot P - \lambda P^{-1}\cdot Id\cdot P) = p_D(\lambda) $$ Por tanto, $D$ sería una matriz diagonal con polinomio característico $\lambda^2$, y tendría que ser la matriz cero, pero entonces $A$ también lo sería y no lo es...
Observa el resultado de buscar los autovalores y autovectores de A: ¿detectas algún problema?
A = [0,1;0,0]
[P, D] = eig(A)
Sin embargo, si añadimos un ruido a la matriz $A$, se vuelve diagonalizable (aunque los autovalores pueden ser un par de números complejos conjugados)
% Un numero real aleatorio en el intervalo [0,1]
rand
% Con esta perturbacion aleatoria,
% siempre obtenemos autovalores reales y distintos
% Uno de ellos es siempre cero
epsilon = 1e-4
A = [epsilon*rand,1; 0,0]
[P, D] = eig(A)
% Con esta otra perturbacion
% siempre obtenemos dos autovalores reales y distintos
% ¿por que? ¿cual es el polinomio caracteristico?
epsilon = 1e-4;
A = [0,1; epsilon*rand,0]
[P, D] = eig(A)
% Un numero aleatorio extraído de una Normal(0,1)
randn
% Con esta otra perturbacion aleatoria,
% a veces obtenemos dos autovalores reales y distintos
% a veces obtenemos dos autovalores complejos imaginarios puros conjugados
% ¿por que?
% Repite varias veces: el resultado es aleatorio
epsilon = 1e-4;
A = [0,1; epsilon*randn,0]
[P, D] = eig(A)
Ejercicio¶
- Explica qué autovalores podemos obtener al introducir un ruido aleatorio en la matriz, de las tres formas descritas:
- $\epsilon$ es un número aleatorio en el intervalo (0,1) $$ A = \left(\begin{array}{cc} \epsilon & 1\\ 0 & 0 \end{array}\right) $$
- $\epsilon$ es un número aleatorio en el intervalo (0,1) $$ A = \left(\begin{array}{cc} 0 & 1\\ \epsilon & 0 \end{array}\right) $$
- $\epsilon$ es un número aleatorio extraído de una Normal(0,1) (puede ser cualquier número real, positivo o negativo) $$ A = \left(\begin{array}{cc} 0 & 1\\ \epsilon & 0 \end{array}\right) $$
Probamos ahora una perturbación general, números positivos o negativos en cualquiera de las cuatro entradas.
%Un array 2x2 de numeros reales en el intervalo [0,1]
randn(2,2)
Ahora los autovalores no son necesariamente reales o imaginarios puros, pero cuando son números complejos, son complejos conjugados: ¿por qué?
% Repite varias veces: el resultado es aleatorio
epsilon = 1e-4;
ruido = epsilon*randn(2,2);
A = [0,1; 0,0] + ruido
[P, D] = eig(A)
Casi todas las matrices son diagonalizables¶
Tras añadir un ruido aleatorio a cualquier matriz $n\times n$ (en todas sus entradas), el polinomio característico tiene $n$ raíces distintas con probabilidad $1$.
Es importante tener ésto en cuenta, porque los errores de redondeo son frecuentes...
Factorización LU y PLU¶
Hay otras muchas formas de diagonalizar matrices. Nos interesa ahora la factorización LU, porque ayuda a resolver sistemas de ecuaciones lineales de forma más eficiente.
Recordatorio del método de Gauss¶
Si no recuerdas el método de Gauss, puedes usar esta página web interactiva que ejecuta el método de Gauss paso a paso, o leer el capítulo 6.1 del libro de Burden y Faires, o tus apuntes de Álgebra Lineal...
Recordamos que el método de Gauss para resolver sistemas de ecuaciones consiste en aplicar de cierta forma sistemática las siguientes operaciones:
- Intercambiar dos filas de una matriz
- Multiplicar una fila por un número
- Sumar a una fila un múltiplo de otra
Para una matriz $n\times n$ general, puede ser necesario hacer la última operación n-1 veces para hacer ceros debajo del primer pivote, después $n-2$ veces para hacer ceros debajo del segundo pivote, etc, lo que arroja un total de $$ (n-1) + (n-2) + \dots + 1 = n(n-1)/2=\frac{n^2}{2} + \text{términos de menor grado}. $$ aplicaciones de la operación "suma a una fila un múltiplo de otra". Sumar a una fila un múltiplo de otra requiere n multiplicaciones y n-1 sumas, aunque si la fila comienza con ceros, nos ahorramos unas cuantas operaciones...
En total, si la matriz no es "especial", hacen falta $2n^3/3$ operaciones para resolver un sistema de ecuaciones con este método, como puedes consultar en la bibliografía.
Sin embargo, si el sistema es triangular inferior, son necesarias muchas menos operaciones.
La solución del sistema: $$\begin{matrix}\ell _{1,1}x_{1}&&&&&&&=&b_{1}\\\ell _{2,1}x_{1}&+&\ell _{2,2}x_{2}&&&&&=&b_{2}\\\vdots &&\vdots &&\ddots &&&&\vdots \\\ell _{m,1}x_{1}&+&\ell _{m,2}x_{2}&+&\dotsb &+&\ell _{m,m}x_{m}&=&b_{m}\\\end{matrix} $$ se puede encontrar despejando las variables sucesivas: $$ \begin{aligned}x_{1}&={\frac {b_{1}}{\ell _{1,1}}},\\x_{2}&={\frac {b_{2}-\ell _{2,1}x_{1}}{\ell _{2,2}}},\\&\ \ \vdots \\x_{m}&={\frac {b_{m}-\sum _{i=1}^{m-1}\ell _{m,i}x_{i}}{\ell _{m,m}}}.\end{aligned} $$ con un total de $n^2$ operaciones (+ términos de menor orden)
Si el sistema es triangular superior, se puede encontrar la solución de forma similar.
Factorización LU¶
Por tanto, sería interesante poder factorizar una matriz arbitraria $M$ como el producto de una matriz triangular inferior y otra triangular superior $$ M = L\cdot U. $$ Una vez hallada la factorización, podríamos resolver el sistema $M\cdot \mathbf{x} = \mathbf{b}$ de la forma siguiente:
- Resuelve primero $L\cdot \mathbf{y} = \mathbf{b}$ => el resultado es un vector $\mathbf{y}$, tras $\approx n^2$.
- Resuelve ahora $U\cdot \mathbf{x} = \mathbf{y}$, tras $\approx n^2$.
El vector $\mathbf{x}$ es solución del problema original: $$ M\cdot \mathbf{x} = L\cdot U \cdot \mathbf{x} = L\cdot \mathbf{y} = \mathbf{b}, $$ y hemos necesitado $\approx 2n^2$ operaciones en vez de $\approx 2n^3/3$, que cuando $n$ es grande es mucho mayor.
La factorización LU no siempre existe¶
Desgraciadamente, no siempre es posible encontrar una factorización LU. Por ejemplo, la matriz $$ A = \left(\begin{array}{cc} 0 & 1\\ 1 & 0 \end{array}\right) $$ no puede ser el producto de una matriz triangular inferior y otra superior. $$ L = \left(\begin{array}{cc} a & 0\\ b & c \end{array}\right),\; U = \left(\begin{array}{cc} d & e\\ 0 & f \end{array}\right) $$
$$ L\cdot U = \left(\begin{array}{cc} ad & ae\\ bd & be+cf \end{array}\right)\neq A $$dado que si $ad=0$, entonces o bien $ae$ o bien $bd$ es cero.
- Las matrices a las que se puede aplicar el método de Gauss sin permutar filas siempre admiten una factorización LU (más detalles en el libro).
Factorización PLU¶
Si añadimos un ruido a una matriz, es casi seguro que admitirá una factorización LU. Sin embargo, la factorización PLU no sólo existe siempre, sino que además es menos sensible a errores numéricos.
El método de Gauss se puede aplicar a cualquier matriz invertible, pero puede ser necesario permutar filas. Siguiendo la pista de estas operaciones, podemos encontrar una factorización un poco distinta: $$ PA = LU $$ ó $$ A = PLU $$ donde $P$ es una matriz de permutación: una matriz que en cada fila tiene todos ceros menos exactamente un uno, y en cada columna tiene todo ceros menos exactamente un uno.
Por ejemplo: las únicas matrices de permutación $2\times 2$ son la identidad y: $$ P = \left(\begin{array}{cc} 0 & 1\\ 1 & 0 \end{array}\right) $$
Las matrices de permutación $3\times 3$ son seis: $$ \left(\begin{array}{ccc} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{array}\right),\; \left(\begin{array}{ccc} 0 & 1 & 0\\ 0 & 0 & 1\\ 1 & 0 & 0 \end{array}\right),\; \left(\begin{array}{ccc} 0 & 0 & 1\\ 1 & 0 & 0\\ 0 & 1 & 0 \end{array}\right),\; $$ $$ \left(\begin{array}{ccc} 0 & 1 & 0\\ 1 & 0 & 0\\ 0 & 0 & 1 \end{array}\right),\; \left(\begin{array}{ccc} 1 & 0 & 0\\ 0 & 0 & 1\\ 0 & 1 & 0 \end{array}\right),\; \left(\begin{array}{ccc} 0 & 0 & 1\\ 0 & 1 & 0\\ 1 & 0 & 0 \end{array}\right).\; $$
Se llaman matrices de permutación porque el efecto de multiplicar una matriz $A$ por una matriz de permutación es permutar (cambiar de orden) las filas de $A$: $$ \left(\begin{array}{ccc} 0 & 1 & 0\\ 1 & 0 & 0\\ 0 & 0 & 1 \end{array}\right)\cdot \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} a_{21} & a_{22} & a_{23}\\ a_{11} & a_{12} & a_{13}\\ a_{31} & a_{32} & a_{33} \end{array}\right) $$
Resolver un sistema de ecuaciones $$ P\cdot \mathbf{x} = \mathbf{b} $$ es también muy sencillo cuando $P$ es una matriz de permutación: $\mathbf{x}$ es igual a $\mathbf{b}$, pero después de permutar sus filas. Otra forma de verlo es usar que las matrices de permutación son ortogonales: $$ P\cdot P^t = P^t\cdot P = Id, $$ luego, multiplicando la ecuación anterior por $P^t$: $$ P^t\cdot P\cdot \mathbf{x} = \mathbf{x} = P^t\mathbf{b}, $$ aunque un producto de matriz por vector requiere $\approx n^2$ operaciones, y una permutación sólo $n$.
Por tanto, resolver un sistema $P\cdot L\cdot U\cdot \mathbf{x} = \mathbf{b}$ es también sencillo:
- Multiplica primero por $P^t$ (o permuta directamente las filas de $\mathbf{b}$) para obtener $\tilde{\mathbf{b}}=P^t\cdot \mathbf{b}$.
- Resuelve después $L\cdot \mathbf{y} = \tilde{\mathbf{b}}$ => el resultado es un vector $\mathbf{y}$.
- Resuelve finalmente $U\cdot \mathbf{x} = \mathbf{y}$, tras $\approx 2n^2$.
El vector $\mathbf{x}$ es solución del problema original: $$ M\cdot \mathbf{x} = P\cdot L\cdot U \cdot \mathbf{x} = P\cdot L\cdot \mathbf{y} = P\cdot\tilde{\mathbf{b}} = \mathbf{b}, $$
Una de las dos matrices de la factorización tiene 1 en la diagonal¶
El número de grados de libertad de las matrices $n\times n$ es por supuesto $n^2$.
¿Cuántos grados de libertad tienen una matriz triangular superior y una triangular inferior? Resulta que son más de $\frac{n(n+1)}{2}$ cada uno, luego de forma conjunta tienen más de $n^2$ grados de libertad.
Sin embargo, si además de pedir a $L$ que sea triangular inferior le pedimos que tenga unos en la diagonal, el número de grados de libertad conjunto es $\frac{n(n-1)}{2} + \frac{n(n+1)}{2} = n^2$ .
¿Cuál debe tener unos en la diagonal, $L$, o $U$? Depende del manual, o de la librería de software (es importante leer la documentación).
Teorema¶
Teorema: toda matriz invertible $A$ admite al menos una factorización $$ A = PLU $$
- P es una matriz de permutación
- L es una matriz triangular con unos en la diagonal
- U es una matriz superior
Atención: toda matriz invertible $A$ admite también una factorización $PA = LU$. Es muy similar, y está relacionada con la anterior, pero no son idénticas. Es importante saber cuál de ellas calcula nuestro software.
Aplicación al polinomio interpolador¶
Buscamos el polinomio interpolador
$$ p(x) = a_0 + a_1\cdot x + \dots + a_n\cdot x^n $$que pasa por los $n+1$ puntos de interpolación $$ p(x_0)=y_0,\:p(x_1)=y_1,\dots,p(x_n)=y_n. $$ Desarrollando esta expresión, los coeficientes del polinomio interpolador son la solución del sistema de ecuaciones $$ V\cdot \mathbf{a} = \mathbf{y}, $$ donde $$ {\displaystyle V={\begin{bmatrix}1&x _{0}&x _{0}^{2}&\dots &x _{0}^{n}\\1&x _{1}&x _{1}^{2}&\dots &x _{1}^{n}\\1&x _{2}&x _{2}^{2}&\dots &x _{2}^{n}\\\vdots &\vdots &\vdots &\ddots &\vdots \\1&x _{n}&x _{n}^{2}&\dots &x _{n}^{n}\end{bmatrix}}.} $$ $$ {\displaystyle \mathbf{a}={\begin{bmatrix}a _{0}\\a_1\\\vdots\\a_n\end{bmatrix}}; \mathbf{y}={\begin{bmatrix}y_{0}\\y_1\\\vdots\\y_n\end{bmatrix}},} $$
Si los puntos $x_0,\dots,x_n$ los conocemos de antemano, pero los puntos $y_0,\dots,y_n$ pueden depender de la aplicación, nos interesa resolver muchos sistemas del tipo $$ V\cdot \mathbf{a} = \mathbf{y} $$ para muchos vectores $\mathbf{y}$ distintos. Lo podemos hacer aprovechando una factorización PLU de V.
function V = vdm(xs)
% Matriz de van der Monde dada una lista de puntos xi
%
% Tcrea una matriz NxN cuyo elemento M[i,j] es xi^j
%
% xs: lista de puntos
%
% Devuelve un array NxN
N = length(xs);
M = zeros(N, N);
for i = [1:N]
for j = [1:N]
V(i,j) = xs(i)^(j-1);
end
end
end
xs = [-1,0,1,2,3];
V = vdm(xs)
%Otra versión más rápida (vectorizada)
function V = vdm2(xs)
% Matriz de van der Monde dada una lista de puntos xi
%
% Tcrea una matriz NxN cuyo elemento M[i,j] es xi^j
%
% xs: lista de puntos
%
% Devuelve un array NxN
N = length(xs);
M = zeros(N, N);
for j = [1:N]
V(:,j) = xs.^(j-1);
end
end
xs = [-1,0,1,2,3];
V = vdm2(xs)
[L, U, P] = lu(V)
P*L*U
Esta llamada a solve usa el método de Gauss => $\approx 2n^3/3$ operaciones.
f = @(x)(sin(x));
ys = f(xs)';
coeficientes = V\ys
% Los coeficientes son los mismos
ys = f(xs)';
vs = L\(P'*ys);
ws = U\vs;
coeficientes = ws
N = length(xs);
x_eval = linspace(min(xs),max(xs),200);
y_eval = polyval(flipud(coeficientes), x_eval);
y_fun = f(x_eval);
hold on;
plot(xs, ys, 'bo')
plot(x_eval, y_eval, 'b-')
plot(x_eval, y_fun, 'g-')
legend('nodos de interpolacion', 'polinomio interpolador', 'funcion')
hold off;
f =@(x)(1./(1+25*x.^2));
ys = f(xs)';
vs = L\(P'*ys);
ws = U\vs;
coeficientes = ws
N = length(xs);
x_eval = linspace(min(xs),max(xs),200);
y_eval = polyval(flipud(coeficientes), x_eval);
y_fun = f(x_eval);
hold on;
plot(xs, ys, 'bo')
plot(x_eval, y_eval, 'b-')
plot(x_eval, y_fun, 'g-')
legend('nodos de interpolacion', 'polinomio interpolador', 'funcion')
hold off;
Tiempo de ejecución¶
Medimos el tiempo de ejecución:
- cpu time es número de segundos*cpu empleados
- wall time es el tiempo total
La dos mediciones pueden ser muy distintas:
- uno o varios cores
- cuellos de botella: la cpu puede estar parada esperando a que le lleguen datos de la RAM
- carga del sistema, requerimientos de cada proceso en el momento de la medición
Observa los tiempos de ejecución: ¿responde a tus expectativas sobre qué proceso es más rápido y cuál más lento?
Ejercicio¶
- Aprende a medir tiempos de ejecución usando la función
timeitde matlab (documentación de timeit).
xs = linspace(-1,1,500);
f = @(x)(sin(x));
ys = f(xs)';
callV = @()(vdm(linspace(-1,1,500)));
timeit(callV)
%%time
V=vdm(linspace(-1,1,500));
%%time
V=vdm2(linspace(-1,1,500));
%%time
coeficientes = V\ys;
%%time
[P, L, U] = lu(V);
%%time
vs = L\(P'*ys);
ws = U\vs;
Ejercicio¶
- Crea una función que recibe como argumento un número
ny devuelve una matriz cuadradaMde dimensionesnxncuya entradaM[i,j]es 2 en la diagonal principal (es decir, si i=j), -1 si se trata de la diagonal encima o por debajo de la diagonal principal. $$ \left(\begin{array}{rrrrr} 2.0 & -1.0 & 0.0 & 0.0 & 0.0 \\ -1.0 & 2.0 & -1.0 & 0.0 & 0.0 \\ 0.0 & -1.0 & 2.0 & -1.0 & 0.0 \\ 0.0 & 0.0 & -1.0 & 2.0 & -1.0 \\ 0.0 & 0.0 & 0.0 & -1.0 & 2.0 \end{array}\right) $$ - Al crear la matriz, ¿has usado un bucle o dos bucles anidados (u otro mecanismo)? Piensa si puedes crear la matriz de forma más eficiente.
- Calcula la factorización PLU de la matriz anterior. ¿Qué clase de matrices son P, L y U? ¿tienen alguna propiedad especial, además de ser, respectivamente, una matriz de permutación, triangular inferior y triangular superior?
- Resuelve el sistema de ecuaciones $M\cdot \mathbf{x}= \mathbf{y}$, para $n=100$, y para el vector $\mathbf{y}$ que tiene entrada $i$ igual a $i(n-i)$. Compara el tiempo necesario para resolver el sistema usando solve, y usando la factorización PLU.