Tema 7 - Resolución de sistemas lineales por métodos directos

Introducción

Sistema de ecuaciones lineal

Definición de sistema de ecuaciones lineal. Un sistema de ecuaciones lineal es un conjunto de ecuaciones del tipo: \[ \left. \begin{align*} E_1&: a_{11}x_1+a_{12}x_2+\cdots +a_{1n}x_n= b_1,\\ E_2& : a_{21}x_1+a_{22}x_2+\cdots +a_{2n}x_n= b_2,\\ & \vdots \\ E_n& : a_{n1}x_1+a_{n2}x_2+\cdots +a_{nn}x_n= b_n, \end{align*} \right\} \] donde \(E_i\) significa la ecuación \(i\)-ésima y recordemos que los coeficientes \((a_{ij})_{i,j=1,\ldots,n}\) son los coeficientes del sistema que se supone conocidos, \((b_1,\ldots,b_n)\) es el vector de términos independientes también conocido y \((x_1,x_2,\ldots,x_n)\) es el vector de incógnitas que tenemos que hallar.

Suponemos también que el sistema anterior tiene solución única, lo que es equivalente a \(\mathrm{det}((a_{ij})_{i,j=1,\ldots,n})\neq 0\).

Sistema de ecuaciones lineal

Un sistema de ecuaciones lineal se suele escribir en forma matricial de la forma siguiente: \[ \mathbf{A}\cdot\mathbf{x}=\mathbf{b}, \] donde \(\mathbf{A}\) es la denominada matriz del sistema, \(\mathbf{x}\) el vector de incógnitas y \(\mathbf{b}\), el vector de términos independientes:

\[ \mathbf{A}=\begin{bmatrix}a_{11}&a_{12}&\ldots&a_{1n} \\a_{21}&a_{22}&\ldots&a_{2n} \\\vdots&\vdots&\ddots&\vdots \\a_{n1}&a_{n2}&\ldots&a_{nn} \\\end{bmatrix},\quad \mathbf{x}=\begin{bmatrix}x_{1} \\x_2 \\\vdots \\x_n \\\end{bmatrix},\quad \mathbf{b}=\begin{bmatrix}b_{1} \\b_2 \\\vdots \\b_n \\\end{bmatrix}. \]

Sistema de ecuaciones lineal

En este curso supondremos que los sistemas de ecuaciones son compatibles determinados, es decir, tienen solución única.

El Teorema de Rouché Frobenius dice que un sistema de ecuaciones lineal es compatible si el rango de la matriz del sistema \(\mathbf{A}\) coincide con el rango de la matriz ampliada que se construye añadiendo como columna el vector de términos independientes a la matriz \(\mathbf{A}\): \(\mathbf{\overline{A}}=(\mathbf{A}|\mathbf{b})\).

Es decir, si \(\mathrm{rank}(\mathbf{A})=\mathrm{rank}(\mathbf{\overline{A}})\), el sistema es compatible. Si además \(\mathrm{rank}(\mathbf{A})=\mathrm{rank}(\mathbf{\overline{A}})=n\), el sistema es compatible determinado, es decir, tiene solución única.

Ejemplo

Consideremos el sistema de ecuaciones siguiente: \[ \left. \begin{align*} E_1 &: 6x_1+5x_2+3x_3= 1,\\ E_2 &:6x_1+4x_2+2x_3= -1,\\ E_3 &: 3x_1+6x_2+4x_3= 2.\\ \end{align*} \right\} \] o escrito en forma matricial: \[ \begin{bmatrix}6&5&3 \\6&4&2 \\3&6&4 \\\end{bmatrix}\cdot \begin{bmatrix}x_{1} \\x_2 \\x_3 \\\end{bmatrix}=\begin{bmatrix}1 \\-1 \\2 \\\end{bmatrix}. \] Puede comprobarse que el sistema anterior es compatible determinado y tiene como solución:

\[ x_1=0.3333333,\quad x_2=-3.5,\quad x_3=5.5. \]

Métodos directos

En este capítulo, vamos a ver métodos directos para resolver un sistema de ecuaciones de \(n\) ecuaciones y \(n\) variables o incógnitas:

Un método directo para resolver un sistema lineal se caracteriza porque lo resuelve en un número finito de pasos.

Para simplificar el sistema lineal podemos realizar las operaciones siguientes, es decir, que cualquier operación de las que sigue transforma el sistema lineal en un sistema lineal equivalente con la misma solución:

Simplificación del sistema

La ecuación \(E_i\) se puede multiplicar por cualquier constante \(k\neq 0\): \[ (kE_i)\rightarrow (E_i). \]
Se puede transformar la ecuación \(E_i\) a una nueva ecuación añadiendo cualquier ecuación \(E_j\) multiplicada por una constante \(k\): \[ (E_i+kE_j)\rightarrow (E_i). \]
Se pueden intercambiar cualquier par de ecuaciones \(E_i\) y \(E_j\): \[ (E_i)\leftrightarrow (E_j). \]

Ejemplo

Consideremos el sistema anterior: \[ \left. \begin{align*} E_1 &: 6x_1+5x_2+3x_3= 1,\\ E_2 &:6x_1+4x_2+2x_3= -1,\\ E_3 &: 3x_1+6x_2+4x_3= 2.\\ \end{align*} \right\} \] Realizamos las operaciones siguientes:

\((E_3)\leftrightarrow (E_1)\): \[ \left. \begin{align*} E_1 &: 3x_1+6x_2+4x_3= 2,\\ E_2 &:6x_1+4x_2+2x_3= -1,\\ E_3 &: 6x_1+5x_2+3x_3= 1.\\ \end{align*} \right\} \]
\((E_2-2E_1)\rightarrow (E_2)\): \[ \left. \begin{align*} E_1 &: 3x_1 & +6x_2+4x_3& = 2,\\ E_2 &: & -8x_2 -6x_3 & = -5,\\ E_3 &: 6x_1 & +5x_2+3x_3& = 1.\\ \end{align*} \right\} \]

Ejemplo (continuación)

\((E_3-2E_1)\rightarrow (E_3)\): \[ \left. \begin{align*} E_1 &: 3x_1 & +6x_2+4x_3& = 2,\\ E_2 &: & -8x_2 -6x_3 & = -5,\\ E_3 &: & -7x_2-5x_3& = -3.\\ \end{align*} \right\} \]
\((E_3-\frac{7}{8}E_2)\rightarrow (E_3)\): \[ \left. \begin{align*} E_1 &: 3x_1 & +6x_2+4x_3& = 2,\\ E_2 &: & -8x_2 -6x_3 & = -5,\\ E_3 &: & 0.25x_3& = 1.375.\\ \end{align*} \right\} \] Fijarse que hemos transformado el sistema de ecuaciones en un sistema de ecuaciones donde la matriz del sistema es triangular superior.

La resolución del sistema para este tipo de matrices es sencilla: basta empezar por la última e ir sustituyendo los valores obtenidos en las demás: \[ x_3=\frac{1.375}{0.25}=5.5,\quad x_2=\frac{-5+6\cdot 5.5}{-8}=-3.5,\quad x_1=\frac{2-4\cdot 5.5-6\cdot (-3.5)}{3}=0.3333333. \]

Ejemplo (continuación)

Escribamos matricialmente las operaciones realizadas:

\[ \begin{bmatrix}6&5&3&1 \\6&4&2&-1 \\3&6&4&2 \\\end{bmatrix}\rightarrow \begin{bmatrix}3&6&4&2 \\6&4&2&-1 \\6&5&3&1 \\\end{bmatrix}\rightarrow \begin{bmatrix}3&6&4&2 \\0&-8&-6&-5 \\0&-7&-5&-3 \\\end{bmatrix}\rightarrow \begin{bmatrix}3&6&4&2 \\0&-8&-6&-5 \\0&0&0.25&1.375 \\\end{bmatrix} \] Si no hubiésemos permutado las filas primera y tercera, las matrices resultantes de reducir la matriz del sistema a una matriz triangular superior serían: \[ \begin{bmatrix}6&5&3&1 \\6&4&2&-1 \\3&6&4&2 \\\end{bmatrix}\rightarrow \begin{bmatrix}6&5&3&1 \\0&-1&-1&-2 \\0&3.5&2.5&1.5 \\\end{bmatrix}\rightarrow \begin{bmatrix}6&5&3&1 \\0&-1&-1&-2 \\0&0&-1&-5.5 \\\end{bmatrix} \] Las soluciones del sistema ahora son fáciles de hallar: \[ x_3=\frac{-5.5}{-1}=5.5,\quad x_2=\frac{-2+5.5}{-1}=-3.5,\quad x_1=\frac{1-3\cdot 5.5-5\cdot (-3.5)}{6}=0.3333333. \] Acabamos de descubrir el método de eliminación de Gauss con sustitución hacia atrás.

Método de eliminación de Gauss

Consideramos un sistema de ecuaciones: \[ \left. \begin{align*} E_1&: a_{11}x_1+a_{12}x_2+\cdots +a_{1n}x_n= b_1,\\ E_2& : a_{21}x_1+a_{22}x_2+\cdots +a_{2n}x_n= b_2,\\ & \vdots \\ E_n& : a_{n1}x_1+a_{n2}x_2+\cdots +a_{nn}x_n= b_n, \end{align*} \right\} \] y denotamos por \(\overline{\mathbf{A}}\) la matriz ampliada del sistema:

\[ \overline{\mathbf{A}}=\begin{bmatrix}a_{11}&a_{12}&\ldots&a_{1n}&b_1 \\a_{12}&a_{22}&\ldots&a_{2n}&b_2 \\\vdots&\vdots&\ddots&\vdots&\vdots \\a_{n1}&a_{n2}&\ldots&a_{nn}&b_n \\\end{bmatrix}. \]

Método de eliminación de Gauss

El método de eliminación de Gauss consiste en realizar operaciones del tipo vistas anteriormente de cara a transformar el sistema lineal en un sistema lineal equivalente cuya matriz del sistema sea triangular superior.

Resolver un sistema de ecuaciones cuya matriz sea triangular superior es relativamente sencillo. La solución del mismo puede hallarse usando la técnica de sustitución hacia atrás:

Consideremos un sistema de ecuaciones con matriz del sistema triangular superior: \[ \left. \begin{align*} E_1&: a_{11}x_1& +a_{12}x_2+\cdots + &a_{1n}x_n & = b_1,\\ E_2& : & a_{22}x_2+\cdots +& a_{2n}x_n & = b_2,\\ & \vdots &\vdots & &\vdots \\ E_n& : && a_{nn}x_n & = b_n, \end{align*} \right\} \]

Sistema lineal triangular superior

La solución sería la siguiente: \[ \begin{align*} x_n & = \frac{b_n}{a_{nn}},\\ x_{n-1} &= \frac{b_{n-1}-a_{nn}x_n}{a_{n-1,n-1}},\\ \vdots & \\ x_i & = \frac{b_i-\displaystyle\sum_{j=i+1}^n a_{ij}x_j}{a_{ii}},\\ \vdots & \\ x_1 & = \frac{b_1-\displaystyle\sum_{j=2}^n a_{1j}x_j}{a_{11}}. \end{align*} \]

Sistema lineal triangular superior

Observación.

Todos los valores \(a_{ii}\neq 0\) son diferentes de \(0\), para \(i=1,2,\ldots,n\) ya que si algún \(a_{ii}=0\) fuese \(0\), el determinante del sistema sería \(0\) y el sistema no tendría solución única o sería incompatible y hemos dicho al principio del tema que suponemos que todos nuestros sistemas son compatibles con solución única.

En el ejemplo realizado anteriormente teníamos un sistema con matriz triangular superior.

Sistema lineal triangular superior

Proposición. El número de operaciones básicas (sumas, restas, multiplicaciones y divisiones) requerida para resolver un sistema lineal con matriz triangular superior vale: \(n^2\).

Demostración

En el paso \(i\)-ésimo para calcular el valor de \(x_i\) realizamos \(n-i\) sumas/restas y \(n-i+1\) multiplicaciones/divisiones.

El número de operaciones total será, pues, \[ \sum_{i=1}^n ((n-i)+(n-i+1))=\sum_{i=1}^n (2(n-i)+1)=n+2\sum_{i=1}^{n-1} i=n+2\cdot\frac{n\cdot (n-1)}{2}=n^2. \]

Método de eliminación de Gauss

Veamos cómo transformar un sistema lineal en otro equivalente con matriz del sistema triangular superior.

Primer paso. Supongamos que \(a_{11}\neq 0\). Realizamos las operaciones siguientes: \[ \left(E_j-\left(\frac{a_{j1}}{a_{11}}\right)E_1\right)\longrightarrow (E_j), \] para \(j=2,3,\ldots,n\).

Método de eliminación de Gauss

Es decir, cambiamos las filas \(2,3,\ldots,n\) haciendo que en la primera columna aparezcan ceros desde la componente \(2\) hasta la \(n\). La nueva matriz ampliada del sistema será:

\[ \overline{\mathbf{A}^{(2)}}=\begin{bmatrix}a_{11}&a_{12}&\ldots&a_{1n}&b_1 \\0&a_{22}^{(2)}&\ldots&a_{2n}^{(2)}&b_2^{(2)} \\\vdots&\vdots&\ddots&\vdots&\vdots \\0&a_{n2}^{(2)}&\ldots&a_{nn}^{(2)}&b_n^{(2)} \\\end{bmatrix}. \]

Método de eliminación de Gauss

Segundo paso. Tenemos que el sistema corresponde a la matriz ampliada \(\overline{\mathbf{A}^{(2)}}\). Supongamos que \(a_{22}^{(2)}\neq 0\). Realizamos las operaciones siguientes: \[ \left(E_j-\left(\frac{a_{j2}^{(2)}}{a_{22}^{(2)}}\right)E_2\right)\longrightarrow (E_j), \] para \(j=3,\ldots,n\).

Método de eliminación de Gauss

Es decir, cambiamos las filas \(3,\ldots,n\) haciendo que en la segunda columna aparezcan ceros desde la componente \(3\) hasta la \(n\). La nueva matriz ampliada del sistema será:

\[ \overline{\mathbf{A}^{(3)}}=\begin{bmatrix}a_{11}&a_{12}&a_{13}&\ldots&a_{1n}&b_1 \\0&a_{22}^{(2)}&a_{23}^{(2)}&\ldots&a_{2n}^{(2)}&b_2^{(2)} \\0&0&a_{33}^{(3)}&\ldots&a_{3n}^{(3)}&b_3^{(3)} \\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots \\0&0&a_{n3}^{(3)}&\ldots&a_{nn}^{(3)}&b_n^{(3)} \\\end{bmatrix}. \]

Método de eliminación de Gauss

Paso \(k\)-ésimo. Supongamos que el sistema ha sido transformado a un sistema con la matriz ampliada siguiente:

\[ \overline{\mathbf{A}^{(k)}}=\begin{bmatrix}a_{11}&\ldots&\ldots&\ldots&\ldots&a_{1n}&b_1 \\\vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\vdots \\0&\ldots&0&a_{kk}^{(k)}&\ldots&a_{kn}^{(k)}&b_{k}^{(k)} \\\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots \\0&\ldots&0&a_{nk}^{(k)}&\ldots&a_{nn}^{(k)}&b_n^{(k)} \\\end{bmatrix}. \]

Método de eliminación de Gauss

Realizamos las operaciones siguientes: \[ \left(E_j-\left(\frac{a_{jk}^{(k)}}{a_{kk}^{(k)}}\right)E_{k}\right)\longrightarrow (E_j), \] para \(j=k+1,\ldots,n\).

Método de eliminación de Gauss

Es decir, cambiamos las filas \(k+1,\ldots,n\) haciendo que en la columna \(k\)-ésima aparezcan ceros desde la componente \(k+1\) hasta la \(n\). La nueva matriz ampliada del sistema será:

\[ \overline{\mathbf{A}^{(k+1)}}=\begin{bmatrix}a_{11}&\ldots&\ldots&\ldots&\ldots&a_{1n}&b_1 \\\vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\vdots \\0&\ldots&0&a_{k+1,k+1}^{(k+1)}&\ldots&a_{k+1,n}^{(k+1)}&b_{k+1}^{(k+1)} \\\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots \\0&\ldots&0&a_{n,k+1}^{(k+1)}&\ldots&a_{nn}^{(k+1)}&b_n^{(k+1)} \\\end{bmatrix}. \]

Método de eliminación de Gauss

En el paso \(n-1\)-ésimo, la matriz del sistema que corresponderá a las \(n\) primeras columnas de la matriz ampliada del sistema \(\overline{\mathbf{A}^{(n)}}\) será triangular superior y el sistema se podrá resolver como hemos indicado anteriormente.

Los valores \(a_{kk}^{(k)}\) se denominan pivotes.

Escribamos cómo se calculan las componentes de la matriz \(\overline{\mathbf{A}^{(k+1)}}=(a_{ij}^{(k+1)})_{i,j=1,\ldots,n}\) en función de las componentes de la matriz \(\overline{\mathbf{A}^{(k)}}=(a_{ij}^{(k)})_{i,j=1,\ldots,n}\): \[ a_{ij}^{(k+1)}= \begin{cases} a_{ij}^{(k)}, & \mbox{si }i=1,\ldots,k,\ j=1,\ldots,n+1,\\ 0, &\mbox{si }i=k+1,\ldots,n,\ j=1,\ldots,k,\\ a_{ij}^{(k)}-\frac{a_{ik}^{(k)}}{a_{kk}^{(k)}}a_{kj}^{(k)}, &\mbox{si }i=k+1,\ldots,n,\ j=k+1,\ldots,n+1. \end{cases} \]

Ejemplo

Consideremos el siguiente sistema de ecuaciones con \(5\) ecuaciones y \(5\) incógnitas: \[ \begin{align*} E_1: & x_1 & + x_2 & -3 x_3 & - x_4 & -2 x_5 & = 2,\\ E_2: & x_1 & +2 x_2 & +2 x_3 & & +3 x_5 & = 2,\\ E_3: & x_1 & - x_2 & +3 x_3& +2 x_4& & = 2,\\ E_4: & & x_2 & & +4 x_4 & - x_5 & = 3,\\ E_5: & x_1 & +3 x_2& + x_3 & & +5 x_5 & = 1. \end{align*} \] Vamos a ver las matrices que resultan aplicando el método de Gauss:

Primer paso. Operaciones realizadas:
- \((E_2-E_1)\rightarrow (E_2)\): \[\begin{bmatrix}1&1&-3&-1&-2&2 \\1&2&2&0&3&2 \\1&-1&3&2&0&2 \\0&1&0&4&-1&3 \\1&3&1&0&5&1 \\\end{bmatrix}\longrightarrow \begin{bmatrix}1&1&-3&-1&-2&2 \\0&1&5&1&5&0 \\1&-1&3&2&0&2 \\0&1&0&4&-1&3 \\1&3&1&0&5&1 \\\end{bmatrix}\]

Ejemplo (continuación)

- \((E_3-E_1)\rightarrow (E_3)\): \[\begin{bmatrix}1&1&-3&-1&-2&2 \\0&1&5&1&5&0 \\1&-1&3&2&0&2 \\0&1&0&4&-1&3 \\1&3&1&0&5&1 \\\end{bmatrix}\longrightarrow \begin{bmatrix}1&1&-3&-1&-2&2 \\0&1&5&1&5&0 \\0&-2&6&3&2&0 \\0&1&0&4&-1&3 \\1&3&1&0&5&1 \\\end{bmatrix}\]
- \((E_5-E_1)\rightarrow (E_5)\): \[\begin{bmatrix}1&1&-3&-1&-2&2 \\0&1&5&1&5&0 \\0&-2&6&3&2&0 \\0&1&0&4&-1&3 \\1&3&1&0&5&1 \\\end{bmatrix}\longrightarrow \begin{bmatrix}1&1&-3&-1&-2&2 \\0&1&5&1&5&0 \\0&-2&6&3&2&0 \\0&1&0&4&-1&3 \\0&2&4&1&7&-1 \\\end{bmatrix}\]

Ejemplo (continuación)

Hemos calculado la matriz \(\overline{\mathbf{A}}^{(2)}\): \[ \overline{\mathbf{A}}^{(2)}=\begin{bmatrix}1&1&-3&-1&-2&2 \\0&1&5&1&5&0 \\0&-2&6&3&2&0 \\0&1&0&4&-1&3 \\0&2&4&1&7&-1 \\\end{bmatrix} \]
- \((E_3+2E_2)\rightarrow (E_3)\): \[\begin{bmatrix}1&1&-3&-1&-2&2 \\0&1&5&1&5&0 \\0&-2&6&3&2&0 \\0&1&0&4&-1&3 \\0&2&4&1&7&-1 \\\end{bmatrix}\longrightarrow \begin{bmatrix}1&1&-3&-1&-2&2 \\0&1&5&1&5&0 \\0&0&16&5&12&0 \\0&1&0&4&-1&3 \\0&2&4&1&7&-1 \\\end{bmatrix}\]

Ejemplo (continuación)

- \((E_4-E_2)\rightarrow (E_4)\): \[\begin{bmatrix}1&1&-3&-1&-2&2 \\0&1&5&1&5&0 \\0&0&16&5&12&0 \\0&1&0&4&-1&3 \\0&2&4&1&7&-1 \\\end{bmatrix}\longrightarrow \begin{bmatrix}1&1&-3&-1&-2&2 \\0&1&5&1&5&0 \\0&0&16&5&12&0 \\0&0&-5&3&-6&3 \\0&2&4&1&7&-1 \\\end{bmatrix}\]
- \((E_5-2E_2)\rightarrow (E_5)\): \[\begin{bmatrix}1&1&-3&-1&-2&2 \\0&1&5&1&5&0 \\0&0&16&5&12&0 \\0&0&-5&3&-6&3 \\0&2&4&1&7&-1 \\\end{bmatrix}\longrightarrow \begin{bmatrix}1&1&-3&-1&-2&2 \\0&1&5&1&5&0 \\0&0&16&5&12&0 \\0&0&-5&3&-6&3 \\0&0&-6&-1&-3&-1 \\\end{bmatrix}\]

Ejemplo (continuación)

Hemos calculado la matriz \(\overline{\mathbf{A}}^{(3)}\): \[ \overline{\mathbf{A}}^{(3)}=\begin{bmatrix}1&1&-3&-1&-2&2 \\0&1&5&1&5&0 \\0&0&16&5&12&0 \\0&0&-5&3&-6&3 \\0&0&-6&-1&-3&-1 \\\end{bmatrix}. \]
- \(\left(E_4+\frac{5}{16}E_3\right)\rightarrow (E_4)\): \[\begin{bmatrix}1&1&-3&-1&-2&2 \\0&1&5&1&5&0 \\0&0&16&5&12&0 \\0&0&-5&3&-6&3 \\0&0&-6&-1&-3&-1 \\\end{bmatrix}\longrightarrow \begin{bmatrix}1&1&-3&-1&-2&2 \\0&1&5&1&5&0 \\0&0&16&5&12&0 \\0&0&0&4.5625&-2.25&3 \\0&0&-6&-1&-3&-1 \\\end{bmatrix}\]

Ejemplo (continuación)

- \(\left(E_5+\frac{6}{16}E_3\right)\rightarrow (E_5)\): \[\begin{bmatrix}1&1&-3&-1&-2&2 \\0&1&5&1&5&0 \\0&0&16&5&12&0 \\0&0&0&4.5625&-2.25&3 \\0&0&-6&-1&-3&-1 \\\end{bmatrix}\longrightarrow \begin{bmatrix}1&1&-3&-1&-2&2 \\0&1&5&1&5&0 \\0&0&16&5&12&0 \\0&0&0&4.5625&-2.25&3 \\0&0&0&0.875&1.5&-1 \\\end{bmatrix}\]
Hemos calculado la matriz \(\overline{\mathbf{A}}^{(4)}\): \[ \overline{\mathbf{A}}^{(4)}=\begin{bmatrix}1&1&-3&-1&-2&2 \\0&1&5&1&5&0 \\0&0&16&5&12&0 \\0&0&0&4.5625&-2.25&3 \\0&0&0&0.875&1.5&-1 \\\end{bmatrix}. \]

Ejemplo (continuación)

- \(\left(E_5-\frac{0.875}{4.5625}E_4\right)\rightarrow (E_5)\): \[\begin{bmatrix}1&1&-3&-1&-2&2 \\0&1&5&1&5&0 \\0&0&16&5&12&0 \\0&0&0&4.5625&-2.25&3 \\0&0&0&0.875&1.5&-1 \\\end{bmatrix}\longrightarrow \begin{bmatrix}1&1&-3&-1&-2&2 \\0&1&5&1&5&0 \\0&0&16&5&12&0 \\0&0&0&4.5625&-2.25&3 \\0&0&0&0&1.931507&-1.575342 \\\end{bmatrix}\]
Hemos calculado la matriz \(\overline{\mathbf{A}}^{(5)}\): \[ \overline{\mathbf{A}}^{(5)}=\begin{bmatrix}1&1&-3&-1&-2&2 \\0&1&5&1&5&0 \\0&0&16&5&12&0 \\0&0&0&4.5625&-2.25&3 \\0&0&0&0&1.931507&-1.575342 \\\end{bmatrix}. \]

Ejemplo (continuación)

La solución del sistema será usando el método de sustitución hacia atrás: \[ \begin{align*} x_5 & =\frac{-1.575342}{1.931507}=-0.815603,\\ x_4 & =\frac{3-(-2.25)\cdot (-0.8156028)}{4.5625}=0.2553191,\\ x_3 & =\frac{0-12\cdot (-0.8156028)-5\cdot 0.2553191}{16}=0.5319149,\\ x_2 & =\frac{0-5\cdot (-0.8156028)-1\cdot 0.2553191-5\cdot 0.5319149}{1}=1.1631206,\\ x_1 & =\frac{2-(-2)\cdot (-0.8156028)-(-1)\cdot 0.2553191-(-3)\cdot 0.5319149-1\cdot 1.1631206}{1}=1.0567376. \end{align*} \]

Método de Gauss. Pseudocódigo

INPUT número de incógnitas y ecuaciones \(n\), matriz ampliada \(\overline{\mathbf{A}}=(a_{ij})_{i=1,\ldots,n,j=1,\ldots,n+1}\).
For k=1,...,n-1 (proceso de eliminación)
- Sea \(p\) el menor entero tal que \(a_{pk}\neq 0\) (buscamos el primer elemento de la columna \(k\)-ésima que sea distinto de \(0\), si dicho elemento no existe, el sistema no tiene solución.)
- If \(p\) no existe then el sistema no tiene solución.
- If \(p\neq k\) then realizar la operación \((E_p)\leftrightarrow (E_k)\) (cambiar las filas \(p\) y \(k\))

Método de Gauss. Pseudocódigo

- For \(j=k+1,\ldots,n\) do
  - Set \(m_{jk}=\frac{a_{jk}}{a_{kk}}\) (calculamos el valor por el que debemos multiplicar la fila \(j\)-ésima);
  - Realizar la operación\((E_j-m_{jk}E_k)\rightarrow (E_j)\). (cambiamos la fila \(j\)-ésima por la “antigua” fila \(j\)-ésima menos la fila \(k\)-ésima multiplicada por \(m_{jk}\))
Set \(x_n=\frac{a_{n,n+1}}{a_{nn}}\) (Resolvemos el sistema triangular del cuya matriz ampliada es \(\overline{\mathbf{A}^{(n)}}\)).
For \(i=n-1,\ldots,1\), set \(x_i=\frac{a_{i,n+1}-\sum_{j=i+1}^n a_{ij}x_j}{a_{ii}}\).
Dar la solución \(x_1,\ldots,x_n\).

Implementación del método de Gauss

El método de Gauss para resolver el sistema: \[ \begin{align*} E_1: & x_1 & + x_2 & -3 x_3 & - x_4 & -2 x_5 & = 2,\\ E_2: & x_1 & +2 x_2 & +2 x_3 & & +3 x_5 & = 2,\\ E_3: & x_1 & - x_2 & +3 x_3& +2 x_4& & = 2,\\ E_4: & & x_2 & & +4 x_4 & - x_5 & = 3,\\ E_5: & x_1 & +3 x_2& + x_3 & & +5 x_5 & = 1. \end{align*} \] se puede encontrar implementado en:

Método de Gauss. Número de operaciones

Vamos a contar el número de operaciones básicas (sumas/restas y multiplicaciones/divisiones) que requiere el método de Gauss para resolver un sistema de ecuaciones \(n\times n\).

Primero contemos el número de operaciones para pasar la matriz del sistema a una matriz triangular superior.

Recordemos que el método de Gauss consta de \(n-1\) pasos. Veamos cuántas operaciones realizamos en el paso \(k\)-ésimo, es decir, para pasar de la matriz \(\overline{\mathbf{A}^{(k)}}\) a la matriz \(\overline{\mathbf{A}^{(k+1)}}\), \(k=1,\ldots,n-1\):

Método de Gauss. Número de operaciones

Número de multiplicaciones/divisiones:
- cálculo de los \(m_{ik}=\frac{a_{ik}}{a_{kk}}\) para \(i=k+1,\ldots,n,\ \Rightarrow n-k\) operaciones,
- cálculo de los nuevos \(a_{ij}^{(k)}\) para \(i=k+1,\ldots,n,\ j=k+1,\ldots,n+1,\ \Rightarrow (n-k)\cdot (n-k+1)\) operaciones,
- número total de operaciones: \[ \begin{align*} & n-k+(n-k)\cdot (n-k+1)\\ & =(n-k)\cdot (1+n-k+1) =(n-k)\cdot (n-k+2). \end{align*} \]

Método de Gauss. Número de operaciones

Número de sumas/restas
- cálculo de los nuevos \(a_{ij}^{(k)}\) para \(i=k+1,\ldots,n,\ j=k+1,\ldots,n+1,\ \Rightarrow (n-k)\cdot (n-k+1)\),
- número total de operaciones: \[ (n-k)\cdot (n-k+1). \]

Método de Gauss. Número de operaciones

En total, el número de operaciones básicas para pasar la matriz del sistema a una matriz triangular superior será: \[ \begin{align*} & \sum_{k=1}^{n-1} \left((n-k)\cdot (n-k+2)+(n-k)\cdot (n-k+1)\right)\\ & =\sum_{k=1}^{n-1} (n-k)(2(n-k)+3) = \sum_{k=1}^{n-1} k\cdot (2k+3)\\ & =2\sum_{k=1}^{n-1}k^2+3\sum_{k=1}^{n-1}k = \frac{1}{3} (n-1) n (2 n-1)+\frac{3}{2} (n-1) n \\ & = \frac{1}{6} (n-1) n (4 n+7) = \frac{2n^3}{3}+\frac{n^2}{2}-\frac{7n}{6}. \end{align*} \]

Método de Gauss. Número de operaciones

A continuación, tenemos que tener en cuenta el número de operaciones requerida para resolver un sistema triangular superior.

Recordemos que dicho valor ya está calculado y valía \(n^2\).

En resumen, el número total de operaciones básicas para resolver un sistema de ecuaciones \(n\times n\) usando el método de Gauss sería: \[ \frac{2n^3}{3}+\frac{n^2}{2}-\frac{7n}{6}+n^2 = \frac{2n^3}{3}+\frac{3n^2}{2}-\frac{7n}{6}. \]

Observación. El orden computacional del método de Gauss vale \(O(n^3)\).

Técnicas de pivotaje

Recordemos que en el paso \(k\)-ésimo del algoritmo de Gauss, de cara a hacer ceros en la columna \(k\)-ésima desde la componente \(j=k+1\) hasta la componente \(j=n\), hay que multiplicar la fila \(k\)-ésima por \(\frac{a_{jk}^{(k)}}{a_{kk}^{(k)}}\).

Los errores de redondeo se amplifican si el valor \(a_{kk}^{(k)}\) es pequeño. Por tanto, podemos modificar el algoritmo de Gauss eligiendo el valor de la columna \(k\)-ésima que maximiza \(|a_{jk}^{(k)}|\) desde \(j=k\) hasta \(j=n\). Llamemos \(j_{\mathrm{max}}\) al índice correspondiente a dicho valor.

A continuación, realizamos un cambio entre las filas \(k\) y \(j_{\mathrm{max}}\): \(E_{k}\leftrightarrow E_{j_{\mathrm{max}}}.\)

La modificación del algoritmo de Gauss usando la técnica anterior se denomina método de Gauss con pivotaje parcial.

Ilustremos la necesidad de realizar el pivotaje parcial con el siguiente ejemplo:

Ejemplo

Consideremos el sistema de ecuaciones siguiente: \[ \left. \begin{align*} 0.010534x_1+0.02x_2= & -1,\\ x_1+2x_2= & 3. \end{align*} \right\} \] Suponemos que trabajamos con \(5\) dígitos significativos.

Resolver el sistema anterior es equivalente desde el punto de vista geométrico a hallar la intersección de las rectas: \(r_1: 0.010534x_1+0.02x_2= -1\) y \(r_2: x_1+2x_2= 3\).

Ahora bien, como el determinante de la matriz del sistema anterior tiene un valor pequeño (\(0.001068\)), las dos rectas anteriores serían “casi paralelas”. Esto implica que se trata de un problema inestable desde el punto de vista numérico.

La solución del sistema anterior usando en principio los dígitos necesarios sería: \[ x_1=-1928.8389513,\quad x_2=965.9194757. \]

Ejemplo (continuación)

Si resolvemos el sistema sin realizar pivotaje parcial y teniendo en cuenta que trabajamos con \(5\) cifras significativas, obtenemos la matriz triangular siguiente: \[ \begin{bmatrix}0.010534&0.02&-1 \\0&0.1014&97.931 \\\end{bmatrix}, \] cuya solución es: \[ x_2=\frac{97.931}{0.1014}=965.79,\quad x_1=\frac{-1-0.02\cdot 965.79}{0.010534}=\frac{-20.316}{0.010534}=-1928.6. \]

Ejemplo (continuación)

Si, en cambio, resolvemos el sistema realizando pivotaje parcial, hay que permutar inicialmente las filas \(1\) y \(2\) obteniendo la matriz triangular siguiente después de aplicar el algoritmo de Gauss a la matriz resultante de la permutación: \[ \begin{bmatrix}1&2&3 \\0&-0.001068&-1.0316 \\\end{bmatrix}, \] cuya solución es: \[ x_2=\frac{-1.0316}{-0.001068}=965.92,\quad x_1=\frac{3-2\cdot 965.92}{1}=\frac{-1928.8}{1}=-1928.8. \] Comprobamos que se comete menor error realizando pivotaje parcial.

Pseudocódigo. Pivotaje parcial

INPUT número de incógnitas y ecuaciones \(n\), matriz ampliada \(\overline{\mathbf{A}}=(a_{ij})_{i=1,\ldots,n,j=1,\ldots,n+1}\).
For k=1,...,n-1 (proceso de eliminación)
- Sea \(p\) el menor entero tal que \(|a_{pk}|=\max_{j=k,\ldots,n}|a_{jk}|\).
- If \(a_{pk}=0\), then el sistema no tiene solución.
- If \(p\neq k\) then realizar la operación \((E_p)\leftrightarrow (E_k)\) (cambiar las filas \(p\) y \(k\))

Pseudocódigo. Pivotaje parcial

- For \(j=k+1,\ldots,n\) do
  - Set \(m_{jk}=\frac{a_{jk}}{a_{kk}}\) (calculamos el valor por el que debemos multiplicar la fila \(j\)-ésima);
  - Realizar la operación\((E_j-m_{jk}E_k)\rightarrow (E_j)\). (cambiamos la fila \(j\)-ésima por la “antigua” fila \(j\)-ésima menos la fila \(k\)-ésima multiplicada por \(m_{jk}\))
Set \(x_n=\frac{a_{n,n+1}}{a_{nn}}\) (Resolvemos el sistema triangular del cuya matriz ampliada es \(\overline{\mathbf{A}^{(n)}}\)).
For \(i=n-1,\ldots,1\), set \(x_i=\frac{a_{i,n+1}-\sum_{j=i+1}^n a_{ij}x_j}{a_{ii}}\).
Dar la solución \(x_1,\ldots,x_n\).

Ejemplo

Resolvamos el sistema propuesto anteriormente usando pivotaje parcial: \[ \begin{align*} E_1: & x_1 & + x_2 & -3 x_3 & - x_4 & -2 x_5 & = 2,\\ E_2: & x_1 & +2 x_2 & +2 x_3 & & +3 x_5 & = 2,\\ E_3: & x_1 & - x_2 & +3 x_3& +2 x_4& & = 2,\\ E_4: & & x_2 & & +4 x_4 & - x_5 & = 3,\\ E_5: & x_1 & +3 x_2& + x_3 & & +5 x_5 & = 1. \end{align*} \]

Primer paso (no se ha realizado ningún cambio de filas): \[ A^{(2)}=\begin{bmatrix}1&1&-3&-1&-2&2 \\0&1&5&1&5&0 \\0&-2&6&3&2&0 \\0&1&0&4&-1&3 \\0&2&4&1&7&-1 \\\end{bmatrix}. \]

Ejemplo (continuación)

Segundo paso (hemos cambiado las filas \(2\) y \(3\)): \[ A^{(3)}=\begin{bmatrix}1&1&-3&-1&-2&2 \\0&-2&6&3&2&0 \\0&0&8&2.5&6&0 \\0&0&3&5.5&0&3 \\0&0&10&4&9&-1 \\\end{bmatrix}. \]
Tercer paso (hemos cambiado las filas \(3\) y \(5\)): \[ A^{(4)}=\begin{bmatrix}1&1&-3&-1&-2&2 \\0&-2&6&3&2&0 \\0&0&10&4&9&-1 \\0&0&0&4.3&-2.7&3.3 \\0&0&0&-0.7&-1.2&0.8 \\\end{bmatrix}. \]

Ejemplo (continuación)

Cuarto paso (no ha habido cambio de filas): \[ A^{(5)}=\begin{bmatrix}1&1&-3&-1&-2&2 \\0&-2&6&3&2&0 \\0&0&10&4&9&-1 \\0&0&0&4.3&-2.7&3.3 \\0&0&0&0&-1.639535&1.337209 \\\end{bmatrix}. \]

La solución del sistema será usando el método de sustitución hacia atrás: \[ \begin{align*} x_5 & =\frac{1.337209}{-1.639535}=-0.815603,\\ x_4 & =\frac{3.3-(-2.7)\cdot (-0.8156028)}{4.3}=0.2553191,\\ x_3 & =\frac{-1-9\cdot (-0.8156028)-4\cdot 0.2553191}{10}=0.5319149,\\ x_2 & =\frac{0-2\cdot (-0.8156028)-3\cdot 0.2553191-6\cdot 0.5319149}{-2}=1.1631206,\\ x_1 & =\frac{2-(-2)\cdot (-0.8156028)-(-1)\cdot 0.2553191-(-3)\cdot 0.5319149-1\cdot 1.1631206}{1}=1.0567376. \end{align*} \]

Ejemplo (continuación)

El ejemplo anterior se encuentra resuelto en:

Pivotaje parcial escalado

Vamos a modificar ligeramente el pivotaje parcial de la forma que explicamos a continuación.

Dicha modificación es una mejora del pivotaje parcial en el sentido que se escalan las filas de la matriz del sistema a resolver por su máximo valor.

Concretamente, en primer lugar, calculamos el valor máximo en valor absoluto de cada fila \(i\)-ésima de la matriz del sistema a resolver \(\mathbf{A}\), \(n\times n\): \[ s_i=\max_{j=1,\ldots,n}|a_{ij}|. \] Dichos valores \(s_i\) sólo se calculan una vez al principio del algoritmo y no aumentan demasiado el coste computacional del mismo.

Pivotaje parcial escalado

A continuación, imaginemos que estamos en el paso \(k\)-ésimo del algoritmo de Gauss. Entonces, en lugar de hallar el máximo de \(|a_{jk}^{(k)}|\) desde \(j=k\) hasta \(j=n\), hallamos el máximo de \(\frac{|a_{jk}^{(k)}|}{s_j}\) desde \(j=k\) hasta \(j=n\): \[ \max_{j=k,\ldots,n}\frac{|a_{jk}^{(k)}|}{s_j} \]

Es decir, “escalamos” cada fila por el valor \(s_j\) y hallamos el máximo de los valores escalados.

Llamemos \(j_{\mathrm{max}}\) a dicho valor.

Igual que en el pivotaje parcial, realizamos un cambio entre las filas \(k\) y \(j_{\mathrm{max}}\): \(E_{k}\leftrightarrow E_{j_{\mathrm{max}}}.\)

Ejemplo

Resolvamos el sistema propuesto anteriormente usando pivotaje parcial escalado: \[ \begin{align*} E_1: & x_1 & + x_2 & -3 x_3 & - x_4 & -2 x_5 & = 2,\\ E_2: & x_1 & +2 x_2 & +2 x_3 & & +3 x_5 & = 2,\\ E_3: & x_1 & - x_2 & +3 x_3& +2 x_4& & = 2,\\ E_4: & & x_2 & & +4 x_4 & - x_5 & = 3,\\ E_5: & x_1 & +3 x_2& + x_3 & & +5 x_5 & = 1. \end{align*} \]

Los valores \(s_i\) son en este caso:

\[ 3, 3, 3, 4, 5. \]

Primer paso (no se ha realizado ningún cambio de filas ya que el máximo valor entre los valores \(\left|\frac{1}{3}\right|,\ \left|\frac{1}{3}\right|,\ \left|\frac{1}{3}\right|,\ \left|\frac{0}{4}\right|,\ \left|\frac{1}{5}\right|\) vale \(\left|\frac{1}{3}\right|\)): \[ A^{(2)}=\begin{bmatrix}1&1&-3&-1&-2&2 \\0&1&5&1&5&0 \\0&-2&6&3&2&0 \\0&1&0&4&-1&3 \\0&2&4&1&7&-1 \\\end{bmatrix}. \]

Ejemplo (continuación)

Segundo paso (hemos cambiado las filas \(2\) y \(3\) ya que el máximo valor entre los valores \(\left|\frac{1}{3}\right|,\ \left|\frac{-2}{3}\right|,\ \left|\frac{1}{4}\right|,\ \left|\frac{2}{5}\right|\) vale \(\left|\frac{-2}{3}\right|\)): \[ A^{(3)}=\begin{bmatrix}1&1&-3&-1&-2&2 \\0&-2&6&3&2&0 \\0&0&8&2.5&6&0 \\0&0&3&5.5&0&3 \\0&0&10&4&9&-1 \\\end{bmatrix}. \]
Tercer paso (no hemos realizado cambio de filas ya que el máximo valor entre los valores \(\left|\frac{8}{3}\right|,\ \left|\frac{3}{4}\right|,\ \left|\frac{10}{5}\right|\) vale \(\left|\frac{8}{3}\right|\)): \[ A^{(4)}=\begin{bmatrix}1&1&-3&-1&-2&2 \\0&-2&6&3&2&0 \\0&0&8&2.5&6&0 \\0&0&0&4.5625&-2.25&3 \\0&0&0&0.875&1.5&-1 \\\end{bmatrix}. \]

Ejemplo (continuación)

Cuarto paso (no ha habido cambio de filas ya que el máximo valor entre los valores \(\left|\frac{4.5625}{4}\right|,\ \left|\frac{0.875}{5}\right|\) vale \(\left|\frac{4.5625}{4}\right|\)): \[ A^{(5)}=\begin{bmatrix}1&1&-3&-1&-2&2 \\0&-2&6&3&2&0 \\0&0&8&2.5&6&0 \\0&0&0&4.5625&-2.25&3 \\0&0&0&0&1.931507&-1.575342 \\\end{bmatrix}. \]

Ejemplo (continuación)

La solución del sistema será usando el método de sustitución hacia atrás: \[ \begin{align*} x_5 & =\frac{-1.575342}{1.931507}=-0.815603,\\ x_4 & =\frac{3-(-2.25)\cdot (-0.8156028)}{4.5625}=0.2553191,\\ x_3 & =\frac{0-6\cdot (-0.8156028)-2.5\cdot 0.2553191}{8}=0.5319149,\\ x_2 & =\frac{0-2\cdot (-0.8156028)-3\cdot 0.2553191-6\cdot 0.5319149}{-2}=1.1631206,\\ x_1 & =\frac{2-(-2)\cdot (-0.8156028)-(-1)\cdot 0.2553191-(-3)\cdot 0.5319149-1\cdot 1.1631206}{1}=1.0567376. \end{align*} \]

Ejercicio

Programar la variante del método de Gauss con pivotaje parcial escalado que hemos visto en el ejemplo anterior.

Pivotaje global o maximal

Una manera de optimizar el error de redondeo cometido cuando aplicamos el algoritmo de Gauss es realizar un pivotaje global:

es decir, modificar el algoritmo de Gauss en el paso \(k\)-ésimo eligiendo el máximo valor en valor absoluto de la submatriz de \(\mathbf{A}^{(k)}\) formada por las filas \(k, k+1,\ldots,n\) y por las columnas \(k,k+1,\ldots,n\).

Concretamente hallamos \(i_{\mathrm{max}},j_{\mathrm{max}}\) tal que: \[ |a_{i_{\mathrm{max}},j_{\mathrm{max}}}^{(k)}|=\max_{i,j=k,\ldots,n}|a_{ij}^{(k)}|. \]

Pivotaje global o maximal

A continuación, realizamos un cambio entre las filas \(k\) y \(i_{\mathrm{max}}\) y un cambio entre las columnas \(k\) y \(j_{\mathrm{max}}\).

Permutar filas no afecta a las soluciones del sistema de ecuaciones al obtener un sistema de ecuaciones equivalente. Sin embargo, permutar columnas equivale a permutar los valores de las variables correspondientes.

Concretamente, cuando cambiamos las columnas \(k\) y \(j_{\mathrm{max}}\), intercambiamos el “papel” de las variables \(x_k\) y \(x_{j_{\mathrm{max}}}\).

Por tanto, necesitamos recordar todos los cambios de columnas para realizar el cambio inverso en la solución final.

Ejemplo

Apliquemos pivotaje global al sistema que hemos ido desarrollando:

\[ \begin{align*} E_1: & x_1 & + x_2 & -3 x_3 & - x_4 & -2 x_5 & = 2,\\ E_2: & x_1 & +2 x_2 & +2 x_3 & & +3 x_5 & = 2,\\ E_3: & x_1 & - x_2 & +3 x_3& +2 x_4& & = 2,\\ E_4: & & x_2 & & +4 x_4 & - x_5 & = 3,\\ E_5: & x_1 & +3 x_2& + x_3 & & +5 x_5 & = 1. \end{align*} \]

Paso 1. El máximo valor de la submatriz del sistema formada por las 5 primeras filas y columnas vale \(5\) que corresponde a la fila \(5\) y a la columna \(5\). Entonces permutamos la fila \(1\) con la \(5\) y a continuación la columna \(1\) con la \(5\) y aplicamos el algoritmo de Gauss obteniendo la matriz siguiente: \[ \begin{bmatrix}5&3&1&0&1&1 \\0&0.2&1.4&0&0.4&1.4 \\0&-1&3&2&1&2 \\0&1.6&0.2&4&0.2&3.2 \\0&2.2&-2.6&-1&1.4&2.4 \\\end{bmatrix}. \] Las variables del sistema anterior son las siguientes \(x_5,x_2,x_3,x_4\) y \(x_1\).

Ejemplo (continuación)

Paso 2. El máximo valor de la submatriz del sistema formada por las 4 últimas filas y columnas vale \(4\) que corresponde a la fila \(4\) y a la columna \(4\). Entonces permutamos la fila \(2\) con la \(4\) y a continuación la columna \(2\) con la \(4\) y aplicamos el algoritmo de Gauss obteniendo la matriz siguiente: \[ \begin{bmatrix}5&0&1&3&1&1 \\0&4&0.2&1.6&0.2&3.2 \\0&0&2.9&-1.8&0.9&0.4 \\0&0&1.4&0.2&0.4&1.4 \\0&0&-2.55&2.6&1.45&3.2 \\\end{bmatrix}. \] Las variables del sistema anterior son las siguientes \(x_5,x_4,x_3,x_2\) y \(x_1\).

Ejemplo (continuación)

Paso 3. El máximo valor de la submatriz del sistema formada por las 3 últimas filas y columnas vale \(2.9\) que corresponde a la fila \(3\) y a la columna \(3\). Entonces no realizamos ninguna permutación y aplicamos el algoritmo de Gauss obteniendo la matriz siguiente: \[ \begin{bmatrix}5&0&1&3&1&1 \\0&4&0.2&1.6&0.2&3.2 \\0&0&2.9&-1.8&0.9&0.4 \\0&0&0&1.068966&-0.034483&1.206897 \\0&0&0&1.017241&2.241379&3.551724 \\\end{bmatrix}. \] Las variables del sistema anterior son las siguientes \(x_5,x_4,x_3,x_2\) y \(x_1\).

Ejemplo (continuación)

Paso 4. El máximo valor de la submatriz del sistema formada por las 2 últimas filas y columnas vale \(2.2413793\) que corresponde a la fila \(5\) y a la columna \(5\). Entonces permutamos la fila \(4\) con la \(5\) y a continuación la columna \(4\) con la \(5\) y aplicamos el algoritmo de Gauss obteniendo la matriz siguiente: \[ \begin{bmatrix}5&0&1&1&3&1 \\0&4&0.2&0.2&1.6&3.2 \\0&0&2.9&0.9&-1.8&0.4 \\0&0&0&2.241379&1.017241&3.551724 \\0&0&0&0&1.084615&1.261538 \\\end{bmatrix}. \] Las variables del sistema anterior son las siguientes \(x_5,x_4,x_3,x_1\) y \(x_2\).

Ejemplo (continuación)

La solución del sistema será usando el método de sustitución hacia atrás: \[ \begin{align*} x_2 & =\frac{1.261538}{1.084615}=1.163121,\\ x_1 & =\frac{3.551724-1.0172414\cdot 1.1631206}{2.2413793}=1.0567376,\\ x_3 & =\frac{0.4-(-1.8)\cdot 1.1631206-0.9\cdot 1.0567376}{2.9}=0.5319149,\\ x_4 & =\frac{3.2-1.6\cdot 1.1631206-0.2\cdot 1.0567376-0.2\cdot 0.5319149}{4}=0.2553191,\\ x_5 & =\frac{1-3\cdot 1.1631206-1\cdot 1.0567376-1\cdot 0.5319149-0\cdot 0.2553191}{5}=-0.8156028. \end{align*} \]

Factorización de matrices

Introducción

Vamos a ver otro método de resolver un sistema de ecuaciones \(n\times n\) de la forma \(\mathbf{A}\mathbf{x}=\mathbf{b}\).

Dicho método se llama descomposición \(LU\). La idea es factorizar la matriz del sistema \(\mathbf{A}\) de la forma siguiente: \(\mathbf{A}=\mathbf{L}\mathbf{U}\), donde la matriz \(\mathbf{L}\) es triangular inferior con unos en la diagonal y \(\mathbf{U}\) es triangular superior.

Si somos capaces de realizar la descomposición anterior, resolver el sistema lineal \(\mathbf{A}\mathbf{x}=\mathbf{b}\) equivale a resolver dos sistemas triangulares. Veamos cómo.

El sistema \(\mathbf{A}\mathbf{x}=\mathbf{b}\) se puede escribir como \(\mathbf{L}\mathbf{U}\mathbf{x}=\mathbf{b}\). Para resolver el sistema anterior, hacemos lo siguiente:

En primer lugar, resolvemos el sistema \(\mathbf{L}\mathbf{y}=\mathbf{b}\), que es un sistema triangular inferior.
En segundo lugar, resolvemos el sistema \(\mathbf{U}\mathbf{x}=\mathbf{y}\), donde \(\mathbf{y}\) es la solución hallada en el primer paso.

Solución de un sistema triangular inferior

La resolución del sistema \(\mathbf{L}\mathbf{y}=\mathbf{b}\) se puede resolver usando la técnica de sustitución hacia adelante.

La matriz \(\mathbf{L}\) será de la forma:

\[ \begin{bmatrix}1&0&0&\ldots&0 \\l_{21}&1&0&\ldots&0 \\\vdots&\vdots&\ddots&\vdots&\vdots \\l_{n1}&l_{n2}&l_{n3}&\ldots&1 \\\end{bmatrix}. \]

Solución de un sistema triangular inferior

El sistema \(\mathbf{L}\mathbf{y}=\mathbf{b}\) será: \[ \left. \begin{align*} E_1&: y_1 & = b_1,\\ E_2& : l_{21}y_1+y_2 & = b_2,\\ & \vdots &\vdots \\ E_n& : l_{n1} y_1+l_{n2}y+\cdots+ y_n & = b_n. \end{align*} \right\} \]

Solución de un sistema triangular inferior

La solución del sistema anterior se puede obtener por sustitución hacia adelante: \[ \begin{align*} y_1 & = b_1,\\ y_2 &= b_{2}-l_{21}y_1,\\ \vdots & \\ y_i & = b_i-\sum_{j=1}^{i-1} l_{ij}y_j,\\ \vdots & \\ y_n & = b_n-\sum_{j=1}^{n-1} l_{nj}y_j \end{align*} \]

Solución del sistema triangular superior

El número de operaciones básicas para resolver el sistema \(\mathbf{L}\mathbf{y}=\mathbf{b}\) será:

Multiplicaciones/divisiones: \(\displaystyle\sum_{i=1}^{n}(i-1)=\sum_{i=1}^{n-1}i=\frac{n(n-1)}{2}\).
Sumas/restas: \(\displaystyle\sum_{i=1}^{n}(i-1)=\frac{n(n-1)}{2}\).
Total: \(2\frac{n(n-1)}{2}=n^2-n.\)

La solución del sistema \(\mathbf{U}{x}=\mathbf{y}\) se hace usando la técnica de sustitución hacia atrás explicada anteriormente donde recordemos que se requerían \(n^2\) operaciones básicas.

En conclusión, una vez se tiene la descomposición \(LU\), resolver el sistema comporta realizar \(n^2-n+n^2=2n^2-n\) operaciones básicas. Es decir, es de orden \(O(n^2)\).

Ejemplo

Expliquemos cómo resolver el sistema considerado anteriormente:

La matriz del sistema es la siguiente: \[ \mathbf{A}=\begin{bmatrix}1&1&-3&-1&-2 \\1&2&2&0&3 \\1&-1&3&2&0 \\0&1&0&4&-1 \\1&3&1&0&5 \\\end{bmatrix}. \]

Ejemplo (continuación)

La descomposición \(LU\) de la matriz anterior es la siguiente:

\[ \begin{bmatrix}1&1&-3&-1&-2 \\1&2&2&0&3 \\1&-1&3&2&0 \\0&1&0&4&-1 \\1&3&1&0&5 \\\end{bmatrix}=\begin{bmatrix}1&0&0&0&0 \\1&1&0&0&0 \\1&-2&1&0&0 \\0&1&-0.3125&1&0 \\1&2&-0.375&0.191781&1 \\\end{bmatrix}\cdot \begin{bmatrix}1&1&-3&-1&-2 \\0&1&5&1&5 \\0&0&16&5&12 \\0&0&0&4.5625&-2.25 \\0&0&0&0&1.931507 \\\end{bmatrix}. \] En primer lugar resolvemos \(\mathbf{L}\mathbf{y}=\mathbf{b}\): \[ \begin{bmatrix}1&0&0&0&0 \\1&1&0&0&0 \\1&-2&1&0&0 \\0&1&-0.3125&1&0 \\1&2&-0.375&0.191781&1 \\\end{bmatrix}\cdot\begin{bmatrix}y_1\\ y_2\\ y_3\\ y_4\\ y_5\end{bmatrix}=\begin{bmatrix}2 \\2 \\2 \\3 \\1 \\\end{bmatrix}. \]

\[ \begin{align*} y_1 & =2,\ y_2=2- y_1=0,\ y_3=2+2 y_2- y_1=0,\ y_4=3+0.3125 y_3- y_2=3,\\ y_5 & =1-0.191781 y_4+0.375 y_3-2 y_2- y_1=-1.5753425. \end{align*} \]

Ejemplo (continuación)

En segundo lugar, resolvemos \(\mathbf{U}\mathbf{x}=\mathbf{y}\): \[ \begin{bmatrix}1&1&-3&-1&-2 \\0&1&5&1&5 \\0&0&16&5&12 \\0&0&0&4.5625&-2.25 \\0&0&0&0&1.931507 \\\end{bmatrix}\cdot\begin{bmatrix}x_1\\ x_2\\ x_3\\ x_4\\ x_5\end{bmatrix}=\begin{bmatrix}2 \\0 \\0 \\3 \\-1.575342 \\\end{bmatrix}. \]

Algoritmo \(LU\)

El resultado siguiente dice cuando es posible factorizar una matriz \(\mathbf{A}\) de la forma \(\mathbf{A}=\mathbf{L}\mathbf{U}\):

Teorema de factorización \(LU\). Consideremos el sistema de ecuaciones lineal \(\mathbf{A}\mathbf{x}=\mathbf{b}\), donde \(\mathbf{A}\) es la matriz del sistema. Si en el algoritmo de eliminación de Gauss no es necesario permutar ninguna fila, esto es, que \(a_{kk}^{(k)}\neq 0\), para \(k=1,\ldots,n-1\), entonces la matriz \(\mathbf{A}\) se puede factorizar de la forma \(\mathbf{A}=\mathbf{L}\mathbf{U}\), con \(\mathbf{L}\) triangular inferior con unos en la diagonal y \(\mathbf{U}\), triangular superior donde las matrices \(\mathbf{L}\) y \(\mathbf{U}\) tienen la expresión siguiente en función de los valores obtenidos en el proceso de eliminación de Gauss:

\[ \mathbf{A}=\mathbf{L}\cdot\mathbf{U}=\begin{bmatrix}1&0&0&\ldots&0 \\m_{21}&1&0&\ldots&0 \\\vdots&\vdots&\ddots&\vdots&\vdots \\m_{n1}&m_{n2}&m_{n3}&\ldots&1 \\\end{bmatrix}\cdot \begin{bmatrix}a_{11}^{(1)}&a_{12}^{(1)}&a_{13}^{(1)}&\ldots&a_{1n}^{(1)} \\0&a_{22}^{(2)}&a_{23}^{(2)}&\ldots&a_{2n}^{(2)} \\\vdots&\vdots&\ddots&\vdots&\vdots \\0&0&0&\ldots&a_{nn}^{(n)} \\\end{bmatrix}, \] donde \(m_{ji}=\frac{a_{ji}^{(i)}}{a_{ii}^{(i)}}\).

Demostración del Teorema

Sea \(\mathbf{A}^{(1)}=\mathbf{A}\).

Recordemos que en el primer paso del proceso de eliminación de Gauss, definíamos \(m_{j1}=\frac{a_{j1}^{(1)}}{a_{11}^{(1)}}\) y realizábamos las operaciones: \[ (E_j-m_{j1}E_1)\rightarrow (E_j), \] para \(j=2,\ldots, n\) obteniendo la matriz \(\mathbf{A}^{(2)}\):

\[ \mathbf{A}^{(2)}=\begin{bmatrix}a_{11}^{(1)}&a_{12}^{(1)}&\ldots&a_{1n}^{(1)} \\0&a_{22}^{(2)}&\ldots&a_{2n}^{(2)} \\\vdots&\vdots&\ddots&\vdots \\0&a_{n2}^{(2)}&\ldots&a_{nn}^{(2)} \\\end{bmatrix}. \]

Demostración del Teorema (continuación)

El paso anterior puede interpretarse desde otro punto de vista.

Si definimos la matriz: \[ \mathbf{M}^{(1)}=\begin{bmatrix}1&0&0&\ldots&0 \\-m_{21}&1&0&\ldots&0 \\\vdots&\vdots&\ddots&\ldots&\vdots \\-m_{n1}&0&0&\ldots&1 \\\end{bmatrix}, \] podemos escribir la matriz \(\mathbf{A}^{(2)}\) como: \(\mathbf{A}^{(2)}=\mathbf{M}^{(1)}\mathbf{A}^{(1)}=\mathbf{M}^{(1)}\mathbf{A}\).

Es sencillo comprobar que la inversa de \(\mathbf{M}^{(1)}\) vale: \[ \left(\mathbf{M}^{(1)}\right)^{-1}=\begin{bmatrix}1&0&0&\ldots&0 \\m_{21}&1&0&\ldots&0 \\\vdots&\vdots&\ddots&\ldots&\vdots \\m_{n1}&0&0&\ldots&1 \\\end{bmatrix}. \] Basta comprobar que \(\mathbf{M}^{(1)}\cdot \left(\mathbf{M}^{(1)}\right)^{-1}=\mathbf{I}\).

Demostración del Teorema (continuación)

En el segundo paso partíamos de la matriz \(\mathbf{A}^{(2)}\), calculábamos los valores \(m_{j2}=\frac{a_{j2}^{(2)}}{a_{22}^{(2)}}\) y hacíamos las operaciones: \[ (E_j-m_{j2}E_2)\rightarrow (E_j), \] para \(j=3,\ldots, n\) obteniendo la matriz \(\mathbf{A}^{(3)}\): (donde las ecuaciones \(E_j\) se refieren al sistema de ecuaciones con matriz \(\mathbf{A}^{(2)}\))

\[ \mathbf{A}^{(3)}=\begin{bmatrix}a_{11}^{(1)}&a_{12}^{(1)}&a_{13}^{(1)}&\ldots&a_{1n}^{(1)} \\0&a_{22}^{(2)}&a_{23}^{(2)}&\ldots&a_{2n}^{(2)} \\0&0&a_{33}^{(3)}&\ldots&a_{3n}^{(3)} \\\vdots&\vdots&\vdots&\ddots&\vdots \\0&0&a_{n3}^{(3)}&\ldots&a_{nn}^{(3)} \\\end{bmatrix}. \]

Demostración del Teorema (continuación)

Igual que antes, interpretemos este paso desde otro punto de vista. Definimos:

\[ \mathbf{M}^{(2)}=\begin{bmatrix}1&0&0&\ldots&0 \\0&1&0&\ldots&0 \\0&-m_{32}&1&\ldots&0 \\\vdots&\vdots&\vdots&\ddots&\vdots \\0&-m_{n2}&0&\ldots&1 \\\end{bmatrix}, \] cuya inversa es: \[ \left(\mathbf{M}^{(2)}\right)^{-1}=\begin{bmatrix}1&0&0&\ldots&0 \\0&1&0&\ldots&0 \\0&m_{32}&1&\ldots&0 \\\vdots&\vdots&\vdots&\ddots&\vdots \\0&m_{n2}&0&\ldots&1 \\\end{bmatrix}. \] Podemos escribir la matriz \(\mathbf{A}^{(3)}\) como \(\mathbf{A}^{(3)}=\mathbf{M}^{(2)}\cdot \mathbf{A}^{(2)}=\mathbf{M}^{(2)}\mathbf{M}^{(1)}\mathbf{A}\).

Demostración del Teorema (continuación)

Imaginemos que estamos en el paso \(k\)-ésimo del proceso de eliminación de Gauss. Entonces tendremos la matriz:

\[ \mathbf{A}^{(k)}=\begin{bmatrix}a_{11}&\ldots&\ldots&\ldots&\ldots&a_{1n} \\\vdots&\vdots&\ddots&\vdots&\vdots&\vdots \\0&\ldots&0&a_{kk}^{(k)}&\ldots&a_{kn}^{(k)} \\\vdots&\vdots&\vdots&\vdots&\ddots&\vdots \\0&\ldots&0&a_{nk}^{(k)}&\ldots&a_{nn}^{(k)} \\\end{bmatrix}. \]

Demostración del Teorema (continuación)

En este paso, calculábamos los valores \(m_{jk}=\frac{a_{jk}^{(k)}}{a_{kk}^{(k)}}\) y hacíamos las operaciones: \[ (E_j-m_{jk}E_k)\rightarrow (E_j), \] para \(j=k+1,\ldots, n\) obteniendo la matriz \(\mathbf{A}^{(k+1)}\): (donde las ecuaciones \(E_j\) se refieren al sistema de ecuaciones con matriz \(\mathbf{A}^{(k)}\))

\[ \mathbf{A}^{(k+1)}=\begin{bmatrix}a_{11}&\ldots&\ldots&\ldots&\ldots&a_{1n} \\\vdots&\vdots&\ddots&\vdots&\vdots&\vdots \\0&\ldots&0&a_{k+1,k+1}^{(k+1)}&\ldots&a_{k+1,n}^{(k+1)} \\\vdots&\vdots&\vdots&\vdots&\vdots&\vdots \\0&\ldots&0&a_{n,k+1}^{(k+1)}&\ldots&a_{n,n}^{(k+1)} \\\end{bmatrix}. \]

Demostración del Teorema (continuación)

Igual que antes, interpretemos este paso desde otro punto de vista. Definimos:

\[ \mathbf{M}^{(k)}=\begin{bmatrix}1&0&0&0&\ldots&0 \\\vdots&\ddots&0&0&\ldots&0 \\0&\ldots&1&0&\ldots&0 \\0&\ldots&-m_{k+1,k}&1&\ldots&0 \\\vdots&\vdots&\vdots&\vdots&\ddots&\vdots \\0&\ldots&-m_{n,k}&0&\ldots&1 \\\end{bmatrix}, \] cuya inversa es: \[ \left(\mathbf{M}^{(k)}\right)^{-1}=\begin{bmatrix}1&0&0&0&\ldots&0 \\\vdots&\ddots&0&0&\ldots&0 \\0&\ldots&1&0&\ldots&0 \\0&\ldots&m_{k+1,k}&1&\ldots&0 \\\vdots&\vdots&\vdots&\vdots&\ddots&\vdots \\0&\ldots&m_{n,k}&0&\ldots&1 \\\end{bmatrix}. \] Podemos escribir la matriz \(\mathbf{A}^{(k+1)}\) como \(\mathbf{A}^{(k+1)}=\mathbf{M}^{(k)}\cdot \mathbf{A}^{(k)}=\mathbf{M}^{(k)}\mathbf{M}^{(k-1)}\cdots \mathbf{M}^{(1)}\mathbf{A}\).

Demostración del Teorema (continuación)

En el último paso tendremos la expresión siguiente: \[ \mathbf{A}^{(n)}=\mathbf{M}^{(n-1)}\mathbf{M}^{(n-2)}\cdots \mathbf{M}^{(1)}\mathbf{A}, \] donde recordemos que la matriz \(\mathbf{A}^{(n)}\) es triangular superior.

La expresión anterior puede escribirse de la siguiente manera despejando la matriz original del sistema \(\mathbf{A}\): \[ \mathbf{A}=\left(\mathbf{M}^{(1)}\right)^{-1}\left(\mathbf{M}^{(2)}\right)^{-1}\cdots \left(\mathbf{M}^{(n-1)}\right)^{-1}\mathbf{A}^{(n)}. \] Recordemos que las matrices \(\left(\mathbf{M}^{(1)}\right)^{-1}\) eran triangulares inferiores con unos en la diagonal. Como el producto de matrices triangulares inferiores con unos en la diagonal sigue siendo una matriz triangular inferior con unos en la diagonal (ejercicio) tenemos que las matrices \(\mathbf{L}\) y \(\mathbf{U}\) serán:

\[ \mathbf{L}=\left(\mathbf{M}^{(1)}\right)^{-1}\left(\mathbf{M}^{(2)}\right)^{-1}\cdots \left(\mathbf{M}^{(n-1)}\right)^{-1},\ \mathbf{U}=\mathbf{A}^{(n)}. \]

Demostración del Teorema (continuación)

Para acabar la demostración faltaría ver que: \[ \mathbf{L}=\begin{bmatrix}1&0&0&\ldots&0 \\m_{21}&1&0&\ldots&0 \\\vdots&\vdots&\ddots&\vdots&\vdots \\m_{n1}&m_{n2}&m_{n3}&\ldots&1 \\\end{bmatrix}. \]

Ejercicio

Demostrar la última expresión de la matriz \(\mathbf{L}\).

Indicación: Multiplicar primero \(\left(\mathbf{M}^{(1)}\right)^{-1}\left(\mathbf{M}^{(2)}\right)^{-1}\) usando las expresiones de \(\left(\mathbf{M}^{(1)}\right)^{-1}\) y \(\left(\mathbf{M}^{(2)}\right)^{-1}\) vistas en esta demostración para ver cómo se va comportando el producto para después generalizar el resultado.

Algoritmo para calcular las matrices \(\mathbf{L}\) y \(\mathbf{U}\)

Una manera de calcular las matrices \(\mathbf{L}\) y \(\mathbf{U}\) es usar el algoritmo de eliminación de Gauss y aplicar las expresiones vistas en el resultado anterior.

Sin embargo, este método es tedioso y requiere realizar todo el proceso de eliminación de Gauss. Veamos un método mucho más directo.

Recordemos que tenemos una matriz \(\mathbf{A}\), \(n\times n\) tal que se pueden realizar todos los pasos del algoritmo de eliminación de Gauss sin permutar filas y queremos hallar las matrices \(\mathbf{L}\) triangular inferior con unos en la diagonal y \(\mathbf{U}\) triangular superior tal que \(\mathbf{A}=\mathbf{L}\mathbf{U}.\)

Algoritmo para calcular las matrices \(\mathbf{L}\) y \(\mathbf{U}\)

Para hallar dichas matrices, realizamos los pasos siguientes:

La primera fila de la matriz \(\mathbf{U}\) coincide con la primera fila de la matriz \(\mathbf{A}\): \(u_{1j}=a_{1j}\), \(j=2,\ldots,n\).
Calculamos la primera columna de la matriz \(\mathbf{L}\): \(l_{j1}=\frac{a_{j1}}{u_{11}}\).
A continuación, calculamos la segunda fila de la matriz \(\mathbf{U}\) y a continuación la segunda columna de la matriz \(\mathbf{L}\) con \(l_{22}=1\): \[ \begin{align*} u_{2j} & =a_{2j}-l_{21}u_{12},\ j=2,\ldots,n,\\ l_{j2} & =\frac{1}{u_{22}}\left(a_{j2}-l_{j1}u_{12}\right),\ j=3,\ldots,n. \end{align*} \]

Algoritmo para calcular las matrices \(\mathbf{L}\) y \(\mathbf{U}\)

Calculamos la tercera fila de la matriz \(\mathbf{U}\) y a continuación la tercera columna de la matriz \(\mathbf{L}\) con \(l_{33}=1\): \[ \begin{align*} u_{3j} & =a_{3j}-l_{31}u_{1j}-l_{32}u_{2j},\ j=3,\ldots,n,\\ l_{j3} & =\frac{1}{u_{33}}\left(a_{j3}-l_{j1}u_{13}-l_{j2}u_{23}\right),\ j=4,\ldots,n. \end{align*} \]

Algoritmo para calcular las matrices \(\mathbf{L}\) y \(\mathbf{U}\)

En general, supongamos que ya hemos hallado las \(i-1\) primeras filas de la matriz \(\mathbf{U}\) y las \(i-1\) primeras columnas de la matriz \(\mathbf{L}\). Entonces, hallamos la fila \(i\)-ésima de la matriz \(\mathbf{U}\) y a continuación la columna \(i\)-ésima de la matriz \(\mathbf{L}\) con \(l_{ii}=1\): \[ \begin{align*} u_{ij} & =a_{ij}-\sum_{k=1}^{i-1}l_{ik}u_{kj},\ j=i,\ldots,n,\\ l_{ji} & =\frac{1}{u_{ii}}\left(a_{ji}-\sum_{k=1}^{i-1}l_{jk}u_{ki}\right),\ j=i+1,\ldots,n. \end{align*} \]
Al final, calculamos la fila \(n\)-ésima de la matriz \(\mathbf{U}\) y a continuación la columna \(n\)-ésima de la matriz \(\mathbf{L}\) con \(l_{nn}=1\): \[ u_{nn} =a_{nn}-\sum_{k=1}^{n-1}l_{nk}u_{kn}. \]

Factorización \(LU\). Pseudocódigo

INPUT matriz del sistema \(\mathbf{A}=(a_{ij})_{i=1,\ldots,n,j=1,\ldots,n}\).
Set \(\mathbf{L}=\mathbf{0}\).
Set \(\mathbf{U}=\mathbf{0}\).
For i=1,...n
- Set \(l_{ii}=1\)
- Set \(u_{1i}=a_{1i}\) (Calculamos la primera fila de \(\mathbf{U}\))
- Set\(l_{i1}=\frac{a_{i1}}{u_{11}}\) (Calculamos la primera columna de \(\mathbf{L}\))

Factorización \(LU\). Pseudocódigo

For i=2,...n (Vamos a calcular la fila \(i\) de la matriz \(\mathbf{U}\) y a continuación la columna \(i\) de la matriz \(\mathbf{L}\))
- Set \(u_{ii}=a_{ii}-\sum_{k=1}^{i-1}l_{ik}u_{ki}\).
- For j=i+1,...,n
  - Set \(u_{ij}=a_{ij}-\sum_{k=1}^{i-1}l_{ik}u_{kj}\) (fila \(i\)-ésima de \(\mathbf{U}\))
  - Set \(l_{ji}=\frac{1}{u_{ii}}\left(a_{ji}-\sum_{k=1}^{i-1}l_{jk}u_{ki}\right)\) (columna \(i\)-ésima de \(\mathbf{L}\))
Print \(\mathbf{L}\), \(\mathbf{U}\) (Damos las matrices \(\mathbf{L}\) y \(\mathbf{U}\))

Ejemplo anterior

Consideremos el sistema lineal siguiente: \[ \begin{align*} E_1: & x_1 & + x_2 & -3 x_3 & - x_4 & -2 x_5 & = 2,\\ E_2: & x_1 & +2 x_2 & +2 x_3 & & +3 x_5 & = 2,\\ E_3: & x_1 & - x_2 & +3 x_3& +2 x_4& & = 2,\\ E_4: & & x_2 & & +4 x_4 & - x_5 & = 3,\\ E_5: & x_1 & +3 x_2& + x_3 & & +5 x_5 & = 1. \end{align*} \] La matriz del sistema es la siguiente: \[ \mathbf{A}=\begin{bmatrix}1&1&-3&-1&-2&2 \\1&2&2&0&3&2 \\1&-1&3&2&0&2 \\0&1&0&4&-1&3 \\1&3&1&0&5&1 \\\end{bmatrix}. \]

Ejemplo anterior (continuación)

Las matrices resultantes de aplicar el algoritmo de eliminación de Gauss junto con los valores \(m_{ik}\) que vamos escribiendo en la matriz \(\mathbf{L}\) son los siguientes:

Paso 1: \[ \mathbf{A}^{(2)}=\begin{bmatrix}1&1&-3&-1&-2 \\0&1&5&1&5 \\0&-2&6&3&2 \\0&1&0&4&-1 \\0&2&4&1&7 \\\end{bmatrix},\quad \mathbf{L}=\begin{bmatrix}1&0&0&0&0 \\1&1&0&0&0 \\1&*&1&0&0 \\0&*&*&1&0 \\1&*&*&*&1 \\\end{bmatrix} \]
Paso 2: \[ \mathbf{A}^{(3)}=\begin{bmatrix}1&1&-3&-1&-2 \\0&1&5&1&5 \\0&0&16&5&12 \\0&0&-5&3&-6 \\0&0&-6&-1&-3 \\\end{bmatrix},\quad \mathbf{L}=\begin{bmatrix}1&0&0&0&0 \\1&1&0&0&0 \\1&-2&1&0&0 \\0&1&*&1&0 \\1&2&*&*&1 \\\end{bmatrix} \]

Ejemplo anterior (continuación)

Paso 3: \[ \mathbf{A}^{(4)}=\begin{bmatrix}1&1&-3&-1&-2 \\0&1&5&1&5 \\0&0&16&5&12 \\0&0&0&4.5625&-2.25 \\0&0&0&0.875&1.5 \\\end{bmatrix},\quad \mathbf{L}=\begin{bmatrix}1&0&0&0&0 \\1&1&0&0&0 \\1&-2&1&0&0 \\0&1&-0.3125&1&0 \\1&2&-0.375&*&1 \\\end{bmatrix} \]
Paso 4: \[ \mathbf{A}^{(5)}=\begin{bmatrix}1&1&-3&-1&-2 \\0&1&5&1&5 \\0&0&16&5&12 \\0&0&0&4.5625&-2.25 \\0&0&0&0&1.931507 \\\end{bmatrix},\quad \mathbf{L}=\begin{bmatrix}1&0&0&0&0 \\1&1&0&0&0 \\1&-2&1&0&0 \\0&1&-0.3125&1&0 \\1&2&-0.375&0.191781&1 \\\end{bmatrix} \]

Ejemplo anterior (continuación)

La matriz \(\mathbf{U}\) será: \[ \mathbf{U}=\mathbf{A}^{(5)}=\begin{bmatrix}1&1&-3&-1&-2 \\0&1&5&1&5 \\0&0&16&5&12 \\0&0&0&4.5625&-2.25 \\0&0&0&0&1.931507 \\\end{bmatrix}. \]

Ejemplo anterior (continuación)

Calculemos las matrices \(\mathbf{L}\) y \(\mathbf{U}\) aplicando el algoritmo.

Primera fila de \(\mathbf{U}\) coincide con la primera fila de \(\mathbf{A}\).

\[ \mathbf{A}=\begin{bmatrix}1&1&-3&-1&-2 \\1&2&2&0&3 \\1&-1&3&2&0 \\0&1&0&4&-1 \\1&3&1&0&5 \\\end{bmatrix},\Rightarrow \mathbf{U}=\begin{bmatrix}1&1&-3&-1&-2 \\0&*&*&*&* \\0&0&*&*&* \\0&0&0&*&* \\0&0&0&0&* \\\end{bmatrix}. \]

Primera columna de \(\mathbf{L}\) coincide con la primera columna \(\mathbf{A}\) dividido por \(u_{11}=1\): \[ \mathbf{A}=\begin{bmatrix}1&1&-3&-1&-2 \\1&2&2&0&3 \\1&-1&3&2&0 \\0&1&0&4&-1 \\1&3&1&0&5 \\\end{bmatrix},\Rightarrow \mathbf{L}=\begin{bmatrix}1&0&0&0&0 \\1&1&0&0&0 \\1&*&1&0&0 \\0&*&*&1&0 \\1&*&*&*&1 \\\end{bmatrix}. \]

Ejemplo anterior (continuación)

Segunda fila de \(\mathbf{U}\): \(u_{2j}=a_{2j}-l_{21}u_{1j}=a_{2j}-1\cdot u_{1j}\): \[ \begin{align*} u_{22} & =a_{22}-l_{21}u_{12}=2-1\cdot 1=1,\\ u_{23} & =a_{23}-l_{21}u_{13}=2-1\cdot (-3)=5,\\ u_{24} & =a_{24}-l_{21}u_{14}=0-1\cdot (-1)=1,\\ u_{25} & =a_{25}-l_{21}u_{15}=3-1\cdot (-2)=5. \end{align*} \]

\[ \mathbf{U}=\begin{bmatrix}1&1&-3&-1&-2 \\0&1&5&1&5 \\0&0&*&*&* \\0&0&0&*&* \\0&0&0&0&* \\\end{bmatrix}. \]

Ejemplo anterior (continuación)

Segunda columna de \(\mathbf{L}\): \(l_{j2}=\frac{1}{u_{22}}(a_{j2}-l_{j1}u_{12})=\frac{1}{1}(a_{j2}-l_{j1}\cdot 1)=a_{j2}-l_{j1}\): \[ \begin{align*} l_{32} & = a_{32}-l_{31}=-1-1=-2,\\ l_{42} & = a_{42}-l_{41}=1-0=1,\\ l_{52} & = a_{52}-l_{51}=3-1=2. \end{align*} \] \[ \mathbf{L}=\begin{bmatrix}1&0&0&0&0 \\1&1&0&0&0 \\1&-2&1&0&0 \\0&1&*&1&0 \\1&2&*&*&1 \\\end{bmatrix}. \]

Ejemplo anterior (continuación)

Tercera fila de \(\mathbf{U}\): \(u_{3j}=a_{3j}-l_{31}u_{1j}-l_{32}u_{2j}=a_{3j}-1\cdot u_{1j}-(-2)u_{2j}\): \[ \begin{align*} u_{33}& =a_{33}-l_{31}u_{13}-l_{32}u_{23}=3-1\cdot (-3)-(-2)\cdot 5=16,\\ u_{34}& =a_{34}-l_{31}u_{14}-l_{32}u_{24}=2-1\cdot (-1)-(-2)\cdot 1=5,\\ u_{35}& =a_{35}-l_{31}u_{15}-l_{32}u_{25}=0-1\cdot (-2)-(-2)\cdot 5=12. \end{align*} \]

\[ \mathbf{U}=\begin{bmatrix}1&1&-3&-1&-2 \\0&1&5&1&5 \\0&0&16&5&12 \\0&0&0&*&* \\0&0&0&0&* \\\end{bmatrix}. \]

Ejemplo anterior (continuación)

Tercera columna de \(\mathbf{L}\): \(l_{j3}=\frac{1}{u_{33}}(a_{j3}-l_{j1}u_{13}-l_{j2}u_{23})=\frac{1}{16}(a_{j3}-l_{j1}\cdot (-3)-l_{j2}\cdot 5)\): \[ \begin{align*} l_{43} & =\frac{1}{u_{33}}(a_{43}-l_{41}u_{13}-l_{42}u_{23}) = \frac{1}{16}(0-0\cdot (-3)-1\cdot 5)=-0.3125,\\ l_{53} & =\frac{1}{u_{33}}(a_{53}-l_{41}u_{13}-l_{42}u_{23}) = \frac{1}{16}(1-1\cdot (-3)-2\cdot 5)=-0.375. \end{align*} \] \[ \mathbf{L}=\begin{bmatrix}1&0&0&0&0 \\1&1&0&0&0 \\1&-2&1&0&0 \\0&1&-0.3125&1&0 \\1&2&-0.375&*&1 \\\end{bmatrix}. \]

Ejemplo anterior (continuación)

Cuarta fila de \(\mathbf{U}\): \(u_{4j}=a_{4j}-l_{41}u_{1j}-l_{42}u_{2j}-l_{43}u_{3j}=a_{4j}-0\cdot u_{1j}-1\cdot u_{2j}-(-0.3125)u_{3j}\): \[ \begin{align*} u_{44} & =a_{44}-l_{41}u_{14}-l_{42}u_{24}-l_{43}u_{34}=4-0\cdot (-1)-1\cdot 1-(-0.3125)\cdot 5= 4.5625,\\ u_{45} & =a_{45}-l_{41}u_{15}-l_{42}u_{25}-l_{43}u_{35}=-1-0\cdot (-2)-1\cdot 5-(-0.3125)\cdot 12= -2.25. \end{align*} \]

\[ \mathbf{U}=\begin{bmatrix}1&1&-3&-1&-2 \\0&1&5&1&5 \\0&0&16&5&12 \\0&0&0&4.5625&-2.25 \\0&0&0&0&* \\\end{bmatrix}. \]

Ejemplo anterior (continuación)

Cuarta columna de \(\mathbf{L}\) o \(l_{54}\): \[ \begin{align*} l_{54} & =\frac{1}{u_{44}}(a_{54}-l_{51}u_{14}-l_{52}u_{24}-l_{53}u_{34})\\ & =\frac{1}{4.5625}(0-1\cdot (-1)-2\cdot 1-(-0.375)\cdot 5)=0.191781. \end{align*} \] \[ \mathbf{L}=\begin{bmatrix}1&0&0&0&0 \\1&1&0&0&0 \\1&-2&1&0&0 \\0&1&-0.3125&1&0 \\1&2&-0.375&0.191781&1 \\\end{bmatrix}. \]

Ejemplo anterior (continuación)

Por último, quinta fila de \(\mathbf{U}\) o \(u_{55}\): \[ \begin{align*} u_{55} & =a_{55}-l_{51}u_{15}-l_{52}u_{25}-l_{53}u_{35}-l_{54}u_{45} \\ & =5-1\cdot (-2)-2\cdot 5-(-0.375)\cdot 12-0.1917808\cdot (-2.25)=1.9315068: \end{align*} \] \[ \mathbf{U}=\begin{bmatrix}1&1&-3&-1&-2 \\0&1&5&1&5 \\0&0&16&5&12 \\0&0&0&4.5625&-2.25 \\0&0&0&0&1.931507 \\\end{bmatrix}. \]

Propiedades de la descomposición \(LU\)

Proposición. Sea \(\mathbf{A}\) una matriz \(n\times n\) que admite descomposición \(LU\). Sean \(\mathbf{L}\) y \(\mathbf{U}\) las matrices triangular inferior con unos en la diagonal y triangular superior, respectivamente tal que \(\mathbf{A}=\mathbf{L}\cdot\mathbf{U}\).

Sea \(k\) un entero entre \(1\) y \(n\). Sea \(\mathbf{A}_k\) la submatriz de \(\mathbf{A}\) formada por las \(k\) primeras filas y las \(k\) primeras columnas. De la misma manera, definimos \(\mathbf{L}_k\) y \(\mathbf{U}_k\). Entonces \(\mathbf{A}_k=\mathbf{L}_k\cdot \mathbf{U}_k\), es decir las submatrices \(\mathbf{L}_k\) y \(\mathbf{U}_k\) serían la descomposición \(LU\) de la submatriz \(\mathbf{A}_k\).

Demostración de la proposición

Escribimos \(\mathbf{A}=(a_{ij})_{i,j=1,\ldots,n}\), \(\mathbf{L}=(l_{ij})_{i,j=1,\ldots,n}\) y \(\mathbf{U}=(u_{ij})_{i,j=1,\ldots,n}\).

Sea un elemento \(a_{ij}\) de la submatriz \(\mathbf{A}_k\). Como \(\mathbf{A}=\mathbf{L}\cdot\mathbf{U}\), podemos escribir: \(\displaystyle a_{ij}=\sum_{p=1}^n l_{ip}\cdot u_{pj}.\)

Como la matriz \(\mathbf{L}\) es triangular inferior \(l_{ip}=0\), \(p=i+1,\ldots,n\) y como \(\mathbf{U}\) es triangular superior, \(u_{pj}=0\), \(p=j+1,\ldots,n\). Sea \(m=\max\{i,j\}\), entonces \(l_{ip}=u_{pj}=0\), para \(p=m+1,\ldots,n\). Por tanto, \(\displaystyle a_{ij}=\sum_{p=1}^m l_{ip}\cdot u_{pj}\). Ahora bien, como \(a_{ij}\) es un elemento de la submatriz \(\mathbf{A}_k\), entonces \(m=\max\{i,j\}\leq k\) y por tanto: \(\displaystyle a_{ij}=\sum_{p=1}^k l_{ip}\cdot u_{pj}\).

En resumen, hemos demostrado que \(\mathbf{A}_k=\mathbf{L}_k\cdot\mathbf{U}_k.\)

Ejemplo de aplicación

En un ejemplo anterior, considerábamos la matriz:

\[ \mathbf{A}=\begin{bmatrix}1&1&-3&-1&-2 \\1&2&2&0&3 \\1&-1&3&2&0 \\0&1&0&4&-1 \\1&3&1&0&5 \\\end{bmatrix}, \] cuya descomposición \(LU\) era la siguiente: \[ \mathbf{L}=\begin{bmatrix}1&0&0&0&0 \\1&1&0&0&0 \\1&-2&1&0&0 \\0&1&-0.3125&1&0 \\1&2&-0.375&0.191781&1 \\\end{bmatrix},\quad \mathbf{U}=\begin{bmatrix}1&1&-3&-1&-2 \\0&1&5&1&5 \\0&0&16&5&12 \\0&0&0&4.5625&-2.25 \\0&0&0&0&1.931507 \\\end{bmatrix}. \]

Ejemplo de aplicación (continuación)

\(k=1\). La submatriz \(A_1\) vale \(A_1=1\) que trivialmente vale: \[ \mathbf{A}_1=\mathbf{L}_1\cdot \mathbf{U}_1,\ \Rightarrow 1=1\cdot 1. \]
\(k=2\). La submatriz \(A_2\) vale \(A_2=\begin{bmatrix}1&1 \\1&2 \\\end{bmatrix}\) y se verifica que: \[ \mathbf{A}_2=\mathbf{L}_2\cdot \mathbf{U}_2,\ \Rightarrow \begin{bmatrix}1&1 \\1&2 \\\end{bmatrix}=\begin{bmatrix}1&0 \\1&1 \\\end{bmatrix}\cdot \begin{bmatrix}1&1 \\0&1 \\\end{bmatrix}. \]
\(k=3\). La submatriz \(A_3\) vale \(A_3=\begin{bmatrix}1&1&-3 \\1&2&2 \\1&-1&3 \\\end{bmatrix}\) y se verifica que: \[ \mathbf{A}_3=\mathbf{L}_3\cdot \mathbf{U}_3,\ \Rightarrow \begin{bmatrix}1&1&-3 \\1&2&2 \\1&-1&3 \\\end{bmatrix}=\begin{bmatrix}1&0&0 \\1&1&0 \\1&-2&1 \\\end{bmatrix}\cdot \begin{bmatrix}1&1&-3 \\0&1&5 \\0&0&16 \\\end{bmatrix}. \]

Ejemplo de aplicación (continuación)

\(k=4\). La submatriz \(A_4\) vale \(A_4=\begin{bmatrix}1&1&-3&-1 \\1&2&2&0 \\1&-1&3&2 \\0&1&0&4 \\\end{bmatrix}\) y se verifica que: \[ \mathbf{A}_4=\mathbf{L}_4\cdot \mathbf{U}_4,\ \Rightarrow \begin{bmatrix}1&1&-3&-1 \\1&2&2&0 \\1&-1&3&2 \\0&1&0&4 \\\end{bmatrix}=\begin{bmatrix}1&0&0&0 \\1&1&0&0 \\1&-2&1&0 \\0&1&-0.3125&1 \\\end{bmatrix}\cdot \begin{bmatrix}1&1&-3&-1 \\0&1&5&1 \\0&0&16&5 \\0&0&0&4.5625 \\\end{bmatrix}. \]
\(k=5\): no hace falta comprobar nada ya que la condición a comprobar es la condición de descomposición \(LU\) de la matriz \(\mathbf{A}\).

Propiedades de la descomposición \(LU\)

Sea \(k\) un entero entre \(1\) y \(n\). Sea \(\mathbf{A}_k\) la submatriz de \(\mathbf{A}\) formada por las \(k\) primeras filas y las \(k\) primeras columnas.

Entonces, el determinante de \(\mathbf{A}_k\) vale \(\mathrm{det}(\mathbf{A}_k)=a_{11}^{(1)}\cdots a_{kk}^{(k)}\), donde \(a_{ii}^{(i)}\) es el pivote del paso \(i\)-ésimo del algoritmo de descomposición de Gauss.

Demostración de la proposición

Demostración

Usando la proposición anterior tenemos que \(\mathbf{A}_k=\mathbf{L}_k\cdot\mathbf{U}_{k}\) donde \(\mathbf{L}_k\) y \(\mathbf{U}_k\) recordemos que son las submatrices formadas por las \(k\) primeras filas y columnas.

Usando un teorema anterior tenemos que las matrices \(\mathbf{L}_k\) y \(\mathbf{U}_k\) son las siguientes:

\[ \mathbf{L}_k=\begin{bmatrix}1&0&0&\ldots&0 \\m_{21}&1&0&\ldots&0 \\\vdots&\vdots&\ddots&\vdots&\vdots \\m_{k1}&m_{k2}&m_{k3}&\ldots&1 \\\end{bmatrix}\quad \mathbf{U}_k = \begin{bmatrix}a_{11}^{(1)}&a_{12}^{(1)}&a_{13}^{(1)}&\ldots&a_{1k}^{(1)} \\0&a_{22}^{(2)}&a_{23}^{(2)}&\ldots&a_{2k}^{(2)} \\\vdots&\vdots&\ddots&\vdots&\vdots \\0&0&0&\ldots&a_{kk}^{(k)} \\\end{bmatrix}. \] Por tanto, \[ \mathrm{det}(\mathbf{A}_k)=\mathrm{det}(\mathbf{L}_k)\cdot\mathrm{det}(\mathbf{U}_k)=a_{11}^{(1)}\cdots a_{kk}^{(k)}, \] tal como queríamos ver.

Matrices de permutación

Hemos visto anteriormente que no siempre existe la descomposición \(\mathbf{L}\mathbf{U}\) para una matriz \(\mathbf{A}\) de un sistema lineal.

Recordemos que para que dicha descomposición exista, es necesario que \(a_{kk}^{(k)}\neq 0\), siendo \(a_{kk}^{(k)}\) el elemento diagonal de la matriz \(\mathbf{A}^{(k)}\) del algoritmo de eliminación de Gauss en el paso \(k\)-ésimo.

Cuando explicamos el algoritmo de eliminación de Gauss, si en el paso \(k\), el elemento \(a_{kk}^{(k)}=0\), necesitábamos realizar una permutación de filas para poder continuar con el algoritmo.

¿A qué equivale realizar una permutación de filas matricialmente hablando?

Para contestar a la pregunta anterior, necesitamos introducir las matrices de permutación.

Matrices de permutación

Definición de matriz de permutación. Diremos que una matriz \(\mathbf{P}_{ij}=(p_{kl})_{k,l=1,\ldots,n}\), \(n\times n\), es una matriz de permutación si \[ p_{kl}=\begin{cases} 1, & \mbox{si }k=l\neq i,j,\\ 1, & \mbox{si }k=i,\ l=j,\\ 1, & \mbox{si }k=j,\ l=i,\\ 0, &\mbox{en caso contrario.} \end{cases} \]

Ejemplo

Consideremos \(n=5\). A continuación mostramos algunas matrices de permutación:

\[ \mathbf{P}_{14}=\begin{bmatrix}0&0&0&1&0 \\0&1&0&0&0 \\0&0&1&0&0 \\1&0&0&0&0 \\0&0&0&0&1 \\\end{bmatrix},\ \mathbf{P}_{23}=\begin{bmatrix}1&0&0&0&0 \\0&0&1&0&0 \\0&1&0&0&0 \\0&0&0&1&0 \\0&0&0&0&1 \\\end{bmatrix},\ \mathbf{P}_{35}=\begin{bmatrix}1&0&0&0&0 \\0&1&0&0&0 \\0&0&0&0&1 \\0&0&0&1&0 \\0&0&1&0&0 \\\end{bmatrix}. \]

Matrices de permutación

Veamos qué relación existe entre las matrices de permutación y la permutación de filas:

Proposición. Sea \(\mathbf{A}\) una matriz \(n\times n\). Entonces la matriz \(\mathbf{A}_{ij}\) resultante de permutar las filas \(i\) y \(j\) vale: \[ \mathbf{A}_{ij}=\mathbf{P}_{ij}\cdot\mathbf{A}. \]

Demostración

Llamamos \(\mathbf{A}=(a_{kl})_{k,l=1,\ldots,n}\), \(\mathbf{P}_{ij}=(p_{kl})_{k,l=1,\ldots,n}\) y \(\mathbf{A}_{ij}=(a_{kl}^{(ij)})_{k,l=1,\ldots,n}\) a las componentes de la matriz \(\mathbf{A}\), \(\mathbf{P}\) y \(\mathbf{A}_{ij}\), respectivamente.

Sea \(k\) una fila de la matriz \(\mathbf{A}_{ij}\) donde \(k\neq i,j\), en este caso, podemos escribir: \[ a_{kl}^{(ij)}=\sum_{m=1}^n p_{km}a_{ml}=a_{kl}, \] ya que \(p_{km}=0\), si \(m\neq k\) y \(p_{kk}=1\) para \(l=1,\ldots,n\). En resumen, la fila \(k\)-ésima de la matriz \(\mathbf{A}_{ij}\) coincide con la fila \(k\)-ésima de la matriz \(\mathbf{A}\).

Matrices de permutación

Demostración (continuación)

Sea \(k=i\) la fila \(i\)-ésima de la matriz \(\mathbf{A}_{ij}\). En este caso \[ a_{il}^{(ij)}=\sum_{m=1}^n p_{im}a_{ml}=a_{jl}, \] ya que \(p_{im}=0\), si \(m\neq j\) y \(p_{ij}=1\) para \(l=1,\ldots,n\). En resumen, la fila \(i\)-ésima de la matriz \(\mathbf{A}^{(ij)}\) coincide con la fila \(j\)-ésima de la matriz \(\mathbf{A}\).

Usando un razonamiento idéntico, podemos ver que la fila \(j\)-ésima de la matriz \(\mathbf{A}_{ij}\) coincide con la fila \(i\)-ésima de la matriz \(\mathbf{A}\).

Acabamos de demostrar que la matriz \(\mathbf{P}_{ij}\cdot\mathbf{A}\) coincide con la matriz \(\mathbf{A}\) en la que hemos permutado las filas \(i\) y \(j\), tal como queríamos demostrar.

Matrices de permutación

Las matrices de permutación tienen las propiedades siguientes:

La matriz de permutación \(\mathbf{P}_{ij}\) es simétrica.
La matriz de permutación \(\mathbf{P}_{ij}\) es invertible y además \(\mathbf{P}_{ij}^{-1}=\mathbf{P}_{ij}\).

Ejercicio

Demostrar las dos propiedades anteriores.

Matrices de permutación y descomposición \(LU\)

Sea \(\mathbf{A}\) la matriz de un sistema lineal de la que no es posible hallar la descomposición \(LU\).

Sin embargo, existe una permutación de filas tal que es posible hallar la descomposición \(LU\) de la matriz resultante de permutar las filas.

Esto es equivalente a decir que existe una matriz de permutación \(\mathbf{P}_{ij}\) tal que existe la descomposición \(LU\) de la matriz \(\mathbf{P}_{ij}\mathbf{A}\).

Recordemos que la descomposición \(LU\) se usaba para resolver un sistema lineal \(\mathbf{A}\mathbf{x}=\mathbf{b}\). Como no podemos descomponer la matriz \(\mathbf{A}\), multiplicando por la matriz de permutación \(\mathbf{P}_{ij}\) la expresión anterior, tenemos que resolver el sistema anterior es equivalente a resolver: \(\mathbf{P}_{ij}\mathbf{A}\mathbf{x}=\mathbf{P}_{ij}\mathbf{b}\).

A continuación, descomponemos la matriz \(\mathbf{P}_{ij}\mathbf{A}=\mathbf{L}\mathbf{U}\) y resolvemos el sistema anterior. Lo único que cambia es que el vector independiente \(\mathbf{b}\) pasa a ser \(\mathbf{P}_{ij}\mathbf{b}\).

Ejemplo

Nos planteamos resolver el sistema anterior pero hemos cambiado el valor \(a_{11}\) por el valor \(0\):

\[ \begin{align*} E_1: & & x_2 & -3 x_3 & - x_4 & -2 x_5 & = 2,\\ E_2: & x_1 & +2 x_2 & +2 x_3 & & +3 x_5 & = 2,\\ E_3: & x_1 & - x_2 & +3 x_3& +2 x_4& & = 2,\\ E_4: & & x_2 & & +4 x_4 & - x_5 & = 3,\\ E_5: & x_1 & +3 x_2& + x_3 & & +5 x_5 & = 1. \end{align*} \] La matriz del sistema será: \[ \mathbf{A}=\begin{bmatrix}0&1&-3&-1&-2 \\1&2&2&0&3 \\1&-1&3&2&0 \\0&1&0&4&-1 \\1&3&1&0&5 \\\end{bmatrix}. \]

Ejemplo (continuación)

Como \(a_{11}=0\), no podemos realizar la descomposición \(LU\).

Entonces nos planteamos hallar la descomposición \(LU\) de la matriz \(\mathbf{P}_{15}\mathbf{A}=\begin{bmatrix}1&3&1&0&5 \\1&2&2&0&3 \\1&-1&3&2&0 \\0&1&0&4&-1 \\0&1&-3&-1&-2 \\\end{bmatrix}:\) que resulta de permutar las filas \(1\) y \(5\): \[ \mathbf{L}=\begin{bmatrix}1&0&0&0&0 \\1&1&0&0&0 \\1&4&1&0&0 \\0&-1&-0.5&1&0 \\0&-1&1&-0.6&1 \\\end{bmatrix},\quad \mathbf{U}=\begin{bmatrix}1&3&1&0&5 \\0&-1&1&0&-2 \\0&0&-2&2&3 \\0&0&0&5&-1.5 \\0&0&0&0&-7.9 \\\end{bmatrix}. \]

Ejemplo (continuación)

El sistema que tenemos que resolver es el siguiente: \[ \begin{bmatrix}1&3&1&0&5 \\1&2&2&0&3 \\1&-1&3&2&0 \\0&1&0&4&-1 \\0&1&-3&-1&-2 \\\end{bmatrix}\begin{bmatrix}x_1\\ x_2\\ x_3\\ x_4\\ x_5\end{bmatrix}=\begin{bmatrix}0&0&0&0&1 \\0&1&0&0&0 \\0&0&1&0&0 \\0&0&0&1&0 \\1&0&0&0&0 \\\end{bmatrix}\cdot \begin{bmatrix}2 \\2 \\2 \\3 \\1 \\\end{bmatrix}=\begin{bmatrix}1 \\2 \\2 \\3 \\2 \\\end{bmatrix}. \]

En primer lugar resolvemos \(\mathbf{L}\mathbf{y}=\mathbf{b}\): \[ \begin{bmatrix}1&0&0&0&0 \\1&1&0&0&0 \\1&4&1&0&0 \\0&-1&-0.5&1&0 \\0&-1&1&-0.6&1 \\\end{bmatrix}\cdot\begin{bmatrix}y_1\\ y_2\\ y_3\\ y_4\\ y_5\end{bmatrix}=\begin{bmatrix}1 \\2 \\2 \\3 \\2 \\\end{bmatrix}. \]

\[ \begin{align*} y_1 & =1,\ y_2=2- y_1=1,\ y_3=2-4 y_2- y_1=-3,\ y_4=3+0.5 y_3+ y_2=2.5,\\ y_5 & =2+0.6 y_4- y_3+ y_2=7.5. \end{align*} \]

Ejemplo (continuación)

En segundo lugar, resolvemos \(\mathbf{U}\mathbf{x}=\mathbf{y}\): \[ \begin{bmatrix}1&3&1&0&5 \\0&-1&1&0&-2 \\0&0&-2&2&3 \\0&0&0&5&-1.5 \\0&0&0&0&-7.9 \\\end{bmatrix}\cdot\begin{bmatrix}x_1\\ x_2\\ x_3\\ x_4\\ x_5\end{bmatrix}=\begin{bmatrix}1 \\1 \\-3 \\2.5 \\7.5 \\\end{bmatrix}. \]

La solución del sistema será usando el método de sustitución hacia atrás: \[ \begin{align*} x_5 & =\frac{7.5}{-7.9}=-0.949367,\\ x_4 & =\frac{2.5-(-1.5)\cdot (-0.9493671)}{5}=0.2151899,\\ x_3 & =\frac{-3-3\cdot (-0.9493671)-2\cdot 0.2151899}{-2}=0.2911392,\\ x_2 & =\frac{1-(-2)\cdot (-0.9493671)-0\cdot 0.2151899-1\cdot 0.2911392}{-1}=1.1898734,\\ x_1 & =\frac{1-5\cdot (-0.9493671)-0\cdot 0.2151899-1\cdot 0.2911392-3\cdot 1.1898734}{1}=1.8860759. \end{align*} \]

Ejemplo (continuación)

La mejora de la descomposición \(LU\) con matrices de permutación está implementada en:

Matrices diagonal dominantes

Veamos que si la matriz del sistema es estrictamente diagonal dominante, el algoritmo de eliminación de Gauss se puede realizar de forma óptima:

Teorema. Eliminación de Gauss para matrices diagonal dominantes. Sea \(\mathbf{A}\) una matriz \(n\times n\) estrictamente diagonal dominante. Entonces \(\mathbf{A}\) es no singular (\(\mathrm{det}(\mathbf{A})\neq 0\)).

Consideremos un sistema lineal de la forma \(\mathbf{A}\mathbf{x}=\mathbf{b}\) cuya matriz del sistema es \(\mathbf{A}\). Entonces el algoritmo de eliminación de Gauss se puede realizar sin permutar filas ni columnas y el método es estable con respecto a los errores de redondeo.

Demostración del Teorema

Veamos por reducción al absurdo que la matriz \(\mathbf{A}\) es no singular.

Consideremos el sistema lineal homogéneo \(\mathbf{A}\mathbf{x}=\mathbf{0}\) y supongamos que tenemos una solución \(\mathbf{x}=(x_i)_{i=1,\ldots,n}\neq\mathbf{0}\). Sea \(k\) el índice tal que: \[ 0 < |x_k|=\max_{i=1,\ldots,n}|x_i|. \] Como \(\displaystyle\sum_{j=1}^n a_{ij}x_j=0\), para todo \(i=1,\ldots,n\), tenemos que para \(i=k\), \[ \sum_{j=1}^n a_{kj}x_j=a_{kk}x_k+\sum_{j=1,j\neq k}^n a_{kj}x_j=0,\ \Rightarrow a_{kk}x_k=-\sum_{j=1,j\neq k}^n a_{kj}x_j. \] Usando la desigualdad triangular, tenemos: \[ |a_{kk}||x_k|\leq \sum_{j=1,j\neq k}^n |a_{kj}||x_j|,\ \Rightarrow |a_{kk}|\leq \sum_{j=1,j\neq k}^n |a_{kj}|\frac{|x_j|}{|x_k|}\leq \sum_{j=1,j\neq k}^n |a_{kj}|. \] La última desigualdad contradice que la matriz \(\mathbf{A}\) sea estrictamente diagonal dominante ya que dicha condición fallaría para la fila \(k\)-ésima. Por tanto, la matriz \(\mathbf{A}\) es no singular.

Demostración del Teorema (continuación)

Para ver que podemos realizar el algoritmo de eliminación de Gauss sin permutación de filas ni de columnas, demostraremos que las matrices obtenidas por el algoritmo \(\mathbf{A}^{(2)},\mathbf{A}^{(3)},\ldots,\mathbf{A}^{(n)}\) son todas estrictamente diagonal dominantes. Este hecho asegurará que el valor \(a_{kk}^{(k)}\neq 0\), para \(k=2,\ldots,n\) y por tanto, podemos continuar sin permutar filas ni columnas.

Como la matriz \(\mathbf{A}\) es estrictamente diagonal dominante, \(a_{11}^{(1)}\neq 0\) y podemos calcular la matriz \(\mathbf{A}^{(2)}\): \[ a_{ij}^{(2)}=a_{ij}^{(1)}-\frac{a_{1j}^{(1)}a_{i1}^{(1)}}{a_{11}^{(1)}}, \] para \(j=2,\ldots,n,\ i=2,\ldots,n\).

Demostración del Teorema (continuación)

Consideremos la fila \(i\) de la matriz \(\mathbf{A}^{(2)}\), \(i=2,\ldots,n\). Veamos que \(\displaystyle |a_{ii}^{(2)}|>\sum_{j=2,j\neq i}^n |a_{ij}^{(2)}|\) (pensemos que \(a_{i1}^{(2)}=0,\ i=2,\ldots,n\)) \[ \begin{align*} \sum_{j=2,j\neq i}^n |a_{ij}^{(2)}| & =\sum_{j=2,j\neq i}^n\left|a_{ij}^{(1)}-\frac{a_{1j}^{(1)}a_{i1}^{(1)}}{a_{11}^{(1)}}\right|\leq \sum_{j=2,j\neq i}^n |a_{ij}^{(1)}|+\sum_{j=2,j\neq i}^n\left|\frac{a_{1j}^{(1)}a_{i1}^{(1)}}{a_{11}^{(1)}}\right|\\ & < |a_{ii}^{(1)}|-|a_{i1}^{(1)}|+\frac{|a_{i1}^{(1)}|}{|a_{11}^{(1)}|}(|a_{11}^{(1)}|-|a_{1i}^{(1)}|). \end{align*} \] En la última desigualdad hemos usado que la matriz \(\mathbf{A}^{(1)}=\mathbf{A}\) es estrictamente diagonal dominante y, por tanto, \[ \begin{align*} \sum_{j=2,j\neq i}^n |a_{ij}^{(1)}| & =\sum_{j=1,j\neq i}^n |a_{ij}^{(1)}|-|a_{i1}^{(1)}|< |a_{ii}^{(1)}|-|a_{i1}^{(1)}|,\\ \sum_{j=2,j\neq i}^n |a_{1j}^{(1)}|& =\sum_{j=2}^n |a_{1j}^{(1)}|-|a_{1i}^{(1)}|< |a_{11}^{(1)}|-|a_{1i}^{(1)}|. \end{align*} \]

Demostración del Teorema (continuación)

Hemos demostrado que: \[ \sum_{j=2,j\neq i}^n |a_{ij}^{(2)}|< |a_{ii}^{(1)}|-|a_{i1}^{(1)}|+\frac{|a_{i1}^{(1)}|}{|a_{11}^{(1)}|}(|a_{11}^{(1)}|-|a_{1i}^{(1)}|)=|a_{ii}^{(1)}|-\frac{|a_{i1}^{(1)}|}{|a_{11}^{(1)}|}|a_{1i}^{(1)}|\leq \left|a_{ii}^{(1)}-\frac{a_{i1}^{(1)}}{a_{11}^{(1)}}a_{1i}^{(1)}\right|=|a_{ii}^{(2)}|, \] como queríamos demostrar. Por tanto, la matriz \(\mathbf{A}^{(2)}\) es estrictamente diagonal dominante.

Usando el mismo razonamiento anterior se puede demostrar por inducción que las matrices \(\mathbf{A}^{(k)}\), para \(k=2,\ldots,n\) que se van obteniendo en el algoritmo de eliminación de Gauss son diagonal dominantes y por tanto el proceso de Gauss se puede realizar sin permutar filas ni columnas.

Matrices simétricas y definidas positivas

Sea \(\mathbf{A}\) una matriz \(n\times n\) simétrica y estrictamente definida positiva.

Por el criterio de Sylvester sabemos que para todo \(k=1,\ldots,n\), \(\mathrm{det}(\mathbf{A}_k)>0\), donde \(\mathbf{A}_k\) recordemos que es la submatriz formada por las \(k\) primeras filas y columnas de la matriz \(\mathbf{A}\).

Usando una proposición anterior, el valor de \(\mathrm{det}(\mathbf{A}_k)>0\) era \(\mathrm{det}(\mathbf{A}_k)=a_{11}^{(1)}\cdots a_{kk}^{(k)}\) donde \(a_{ii}^{(i)}\) es el pivote del paso \(i\)-ésimo del algoritmo de eliminación de Gauss.

De las dos afirmaciones anteriores, se deduce que, como \(\mathrm{det}(\mathbf{A}_k)=a_{11}^{(1)}\cdots a_{kk}^{(k)} >0\), para \(k=1,\ldots,n\), los pivotes \(a_{kk}^{(k)}\) serán positivos \(a_{kk}^{(k)}>0\) y se puede aplicar por tanto, el algoritmo de eliminación de Gauss sin permutar filas ni columnas.

Matrices simétricas y definidas positivas

En el caso en que la matriz del sistema \(\mathbf{A}\) sea simétrica y definida positiva, tenemos el siguiente resultado:

Proposición: descomposición de una matriz simétrica y definida positiva. Sea \(\mathbf{A}\) una matriz del sistema lineal simétrica y estrictamente definida positiva. Entonces existen matrices \(\mathbf{L}\) triangular inferior con unos en la diagonal y \(\mathbf{D}\) matriz diagonal con todos sus elementos positivos tal que \(\mathbf{A}=\mathbf{L}\cdot\mathbf{D}\cdot\mathbf{L}^\top\).

Demostración de la proposición

Hemos visto anteriormente que la matriz \(\mathbf{A}\) admite descomposición \(LU\) ya que puede realizarse el algoritmo de eliminación de Gauss sin permutar filas ni columnas.

Entonces \(\mathbf{A}=\mathbf{L}\mathbf{U}\). Como \(\mathbf{A}^\top =\mathbf{A}\), al ser la matriz \(\mathbf{A}\) simétrica, tenemos que: \[ \mathbf{A}=\mathbf{A}^\top = \mathbf{U}^{\top}\cdot \mathbf{L}^\top. \] Ahora escribimos la matriz \(\mathbf{U}^\top\) como \(\mathbf{U}^\top =\mathbf{\tilde{U}}^\top\cdot\mathbf{D}\), donde \(\mathbf{\tilde{U}}^\top\) es una matriz triangular inferior con unos en la diagonal y \(\mathbf{D}\) es la matriz diagonal siguiente:

\[ \mathbf{D}=\begin{bmatrix}u_{11}&0&\ldots&0 \\0&u_{22}&\ldots&0 \\\vdots&\vdots&\ddots&\vdots \\0&\ldots&0&u_{nn} \\\end{bmatrix}. \] Como los elementos de la diagonal de \(\mathbf{D}\) son los elementos diagonales de la matriz \(\mathbf{U}\) y éstos, a su vez, son los pivotes en el algoritmo de eliminación de Gauss, serán positivos.

Entonces tenemos que \(\mathbf{A}=\mathbf{A}^\top =\mathbf{\tilde{U}}^\top\mathbf{D}\mathbf{L}^\top =\mathbf{L}\mathbf{U}.\)

Por la unicidad de la descomposición \(LU\) tenemos que: \(\mathbf{\tilde{U}}^\top =\mathbf{L},\quad \mathbf{D}\mathbf{L}^\top =\mathbf{U}\,\Rightarrow \mathbf{A}=\mathbf{L}\mathbf{D}\mathbf{L}^\top,\) tal como queríamos demostrar.

Factorización \(\mathbf{L}\mathbf{D}\mathbf{L}^\top\)

Sea \(\mathbf{A}\) una matriz simétrica y definida positiva. Vamos a dar el algoritmo para la descomposición \(\mathbf{A}=\mathbf{L}\mathbf{D}\mathbf{L}^\top\).

La idea es ir calculando las columnas de la matriz \(\mathbf{L}\) y los elementos diagonales de la matriz \(\mathbf{D}\).

Primer paso: \(d_{11}=a_{11}\) y a continuación primera columna de \(\mathbf{L}\): \[l_{j1}=\frac{a_{j1}}{d_{11}},\ j=2,\ldots,n.\]
Segundo paso: \(d_{22}=a_{22}-l_{21}^2 d_{11}\) y a continuación segunda columna de \(\mathbf{L}\): \[l_{j2}=\frac{1}{d_{22}}\left(a_{j2}-l_{j1}l_{21} d_{11}\right),\ j=3,\ldots,n.\]
\(\vdots\)

Factorización \(\mathbf{L}\mathbf{D}\mathbf{L}^\top\)

Paso \(i\)-ésimo: \(\displaystyle d_{ii}=a_{ii}-\sum_{j=1}^{i-1}l_{ij}^2 d_{jj}\) y a continuación la columna \(i\)-ésima de \(\mathbf{L}\): \[l_{ji}=\frac{1}{d_{ii}}\left(a_{ji}-\sum_{k=1}^{i-1}l_{jk}l_{ik} d_{kk}\right),\ j=i+1,\ldots,n.\]
\(\vdots\)
Hasta llegar al paso \(n\)-ésimo que calculamos \(\displaystyle d_{nn}=a_{nn}-\sum_{j=1}^{n-1}l_{nj}^2 d_{jj}\).

Factorización \(\mathbf{L}\mathbf{D}\mathbf{L}^\top\)

Observación. Para calcular la columna \(i\)-ésima de la matriz \(\mathbf{L}\) en el paso \(i\)-ésimo, \(l_{ji}\), \(j=i+1,\ldots,n\) se usan las \(i-1\) primeras columnas de la matriz \(\mathbf{L}\) que ya se suponen calculadas en pasos anteriores, fijarse que en el sumatorio para \(k\) salen los valores \(l_{jk}\), \(l_{ik}\), los cuales son elementos de la columna \(k\)-ésima, para \(k\) desde \(1\) hasta \(i-1\).

Factorización \(\mathbf{L}\mathbf{D}\mathbf{L}^\top\). Pseudocódigo

INPUT matriz del sistema \(\mathbf{A}=(a_{ij})_{i=1,\ldots,n,j=1,\ldots,n}\).
Set \(\mathbf{L}=\mathbf{Id}_n\).
Set \(d_{11}=\cdots = d_{nn}=0\).
For i=1,...n
- Set \(d_{ii}=a_{ii}-\sum_{j=1}^{i-1} l_{ij}^2 d_{jj}\).
- For j=i+1,...,n
- Set \(l_{ji}=\frac{1}{d_{ii}}\left(a_{ji}-\sum_{k=1}^{i-1}l_{jk}l_{ik} d_{kk}\right)\)
Print \(\mathbf{L}\), \(d_{11},\ldots,d_{nn}\) (Damos la matriz \(\mathbf{L}\) y los valores \(d_{ii}\))

Ejemplo

Consideremos la matriz simétrica y definida positiva siguiente: \[ \mathbf{A}=\begin{bmatrix}5&2&-2&-1&-1 \\2&7&1&1&6 \\-2&1&9&2&1 \\-1&1&2&11&-1 \\-1&6&1&-1&13 \\\end{bmatrix}. \] Vemos que la matriz anterior es simétrica. Comprobemos, aplicando el criterio de Sylvester que es estrictamente definida positiva: \[ \begin{align*} \mathrm{det}(\mathbf{A}_1) & =5 >0,\quad \mathrm{det}(\mathbf{A}_2) =\mathrm{det}\left(\begin{bmatrix}5&2 \\2&7 \\\end{bmatrix}\right)=31 >0,\\ \mathrm{det}(\mathbf{A}_3) & =\mathrm{det}\left(\begin{bmatrix}5&2&-2 \\2&7&1 \\-2&1&9 \\\end{bmatrix}\right)=238 >0. \end{align*} \]

Ejemplo (continuación)

\[ \begin{align*} \mathrm{det}(\mathbf{A}_4) & =\mathrm{det}\left(\begin{bmatrix}5&2&-2&-1 \\2&7&1&1 \\-2&1&9&2 \\-1&1&2&11 \\\end{bmatrix}\right)=2451 >0,\\ \mathrm{det}(\mathbf{A}_5) & =\mathrm{det}\left(\begin{bmatrix}5&2&-2&-1&-1 \\2&7&1&1&6 \\-2&1&9&2&1 \\-1&1&2&11&-1 \\-1&6&1&-1&13 \\\end{bmatrix}\right)=13247 >0. \end{align*} \]

A continuación apliquemos el algoritmo anterior:

Paso 1: \(d_{11}=a_{11}=5\). Primer columna de \(\mathbf{L}\): \[ \begin{align*} l_{21} & =\frac{a_{21}}{d_{11}}=\frac{2}{5}=0.4,\ l_{31}=\frac{a_{31}}{d_{11}}=\frac{-2}{5}=-0.4,\ l_{41}=\frac{a_{41}}{d_{11}}=\frac{-1}{5}=-0.2,\\ l_{51} & =\frac{a_{51}}{d_{11}}=\frac{-1}{5}=-0.2. \end{align*} \]

Ejemplo (continuación)

Paso 2: \(d_{22}=a_{22}-l_{21}^2 d_{11}=7-0.4^2\cdot 5=6.2.\) Segunda columna de \(L\): \[ \begin{align*} l_{32} & =\frac{1}{d_{22}}\left(a_{32}-l_{31}l_{21}d_{11}\right)=\frac{1}{6.2}(1-(-0.4)\cdot 0.4\cdot 5)=0.2903226,\\ l_{42} & =\frac{1}{d_{22}}\left(a_{42}-l_{41}l_{21}d_{11}\right)=\frac{1}{6.2}(1-(-0.2)\cdot 0.4\cdot 5)=0.2258065,\\ l_{52} & =\frac{1}{d_{22}}\left(a_{52}-l_{51}l_{21}d_{11}\right)=\frac{1}{6.2}(6-(-0.2)\cdot 0.4\cdot 5)=1.0322581. \end{align*} \]
Paso 3: \(d_{33}=a_{33}-l_{31}^2 d_{11}-l_{32}^2 d_{22}=9-(-0.4)^2\cdot 5-0.2903226^2\cdot 6.2=7.6774194.\) Tercera columna de \(L\): \[ \begin{align*} l_{43} & =\frac{1}{d_{33}}\left(a_{43}-l_{41}l_{31}d_{11}-l_{42}l_{32}d_{22}\right)\\ & =\frac{1}{7.6774194}(2-(-0.2)\cdot (-0.4)\cdot 5-0.2258065\cdot 0.2903226\cdot 6.2)=0.1554622,\\ l_{53} & =\frac{1}{d_{33}}\left(a_{53}-l_{51}l_{31}d_{11}-l_{52}l_{32}d_{22}\right)\\ & =\frac{1}{7.6774194}(1-(-0.2)\cdot (-0.4)\cdot 5-1.0322581\cdot 0.2903226\cdot 6.2)=-0.1638655. \end{align*} \]

Ejemplo (continuación)

Paso 4: \[ \begin{align*} d_{44} & =a_{44}-l_{41}^2 d_{11}-l_{42}^2 d_{22}-l_{43}^2 d_{33}\\ & =11-(-0.2)^2\cdot 5-0.2258065^2\cdot 6.2-0.1554622^2\cdot 7.6774194=10.2983193. \end{align*} \] Cuarta columna de \(L\): \[ \begin{align*} l_{54} & =\frac{1}{d_{44}}\left(a_{54}-l_{51}l_{41}d_{11}-l_{52}l_{42}d_{22}-l_{53}l_{43}d_{33}\right)\\ & =\frac{1}{10.2983193}(-1-(-0.2)\cdot (-0.2)\cdot 5-1.0322581\cdot 0.2258065\cdot 6.2\\ & \quad -(-0.1638655)\cdot 0.1554622\cdot 7.6774194)=-0.2378621. \end{align*} \]
Paso 5 y último: \[ \begin{align*} d_{55} & =a_{55}-l_{51}^2 d_{11}-l_{52}^2 d_{22}-l_{53}^2 d_{33}-l_{54}^2 d_{44}\\ & =13-(-0.2)^2\cdot 5-1.0322581^2\cdot 6.2\\ &\quad -(-0.1638655)^2\cdot 7.6774194-(-0.2378621)^2\cdot 10.2983193 =5.4047328. \end{align*} \]

Ejemplo (continuación)

Las matrices \(\mathbf{L}\) y \(\mathbf{D}\) serán: \[ \begin{align*} \mathbf{L} & =\begin{bmatrix}1&0&0&0&0 \\0.4&1&0&0&0 \\-0.4&0.290323&1&0&0 \\-0.2&0.225806&0.155462&1&0 \\-0.2&1.032258&-0.163866&-0.237862&1 \\\end{bmatrix},\\ \mathbf{D} & =\begin{bmatrix}5&0&0&0&0 \\0&6.2&0&0&0 \\0&0&7.677419&0&0 \\0&0&0&10.298319&0 \\0&0&0&0&5.404733 \\\end{bmatrix}. \end{align*} \]

Número de operaciones

Para resolver un sistema lineal \(\mathbf{A}\mathbf{x}=\mathbf{b}\) donde la matriz \(\mathbf{A}\) es simétrica y estrictamente definida positiva necesitamos realizar los pasos siguentes:

Paso 1: Descomponer la matriz \(\mathbf{A}\) usando la descomposición: \(\mathbf{A}=\mathbf{L}\mathbf{D}\mathbf{L}^\top\). El sistema a resolver quedará: \(\mathbf{L}\mathbf{D}\mathbf{L}^\top\mathbf{x}=\mathbf{b}\).
Paso 2: Resolver el sistema \(\mathbf{L}\mathbf{y}=\mathbf{b}\).
Paso 3: Resolver el sistema \(\mathbf{D}\mathbf{z}=\mathbf{y}\).
Paso 4: Resolver el sistema \(\mathbf{L}^\top\mathbf{x}=\mathbf{z}\).

Los pasos 2 y 4 consisten en resolver un sistema lineal triangular inferior y triangular superior. Vimos en su momento que resolver un sistema triangular superior requería \(n^2-n\) operaciones básicas.

Se puede ver de la misma manera que resolver un sistema lineal triangular inferior requiere también \(n^2-n\) operaciones.

Número de operaciones

El número de operaciones básicas para resolver un sistema diagonal (paso 3) vale \(n\) ya que el sistema a resolver es el siguiente:

\[ \mathbf{D}\mathbf{z}=\mathbf{y},\Rightarrow \begin{bmatrix}d_{11}&0&\ldots&0 \\0&d_{22}&\ldots&0 \\\vdots&\vdots&\ddots&\vdots \\0&\ldots&0&d_{nn} \\\end{bmatrix}\begin{bmatrix}z_1\\\vdots\\ z_n\end{bmatrix}=\begin{bmatrix}y_1\\\vdots\\ y_n\end{bmatrix}. \] La solución del mismo es: \[ z_1=\frac{y_1}{d_{11}},\ldots, z_n=\frac{y_n}{d_{nn}}. \] En total, hemos realizado \(n\) divisiones, lo que equivale a \(n\) operaciones básicas.

Número de operaciones. Paso 1

En el paso 1 o en la factorización \(\mathbf{L}\mathbf{D}\mathbf{L}^\top\) realizamos las operaciones siguientes en el paso \(i\)-ésimo:

Cálculo de \(d_{ii}\): \(i\) sumas/restas y \(2(i-1)\) multiplicaciones (en la suma de \(j\) desde \(1\) hasta \(i-1\) hay que calcular \(l_{ij}^2\) y seguidamente multiplicar por \(d_{jj}\)). En total, \(2(i-1)+i=3i-2\) operaciones básicas.
Cálculo de \(l_{ji}\), para \(j=i+1,\ldots,n\):
- \(i\) sumas/restas y \(2(i-1)+1\) multiplicaciones divisiones. En total, \(3i-1\) operaciones básicas.
- En total, para calcular todos los \(l_{ji}\), tenemos que realizar \(\displaystyle\sum_{j=i+1}^n (3i-1)=(3i-1)(n-i)\) operaciones básicas.

Número de operaciones. Paso 1

En el paso \(i\)-ésimo, tenemos que realizar un total de \(3i-2+(3i-1)(n-i)\) operaciones básicas.

Como el índice \(i\) va desde \(1\) hasta \(n\), el número total de operaciones básicas a realizar en el paso 1 será: \[ \sum_{i=1}^n (3i-2+(3i-1)(n-i))=\frac{1}{2} n (n (n+2)-1)=\frac{n^3}{2}+n^2-\frac{n}{2}. \]

Número de operaciones

En resumen, el número total de operaciones para resolver un sistema lineal \(\mathbf{A}\mathbf{x}=\mathbf{b}\), donde la matriz del sistema \(\mathbf{A}\) es simétrica y definida positiva usando la descomposición \(\mathbf{A}=\mathbf{L}\mathbf{D}\mathbf{L}^\top\) viene dado en la tabla siguiente:

Pasos	Número de operaciones
\(1\)	\(\frac{n^3}{2}+n^2-\frac{n}{2}\)
\(2\)	\(n^2-n\)
\(3\)	\(n\)
\(4\)	\(n^2-n\)
Total	\(\frac{n^3}{2}+3n^2-\frac{3n}{2}\)

Descomposición de Choleski

Hemos visto que si una matriz \(\mathbf{A}\) es simétrica y definida positiva, se puede descomponer como: \(\mathbf{A}=\mathbf{L}\mathbf{D}\mathbf{L}^\top\), donde \(\mathbf{L}\) es una matriz triangular inferior con unos en la diagonal y \(\mathbf{D}\) es una matriz diagonal con todas sus componentes diagonales positivas.

Usando la descomposición anterior, podemos escribir: \[ \mathbf{A}=\mathbf{L}\mathbf{D}\mathbf{L}^\top=\mathbf{L}\mathbf{D}^{\frac{1}{2}}\mathbf{D}^{\frac{1}{2}}\mathbf{L}^\top, \] donde

\[ \mathbf{D}^{\frac{1}{2}}=\begin{bmatrix}\sqrt{d_{11}}&0&\ldots&0 \\0&\sqrt{d_{22}}&\ldots&0 \\\vdots&\vdots&\ddots&\vdots \\0&\ldots&0&\sqrt{d_{nn}} \\\end{bmatrix}, \]

Descomposición de Choleski

suponiendo que la matriz diagonal \(\mathbf{D}\) es la siguiente: \[ \mathbf{D}=\begin{bmatrix}d_{11}&0&\ldots&0 \\0&d_{22}&\ldots&0 \\\vdots&\vdots&\ddots&\vdots \\0&\ldots&0&d_{nn} \\\end{bmatrix}. \]

Si llamamos \(\tilde{\mathbf{L}}=\mathbf{L}\mathbf{D}^{\frac{1}{2}}\), tenemos que la matriz \(\mathbf{A}\) puede descomponerse como \(\mathbf{A}=\tilde{\mathbf{L}}\tilde{\mathbf{L}}^\top\). A dicha descomposición se le llama descomposición de Choleski.

Descomposición de Choleski

Para realizar la descomposición de Choleski usamos un procedimiento similar cuando realizábamos la descomposición \(\mathbf{L}\mathbf{D}\mathbf{L}^\top\).

Es decir, vamos calculando las columnas de \(\tilde{\mathbf{L}}\):

Primer paso: \(\tilde{l}_{11}=\sqrt{a_{11}}\) y a continuación la primera columna de \(\tilde{\mathbf{L}}\): \(\tilde{l}_{j1}=\frac{a_{j1}}{\tilde{l}_{11}}\).
Segundo paso: \(\tilde{l}_{22}=\sqrt{a_{22}-\tilde{l}_{21}}\) y a continuación la segunda columna de \(\tilde{\mathbf{L}}\): \(\tilde{l}_{j2}=\frac{a_{j2}-\tilde{l}_{j1}\tilde{l}_{21}}{\tilde{l}_{22}}\), para \(j=3,\ldots,n\)
\(\vdots\)

Descomposición de Choleski

Paso \(i\)-ésimo: \(\tilde{l}_{ii}=\left(a_{ii}-\sum_{k=1}^{i-1}\tilde{l}_{ik}^2\right)^{\frac{1}{2}}\) y a continuación la columna \(i\)-ésima de \(\tilde{\mathbf{L}}\): \(\tilde{l}_{ji}=\frac{\left(a_{ji}-\sum_{k=1}^{i-1}\tilde{l}_{jk}\tilde{l}_{ik}\right)}{\tilde{l}_{ii}}\), para \(j=i+1,\ldots,n\).

*\(\vdots\)

Paso \(n\)-ésimo: \(\tilde{l}_{nn}=\left(a_{nn}-\sum_{k=1}^{n-1}\tilde{l}_{nk}^2\right)^{\frac{1}{2}}\).

Descomposición de Choleski. Pseudocódigo

INPUT matriz del sistema \(\mathbf{A}=(a_{ij})_{i=1,\ldots,n,j=1,\ldots,n}\).
Set \(\tilde{\mathbf{L}}=\mathbf{0}\).
Set \(\tilde{l}_{11}=\sqrt{a_{11}}\).
For i=2,...n
- Set \(\tilde{l}_{j1}=\frac{a_{j1}}{\tilde{l}_{11}}\).

Descomposición de Choleski. Pseudocódigo

For i=2,...n
- Set \(\tilde{l}_{ii}=\left(a_{ii}-\sum_{k=1}^{i-1}\tilde{l}_{ik}^2\right)^{\frac{1}{2}}\).
- For j=i+1,...,n
- Set \(\tilde{l}_{ji}=\frac{\left(a_{ji}-\sum_{k=1}^{i-1}\tilde{l}_{jk}\tilde{l}_{ik}\right)}{\tilde{l}_{ii}}\).
Set \(\tilde{l}_{nn}=\left(a_{nn}-\sum_{k=1}^{n-1}\tilde{l}_{nk}^2\right)^{\frac{1}{2}}\).
Print \(\tilde{\mathbf{L}}\) (Damos la matriz \(\tilde{\mathbf{L}}\)).

Ejemplo

Vamos a hallar la descomposición de Choleski de la matriz siguiente vista anteriormente: \[ \mathbf{A}=\begin{bmatrix}5&2&-2&-1&-1 \\2&7&1&1&6 \\-2&1&9&2&1 \\-1&1&2&11&-1 \\-1&6&1&-1&13 \\\end{bmatrix}. \]

Paso 1: \(\tilde{l}_{11}=\sqrt{a_{11}}=\sqrt{5}=2.236068\) y a continuación primera columna de \(\tilde{\mathbf{L}}\): \[ \begin{align*} \tilde{l}_{21} & =\frac{a_{21}}{\tilde{l}_{11}}=\frac{2}{2.236068}=0.894427,\\ \tilde{l}_{31} & =\frac{a_{31}}{\tilde{l}_{11}}=\frac{-2}{2.236068}=-0.894427,\\ \tilde{l}_{41} & =\frac{a_{41}}{\tilde{l}_{11}}=\frac{-1}{2.236068}=-0.447214,\\ \tilde{l}_{51} & =\frac{a_{51}}{\tilde{l}_{11}}=\frac{-1}{2.236068}=-0.447214. \end{align*} \]

Ejemplo (continuación)

\[ \tilde{\mathbf{L}}=\begin{bmatrix}2.236068&0&0&0&0 \\0.894427&*&0&0&0 \\-0.894427&*&*&0&0 \\-0.447214&*&*&*&0 \\-0.447214&*&*&*&* \\\end{bmatrix}. \]

Paso 2: \(\tilde{l}_{22}=\sqrt{a_{22}-\tilde{l}_{21}^2}=\sqrt{7-0.894427^2}=2.4899799\) y a continuación segunda columna de \(\tilde{\mathbf{L}}\): \[ \begin{align*} \tilde{l}_{32} & =\frac{a_{32}-\tilde{l}_{31}\tilde{l}_{21}}{\tilde{l}_{22}}=\frac{1-(-0.894427)\cdot 0.894427}{2.48998}=0.722897,\\ \tilde{l}_{42} & =\frac{a_{42}-\tilde{l}_{41}\tilde{l}_{21}}{\tilde{l}_{22}}=\frac{1-(-0.447214)\cdot 0.894427}{2.48998}=0.562254,\\ \tilde{l}_{52} & =\frac{a_{52}-\tilde{l}_{51}\tilde{l}_{21}}{\tilde{l}_{22}}=\frac{6-(-0.447214)\cdot 0.894427}{2.48998}=2.570302. \end{align*} \]

Ejemplo (continuación)

\[ \tilde{\mathbf{L}}=\begin{bmatrix}2.236068&0&0&0&0 \\0.894427&2.48998&0&0&0 \\-0.894427&0.722897&*&0&0 \\-0.447214&0.562254&*&*&0 \\-0.447214&2.570302&*&*&* \\\end{bmatrix}. \]

Paso 3: \(\tilde{l}_{33}=\sqrt{a_{33}-\tilde{l}_{31}^2-\tilde{l}_{32}^2}=\sqrt{9-(-0.894427)^2-0.722897^2}=2.7708156\) y a continuación tercera columna de \(\tilde{\mathbf{L}}\): \[ \begin{align*} \tilde{l}_{43} & =\frac{a_{43}-\tilde{l}_{41}\tilde{l}_{31}-\tilde{l}_{42}\tilde{l}_{32}}{\tilde{l}_{33}}=\frac{2-(-0.447214)\cdot (-0.894427)-0.562254\cdot 0.722897}{2.770816}\\ & =0.430757,\\ \tilde{l}_{53} & =\frac{a_{53}-\tilde{l}_{51}\tilde{l}_{31}-\tilde{l}_{52}\tilde{l}_{32}}{\tilde{l}_{33}}=\frac{1-(-0.447214)\cdot (-0.894427)-2.570302\cdot 0.722897}{2.770816}\\ & =-0.454041. \end{align*} \]

Ejemplo (continuación)

\[ \tilde{\mathbf{L}}=\begin{bmatrix}2.236068&0&0&0&0 \\0.894427&2.48998&0&0&0 \\-0.894427&0.722897&2.770816&0&0 \\-0.447214&0.562254&0.430757&*&0 \\-0.447214&2.570302&-0.454041&*&* \\\end{bmatrix}. \]

Paso 4: \(\tilde{l}_{44}=\sqrt{a_{44}-\tilde{l}_{41}^2-\tilde{l}_{42}^2-\tilde{l}_{43}^2}=\sqrt{11-(-0.447214)^2-0.562254^2-0.430757^2}=3.2090995\) y a continuación cuarta columna de \(\tilde{\mathbf{L}}\): \[ \begin{align*} \tilde{l}_{54} & =\frac{a_{54}-\tilde{l}_{51}\tilde{l}_{41}-\tilde{l}_{52}\tilde{l}_{42}-\tilde{l}_{53}\tilde{l}_{43}}{\tilde{l}_{44}}\\ & =\frac{-1-(-0.447214)\cdot (-0.447214)-2.570302\cdot 0.562254-(-0.454041)\cdot 0.430757}{3.209099} \\ & =-0.763323. \end{align*} \]

Ejemplo (continuación)

\[ \tilde{\mathbf{L}}=\begin{bmatrix}2.236068&0&0&0&0 \\0.894427&2.48998&0&0&0 \\-0.894427&0.722897&2.770816&0&0 \\-0.447214&0.562254&0.430757&3.209099&0 \\-0.447214&2.570302&-0.454041&-0.763323&* \\\end{bmatrix}. \]

Paso 5 y último: \[ \begin{align*} \tilde{l}_{55} & =\sqrt{a_{55}-\tilde{l}_{51}^2-\tilde{l}_{52}^2-\tilde{l}_{53}^2-\tilde{l}_{54}^2}\\ & =\sqrt{13-(-0.447214)^2-2.570302^2-(-0.454041)^2-(-0.763323)^2}\\ & =2.3248081. \end{align*} \]

Ejemplo (continuación)

Veamos el ejemplo implementado:

Matrices tridiagonales. Factorización de Crout

En muchos problemas de ingeniería y de machine learning, nos aparecen sistemas lineales cuya matriz del sistema es tridiagonal:

\[\mathbf{A}=\begin{bmatrix}a_{11}&a_{12}&0&0&\ldots&0 \\a_{21}&a_{22}&a_{23}&0&\ldots&0 \\0&a_{32}&a_{33}&a_{34}&\ldots&0 \\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots \\0&\ldots&0&a_{n-1,n-2}&a_{n-1,n-1}&a_{n-1,n} \\0&\ldots&\ldots&0&a_{n,n-1}&a_{n,n} \\\end{bmatrix}.\]

Matrices tridiagonales. Factorización de Crout

En este caso, vamos a factorizar la matriz \(\mathbf{A}\) de la forma \(\mathbf{A}=\mathbf{L}\mathbf{U}\) donde las matrices \(\mathbf{L}\) y \(\mathbf{U}\) tienen la forma siguiente:

\[ \mathbf{L}=\begin{bmatrix}l_{11}&0&0&0&\ldots&0 \\l_{21}&l_{22}&0&0&\ldots&0 \\0&l_{32}&l_{33}&0&\ldots&0 \\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots \\0&\ldots&0&l_{n-1,n-2}&a_{n-1,n-1}&0 \\0&\ldots&\ldots&0&l_{n,n-1}&l_{n,n} \\\end{bmatrix},\ \mathbf{U}=\begin{bmatrix}1&u_{12}&0&0&\ldots&0 \\0&1&u_{23}&0&\ldots&0 \\0&0&1&u_{34}&\ldots&0 \\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots \\0&\ldots&0&0&1&u_{n-1,n} \\0&\ldots&\ldots&0&0&1 \\\end{bmatrix}. \]

Matrices tridiagonales. Factorización de Crout

Usando que \(\mathbf{A}=\mathbf{L}\mathbf{U}\), podemos hallar las componentes de las matrices \(\mathbf{L}\) y \(\mathbf{U}\).

Consideremos la fila \(i\)-ésima de la matriz \(\mathbf{A}\): \[ (0,\ldots,0,a_{i,i-1},a_{i,i},a_{i,i+1},0,\ldots,0), \]

Matrices tridiagonales. Factorización de Crout

entonces:

el valor \(a_{i,i-1}\) se obtiene multiplicando la fila \(i\)-ésima de la matriz \(\mathbf{L}\)

\[(0,\ldots,0,\overbrace{l_{i,i-1}}^{i-1},\overbrace{l_{ii}}^{i},0,\ldots,0)\]

por la columna \(i-1\)-ésima de la matriz \(\mathbf{U}\) \[ (0,\ldots,0,\overbrace{u_{i-1,i-2}}^{i-2},\overbrace{1}^{i-1},0,\ldots,0), \]

donde hemos indicado entre llaves las componentes de cada valor: \[ a_{i,i-1}=0\cdot u_{i-1,i-2}+l_{i,i-1}\cdot 1+l_{ii}\cdot 0=l_{i,i-1}. \]

Matrices tridiagonales. Factorización de Crout

el valor \(a_{i,i}\) se obtiene multiplicando la fila \(i\)-ésima de la matriz \(\mathbf{L}\) \[ (0,\ldots,0,\overbrace{l_{i,i-1}}^{i-1},\overbrace{l_{ii}}^{i},0,\ldots,0) \] por la columna \(i\)-ésima de la matriz \(\mathbf{U}\) \[ (0,\ldots,0,\overbrace{u_{i-1,i}}^{i-1},\overbrace{1}^{i},0,\ldots,0): \] \[ a_{i,i}=l_{i,i-1}\cdot u_{i-1,i}+l_{ii}\cdot 1=l_{i,i-1}\cdot u_{i-1,i}+l_{ii}. \]

Matrices tridiagonales. Factorización de Crout

el valor \(a_{i,i+1}\) se obtiene multiplicando la fila \(i\)-ésima de la matriz \(\mathbf{L}\) \[ (0,\ldots,0,\overbrace{l_{i,i-1}}^{i-1},\overbrace{l_{ii}}^{i},0,\ldots,0) \] por la columna \(i+1\)-ésima de la matriz \(\mathbf{U}\) \[ (0,\ldots,0,\overbrace{u_{i,i+1}}^{i},\overbrace{1}^{i+1},0,\ldots,0): \] \[ a_{i,i}=l_{i,i-1}\cdot 0+l_{i,i}\cdot u_{i,i+1}+0\cdot 1=l_{i,i}\cdot u_{i,i+1}. \]

Matrices tridiagonales. Factorización de Crout

A partir del primero de los tres resultados anteriores, \(a_{i,i-1}=l_{i,i-1},\) podemos calcular los valores \(l_{i,i-1}\) que van por debajo de la diagonal de la matriz \(\mathbf{L}\).

Para \(i=1\) tenemos que \(l_{11}=a_{11}\). Entonces usando los dos resultados últimos \[ \begin{align*} a_{i,i} & =l_{i,i-1}\cdot u_{i-1,i}+l_{i,i},\\ a_{i,i+1} &=l_{i,i}\cdot u_{i,i+1}, \end{align*} \] vamos calculando los valores \(u_{i,i+1}\) y \(l_{ii}\) alternativamente.

Matrices tridiagonales. Factorización de Crout

Es decir,

\(i=1\): \(l_{11}=a_{11}\), \[ a_{12}=l_{11}\cdot u_{12},\ \Rightarrow u_{12}=\frac{a_{12}}{l_{11}}. \]
\(i=2\): \[ \begin{align*} a_{22} & =l_{21}\cdot u_{12}+l_{22},\ \Rightarrow l_{22}=a_{22}-l_{21}\cdot u_{12},\\ a_{23} & =l_{22}\cdot u_{23},\ \Rightarrow u_{23}=\frac{a_{23}}{l_{22}}, \end{align*} \]
y así sucesivamente hasta llegar a \(i=n\).

Factorización de Crout. Pseudocódigo

INPUT matriz del sistema \(\mathbf{A}=(a_{ij})_{i=1,\ldots,n,j=1,\ldots,n}\).
Set \(\mathbf{L}=\mathbf{0}\).
Set \(\mathbf{U}=\mathbf{Id}_n\).
Set \(l_{11}=a_{11}\).
Set \(u_{12}=\frac{a_{12}}{l_{11}}\).
For i=2,...n-1
- Set \(l_{i,i-1}=a_{i,i-1}\). (fila \(i\)-ésima de \(\mathbf{L}\))
- Set \(l_{ii}=a_{ii}-l_{i,i-1}u_{i-1,i}\).
- Set \(u_{i,i+1}=\frac{a_{i,i+1}}{l_{i,i}}\). (columna \(i+1\)-ésima de \(\mathbf{U}\))

Factorización de Crout. Pseudocódigo

Set \(l_{n,n-1}=a_{n,n-1}\).
Set \(l_{n,n}=a_{n,n}-l_{n,n-1}u_{n-1,n}\)
Print \(\mathbf{L}\), \(\mathbf{U}\).

Ejemplo

Nos piden resolver el sistema de ecuaciones siguiente: \[ \begin{align*} E_1: & 5x_1 & +8 x_2 & & & & =& 1,\\ E_2: & 6x_1 & +6 x_2 & +2 x_3 & & & = & -1,\\ E_3: & & 6x_2 & +6 x_3 & +7 x_4 & & = & 4,\\ E_4: &&&+4 x_3&+5 x_4&+5 x_5&=&2,\\ E_5: & & & &6x_4 & +6 x_5 & = & 3. \end{align*} \] La matrix del sistema es tridiagonal: \[ \mathbf{A}=\begin{bmatrix}5&8&0&0&0 \\6&6&2&0&0 \\0&6&6&7&0 \\0&0&4&5&5 \\0&0&0&6&6 \\\end{bmatrix}. \]

Ejemplo (continuación)

Vamos a usar la descomposición de Crout para resolver el sistema:

Paso 1: \(l_{11}=a_{11}=5\), \(u_{12}=\frac{a_{12}}{l_{11}}=\frac{8}{5}=1.6.\)

\[ \mathbf{L}=\begin{bmatrix}5&0&0&0&0 \\*&*&0&0&0 \\0&*&*&0&0 \\0&0&*&*&0 \\0&0&0&*&* \\\end{bmatrix},\quad \mathbf{U}=\begin{bmatrix}1&1.6&0&0&0 \\0&1&*&0&0 \\0&0&1&*&0 \\0&0&0&1&* \\0&0&0&0&1 \\\end{bmatrix}. \]

Ejemplo (continuación)

Paso 2: \[ \begin{align*} l_{21} & =a_{21}=6,\\ l_{22} & =a_{22}-l_{21}u_{12}=6-6\cdot 1.6=-3.6,\\ u_{23} & =\frac{a_{23}}{l_{22}}=\frac{2}{-3.6}=-0.5555556. \end{align*} \]

\[ \mathbf{L}=\begin{bmatrix}5&0&0&0&0 \\6&-3.6&0&0&0 \\0&*&*&0&0 \\0&0&*&*&0 \\0&0&0&*&* \\\end{bmatrix},\quad \mathbf{U}=\begin{bmatrix}1&1.6&0&0&0 \\0&1&-0.555556&0&0 \\0&0&1&*&0 \\0&0&0&1&* \\0&0&0&0&1 \\\end{bmatrix}. \]

Ejemplo (continuación)

Paso 3: \[ \begin{align*} l_{32} & =a_{32}=6,\\ l_{33} & =a_{33}-l_{32}u_{23}=6-6\cdot (-0.5555556)=9.3333333,\\ u_{34} & =\frac{a_{34}}{l_{33}}=\frac{7}{9.3333333}=0.75. \end{align*} \]

\[ \mathbf{L}=\begin{bmatrix}5&0&0&0&0 \\6&-3.6&0&0&0 \\0&6&9.333333&0&0 \\0&0&*&*&0 \\0&0&0&*&* \\\end{bmatrix},\quad \mathbf{U}=\begin{bmatrix}1&1.6&0&0&0 \\0&1&-0.555556&0&0 \\0&0&1&0.75&0 \\0&0&0&1&* \\0&0&0&0&1 \\\end{bmatrix}. \]

Ejemplo (continuación)

Paso 4: \[ \begin{align*} l_{43} & =a_{43}=4,\\ l_{44} & =a_{44}-l_{43}u_{34}=5-4\cdot (0.75)=2,\\ u_{45} & =\frac{a_{45}}{l_{44}}=\frac{5}{2}=2.5. \end{align*} \]

\[ \mathbf{L}=\begin{bmatrix}5&0&0&0&0 \\6&-3.6&0&0&0 \\0&6&9.333333&0&0 \\0&0&4&2&0 \\0&0&0&*&* \\\end{bmatrix},\quad \mathbf{U}=\begin{bmatrix}1&1.6&0&0&0 \\0&1&-0.555556&0&0 \\0&0&1&0.75&0 \\0&0&0&1&2.5 \\0&0&0&0&1 \\\end{bmatrix}. \]

Ejemplo (continuación)

Paso 5 y último:

\[ \begin{align*} l_{54} & =a_{54}=6,\\ l_{55} & =a_{55}-l_{54}u_{45}=6-6\cdot (2.5)=-9. \end{align*} \]

\[ \mathbf{L}=\begin{bmatrix}5&0&0&0&0 \\6&-3.6&0&0&0 \\0&6&9.333333&0&0 \\0&0&4&2&0 \\0&0&0&6&-9 \\\end{bmatrix},\quad \mathbf{U}=\begin{bmatrix}1&1.6&0&0&0 \\0&1&-0.555556&0&0 \\0&0&1&0.75&0 \\0&0&0&1&2.5 \\0&0&0&0&1 \\\end{bmatrix}. \]

Ejemplo (continuación)

A continuación, vamos a resolver el sistema.

Como hemos descompuesto \(\mathbf{A}\) de la forma \(\mathbf{A}=\mathbf{L}\mathbf{U}\), resolver el sistema \(\mathbf{A}\mathbf{x}=\mathbf{L}\mathbf{U}\mathbf{x}=\mathbf{b}\) con:

\[ \mathbf{A}=\begin{bmatrix}5&8&0&0&0 \\6&6&2&0&0 \\0&6&6&7&0 \\0&0&4&5&5 \\0&0&0&6&6 \\\end{bmatrix},\quad \mathbf{b}=\begin{bmatrix}1 \\-1 \\4 \\2 \\3 \\\end{bmatrix}, \]

recordemos que primero tenemos que resolver el sistema \(\mathbf{L}\mathbf{y}=\mathbf{b}\), y, una vez hallado \(\mathbf{y}\), resolver el sistema \(\mathbf{U}\mathbf{x}=\mathbf{y}\) y hallar la solución \(\mathbf{x}\).

Ejemplo (continuación)

Resolución de \(\mathbf{L}\mathbf{y}=\mathbf{b}\):

\[ \begin{align*} E_1: & 5y_1 & & & & & = & 1,\\ E_2: & 6y_1 & -3.6 y_2 & & & & = & -1,\\ E_3: & & 6y_2 & +9.333333 y_3 & & & = & 4,\\ E_4: &&&+4 y_3&+2 y_4&&=&2,\\ E_5: & & & &6y_4 & -9 y_5 & = & 3. \end{align*} \]

Usamos el método de sustitución hacia adelante:

\[ \begin{align*} y_1 & =\frac{1}{5}=0.2,\\ y_2 & =\frac{-1-6\cdot 0.2}{-3.6}=0.6111111,\\ y_3 & =\frac{4-6\cdot 0.6111111}{9.3333333}=0.0357143,\\ y_4 & =\frac{2-4\cdot 0.0357143}{2}=0.9285714,\\ y_5 & =\frac{3-6\cdot 0.9285714}{-9}=0.2857143. \end{align*} \]

Ejemplo (continuación)

Resolución de \(\mathbf{U}\mathbf{x}=\mathbf{y}\):

\[ \begin{align*} E_1: & x_1 & +1.6 x_2 & & & & =& 0.2,\\ E_2: & & x_2 & -0.555556 x_3 & & & = & 0.6111111,\\ E_3: & & & x_3 & +0.75 x_4 & & = & 0.0357143,\\ E_4: & & & & x_4 & +2.5 x_5&=&0.9285714,\\ E_5: & & & & & x_5 & = & 0.2857143. \end{align*} \]

Usamos el método de sustitución hacia atrás:

\[ \begin{align*} x_5 & =0.2857143,\\ x_4 &=0.9285714-2.5\cdot 0.2857143=0.2142857,\\ x_3 & =0.0357143-0.75\cdot 0.2142857=-0.125,\\ x_2 & =0.6111111-(-0.5555556)\cdot (-0.125)=0.5416667,\\ x_1 & =0.2-1.6\cdot 0.5416667=-0.6666667. \end{align*} \]

Métodos de ortogonalización

Introducción

Los métodos de ortogonalización para resolver un sistema lineal \(\mathbf{A}\mathbf{x}=\mathbf{b}\) consisten en descomponer la matriz del sistema \(\mathbf{A}\), regular \(n\times n\), de la forma \(\mathbf{A}=\mathbf{Q}\mathbf{R}\), donde la matriz \(\mathbf{Q}\) es ortogonal, (\(\mathbf{Q}^\top =\mathbf{Q}^{-1}\)) y la matriz \(\mathbf{R}\) es triangular superior.

Una vez obtenida dicha descomposición, resolver el sistema es sencillo: \[ \mathbf{A}\mathbf{x}=\mathbf{b},\ \Rightarrow \mathbf{Q}\mathbf{R}\mathbf{x}=\mathbf{b},\ \Rightarrow \mathbf{R}\mathbf{x}=\mathbf{Q}^\top\mathbf{b}, \] y solamente tenemos que resolver un sistema triangular superior usando el método de sustitución hacia atrás explicado anteriormente.

Ejemplo

Consideremos el siguiente sistema de ecuaciones con \(5\) ecuaciones y \(5\) incógnitas: \[ \begin{align*} E_1: & 5x_1 & + x_2 & -3 x_3 & - x_4 & -2 x_5 & = 2,\\ E_2: & x_1 & +2 x_2 & +2 x_3 & & +3 x_5 & = 2,\\ E_3: & x_1 & - x_2 & +3 x_3& +2 x_4& & = 2,\\ E_4: & & x_2 & & +4 x_4 & - x_5 & = 3,\\ E_5: & x_1 & +3 x_2& + x_3 & & +5 x_5 & = 1. \end{align*} \] La matriz del sistema es la siguiente: \(\begin{bmatrix}5&1&-3&-1&-2 \\1&2&2&0&3 \\1&-1&3&2&0 \\0&1&0&4&-1 \\1&3&1&0&5 \\\end{bmatrix}\) con la siguiente descomposición \(QR\):

Ejemplo (continuación)

\[ \begin{align*} \mathbf{A} & = \mathbf{Q}\cdot \mathbf{R},\\ \begin{bmatrix}5&1&-3&-1&-2 \\1&2&2&0&3 \\1&-1&3&2&0 \\0&1&0&4&-1 \\1&3&1&0&5 \\\end{bmatrix}& =\begin{bmatrix}0.944911&-0.167701&-0.278534&0.024779&-0.028702 \\0.188982&0.463645&0.41874&-0.206152&-0.729025 \\0.188982&-0.364997&0.853128&0.202194&0.249705 \\0&0.276214&-0.068225&0.949537&-0.132028 \\0.188982&0.739859&0.120802&-0.11994&0.622828 \\\end{bmatrix}\cdot \\ & \begin{bmatrix}5.291503&1.70084&-1.70084&-0.566947&-0.377964 \\0&3.620379&1.075262&0.542564&5.149422 \\0&0&4.353269&1.71189&2.485526 \\0&0&0&4.177754&-2.217249 \\0&0&0&0&1.116499 \\\end{bmatrix}. \end{align*} \]

Ejemplo (continuación)

El sistema a resolver será: \[ \begin{align*} & \mathbf{R}\mathbf{x} =\mathbf{Q}^\top\cdot\mathbf{b},\\ & \begin{bmatrix}5.291503&1.70084&-1.70084&-0.566947&-0.377964 \\0&3.620379&1.075262&0.542564&5.149422 \\0&0&4.353269&1.71189&2.485526 \\0&0&0&4.177754&-2.217249 \\0&0&0&0&1.116499 \\\end{bmatrix} \begin{bmatrix}x_1\\ x_2\\ x_3\\ x_4\\ x_5\end{bmatrix} = \\ & \begin{bmatrix}0.944911&0.188982&0.188982&0&0.188982 \\-0.167701&0.463645&-0.364997&0.276214&0.739859 \\-0.278534&0.41874&0.853128&-0.068225&0.120802 \\0.024779&-0.206152&0.202194&0.949537&-0.11994 \\-0.028702&-0.729025&0.249705&-0.132028&0.622828 \\\end{bmatrix}\cdot\begin{bmatrix}2 \\2 \\2 \\3 \\1 \\\end{bmatrix}=\begin{bmatrix}2.834734 \\1.430395 \\1.902795 \\2.770313 \\-0.789298 \\\end{bmatrix}, \end{align*} \]

Ejemplo (continuación)

La solución del sistema será usando el método de sustitución hacia atrás: \[ \begin{align*} x_5 & =\frac{-0.789298}{1.116499}=-0.706941,\\ x_4 & =\frac{2.770313-(-2.2172487)\cdot (-0.7069409)}{4.1777542}=0.2879177,\\ x_3 & =\frac{1.9027949-(2.4855259)\cdot (-0.7069409)-1.7118895\cdot 0.2879177}{4.3532693}=0.7275064,\\ x_2 & =\frac{1.4303949-5.1494217\cdot (-0.7069409)-0.5425636\cdot 0.2879177-1.0752624\cdot 0.7275064}{3.6203788}\\ & =1.1413882,\\ x_1 & =\frac{1}{5.2915026}(2.8347335-(-0.3779645)\cdot (-0.7069409)-(-0.5669467)\cdot 0.2879177\\ & \quad -(-1.7008401)\cdot 0.7275064 -1.7008401\cdot 1.1413882)=0.3830334. \end{align*} \]

Método de ortogonalización de Gram-Schmidt

El método de descomposición \(QR\) se basa en el método de ortogonalización de Gram-Schmidt. Dicho método resuelve el problema siguiente:

Problema:

Sean \(\mathbf{v}_1,\ldots,\mathbf{v}_{k}\in\mathbb{R}^n\), \(k\) vectores linealmente independientes en el espacio vectorial \(\mathbb{R}^n\) (\(n\geq k\)). Queremos hallar \(\mathbf{u}_1,\ldots,\mathbf{u}_k\in\mathbb{R}^n\), que cumplen lo siguiente:

Son ortogonales dos a dos: \(\mathbf{u}_i^\top\cdot\mathbf{u}_j =0\), si \(i\neq j\).
Son ortonormales: \(\|\mathbf{u}_i\|=1\), \(i=1,\ldots,k\).
El subespacio vectorial generado por \(\mathbf{u}_1,\ldots,\mathbf{u}_i\) es el mismo que el subespacio vectorial generado por \(\mathbf{v}_1,\ldots,\mathbf{v}_i\), para \(i=1,\ldots,k\): \[ \begin{align*} & <\mathbf{u}_1,\ldots,\mathbf{u}_i> =\{\mathbf{x}=\alpha_1\mathbf{u}_1+\cdots +\alpha_i\mathbf{u}_i,\ \alpha_1,\ldots,\alpha_i\in\mathbb{R}\}\\ & =<\mathbf{v}_1,\ldots,\mathbf{v}_i>=\{\mathbf{x}=\beta_1\mathbf{v}_1+\cdots +\beta_i\mathbf{v}_i,\ \beta_1,\ldots,\beta_i\in\mathbb{R}\}. \end{align*} \]

Método de ortogonalización de Gram-Schmidt

Paso 1: Elegimos \(\mathbf{u}_1'=\mathbf{v}_1\) y \(\mathbf{u}_1=\frac{\mathbf{u}_1'}{\|\mathbf{u}_1'\|}\).
Paso 2: como \(<\mathbf{u}_1,\mathbf{u}_2>=<\mathbf{v}_1,\mathbf{v}_2>=<\mathbf{u_1},\mathbf{v}_2>\), definimos: \(\mathbf{u}_2'=\mathbf{v}_2+\lambda_{12}\mathbf{u}_1\) y hallamos \(\lambda_{12}\) imponiendo que \(\mathbf{u}_1\) y \(\mathbf{u}_2'\) deben ser ortogonales: \[ \begin{align*} & \mathbf{u}_1^\top\cdot\mathbf{u}_2' =0,\ \Rightarrow \mathbf{u}_1^\top\cdot (\mathbf{v}_2+\lambda_{12}\mathbf{u}_1)=\mathbf{u}_1^\top\cdot\mathbf{v}_2 +\lambda_{12}\cdot\mathbf{u}_1^\top\cdot\mathbf{u}_1 =0,\ \Rightarrow\\ & \lambda_{12} =-\frac{\mathbf{u}_1^\top\cdot\mathbf{v}_2}{\mathbf{u}_1^\top\cdot\mathbf{u}_1}=-\frac{\mathbf{u}_1^\top\cdot\mathbf{v}_2}{\|\mathbf{u}_1\|^2}=-\mathbf{u}_1^\top\cdot\mathbf{v}_2. \end{align*} \]

Entonces, \[ \mathbf{u}_2'=\mathbf{v}_2-(\mathbf{u}_1^\top\cdot\mathbf{v}_2)\cdot\mathbf{u}_1, \]

Método de ortogonalización de Gram-Schmidt

y por último, para hacer \(\mathbf{u}_2\) ortonormal, o que su norma euclídea sea \(1\), dividimos \(\mathbf{u}_2'\) por su norma: \[ \mathbf{u}_2 =\frac{\mathbf{u}_2'}{\|\mathbf{u}_2'\|}=\frac{\mathbf{v}_2-(\mathbf{u}_1^\top\cdot\mathbf{v}_2)\cdot\mathbf{u}_1}{\|\mathbf{v}_2-(\mathbf{u}_1^\top\cdot\mathbf{v}_2)\cdot\mathbf{u}_1\|} \]

Paso \(i\)-ésimo. Supongamos que hemos hallado \(\mathbf{u}_1,\ldots,\mathbf{u}_{i-1}\) que cumplen las condiciones anteriores. Vamos a hallar el valor siguiente \(\mathbf{u}_i\).

Como \(<\mathbf{u}_1,\ldots,\mathbf{u}_{i-1},\mathbf{u}_i>=<\mathbf{v}_1,\ldots,\mathbf{v}_{i-1},\mathbf{v}_i>=<\mathbf{u}_1,\ldots,\mathbf{u}_{i-1},\mathbf{v}_i>\), definimos: \[ \mathbf{u}_i' =\mathbf{v_i}+\sum_{j=1}^{i-1}\lambda_{ji}\cdot\mathbf{u}_j. \] Vamos a hallar los valores \(\lambda_{ji}\).

Método de ortogonalización de Gram-Schmidt

Como \(\mathbf{u}_l^\top\cdot\mathbf{u}_i' =0\), para \(l=1,\ldots,i-1\), tenemos que: \[ \mathbf{u}_l^\top\cdot\mathbf{u}_i' =\mathbf{u}_l^\top\cdot\left(\mathbf{v_i}+\sum_{j=1}^{i-1}\lambda_{ji}\cdot\mathbf{u}_j\right)=\mathbf{u}_l^\top\cdot\mathbf{v_i} +\lambda_{li}\cdot(\mathbf{u}_l^\top\cdot\mathbf{u}_l) =0, \] y, por tanto, el valor de \(\lambda_{li}\) vale: \[ \lambda_{li} = -\frac{\mathbf{u}_l^\top\cdot\mathbf{v}_i}{\mathbf{u}_l^\top\cdot\mathbf{u}^l}=-\frac{\mathbf{u}_l^\top\cdot\mathbf{v}_i}{\|\mathbf{u}_l\|^2}=-\mathbf{u}_l^\top\cdot\mathbf{v}_i. \] Por tanto, \[ \mathbf{u}_i' =\mathbf{v_i}-\sum_{j=1}^{i-1}(\mathbf{u}_j^\top\cdot\mathbf{v}_i)\cdot\mathbf{u}_j. \]

Método de ortogonalización de Gram-Schmidt

y por último, para hacer \(\mathbf{u}_i\) ortonormal, o que su norma euclídea sea \(1\), lo dividimos por su norma: \[ \mathbf{u}_i =\frac{\mathbf{u}_i'}{\|\mathbf{u}_i'\|}=\frac{\mathbf{v_i}-\sum_{j=1}^{i-1}(\mathbf{u}_j^\top\cdot\mathbf{v}_i)\cdot\mathbf{u}_j}{\|\mathbf{v_i}-\sum_{j=1}^{i-1}(\mathbf{u}_j^\top\cdot\mathbf{v}_i)\cdot\mathbf{u}_j\|}. \]

Ejemplo

Hallemos los vectores resultantes de aplicar el método de ortogonalización de Gram-Schmidt a las columnas de la matriz del sistema de ecuaciones del ejemplo anterior: \[ \mathbf{A}=\begin{bmatrix}5&1&-3&-1&-2 \\1&2&2&0&3 \\1&-1&3&2&0 \\0&1&0&4&-1 \\1&3&1&0&5 \\\end{bmatrix} \] Los vectores \(\mathbf{v}_1,\mathbf{v}_2,\mathbf{v}_3,\mathbf{v}_4,\mathbf{v}_5\) serían los siguientes: \[ \mathbf{v}_1=\begin{bmatrix}5 \\1 \\1 \\0 \\1 \\\end{bmatrix},\ \mathbf{v}_2=\begin{bmatrix}1 \\2 \\-1 \\1 \\3 \\\end{bmatrix},\ \mathbf{v}_3=\begin{bmatrix}-3 \\2 \\3 \\0 \\1 \\\end{bmatrix},\ \mathbf{v}_4=\begin{bmatrix}-1 \\0 \\2 \\4 \\0 \\\end{bmatrix},\ \mathbf{v}_5=\begin{bmatrix}-2 \\3 \\0 \\-1 \\5 \\\end{bmatrix}. \]

Ejemplo (continuación)

El vector \(\mathbf{u}_1'\) sería \(\mathbf{u}_1'=\mathbf{v}_1\) y el vector \(\mathbf{u}_1\) sería: \[ \mathbf{u}_1 =\frac{\mathbf{u}_1'}{\|\mathbf{u}_1'\|}=\frac{1}{\sqrt{5^2+ 1^2+1^2+0^2+1^2}}\begin{bmatrix}5 \\1 \\1 \\0 \\1 \\\end{bmatrix}=\frac{1}{5.2915026}\begin{bmatrix}5 \\1 \\1 \\0 \\1 \\\end{bmatrix}=\begin{bmatrix}0.944911 \\0.188982 \\0.188982 \\0 \\0.188982 \\\end{bmatrix}. \]
Para calcular \(\mathbf{u}_2\), en primer lugar calculamos \(\lambda_{12}=-\mathbf{u}_1^\top\cdot \mathbf{v}_2\):

\[ \begin{align*} -\mathbf{u}_1^\top\cdot \mathbf{v}_2 & =-\begin{bmatrix}0.944911&0.188982&0.188982&0&0.188982 \\\end{bmatrix}\cdot \begin{bmatrix}1 \\2 \\-1 \\1 \\3 \\\end{bmatrix}\\ & =-(0.944911\cdot1+0.188982\cdot2+0.188982\cdot (-1)+0\cdot1+0.188982\cdot3)=-1.70084. \end{align*} \]

Ejemplo (continuación)

El valor de \(\mathbf{u}_2'\) será: \[ \mathbf{u}_2'=\mathbf{v}_2+\lambda_{12}\cdot\mathbf{u}_1 =\begin{bmatrix}1 \\2 \\-1 \\1 \\3 \\\end{bmatrix}-1.70084\cdot \begin{bmatrix}0.944911 \\0.188982 \\0.188982 \\0 \\0.188982 \\\end{bmatrix}=\begin{bmatrix}-0.607143 \\1.678571 \\-1.321429 \\1 \\2.678571 \\\end{bmatrix}. \] Por último, hacemos que el vector \(\mathbf{u}_2\) sea ortonormal o tenga norma euclídea \(1\):

\[ \begin{align*} \mathbf{u}_2 & = \frac{\mathbf{u}_2'}{\|\mathbf{u}_2'\|}=\frac{1}{\|\begin{bmatrix}-0.607143&1.678571&-1.321429&1&2.678571 \\\end{bmatrix}\|}\cdot \begin{bmatrix}-0.607143 \\1.678571 \\-1.321429 \\1 \\2.678571 \\\end{bmatrix}\\ & =\frac{1}{3.620379}\cdot \begin{bmatrix}-0.607143 \\1.678571 \\-1.321429 \\1 \\2.678571 \\\end{bmatrix}=\begin{bmatrix}-0.167701 \\0.463645 \\-0.364997 \\0.276214 \\0.739859 \\\end{bmatrix} \end{align*} \]

Ejemplo (continuación)

Recordemos que para calcular \(\mathbf{u}_3\), primero calculamos \(\mathbf{u}_3'\): \[ \mathbf{u}_3'=\mathbf{v}_3+\lambda_{13}\cdot\mathbf{u}_{13}+\lambda_{23}\cdot\mathbf{u}_2. \] Los valores de \(\lambda_{13}\) y \(\lambda_{23}\) son los siguientes: \[ \begin{align*} \lambda_{13} & = -\mathbf{u}_1^\top\cdot\mathbf{v}_3=-\begin{bmatrix}0.944911&0.188982&0.188982&0&0.188982 \\\end{bmatrix}\cdot \begin{bmatrix}-3 \\2 \\3 \\0 \\1 \\\end{bmatrix}\\ & =-(0.944911\cdot (-3)+0.188982\cdot2+0.188982\cdot 3+0\cdot0+0.188982\cdot1)=1.70084,\\ \lambda_{23} & = -\mathbf{u}_2^\top\cdot\mathbf{v}_3=-\begin{bmatrix}-0.167701&0.463645&-0.364997&0.276214&0.739859 \\\end{bmatrix}\cdot \begin{bmatrix}-3 \\2 \\3 \\0 \\1 \\\end{bmatrix}\\ & =-(-0.167701\cdot (-3)+0.463645\cdot2+(-0.364997)\cdot 3+0.276214\cdot0+0.739859\cdot1)\\ & =-1.075262. \end{align*} \]

Ejemplo (continuación)

El valor de \(\mathbf{u}_3'\) será: \[ \begin{align*} \mathbf{u}_3' & =\mathbf{v}_3+\lambda_{13}\cdot\mathbf{u}_1+\lambda_{23}\cdot\mathbf{u}_2 =\begin{bmatrix}-3 \\2 \\3 \\0 \\1 \\\end{bmatrix}+1.70084\cdot \begin{bmatrix}0.944911 \\0.188982 \\0.188982 \\0 \\0.188982 \\\end{bmatrix}+(-1.075262)\cdot \begin{bmatrix}-0.167701 \\0.463645 \\-0.364997 \\0.276214 \\0.739859 \\\end{bmatrix}\\ & =\begin{bmatrix}-1.212534 \\1.822888 \\3.713896 \\-0.297003 \\0.525886 \\\end{bmatrix}. \end{align*} \]

Ejemplo (continuación)

Por último, hacemos que el vector \(\mathbf{u}_3\) sea ortonormal o tenga norma euclídea \(1\): \[ \begin{align*} \mathbf{u}_3 & = \frac{\mathbf{u}_3'}{\|\mathbf{u}_3'\|}=\frac{1}{\|\begin{bmatrix}-1.212534&1.822888&3.713896&-0.297003&0.525886 \\\end{bmatrix}\|}\cdot \begin{bmatrix}-1.212534 \\1.822888 \\3.713896 \\-0.297003 \\0.525886 \\\end{bmatrix}\\ & =\frac{1}{4.353269}\cdot \begin{bmatrix}-1.212534 \\1.822888 \\3.713896 \\-0.297003 \\0.525886 \\\end{bmatrix}=\begin{bmatrix}-0.278534 \\0.41874 \\0.853128 \\-0.068225 \\0.120802 \\\end{bmatrix} \end{align*} \]

Ejemplo (continuación)

Recordemos que para calcular \(\mathbf{u}_4\), primero calculamos \(\mathbf{u}_4'\): \[ \mathbf{u}_4'=\mathbf{v}_4+\lambda_{14}\cdot\mathbf{u}_1+\lambda_{24}\cdot\mathbf{u}_2+\lambda_{34}\cdot\mathbf{u}_3. \] Los valores de \(\lambda_{14}\), \(\lambda_{24}\) y \(\lambda_{34}\) son los siguientes: \[ \begin{align*} \lambda_{14} & = -\mathbf{u}_1^\top\cdot\mathbf{v}_4=-\begin{bmatrix}0.944911&0.188982&0.188982&0&0.188982 \\\end{bmatrix}\cdot \begin{bmatrix}-1 \\0 \\2 \\4 \\0 \\\end{bmatrix}\\ & =-(0.944911\cdot (-1)+0.188982\cdot0+0.188982\cdot 2+0\cdot4+0.188982\cdot0)=0.566947,\\ \lambda_{24} & = -\mathbf{u}_2^\top\cdot\mathbf{v}_4=-\begin{bmatrix}-0.167701&0.463645&-0.364997&0.276214&0.739859 \\\end{bmatrix}\cdot \begin{bmatrix}-1 \\0 \\2 \\4 \\0 \\\end{bmatrix}\\ & =-(-0.167701\cdot (-1)+0.463645\cdot0+(-0.364997)\cdot 2+0.276214\cdot4+0.739859\cdot0) =-0.542564,\\ \lambda_{34} & = -\mathbf{u}_3^\top\cdot\mathbf{v}_4=-\begin{bmatrix}-0.278534&0.41874&0.853128&-0.068225&0.120802 \\\end{bmatrix}\cdot \begin{bmatrix}-1 \\0 \\2 \\4 \\0 \\\end{bmatrix}\\ & =-(-0.278534\cdot (-1)+0.41874\cdot0+0.853128\cdot 2+(-0.068225)\cdot4+0.120802\cdot0) =-1.71189. \end{align*} \]

Ejemplo (continuación)

El valor de \(\mathbf{u}_4'\) será: \[ \begin{align*} \mathbf{u}_4' & =\mathbf{v}_4+\lambda_{14}\cdot\mathbf{u}_1+\lambda_{24}\cdot\mathbf{u}_2+\lambda_{34}\cdot\mathbf{u}_3 =\begin{bmatrix}-1 \\0 \\2 \\4 \\0 \\\end{bmatrix}+0.566947\cdot \begin{bmatrix}0.944911 \\0.188982 \\0.188982 \\0 \\0.188982 \\\end{bmatrix}+(-0.542564)\cdot \begin{bmatrix}-0.167701 \\0.463645 \\-0.364997 \\0.276214 \\0.739859 \\\end{bmatrix}\\ & \qquad +(-1.71189)\cdot \begin{bmatrix}-0.278534 \\0.41874 \\0.853128 \\-0.068225 \\0.120802 \\\end{bmatrix} =\begin{bmatrix}0.103523 \\-0.861251 \\0.844716 \\3.96693 \\-0.501078 \\\end{bmatrix}. \end{align*} \]

Ejemplo (continuación)

Por último, hacemos que el vector \(\mathbf{u}_4\) sea ortonormal o tenga norma euclídea \(1\): \[ \begin{align*} \mathbf{u}_4 & = \frac{\mathbf{u}_4'}{\|\mathbf{u}_4'\|}=\frac{1}{\|\begin{bmatrix}0.103523&-0.861251&0.844716&3.96693&-0.501078 \\\end{bmatrix}\|}\cdot \begin{bmatrix}0.103523 \\-0.861251 \\0.844716 \\3.96693 \\-0.501078 \\\end{bmatrix}\\ & =\frac{1}{4.177754}\cdot \begin{bmatrix}0.103523 \\-0.861251 \\0.844716 \\3.96693 \\-0.501078 \\\end{bmatrix}=\begin{bmatrix}0.024779 \\-0.206152 \\0.202194 \\0.949537 \\-0.11994 \\\end{bmatrix}. \end{align*} \]

Ejemplo (continuación)

Recordemos que para calcular \(\mathbf{u}_5\), primero calculamos \(\mathbf{u}_5'\): \[ \mathbf{u}_5'=\mathbf{v}_5+\lambda_{15}\cdot\mathbf{u}_1+\lambda_{25}\cdot\mathbf{u}_2+\lambda_{35}\cdot\mathbf{u}_3+\lambda_{45}\cdot\mathbf{u}_4. \] Los valores de \(\lambda_{15}\), \(\lambda_{25}\), \(\lambda_{35}\) y \(\lambda_{45}\) son los siguientes: \[ \begin{align*} \lambda_{15} & = -\mathbf{u}_1^\top\cdot\mathbf{v}_5=-\begin{bmatrix}0.944911&0.188982&0.188982&0&0.188982 \\\end{bmatrix}\cdot \begin{bmatrix}-2 \\3 \\0 \\-1 \\5 \\\end{bmatrix}\\ & =-(0.944911\cdot (-2)+0.188982\cdot3+0.188982\cdot 0+0\cdot (-1)+0.188982\cdot5)=0.377964,\\ \lambda_{25} & = -\mathbf{u}_2^\top\cdot\mathbf{v}_5=-\begin{bmatrix}-0.167701&0.463645&-0.364997&0.276214&0.739859 \\\end{bmatrix}\cdot \begin{bmatrix}-2 \\3 \\0 \\-1 \\5 \\\end{bmatrix}\\ & =-(-0.167701\cdot (-2)+0.463645\cdot3+(-0.364997)\cdot 0+0.276214\cdot (-1)+0.739859\cdot5) =-5.149422,\\ \lambda_{35} & = -\mathbf{u}_3^\top\cdot\mathbf{v}_5=-\begin{bmatrix}-0.278534&0.41874&0.853128&-0.068225&0.120802 \\\end{bmatrix}\cdot \begin{bmatrix}-2 \\3 \\0 \\-1 \\5 \\\end{bmatrix}\\ & =-(-0.278534\cdot (-2)+0.41874\cdot3+0.853128\cdot 0+(-0.068225)\cdot (-1)+0.120802\cdot5) =-2.485526, \end{align*} \]

Ejemplo (continuación)

\[ \begin{align*} \lambda_{45} & = -\mathbf{u}_4^\top\cdot\mathbf{v}_5=-\begin{bmatrix}0.024779&-0.206152&0.202194&0.949537&-0.11994 \\\end{bmatrix}\cdot \begin{bmatrix}-2 \\3 \\0 \\-1 \\5 \\\end{bmatrix}\\ & =-(0.024779\cdot (-2)+(-0.206152)\cdot3+0.202194\cdot 0+0.949537\cdot (-1)+(-0.11994)\cdot5) =2.217249. \end{align*} \]

El valor de \(\mathbf{u}_5'\) será: \[ \begin{align*} \mathbf{u}_5' & =\mathbf{v}_5+\lambda_{15}\cdot\mathbf{u}_1+\lambda_{25}\cdot\mathbf{u}_2+\lambda_{35}\cdot\mathbf{u}_3+\lambda_{45}\cdot\mathbf{u}_4 =\begin{bmatrix}-2 \\3 \\0 \\-1 \\5 \\\end{bmatrix}+0.377964\cdot \begin{bmatrix}0.944911 \\0.188982 \\0.188982 \\0 \\0.188982 \\\end{bmatrix}\\ & \qquad +(-5.149422)\cdot \begin{bmatrix}-0.167701 \\0.463645 \\-0.364997 \\0.276214 \\0.739859 \\\end{bmatrix} +(-2.485526)\cdot \begin{bmatrix}-0.278534 \\0.41874 \\0.853128 \\-0.068225 \\0.120802 \\\end{bmatrix}+2.217249\cdot \begin{bmatrix}0.024779 \\-0.206152 \\0.202194 \\0.949537 \\-0.11994 \\\end{bmatrix} =\begin{bmatrix}-0.032045 \\-0.813955 \\0.278796 \\-0.147409 \\0.695387 \\\end{bmatrix}. \end{align*} \]

Ejemplo (continuación)

Por último, hacemos que el vector \(\mathbf{u}_5\) sea ortonormal o tenga norma euclídea \(1\): \[ \begin{align*} \mathbf{u}_5 & = \frac{\mathbf{u}_5'}{\|\mathbf{u}_5'\|}=\frac{1}{\|\begin{bmatrix}-0.032045&-0.813955&0.278796&-0.147409&0.695387 \\\end{bmatrix}\|}\cdot \begin{bmatrix}-0.032045 \\-0.813955 \\0.278796 \\-0.147409 \\0.695387 \\\end{bmatrix}\\ & =\frac{1}{1.116499}\cdot \begin{bmatrix}-0.032045 \\-0.813955 \\0.278796 \\-0.147409 \\0.695387 \\\end{bmatrix}=\begin{bmatrix}-0.028702 \\-0.729025 \\0.249705 \\-0.132028 \\0.622828 \\\end{bmatrix}. \end{align*} \]

Ejemplo (continuación)

En resumen, hemos convertido las columnas de la matriz \[ \mathbf{A}=\begin{bmatrix}5&1&-3&-1&-2 \\1&2&2&0&3 \\1&-1&3&2&0 \\0&1&0&4&-1 \\1&3&1&0&5 \\\end{bmatrix}, \] en las columnas de la matriz \[ \mathbf{Q}=\begin{bmatrix}0.944911&-0.167701&-0.278534&0.024779&-0.028702 \\0.188982&0.463645&0.41874&-0.206152&-0.729025 \\0.188982&-0.364997&0.853128&0.202194&0.249705 \\0&0.276214&-0.068225&0.949537&-0.132028 \\0.188982&0.739859&0.120802&-0.11994&0.622828 \\\end{bmatrix}, \] de tal forma que si construimos la submatriz formada por las \(k\) primeras columnas de la matriz \(\mathbf{A}\) y las \(k\) primeras columnas de la matriz \(\mathbf{Q}\), es decir, dicha submatriz tendrá \(5\) filas y \(2k\) columnas, entonces el rango de dicha submatriz será siempre \(k\), con \(k=1,2,3,4,5\).

Ejemplo (continuación)

Ejercicio

Verificar la afirmación anterior para \(k=3\), es decir, demostrar que el rango de la submatriz siguiente es \(3\): \[ \begin{bmatrix}5&1&-3&0.944911&-0.167701&-0.278534 \\1&2&2&0.188982&0.463645&0.41874 \\1&-1&3&0.188982&-0.364997&0.853128 \\0&1&0&0&0.276214&-0.068225 \\1&3&1&0.188982&0.739859&0.120802 \\\end{bmatrix}. \]

Descomposición \(QR\)

El método de ortogonalización de Gram-Schmidt es la clave para hallar la descomposición \(QR\) de una matriz \(\mathbf{A}\):

Proposición. Descomposición \(QR\).

Sea \(\mathbf{A}\) una matriz cuadrada \(n\times n\) no singular. Aplicamos la ortogonalización de Gram-Schmidt a las columnas de la matriz \(\mathbf{A}\).

La matrices \(\mathbf{Q}\) y \(\mathbf{R}\) de la descomposición \(QR\) de la matriz \(\mathbf{A}\) son las siguientes:

La matriz \(\mathbf{Q}\) es la matriz cuyas columnas están formadas por los vectores ortonormales \(\mathbf{u}_1,\ldots,\mathbf{u}_n\) obtenidos mediante el proceso de Gram-Shmidt y

Descomposición \(QR\)

Proposición (continuación). Descomposición \(QR\).

la matriz \(\mathbf{R}\) es la siguiente:

\[ \mathbf{R}=\begin{bmatrix}\|\mathbf{u}_{1}'\|&-\lambda_{12}&\ldots&-\lambda_{1n} \\0&\|\mathbf{u}_{2}'\|&\ldots&-\lambda_{2n} \\\vdots&\vdots&\ddots&\vdots \\0&\ldots&0&\|\mathbf{u}_{n}'\| \\\end{bmatrix}. \] Es decir, la matriz \(\mathbf{R}\) tiene como diagonal la norma euclídea de los vectores \(\mathbf{u}_i'\), \(i=1,\ldots,n\) de la ortogonalización de Gram-Schmidt y \(r_{ij}=-\lambda_{ij}\), para los elementos fuera de la diagonal con \(j>i\).

Demostración de la proposición

Sean \(\mathbf{v}_1,\ldots,\mathbf{v}_n\) la columnas de la matriz \(\mathbf{A}\).

Recordemos la relación que hay entre los vectores \(\mathbf{v}_i\) y los vectores \(\mathbf{u}_i\): \[ \mathbf{u}_i'=\mathbf{v}_i+\sum_{j=1}^{i-1}\lambda_{ji}\mathbf{u}_j, \] para \(i=1,\ldots,n\).

Si despejamos \(\mathbf{v}_i\) de la expresión anterior, obtenemos: \[ \mathbf{v}_i=\mathbf{u}_i'-\sum_{j=1}^{i-1}\lambda_{ji}\mathbf{u}_j=\|\mathbf{u}_i'\|\cdot \mathbf{u}_i -\sum_{j=1}^{i-1}\lambda_{ji}\mathbf{u}_j. \] Si escribimos la expresión anterior en componentes, obtenemos: \[ v_{ik} = \|\mathbf{u}_i'\|\cdot u_{ik} -\sum_{j=1}^{i-1}\lambda_{ji}u_{jk}, \] para \(k=1,\ldots,n\), e \(i=1,\ldots,n\), donde \(v_{ik}\) es la componente \(k\)-ésima del vector \(\mathbf{v}_i\) y \(u_{jk}\) es la componente \(k\)-ésima del vector \(\mathbf{u}_j\).

Demostración de la proposición (continuación)

Recordemos que \(\mathbf{v}_i\) era la columna \(i\)-ésima de la matriz \(\mathbf{A}\). Entonces si escribimos \(\mathbf{A}=(a_{ij})_{i,j=1,\ldots,n}\), \(a_{ij}=v_{ji}\).

Sea la matriz \(\mathbf{Q}\) ortogonal cuyas columnas son los vectores ortonormales \(\mathbf{u}_1,\ldots,\mathbf{u}_n\). Entonces si escribimos \(\mathbf{Q}=(q_{ij})_{i,j=1,\ldots,n}\), \(q_{ij}=u_{ji}\).

Por último definimos la matriz \(\mathbf{R}=(r_{ij})_{i,j=1,\ldots,n}\) tal como está definida en la proposición.

Entonces, la expresión anterior puede escribirse de la forma siguiente: \[ a_{ki} = \|\mathbf{u}_i'\|\cdot q_{ki} -\sum_{j=1}^{i-1}\lambda_{ji}q_{kj}=\sum_{j=1}^{i-1}q_{kj}r_{ji}+q_{ki}\cdot r_{ii}=\sum_{j=1}^i q_{kj}r_{ji}. \] La expresión anterior es equivalente a afirmar que \(\mathbf{A}=\mathbf{Q}\mathbf{R}\), tal como queríamos demostrar.

Método \(QR\). Pseudocódigo

INPUT número de incógnitas y ecuaciones \(n\), matriz \(\mathbf{A}=(a_{ij})_{i=1,\ldots,n,j=1,\ldots,n}\), vector de términos independientes \(\mathbf{b}=(b_i)_{i=1,\ldots,n}\).
Set \(\mathbf{Q}=\mathbf{0}\) (inicializamos la matriz \(\mathbf{Q}\) a \(\mathbf{0}\))
Set \(\mathbf{R}=\mathbf{0}\) (inicializamos la matriz \(\mathbf{R}\) a \(\mathbf{0}\))
Set \(\mathbf{Q}[:,1]=\mathbf{A}[:,1]/\|\mathbf{A}[:,1]\|_2\). (la primera columna de la matriz \(\mathbf{Q}\) es la primera columna de la matriz \(\mathbf{A}\) normalizada)
Set \(\mathbf{R}[1,1]=\|\mathbf{A}[:,1]\|_2\) (definimos \(R[1,1]\) como la norma euclídea de la primera columna de la matriz \(\mathbf{A}\))

Método \(QR\). Pseudocódigo

For i=2,...,n
- For j=1,...,i-1
  - Set \(\mathbf{R}[j,i]=\mathbf{Q}[:,j]\cdot \mathbf{A}[:,i]\) (la componente \((j,i)\) de la matriz \(R\) es el producto escalar de la columna \(j\)-ésima de la matriz \(\mathbf{Q}\) por la columna \(i\)-ésima de la matriz \(\mathbf{A}\) cambiado de signo.)
- Set \(\mathbf{u}_i'=\mathbf{A}[:,i]-\sum_{j=1}^{i-1}\mathbf{R}[j,i]\cdot\mathbf{Q}[:,j]\) (calculamos el vector \(\mathbf{u}_i'\))
- Set \(\mathbf{R}[i,i]=\|\mathbf{u}_i'\|_2\) (calculamos \(R[i,i]\))
- Set \(\mathbf{Q}[:,i]=\frac{\mathbf{u}_i'}{\|\mathbf{u}_i'\|_2}\) (definimos la columna \(i\)-ésima de la matriz \(\mathbf{Q}\))
Resolvemos \(\mathbf{R}\mathbf{x}=\mathbf{Q}^\top\mathbf{b}\) por sustitución hacia atrás.
Print \(\mathbf{x}\) (damos la solución del sistema)

Ejemplo anterior

Para la matriz del ejemplo anterior: \[ \mathbf{A}=\begin{bmatrix}5&1&-3&-1&-2 \\1&2&2&0&3 \\1&-1&3&2&0 \\0&1&0&4&-1 \\1&3&1&0&5 \\\end{bmatrix}, \] la matriz \(\mathbf{Q}\) sería la matriz cuyas columnas son los vectores \(\mathbf{u}_1,\mathbf{u}_2,\mathbf{u}_3,\mathbf{u}_4\) y \(\mathbf{u}_5\): \[ \mathbf{Q}=\begin{bmatrix}0.944911&-0.167701&-0.278534&0.024779&-0.028702 \\0.188982&0.463645&0.41874&-0.206152&-0.729025 \\0.188982&-0.364997&0.853128&0.202194&0.249705 \\0&0.276214&-0.068225&0.949537&-0.132028 \\0.188982&0.739859&0.120802&-0.11994&0.622828 \\\end{bmatrix}. \]

Ejemplo anterior (continuación)

Los valores \(\lambda_{ij}\) ya fueron calculados y eran: \[ \begin{align*} \lambda_{12} & = -1.70084,\ \lambda_{13}=1.70084,\ \lambda_{14}=0.566947,\ \lambda_{15}=0.377964,\\ \lambda_{23} & = -1.075262,\ \lambda_{24}=-0.542564,\ \lambda_{25}=-5.149422,\\ \lambda_{34} & = -1.71189,\ \lambda_{35}=-2.485526,\\ \lambda_{45} & = 2.217249. \end{align*} \]

Los valores de las normas euclídeas de los vectores \(\mathbf{u}_1',\mathbf{u}_2',\mathbf{u}_3',\mathbf{u}_4'\) y \(\mathbf{u}_5'\) valen:

\[ \|\mathbf{u}_1'\| =\|\mathbf{v}_1\| =5.291503,\ \|\mathbf{u}_2'\| =3.620379,\ \|\mathbf{u}_3'\| =4.353269,\ \|\mathbf{u}_4'\| =4.177754,\ \|\mathbf{u}_5'\| =1.116499. \] La matriz \(\mathbf{R}\) será, entonces, \[ \mathbf{R}=\begin{bmatrix}5.291503&1.70084&-1.70084&-0.566947&-0.377964 \\0&3.620379&1.075262&0.542564&5.149422 \\0&0&4.353269&1.71189&2.485526 \\0&0&0&4.177754&-2.217249 \\0&0&0&0&1.116499 \\\end{bmatrix}. \]

Aplicaciones

Introducción

Vamos a ver cómo podemos usar las técnicas de resolución de sistemas lineales por métodos directos para calcular:

la inversa de una matriz regular \(\mathbf{A}\) cuadrada \(n\times n\),
el determinante de una matriz regular \(\mathbf{A}\) cuadrada \(n\times n\).

Inversa de una matriz \(\mathbf{A}\)

Dada una matriz cuadrada \(n\times n\) y regular \(\mathbf{A}\), hallar la matriz inversa \(\mathbf{A}^{-1}\) equivale a hallar una matriz \(n\times n\), \(\mathbf{A}^{-1}\) tal que: \[ \mathbf{A}\cdot \mathbf{A}^{-1}=\mathbf{I}_n, \] donde \(\mathbf{I}_n\) es la matriz identidad \(n\times n\):

\[ \mathbf{I}_n=\begin{bmatrix}1&0&\ldots&0 \\0&1&\ldots&0 \\\vdots&\vdots&\ddots&\vdots \\0&0&\ldots&1 \\\end{bmatrix}. \]

Inversa de una matriz \(\mathbf{A}\)

Sea \(\mathbf{c}_i\) la columna \(i\)-ésima de la matriz inversa \(\mathbf{A}^{-1}\). Usando la condición anterior, tenemos que dicha columna \(\mathbf{c}_i\) verifica el siguiente sistema lineal: \[ \mathbf{A}\cdot\mathbf{c}_i =\mathbf{e}_i:=\begin{bmatrix}0\\\vdots\\\overbrace{1}^{i)}\\\vdots\\ 0\end{bmatrix}. \]

Inversa de una matriz \(\mathbf{A}\)

Dicho de otra manera, hallar la matriz inversa \(\mathbf{A}^{-1}\) es equivalente a resolver \(n\) sistemas lineales todos con la misma matriz del sistema \(\mathbf{A}\) pero con diferentes vectores de términos independientes.

La resolución del sistema anterior corresponde a la columna \(i\)-ésima de la matriz inversa \(\mathbf{A}^{-1}\), \(\mathbf{c}_i\).

La matriz inversa \(\mathbf{A}^{-1}\) será pues: \(\mathbf{A}^{-1}=(\mathbf{c}_1,\ldots,\mathbf{c}_n)\).

Inversa de una matriz \(\mathbf{A}\)

Como la matriz del sistema es la misma para los \(n\) sistemas lineales anteriores, podemos hallar la inversa usando los métodos vistos anteriormente:

Descomposición \(LU\). Si descomponemos la matriz \(\mathbf{A}\) de la forma \(\mathbf{A}=\mathbf{L}\mathbf{U}\), resolvemos \(\mathbf{L}\mathbf{U}\mathbf{c}_i = \mathbf{e}_i\), para \(i=1,\ldots,n\), usando el método de resolución \(LU\) visto anteriormente, es decir,
- hallamos primero \(\mathbf{y}_i\) resolviendo \(\mathbf{L}\mathbf{y}_i=\mathbf{e}_i\) y seguidamente,
- hallamos \(\mathbf{c}_i\) resolviendo \(\mathbf{U}\mathbf{c}_i=\mathbf{y}_i\).

Inversa de una matriz \(\mathbf{A}\)

Descomposición \(QR\). Si descomponemos la matriz \(\mathbf{A}\) de la forma \(\mathbf{A}=\mathbf{Q}\mathbf{R}\), hacemos lo siguiente para hallar \(\mathbf{A}^{-1}\): \[ \mathbf{A}\cdot\mathbf{A}^{-1}=\mathbf{Id}_n,\ \Rightarrow \mathbf{Q}\cdot\mathbf{R}\cdot\mathbf{A}^{-1}=\mathbf{Id}_n,\ \Rightarrow \mathbf{R}\cdot\mathbf{A}^{-1}= \mathbf{Q}^\top\cdot\mathbf{Id}_n=\mathbf{Q}^\top. \] En este segundo caso, para hallar la columna \(i\)-ésima de la matriz inversa \(\mathbf{c}_i\), hemos de resolver el sistema de ecuaciones \[ \mathbf{R}\cdot\mathbf{c}_i=\mathbf{qt}_i, \] donde \(\mathbf{qt}_i\) es la columna \(i\)-ésima de la matriz \(\mathbf{Q}^\top\) que correspondería a la fila \(i\)-ésima de la matriz \(\mathbf{Q}\). Fijémonos que la matriz del sistema \(\mathbf{R}\) es triangular superior. Por tanto, se puede resolver usando la técnica de sustitución hacia atrás.

Ejemplo

Vamos a hallar la inversa de la matriz vista anteriormente \[ \mathbf{A}=\begin{bmatrix}5&1&-3&-1&-2 \\1&2&2&0&3 \\1&-1&3&2&0 \\0&1&0&4&-1 \\1&3&1&0&5 \\\end{bmatrix}, \] usando la descomposición \(LU\). Dicha descomposición es la siguiente: \[ \mathbf{L}=\begin{bmatrix}1&0&0&0&0 \\0.2&1&0&0&0 \\0.2&-0.666667&1&0&0 \\0&0.555556&-0.270833&1&0 \\0.2&1.555556&-0.458333&0.211982&1 \\\end{bmatrix},\quad\mathbf{U}=\begin{bmatrix}5&1&-3&-1&-2 \\0&1.8&2.6&0.2&3.4 \\0&0&5.333333&2.333333&2.666667 \\0&0&0&4.520833&-2.166667 \\0&0&0&0&1.792627 \\\end{bmatrix}. \]

Ejemplo (continuación)

Primera columna de \(\mathbf{A}^{-1}\). Primero resolvemos: \[ \mathbf{L}\cdot\mathbf{y}_1 =\mathbf{e}_1,\ \Rightarrow \begin{bmatrix}1&0&0&0&0 \\0.2&1&0&0&0 \\0.2&-0.666667&1&0&0 \\0&0.555556&-0.270833&1&0 \\0.2&1.555556&-0.458333&0.211982&1 \\\end{bmatrix}\cdot \begin{bmatrix}y_{11}\\ y_{21}\\ y_{31}\\ y_{41}\\ y_{51}\end{bmatrix}=\begin{bmatrix}1\\0\\0\\0\\0\end{bmatrix}. \] Hallamos el vector \(\mathbf{y}_1\) usando la técnica de sustitución hacia adelante:

\[ \begin{align*} y_{11} & =1,\ y_{21}=0-0.2 y_{11}=-0.2,\ y_{31}=0+0.666667 y_{21}-0.2 y_{11}=-0.333333,\\ y_{41} & =0+0.270833 y_{31}-0.555556 y_{21}=0.0208333,\\ y_{51} & =0-0.211982 y_{41}+0.458333 y_{31}-1.555556 y_{21}-0.2 y_{11}=-0.0460829. \end{align*} \] El vector \(\mathbf{y}_1\) vale: \(\begin{bmatrix}1 \\-0.2 \\-0.333333 \\0.020833 \\-0.046083 \\\end{bmatrix}.\)

Ejemplo (continuación)

A continuación, calculamos la primera columna \(\mathbf{c}_1\) de la matriz \(\mathbf{A}^{-1}\) resolviendo: \[ \mathbf{U}\cdot\mathbf{c}_1=\mathbf{y}_1,\quad \begin{bmatrix}5&1&-3&-1&-2 \\0&1.8&2.6&0.2&3.4 \\0&0&5.333333&2.333333&2.666667 \\0&0&0&4.520833&-2.166667 \\0&0&0&0&1.792627 \\\end{bmatrix}\cdot\begin{bmatrix}c_{11}\\c_{21}\\c_{31}\\c_{41}\\c_{51}\end{bmatrix}=\begin{bmatrix}1 \\-0.2 \\-0.333333 \\0.020833 \\-0.046083 \\\end{bmatrix}. \]

\[ \begin{align*} c_{51} & =\frac{-0.046083}{1.792627}=-0.025707,\\ c_{41} & =\frac{0.020833-(-2.1666667)\cdot c_{51}}{4.5208333}=-0.0077121,\\ c_{31} & =\frac{-0.3333333-2.6666667\cdot c_{51}-2.3333333\cdot c_{41}}{5.3333333}=-0.0462725,\\ c_{21} & =\frac{-0.2-3.4\cdot c_{51}-0.2\cdot c_{41}-2.6\cdot c_{31}}{1.8}=0.0051414,\\ c_{11} & =\frac{1-(-2)\cdot c_{51}-(-1)\cdot c_{41}-(-3)\cdot c_{31}-1\cdot c_{21}}{5}=0.159383. \end{align*} \]

Ejemplo (continuación)

Segunda columna de \(\mathbf{A}^{-1}\). Primero resolvemos: \[ \mathbf{L}\cdot\mathbf{y}_2 =\mathbf{e}_2,\ \Rightarrow \begin{bmatrix}1&0&0&0&0 \\0.2&1&0&0&0 \\0.2&-0.666667&1&0&0 \\0&0.555556&-0.270833&1&0 \\0.2&1.555556&-0.458333&0.211982&1 \\\end{bmatrix}\cdot \begin{bmatrix}y_{12}\\ y_{22}\\ y_{32}\\ y_{42}\\ y_{52}\end{bmatrix}=\begin{bmatrix}0\\1\\0\\0\\0\end{bmatrix}. \] Hallamos el vector \(\mathbf{y}_2\) usando la técnica de sustitución hacia adelante:

\[ \begin{align*} y_{12} & =0,\ y_{22}=1-0.2 y_{12}=1,\ y_{32}=0+0.666667 y_{22}-0.2 y_{12}=0.666667,\\ y_{42} & =0+0.270833 y_{32}-0.555556 y_{22}=-0.375,\\ y_{52} & =0-0.211982 y_{42}+0.458333 y_{32}-1.555556 y_{22}-0.2 y_{12}=-1.1705069. \end{align*} \] El vector \(\mathbf{y}_2\) vale: \(\begin{bmatrix}0 \\1 \\0.666667 \\-0.375 \\-1.170507 \\\end{bmatrix}.\)

Ejemplo (continuación)

A continuación, calculamos la segunda columna \(\mathbf{c}_2\) de la matriz \(\mathbf{A}^{-1}\) resolviendo: \[ \mathbf{U}\cdot\mathbf{c}_2=\mathbf{y}_2,\quad \begin{bmatrix}5&1&-3&-1&-2 \\0&1.8&2.6&0.2&3.4 \\0&0&5.333333&2.333333&2.666667 \\0&0&0&4.520833&-2.166667 \\0&0&0&0&1.792627 \\\end{bmatrix}\cdot\begin{bmatrix}c_{12}\\c_{22}\\c_{32}\\c_{42}\\c_{52}\end{bmatrix}=\begin{bmatrix}0 \\1 \\0.666667 \\-0.375 \\-1.170507 \\\end{bmatrix}. \]

\[ \begin{align*} c_{52} & =\frac{-1.170507}{1.792627}=-0.652956,\\ c_{42} & =\frac{-0.375-(-2.1666667)\cdot c_{52}}{4.5208333}=-0.3958869,\\ c_{32} & =\frac{0.6666667-2.6666667\cdot c_{52}-2.3333333\cdot c_{42}}{5.3333333}=0.6246787,\\ c_{22} & =\frac{1-3.4\cdot c_{52}-0.2\cdot c_{42}-2.6\cdot c_{32}}{1.8}=0.9305913,\\ c_{12} & =\frac{0-(-2)\cdot c_{52}-(-1)\cdot c_{42}-(-3)\cdot c_{32}-1\cdot c_{22}}{5}=-0.151671. \end{align*} \]

Ejemplo (continuación)

Tercera columna de \(\mathbf{A}^{-1}\). Primero resolvemos: \[ \mathbf{L}\cdot\mathbf{y}_3 =\mathbf{e}_3,\ \Rightarrow \begin{bmatrix}1&0&0&0&0 \\0.2&1&0&0&0 \\0.2&-0.666667&1&0&0 \\0&0.555556&-0.270833&1&0 \\0.2&1.555556&-0.458333&0.211982&1 \\\end{bmatrix}\cdot \begin{bmatrix}y_{13}\\ y_{23}\\ y_{33}\\ y_{43}\\ y_{53}\end{bmatrix}=\begin{bmatrix}0\\0\\1\\0\\0\end{bmatrix}. \] Hallamos el vector \(\mathbf{y}_3\) usando la técnica de sustitución hacia adelante:

\[ \begin{align*} y_{13} & =0,\ y_{23}=0-0.2 y_{13}=0,\ y_{33}=1+0.666667 y_{23}-0.2 y_{13}=1,\\ y_{43} & =0+0.270833 y_{33}-0.555556 y_{23}=0.2708333,\\ y_{53} & =0-0.211982 y_{43}+0.458333 y_{33}-1.555556 y_{23}-0.2 y_{13}=0.4009217. \end{align*} \] El vector \(\mathbf{y}_3\) vale: \(\begin{bmatrix}0 \\0 \\1 \\0.270833 \\0.400922 \\\end{bmatrix}.\)

Ejemplo (continuación)

A continuación, calculamos la tercera columna \(\mathbf{c}_3\) de la matriz \(\mathbf{A}^{-1}\) resolviendo: \[ \mathbf{U}\cdot\mathbf{c}_3=\mathbf{y}_3,\quad \begin{bmatrix}5&1&-3&-1&-2 \\0&1.8&2.6&0.2&3.4 \\0&0&5.333333&2.333333&2.666667 \\0&0&0&4.520833&-2.166667 \\0&0&0&0&1.792627 \\\end{bmatrix}\cdot\begin{bmatrix}c_{13}\\c_{23}\\c_{33}\\c_{43}\\c_{53}\end{bmatrix}=\begin{bmatrix}0 \\0 \\1 \\0.270833 \\0.400922 \\\end{bmatrix}. \]

\[ \begin{align*} c_{53} & =\frac{0.400922}{1.792627}=0.22365,\\ c_{43} & =\frac{0.270833-(-2.1666667)\cdot c_{53}}{4.5208333}=0.1670951,\\ c_{33} & =\frac{1-2.6666667\cdot c_{53}-2.3333333\cdot c_{43}}{5.3333333}=0.0025707,\\ c_{23} & =\frac{0-3.4\cdot c_{53}-0.2\cdot c_{43}-2.6\cdot c_{33}}{1.8}=-0.4447301,\\ c_{13} & =\frac{0-(-2)\cdot c_{53}-(-1)\cdot c_{43}-(-3)\cdot c_{33}-1\cdot c_{23}}{5}=0.2133676. \end{align*} \]

Ejemplo (continuación)

Cuarta columna de \(\mathbf{A}^{-1}\). Primero resolvemos: \[ \mathbf{L}\cdot\mathbf{y}_4 =\mathbf{e}_4,\ \Rightarrow \begin{bmatrix}1&0&0&0&0 \\0.2&1&0&0&0 \\0.2&-0.666667&1&0&0 \\0&0.555556&-0.270833&1&0 \\0.2&1.555556&-0.458333&0.211982&1 \\\end{bmatrix}\cdot \begin{bmatrix}y_{14}\\ y_{24}\\ y_{34}\\ y_{44}\\ y_{54}\end{bmatrix}=\begin{bmatrix}0\\0\\0\\1\\0\end{bmatrix}. \] Hallamos el vector \(\mathbf{y}_4\) usando la técnica de sustitución hacia adelante:

\[ \begin{align*} y_{14} & =0,\ y_{24}=0-0.2 y_{14}=0,\ y_{34}=0+0.666667 y_{24}-0.2 y_{14}=0,\\ y_{44} & =1+0.270833 y_{34}-0.555556 y_{24}=1,\\ y_{54} & =0-0.211982 y_{44}+0.458333 y_{34}-1.555556 y_{24}-0.2 y_{14}=-0.2119816. \end{align*} \] El vector \(\mathbf{y}_4\) vale: \(\begin{bmatrix}0 \\0 \\0 \\1 \\-0.211982 \\\end{bmatrix}.\)

Ejemplo (continuación)

A continuación, calculamos la cuarta columna \(\mathbf{c}_4\) de la matriz \(\mathbf{A}^{-1}\) resolviendo: \[ \mathbf{U}\cdot\mathbf{c}_4=\mathbf{y}_4,\quad \begin{bmatrix}5&1&-3&-1&-2 \\0&1.8&2.6&0.2&3.4 \\0&0&5.333333&2.333333&2.666667 \\0&0&0&4.520833&-2.166667 \\0&0&0&0&1.792627 \\\end{bmatrix}\cdot\begin{bmatrix}c_{14}\\c_{24}\\c_{34}\\c_{44}\\c_{54}\end{bmatrix}=\begin{bmatrix}0 \\0 \\0 \\1 \\-0.211982 \\\end{bmatrix}. \]

\[ \begin{align*} c_{54} & =\frac{-0.211982}{1.792627}=-0.118252,\\ c_{44} & =\frac{1-(-2.1666667)\cdot c_{54}}{4.5208333}=0.1645244,\\ c_{34} & =\frac{0-2.6666667\cdot c_{54}-2.3333333\cdot c_{44}}{5.3333333}=-0.0128535,\\ c_{24} & =\frac{0-3.4\cdot c_{54}-0.2\cdot c_{44}-2.6\cdot c_{34}}{1.8}=0.2236504,\\ c_{14} & =\frac{0-(-2)\cdot c_{54}-(-1)\cdot c_{44}-(-3)\cdot c_{34}-1\cdot c_{24}}{5}=-0.066838. \end{align*} \]

Ejemplo (continuación)

Quinta columna de \(\mathbf{A}^{-1}\). Primero resolvemos: \[ \mathbf{L}\cdot\mathbf{y}_5 =\mathbf{e}_5,\ \Rightarrow \begin{bmatrix}1&0&0&0&0 \\0.2&1&0&0&0 \\0.2&-0.666667&1&0&0 \\0&0.555556&-0.270833&1&0 \\0.2&1.555556&-0.458333&0.211982&1 \\\end{bmatrix}\cdot \begin{bmatrix}y_{15}\\ y_{25}\\ y_{35}\\ y_{45}\\ y_{55}\end{bmatrix}=\begin{bmatrix}0\\0\\0\\0\\1\end{bmatrix}. \] Hallamos el vector \(\mathbf{y}_5\) usando la técnica de sustitución hacia adelante:

\[ \begin{align*} y_{15} & =0,\ y_{25}=0-0.2 y_{15}=0,\ y_{35}=0+0.666667 y_{25}-0.2 y_{15}=0,\\ y_{45} & =0+0.270833 y_{35}-0.555556 y_{25}=0,\\ y_{55} & =1-0.211982 y_{45}+0.458333 y_{35}-1.555556 y_{25}-0.2 y_{15}=1. \end{align*} \] El vector \(\mathbf{y}_5\) vale: \(\begin{bmatrix}0 \\0 \\0 \\0 \\1 \\\end{bmatrix}.\)

Ejemplo (continuación)

A continuación, calculamos la quinta columna \(\mathbf{c}_5\) de la matriz \(\mathbf{A}^{-1}\) resolviendo: \[ \mathbf{U}\cdot\mathbf{c}_5=\mathbf{y}_5,\quad \begin{bmatrix}5&1&-3&-1&-2 \\0&1.8&2.6&0.2&3.4 \\0&0&5.333333&2.333333&2.666667 \\0&0&0&4.520833&-2.166667 \\0&0&0&0&1.792627 \\\end{bmatrix}\cdot\begin{bmatrix}c_{15}\\c_{25}\\c_{35}\\c_{45}\\c_{55}\end{bmatrix}=\begin{bmatrix}0 \\0 \\0 \\0 \\1 \\\end{bmatrix}. \]

\[ \begin{align*} c_{55} & =\frac{1}{1.792627}=0.557841,\\ c_{45} & =\frac{0-(-2.1666667)\cdot c_{55}}{4.5208333}=0.2673522,\\ c_{35} & =\frac{0-2.6666667\cdot c_{55}-2.3333333\cdot c_{45}}{5.3333333}=-0.3958869,\\ c_{25} & =\frac{0-3.4\cdot c_{55}-0.2\cdot c_{45}-2.6\cdot c_{35}}{1.8}=-0.5115681,\\ c_{15} & =\frac{0-(-2)\cdot c_{55}-(-1)\cdot c_{45}-(-3)\cdot c_{35}-1\cdot c_{25}}{5}=0.1413882. \end{align*} \]

Ejemplo (continuación)

La matriz inversa \(\mathbf{A}^{-1}\) será: \[ \mathbf{A}^{-1}=\begin{bmatrix}0.159383&-0.151671&0.213368&-0.066838&0.141388 \\0.005141&0.930591&-0.44473&0.22365&-0.511568 \\-0.046272&0.624679&0.002571&-0.012853&-0.395887 \\-0.007712&-0.395887&0.167095&0.164524&0.267352 \\-0.025707&-0.652956&0.22365&-0.118252&0.557841 \\\end{bmatrix}. \]

Inversa de una matriz usando \(LU\)

Observación:

El vector \(\mathbf{y}_i\) es de la forma: \[ \mathbf{y}_i =\begin{bmatrix}0\\\vdots\\0\\\overbrace{1}^{i)}\\y_{i+1,i}\\\vdots\\y_{n,i} \end{bmatrix}. \] Por tanto, sólo hace falta hallar \(n-i\) componentes de dicho vector.

Ejemplo (continuación)

Hallemos la inversa de la matriz \(\mathbf{A}=\begin{bmatrix}5&1&-3&-1&-2 \\1&2&2&0&3 \\1&-1&3&2&0 \\0&1&0&4&-1 \\1&3&1&0&5 \\\end{bmatrix}\) usando la descomposición \(QR\).

La descomposición \(QR\) de la matriz \(\mathbf{A}\) es la siguiente: \[ \begin{align*} \mathbf{Q} & =\begin{bmatrix}0.944911&-0.167701&-0.278534&0.024779&-0.028702 \\0.188982&0.463645&0.41874&-0.206152&-0.729025 \\0.188982&-0.364997&0.853128&0.202194&0.249705 \\0&0.276214&-0.068225&0.949537&-0.132028 \\0.188982&0.739859&0.120802&-0.11994&0.622828 \\\end{bmatrix},\\ \mathbf{R} & =\begin{bmatrix}5.291503&1.70084&-1.70084&-0.566947&-0.377964 \\0&3.620379&1.075262&0.542564&5.149422 \\0&0&4.353269&1.71189&2.485526 \\0&0&0&4.177754&-2.217249 \\0&0&0&0&1.116499 \\\end{bmatrix}. \end{align*} \]

Ejemplo (continuación)

Primera columna \(\mathbf{c}_1\) de \(\mathbf{A}^{-1}\). Hemos de resolver el sistema siguiente \(\mathbf{R}\cdot \mathbf{c}_1=\mathbf{qt}_1\), donde recordemos que \(\mathbf{qt}_1\) sería la primera fila de la matriz \(\mathbf{Q}\): \[ \begin{bmatrix}5.291503&1.70084&-1.70084&-0.566947&-0.377964 \\0&3.620379&1.075262&0.542564&5.149422 \\0&0&4.353269&1.71189&2.485526 \\0&0&0&4.177754&-2.217249 \\0&0&0&0&1.116499 \\\end{bmatrix}\cdot\begin{bmatrix}c_{11}\\c_{21}\\c_{31}\\c_{41}\\c_{51}\end{bmatrix}=\begin{bmatrix}0.944911 \\-0.167701 \\-0.278534 \\0.024779 \\-0.028702 \\\end{bmatrix}. \]

Usamos la técnica de sustitución hacia atrás: \[ \begin{align*} c_{51} & =\frac{-0.028702}{1.116499}=-0.025707,\\ c_{41} & =\frac{0.024779-(-2.2172487)\cdot c_{51}}{4.1777542}=-0.0077121,\\ c_{31} & =\frac{-0.2785341-2.4855259\cdot c_{51}-1.7118895\cdot c_{41}}{4.3532693}=-0.0462725,\\ c_{21} & =\frac{-0.1677015-5.1494217\cdot c_{51}-0.5425636\cdot c_{41}-1.0752624\cdot c_{31}}{3.6203788}=0.0051414,\\ c_{11} & =\frac{0.9449112-(-0.3779645)\cdot c_{51}-(-0.5669467)\cdot c_{41}-(-1.7008401)\cdot c_{31}-1.7008401\cdot c_{21}}{5.2915026}=0.159383. \end{align*} \]

Ejemplo (continuación)

Segunda columna \(\mathbf{c}_2\) de \(\mathbf{A}^{-1}\). Hemos de resolver el sistema siguiente \(\mathbf{R}\cdot \mathbf{c}_2=\mathbf{qt}_2\), donde recordemos que \(\mathbf{qt}_2\) sería la segunda fila de la matriz \(\mathbf{Q}\): \[ \begin{bmatrix}5.291503&1.70084&-1.70084&-0.566947&-0.377964 \\0&3.620379&1.075262&0.542564&5.149422 \\0&0&4.353269&1.71189&2.485526 \\0&0&0&4.177754&-2.217249 \\0&0&0&0&1.116499 \\\end{bmatrix}\cdot\begin{bmatrix}c_{12}\\c_{22}\\c_{32}\\c_{42}\\c_{52}\end{bmatrix}=\begin{bmatrix}0.188982 \\0.463645 \\0.41874 \\-0.206152 \\-0.729025 \\\end{bmatrix}. \]

Usamos la técnica de sustitución hacia atrás: \[ \begin{align*} c_{52} & =\frac{-0.729025}{1.116499}=-0.652956,\\ c_{42} & =\frac{-0.206152-(-2.2172487)\cdot c_{52}}{4.1777542}=-0.3958869,\\ c_{32} & =\frac{0.4187401-2.4855259\cdot c_{52}-1.7118895\cdot c_{42}}{4.3532693}=0.6246787,\\ c_{22} & =\frac{0.4636452-5.1494217\cdot c_{52}-0.5425636\cdot c_{42}-1.0752624\cdot c_{32}}{3.6203788}=0.9305913,\\ c_{12} & =\frac{0.1889822-(-0.3779645)\cdot c_{52}-(-0.5669467)\cdot c_{42}-(-1.7008401)\cdot c_{32}-1.7008401\cdot c_{22}}{5.2915026}=-0.151671. \end{align*} \]

Ejemplo (continuación)

Tercera columna \(\mathbf{c}_3\) de \(\mathbf{A}^{-1}\). Hemos de resolver el sistema siguiente \(\mathbf{R}\cdot \mathbf{c}_3=\mathbf{qt}_3\), donde recordemos que \(\mathbf{qt}_3\) sería la tercera fila de la matriz \(\mathbf{Q}\): \[ \begin{bmatrix}5.291503&1.70084&-1.70084&-0.566947&-0.377964 \\0&3.620379&1.075262&0.542564&5.149422 \\0&0&4.353269&1.71189&2.485526 \\0&0&0&4.177754&-2.217249 \\0&0&0&0&1.116499 \\\end{bmatrix}\cdot\begin{bmatrix}c_{13}\\c_{23}\\c_{33}\\c_{43}\\c_{53}\end{bmatrix}=\begin{bmatrix}0.188982 \\-0.364997 \\0.853128 \\0.202194 \\0.249705 \\\end{bmatrix}. \]

Usamos la técnica de sustitución hacia atrás: \[ \begin{align*} c_{53} & =\frac{0.249705}{1.116499}=0.22365,\\ c_{43} & =\frac{0.202194-(-2.2172487)\cdot c_{53}}{4.1777542}=0.1670951,\\ c_{33} & =\frac{0.8531281-2.4855259\cdot c_{53}-1.7118895\cdot c_{43}}{4.3532693}=0.0025707,\\ c_{23} & =\frac{-0.3649973-5.1494217\cdot c_{53}-0.5425636\cdot c_{43}-1.0752624\cdot c_{33}}{3.6203788}=-0.4447301,\\ c_{13} & =\frac{0.1889822-(-0.3779645)\cdot c_{53}-(-0.5669467)\cdot c_{43}-(-1.7008401)\cdot c_{33}-1.7008401\cdot c_{23}}{5.2915026}=0.2133676. \end{align*} \]

Ejemplo (continuación)

Cuarta columna \(\mathbf{c}_4\) de \(\mathbf{A}^{-1}\). Hemos de resolver el sistema siguiente \(\mathbf{R}\cdot \mathbf{c}_4=\mathbf{qt}_4\), donde recordemos que \(\mathbf{qt}_4\) sería la cuarta fila de la matriz \(\mathbf{Q}\): \[ \begin{bmatrix}5.291503&1.70084&-1.70084&-0.566947&-0.377964 \\0&3.620379&1.075262&0.542564&5.149422 \\0&0&4.353269&1.71189&2.485526 \\0&0&0&4.177754&-2.217249 \\0&0&0&0&1.116499 \\\end{bmatrix}\cdot\begin{bmatrix}c_{14}\\c_{24}\\c_{34}\\c_{44}\\c_{54}\end{bmatrix}=\begin{bmatrix}0 \\0.276214 \\-0.068225 \\0.949537 \\-0.132028 \\\end{bmatrix}. \]

Usamos la técnica de sustitución hacia atrás: \[ \begin{align*} c_{54} & =\frac{-0.132028}{1.116499}=-0.118252,\\ c_{44} & =\frac{0.949537-(-2.2172487)\cdot c_{54}}{4.1777542}=0.1645244,\\ c_{34} & =\frac{-0.0682252-2.4855259\cdot c_{54}-1.7118895\cdot c_{44}}{4.3532693}=-0.0128535,\\ c_{24} & =\frac{0.2762142-5.1494217\cdot c_{54}-0.5425636\cdot c_{44}-1.0752624\cdot c_{34}}{3.6203788}=0.2236504,\\ c_{14} & =\frac{0-(-0.3779645)\cdot c_{54}-(-0.5669467)\cdot c_{44}-(-1.7008401)\cdot c_{34}-1.7008401\cdot c_{24}}{5.2915026}=-0.066838. \end{align*} \]

Ejemplo (continuación)

Quinta columna \(\mathbf{c}_5\) de \(\mathbf{A}^{-1}\). Hemos de resolver el sistema siguiente \(\mathbf{R}\cdot \mathbf{c}_5=\mathbf{qt}_5\), donde recordemos que \(\mathbf{qt}_5\) sería la quinta fila de la matriz \(\mathbf{Q}\): \[ \begin{bmatrix}5.291503&1.70084&-1.70084&-0.566947&-0.377964 \\0&3.620379&1.075262&0.542564&5.149422 \\0&0&4.353269&1.71189&2.485526 \\0&0&0&4.177754&-2.217249 \\0&0&0&0&1.116499 \\\end{bmatrix}\cdot\begin{bmatrix}c_{15}\\c_{25}\\c_{35}\\c_{45}\\c_{55}\end{bmatrix}=\begin{bmatrix}0.188982 \\0.739859 \\0.120802 \\-0.11994 \\0.622828 \\\end{bmatrix}. \]

Usamos la técnica de sustitución hacia atrás: \[ \begin{align*} c_{55} & =\frac{0.622828}{1.116499}=0.557841,\\ c_{45} & =\frac{-0.11994-(-2.2172487)\cdot c_{55}}{4.1777542}=0.2673522,\\ c_{35} & =\frac{0.1208024-2.4855259\cdot c_{55}-1.7118895\cdot c_{45}}{4.3532693}=-0.3958869,\\ c_{25} & =\frac{0.7398594-5.1494217\cdot c_{55}-0.5425636\cdot c_{45}-1.0752624\cdot c_{35}}{3.6203788}=-0.5115681,\\ c_{15} & =\frac{0.1889822-(-0.3779645)\cdot c_{55}-(-0.5669467)\cdot c_{45}-(-1.7008401)\cdot c_{35}-1.7008401\cdot c_{25}}{5.2915026}=0.1413882. \end{align*} \]

Ejemplo (continuación)

La inversa está calculada en el siguiente enlace:

Determinante de una matriz \(\mathbf{A}\)

Para calcular el determinante de una matriz \(\mathbf{A}\) cuadrada \(n\times n\), podemos aprovechar la descomposición \(LU\) de dicha matriz y afirmar que \(\mathrm{det}(\mathbf{A})=\mathrm{det}(\mathbf{L})\cdot\mathrm{det}(\mathbf{U})\).

Ahora bien, como la matriz \(\mathbf{L}\) es triangular inferior con \(1\) en la diagonal, \(\mathrm{det}(\mathbf{L})=1\) y como la matriz \(\mathbf{U}=(u_{ij})_{i,j=1,\ldots,n}\) es triangular superior, \(\mathrm{det}(\mathbf{U})=u_{11}\cdots u_{nn}\).

En resumen, \(\mathrm{det}(\mathrm{A})=u_{11}\cdots u_{nn}\).

Observación.

Si quisiéramos usar la descomposición \(QR\) de la matriz \(\mathbf{A}\), sólo podemos conocer \(|\mathrm{det}(\mathbf{A})|\) ya que \(\mathrm{det}(\mathbf{A})=\mathrm{det}(\mathbf{Q})\cdot\mathrm{det}(\mathbf{R})\).

La matriz \(\mathbf{Q}\), al ser ortogonal, verifica que \(\mathrm{det}(\mathbf{Q})=\pm 1\) y la matriz \(\mathbf{R}=(r_{ij})_{i,j=1,\ldots,n}\) al ser triangular superior, verifica que \(\mathrm{det}(\mathbf{R})=r_{11}\cdots r_{nn}\).

En resumen, \(\mathrm{det}(\mathrm{A})=\pm r_{11}\cdots r_{nn}\). Por tanto, sólo podemos conocer \(|\mathrm{det}(\mathbf{A})|\).

Ejemplo

Vamos a hallar el determinante de la matriz vista anteriormente \[ \mathbf{A}=\begin{bmatrix}5&1&-3&-1&-2 \\1&2&2&0&3 \\1&-1&3&2&0 \\0&1&0&4&-1 \\1&3&1&0&5 \\\end{bmatrix}, \] usando la descomposición \(LU\). Recordemos que dicha descomposición era la siguiente: \[ \mathbf{L}=\begin{bmatrix}1&0&0&0&0 \\0.2&1&0&0&0 \\0.2&-0.666667&1&0&0 \\0&0.555556&-0.270833&1&0 \\0.2&1.555556&-0.458333&0.211982&1 \\\end{bmatrix},\quad\mathbf{U}=\begin{bmatrix}5&1&-3&-1&-2 \\0&1.8&2.6&0.2&3.4 \\0&0&5.333333&2.333333&2.666667 \\0&0&0&4.520833&-2.166667 \\0&0&0&0&1.792627 \\\end{bmatrix}. \]

El determinante de la matriz \(\mathbf{A}\) vale: \[ \mathrm{det}(\mathbf{A})=u_{11}\cdot u_{22}\cdot u_{33}\cdot u_{44}\cdot u_{55}=5\cdot 1.8\cdot 5.333333\cdot 4.520833\cdot 1.792627=389. \]

Análisis del error

Introducción

Cuando resolvemos un sistema lineal \(\mathbf{A}\mathbf{x}=\mathbf{b}\), la matriz del sistema \(\mathbf{A}\) y el vector de términos independientes \(\mathbf{b}\) pueden estar afectados por un error en los datos.

Entonces, en realidad, resolvemos el sistema \[ (\mathbf{A}+\delta(\mathbf{A}))(\mathbf{x}+\delta(\mathbf{x}))=\mathbf{b}+\delta(\mathbf{b}), \] donde \(\mathbf{A}+\delta(\mathbf{A})\) sería la matriz del sistema perturbada o con error y \(\mathbf{b}+\delta(\mathbf{b})\) sería el vector de términos independientes perturbado o con error.

En la expresión anterior, indicamos por \(\mathbf{x}\) la solución exacta o sin errores del sistema lineal y \(\mathbf{x}+\delta(\mathbf{x})\) la solución del sistema lineal perturbado o con error.

La cuestión que nos planteamos es cómo afectan dichos errores a la solución del sistema \(\mathbf{x}+\delta(\mathbf{x})\) o, dicho de otra forma, ¿existe alguna forma de controlar el error de la solución del sistema, \(\delta(\mathbf{x})\) en función de alguna característica de la matriz del sistema \(\mathbf{A}\) y el vector de términos independientes \(\mathbf{b}\)?

Error en la matriz del sistema

Supongamos de cara a simplificar nuestro problema que sólo existe error en los datos de la matriz del sistema \(\mathbf{A}\). En este caso, resolvemos el sistema siguiente: \[ (\mathbf{A}+\delta(\mathbf{A}))(\mathbf{x}+\delta(\mathbf{x}))=\mathbf{b}. \] Si desarrollamos la expresión anterior, obtenemos: \[ \mathbf{A}\mathbf{x}+\mathbf{A}\delta(\mathbf{x})+\mathbf{\delta(A)}\mathbf{x}+\delta(\mathbf{A})\delta(\mathbf{x})=\mathbf{b}. \] Usando que \(\mathbf{A}\mathbf{x}=\mathbf{b}\) y despreciando el término \(\delta(\mathbf{A})\delta(\mathbf{x})\) ya que es el producto de dos errores y su valor será muy pequeño en comparación con los demás, nos queda: \[ \mathbf{A}\delta(\mathbf{x})+\mathbf{\delta(A)}\mathbf{x}\approx \mathbf{0},\ \Rightarrow \mathbf{A}\delta(\mathbf{x})\approx -\mathbf{\delta(A)}\mathbf{x},\ \Rightarrow \delta(\mathbf{x})\approx -\mathbf{A}^{-1}\mathbf{\delta(A)}\mathbf{x}. \]

Error en la matriz del sistema

Tomando una norma vectorial \(\|\cdot\|\) y su normal matricial correspondiente, obtenemos que: \[ \begin{align*} \|\delta(\mathbf{x})\| & \approx \|\mathbf{A}^{-1}\delta(\mathbf{A})\mathbf{x}\|\leq \|\mathbf{A}^{-1}\|\cdot\|\delta(\mathbf{A})\|\cdot\|\mathbf{x}\|\\ & =\|\mathbf{A}^{-1}\|\cdot\|\mathbf{A}\|\frac{\|\delta(\mathbf{A})\|}{\|\mathbf{A}\|}\cdot\|\mathbf{x}\|. \end{align*} \] Es decir, \[ \frac{\|\delta(\mathbf{x})\|}{\|\mathbf{x}\|}\leq \|\mathbf{A}^{-1}\|\cdot\|\mathbf{A}\|\cdot \frac{\|\delta(\mathbf{A})\|}{\|\mathbf{A}\|}. \] La expresión anterior nos dice que el error relativo cometido en la solución del sistema lineal está acotado por la cantidad \(\|\mathbf{A}^{-1}\|\cdot\|\mathbf{A}\|\) por el error relativo de la matriz del sistema.

Error en la matriz del sistema

A la cantidad \(\mu(\mathbf{A}):=\|\mathbf{A}^{-1}\|\cdot\|\mathbf{A}\|\) se le llama número de condición de la matriz \(\mathbf{A}\) y, como hemos visto, juega un papel fundamental a la hora de controlar el error relativo cometido en la solución del sistema.

Es decir, si el número de condición de una matriz es grande, el sistema lineal correspondiente será inestable en el sentido de que una pequeña perturbación en la matriz del sistema dará lugar a una solución con un error relativo grande.

Error en el vector de términos independientes

Supongamos ahora que sólo existe error en los datos del vector de términos independientes \(\mathbf{b}\). En este caso, resolvemos el sistema siguiente: \[ \mathbf{A}(\mathbf{x}+\delta(\mathbf{x}))=\mathbf{b}+\delta(\mathbf{b}). \] Si desarrollamos la expresión anterior, obtenemos: \[ \mathbf{A}\mathbf{x}+\mathbf{A}\delta(\mathbf{x})=\mathbf{b}+\delta(\mathbf{b}). \] Usando que \(\mathbf{A}\mathbf{x}=\mathbf{b}\), nos queda: \[ \mathbf{A}\delta(\mathbf{x})=\delta(\mathbf{b}), \Rightarrow \delta(\mathbf{x})= \mathbf{A}^{-1}\delta(\mathbf{b}). \]

Error en el vector de términos independientes

Tomando igual que antes una norma vectorial \(\|\cdot\|\) y su normal matricial correspondiente, obtenemos que: \[ \frac{\|\delta(\mathbf{x})\|}{\|\mathbf{x}\|}=\frac{\|\mathbf{A}^{-1}\delta(\mathbf{b})\|}{\|\mathbf{A}^{-1}\mathbf{b}\|}. \] A continuación veamos que \(\frac{1}{\|\mathbf{A}^{-1}b\|}\leq\frac{\|\mathbf{A}\|}{\|\mathbf{b}\|}\): \[ \begin{align*} \|\mathbf{b}\| & =\|\mathbf{A}\mathbf{A}^{-1}\mathbf{b}\|\leq\|\mathbf{A}\|\cdot\|\mathbf{A}^{-1}\mathbf{b}\|,\ \Rightarrow \frac{\|\mathbf{b}\|}{\|\mathbf{A}\|}\leq \|\mathbf{A}^{-1}\mathbf{b}\|,\ \Rightarrow\\ \frac{\|\mathbf{A}\|}{\|\mathbf{b}\|} & \geq\frac{1}{\|\mathbf{A}^{-1}\mathbf{b}\|}. \end{align*} \]

Error en el vector de términos independientes

Por tanto, \[ \frac{\|\delta(\mathbf{x})\|}{\|\mathbf{x}\|}=\frac{\|\mathbf{A}^{-1}\delta(\mathbf{b})\|}{\|\mathbf{A}^{-1}\mathbf{b}\|}\leq \|\mathbf{A}^{-1}\|\cdot\|\mathbf{A}\|\cdot\frac{\|\delta(\mathbf{b})\|}{\|\mathbf{b}\|}=\mu(\mathbf{A})\cdot\frac{\|\delta(\mathbf{b})\|}{\|\mathbf{b}\|}. \] Dicho en otras palabras, el error relativo de la solución del sistema es menor que el número de condición de la matriz del sistema por el error relativo del vector de términos independientes.

En resumen, si el número de condición de una matriz es grande, el sistema lineal correspondiente será inestable en el sentido de que una pequeña perturbación en el vector de términos independientes dará lugar a una solución con un error relativo grande.

Error en ambos

Supongamos ahora que tenemos error en los datos de la matriz del sistema \(\mathbf{A}\) y en el vector de términos independientes \(\mathbf{b}\). En este caso, resolvemos el sistema siguiente: \[ (\mathbf{A}+\delta(\mathbf{A}))(\mathbf{x}+\delta(\mathbf{x}))=\mathbf{b}+\delta(\mathbf{b}). \] En este caso, puede demostrarse la desigualdad siguiente fijada una norma vectorial junto con la correspondiente norma matricial: \[ \frac{\|\delta(\mathbf{x})\|}{\|\mathbf{x}\|}\leq \mu(\mathbf{A})\left(\frac{\|\delta(\mathbf{b})\|}{\|\mathbf{b}\|}+\frac{\|\delta(\mathbf{A})\|}{\|\mathbf{A}\|}\cdot \frac{\|\mathbf{x}+\delta(\mathbf{x})\|}{\|\mathbf{x}\|}\right). \]

Error en ambos

Observaciones:

El número de condición de la matriz del sistema es fundamental para estudiar la estabilidad de las soluciones, es decir, si el error relativo de la solución del sistema es grande o no.
En la práctica \(\|\mathbf{x}+\delta(\mathbf{x})\|\approx \|\mathbf{x}\|.\)
Si \(\mathbf{U}\) es una matriz ortogonal, y la norma considerada es la norma euclídea, entonces en general \(\mu_2(\mathbf{U}\mathbf{A})=\|\mathbf{U}\mathbf{A}\|_2\cdot\|(\mathbf{U}\mathbf{A})^{-1}\|_2=\mu_2(\mathbf{A})\) ya que:

Error en ambos

\[ \begin{align*} \mu_2(\mathbf{U}\mathbf{A}) & =\|\mathbf{U}\mathbf{A}\|_2\cdot\|(\mathbf{U}\mathbf{A})^{-1}\|_2\\ &=\sqrt{\rho((\mathbf{U}\mathbf{A})^\top\mathbf{U}\mathbf{A})}\cdot \sqrt{\rho(((\mathbf{U}\mathbf{A})^{-1})^\top(\mathbf{U}\mathbf{A})^{-1})}\\ & = \sqrt{\rho(\mathbf{A}^\top\mathbf{U}^\top\mathbf{U}\mathbf{A})}\cdot\sqrt{\rho((\mathbf{A}^{-1}\mathbf{U}^{-1})^\top \mathbf{A}^{-1}\mathbf{U}^{-1})}\\ & = \sqrt{\rho(\mathbf{A}^\top\mathbf{A})}\cdot\sqrt{\rho((\mathbf{A}^{-1}\mathbf{U}^\top)^\top \mathbf{A}^{-1}\mathbf{U}^\top)}\\ & = \|\mathbf{A}\|_2\cdot\sqrt{\rho(\mathbf{U}(\mathbf{A}^{-1})^\top\mathbf{A}^{-1}\mathbf{U}^\top)}= \|\mathbf{A}\|_2\cdot\sqrt{\rho((\mathbf{A}^{-1})^\top\mathbf{A}^{-1})}\\ & = \|\mathbf{A}\|_2\cdot\|\mathbf{A}^{-1}\|_2=\mu_2(\mathbf{A}), \end{align*} \] donde \(\rho(\cdot)\) representa el radio espectral de una matriz y hemos usado que \(\mathbf{U}^{-1}=\mathbf{U}^\top\).

En el penúltimo paso hemos utilizado que los valores propios de la matriz \(\mathbf{U}(\mathbf{A}^{-1})^\top\mathbf{A}^{-1}\mathbf{U}^\top\) son los mismos que los valores propios de la matriz \(\mathbf{A}^{-1})^\top\mathbf{A}^{-1}\) ya que la primera representa un cambio de base respecto de la segunda matriz. Por tanto, tendrán el mismo radio espectral, como vimos en el capítulo de Preliminares.

Error en ambos

Si resolvemos el sistema lineal \(\mathbf{A}\mathbf{x}=\mathbf{b}\) usando la descomposición \(QR\), tenemos que, usando la propiedad anterior, \(\mu_2(\mathbf{R})=\mu_2(\mathbf{Q}^\top\mathbf{A})=\mu_2(\mathbf{A})\) al ser la matriz \(\mathbf{Q}^\top\) ortogonal.

Entonces, cuando resolvemos el sistema triangular superior \(\mathbf{R}\mathbf{x}=\mathbf{Q}^\top\mathbf{b}\), la propagación de los errores de la solución \(\mathbf{x}\) no empeoran con respecto a la resolución del sistema original \(\mathbf{A}\mathbf{x}=\mathbf{b}\) al tener las dos matrices del sistema el mismo número de condición.

Ejemplo

Consideremos el sistema de ecuaciones siguiente visto anteriormente:

\[ \left. \begin{align*} 0.010534x_1+0.02x_2= & -1,\\ x_1+2x_2= & 3. \end{align*} \right\} \] La matriz del sistema es \(\mathbf{A}=\begin{bmatrix}0.010534&0.02 \\1&2 \\\end{bmatrix}\) cuya norma euclídea vale \(\|\mathbf{A}\|_2=2.2361822\).

La inversa de la matriz anterior vale \(\mathbf{A}^{-1}=\begin{bmatrix}1872.659176&-18.726592 \\-936.329588&9.863296 \\\end{bmatrix}\) cuya norma euclídea vale \(\|\mathbf{A}^{-1}\|_2=2093.803538.\)

El número de condición de la matriz del sistema vale: \[ \mu_2(\mathbf{A})=\|\mathbf{A}\|_2\cdot \|\mathbf{A}^{-1}\|_2=2.2361822\cdot2093.803538=4682.126158. \] El número de condición anterior es grande. Por tanto, la solución del sistema será inestable y sensible a los errores en la matriz del sistema y al vector de términos independientes.