Cosas básicas sobre errores en las medidas. Notaciones y tratamiento estadístico

En mi blog de Física, estoy escribiendo un artículito sobre el tratamiento básico del error en la medida de magnitudes físicas. Os invito a echar un vistazo siguiendo este enlace.

Un ejercicio de aplicación de la factorización de Crout a la resolución de un sistema de ecuaciones lineales compatible determinado

Se quiere resolver el sistema de ecuaciones lineales $$\left\{ \begin{matrix}x_1&-&x_2&&&=&4\\ x_1&+&x_2&+&x_3&=&3 \\ x_1&&&-&x_3&=&2 \end{matrix}\right.$$, que en otro artículo ya había sido resuelto por el método de Doolittle. Ahora vamos a resolverlo por el método de Crout.

Escribamos el sistema en forma matricial $$\begin{pmatrix}1&-1&0\\1&1&1\\1&0&-1\end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}=\begin{pmatrix}4\\3\\2\end{pmatrix}$$ Expresado así, $AX=B$, donde $X=(x_1\,x_2\,x_3)^\top$, $b=(4\,3\,2)^\top$ y $A=\begin{pmatrix}1&-1&0\\1&1&1\\1&0&-1\end{pmatrix}$, vamos a factorizar la matriz $A$ (que es regular) de la forma $A=LU$ por el método de Crout, donde $L$ es una matriz triangular inferior, y $U$ es una matriz triangular superior con unos en la diagonal principal —recordemos que los unos en la digonal principal estaban en la matriz $L$ en el método de Doolittle—.

Al obtener $A=LU$, el sistema de ecuacione puede escribirse de la forma $LUX=B$, esto es, $L(UX)=b$. Denotando $UX=Y$, la resolución constará de los siguientes pasos:

1. Resolveremos $LY=B$ para determinar el vector $Y=(y_1\,y_2\,y_3)^\top$
2. Una vez conocido $Y$, resolveremos finalmente $UX=Y$ para determinar el vector $X=(x_1\,x_2\,x_3)^\top$

Cálculo de las matrices $L$ y $U$

De acuerdo con el m. de Crout, la matriz triangular inferior es de la forma $L=\begin{pmatrix}\ell_{11}&0&0\\ \ell_{21}&\ell_{22}&0\\\ell_{31}&\ell_{32}&\ell_{33}\end{pmatrix}$ y la matriz triangular superior es $U=\begin{pmatrix}1&u_{12}&u_{13}\\ 0&1&u_{23}\\0&0&1\end{pmatrix}$. Entonces, como $$\begin{pmatrix}\ell_{11}&0&0\\ \ell_{21}&\ell_{22}&0\\\ell_{31}&\ell_{32}&\ell_{33}\end{pmatrix}\begin{pmatrix}1&u_{12}&u_{13}\\ 0&1&u_{23}\\0&0&1\end{pmatrix}=\begin{pmatrix}1&-1&0\\1&1&1\\1&0&-1\end{pmatrix}$$ por la definición de producto de matrices se tiene que
  $a_{11}=1=\ell_{11}$
  $a_{21}=1=\ell_{21}$
  $a_{31}=1=\ell_{31}$

  $a_{12}=-1=\ell_{11}\,u_{12}=1\cdot u_{12}\Rightarrow u_{12}=-1$
  $a_{22}=1=\ell_{21}\,u_{12}+1\cdot \ell_{22}=1\cdot u_{12}+\ell_{22}=-1+\ell_{22}\Rightarrow \ell_{22}=2$
  $a_{32}=0=\ell_{31}\,u_{12}=\ell_{32}\,u_{12}+1\cdot \ell_{32}=1\cdot (-1)+\ell_{32}\Rightarrow \ell_{32}=1$

  $a_{13}=0=\ell_{11}\,u_{13}=1\cdot u_{13}\Rightarrow u_{13}=0$
  $a_{23}=1=\ell_{21}\,u_{13}+\ell_{22}\,u_{23}=0+\ell_{22}\,u_{23}=0+2u_{23}\Rightarrow u_{33}=\frac{1}{2}$
  $a_{33}=-1=\ell_{31}\,u_{13}+\ell_{32}\,u_{23}+1 \cdot \ell_{33}=1 \cdot 0+1\cdot u_{23}+\ell_{33}=\frac{1}{2}+\ell_{33} \Rightarrow \ell_{33}=-\frac{3}{2}$

Por tanto $L=\begin{pmatrix}1&0&0\\ 1&2&0\\1&1&-\frac{3}{2}\end{pmatrix}$ y $U=\begin{pmatrix}1&-1&0\\ 0&1&\frac{1}{2}\\0&0&1\end{pmatrix}$

Abordamos ahora el primer paso, resolviendo $LY=B$

Escribiendo el sistema de ecuaciones, $\left\{ \begin{matrix}y_1&&&&&=&4\\ y_1&+&y_2&&&=&3 \\ y_1&+&\frac{1}{2}\,y_2&+&y_3&=&2 \end{matrix}\right.$, encontramos fácilmente $y_1=4$, $y_2=-1$ y $y_3=-\frac{3}{2}$

Abordamos finalmente el segundo paso, resolviendo $UX=Y$

Escribiendo el sistema de ecuaciones, $\left\{ \begin{matrix}x_1&-&x_2&&&=&4\\ &&2x_2&+&x_3&=&-1 \\ &&&&-\frac{3}{2}x_3&=&-\frac{3}{2}\end{matrix}\right.$, de donde $x_1=3$, $x_2=-1$ y $x_3=1$.

Algoritmo de Crout

Arriba hemos realizado las operaciones paso a paso, si bien con un poco de paciencia podemos inducir las expresiones matemáticas que dan valor a los elementos de $L$ y $U$ en el caso general de tener que factorizar una matriz $A$ de orden $n$; esto lo podemos hacer partiendo de las regularidades que encontraremos para matrices de orden $3$. Obtendremos así lo que podemos entender como el algoritmo que cómodamente implementaremos mediante un lenguaje de programación, y, así, automatizaremos los cálculos. Se puede comprobar que:

• $\ell_{i1}=a_{i1}$ si $j=1$ para $i=1,\ldots,n$
• $u_{1j}=\dfrac{a_{1j}}{a_{11}}$ si $i=1$ para $j=2\,\ldots,n$
• $u_{ij}=\dfrac{a_{ij}-\displaystyle \sum_{k=1}^{j-1}\,\ell_{ik}\,u_{kj}}{\ell_{ii}}$ si $i\lt j$ para $i=2,\ldots,n$
• $\ell_{ij}=\dfrac{a_{ij}-\displaystyle \sum_{k=1}^{i-1}\,\ell_{ik}\,u_{kj}}{u_{jj}}$ si $i\ge j$ para $j=2,\ldots,n$

$\diamond$

Las matrices de permutación y los métodos de reducción con pivotaje

Las matrices de permutación permiten formalizar los procedimientos de pivotaje en la resolución de sistemas de ecuaciones lineales. Tal como ya se ha comentado en otros artículos, el pivotaje puede ser parcial (p. maximal por columnas) —en el que solamente se permutan filas de una matriz o de un sistema de ecuaciones lineales— o bien p. total (o p. completo) —en el que se permutan filas (entre el conjunto de las mismas), y también columnas (entre el conjunto de columnas)—. En este breve artículo hablaré de las operaciones de permutación maximal por columnas, que ilustraré con un sencillo ejemplo, aprovechándolo para mostrar cómo podemos utilizar esta operación para formalizar la resolución de un sistema de ecuaciones lineales en el que intercambiaremos (permutaremos) dos de sus filas.

Matrices de permutación en la pivotación maximal por columnas

Denominamos matriz de permutación a una matriz cuadrada $P$ formada a partir de la alteración del orden de las filas de la matriz identidad $I$ tal que al multiplicar por la izquierda a una matriz cuadrada del mismo orden $A$ permuta sus filas convenientemente. Así, por ejemplo, dada la matriz $A=\begin{pmatrix}0&1&2\\1&2&3\\4&-1&1\end{pmatrix}$ para intercambiar la primera y la tecera filas —siendo $I=\begin{pmatrix}1&0&0\\0&1&0\\0&0&1\end{pmatrix}$ la matriz identidad del orden de la matriz dada— realizaremos el producto de $P=\begin{pmatrix}0&0&1\\0&1&0\\1&0&0\end{pmatrix}$ —notemos que intercambiamos las filas primera y tercera de la matriz identidad; las mismas que queremos inteercambiar en la matriz $A$— por la matriz $A$ por la izquierda: $PA=\begin{pmatrix}0&0&1\\0&1&0\\1&0&0\end{pmatrix}\begin{pmatrix}0&1&2\\1&2&3\\4&-1&1\end{pmatrix}$ obteniendo como resultado la matriz pedida (con las filas primera y tercera permutadas): $\begin{pmatrix}4&-1&1\\1&2&3\\0&1&2\end{pmatrix}$.

Aplicación de las matrices de permutación para la técnica de pivotaje maximal por columnas (p. parcial) en la resolución gaussiana de un sistema de ecuaciones lineales

Consideremos ahora el sistema de ecuacions lineales $$\left\{\begin{matrix} &&x_2&+&2x_3&=&1 \\ x_1&+&2x_2&+&3x_3&=&0 \\ 4x_1&-&x_2&+&x_3&=&0 \\ \end{matrix}\right.$$ que podemos expresar matricialmente de la forma $AX=B \quad (1)$, donde $A=\begin{pmatrix}0&1&2\\1&2&3\\4&-1&1\end{pmatrix}$ es la matriz de los coeficientes del sistema; $X=(x_1\,x_2\,x_3)^\top$ es la matriz columna de las incógnitas, y $B=(1\,0\,0)^\top$ es la matriz columna de los términos independientes.

Para reducir el sistema por el método de Gauss, nos interesa permutar la primera ecuación con alguna de las otras dos, habida cuenta de que el coeficiente de la primera incógnita de la primera ecuación es $0$. Para ello, debemos elegir el elemento pivote de la primera columna, que, como norma general, escogemos el elemento máximo de la misma; por lo tanto, en el caso que nos ocupa, intercambiaremos (permutaremos) la primera y la tercera filas. Expresado en forma matricial esto se va a traducir en multiplicar por la matriz de permutación correspondiente por la izquierda a ambos miembros de (1), así que escribiremos el sistema equivalente $PAX=PB$, esto es $$\begin{pmatrix}0&0&1\\0&1&0\\1&0&0\end{pmatrix}\begin{pmatrix}0&1&2\\1&2&3\\4&-1&1\end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}=\begin{pmatrix}0&0&1\\0&1&0\\1&0&0\end{pmatrix}\begin{pmatrix}1\\0\\0\end{pmatrix}$$ y realizando las operaciones de multiplicación de matrices llegamos a $$\begin{pmatrix}4&-1&1\\1&2&3\\0&1&2\end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}=\begin{pmatrix}0\\0\\1\end{pmatrix}$$ Es decir, las ecuaciones del sistema equivalente quedan ahora en el siguiente orden $$\left\{\begin{matrix} 4x_1&-&x_2&+&x_3&=&0 \\ x_1&+&2x_2&+&3x_3&=&0 \\ &&x_2&+&2x_3&=&1 \\ \end{matrix}\right.$$ Podemos ahora realizar ya los pasos de reducción por Gauss: de $-\frac{1}{4}\,e_1+e_2 \rightarrow e_2$ obtenemos el sistema equivalente $$\left\{\begin{matrix} 4x_1&-&x_2&+&x_3&=&0 \\ &&9x_2&+&11x_3&=&0 \\ &&x_2&+&2x_3&=&1 \\ \end{matrix}\right.$$ y, finalmente, mediante la operación elemental (entre la segunda fila y la tercera del paso anterior) $-\frac{1}{9}\,e_1+e_3 \rightarrow e_3$ llegamos a $$\left\{\begin{matrix} 4x_1&-&x_2&+&x_3&=&0 \\ &&9x_2&+&11x_3&=&0 \\ &&&&\frac{7}{9}x_3&=&1 \\ \end{matrix}\right.$$ por lo que, despejando la tercera incóngita de la tercera ecuación, y sustituyendo regresivamente a las ecuaciones segunda y primera, obtenemos la solución: $x_1=-\dfrac{5}{7}$, $x_2=-\dfrac{11}{7}$ y $x_3=\dfrac{9}{7}$.

En ulteriores artículos mostraré cómo intervienen las matrices de permutación en la resolución de sistemas de ecuaciones lineales en el caso de optar por los métodos de factorización $LU$ de la matriz de los coeficientes del sistema.

Conviene tener en cuenta que toda matriz de permutación $P$ es simétrica: $P=P^\top$; y ortogonal: $P^{-1}=P^\top$, como puede demostrarse sin dificultad. Estas propiedades de las matrices de permutación facilitan la realización de determinados cálculos, tal y como veremos más adelante en este contexto de los métodos de reducción con pivotaje.$\diamond$

Un ejercicio de aplicación de la factorización de Doolittle a la resolución de un sistema de ecuaciones lineales compatible determinado

Se quiere resolver el sistema de ecuaciones lineales $$\left\{ \begin{matrix}x_1&-&x_2&&&=&4\\ x_1&+&x_2&+&x_3&=&3 \\ x_1&&&-&x_3&=&2 \end{matrix}\right.$$

Escribamos el sistema en forma matricial $$\begin{pmatrix}1&-1&0\\1&1&1\\1&0&-1\end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}=\begin{pmatrix}4\\3\\2\end{pmatrix}$$ Expresado así, $Ax=b$, donde $x=(x_1,x_2,x_3)^\top$, $b=(4,3,2)^\top$ y $A=\begin{pmatrix}1&-1&0\\1&1&1\\1&0&-1\end{pmatrix}$, vamos a factorizar la matriz $A$ (que es regular) de la forma $A=LU$ por el método de Doolittle, donde $L$ es una matriz triangular inferior con unos en la diagonal principal, y $U$ es una matriz triangular superior.

Al obtener $A=LU$, el sistema de ecuacione puede escribirse de la forma $LUX=B$, esto es, $L(UX)=b$. Denotando $UX=Y$, la resolución constará de los siguientes pasos:

1. Resolveremos $LY=b$ para determinar el vector $Y=(y_1\,y_2\,y_3)^\top$
2. Una vez conocido $Y$, resolveremos $Ux=X$ para determinar $X=(x_1\,x_2\,x_3)^\top$

Cálculo de las matrices $L$ y $U$

De acuerdo con el m. de Doolittle, la matriz triangular inferior es de la forma $L=\begin{pmatrix}1&0&0\\ \ell_{21}&1&0\\\ell_{31}&\ell_{32}&1\end{pmatrix}$ y la matriz triangular superior es $U=\begin{pmatrix}u_{11}&u_{12}&u_{13}\\ 0&u_{22}&u_{23}\\0&0&u_{33}\end{pmatrix}$. Entonces, como $$\begin{pmatrix}1&0&0\\ \ell_{21}&1&0\\\ell_{31}&\ell_{32}&1\end{pmatrix}\begin{pmatrix}u_{11}&u_{12}&u_{13}\\ 0&u_{22}&u_{23}\\0&0&u_{33}\end{pmatrix}=\begin{pmatrix}1&-1&0\\1&1&1\\1&0&-1\end{pmatrix}$$ por la definición de producto de matrices se tiene que
  $a_{11}=1=u_{11}$
  $a_{21}=1=\ell_{21}\,u_{11}=\ell_{21}\cdot 1\Rightarrow \ell_{21}=1$
  $a_{31}=1=\ell_{31}\,u_{11}=\ell_{31}\cdot 1\Rightarrow \ell_{31}=1$

  $a_{12}=-1=u_{12}$
  $a_{22}=1=\ell_{21}\,u_{12}+u_{22}=1\cdot (-1)+u_{22}\Rightarrow u_{22}=2$
  $a_{32}=0=\ell_{31}\,u_{12}=\ell_{32}\,u_{22}=1\cdot (-1)+\ell_{32}\cdot 2\Rightarrow \ell_{32}=\frac{1}{2}$

  $a_{13}=0=u_{13}$
  $a_{23}=1=\ell_{21}\,u_{13}+u_{23}=0+u_{23}\Rightarrow u_{23}=1$
  $a_{33}=-1=\ell_{31}\,u_{13}+\ell_{32}\,u_{23}+u_{33} =0+\ell_{32}\cdot 1+u_{33} \Rightarrow u_{33}=-\frac{3}{2}$

Por tanto $L=\begin{pmatrix}1&0&0\\ 1&1&0\\1&\frac{1}{2}&1\end{pmatrix}$ y $U=\begin{pmatrix}1&-1&0\\ 0&2&1\\0&0&-\frac{3}{2}\end{pmatrix}$

Abordamos ahora el primer paso, resolviendo $LY=B$

Escribiendo el sistema de ecuaciones, $\left\{ \begin{matrix}y_1&&&&&=&4\\ y_1&+&y_2&&&=&3 \\ y_1&+&\frac{1}{2}\,y_2&+&y_3&=&2 \end{matrix}\right.$, encontramos fácilmente $y_1=4$, $y_2=-1$ y $y_3=-\frac{3}{2}$

Abordamos finalmente el segundo paso, resolviendo $UX=Y$

Escribiendo el sistema de ecuaciones, $\left\{ \begin{matrix}x_1&-&x_2&&&=&4\\ &&2x_2&+&x_3&=&-1 \\ &&&&-\frac{3}{2}x_3&=&-\frac{3}{2}\end{matrix}\right.$, de donde $x_1=3$, $x_2=-1$ y $x_3=1$.

Algoritmo de Doolittle

Arriba hemos realizado las operaciones paso a paso, si bien con un poco de paciencia podemos inducir las expresiones matemáticas que dan valor a los elementos de $L$ y $U$ en el caso general de tener que factorizar una matriz $A$ de orden $n$; esto lo podemos hacer partiendo de las regularidades que encontraremos para matrices de orden $3$. Obtendremos así lo que podemos entender como el algoritmo que cómodamente implementaremos mediante un lenguaje de programación, y, así, automatizaremos los cálculos. Se puede comprobar que:

• $u_{1j}=a_{1j}$ si $i=1$ y para $j=1,\ldots,n$
• $\ell_{i1}=\dfrac{a_{i1}}{u_{11}}$ si $j=1$ y para $i=2,\ldots,n$
• $\ell_{ij}=\dfrac{a_{ij}-\displaystyle \sum_{k=1}^{j-1}\,\ell_{ik}\,u_{kj}}{u_{jj}}$ si $i\gt j$ para $i=2,\ldots,n$
• $u_{ij}=\dfrac{a_{ij}-\displaystyle \sum_{k=1}^{i-1}\,\ell_{ik}\,u_{kj}}{u_{jj}}$ si $i\le j$ para $j=2,\ldots,n$

$\diamond$

Factorización gaussiana LU de una matriz regular

Preliminares

Toda matriz de orden $n$, regular, $A$, cuyos determinantes principales sean distintos de cero, puede expresarse de la forma $A=LU$ mediante operaciones de reducción gaussiana sin necesidad de pivotaje, donde $L$ es una matriz triangular inferior de orden $n$ con unos en la diagonal principal y $U$ es una matriz triangular superior de orden $n$, siendo ambas matrices regulares (sus determinantes son no nulos).

En estas condiciones, las matrices $L$ y $U$ son únicas. En efecto, supongamos que no lo sean, esto es, si podemos escribir $A=L_1\,U_1=L_2\,U_2$; entonces, como todas estas matrices son inversibles, podemos escribir $L_{2}^{-1}\,L_1\,U_1=L_{2}^{-1}\,L_2\,U_2$, luego $L_{2}^{-1}\,L_1\,U_1=I\,U_2$ y por tanto $L_{2}^{-1}\,L_1\,U_1=U_2$; multiplicando ahora por $U_{2}^{-1}$ por la derecha en ambos miembros, tenemos $L_{2}^{-1}\,L_1\,U_1\,U_{2}^{-1}=U_2\,U_{1}^{-1}$, esto es $L_{2}^{-1}\,L_1\,U_1\,U_{1}^{-1}=U_2\,U_{1}^{-1}$ y por tanto, $L_{2}^{-1}\,L_1\,I=U_2\,U_{1}^{-1}$, que es lo mismo que $L_{2}^{-1}\,L_1=U_2\,U_{1}^{-1}$. La matriz inversa de una matriz triangular inferior es otra matriz triangular inferior, y lo mismo sucede con la inversa de una matriz triangular superior, así que, como el producto de dos matrices triangulares inferiores es otra matriz triangular inferior (primer miembro) y el producto de dos matrices triangulares superiores es una matriz triangular superior, la única manera de que esto sea posible es que ambos miembros de la igualdad sean la matriz identidad, por consiguiente: $L_{2}^{-1}\,L_1=I \Leftrightarrow L_2=L_1$ y $U_2\,U_{1}^{-1}\Leftrightarrow U_2=U_1$. $\square$

Observaciones:

• En una futura entrada trataré el caso en que sea necesario el pivotaje por el hecho de encontrar elementos nulos en la diagonal de alguno de los pasos de reducción y que supondrían divisiones por cero —sería el caso en el que no todos los determiantes principales sea distinto de cero, incumpliendo la condición dada arriba—, o bien números muy pequeños en algún elemento de la diagonal que harían que los errores de redondeo se amplificasen enormemente. Trataremos pues más adelante los casos de pivotaje, bien sea pivotaje parcial (también llamado maximal por columnas) en el que hay que permutar filas de la matriz; o bien se trate de pivotaje total (también llamado pivotaje completo), en el que no sólo hay que permutar filas sino también columnas. En los casos de pivotaje veremos que intervendran las llamadas matrices de permutación.
• Si no son necesarias operaciones de pivotaje, hay otros métodos para expresar una matriz $A$ como producto de una matriz triangular inferior por una matriz triangular superior, pero que no se basan en la reducción gaussiana:
  • El método de Doolittle permite una factorización de este tipo, obteniendo la matriz triangular inferior $L$ con unos en la diagonal principal, sin utilizar algoritmos de reducción (véase este ejemplo práctico).
  • El método de Crout permite también una factorización similar, obteniendo la matriz triangular superior $U$ con unos en la diagonal principal, sin utilizar algoritmos de reducción (véase este ejemplo práctico)
• Además de la factorización de $A$ mediante el producto de una matriz triangular inferior por una matriz triangular superior hay otras formas de factorización:
  • Si $A$ es una matriz simétrica, puede realizarse una factorización más conveniente de la forma $A=LDL^\top$, donde $L$ es una matriz triangular inferior con unos en la diagonal y $D$ es una matriz diagonal.
  • Si $A$ es una matriz definida positiva, el criterio de Sylvester asegura que los elementos de la diagonal, $d_{kk}$, son positivos. Y definiendo $\mathcal{L}:=L\,D^{\frac{1}{2}}$, con $D^{\frac{1}{2}}=\text{diagonal}(d_{11}^{\frac{1}{2}},\ldots,d_{nn}^{\frac{1}{2}})$, entonces la matriz $A$ puede expresarse de la forma $A=\mathcal{L}\,\mathcal{L}^\top$, que se denomina factorización de Cholesky.
    • Si $A$ es simétrica y definida positiva, entonces la factorización de Cholesky de $A$ que se acaba de exponer es más sencilla, pues $\mathcal{L}=L$. Así, $A=L\,L^\top$; en otras palabras, en tales condiciones la matriz triangular superior $U$ de la factorización $LU$ es $U=L^\top$

-oOo-
Ensayemos ahora un poco con una matriz regular como, por ejemplo, $A=\begin{pmatrix}3&-1&1\\1&1&1\\2&1&0\end{pmatrix}$, cuya reducción no requiere pivotamiento parcial pues no nos vamos a encontrar con ceros en la diagonal principal de la matrices que iremos obteniendo.

Más adelante, extenderemos lo que vamos a deducir ahora al caso en que sí se requiere pivotamiento.

Empezaremos pues con el proceso de reducción. Siendo esta matriz de orden $3$, tendremos que llevar a cabo dos etapas.

Obtenemos ceros por debajo de la columna del elemento $a_{11}$ haciendo las siguientes transformaciones $\begin{matrix}f_1 \rightarrow f_1 \\ -\frac{1}{3}f_1+f_2\rightarrow f_2\\-\frac{2}{3}f_1+f_3\rightarrow f_3\end{matrix}$, obteniendo la matriz $\begin{pmatrix}3&-1&1\\0&\frac{4}{3}&\frac{2}{3}\\0&\frac{5}{3}&-\frac{2}{3}\end{pmatrix}$. Finalmente, anularemos el elemento de la segunda fila y la segunda columna mendiante $\begin{matrix}f_1 \rightarrow f_1 \\f_2\rightarrow f_2\\-\frac{5}{4}f_2+f_3\rightarrow f_3\end{matrix}$ lo cual da lugar a la matriz reducida por Gauss $\begin{pmatrix}3&-1&1\\0&\frac{4}{3}&\frac{2}{3}\\0&0&-\frac{3}{2}\end{pmatrix}$, que, obviamente es una matriz triangular superior, y a la que denominaremos $U$.

Utilidades de la factorización $A=LU$

La factorización de una matriz regular como producto de una matriz triangular inferior y una matriz triangular superior facilita la resolución de un sistema de ecuaciones lineales, así como diversos cálculos matriciales. Veamos tres utilidades importantes:

• Dado un sistema de ecuaciones lineales $Ax=b$, tenemos que, al factorizar $A$ de este modo podemos escribir $LUx=b$, esto es $L(Ux)=b$, con lo cual podemos resolver el sistema en dos pasos: $\left\{\begin{matrix}\text{Primer paso: resolvemos}\quad Ly=b\\\text{Segundo paso: conocido}\quad y \quad \text{resolvemos:} \quad Ux=y \quad \text{y acabamos calculando} \quad x\end{matrix}\right.$
• Cálculo de la matriz inversa de $A$: Como $A=LU$, se tiene que $A^{-1}=(LU)^{-1}=U^{-1}\,L^{-1}$
• Cálculo del determinante de $A$: Como $A=LU$, y teniendo en cuenta que el determinante de un producto de matrices es igual al producto de sus determinantes se tiene que $\text{det}(A)=\text{det}(L)\cdot \text{det}(U)$. Entonces, si los elementos de la diagonal de $L$ son unos (factorización de Doolittle), resulta que $\text{det}(A)=1\cdot \displaystyle \prod_{i=1}^{n}\,u_{ii}=\displaystyle \prod_{i=1}^{n}\,u_{ii}$; y en el caso de que los unos estén en la diagonal de $U$ (factorización de Crout), tendremos que $\text{det}(A)=\displaystyle \prod_{i=1}^{n}\,\ell_{ii}\cdot 1=\displaystyle \prod_{i=1}^{n}\,\ell_{ii}$

Formalicemos un poco lo que acabamos de hacer acerca de cómo calcular la factorización $A=LU$:

Notemos ahora que, de los multiplicadores de la transformación, podemos escribir esta otra matriz triangular inferior y a la que denominaremos $L$ y escribiremos $L=\begin{pmatrix}1&0&0\\\frac{1}{3}&1&0\\\frac{2}{3}&\frac{5}{4}&1\end{pmatrix}$, donde hemos cambiado el signo de los números distintos de cero y de uno que hemos utilizado para efectuar las transformaciones. Pues bien, puede comprobarse que haciendo el producto $LU$ se obtiene la matriz $A$; en efecto $$\begin{pmatrix}1&0&0\\\frac{1}{3}&1&0\\\frac{2}{3}&\frac{5}{4}&1\end{pmatrix}\begin{pmatrix}3&-1&1\\0&\frac{4}{3}&\frac{2}{3}\\0&0&-\frac{3}{2}\end{pmatrix}=\begin{pmatrix}3&-1&1\\1&1&1\\2&1&0\end{pmatrix}$$ Esta es pues una manera sencilla de efectuar la factorización de una matriz regular dada $A$ como producto de una matriz triangular inferior por una matriz triangular superior, $A=LU$, siempre que no tengamos que utilizar técnicas de pivotamiento parcial —por encontrarnos con ceros en los elementos de la diagonal principal en las etapas del proceso de reducción, en cuyo caso deberemos extender esta técnica, tal y como ya trataré más adelante artícul—. Veremos también las utilidades de dicha factorización en el próximo artículo.

Cálculo de la matriz triangular inferior $L$

Démonos cuenta de que la matriz obtenida $E$ (la matriz de multiplicadores) puede a su vez descomponerse como un producto de dos matrices que corresponden a las dos etapas de que consta el proceso de reducción de esta matriz de orden $3$ y que vamos a denotar como $E(1)=\begin{pmatrix}1&0&0\\-\frac{1}{3}&1&0\\-\frac{2}{3}&0&1\end{pmatrix}$ y $E(2)=\begin{pmatrix}1&0&0\\0&1&0\\0& -\frac{5}{4}&1\end{pmatrix}$, y sus inversas, $E(1)^{-1}=\begin{pmatrix}1&0&0\\\frac{1}{3}&1&0\\\frac{2}{3}&0&1\end{pmatrix}$ y $E(2)^{-1}=\begin{pmatrix}1&0&0\\0&1&0\\0&\frac{5}{4}&1\end{pmatrix}$ ;donde, si recordemos lo que hemos hecho, y con vistas a formalizar un método para poder hacer descomposiciones de matrices de orden cualquiera, $E(1)=\begin{pmatrix}1&0&0\\-\frac{a_{21}^{(1)}}{a_{11}^{(1)}}&1&0\\-\frac{a_{31}^{(1)}}{a_{11}^{(1)}}&0&1\end{pmatrix}$ y $E(2)=\begin{pmatrix}1&0&0\\0&1&0\\0& -\frac{a_{32}^{(2)}}{a_{22}^{(2)}}&1\end{pmatrix}$, y sus inversas —se demuestra fácilmente que la matriz inversa $E(j)^{-1}$ se obtiene cambiando simplemente el signo de los elementos que corresponden a los multiplicadores de la matriz $E(j)$—, esto es, $E(1)^{-1}=\begin{pmatrix}1&0&0\\ \frac{a_{21}^{(1)}}{a_{11}^{(1)}}&1&0\\ \frac{a_{31}^{(1)}}{a_{11}^{(1)}}&0&1\end{pmatrix}$ y $E(2)^{-1}=\begin{pmatrix}1&0&0\\0&1&0\\0& \frac{a_{32}^{(2)}}{a_{22}^{(2)}}&1\end{pmatrix}$ (los superíndices —entre paréntesis— indican el orden de la etapa de reducción), con lo cual, $L=E(1)^{-1}\,E(2)^{-1}=\begin{pmatrix}1&0&0\\\frac{1}{3}&1&0\\\frac{2}{3}&0&1\end{pmatrix}\begin{pmatrix}1&0&0\\0&1&0\\0&\frac{5}{4}&1\end{pmatrix}=\begin{pmatrix}1&0&0\\\frac{1}{3}&1&0\\\frac{2}{3}&\frac{5}{4}&1\end{pmatrix}$.

Cálculo de la matriz triangular superior $U$

Observemos que, denotando por $A(1)\equiv A$, al ser $A$ una matriz de orden $3$, nos encontramos con que terminamos obteniendo la matriz triangular superior $U$ en dos etapas $\left\{\begin{matrix}A(2)=E(1)\,A(1)\\ U \equiv A(3)=E(2)\,A(2)\end{matrix}\right\}$, y terminamos. Si en lugar de simplemente $2$ etapas, tuviésemos $n-1$ (caso de una matriz de orden $n$), seguieremos la recurrencia $A(j+1)=E(j)\,A(j)$ para $j=1,2,\ldots,n-1$, siendo —recordémoslo— $A(1)\equiv A$ y $A(n)\equiv U$.

-oOo-

Observación: Conociendo $L$, para calcular $U$, fijémonos en que de la factorización que buscamos, $A=LU$, tenemos que, multiplicando ambos miembros por $L^{-1}$ por la izquierda, se tiene que $U=L^{-1}\,A=(E(1)^{-1}\,E(2)^{-1})^{-1}\,A=(E(2)^{-1})^{-1}\,(E(1)^{-1})^{-1}\,A=E(2)\,E(1)\,A$, con lo cual la matriz triangular superior $U$ es, como debe ser —valga como comprobación—, $$U=\begin{pmatrix}1&0&0\\0&1&0\\0&-\frac{5}{4}&1\end{pmatrix}\begin{pmatrix}1&0&0\\-\frac{1}{3}&1&0\\-\frac{2}{3}&0&1\end{pmatrix}\begin{pmatrix}3&-1&1\\1&1&1\\2&1&0\end{pmatrix}=\begin{pmatrix}3&-1&1\\0&\frac{4}{3}&\frac{2}{3}\\0&0&-\frac{3}{2}\end{pmatrix}$$

-oOo-

Así, en el caso de factorizar una matriz de orden $4$, la matriz triangular inferior haremos el siguiente cálculo, $L=E(1)\,E(2)\,E(3)$; y para calcular la matriz triangular superior $U$ tendríamos tres pasos en lugar de dos: $\left\{\begin{matrix} A(2)=E(1)\,A(1)\\ A(3)=E(2)\,A(2) \\ U\equiv A(4)=E(3)\,A(3) \end{matrix}\right\}$ donde $A(1) \equiv A$, con las respectivas matrices de multiplicadores (ahora hay tres): $E(1)=\begin{pmatrix}1&0&0&0 \\ -\frac{a_{21}^{(1)}}{a_{11}^{(1)}} & 1 & 0 & 0 \\ -\frac{a_{31}^{(1)}}{a_{11}^{(1)}} & 0 & 1 & 0 \\ -\frac{a_{41}^{(1)}}{a_{11}^{(1)}} & 0 & 0 & 1 \end{pmatrix}$, $E(2)=\begin{pmatrix}1& 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & -\frac{a_{32}^{(2)}}{a_{22}^{(2)}} & 1 & 0 \\ 0 & -\frac{a_{42}^{(2)}}{a_{22}^{(2)}} & 0 & 1 \end{pmatrix}$ y $E(3)=\begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & -\frac{a_{43}^{(3)}}{a_{33}^{(3)}} & 1\end{pmatrix}$, así como sus respectivas inversas: $E(1)^{-1}=\begin{pmatrix}1&0&0&0 \\ \frac{a_{21}^{(1)}}{a_{11}^{(1)}} & 1 & 0 & 0 \\ \frac{a_{31}^{(1)}}{a_{11}^{(1)}} & 0 & 1 & 0 \\ \frac{a_{41}^{(1)}}{a_{11}^{(1)}} & 0 & 0 & 1 \end{pmatrix}$, $E(2)^{-1}=\begin{pmatrix}1& 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & \frac{a_{32}^{(2)}}{a_{22}^{(2)}} & 1 & 0 \\ 0 & \frac{a_{42}^{(2)}}{a_{22}^{(2)}} & 0 & 1 \end{pmatrix}$ y $E(3)^{-1}=\begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & \frac{a_{43}^{(3)}}{a_{33}^{(3)}} & 1\end{pmatrix}$

Generalización. Factorización gaussiana $LU$ de una matriz de orden $n$ (sin que sea necesario el pivotaje)

A partir de lo que se acaba de decir, resulta ya muy fáoil generalizar lo que se ha encontrado para una matriz $A$ de orden $n$ —siempre que no haya que recurrir a permutar filas (pivotamiento) por encontrarnos con ceros o bien con elementos próximos a cero en las posiciones diagonales—, podemos escribir:

1. $L=E(1)^{-1}\,E(2)^{-1}\,\ldots\,E(n-1)^{-1}$
2. $U\equiv A(n)=E(n-1)\,A(n-1)=E(n-1)\,E(n-2)\,A(n-2)$
   $=\ldots=E(n-1)\,E(n-2)\,\ldots\,E(2)\,\,E(1)\,A(1)$, donde $A(1)\equiv A$

donde $$E(j)=\begin{pmatrix} 1&0&\ldots &0&\ldots &0&0&0\\ 0&1&\ldots &0&\ldots &0&0&0\\ \vdots & \vdots & \ddots &\vdots &\vdots &\vdots &\vdots &\vdots\\ \vdots & \vdots & \vdots &1 &\vdots &\vdots &\vdots &\vdots\\ 0 & 0 & \ldots &-\frac{a_{j+1,j}^{(j)}}{a_{jj}^{(j)}} &1 &\ldots &\ldots &0\\ \vdots & \vdots & \vdots &\vdots &\vdots &\ddots &\vdots &\vdots\\ 0 & 0 & \ldots &-\frac{a_{n-1,j}^{(j)}}{a_{jj}^{(j)}} &\ldots &\ldots &1 &0\\ 0 & 0 & \ldots &-\frac{a_{n,j}^{(j)}}{a_{jj}^{(j)}} &\ldots &\ldots &0 &1\\ \end{pmatrix}$$ y su respectiva matriz inversa $$E(j)^{-1}=\begin{pmatrix} 1&0&\ldots &0&\ldots &0&0&0\\ 0&1&\ldots &0&\ldots &0&0&0\\ \vdots & \vdots & \ddots &\vdots &\vdots &\vdots &\vdots &\vdots\\ \vdots & \vdots & \vdots &1 &\vdots &\vdots &\vdots &\vdots\\ 0 & 0 & \ldots &\frac{a_{j+1,j}^{(j)}}{a_{jj}^{(j)}} &1 &\ldots &\ldots &0\\ \vdots & \vdots & \vdots &\vdots &\vdots &\ddots &\vdots &\vdots\\ 0 & 0 & \ldots &\frac{a_{n-1,j}^{(j)}}{a_{jj}^{(j)}} &\ldots &\ldots &1 &0\\ 0 & 0 & \ldots &\frac{a_{n,j}^{(j)}}{a_{jj}^{(j)}} &\ldots &\ldots &0 &1\\ \end{pmatrix}$$ de donde se obtiene la matriz triangular inferior $$L=\begin{pmatrix} 1&0&\ldots & 0 &\ldots & 0 & 0 & 0 \\ \frac{a_{21}^{(1)}}{a_{11}^{(1)}}&1&\ldots & 0 &\ldots & 0 & 0 & 0 \\ \vdots&\vdots&\ddots & \vdots &\vdots & \vdots & \vdots & \vdots \\ \vdots&\vdots&\vdots & 1 &\vdots & \vdots & \vdots & \vdots \\ \frac{a_{j1}^{(1)}}{a_{11}^{(1)}}& \frac{a_{j2}^{(2)}}{a_{22}^{(2)}}&\ldots & \frac{a_{j+1,j}^{(j)}}{a_{jj}^{(j)}} & 1 & \ldots & \ldots & 0 \\ \vdots&\vdots&\vdots & \vdots &\vdots & \ddots & \vdots & \vdots \\ \frac{a_{n-1,1}^{(1)}}{a_{11}^{(1)}}& \frac{a_{n-1,2}^{(2)}}{a_{22}^{(2)}}&\ldots & \frac{a_{n-1,j}^{(j)}}{a_{jj}^{(j)}} & \ldots & \ldots & 1 & 0 \\ \frac{a_{n,1}^{(1)}}{a_{11}^{(1)}}& \frac{a_{n,2}^{(2)}}{a_{22}^{(2)}}&\ldots & \frac{a_{n,j}^{(j)}}{a_{jj}^{(j)}} & \ldots & \ldots & \frac{a_{n,n-1}^{(n-1)}}{a_{n-1,n-1}^{(n-1)}} & 1 \\ \end{pmatrix}$$ $\diamond$

-oOo-

Referencias:
[1] A. Aubanell, A. Benseny y A. Deslshams, Útiles de cálculo numérico (Labor, Barcelona, 1993).
[2] C. Moreno, Introducción al cálculo numérico (UNED, Madrid, 2011).

Métodos gaussianos de resolución de sistemas de ecuaciones lineales $AX=B$. Acerca de cómo minimizar errores de redondeo: técnicas de pivotaje parcial y total

Al utilizar los métodos de reducción (m. gaussianos) —m. de Gauss, y m. de Gauss-Jordan— y se llegue a obtener algún elemento pivote nulo, $a_{ii}^{\ast}=0$, es claro que será necesario intercambiar ecuaciones (véase [1], pp. 165-166). También deberá hacerse eso cuando, los elementos de la diagonal sean números próximos a cero, para evitar la amplificación de los múltiples errores de redondeo —utilizando, claro está, un número de dígitos limitado— (véase [2], pp. 79-81). Para todos estos casos, habrá que elegir siempre como pivote el mayor elemento en valor absoluto de la columna que va a reducirse. Hay dos técnicas de pivotaje:

• Pivotaje parcial (o pivotaje maximal por columnas): Esta técnica consiste en seleccionar como elemento pivote el elemento de la columna implicada que esté por debajo de la diagonal y tenga el mayor valor absoluto, para, en el caso de ser necesario, intercambiar las ecuaciones antes de efectuar la operación de reducción.
• Pivotaje total (o pivotaje completo): Esta otra técnica consiste en seleccionar como elemento pivote el elemento máximo, intercambiando no sólo las ecuaciones sino también las incógnitas que sean necesarias antes de efectuar la operación de reducción.$\diamond$

-oOo-

Referencias:
[1] F. García y A. Nevot, Ejercicios resueltos de cálculo numérico (Paraninfo, Madrid, 1992).
[2] J.C. Mason, Métodos matriciales (Anaya, Madrid, 1986).

Las matrices y el cálculo numérico

Introducción

Las matrices tienen un relevante papel en el cálculo (y análisis) numérico, pues muchos algoritmos se fundamental en ellas, especialmente los algoritmos de resolución numérica de sistemas de ecuaciones lineales que podemos escribir de la forma $Ax=b$, donde en general, $A$ es una matriz de una aplicación lineal de un endomorfismo $f:V \rightarrow V$ sobre un cuerpo $\mathbb{K}$ —que puede ser $\mathbb{R}$ o bien $\mathbb{C}$— de tamaño $m \times n$; $x$ es el vector columna de las incógnitas (de tamaño $n\times 1$), y $b$ es el vector columna de los términos independientes (de tamaño $m\times 1$). Existen dos clases de métodos numéricos para resolver sistemas de ecuaciones lineales: los métodos directos, entre los que son muy importantes los derivados del método de Gauss —el que se estudia en el bachillerato, junto con el método de Cramer o de los determinantes— y a los que se les da el nombre de métodos gaussianos, como los basados en la factorización $A=LU$ (m. de Doolittle, m. de Crout, y m. de Cholesky); los que se basan en la factorización $A=QR$ (m. de Gram-Schmidt, m. de Householder), y los iterativos (por contraposición a los métodos directos) como son el m. de Gauss-Seidel, el m. de Jacobi, y el m. de sobrerrelajación. La razón de ser de estos métodos numéricos es la de procurar la eficiencia —claramente, por ejemplo, el método de Cramer no es eficiente— así como la estabilidad del mismo en la resolución de un problema, especialmente cuando la matriz $A$ es de gran tamaño.

Definiciones preliminares

Veamos primero algunos tipos de matrices para afianzar los conceptos a la hora de ponernos a estudiar métodos numéricos matriciales, a la vez que iremos haciendo diversos comentarios relacionados con la aplicación de las matrices a los métodos numéricos de resolución de sistemas de ecuaciones lineales, cálculo de determinantes, cálculo de la matriz inversa asociada a una matriz regular, y cálculo de valores.

Los elementos de una matriz $A\in \mathcal{M}_{m,n}$ (espacio vectorial de las matrices en un cuerpo, que, en general, será $\mathbb{C}$, de $m$ filas y $n$ columnas) se denotan por $m_{ij}$; para $i=1,\ldots,m$ y $j=1,\ldots,n$

Dada una matriz $A\in \mathcal{M}_{m,n}$ compleja, designamos por $A^{\top}$ a su matriz traspuesta; por $\bar{A}$ su matriz conjugada, y por $A^{\ast}$ (en algunos libros se denota por $A^\dagger$) a la matriz adjunta de $A$, que es la traspuesta de la matriz conjugada, esto es $A^{\ast}\equiv(\bar{A})^{\top}=\overline{A^{\top}}$.

Propiedades muy útiles en los cálculos con matrices rectangulares

1. $(A+B)^\top=A^\top+B^\top$
2. $\overline{A+B}=\bar{A}+\bar{B}$
3. $(AB)^\top=B^\top\,A^\top$
4. $\overline{AB}=\bar{A}\,\bar{B}$

Hablemos de las matrices cuadradas de orden $n$

Para matrices cuadradas, $m=n$, designamos por $I$ a la matriz identidad, donde sus elementos son $\delta_{ij}$ (delta de Kronecker), es decir, $\delta_{ij}=\left\{\begin{matrix}1&\text{si}&i=j\\0&\text{si}&i\neq j \end{matrix}\right.$

Dada una matriz cuadrada $A$ de tamaño $m \times m$, llamamos matrices principales $(A)_k$ de dicha matriz a las submatrices $k\times k$ formadas por las primeras $k$ filas y las primeras $k$ columnas. Para designar a los determinantes de esas matrices principales (menores principales) escribiremos $\Delta_k \equiv \text{det}(A)_k$

Es importante distinguir entre los diversos tipos que puede presentar una matriz cuadrada $A$:

• $A$ es singular si $\text{det}\,A=0$
• $A$ es regular si $\text{det}\,A\neq 0$ (si una matriz es regular, entonces tiene asociada una matriz inversa $A^{-1}$ que es única, de manera que $AA^{-1}=A^{-1}A=I$), y lo recíproco también es cierto. Algunas proposiciones importantes en el cálculo numérico en relación a las matrices invesibles son las siguientes:
  • Toda matriz inversible admite una factorización $A=QR$, donde $Q$ es una matriz ortogonal ($Q=Q^\top$) y $R$ es una matriz triangular superior. Para obtener la matriz $Q$ deberemos obtener una base ortogonal. Dos métodos importantes para ortogonalizar son los de Gram-Schmidt y el de Householder.
• $A$ es hermitiana si $A^{\ast}=A$. En el caso de que el cuerpo sea el de los números reales, esto se traduce en $A^{\top}=A$ y decimos que $A$ es una matriz simétrica.
• $A$ es unitaria si $A^{-1}=A^{\ast}$. En el caso de que el cuerpo sea el de los números reales, esto se traduce en $A^{-1}=A^{\top}$ y decimos que $A$ es una matriz ortogonal.
• $A$ es una matriz banda $(p,q)$ si $a_{ij}=0$ para $i\ge j+p$ ó $j \ge i+q$ y puede ser:
  • $A$ Hessenberg superior si $p=2$ y $q=n$.
  • $A$ Hessenberg inferior si $p=n$ y $q=2$.
  • $A$ triangular superior si $p=1$ y $q=n$.
  • $A$ triangular inferior si $p=n$ y $q=1$.
  • $A$ es: diagonal, tridiagonal, cuatridiagonal, pentadiagonal, ... si $p=q=1,2,3,5,\ldots$, respectivamente.
• $A$ es definida positiva si es hermitiana ($A^{\ast}=A$) y $x^{\ast}Ax\gt 0$ para todo vector $x\neq 0$.Algunas propiedades muy útilies de las matrices definidas positivas son las siguientes:
  • Recordemos el importante criterio de Sylvester que proporciona una condición necesaria y suficiente para que una matriz sea definida positiva (véase [5]), que es válido únicamente para matrices simétricas (en el caso de que $\mathbb{K}=\mathbb{R}$) o hermitianas (en el caso de que $\mathbb{K}=\mathbb{C}$). Así, para una matriz simétrica (definida sobre el cuerpo de los números reales), $A^\top=A$, decimos que es definida positiva si y sólo si todos sus menores principales —los menores principales de una matriz cuadrada son los determinantes de las submatrices principales de orden $k\le n$, esto es, los menores de orden inferior o igual a $n$ cuyos elementos son simétricos con respecto de la diagonal principal de la matriz $A$— son positivos: $\Delta_k\gt 0$, $\forall\, k=1,\ldots,n$.
  • Si una matriz $A$ es definida positiva, entonces $A$ es regular, y por tanto inversible.
  • Si $A$ es una matriz definida positiva, entonces puede resolverse un sistema de ecuaciones lineales $Ax=b$ por el método de Gauss sin intercambiar ni ecuaciones ni incógnitas; una consecuencia de ello es que $A$ admite una descomposición como el producto de una matriz triangular inferior $L$ y una matriz triangular superior $U$, de la forma $A=LU$.
  • Una matriz definida positiva y simétrica admite una factorización $A=LL^\top$ (factorización de Cholesky), donde $L$ es una matriz triangular inversible. Lo recíproco también es cierto.
• $A$ es estrictamente diagonal dominante si $|a_{ii}|\gt \displaystyle \sum_{j=i+1}^{n}\,|a_{ij}|$, para $i=1,\dots,n$. Algunas de las propiedades de estas matrices son las siguientes:
  1. Si $A$ es una matriz estrictamente diagonal dominante, entonces $A$ es regular, y por tanto inversible.
  2. Si $A$ es una matriz estrictamente diagonal dominante, entonces puede resolverse un sistema lineal $Ax=b$ por Gauss, sin intercambiar ni ecuaciones ni incóngitas, con lo cual $A$ admite una descomposición $A=LU$
• Si los determinantes principales (menores principales) $\Delta_k$ de $A$ no se anulan, entonces $A$ admite una factorización $A=LU$
• Diremos que dos matrices $A$ y $B$ son semejantes si existe una matriz inversible $C$ —a la que denominaremos matriz de transformación de semejanza de $A$ y $B$— tal que se verifique $B=C^{-1}AC$. En el caso de que $A$ sea semejante a una matriz diagonal, diremos que $A$ es diagonalizable; y si se puede encontrar una por cajas de Jordan (por cajas) que sea semejante a $A$, diremos que se puede descomponer por Jordan.
  • Lema de Shur. Toda matriz $A$ es semejante a una matriz triangular superior mediante una transformación de semejanza con una matriz unitaria.

Subrayemos también estas propiedades referidas al determinante y a la matriz inversa

1. $\text{det}\,A^\top=\text{det}\,A$
2. $\text{det}\,AB=\text{det}\,A\cdot\text{det}\,B$
3. En general, $\text{det}\,A+B\neq\text{det}\,A+\text{det}\,B$
4. $(AB)^{-1}=B^{-1}\,A^{-1}$

Acerca de los valores y vectores propios de una matriz cuadrada $A$ de orden $n$

El cálculo de los vectores propios de una matriz es muy importante en los métodos numéricos matriciales —los valors propios son los ceros del polinomio característico de una matriz $A$, y los vectores propios son los vectores asociados a los valores propios&mdash. Para calcularlos, se utilizan métodos como el de Krylov y el método de Givens, entre otros que derivan de los métodos para la resolución de sistemas de ecuaciones lineales.

Dada una matriz $A$ de orden $n$, decimos que $v\in V\equiv \mathbb{R}^n$ (espacio vectorial con el que se construye la aplicación lineal) es un vector propio de $A$ de valor propio $\lambda$ cuando se cumple que $A\,v=\lambda\,v$, condición que puede escribirse de la forma $(A-\lambda\,I)\,v=0_V$, y tal cosa es cierta si y sólo si se cumple que su valor propio $\lambda$ satiface la ecuación característica $p_{A}(\lambda)\equiv \text{det}\,(A-\lambda\,I)=0$. El primer miembro de esa igualdad es un polinomio al cual denominamos polinomio característico de $A$. Un concepto importante en el análisis numérico con métodos matriciales es el de radio espectral de $A$, que se suele designar como $\rho(A):=\text{máx}\,\{|\lambda_j|\}$, con $j=1, \ldots,r$ valores propios de $A$, es decir, se define como el módulo máximo de todos los valores propios de $A$. Se demuestra que el radio espectral es menor o igual que el valor de cualquier norma natural de $A$, $\rho(A)\le \left\|A\right\|_{.}$. Es sabido que los métodos iterativos para el cálculo de los valores y vectores propios, así como para resolver sistemas de ecuaciones lineales, y —tal como se demuestra a partir de resultados que todavía no se han expuesto— convergen si $\rho(A)\lt 1$, lo cual equivale a $\displaystyle \lim_{n\,\rightarrow\, \infty}\,A^n=0_{matriz}$, así como a $\displaystyle \lim_{n\,\rightarrow\, \infty}\,\left\|A^n\right\|_{.}=0$. Acerca de los vectores y valores pripios de $A$ podemos destacar las siguientes propiedades:

1. Los vectores propios asociados a valores propios distintos son linealmente independientes, luego si los valores propios de una matriz $A$ son todos distintos ésta es diagonalizable ($A$ es semejante a una matriz diagonal).
2. Dos matrices semejantes $A$ y $B$ tienen los mismos valores propios, y se cumple que $v$ es un vector propio de valor propio $\lambda$ de $A$ si y sólo si $S^{-1}v$ es un vector propio de valor propio $\lambda$ de $B$, donde $S$ es la matriz de la transformación de semejanza $B=S^{-1}AS$.
3. Si $A$ es semejante a una matriz diagonal $D$ por una transformación de semejanza $S$, esto es $D=S^{-1}AS$, entonces los elementos de la diagonal de $D$ son los valores propios de $A$, y los vectores asociados a las columnas de la matriz de la transformación de semejanza $S$ constituyen una base de vectores propios de $A$ (vectores propios por la derecha), y los vectores asociados a las filas de la matriz $P^{-1}$ constituyen una base de vectores propios por la izquierda de $A$.
4. Las matrices $A$ y $A^\top$ tienen los mismos valores propios.
5. Los vectores propios por la derecha (los vectores propios de $A$) son ortogonales a los vectores propios por la izquierda (los vectores propios de $A^\top$) para valores propios diferentes.
6. Una matriz $A$ es regular (inversible) si y sólo si todos sus valores propios son distintos de cero; y, en ese caso, los valores propios de $A^{-1}$ son los recíprocos de los de $A$; si $v$ es un vector propio de valor propio $\lambda$ de $A$ si y sólo si $v$ es un vector propio de valor propio $\dfrac{1}{\lambda}$ de $A^{-1}$.
7. Sea $p$ es un polinomio no nulo, entonces $\lambda$ es un valor propio de $A$ si y sólo si $p(\lambda)$ es un valor propio de $p(A)$. Y $v$ es un vector propio de valor propio $\lambda$ de $A$ si y sólo si es también vector propio de valor propio de $p(\lambda)$ de $p(A)$.
8. Sea una matriz $A$ por bloques de matrices cuadradas, de la forma $A=\begin{pmatrix}A_{11}&A_{12}\\0&A_{22}\end{pmatrix}$. Entonces, el conjunto de valores propios de $A$ es la unión de los conjuntos de valores propios de las submatrices cuadradas de la diagonal por bloques, $A_{11}$ y $A_{22}$. Luego, los valores propios de una matriz triangular son los elementos de la diagonal principal.
9. Sean $\lambda_1,\ldots,\lambda_n$ los valores propios de una matriz $A$ de orden $n$, que aprecen escritos tantas veces repetidos como indique su multiplicidad). Entonces:
   • $\displaystyle \sum_{i=1}^{n}\,\lambda_i=\text{traza}\,A\equiv \sum_{i=1}^{n}\,a_{ii}$
   • $\displaystyle \prod_{i=1}^{n}\,\lambda_i=\text{det}\,A$
10. El determinante de una matriz triangular es igual al producto de los elementos de la diagonal principal (consecuencia de propiedades anteriores).
11. Sea $A$ una matriz hermítica, esto es, una matriz tal que $A^\ast\equiv (\bar{A})^\top= A$. Entonces, todos sus valores propios son números reales.
12. Toda matriz hermítica, $A=A^\ast$, es semejante a una matriz diagonal $D$ mediante una transformación de semejanza cuya matriz de la transformación es una matriz unitaria $U$ (es decir $D=U^{\ast}AU$) —recordemos $U$ es unitaria si $U^{\ast}=U$—, entonces las columnas de $U$ forman una base de vectores propios ortonormales de $A$. Y en el caso de que el cuerpo sobre el que esté definida $A$ sea el de los números reales, $A$ es diagonalizable mediante una transformación de semejanza cuya matriz de la transformación es una matriz ortogonal $S=S^\top$, estando sus columnas asociadas a los respectivos vectores propios ortonormales de $A$.

Sobre las normas matriciales

Es sabido que podemos definir un espacio normado a partir de un espacio vectorial real o complejo. También podemos hablar de normas matriciales. Existe una definición genérica para la norma matricial de una matriz si bien en casos relativamente sencillos, y para matrices cuadradas, podemos trabajar con normas matriciales que se definen como herencia de las normas vectoriales y constituyen herramientas importante a la hora de analizar el comportamiento (convergencia, estabilidad,...) de algoritmos que hagan uso de matrices.

Definición genérica de una norma matricial de una matriz cuadrada

En general, siempre podemos definir una norma matricial $\left\|\,\right\|$ consistente con una norma vectorial dada (que notamos de la misma manera) mediante $$\left\| A \right\|:= \text{máx}\, \left\{ \dfrac{\left\|Ax\right\|}{\left\|x\right\|}: x \neq 0_V\right\}$$ Nota: en algunos libros (p.e. en [1]), se da esta otra condición: $$\left\| A \right\|:= \text{máx}\, \left\{ \left\| Ax \right\|: \left\| x \right\|=1 \right\}$$

Revisemos, a continuación, lo que sabemos sobre un espacio vectorial normado, para a continuación ir adentrándonos en la noción de normas matriciales.

Espacio vectorial normado

Consideremos un espacio vectorial $V$ sobre un cuerpo $\mathbb{K}$ (ya sea éste $\mathbb{R}$ o bien $\mathbb{C}^n$). Se define una norma en $V\equiv \mathbb{R}^n$ (o en $V\equiv \mathbb{C}^n$) como una aplicación $\left\|\,\right\|: V \rightarrow \mathbb{R}^{+}\cup\{0\};\,x\mapsto \left||x\right\|$, que cumple las siguiente propiedades:

1. Sea $x\in V$, entonces $\left||x\right\|=0_{\mathbb{R}}\Leftrightarrow x=0_V$
2. $\left||kx\right\|=|k|\,\left||x\right\|\;\forall\,\text{escalar}\,c\,\text{y}\,\forall\, x\in V$
3. $\left||x+y\right\|\le \left||x\right\|+\left||y\right\|\;\forall\, x,y\in V$

Nota: Recordemos que una norma en un espacio vectorial induce una métrica, que nos permite manejar distancias.

Unas normas vectoriales muy empleadas con las p-normas de Hölder: $\left\|x\right\|_{p}:=\left(\displaystyle \sum_{i=1}^{n}\,|x_i|^p\right)^{\frac{1}{p}}\;\text{donde}\;p\ge 1$, donde $|x|$ indica el módulo de un vector. Así, por ejemplo, se utilizan las siguientes normas particulares para todo conjunto de vectores $\{x_i\}\,,\,i=1\ldots,n$ de $V$ :

1. La norma suma de módulos se define haciendo $p=1$ en la norma genérica de Hölder: $\left\|x\right\|_1=\displaystyle \sum_{i=1}^{n}\,|x_i|$
2. La norma euclídea corresponde a poner $p=2$ en la norma genérica de Hölder: $\left\|x\right\|_1=\displaystyle \left(\sum_{i=1}^{n}\,|x_i|^2\right)^{\frac{1}{2}}$
3. La norma del máximo se define como el límite de la norma genérica de Hölder al hacer tender $p\rightarrow +\infty$, esto es $\left\|x\right\|_{\infty}:=\displaystyle \lim_{p\rightarrow +\infty}\,\left( \sum_{i=1}^{n}\,|x_i|^p\right)^{\frac{1}{p}}=\text{máx}_{i=1,\ldots\,n}\,\{|x_i|\}$

Normas matriciales subordinadas a alguna norma vectorial

Dicho esto, entedemos una norma matricial como una norma en el espacio vectorial de las matrices cuadradas $\mathcal{M}_{n,n}$ que cumpla el que sea multiplicativa, esto es, $\left\|AB\right\| \le \left\|A \right\|\,\left\|B \right\|\;\forall A,B\in \mathcal{M}_{n,n}$, y diremos que dicha norma matrical $\left\|\,\right\|$ es consistente con una cierta norma vectorial (está subordianda a una norma vectorial) si y sólo si $\left\|Ax\right\|\le \left\|A\right\|\,\left\|x\right\|\,\forall \,A\in \mathcal{M}_{n,n}\,\text{y}\,\forall x\in \mathbb{K}^n$

Ya se ha comentado más arriba que, en casos sencillos, pueden determinarse expresiones que proporcionen la norma matricial en función de los coeficientes de una matriz. Así, en particular, podemos hablar de las siguientes normas subordinadas a las normas vectoriales que hemos comentado arriba:

1. Norma matricial subordinada 1. Es la norma matricial subordinada a la norma vectorial suma de módulos (normal subordinada 1): $\left\|A\right\|_{1}:=\text{máx}_{1\le j \le n}\,\left\{\displaystyle \sum_{i=1}^{n}\,|a_{ij}|\right\}$ (esto es, el máximo por columnas de la suma de módulos de los elementos de todas y cada una de las columnas de la matriz)
2. Norma matricial subordinada 2. Es la norma subordinada a la norma vectorial euclídea: $\left\|A\right\|_{2}:=\rho(A^\top\,A)^{\frac{1}{2}}$ —que también puede escribirse como $\left\|A\right\|_{2}:=\sqrt{\rho(A^\top\,A)}$—, donde recordemos que $\rho$ denota el radio espectral (en este caso, de la matriz $A^\top\,A$); esto es, $\rho(A^\top\,A)=\text{máx}\,\{\text{valores propios en valor absoluto de}\,A^\top\,A\}$.
   
   En el caso de que $A$ sea simétrica ($A=A^\top$), esto se reduce a $\rho(A)$.
   
   Cabe también comentar el siguiente teorema: $\rho(A)=\text{ínfimo}\,\left\{\left\|A\right\|:\, \text{para cualquier norma matricial subordinada}\,\left\|.\right\|\right\}$. Como corolario de este teorema, se demuestra un resultado que ya había comentado más rarriba: (i) $\displaystyle \lim_{n\,\rightarrow\, \infty}\,A^n=0_{matriz}$, (ii) $\displaystyle \lim_{n\,\rightarrow\, \infty}\,\left\|A^n\right\|=0$, y (iii) $\rho(A)\lt 1$.
3. Norma matricial subordinada $\infty$. Es la norma matricial subordinada a la norma vectorial infinito: $\left\|A\right\|_{\infty}:=\text{máx}_{1\le i \le n}\,\left\{\displaystyle \sum_{j=1}^{n}\,|a_{ij}|\right\}$ (esto es, el máximo por filas de la suma de módulos de los elementos de todas y cada una de las filas de la matriz)

Otra norma matricial: la norma de Fröbenius

Además de las normas matriciales subordinadas que acabamos de comentar, también manejaremos esta otra norma específica, referida al producto escalar de matrices cuadradas, producto que se define como $A:B\,:=\text{traza}\,(A^\top\,B)$. Dicha norma es la que llamamos norma de Fröbenius y se define así: $$\left\|A\right\|_{F}:=(A:A)^{\frac{1}{2}}=\displaystyle \left( \sum_{i,j}^{n}\,a_{ij}^2 \right)^{\frac{1}{2}}$$ Téngase en cuenta que todas las normas matriciales toman el valor $1$ para la matriz identidad $A\equiv I$, pero la norma de Fröbenius toma el valor $\sqrt{n}$ para dicha matriz.

Acerca de la inestabilidad de los sitemas de eucaciones lineales $Ax=b$

En relación con la norma matricial (de la matriz de los coeficientes del sistema) es relevante también mencionar que existe una cantidad que denominamos número condición, $\text{Cond}(A)$ o $K(A)$ y que definimos de la forma $\text{Cond}(A):=\left\|A\right\|\cdot \left\|A^{-1}\right\|$. Se espera que la matriz se comporte bien si $\text{Cond}(A)\approx 1$, y un comportamiento defectuoso (por mal condicionamiento) si $\text{Cond}(A)$ es significativamente mayor que $1$.
$\diamond$

¿Y los sitemas de ecuaciones no lineales? ¿Cómo pueden resolverse numéricamente?

Solamente comentaré aquí, sin entrar en detalles, que existen métodos iterativos, como el método de Newton, que consisten en encontrar matrices apropiadas para construir una matriz iteradora. En el caso del método de Newton, dicha matriz se construye a partir de la matriz jacobiana (matriz de las derivadas parciales $\dfrac{\partial\,f_i}{\partial\,x_j}$ de las funciones $f_{i}(x_1,\ldots,x_n)=0$ ($i,j=1\ldots,n) de varias variables que corresponden a cada una de las ecuaciones del sistema (véase, p.e., [4] pp. 176-177, y un ejemplo de resolución de un sistema no lineal en las pp. 225-226).

-oOo-

Referencias:
[1] C. Moreno, Introducción al cálculo numérico (UNED, Madrid, 2011).
[2] A. Aubanell, A. Benseny y A. Deslshams, Útiles de cálculo numérico (Labor, Barcelona, 1993).
[3] B.P., Demidovich y I.A. Maron, Cálculo numérico fundamental (Paraninfo, Madrid, 1988).
[4] F. García y A. Nevot, Ejercicios resueltos de cálculo numérico (Paraninfo, Madrid, 1992).
[5] vv.aa., https://ca.wikipedia.org/wiki/Criteri_de_Sylvester (Vikipèdia).