versión pdf. - Universidad de Alcalá

Tema02: Sistemas de Ecuaciones.
Fernando San Segundo
Curso 2014-2015. Cálculo Numérico y Estadística. Grado en Química. Universidad de Alcalá.
Contents
Introducción: sistemas de ecuaciones. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Eliminación Gaussiana en sistemas de ecuaciones lineales.
1
2
Sistemas incompatibles e indeterminados. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
Matrices escalonadas, eliminación Gaussiana y rango. . . . . . . . . . . . . . . . . . . . . . . . . .
5
Teorema de Rouche-Frobenius. Matrices no singulares. . . . . . . . . . . . . . . . . . . . . . . . . .
8
Inversa de una matriz cuadrada. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
Determinantes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11
Complejidad computacional del método de eliminación. Número de operaciones. . . . . . . . . . .
12
Factorización de matrices.
13
Factorización LU. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
Diagonalización. Autovalores y autovectores. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
21
Descomposición en valores singulares (SVD). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
27
Factorización QR. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
32
Sistemas de ecuaciones no lineales.
33
Sistemas de ecuaciones polinómicas. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
34
Métodos numéricos para sistemas de ecuaciones no lineales. . . . . . . . . . . . . . . . . . . . . . .
35
Bibliografía.
36
Advertencia preliminar.
En este tema vamos a hacer uso de las matrices de R y de las operaciones con ellas. Si no recuerdas como se
llevan a cabo esas operaciones, te recomiendo que releas los tutoriales pertinentes del curso de Estadística,
muy especialmente la Sección 2 del Tutorial 3.
Introducción: sistemas de ecuaciones.
En el tema anterior hemos visto algunos métodos para localizar las soluciones de una ecuación de la forma
f (x) = 0,
1
donde f es una función de una variable. En este tema vamos a generalizar estas ideas para estudiar las
soluciones de sistemas de ecuaciones. Un sistema puede ser algo tan sencillo como
(
x − 3y = 1
3x − 5y = 7
o algo más complicado como este sistema:

(1 − x1 ) + 100 · (−x1 x2 ) = 0



(2 − x ) + 100 · (−x x − x x ) = 0
2
1 2
3 2

−x
+
100
·
(x
x
−
x
x
)
=
0
3
1
2
3
2



−x4 + 100 · (x3 x2 ) = 0
que aparece al estudiar el equilibrio en un reactor químico (ver pág. 71 de Numerical Methods for Chemical
Engineering, K.J. Beers, Cambridge University Press 2007).
En este tema nos vamos a centrar en sistemas como el primero, aunque al final diremos algo sobre esos otros
sistemas más complicados.
Eliminación Gaussiana en sistemas de ecuaciones lineales.
Un sistema de m ecuaciones lineales con n incógnitas es un sistema que se puede escribir en la forma:

a11 x1 + a12 x2 + · · · + a1n x2 = b1




a21 x1 + a22 x2 + · · · + a2n xn = b2
..


.



am1 x1 + am2 x2 + · · · + amn xn = bn
o en la forma matricial equivalente:
A·X =B
donde A es la matriz de coeficientes de dimensiones m × n,



A=

a11
a21
..
.
a12
a22
..
.
···
···
..
.
a1n
a2n
..
.
am1
am2
···
amn



,

B es la matriz columna de términos independientes, con dimensiones m × 1:



B=

b1
b2
..
.





bm
y X es la matriz columna de incógnitas, de dimensiones n × 1:
2

x1
x2
..
.


X=






xn
Como ilustración, el sistema de nuestro primer ejemplo es lineal, y su forma matricial es:
1 −3
x
1
=
.
3 −5
y
7
A menudo, para representar un sistema de ecuaciones

a11 · · ·
 a21 · · ·

 ..
..
 .
.
am1
se usa también la notación de matrices ampliadas:

a1n b1
a2n b2 

.. 
..
. 
.
···
amn
bm
En esta notación incluimos en la misma matriz los coeficientes aij y, como última columna añadimos los
términos independientes bi , separando estos últimos con una línea vertical en la notación que usamos aquí.
Sistemas escalonados, eliminación de variables, transformaciones elementales y sistemas equivalentes.
El siguiente sistema es especialmente fácil de resolver:

 2x1 + 3x2 + x3 = 7
−x2 + 5x3 = −7

x3 = −1
porque es un sistema escalonado: el valor de x3 se obtiene de la última ecuación, a continuación se sustituye
ese valor en la segunda ecuación y se despeja x2 y con esos valores despejamos x1 de la primera ecuación.
El siguiente sistema no es escalonado:
2x1 + x2 + 4x3 = 1
2x1 + 2x2 + 7x3 = 0

x1 − x2 − 3x3 = 1


Pero se puede transformar en uno escalonado. Empezamos en la primera columna restando la primera fila de
la segunda. Representamos así esta transformación: f2 − f1 → f2 . También reemplazamos la tercera fila con
la transformación: −2f3 + f1 → f3 .

 2x1 + x2 + 4x3 = 1
x2 + 3x3 = −1

3x2 + 10x3 = −1
Con esto hemos eliminado x1 de las dos últimas ecuaciones. Ahora hacemos la transformación f3 − 3f2 → f3 :

 2x1 + x2 + 4x3 = 1
x2 + 3x3 = −1

x3 = 2
Y así hemos obtenido un sistema equivalente. ¿Pero qué tienen que ver las soluciones de este sistema con las
del sistema original?
• Las transformaciones elementales de un sistema son:
3
– El intercambio de posiciones de dos ecuaciones del sistema: fi ↔ fj .
– La sustitución de una ecuación por una combinación lineal de ecuaciones, en la que esa ecuación
tiene un coeficiente no nulo:
a1 f1 + · · · + ai fi + · · · + am fm → fi
siempre que
ai 6= 0
• Dos sistemas son equivalentes si tienen exactamente el mismo conjunto de soluciones.
• Teoremas:
– Los sistemas conectados por transformaciones elementales son siempre equivalentes.
– Cualquier sistema se puede transformar en un sistema escalonado mediante transformaciones
elementales.
La Eliminación Gaussiana, que veremos más abajo, es un algoritmo que lleva a cabo de manera sistemática
esa transformación de un sistema a la forma escalonada. Lo más habitual es llevar a cabo las transformaciones
elementales directamente sobre la expresión matricial del sistema.
Sistemas incompatibles e indeterminados.
Los siguientes ejemplos muestran algunas situaciones típicas con las que nos encontramos al tratar con
sistemas de ecuaciones lineales.
• Aplicando transformaciones elementales para escalonar el sistema

 x1 + 2x2 − 3x3 = −1
3x1 − 2x2 + 2x3 = 7

5x1 + 3x2 − 4x3 = 2
obtenemos este sistema equivalente:

 x1 + 2x2 − 3x3 = −1
−7x2 + 11x3 = 10

0x3 = 21
La última ecuación es imposible de satisfacer, sea cual sea el valor de x3 . Así que este sistema (y por lo
tanto el sistema incial) no tiene soluciones: es un sistema incompatible.
• Aplicando transformaciones elementales para escalonar el sistema

x1 + 2x2 − 3x3 = 6

2x1 − x2 + 4x3 = 2

4x1 + 3x2 − 2x3 = 14
obtenemos este sistema equivalente:

 x1 + 2x2 − 3x3 = 6
−5x2 + 10x3 = −10

0x3 = 0
A diferencia de lo que sucedía en el ejemplo anterior, cualquier valor de x3 satisface la última ecuación.
Una vez elegido ese valor podemos usarlo como antes, sustituyendo una por una en las ecuaciones
precedentes, para obtener los valores de x2 y x1 . Pero como el valor de x3 es arbitrario, eso significa
que el sistema tiene infinitas soluciones: es un sistema indeterminado.
• Los ejemplos más sencillos que hemos visto antes corresponden a sistemas con solución única. Esos
sistemas se llaman sistemas compatibles determinados.
4
Todos los sistemas de ecuaciones lineales, con independencia del número de ecuaciones y variables, pertenecen
a una de esas tres categorías. Podemos, no obstante, hacer algunas observaciones:
• Un sistema con menos ecuaciones que incógnitas no puede tener solución única. Pero puede ser
incompatible y no tener soluciones.
• Un sistema con más ecuaciones que incógnitas es, a menudo, un sistema incompatible. Pero puede que
al escalonarlo algunas de esas ecuaciones extra resulten ser ecuaciones superfluas, y el sistema resulte
compatible.
Matrices escalonadas, eliminación Gaussiana y rango.
Una matriz en forma escalonada por filas (row echelon form) se caracteriza por esta propiedad:
Si el primer elemento no nulo de la fila i está en la columna k, entonces en la siguiente fila, la fila i + 1,
todos los elementos de las columnas 1 a k (inclusive) son nulos. Y si todos los elementos de la fila i son 0,
entonces todas las filas a partir de la i + 1 son también filas nulas.
Esta propiedad nos garantiza que el número de elementos nulos aumenta al pasar de cada fila a la siguiente,
hasta llegar a la posible presencia de una cierta cantidad filas nulas en la parte inferior de la matriz.
Eliminación Gaussiana y eliminación de Gauss-Jordan.
Para obtener una matriz escalonada equivalente a la matriz dada podemos usar el Algoritmo de Eliminación
Gaussiana. La idea general es avanzar columna a columna, de izquierda a derecha, y descendiendo un
peldaño cada vez. Usamos el primer elemento no nulo de cada fila para hacer cero los elementos que aparecen
debajo en esa misma columna (si es preciso, intercambiamos filas para conseguir un elemento no nulo en la
posición adecuada). Puedes ver una descripción más detallada y ejemplos en este enlace de la Wikipedia:
http://es.wikipedia.org/wiki/Eliminación_de_Gauss-Jordan#Algoritmo_de_eliminación_de_Gauss-Jordan
y si quieres referencias más completas puedes consultar alguno de los libros que aparecen en las referencias al
final de estas notas:
Aquí no nos entretendremos en los detalles porque vamos a usar el ordenador para aplicar el metodo (aunque
te recomendamos que practiques resolviendo a mano algunos ejemplos para asimilar esos detalles).
Rango.
Sea A una matriz cualquiera y sea A0 una matriz escalonada equivalente (que se obtiene de A mediante
transformaciones elementales). El rango de A es el número de filas con elementos no nulos que contiene A0 .
Por ejemplo, dado el sistema:

2x1 +2x2
−x3
+x4 = 4



4x1 +3x2
−x3 +2x4 = 6
8x1 +5x2 −3x3 +4x4 = 12



2x1
−x3
+x4 = 3
Su forma matricial es:


2 2 −1 1 4
 4 3 −1 2 6 


 8 5 −3 4 12 
2 0 −1 1 3
y una forma escalonada equivalente es:


2
2 −1 1
4
 0 −1
1 0 −2 


 0
0 −2 0
2 
0
0
0 0 −1
5
Así que el rango de la matriz es 4 porque hay cuatro peldaños no nulos. Fíjate en que este sistema es
compatible determinado. En cambio para el sistema indeterminado que hemos visto antes:

 x1 + 2x2 − 3x3 = 6
2x1 − x2 + 4x3 = 2

4x1 + 3x2 − 2x3 = 14
obtenemos esta matriz ampliada escalonada:

1
 0
0

6
2 −3
−5 10 −10 
0
0
0
Fíjate en que en este caso tanto la matriz ampliada como la matriz de coeficientes tienen rango 2.
En el ejemplo de sistema incompatible:

 x1 + 2x2 − 3x3 = −1
3x1 − 2x2 + 2x3 = 7

5x1 + 3x2 − 4x3 = 2
la matriz ampliada escalonada es:
1
 0
0


2 −3 −1
−7 11 10 
0
0 21
Y en este caso la observación esencial es que el rango de la matriz ampliada (que es 3) es mayor que el rango
de la matriz de coeficientes (que es 2).
Eliminación de Gauss-Jordan.
Pero sí queremos comentar un aspecto que hasta ahora no hemos destacado en los ejemplos. El procedimiento
de eliminación Gaussiana, tal como lo hemos descrito y usado, procede de arriba a abajo en la matriz, y se
detiene al alcanzar la última fila. Pero es posible hacer un procedimiento de eliminación más completo, la
eliminación de Gauss-Jordan, en el que al llegar a la última fila de la matriz “volvemos hacia arriba” usndo
el primer elemento no nulo de cada peldaño para eliminar todos los elementos que tiene por encima.
Por ejemplo, en el sistema escalonado

 2x1 + x2 + 4x3 = 1
x2 + 3x3 = −1

x3 = 2
que hemos analizado antes, la matriz ampliada escalonada es:


2 1 4
1
 0 1 3 −1 
0 0 1
2
Y para completar “hacaia arriba” la eliminación ahora podemos hacer transformaciones adicionales. Empezamos trabajando en la tercer columna, con: f1 − 4f3 → f1 seguido de f2 − 3f3 → f2 :


2 1 0 −7
 0 1 0 −7 
0 0 1
2
Ahora la tercera columna sólo contiene un elemento no nulo. Seguimos en la segunda columna con f1 −f2 → f1 :


2 0 0
0
 0 1 0 −7 
0 0 1
2
Y ahora es inmediato que la solución única del sistema es:
x1 = 0, x2 = −7, x3 = 2.
6
Eliminación gaussiana en R. Forma escalonada. Solución de un sistema.
En R, usando la librería pracma disponemos de la función rref que permite, dada una matriz A, obtener una
matriz A0 escalonada equivalente (en la forma de Gauss-Jordan, escalonada “hacia abajo y hacia arriba”).
Recuerda instalar la librería para usarla si no lo has hecho previamente.
Por ejemplo, dada la matriz de uno de nuestros ejemplos previos,


1
2 −3
6
 2 −1
4 2 
4
3 −2 14
la forma escalonada se obtiene así:
library(pracma)
(A = matrix(c(1, 2, -3, 6, 2, -1, 4, 2, 4, 3, -2, 14), nrow = 3, byrow = TRUE))
##
[,1] [,2] [,3] [,4]
## [1,]
1
2
-3
6
## [2,]
2
-1
4
2
## [3,]
4
3
-2
14
rref(A)
##
[,1] [,2] [,3] [,4]
## [1,]
1
0
1
2
## [2,]
0
1
-2
2
## [3,]
0
0
0
0
En particular, con el truco habitual de contar los booleanos como 1/0, esta función puede usarse para calcular
el rango de la matriz así:
rango = function(A){
library(pracma)
return(sum(rowSums(rref(A)) > 0))
}
rango(A)
## [1] 2
En cualquier caso, la función rref está pensada para usos básicamente pedagógicos. Y aunque es suficiente
para los problemas que vamos a encontrar en este curso, no es la más adecuada para un uso profesional. Más
adelante hablaremos de una forma más robusta y numéricamente eficiente de obtener el rango.
Si lo que deseamos es simplemente obtener las soluciones de un sistema, entonces podemos usar la función
solve, que usa como argumentos la matriz A de coeficientes y la matriz B de términos independientes. Por
ejemplo, para resolver el siguiente sistema (que ya hemos visto)

 2x1 + x2 + 4x3 = 1
2x1 + 2x2 + 7x3 = 0

x1 − x2 − 3x3 = 1
haríamos esto:
7
(A = matrix(c(2, 1, 4, 2, 2, 7, 1, -1, -3), nrow = 3, byrow = TRUE))
##
[,1] [,2] [,3]
## [1,]
2
1
4
## [2,]
2
2
7
## [3,]
1
-1
-3
(B = c(1, 0, 1)) # Fijate en que os introducimos como un vector.
## [1] 1 0 1
(solucion = solve(A,
B))
## [1] -8.881784e-16 -7.000000e+00
2.000000e+00
Es muy importante que te fijes en que el valor de la primera variable en la solución. La solución exacta, que
hemos obtenido antes, es 0. R nos devuelve aquí un un valor extremadamente pequeño porque utiliza métodos
numéricos (aproximados) para resolver este sistema. En casos como este puede resultar interesante la función
zapsmall que elimina estas diferencias extremadamente pequeñas:
zapsmall(solucion)
## [1]
0 -7
2
y, como ves, nos devuelve en este caso la solución exacta del sistema.
Teorema de Rouche-Frobenius. Matrices no singulares.
La notación de la matriz ampliada nos permite distinguir entre el rango de la matriz de coeficientes A y el
rango de la matriz ampliada (A|B). Los ejemplos anteriores apuntan a que la relación entre esos dos números
es la clave para la clasificación de los sistemas de ecuaciones lineales.
Teorema de Rouché-Frobenius. Dado un sistema de m ecuaciones lineales con n incógnitas:
• Si rango(A|B) > rango(A) el sistema es incompatible (no tiene soluciones).
• Si rango(A|B) = rango(A) = n entonces el sistema es compatible y determinado (tiene solución única).
• Si rango(A|B) = rango(A) < n entonces el sistema es compatible e indeterminado (tiene infinitas
soluciones).
A la vista de esta tipología de sistemas de ecuaciones lineales, parece claro que el modelo ideal de sistema con
solución única es un sistema
AX = B
de n ecuaciones con n incógnitas, en el que la matriz A es de rango máximo n. Si hay menos ecuaciones
que incógnitas el sistema no puede tener solución única, y si hay más soluciones que incógnitas, para el
que sistema admita soluciones tiene que haber ecuaciones redundantes. Estas observaciones justifican la
importancia del tipo especial de matrices que vamos a definir aquí:
Una matriz cuadrada A de n filas es no singular si su rango es igual a n. Y es singular si ese rango es menor.
8
Otra forma de calcular el rango en R. Aunque, como hemos dicho, se puede obtener el rango a partir
del resultado de rref, e incluso es posible programar una función que haga esta tarea por nosotros, a veces
queremos una forma rápida de llegar al rango de una matriz. Una forma de hacerlo es esta:
A
##
[,1] [,2] [,3]
## [1,]
2
1
4
## [2,]
2
2
7
## [3,]
1
-1
-3
(rango = qr(A)$rank)
## [1] 3
Más adelante explicaremos qué es esa función qr que hemos utilizado, pero no es necesario entenderlo para
obtener el rango.
El caso especial de los sistemas homogéneos.
Un sistema lineal es homogéneo cuando B = 0. Es decir, cuando todos los términos independientes son
0. Es fácil ver que en ese caso el sistema siempre es compatible, porque siempre admite la solución trivial
x1 = x2 = · · · = xn = 0. Así que en los sistemas homogéneos la discusión se limita a analizar si son
determinados (sólo la solución trivial) o indeterminados (soluciones no triviales) y a comprender la estructura
de las soluciones en el caso de soluciones no triviales. Pero para eso es necesario el lenguaje de la teoría de
espacios vectoriales, en el que no vamos a entrar aquí. Nos conformaremos con decir esto:
Teorema: Si la matriz del sistema homogéneo
A·X =0
es de dimensiones m × n y el rango de A es r, entonces el conjunto de soluciones del sistema es un subespacio
vectorial de Rn de dimensión d = n − r.
Inversa de una matriz cuadrada.
El sistema de ecuaciones lineales más sencillo es el que tiene una única variable y una única ecuación:
ax = b
Incluso en este caso tan sencillo ya aparece la tipología completa de sistemas de ecuaciones lineales que hemos
descrito. Por ejemplo la ecuación
0x = 1
es incompatible, mientras que
0x = 0
es compatible indeterminada. Los casos compatibles determinados, como
3x = 7
se resuelven, como aprendemos todos en la escuela, multiplicando los dos miembros de la ecuación por
3−1 =
9
1
3
así:
y aprovechando el hecho de que
7
3
x=
3
3
3
=1
3
Generalizando, este método se puede resumir así: la solución de
ax = b
es
x = a−1 b siempre que sea
b 6= 0
La idea de inversa de una matriz es una generalización de este resultado. Ante un sistema lineal como
AX = B
podemos pensar en “multiplicar por A−1 ”, donde A−1 es una matriz que multiplicada por A produce el
equivalente matricial del 1, que es la matriz identidad de dimensión n


1 0 ··· 0
 0 1 ··· 0 


I= . . .
. . ... 
 .. ..

0 0 ··· 1
Cuando sea necesario recordar la dimensión, escribiremos In en lugar de I.
Resulta que, para que estas ideas puedan aplicarse de forma sencilla, tenemos que empezar limitando nuestro
trabajo al caso de las matrices cuadradas.
Dada una matriz cuadrada A de dimensiones n × n, la matriz inversa de A, si existe, es una matriz A−1 ,
también de dimensiones n × n, que cumple:
A−1 · A = In
¿Cuándo existe la inversa? El rango nos permite responder a esa pregunta.
Teorema: La inversa de la matriz A existe si y sólo si la matriz A es de rango máximo, rango(A) = n.
EN R, para obtener la inversa podemos usar de nuevo solve, en este caso sin términos independientes.
(A = matrix(c(2, 1, 4, 2, 2, 7, 1, -1, -3), nrow = 3, byrow = TRUE))
##
[,1] [,2] [,3]
## [1,]
2
1
4
## [2,]
2
2
7
## [3,]
1
-1
-3
(invA = solve(A))
##
[,1] [,2] [,3]
## [1,]
-1
1
1
## [2,] -13
10
6
## [3,]
4
-3
-2
y podemos comprobar que el resultado es, de hecho, la inversa, multiplicando ambas matrices (recuerda que
en R el producto matricial es %*%)
10
invA %*% A
##
[,1]
[,2]
[,3]
## [1,] 1.000000e+00 4.440892e-16 1.332268e-15
## [2,] -8.881784e-16 1.000000e+00 -1.065814e-14
## [3,] 0.000000e+00 0.000000e+00 1.000000e+00
De nuevo aparecen los problemas asociados al redondeo, que podemos esquivar usando zapsmall
zapsmall(invA %*% A)
##
[,1] [,2] [,3]
## [1,]
1
0
0
## [2,]
0
1
0
## [3,]
0
0
1
Volviendo a la idea inicial, cuando existe la inversa podemos usarla para resolver el sistema de ecuaciones.
Teorema: Dado el sistema de ecuaciones lineales
AX = B
si A es una matriz cuadrada de dimensiones n × n y rango máximo n, entonces la solución del sistema viene
dada por el producto matricial:
X = A−1 B
Veámoslo en el ejemplo del sistema
2x1 + x2 + 4x3 = 1
2x1 + 2x2 + 7x3 = 0

x1 − x2 − 3x3 = 1


cuya matriz es la matriz A de la que hemos calculado la inversa con R. Para obtener la solución basta con
hacer (de nuevo, atención al redondeo):
zapsmall(invA %*% c(1, 0, 1))
##
[,1]
## [1,]
0
## [2,]
-7
## [3,]
2
Determinantes.
Seguramente habrás aprendido, en cursos previos de matemáticas, que se puede usar determinantes para
resolver sistemas de ecuaciones lineales. Dada una matriz A, representaremos su determinante con det(A) o
con |A|. Aquí no nos vamos a ocupar en detalle de las propiedades de los determinantes, porque en este curso
es prioritaria la eficiencia de los métodos. Y los métodos basados en el cálculo elemental de determinantes,
tal y como se explica en las matemáticas elementales, son extremadamente ineficientes. Nos limitaremos a
algunas observaciones que pueden resultar útiles:
• La inversa de una matriz cuadrada A existe si y sólo si el determinante de A es distinto de 0.
11
• Esa condición, det(A) 6= 0 es equivalente al hecho de que A sea no songular; es decir, de rango n.
• Y, para quienes conozcan el lenguaje del Álgebra Lineal, también es equivalente al hecho de que los
vectores que forman las columnas de A son linealmente independientes.
• En R existe una función det para calcular el determinante de una matriz:
A
##
[,1] [,2] [,3]
## [1,]
2
1
4
## [2,]
2
2
7
## [3,]
1
-1
-3
det(A)
## [1] -1
Aunque en la matemática eleemntal se explican fórmulas para calcular determinantes de matrices 2 × 2 o 3 × 3,
e incluso otras basadas en la idea de desarrollar el determinante por filas o columnas, es importante saber
que la forma más eficiente de calcular el determinante de una matriz (salvo en las dimensiones más bajas, 2
o 3) es usando algún procedimiento relacionado con la eliminación gaussiana (como los procedimientos de
factorización que se explican más adelante).
Complejidad computacional del método de eliminación. Número de operaciones.
Cuando disponemos de varios métodos para resolver un mismo problema tiene sentido compararlos, para
tratar de decidir cuál es, en algún sentido, el mejor método. Naturalmente, la elección depende del criterio
que se utilice para comparar los métodos. En los métodos computacionales en los que se emplea el ordenador
para hacer cálculos, hay varios criterios que a menudo se tienen en cuenta:
• El primero en importancia suele ser la rapidez del método y es al que más espacio le vamos a dedicar
aquí.
• A menudo también es importante el uso eficiente de la memoria. Aunque los ordenadores modernos
disponen de gran capacidad, no es sensato actuar como si esa capacidad fuese ilimitada.
• Otro factor esencial, más importante aún que la memoria, es la estabilidad numérica del método. A
veces el resultado de un método numérico puede verse alterado de forma drástica por un pequeño error
de redondeen los valores de entrada (como hemos visto al hablar de caos determinista en el método de
Newton). Es preferible utilizar métodos que no sean tan sensibles a pequeñas diferencias o errores de
redondeo, siempre que sea posible.
• Finalmente, no queremos olvidarlo, la sencillez del método es también importante. Un método muy
complicado es más difícil de utilizar y más proclive a los errores.
Vamos a centrarnos, para no complicar la discusión, en el primer factor, la velocidad del método. El factor
más determinante de esa velocidad es el número de operaciones aritméticas que realiza el ordenador. Y,
dentro de esas operaciones, las multiplicaciones y divisiones son mucho más importantes que las sumas y
restas, porque el ordenador suma a una velocidad mucho mayor que la que alcanza al multiplicar.
En resumen, vamos a suponer que queremos resolver un sistema de ecuaciones lineal
AX = B
12
en el que, para simplificar, supondremos que la matriz es cuadrada, de dimensiones n × n. Entonces resolver
el sistema por el método de eliminación gaussiana requiere de aproximadamente
3
n
O
3
multiplicaciones (contamos las divisiones como multiplicaciones). La notación O() que hemos usado aquí
se lee “del orden de”. Es la notación estándar cuando queremos describir la eficiencia de un algoritmo en
función del tamaño del problema; en este caso, en función del tamaño de la matriz.
Factorización de matrices.
Factorización LU.
A menudo sucede que tenemos que resolver varios sistemas de ecuaciones lineales que tienen la misma matriz
de coeficientes A y en los que lo único que cambia son los términos independientes:
AX = B1 , AX = B2 , AX = B3 , . . . , AX = Bk .
En casos como estos, es interesante tratar de descomponer la matriz como producto de una matriz triangular
inferior L (del inglés lower) y una matriz triangular superior U (del inglés upper).
A = LU
Y la razón por la que es interesante es porque, una vez hecho esto, para resolver un sistema de la forma:
AX = B, es decir LU X = B
podemos proceder así:
1. Hacemos un cambio de variable, llamando Y = U X (los valores de Y son desconocidos, los averiguaremos
en el siguiente paso).
2. Resolvemos el sistema LY = B, y obtenemos los valores de Y .
3. Ahora usamos esos valores de Y como términos independientes para resolver el sistema Y = U X,
obteniendo los valores de X.
¿Qué se gana con esto? Pues, inicialmente, lo que se gana es que la solución de un sistema con matriz
triangular (ya sea superior o inferior) es muy fácil de obtener, como hemos visto antes: si el sistema es
triangular superior basta con sustituir de abajo hacia arriba, como hacíamos en la eliminación Gaussiana. Y
si es triangular inferior, sustituimos de arriba hacia abajo. Así que los dos pasos 2 y 3 del método anterior
son extremadamente sencillos. Técnicamente, son de orden O(n2 ) para una matriz cuadrada de n filas. Y eso
resulta muy ventajoso si, como decíamos antes, tenemos que resolver varios sistemas con la misma matriz A.
Hacemos la descomposición A = LU una sola vez, y luego la usamos para resolver todos los sistemas que
sean necesarios. Naturalmente, el proceso de obtener la factorización LU tiene un coste, que es esencialmente
el de una eliminación gaussiana (del orden O(n3 ), como sabemos). Pero, insistimos, sólo es necesario hacer
esa factorización una vez.
El principal problema de este enfoque es que la factorización LU no siempre es posible. Para entender lo
que sucede pueder ser bueno pensar en las transformaciones elementales de la eliminación gaussiana de otra
manera.
13
Matrices de transformación.
La mejor forma de entender lo que vamos a hacer es con un ejemplo. Supongamos que queremos aplicar
eliminación gaussiana a esta matriz A, para eliminar todos los elementos por debajo de la diagonal:
(A = matrix(c(2, 3, 5, -1, 4, 3, 2, 1, 4), nrow = 3))
##
[,1] [,2] [,3]
## [1,]
2
-1
2
## [2,]
3
4
1
## [3,]
5
3
4
El primer paso podría ser la transformación elemental 2f2 − 3f1 → f2 , para conseguir un 0 como primer
elemento de la fila 2. Vamos a tomar la matriz identidad de tamaño 3, que se fabrica en R con:
(I3 = diag(3))
##
[,1] [,2] [,3]
## [1,]
1
0
0
## [2,]
0
1
0
## [3,]
0
0
1
y vamos a aplicarle a esta matriz identidad la misma transformación 2f2 − 3f1 → f2 que queremos a hacer
con A:
T1 = I3
T1[2, ] = 2 * T1[2, ] - 3 * T1[1, ]
T1
##
[,1] [,2] [,3]
## [1,]
1
0
0
## [2,]
-3
2
0
## [3,]
0
0
1
Llamaremos T1 a esta matriz, donde la T viene de la palabra transformación por razones que vamos a ver
enseguida. Porque lo más interesante viene ahora. Multipliquemos esta matriz a la izquierda por A:
T1 %*% A
##
[,1] [,2] [,3]
## [1,]
2
-1
2
## [2,]
0
11
-4
## [3,]
5
3
4
¿Ves lo que ha ocurrido? El efecto de multiplicar esta matriz por la izquierda con A es que hacemos en
las filas de A las mismas transformaciones que hemos hecho en las filas de la matriz identidad. Vamos a
probar otra vez. Ahora, para conseguir un 0 en la fila 3, columna 1, usaremos la transformación elemental
2f3 − 5f1 → f3 . Hagamos esto en la matriz identidad:
14
T2 = I3
T2[3, ] = 2 * T2[3, ] - 5 * T2[1, ]
T2
##
[,1] [,2] [,3]
## [1,]
1
0
0
## [2,]
0
1
0
## [3,]
-5
0
2
Y ahora multiplicamos por T2 para seguir avanzando en la eliminación:
T2 %*% T1 %*% A
##
[,1] [,2] [,3]
## [1,]
2
-1
2
## [2,]
0
11
-4
## [3,]
0
11
-2
Finalmente, para completar la eliminación usamos una matriz T3 (¿cuál es la transformación?):
T3 = I3
T3[3, ] = T3[3, ] - T3[2, ]
T3
##
[,1] [,2] [,3]
## [1,]
1
0
0
## [2,]
0
1
0
## [3,]
0
-1
1
(U = T3 %*% T2 %*% T1 %*% A)
##
[,1] [,2] [,3]
## [1,]
2
-1
2
## [2,]
0
11
-4
## [3,]
0
0
2
La matriz que hemos obtenido es una matriz triangular superior, asi que si, como hemos hecho en el código,
la llamamos U tenemos:
T3 T2 T1 A = U
Para cualquiera con un poco de experiencia con matrices rsulta natural,al ver esta ecuación, pensar en
despejar A así:
A = (T3 T2 T1 )−1 U
y usando la relación entre inversas y productos, esto es:
A = T1−1 T2−1 T3−1 U
Esto nos lleva a pensar en las inversas de las matrices de transformación Ti . Y es bastante fácil ver que la
inversa de cada una de ellas es una nueva matriz de transformación, la matriz de la transformación inversa.
Por ejemplo, si antes hemos hecho la transformación T 2, que era
f30 = 2f3 − 5f1
15
(aquí f30 representa la fila 3 después de la transformación), despejamos obteniendo
f3 =
f30 + 5f1
2
Está claro que para volver al estado anterior tenemos que hacer esta transformación inversa. Así que la
inversa de T 2 se obtiene así:
T2inv = I3
T2inv[3, ] = (T2inv[3, ] + 5 * T2inv[1, ]) / 2
T2inv
##
[,1] [,2] [,3]
## [1,] 1.0
0 0.0
## [2,] 0.0
1 0.0
## [3,] 2.5
0 0.5
Comprobémoslo:
T2inv %*% T2
##
[,1] [,2] [,3]
## [1,]
1
0
0
## [2,]
0
1
0
## [3,]
0
0
1
Ahora hay que fijarse en un par de detalles importantes:
• Las matrices T1 y sus inversas son matrices triangulares inferiores.
• El producto de matrices triangulares inferiores es una matriz triangular inferior.
Eso significa, volviendo a nuestro ejemplo, que la matriz el producto
L = T1−1 T2−1 T3−1
es trinagular inferior. Esta matriz L aparecía en la ecuación
A = T1−1 T2−1 T3−1 U
que ahora podemos escribir:
A = LU
con la factorización que queríamos.
¿Qué hace falta para que ese método funcione? Si lo examinas atentamente, verás que lo único importante es
que podamos hacer la eliminación gaussiana sin necesidad de intercambios de posiciones entre filas.
¿Qué sucede cuando hay que intercambiar filas? Veámoslo de nuevo con un ejemplo:
(A = matrix (c(0, 1, 1, -1, 0, 3, 4, 1, 6), nrow=3))
##
[,1] [,2] [,3]
## [1,]
0
-1
4
## [2,]
1
0
1
## [3,]
1
3
6
16
Para esta matriz, en el primer paso de la eliminación gaussiana es necesario intercambiar la primera fila con,
por ejemplo, la segunda. Afortunadamente, también es posible expresar una transformación como f1 ↔ f2
usando exactamente las mismas ideas que antes. Concretamente, hacemos esa misma transformación en la
matriz identidad:
P1 = I3
temporal = P1[1, ]
P1[1, ] = P1[2, ]
P1[2, ] = temporal
P1
##
[,1] [,2] [,3]
## [1,]
0
1
0
## [2,]
1
0
0
## [3,]
0
0
1
Y ahora multiplicamos P 1 (por la izquierda) por A:
P1 %*% A
##
[,1] [,2] [,3]
## [1,]
1
0
1
## [2,]
0
-1
4
## [3,]
1
3
6
Como ves, el resultado es que hemos obtenido la transformación deseada en A. ¡Pero atención! La matriz P 1
no es triangular (ni superior ni inferior). Así que (tras un poco de reflexión) es fácil ver que si la eliminación
gaussiana requiere intercambios de posiciones de las filas, entonces la factorización LU no es posible. Pero
todavía podemos conservar buena parte de sus ventajas, introduciendo las transformaciones de reordenación
de filas en el proceso. De hecho, es posible analizar los intercambios de filas que van a ser necesarios y llevarlos
a cabo todos antes de empezar con la eliminación. Veámoslo en un nuevo ejemplo.
(A = matrix (c(1, 2, 1, 3, 6, 4, 5, 13, 7), nrow=3))
##
[,1] [,2] [,3]
## [1,]
1
3
5
## [2,]
2
6
13
## [3,]
1
4
7
Empezamos la eliminación con las dos transformaciones:
f2 − 2f1 → f2 ,
f3 − f1 → f3
Traduciendo matricialmente estas transformaciones:
T1 = I3
T1[2, ] = T1[2, ] - 2 * T1[1, ]
T1
##
[,1] [,2] [,3]
## [1,]
1
0
0
## [2,]
-2
1
0
## [3,]
0
0
1
17
T2 = I3
T2[3, ] = T2[3, ] - T2[1, ]
T2
##
[,1] [,2] [,3]
## [1,]
1
0
0
## [2,]
0
1
0
## [3,]
-1
0
1
T2 %*% T1 %*% A
##
[,1] [,2] [,3]
## [1,]
1
3
5
## [2,]
0
0
3
## [3,]
0
1
2
Hemos conseguido nuestro objetivo de hacer dos ceros en la primera columna, pero como efecto colateral el
elemento (2, 2) de la matriz resultante es cero. Podemos salir del atolladero cambiando la segunda fila con la
tercera, en la transformación f2 ↔ f3 . Matricialmente:
P1 = I3
temporal = P1[2, ]
P1[2, ] = P1[3, ]
P1[3, ] = temporal
P1
##
[,1] [,2] [,3]
## [1,]
1
0
0
## [2,]
0
0
1
## [3,]
0
1
0
y ahora, multiplicando:
P1 %*% T2 %*% T1 %*% A
##
[,1] [,2] [,3]
## [1,]
1
3
5
## [2,]
0
1
2
## [3,]
0
0
3
obtenemos una matriz escalonada (triangular inferior). Jugando con las ideas de este ejemplo es posible
demostrar que si una matriz cuadrada A es de rango máximo entonces existe una factorización de la forma:
A = P LU
en la que:
• P es el producto de varias de estas matrices de intercambio de posiciones de filas. Además, la inversa
de P es muy fácil de obtener, porque coincide con su traspuesta: P t = P −1 .
• Como antes L es triangular inferior y U es triangular superior.
18
Esta descomposición permite resolver sistemas de la forma AX = B con una estrategia similar a la que
hemos descrito en el caso de la descomposición LU . En R podemos usar la librería Matrix para obtener estas
descomposiciones. Por ejemplo
library(Matrix)
## Warning: package 'Matrix' was built under R version 3.1.3
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:pracma':
##
##
expm, lu, tril, triu
##
## The following objects are masked from 'package:base':
##
##
crossprod, tcrossprod
(A = matrix(c(2, 3, 5, -1, 4, 3, 2, 1, 4), nrow = 3))
##
[,1] [,2] [,3]
## [1,]
2
-1
2
## [2,]
3
4
1
## [3,]
5
3
4
LuA
L =
U =
P =
= expand(lu(A))
LuA$L
LuA$U
LuA$P
Como ves, la función lu (y la función auxiliar expand que interviene en el proceso) usa una notación especial,
en la que los puntos representan ceros. Para comprobar el resultado hacemos:
P %*% L %*% U
##
##
##
##
##
3 x 3 Matrix of class "dgeMatrix"
[,1] [,2] [,3]
[1,]
2
-1
2
[2,]
3
4
1
[3,]
5
3
4
Factorización de Cholesky. No vamos a profundizar más en el tema de la factorización de matrices, pero
es conveniente saber que, cuando se trabaja con matrices que tienen características especiales, existen otros
tipos de factorizaciones similares que pueden ser muy útiles. Por ejemplo, las matrices definidas positivas son
especialmente importantes en muchos problemas que involucran el cálculo de distancias. Sin ir más lejos, en
el método de mínimos cuadrados para la regresión lineal que hemos visto en la parte de Estadística del curso.
Una matriz cuadrada A de n filas es definida positiva si
X t AX > 0
19
para cualquier vector X 6= 0 (pensado aquí como una matriz columna n × 1). La matriz A es semidefinida
positiva si
X t AX ≥ 0
para cualquier vector X. lo interesante aquí es que permitimos que algún vector no nulo X pueda dar como
resultado $ Xˆt A X = 0$.
Para el caso de matrices definidas positivas se disponde de la factorización de Cholesky que se describe
en este teorema:
Teorema: Una matriz A es definida positiva si y sólo si se puede factorizar
A = U tU
donde U es una matriz triangular superior con elementos no nulos en su diagonal.
Detalles sobre las matrices definidas y semidefinidas positivas.
es definida positiva:
Por ejemplo, la siguiente matriz A
(A = matrix(c(4, -2, 14, -2, 10, -10, 14, -10, 66), nrow = 3))
##
[,1] [,2] [,3]
## [1,]
4
-2
14
## [2,]
-2
10 -10
## [3,]
14 -10
66
Una forma de comprobar si la matriz A es definida positiva, que no es muy eficiente, pero que sirve para
matrices pequeñas, es calcular los determinantes de los menores principales de A, que son las matrices que se
apoyan en la diagonal principal (en el primero hemos usado as.matrix porque al ser de dimensión 1 R lo
convierte en un escalar, y det produciría un mensaje de error):
A1 = as.matrix(A[1, 1])
det(A1)
## [1] 4
A2 = A[1:2, 1:2]
det(A2)
## [1] 36
A3 = A[1:3, 1:3]
det(A3)
## [1] 576
Cuando todos esos determinantes son positivos, la matriz es definida positiva.
En ese caso, R nos permite obtener la factorización de Cholesky usando la función chol así:
(U = chol(A))
20
##
[,1] [,2] [,3]
## [1,]
2
-1
7
## [2,]
0
3
-1
## [3,]
0
0
4
Y la comprobación de la factorización es:
t(U) %*% U
##
[,1] [,2] [,3]
## [1,]
4
-2
14
## [2,]
-2
10 -10
## [3,]
14 -10
66
que, como puedes comprobar, es la matriz A original. Como hemos dicho, se puede profundizar mucho más
en este tema de la factorización de matrices de tipo especial y sus aplicaciones. Remitimos al lector a la
bibliografía.
Más adelante vamos a usar esta otra propiedad, que permite obtener fácilmente matrices semidefinidas
positivas. Dada una matriz A cualquiera (¡cuadrada o rectangular!) las matrices AAt y At A son simétricas y
semidefinidas positivas.
Diagonalización. Autovalores y autovectores.
La factorización LU es útil, entre otras cosas, porque nos permite resolver sistemas lineales AX = B usando
la estructura de la matriz A. En concreto, reducimos sistemas de ecuaciones a una forma triangular, más
sencilla. Pero existe una forma todavía más sencilla de sistemas de ecuaciones, los que tienen matrices de
coeficientes diagonales.
Una matriz A es diagonal si todos los elementos aij en los que i 6= j son iguales a 0. Por ejemplo, la siguiente
matriz A es diagonal:


2
0 0 0
 0 −3 0 0 
 = 0.
A=
 0
0 1 0 
0
0 0 7
En R es fácil crear este tipo de matrices con la función diag (que ya hemos usado para crear matrices
identidad):
(A = diag(c(2, -3, 1, 7)))
##
##
##
##
##
[1,]
[2,]
[3,]
[4,]
[,1] [,2] [,3] [,4]
2
0
0
0
0
-3
0
0
0
0
1
0
0
0
0
7
Un sistema de ecuaciones con matriz diagonal es extremadamente sencillo, claro. Por ejemplo, con esta misma
matriz A y matriz de términos independientes


5
 −4 

B=
 0 ,
1
21
el sistema de ecuaciones lineales AX = B se reduce a:

2x1 = 5



−3x1 = −4
1x1 = 0



7x1 = 1
cuyo análisis (y solución) es inmediato.
Teniendo esto en cuenta, es natural preguntarse si será posible hacer alguna factorización de una matriz A
que la convierta o relacione de alguna manera con una matriz diagonal. Ese es el sentido de esta definición:
Una matriz cuadrada A es diagonalizable si existen una matriz P no singular y una matriz diagonal D
tales que
A = P DP −1
En un sentido amplio, la diagonalización de matrices se refiere al conjunto de técnicas disponibles para
conseguir este objetivo. La diagonalización es, en muchos problemas computacionales, un ingrediente esencial
para entender la estructura de las matrices que intervienen, la geometría subyacente al problema. Y un
ingrediente clave de la diagonalización son las nociones de autovalor y autovector, de las que nos vamos a
ocupar a continuación.
Autovalores y autovectores.
Sea A una matriz cuadrada de n filas. Un autovalor λ de una matriz es un número (real o complejo) que
tiene la propiedad de que:
rango(A − λI) < n
o, de forma equivalente,
det(A − λI) = 0,
donde I es la matriz identidad de dimensión n. Esta
a11 − λ
a12
a21
a
−λ
22
..
..
.
.
an1
an2
última ecuación se puede visualizar así:
···
a1n
···
a2n
= 0.
..
..
.
.
· · · ann − λ Y al desarrollar el determinante se obtiene una ecuación polinómica de grado n en λ. El polinomio det(A − λI)
se denomina polinomio característico de la matriz A. Los autovalores se llaman también valores propios
o eigenvalores (en inglés, eigenvalues).
Por ejemplo, dada la matriz A
A=
−12 −21
10
17
para encontrar sus autovalores planteamos la ecuación:
−12 − λ
|A − λI| = 10
Desarrollando el determinante esto es:
−21 =0
17 − λ λ2 − 5λ + 6 = 0
Y esta ecuación tiene dos soluciones, λ1 = 2 y λ2 = 3, que son los dos autovalores de esta matriz A.
22
Vamos ahora a definir los autovectores (o vectores propios o eigenvectores, en inglés eigenvectors) asociados
con el autovalor λ. Puesto que el sistema
(A − λI)X = 0
es homogéneo siempre tiene, al menos, la solución nula X = 0, sea cual sea el valor de λ. Pero cuando λ es
un autovalor de A, entonces el rango de la matriz (A − λI) es menor que n y eso nos garantiza la existencia
de soluciones no triviales del sistema. Eso justifica esta definición:
Un autovector v de la matriz A asociado con el autovalor λ es un vector no nulo que cumple
(A − λI)v = 0
o, de forma equivalente,
Av = λv.
Por ejemplo, para la matriz A que hemos visto antes, si se usa el autovector λ1 = 2, un autovector v = (v1 , v2 )t
(vector columna) tiene que ser una solución no nula de este sistema de ecuaciones:
−14 −21
v1
0
(A − 2I)v =
=
10
15
v2
0
Esto es:
−14v1 − 21v2 = 0
10v1 + 15v2 = 0
Dividiendo la primera ecuación por −7 y la segunda por 5 se obtiene este sistema equivalente:
2v1 + 3v2 = 0
2v1 + 3v2 = 0
que muestra de forma evidente que el rango del sistema es 1. ¡Naturalmente! Para eso hemos sustituido el
autovalor λ1 = 2. Precisamente lo que buscábamos era esta caída del rango, para hacer posible la existencia
de soluciones no nulas de este sistema. Un autovector es, entonces, una de esas soluciones no nulas. Por
ejemplo, podemos tomar v = (−3, 2) como autovector asociado a λ1 .
Para el segundo autovalor λ2 = 3 tenemos que resolver el sistema:
−15v1 − 21v2 = 0
10v1 + 14v2 = 0
Y, de nuevo, dividiendo la primera ecuación por −3 y la segunda por 2 se llega a
5v1 + 7v2 = 0
5v1 + 7v2 = 0
Con lo que podemos tomar como autovector v 0 = c(−7, 5), asociado con el autovalor λ2 = 3.
Autovectores y diagonalización.
Los autovalores y autovectores están relacionados con el problema de la diagonalización. Veámoslo en el
ejemplo anterior. Si creamos una matriz diagonal D con los autovalores de A que hemos obtenido:
D = diag(c(2, 3))
y una matriz P que tenga como columnas los autovectores de A que hemos calculado (en el mismo orden que
los correspondientes autovalores):
23
(P = matrix(c(-3, 2, -7, 5), nrow = 2))
##
[,1] [,2]
## [1,]
-3
-7
## [2,]
2
5
Entonces al hacer el producto P −1 DP se obtiene:
P %*% D %*% solve(P)
##
[,1] [,2]
## [1,] -12 -21
## [2,]
10
17
que precisamente la matriz A. En resumen, hemos comprobado que se cumple esta factorización de la matriz
A:
A = P DP −1
Y esa es, precisamente, nuestra definición de diagonalización.
la diagonalización requiere muchas veces que se tengan en cuenta los posibles autovalores complejos de la
matriz. Por ejemplo, dada la matriz A:
3 −4
A=
4
3
sus autovalores y los autovectores asociados se calculan con R así:
(A = matrix(c(3, 4, -4, 3), nrow = 2))
##
[,1] [,2]
## [1,]
3
-4
## [2,]
4
3
(autoA = eigen(A))
##
##
##
##
##
##
##
$values
[1] 3+4i 3-4i
$vectors
[,1]
[,2]
[1,] 0.7071068+0.0000000i 0.7071068+0.0000000i
[2,] 0.0000000-0.7071068i 0.0000000+0.7071068i
Fíjate en particular en que los autovalores son los números complejos (conjugados) 3 + 4i y 3 − 4i. Si
procedemos como en el ejemplo anterior, podemos definir una matriz diagonal con estos autovalores:
(D = diag(autoA$values))
##
[,1] [,2]
## [1,] 3+4i 0+0i
## [2,] 0+0i 3-4i
y una matriz P cuyas columnas son los autovectores:
24
(P = autoA$vectors)
##
[,1]
[,2]
## [1,] 0.7071068+0.0000000i 0.7071068+0.0000000i
## [2,] 0.0000000-0.7071068i 0.0000000+0.7071068i
Entonces, de nuevo:
P %*% D %*% solve(P)
##
[,1] [,2]
## [1,] 3+0i -4+0i
## [2,] 4+0i 3+0i
es la matriz A. A primera vista puede confundirte el aspecto complejo, pero fíjate en que las partes imaginarias
de todos los elementos son 0. En R usamos las funciones Re e Im para obtener, respectivamente, la parte real
e imaginaria de un número complejo. Así que aquí podemos hacer más evidente la diagonalización así:
Re(P %*% D %*% solve(P))
##
[,1] [,2]
## [1,]
3
-4
## [2,]
4
3
Im(P %*% D %*% solve(P))
##
[,1] [,2]
## [1,]
0
0
## [2,]
0
0
Condiciones para la posible diagonalización de una matriz.
¿Qué matrices son diagonalizables? Empecemos por decir que no todas, desde luego, ni siquiera teniendoen
cuenta los autovalores y autovectores complejos. Los ejemplos anteriores indican que la matriz formada por
los autovectores juega un papel importante en la diagonalización. En particular hemos necesitado construir
una matriz no singular de autovectores. Y, con un poco de observación, es fácil ver que esa es una condición
necesaria y suficiente para que la matriz A sea diagonalizable. Podemos decir entonces que una matriz no
es diagonalizable cuando no tiene suficientes autovectores suficientemente distintos. Técnicamente, esto se
puede enunciar así:
Teorema: Una matriz cuadrada A de dimensión n es diagonalizable si y sólo existe una matriz no singular
P cuyas columnas son todas autovectores de A.
Para ver un ejemplo de una matriz no diagonalizable, tomemos esta matriz A.
(A = matrix(c(3, 1, -1, 1), nrow = 2))
##
[,1] [,2]
## [1,]
3
-1
## [2,]
1
1
Si le preguntamos a R por sus autovalores y autovectores obtenemos:
25
eigen(A)
##
##
##
##
##
##
##
$values
[1] 2 2
$vectors
[,1]
[,2]
[1,] 0.7071068 0.7071068
[2,] 0.7071068 0.7071068
Sólo hay un autovalor, λ = 2. Eso no es necesariamente malo. Pero lo malo es que cuando miramos los
autovectores, vemos que R repite el mismo en las dos columnas de la matriz. Y la razón por la que R hace eso
es porque le ha sido imposible encontrar un segundo autovector que no fuera, esencialmente, una repetición
del primero. Con más detalle, en este caso la matriz A − λI = A − 2I es:
A - 2 * diag(2)
##
[,1] [,2]
## [1,]
1
-1
## [2,]
1
-1
Y por lo tanto el sistema de ecuaciones para encontrar autovectores v es:
v1 − v2 = 0
v1 − v2 = 0
Podemos tomar como autovector v = (1, 1). Pero si tratamos de encontrar un segundo autovector, como
v 0 = (2, 2), descubriremos que la matriz cuyas columnas son v y v 0 tiene rango 1. Es imposible encontrar el
segundo autovector linealmente independiente de (1, 1) que necesitaríamos para la diagonalización.
Este último ejemplo hace evidente la pregunta de si las matrices diagonalizables son muy raras. Vamos a ver
una familia especial, pero muy importante, de matrices diagonalizables.
Una matriz A es simétrica si At = A. Es decir, si A permanece invariante al intercambiar filas por columnas.
Las matrices simétricas aparecen de manera natural en muchos problemas importantes, así que el siguiente
resultado es especialmente útil.
Teorema: 1. Una matriz real simétrica siempre es diagonalizable. Además los autovalores de A son todos
números reales y se puede elegir la matriz P de forma que se cumpla P −1 = P t . 2. Para una matriz simétrica
A, ser definida positiva es equivalente a que todos los autovalores de la matriz sean números reales positivos.
Y ser semidefinida negativa es equivalente a que todos los autovalores sean no negativos (alguno puede ser 0).
La propiedad P −1 = P t significa que la matriz P es una matriz ortogonal, porque su traspuesta coincide
con su inversa.
Sistema de ecuaciones lineales y diagonalización. Cuando la matriz A es diagonalizable, para resolver
un sistema de ecuaciones lineales
AX = B
escribimos
y entonces el sistema se convierte en:
A = P DP −1 con P −1 = P t
P DP −1 X = B
26
Si usamos el cambio de variables Y = P −1 X entonces esto es:
P DY = B
o, lo que es lo mismo,
DY = P t B
Resolver este sistema es inmediato, porque D es diagonal. Y una vez que sabemos cuál es el vector Y , obtener
X es tan sencillo como multiplicar por P :
X = PY
Forma canónica de Jordan.
¿Y qué sucede cuando la matriz A no es diagonalizable, porque no hay suficientes autovectores independientes?
Una primera respuesta es que en ese caso, sea cual sea la matriz cuadrada A, todavía podemos obtener una
factorización, que es una suerte de generalización de la diagonalización. Concretamente, podemos escribir
A = P JP −1
donde P es una matriz no singular, y J es la denominada forma canónica de Jordan asociada a la matriz A.
Esta matriz J es una matriz constituida por bloque, situados a lo largo de la diagonal, donde cada bloque es
o bien una matriz de forma diagonal:


λ 0 ··· 0
 0 λ ··· 0 


 ..
..
.. 
..
 .
.
. 
.
0 ···
0 λ
o bien una matriz con unos encima de la diagonal como esta:

λ 1 0 ···
 0 λ 1 ···

 .. .. . .
..
 . .
.
.
0 0 ···
0
0
0
..
.





λ
Tanto la diagonalización como la forma canónica de Jordan juegan un papel importante en el análisis teórico
del comportamiento de los sistemas de ecuaciones diferenciales lineales, de los que hablaremos más adelante
en el curso. Pero desde el punto de vista numérico, es preferible utilizar la factorización que vamos a describir
en el último apartado.
Descomposición en valores singulares (SVD).
Cuando una matriz no es diagonalizables, como hemos dicho, podemos construir la forma canónica de Jordan
asociada. Pero es un hecho conocido que el cálculo de la forma de Jordan no tiene buenas propiedades,
desde el punto de vista de la estabilidad numérica. Por esa razón es más adecuado utilizar factorizaciones
alternativas, como la que vamos a ver aquí:
Teorema: Dada una matriz A cualquiera, existen dos matrices ortogonales U y V y una matriz diagonal Σ
con valores no negativos en la diagonal tales que se tiene esta factorización:
A = U ΣV t
Aviso: no te dejes confundir por la notación. U no es triangular superior. A pesar de que pueda resultar
confusa la notación U, V es la más usada aquí y, como veremos, también es la que usa R.
La factorización que hemos descrito es la descomposición en valores singulares (en inglés, singular value
decomposition, a menudo abreviado como SVD)
27
El punto de partida es la observación de que las matrices simétricas son todas diagonalizables, y su
diagonalización tiene buenas propiedades. Si empezamos con una matriz A que no es simétrica, podemos
construir a partir de ella la matriz AAt que sí es simétrica. Veámoslo en un ejemplo. Empezamos con:
(A = matrix(c(-2, 0, 2, 0, 2, 0, -2, 0, 0), nrow = 3))
##
[,1] [,2] [,3]
## [1,]
-2
0
-2
## [2,]
0
2
0
## [3,]
2
0
0
y primero calculamos la matriz S = AAt :
(S = A %*% t(A))
##
[,1] [,2] [,3]
## [1,]
8
0
-4
## [2,]
0
4
0
## [3,]
-4
0
4
Esta matriz S sí que es simétrica. Como hemos visto,eso significa que S es diagonalizable y semidefinida
positiva. En particular, sus autovalores son números reales no negativos (puede que alguno sea 0). En este
ejemplo:
eigen(S)
##
##
##
##
##
##
##
##
$values
[1] 10.472136
4.000000
1.527864
$vectors
[,1] [,2]
[,3]
[1,] 0.8506508
0 0.5257311
[2,] 0.0000000
1 0.0000000
[3,] -0.5257311
0 0.8506508
Si diagonalizamos esa matriz S tendremos S = P DP −1 , donde:
(D = diag(eigen(S)$values))
##
[,1] [,2]
[,3]
## [1,] 10.47214
0 0.000000
## [2,] 0.00000
4 0.000000
## [3,] 0.00000
0 1.527864
(P = eigen(S)$vectors)
##
[,1] [,2]
[,3]
## [1,] 0.8506508
0 0.5257311
## [2,] 0.0000000
1 0.0000000
## [3,] -0.5257311
0 0.8506508
28
Pero puesto que los autovalores son números no negativos, podemos tomar la raíz cuadrada de la matriz D, a
la que llamaremos Σ (es la notación habitual, qué se le va a hacer), de manera que Σ es una matriz diagonal
que cumple:
Σ2 = D
√
Es decir, que si λ > 0 ocupa la posición i en la diagonal de D, entonces λ ocupa la posición i en la diagonal
de Σ (los elementos que no aparecen en las siguientes matrices son todos cero):


 √

λ1 √
λ1




λ2
λ2








..
..




.
.




√



.
λr
D=
Σ=
λr

,




0
0








.
.
.
.




.
.
0
0
El número r coincide con el rango de A. Vamos a usar también esta matriz:
 1

√
 λ1



1


√


λ2




..


.


K=
.
1


√


λr




0




.
..


0
Aunque no podemos entrar en detalle a justificar el resultado, la idea es que si tomamos la matriz
V = P,
y usamos como U el producto:
U = AV K
entonces se cumplirá A = U ΣV t . Veámoslo en un ejemplo. Empecemos con la matriz:
(A = matrix(c(2, 1, 1, 1, 2, 1, -2, -2, -1), nrow = 3))
##
[,1] [,2] [,3]
## [1,]
2
1
-2
## [2,]
1
2
-2
## [3,]
1
1
-1
Como primer paso construimos la matriz At A, que vamos a guardar en la variable H de R:
(H = (t(A) %*% A))
##
[,1] [,2] [,3]
## [1,]
6
5
-7
## [2,]
5
6
-7
## [3,]
-7
-7
9
29
Para diagonalizarla usamos la función eig de R, que nos proporcionan los autovalores de H y una matriz
ortogonal de autovectores.
eigH = eigen(H)
(D = diag(eigH$values))
##
[,1] [,2]
[,3]
## [1,] 19.94987
0 0.00000000
## [2,] 0.00000
1 0.00000000
## [3,] 0.00000
0 0.05012563
(P = eigH$vectors)
##
[,1]
[,2]
[,3]
## [1,] -0.5245245 7.071068e-01 0.4742089
## [2,] -0.5245245 -7.071068e-01 0.4742089
## [3,] 0.6706326 -5.551115e-16 0.7417897
Para comprobar que tenemos la diagonalización de H hagamos el producto:
P %*% D %*% solve(P)
##
[,1] [,2] [,3]
## [1,]
6
5
-7
## [2,]
5
6
-7
## [3,]
-7
-7
9
Una vez comprobado esto, vamos a construir Σ y K:
(Sigma = sqrt(D))
##
[,1] [,2]
[,3]
## [1,] 4.466528
0 0.0000000
## [2,] 0.000000
1 0.0000000
## [3,] 0.000000
0 0.2238875
(K = diag(1/diag(Sigma)))
##
[,1] [,2]
[,3]
## [1,] 0.2238875
0 0.000000
## [2,] 0.0000000
1 0.000000
## [3,] 0.0000000
0 4.466528
Y ahora usamos P como matriz V , a partir de la cual obtenemos también U :
V = P
(U = A %*% P %*% K)
30
##
[,1]
[,2]
[,3]
## [1,] -0.6525961 0.7071068 -0.2722469
## [2,] -0.6525961 -0.7071068 -0.2722469
## [3,] -0.3850153 0.0000000 0.9229102
Nos falta simplemente comprobar la identidad A = U ΣV t :
U %*% Sigma %*% t(V)
##
[,1] [,2] [,3]
## [1,]
2
1
-2
## [2,]
1
2
-2
## [3,]
1
1
-1
Afortunadamente, R nos proporciona una forma mucho más sencilla de hacer este trabajo mediante la función
svd:
(svdA = svd(A))
##
##
##
##
##
##
##
##
##
##
##
##
##
##
$d
[1] 4.4665282 1.0000000 0.2238875
$u
[,1]
[,2]
[,3]
[1,] -0.6525961 7.071068e-01 -0.2722469
[2,] -0.6525961 -7.071068e-01 -0.2722469
[3,] -0.3850153 8.326673e-17 0.9229102
$v
[,1]
[,2]
[,3]
[1,] -0.5245245 7.071068e-01 0.4742089
[2,] -0.5245245 -7.071068e-01 0.4742089
[3,] 0.6706326 2.220446e-16 0.7417897
Comentarios sobre la SVD.
• Como hemos advertido, la SVD es una técnica factorización numéricamente más eficiente y se usa a
menudo para trabajar con grandes matrices o en problemas en los que los condicionantes numéricos son
importantes.
• Si la matriz A es simétrica definida positiva, entonces la SVD coincide con la diagonalización de A.
• Aunque en las secciones anteriores nos hemos centrado en el caso de matrices cuadradas, una de las
ventajas de la SVD es que se aplica sin ninguna modificación a matrices rectangulares.
Veamos un ejemplo de este última afirmación.
(A = matrix(c(1:12), nrow = 3))
##
[,1] [,2] [,3] [,4]
## [1,]
1
4
7
10
## [2,]
2
5
8
11
## [3,]
3
6
9
12
31
svdA = svd(A)
svdA$u %*% diag(svdA$d) %*% t(svdA$v)
##
[,1] [,2] [,3] [,4]
## [1,]
1
4
7
10
## [2,]
2
5
8
11
## [3,]
3
6
9
12
Como ves, el algoritmo se aplica exactamente igual para dar la factorización que buscamos. La diferencia en
este caso es que U es una matriz (3, 3), mientras que V es (3, 4) y la matriz diagonal Σ es (3, 3).
• Las aplicaciones de la SVD son numerosísimas. Recomendamos al lector la consulta de la página de la
Wikipedia (en inglés) sobre la SVD, en la que encontrará información adicional:
http://en.wikipedia.org/wiki/Singular_value_decomposition
Factorización QR.
La factorización QR es, junto con la descomposición SVD que hemos visto, uno de los algoritmos matriciales
numéricos más importantes. En esencia se trata de que cualquier matriz cuadrada A (de números reales) se
puede factorizar en la forma:
A = QR
donde Q es una matriz orotogonal (Q−1 = Qt ) mientras que R es una matriz triangular superior. En R
disponemos de la función qr (que ya hemos usado antes para calcular un rango) que permite obtener esa
factorización. Por ejemplo, con la matriz (singular)
(A = matrix(c(1, 2, 4, 2, -1, 3, -3, 4, -2), nrow = 3))
##
[,1] [,2] [,3]
## [1,]
1
2
-3
## [2,]
2
-1
4
## [3,]
4
3
-2
se obtiene
(qrA = qr(A))
##
##
##
##
##
##
##
##
##
##
##
##
##
##
$qr
[,1]
[,2]
[,3]
[1,] -4.5825757 -2.6186147 6.546537e-01
[2,] 0.4364358 2.6726124 -5.345225e+00
[3,] 0.8728716 0.1157322 -1.332268e-15
$rank
[1] 2
$qraux
[1] 1.218218e+00 1.993280e+00 1.332268e-15
$pivot
[1] 1 2 3
32
##
## attr(,"class")
## [1] "qr"
Este resultado es, tal como aparece, difícil de interpretar. Para ahorrar espacio, R está incluyendo en una
única matriz la información necesaria para recuperar las dos matrices Q y R de la descomposición. No hay
que preocuparse. Tenemos dos funciones auxiliares que facilitan el trabajo:
(Q = qr.Q(qrA))
##
[,1]
[,2]
[,3]
## [1,] -0.2182179 0.5345225 -0.8164966
## [2,] -0.4364358 -0.8017837 -0.4082483
## [3,] -0.8728716 0.2672612 0.4082483
(R = qr.R(qrA))
##
[,1]
[,2]
[,3]
## [1,] -4.582576 -2.618615 6.546537e-01
## [2,] 0.000000 2.672612 -5.345225e+00
## [3,] 0.000000 0.000000 -1.332268e-15
y para comprobar que esta es la descomposición que esperábamos, multipliquemos:
Q %*% R
##
[,1] [,2] [,3]
## [1,]
1
2
-3
## [2,]
2
-1
4
## [3,]
4
3
-2
Fíjate además en que uno de los componentes de la respuesta que se obtiene con qr es rank, el rango de la
matriz. El algoritmo QR es un algoritmo muy adecuado, desde el punto de vista numérico, para el cálculo
del rango de una matriz. Y es, sin duda, el algoritmo más recomendado para el cálculo de los autovalores de
una matriz. La descomposición QR juega un papel importante también en algunos métodos para la solución
del problema de ajuste por mínimos cuadrados de la Estadística. Para obtener más información consulta la
Bibliografía y la página de la Wikipedia (en inglés):
http://en.wikipedia.org/wiki/QR_decomposition
Sistemas de ecuaciones no lineales.
Comparado con lo que hemos visto en el caso de los sistemas lineales, el mundo de los sistemas no lineales es
extremadamente complicado, y gran parte de los problemas que surgen son muy difíciles.
Un primer paso para adentrarnos en ese territorio puede llevarnos a considerar los sistemas de ecuaciones
polinómicas.
33
Sistemas de ecuaciones polinómicas.
Un sistema de ecuaciones polinómico es de la forma:

f1 (x1 , x2 , . . . , xn ) = 0



 f2 (x1 , x2 , . . . , xn ) = 0
..

.



fk (x1 , x2 , . . . , xn ) = 0
donde f1 , f2 , . . . , fk son polinomios. Por ejemplo, el siguiente sistema
3x2 + 2xy + y 2 − x − 2 = 0
2x2 + xy − y 2 − 1 = 0
aparece al estudiar los puntos de corte de la elipse y la hipérbola de la siguiente figura:
Existe un procedimiento general para abordar este tipo de problemas, que sepuede ver como una generalización
de la Eliminación Gaussiana y que consiste en el cálculo de lo que se conoce como una base de Groebner
para el sistema (las bases de Groebner son un hallazgo reciente, hecho en 1965 por Buchberger). Este tipo
de métodos pertenecen al ámbito del Álgebra Computacional, y son métodos simbólicos. Por eso no es de
extrañar que R, que es un programa esencialmente numérico, no incluya herramientas para el cálculo de
bases de Groebner. Pero podemos usar GeoGebra o Wolfram Alpha para esto. Por ejemplo, en la vista CAS
de GeoGebra puedes teclear uno tras otro estos comandos:
f(x, y) := 2x2 + x y - y2 - 1
g(x, y) :=3x2 + 2x y + y2 - x - 2
GroebnerLex[{f, g}, {x, y}]
y obtendrás un nuevo sistema de dos ecuaciones.
−22x4 + 13x3 + 29x2 − 6x − 9 = 0
−9y − 22x3 + 13x2 + 14x − 3 = 0
con dos propiedades importantes:
• El nuevo sistema es equivalente al primero, en el mismo sentido que vimos para la eliminación gaussiana:
sus conjuntos de soluciones coinciden.
• La primera ecuación del sistema sólo contiene la variable x. Hemos eliminado la y. Fíjate además que
la segunda ecuación permite despejar la y en función de la x de maera muy sencilla.
Por esta última razón decimos que es un sistema de eliminación. Una vez que disponemos de este sistema
podemos utilizar algún procedimiento de los que hemos estudiado (como el método de Newton) para encontrar
las soluciones x de la primera ecuación. Sustituimos en la segunda y obtenemos los correspondientes valores
de la y.
Si quieres pedirle a GeoGebra que haga estos pasos de forma más o menos automática por nosotros, prueba a
ejecutar:
Resuelve[{f(x, y)=0, g(x, y)=0}, {x, y}]
y obtendrás como respuesta:
34
{{x = 0.6920811632169, y = 0.6247819618866}, {x = 1.260739706209, y = (-0.9747143706442)}}
El mayor inconveniente de estos métodos de eliminación simbólicos es que los algoritmos que se conocen son
muy poco eficientes comparados con los métodos numéricos. En cualquier caso, muchas veces se utiliza una
cierta combinación híbrida de métodos simbólicos (que son exactos) para localizar las raíces, seguidos de
métodos numéricos para estimar esas raíces de forma eficiente.
Para ver ejemplos con Wolfram Alpha prueba a ejecutar:
GroebnerBasis({2 * x^2 + x * y - y^2 - 1, 3 * x^2 + 2 * x * y + y^2 - x - 2}, {y, x} )
y también
GroebnerBasis({2 * x^2 + x * y - y^2 - 1, 3 * x^2 + 2 * x * y + y^2 - x - 2}, {x, y} )
para ver la diferencia. Por último, ejecuta
GroebnerBasis({2 * x^2 + x * y - y^2 - 1, 3 * x^2 + 2 * x * y + y^2 - x - 2}, {x}, {y} )
para una tercera versión de este método (esta última es interesante desde el punto de vista numérico). PAra
profundizar en este tema recomendamos el libro de Cox, Little y O’Shea que aparece en la Bibliografía.
Métodos numéricos para sistemas de ecuaciones no lineales.
Existen otros métodos para sistemas no polinómicos (o incluso para los polinómicos cuando la complejidad de
los métodos simbólicos lo haga aconsejable), entre los que se incluyen generalizaciones multidimensionales del
método de Newton o de los métodos iterativos de los que hemos hablado brevemente en el tema anterior. EN
este curso no disponemos de tiempo para extendernos en el estudio detallado de esos métodos (nos remitimos
a la bibliografía). Pero sí queremos mostrar un par de ejemplos de cómo se puede abordar la estimaión de
las soluciones de estos sistemas usando R. En particular, vamos a usar la función multiroot de la librería
rootSolve. Para empezar vamos a ver como encontrar aproximadamente los puntos de corte de la elipse y la
hipérbola de nuestro ejemplo previo. Empezamos por cargar la librería y definir el sistema de ecuaciones:
library(rootSolve)
##
## Attaching package: 'rootSolve'
##
## The following objects are masked from 'package:pracma':
##
##
gradient, hessian
F = function(x){c(2 * x[1]^2 + x[1] * x[2] - x[2]^2 - 1, 3 * x[1]^2 + 2 * x[1] * x[2] + x[2]^2 - x[1] Como ves, hemos definido el sistema como una función F de las variables (x1 , x2 ), y el valor de esa función
es un vector con dos componentes, que se corresponden con las ecuaciones del sistema (F : R2 → R2 ).
Una exploración gráfica preliminar (usando por ejemplo GeoGebra) nos permite estimar que elpunto de corte
del primer cuadrante está cerca de (0.7, 0.6). Así que tomamos ese valor como estimación inicial para lanzar
el método:
35
options(digits=15)
multiroot(F, start = c(0.7, 0.6))
##
##
##
##
##
##
##
##
##
##
##
$root
[1] 0.692081127846604 0.624782061427747
$f.root
[1] -1.75507809618836e-07
1.06462065385671e-07
$iter
[1] 3
$estim.precis
[1] 1.40984937502253e-07
Para el otro punto de intersección tomamos (1.2, −1) como estimación inicial:
options(digits=15)
multiroot(F, start = c(1.2, -1))
##
##
##
##
##
##
##
##
##
##
##
$root
[1] 1.260739706235216 -0.974714370658071
$f.root
[1] 6.20581364074724e-11 1.12921227923835e-10
$iter
[1] 4
$estim.precis
[1] 8.74896821656534e-11
Puedes compararlas con las soluciones que se han obtenido antes con Wolfram Alpha.
Apenas hemos rozado la superficie de este tipo de problemas. Recomendamos al lector que, si ha de hacer
uso de ellos, más allá de lo ocasional, profundice en su estudio usando las fuentes incluidas en las referencias
que aparecen más abajo.
Bibliografía.
Las mismas del tema anterior y además:
• Como texto de referencia para Álgebra Lineal, Linear Algebra and Its Applications de G. Strang, ver
http://math.mit.edu/~gs/linearalgebra/
• Como introducción a los sisemas de ecuaciones polinómicas: Ideals, Varieties and Algorithms, de Cox,
Little, O’Shea. Ed. Springer.
• Numerical Recipes 3rd Edition, Press, Teukolsky, Vetterling y Flannery. Cambridge University Press
(2007).
36

Download Report