Àlgebra Lineal Numèrica

Adrián Rodrı́guez Muñoz

11 de febrero de 2019
Capı́tulo 1

Errores

1.1. Ejemplos de errores numéricos

1.1.1. Ejemplo 1: Errores en los datos

Consideramos el sistema de ecuaciones:

1x + 2y = 3
0,4999x + 1,001y = 1,5

Con solución x = 1, y = 1
Si cambiamos 0,499 → 0,500 para obtener el sistema

1x + 2y = 3
0,500x + 1,001y = 1,5

Las solución ahora es x = 3, y = 0. Un cambio de 10− 3 en uno de los parámetros ha llevado a

cabo un cambio del orden del 101 en la solución; este tipo de error recibe el nombre de error en
los datos.

1.1.2. Ejemplo 2: Error en las operaciones

Consideramos el calculo de las raı́ces del polinomio x2 − 18x + 1 usando aritmética de 4 cifras
√
decimales en el calculo de :
√
18 ± 182 − 4 · 1 · 1 √
x= = 9 ± 80
2·1
√
Si tenemos como valor para 80 ' 8,9442, obtenemos las soluciones:

2
Algebra Lineal Numérica 3

x+ = 9 + 8,9442 = 17,9442
x− = 9 − 8,9442 = 0,0558

Hemos obtenido la primera solución con 6 dı́gitos significativos, pero la segunda solo con 3;
√
esto es una perdida importante de información debida al error cometido en la operación al
truncar a 4 dı́gitos.

1.1.3. Ejemplo 3: Error en el algoritmo

R1
Consideramos el calculo de In := 0 xn ex−1 dx, para n = 0, 1, 2, ...
Algoritmo: Obtenemos una recurrencia a partir de la formula de integración por partes
Z 1 Z 1
n x−1
In = x e dx = [xn ex−1 ]10 − nxx−1 ex−1 dx = 1 − nIn−1
0 0

R1 1
Con el valor inicial I0 = 0 x0 ex−1 dx = [ex−1 ]10 = 1 − e ' 0,63212055
Ası́ calculamos I0 = 0,63212055, I1 = 0,367899, ..., I9 = −0,068480; lo cual es un resultado
absurdo ya que xn ex−1 es no negativo para x ∈ [0, 1], y por lo tanto su integral en ese intervalo
también. Aunque la recurrencia sea cierta matemáticamente, el algoritmo empleado es malo
numéricamente (error en el algoritmo).

1.1.4. Ejemplo 4: Error en el algoritmo (2)

Consideramos el cálculo de la derivada de una función f en el punto x0

f (x0 + h) − f (x0 )
f 0 (x0 ) = lim
h→0 h

Sin embargo, un lı́mite es imposible de calcular numéricamente, y por lo tanto estamos obligados
a tomar una aproximación:

f (x0 + h) − f (x0 ) ∆f (x0 )

f 0 (x0 ) ' =
h h

Para h pequeña; la teorı́a sobre lı́mites nos dice que a medida que h se acerca a 0, el valor de la
aproximación se acercará cada vez más al valor real.

Notación. ∆f (x0 ) := f (x0 + h) − f (x0 ) operador diferencia hacia adelante

Desafortunadamente la teorı́a vuelve a fallar en la practica: si calculamos el cociente para ∆fh(x0 )

para h sucesivamente pequeñas, el valor que obtenemos mejorara al principio pero a partir de
un momento volverá a empeorar.
4 Rodriguez Muñoz, Adrián

1.2. Tipos de errores

Tipo “redondeo”: eR
→ Errores en los datos (ej 1)
→ Error en las operaciones (ej 2)
Tipo “truncado”: eT
→ Error en el algoritmo (en el paso de continua a aproximación discreta/numérica, ej 3
y 4)

Error total := eR + eT
En general, el error de tipo redondeo sera alto con un paso pequeño y bajo con un paso más
grande, mientras que el error de tipo truncado sera bajo con pasos pequeños y grande con pasos
más grandes. Debido a que el error total es la suma de los dos, no existe una solución trivial
para el paso óptimo h∗ que deberı́amos emplear en los cálculos; es decir, en la práctica suele ser
inviable encontrar con total exactitud el minimizador del error total que la teorı́a asegura que
existe.

1.3. Representación digital de un número real

Idea. Todo x ∈ R admite una representación única en base b ∈ N>1 con dı́gitos (aj )j∈Z

x = · · · + ak bk + · · · a0 b0 + a−1 b−1 + · · ·

En general la secuencia de dı́gitos aj tendrá un numero infinito de valores distintos de cero; al

ser la capacidad de los ordenadores finita, es imposible representar todos los reales con precisión
exacta lo cual generará errores.
Definición 1.3.1. La representación normalizada (única) en coma flotante de x ∈ R en base b
con dı́gitos (dn )n∈Z + y exponente e ∈ Z es:

x = ±0.d1 d2 ... · be

Con la restricción de que d1 6= 0

Ejemplo 1.3.1. La representación normalizada de 172,456 en base 10 es 0,172456 · 103
Definición 1.3.2. La representación en coma flotante con n dı́gitos de un numero real x con
representación normalizada x = ±0, d1 ..dn ... · be es:

fl por truncado:
flT (x) = ±0.d1 ..dn · be

fl por redondeo: (
±0.d1 ..dn · be 0 ≤ dn < 2b
fl(x) =
±(0.d1 ..dn + b−n ) · be b
2 ≤x<b
Algebra Lineal Numérica 5

Supongamos que queremos almacenar fl(x) con la siguiente estructura (en base b = 2)

[∗] [ ∗ | ∗ | ∗ | ∗ | ∗ | ∗ | ∗ ] [∗] [ ∗ | ∗ | ∗ | ∗ | ∗ ]
|{z} | {z } |{z} | {z }
s1 d1 ,...,dn s2 e1 ,..,ek

Donde s1 es el signo de x, d1 , ..., dn son los dı́gitos de fl(x)1 , s2 és el signo del exponente y
e1 , ..., ek los dı́gitos del exponente de fl(x).
Sin embargo, este mecanismo de representación no es del todo eficiente:

Idea. Una primera mejora se basarı́a en darnos cuenta de que al ser la base 2 y d1 6= 0, que
necesariamente d1 = 1. Ası́ pues, no hay razón para guardarlo y malgastar memoria; esto nos
da un dı́gito más de precisión con el mismo espacio.

Idea. Una segunda mejora serı́a cambiar la manera en el que guardamos el exponente e. Defini-
mos e0 = e + no , con un no fijo, definido como el mı́nimo natural tal que e0 > 0 para todos los e
posibles en la representación original; observamos que la definición de e0 establece una biyección,
ası́ que guardando e0 obtenemos unı́vocamente e. Por ser e0 > 0 siempre, al guardarlo podemos
prescindir el bit asignado al signo de e y utilizarlo para guardar un dı́gito más de e0 /e. El valor
de no para el tipo float (32 bits, 23 para la mantisa, 8 para el exponente) es 127, y 1023 para el
tipo double (64 bits, 52 para la mantisa, 11 para el exponente).

Definición 1.3.3. El ε-máquina de un ordenador es el mı́nimo ε > 0 tal que fl(1 + ε) > 1

1.4. Errores en la representación de los números

Error absoluto y relativo

Definición 1.4.1. Sea x ∈ R un número real y x una aproximación de x

Definimos el error absoluto cometido en aproximar x con x como: ea (x) = x − x

ea (x) ea (x)
Definimos el error relativo cometido en aproximar x con x como: er (x) = x ' x

Nota. Los errores no son conocidos; manejaremos cotas superiores de ellos:

x = x + ea (x) |ea (x)| ≤ εa (x) Cota del error absoluto

x = x(1 + δ) |δ| ≤ εr (x) Cota del error relativo
|{z}
|er (x)|

Observación 1.4.1. x = x(1 + δ) = x + |{z}

δx
ea (x)

1
o, equivalentemente, la mantisa
6 Rodriguez Muñoz, Adrián

Errores cometidos al representar x en punto flotante

Sea x = ±0.d1 ...dn dn+1 ... · be .

Error cometido en punto flotante por trucando: La representación de x en coma flotante con n
dı́gitos por truncado es fl(x) = ±0.d1 ..dn . Entonces por definición tenemos que:

|ea (flT (x))| = |fl(x) − x| = dn+1 ... · b−(n+1) · be ≤ b · b−(n+1) · be = b−n+e = εa (flT (x))
ea (x) b−n+e
|er (flT (x))| = | | ≤ −1+e = b−n+1 = εr (flT (x))
x b

Error cometido en punto flotante por redondeo: Con un argumento similar al anterior vemos
que las cotas para el error absoluto y relativo al hacer punto flotante por redondeo son:

1
|ea (flR (x))| = b−n+e = εa (flR (x))
2
1 −n+1
|er (flR (x))| = b = εr (flR (x))
2
Definición 1.4.2. Usualmente consideraremos b = 10.

Si ea (x) ≤ 21 10−k , diremos que x tiene k cifras decimales correctas.

Si er (x) ≤ 12 10−k+1 , diremos que x tiene k dı́gitos significativos correctos.

En esta definición aparecen los 21 ya que por defecto consideramos coma flotante por redondeo
(al ser el error menor); en caso de utilizar coma flotante por truncado, las definiciones serı́an
idénticas pero no aparecerı́a el 12 .

Propagación de errores en los datos

En principio tenemos que ea (x ± y) = ea (x) ± ea (y). Sin embargo al manejar cotas, estamos
obligados a usar:

Definición 1.4.3.
εa (x ± y) = εa (x) + εa (y)

En el caso del producto, tenemos que er (xy) = xy(1 + δx )(1 + δy ) = xy(1 + δx + δy + δx δy ) '
xy(1 + δx + δy ); nos quedamos con la aproximación de orden uno de esta expresión, ya que si
δx , δy << 1,entonces δx δy << δx + δy . Ası́ pues tenemos que:

Definición 1.4.4.
εr (x · y) / εr (x) + εr (y)

Nota. Ojo con restar variables “cercanas”; se generaran errores relativos grandes:
ea (x − y) x y
er (x − y) = = er (x) − er (y)
x−y x−y x−y
Algebra Lineal Numérica 7

Si ahora tomamos el valor absoluto:

x y
|er (x − y)| ≤ | | · |er (x)| + | | · |er (y)|
x−y x−y

Vemos que si x,y son cercanos, x − y será muy pequeño y al dividir la cota del error incrementa
muchı́simo.
√
Ejemplo 1.4.1. Calcular las raı́zes del polinomio x2 − 18x + 1, con calculo de de 4 cifras
decimales (por truncado) Obtenı́amos:
√
x+ = 9 + 80 = 17,94442 6 dı́gitos significativos
√
x− = 9 − 80 = 0,0558 3 dı́gitos significativos (!)

Analizamos√este resultado en términos de la cota del error relativo en el calculo de x− : Sea

x = 9, y = 80, x = 9, y = 8,4492. Entonces tenemos que:

y 8,9442
er (x − y) ' | | · |er (y)| = · 10−5+1 ' 0.16 × 10−1 < 10−2+1
x−y 0,0588 | {z }
5 dı́g sig

Es decir, obtenemos con exactitud sólo ≈ 2 dı́gitos significativos.

R1
Ejemplo 1.4.2. Consideramos el cálculo de In = 0 xn ex−1 dx, n = 0, 1, 2, ..., y, para ese fin,
haremos uso de la recurrencia In = 1 − nIn−1 ,I0 = 1 − 1e ≈ 0, 63212055. Recordemos que para
I9 obtenemos un valor negativo (cuando la función es positiva en el intervalo de integración). Si
analizamos este problema en términos del error absoluto: ea (In ) / ea (1) + nea (In−1 ), es decir
ea (In ) ' n!ea (I0 ). El error absoluto aumenta factorialmente: este procedimiento no es bueno.

1.5. Fórmula de la propagación del error en los datos

Tenemos x valor real, x aproximación. Queremos calcular f (x), pero solo podemos tener una
aproximación f (x) = f (x).

1.5.1. Error absoluto y relativo de f (x) con cálculo perfecto de f

Supongamos que todo el error viene de los datos y que el cálculo de f es perfecto:

Error absoluto

f (x) = f (x + ea (x)) = f (x) + f 0 (x)ea (x) + θ2 (ea (x)) ' f (x) + f 0 (x)ea (x)
| {z }
ord. dos en ea (x)

∴ |ea (f (x))| ' |f 0 (x)||ea (x)|

8 Rodriguez Muñoz, Adrián

Error relativo
ea (f (x)) f 0 (x)ea (x) xf 0 (x)ea (x) xf 0 (x)
er (f (x)) = ' = = er (x)
f (x) f (x) f (x)x f (x)
0 0
∴ |er (f (x))| ≈ | xff (x)
(x)
||er (x)|, dónde el término | xff (x)
(x)
| recibe el nombre de factor de propaga-
ción/amplificación del error (fdp).
Ejemplo 1.5.1 (Seguridad de la operación raı́z cuadrada). Consideramos el cálculo de f (x) =
√ 1
x 2√ √
x: fdp = √ x
x
= 21 . ∴ |er ( x)| ' 12 |er (x)| Calcular la raı́z cuadrada es una operación segura.
1 1
Ejemplo 1.5.2 (Seguridad de la operación 1−x2
). Consideramos el cálculo de f (x) = 1−x2
: fdp
x 4x2 2
(1−x ) 4x2 1
= 1 = 1−x2
= 4(1 − 1−x2
); vemos que el fdp serı́a muy grande cuando x ' 1; en ese
1−x2
caso, calcular f de esta manera serı́a una operación insegura. En vez de calcularlo dividiendo,
podemos expandir en Taylor: f (x) ≡ 1 + x2 + x4 + x6 + ...
Ejemplo 1.5.3 (Efecto en el cambio de los coeficientes de un sistema lineal en la solución).
Consideramos el sistema:
ax + by = 3
cx + dy = 1,5

Con soluciones x = 3d−1,5b 1,5a−3c

ad−bc ,y = ad−bc . Supongamos que sabemos a = 1, b = 2, d = 1,001 de
manera exacta, pero solo tenemos una aproximación de c dada por c = 0,499; con error absoluto
er (c) ' 1 × 10−3 . Al resolver este sistema, nosotros calculamos x = 3d−1,5b ∼
ad−bc = f (c). El error
absoluto cometido en el cálculo de x esta dado por:
|ea (x)| = |ea (f (c))| ' |f 0 (c)||ea (c)|
(3 · 1,001 − 1,5 · 2) · 2
|f 0 (c)| = ' 6 × 103
(1 · 1,001 − 2 · c)2 |c=0,499
∴ |ea (x)| ' 6
Nota. Supongamos que queremos calcular f (x) ≡ f (x1 , ..., xn ), pero solo tenemos x = (x1 , ..., xn ).
Entonces el error cometido (suponiendo que el cálculo de f es perfecto) viene dado por:
f (x) = f (x + ea (x)) = f (x) + ∇f (x)t ea (x) + θ2 (ea (x))
Tal que:
n
X ∂f
ea (f (x)) ' ∇f (x)t ea (x) = (x)ea (xk )
∂xk
k=1

1.5.2. Error absoluto y relativo de f (x) con cálculo aproximado de f

Supongamos que ademas del error en los datos x = x + ea (x), obtenemos un error relativo
al calcular f acotado por δf . Entonces tenemos la siguiente expresión para el error abosluto
cometido:
0
ea (f (x)) ' |f (x)||ea (x)| + δf |f (x)|
Algebra Lineal Numérica 9

1.6. Análisis del error en el cálculo de raı́zes de un polinomio

Consideramos las raı́ces de un polinomio como una función implı́cita de sus coeficientes. Veremos
que es una relación muy sensible y mal condicionada.

1.6.1. Análisis de un polinomio general

Consideramos el polinomio general:

p(x) = an xn + · · · ak xk + · · · a1 x + a0

(•) Sabemos que las raı́ces de p(x) son α1 , ..., αn

(•) También sabemos que cada raı́z αl es función de los coeficientes: αl = αl (an , ..., a0 )

(•) Queremos estudiar ea (αl (an , ..., a0 )), er ((αl (an , ..., a0 )) l = 1 ÷ n

En general sabemos, de la fórmula del factor de propagación del error:

n
X dαl
ea (αl (an , ..., a0 )) ' ea (aj )
daj
j=1

En nuestro ejemplo supondrems que solo hay error en el coeficiente ak . Por lo tanto:
∂αl
ea (αl (ak )) ' ea (ak )
∂ak

∂αl
Necesitamos calcular ∂ak . Utilizaremos derivación implı́cita respecto ak sobre la siguiente igual-
dad:
an αln + · · · ak αlk + · · · a0 = p(αl ) = p(αl (ak ))

Para obtener:
∂αl ∂αl αk
αlk = p0 (αl (ak )) =⇒ = 0 l
∂ak ∂ak p (αl (ak ))
El numerador de la expresión obtenida es fácil de calcular; para calcular el denominador haremos
uso de la representación de un polinomio como producto de sus factores simples:
n
Y
p(x) = an (x − αl ) (x − αj )
j=1,j6=l
 
n n
Y d  Y
=⇒ p0 (x) = an (x − αj ) + an (x − αl ) (x − αj )
dx
j=1,j6=l j=1,j6=l
n
Y
∴ p0 (αl ) = an (αl − αj )
j=1,j6=l
10 Rodriguez Muñoz, Adrián

De manera que obtenemos la expresión deseada para el error absoluto del cálculo de αl debido
al error en ak :
αlk
∴ |ea (αl (ak ))| ' Qn |ea (ak )|
an j=1,j6=l (αl − αj )

Con respecto al error relativo:

ea (αl (ak )) αk−1 ea (ak )

er (αl (ak )) = ' Qn l |ak |
αl (ak ) an j=1,j6=l (αl − αj ) ak
ak αlk−1
= Qn |er (ak )|
an j=1,j6=l (αl − αj )

Hasta aquı́ en general. En nuestro caso particular de Wilkinson, tenemos n = 20, α1 = 1,

α2 = 2,. . . ,α20 = 20, k = 19 ⇐⇒ a19 = −210, l = 16 ⇐⇒ αl = 16. Calculamos el error
relativo:

(−210) · (1619−1 )
er (αl (ak )) = |er (a19 )| ≈ 0.31 × 1011 er (a19 )
1 · (15!4!(−1)4 )

Consecuencia: Suponemos que queremos calcular el nuevo αl (ak ) (ex. α16 (a19 )) con t dı́gitos
significativos. ¿Cuantos dig. sig. correctos deberı́a tener, aproximadamente, a19 ?

1 1
· 101−t ∼ 0.31 × 1011 · 101−s
2 2
∴1 − t = 12 − s ⇐⇒
s ∼ 11 + t

Es decir que si queremos un valor para la raı́z en el que t dı́gitos sean, con seguridad, correctos,
entonces necesitamos saber 11+t dı́gitos del coeficiente a19 (suponiendo que conocemos los otros
de manera exacta).

1.7. Problema: “No asociatividad de la suma en numérica”

Suponemos δ cota del error relativo en el almacenamiento (δ ∼ 12 101−t ≡ t cifras significativas).
PN 1 P1 1
Queremos analizar las diferencias entre: n=1 n2 vs n=N n2 . Calculamos primero el error
cometido al sumar una secuencia general de n términos:

Valor real Aproximación Cota del error absoluto

s1 = x1 s1 = x1 ε(s1 ) = εa (x1 ) = x1 δ
s2 = s1 + x2 s2 = s1 + x2 εa (s2 ) ' (2x1 + 2x2 )δ
s3 = s2 + x3 s3 = s2 + x3 εa (s3 ) ' (3x1 + 3x2 + 2x3 )δ
s4 = s3 + x4 s4 = s3 + x4 εa (s4 ) ' (4x1 + 4x2 + 3x3 + 2x4 )δ
sN = sN −1 + xN sN = sN −1 + xN εa (sN ) ' (N x1 + N x2 + (N − 1)x3 + (N − 2)x4 + · · · + 2xN )δ
Algebra Lineal Numérica 11

Derivación de las cotas (casos s1 ,s2 ; las demás se obtienen análogamente):

(1) :εa (s2 ) = εa (s1 ) + εa (x2 ) + (s1 + x2 )δ = x1 δ + x2 δ + (x1 + x2 )δ ' (2x1 + 2x2 )δ
(2) :εa (s3 ) = εa (s2 ) + εa (x3 ) + (s2 + x3 )δ ' (2x1 + 2x2 )δ + x3 δ + (x1 + x2 + x3 )δ = (3x1 + 3x2 + 2x3 )δ

Para minimizar el error, querrı́amos: x1 ≤ x2 ≤ x3 ≤ · · ·P≤ xN . En nuestro caso particular con

xj := j12 , este resultado nos indica que la suma inversa ( 1j=N j12 ) darı́a un resultado mejor.
Capı́tulo 2

Sistemas lineales

2.1. Conceptos básicos

2.1.1. Objetivo:

Queremos resolver el sistema lineal A~x = ~b con A ∈ Mn (R); ~b, ~x ∈ Rn ; con A y ~b dados y tales
que det A 6= 0 ∴ ∃! solución ~x ∈ Rn

2.1.2. Métodos basados en A−1

Si n = 1 : ax = b

1. Calculamos a−1 = 1
a

2. x = a−1 b = b
a

Si n ∈ N : A~x = ~b

1. Calculamos A−1
2. ~x = A−1~b
Nota. El único cálculo eficiente de A−1 involucra resolver el sistema matricial AX = In , o
equivalentemente los n sistemas lineales Ax~i = e~i , i = 1 ÷ n. Necesitamos cambiar de estra-
tegia (utilizar la fórmula A−1 = det1 A C T , dónde C es la matriz de cofactores, es terriblemente
ineficiente)

2.1.3. Método de Cramer / Determinantes

a1,1 . a1,i−1 b1 a1,i+1 . a1,n

.. .. .. .. .. .. ..
. . . . . . .
an,1 . an,i−1 bn an,i+1 . a,n
A~x = ~b ⇐⇒ ~xi = para i = 1 ÷ n
|A|

12
Algebra Lineal Numérica 13

Nota. El cálculo del determinante con la fórmula de Leibniz es muy costoso.

Teorema 2.1.1. Fórmula de Leibniz
X n
Y
det A = sgn(σ) aiσ(i)
σ∈Sn i=1

Numero de operaciones de la Fórmula de Leibniz:

n! sumandos =⇒ n! − 1 sumas
Para cada sumando tenemos n + 1 términos multiplicando =⇒ n productos
∴ Cada determinante: n! − 1 sumas, n!n productos

Número de operaciones en el Método de Cramer:

(n + 1)(n! − 1) sumas
(n + 1)n!n productos
n divisiones
∴ (n + 1)2 n! − 1 operaciones

Tiempo que ejecución de la resolución del sistema por el Método de Cramer: (en una máquina
de 1 GFlops)

n Número de operaciones Tiempo

5 4319 0.43 × 10−5 segundos
15 3.3 × 1014 3.9 dias
20 1.07 × 1021 3.4 × 104 años
100 0.95 × 10162 3.02 × 10145 años

2.2. Sistemas con solución trivial

2.2.1. Sistema diagonal

Un sistema diagonal tiene matriz A ≡ D = diag(d11 , ..., dnn ), es decir:

    
d11 x1 b1
.  .   . 
D~x = ~b ⇐⇒ 

..   ..  =  .. 
    
dnn xn bn

Algoritmo: xi ← dbiii , i = 1 ÷ n
Requisito: dii 6= 0 ∀i ∈ [n] ⇐⇒ det D = ni=1 dii 6= 0
Q
Número de operaciones. n divisiones
14 Rodriguez Muñoz, Adrián

2.2.2. Sistema triangular superior

Un sistema triangular superior tiene matriz A ≡ U = (uij ) con uij = 0 si i > j, es decir:
    
u11 . . u1n x1 b1
~
 . .. .
..
 .   . 
. .
U~x = b ⇐⇒ 

 .  =  . 
   

unn xn bn
Pn
bn bi − k=i+1 uik xk
Algoritmo: xn ← dnn , xi ← per a i = (n − 1) ÷ 1
uii Q
Requisito: uii 6= 0 ∀i ∈ [n] ⇐⇒ det U = ni=1 uii 6= 0
Número de operaciones. n2 operaciones

n divisiones
n2 −n
2 sumas
n2 −n
2 multiplicaciones

2.2.3. Sistema lineal con matriz llena general

Para resolverA~x = ~b distinguiremos dos tipos de métodos:

1. Métodos directos: Obtienen la solución exacta en un número finito de pasos, suponiendo

que en el cálculo no se produce error alguno, ej. Método de Gauss.

2. Métodos iterativos: Parten de una aproximación inicial de la solución y calculan una su-
cesión de aproximaciones sucesivas mejoradas, el lı́mite de la cual es la solución buscada.

No hay ningún criterio general para saber qué método usar, pero en la práctica se observa que el
método directo es útil cuando la matriz es densa (pocos ceros) y su dimensión es más pequeña
que 100, y el método iterativo va bien cuando la matriz es escasa (muchos ceros) y su dimensión
es más grande que 100. En casos mixtos (ej. pocos ceros y dimensión alta) cualquiera de los dos
puede ser mejor que el otro caso a caso, o incluso se podrı́a utilizar (y se utiliza en la práctica)
métodos especializados si el sistema es de un cierto tipo especial.

2.3. Métodos directos

2.3.1. Método de Gauss

     
 a11 x1

 + ··· + a1n xn = b1 a1,1 . . a1,n x1 b1
.. .. .. ..  . .. . .
..   ..  =  ... 
A~x = ~b ⇐⇒
    
. . . . ⇐⇒  . .
  .    

an1 x1 + ··· + ann xn = bn an,1 . . an,n xn bn

Algebra Lineal Numérica 15

Paso a SL triangular superior

Denotamos A(0) := A, b(0) := ~b tal que tenemos A(0) ~x = b(0) , y suponemos a11 6= 0. En el primer
(0)
ai1
paso, queremos eliminar x1 de les ecuaciones 2, ..., n. Definimos mi1 := a011
y hacemos:

(ecuación i) := (ecuación i) − mi1 (ecuación 1)

Para i = 2 ÷ n
o equivalentemente:
(1) (0) (0)
aij := aij − mi1 a1j
(1)
bi := b(0) (0)
i − mi1 b1

Para i = 2 ÷ n, j = 1 ÷ n, tal que obtenemos el siguiente sistema equivalente:


(0) (0) (0) (0)


 a11 x1 + a12 x1 + ··· + a1n xn = b1
 (1) (1) (1)
 0 + a22 x2 + ··· + a2n xn = b2

.. .. .. .. ..


 . . . . .
(1) (1) (1)

+ ··· +

 0 + an2 x2 ann xn = bn

(1)
En el segundo paso suponemos a22 6= 0, y queremos eliminar x2 de la ecuación 3, ..., n. Definimos
(1)
ai2
mi2 := a122
y hacemos:

(ecuación i) := (ecuación i) − mi2 (ecuación 2)

Para i = 3 ÷ n
o equivalentemente:
(2) (1) (1)
aij := aij − mi2 a2j
(2)
bi := b(1) (1)
i − mi2 b2

Para i = 2 ÷ n, j = 1 ÷ n.
(0) (1) (2) (k−1)
Seguimos este procedimiento para 3, ..., n − 1. Los elementos a11 , a22 , a33 , ..., akk se llaman
pivotes. Si todos los pivotes son no nulos, entonces el en paso n − 1 tenemos el siguiente SL
equivalente.

(0) (0) (0) (0)


a11 x1 + a12 x1 + · · · + a1n xn = b1
 (1) (1) (1)
a22 x2 + · · · + a2n xn = b2


A (n−1)
~x = ~b(n−1)
⇐⇒ .. .. ..


 . . .
(n−1) (n−1)


 ann xn = bn

Que es un sistema triangular superior.

16 Rodriguez Muñoz, Adrián

Nota.  
(0) (0)  (0) 
a1,1 x1 . . . . . . . . . . . . . . . . . . . . a1,n xn b1
 .. (1)
  . 
 .. 

 . . . . . . . . . . . . . . . . . . a2,n xn

  

(k−1) ..
  k−1 
bk 
ak,k xk . . . . . . . . . . . .
 
A(k) = b(k) =
  
(k) (k)
 b(k) 

 ak+1,k+1 xk+1 . . . ak+1,n xn 
  k+1 
 . 
 .. ..   . 

 . . 
  . 
(k) (k) (k)
an,k+1 xk+1 . . . an,n xn bn

Resolución del sistema triangular superior

bn
xn = dnn
(i−1) Pn (i−1)
bi − k=i+1 aik xk
xi = (i−1)
per a i = (n − 1) ÷ 1
aii

Pseudocódigo del método de Gauss (forma esquemática)

Algorithm 1 Método de Gauss

1: procedure Gauss(A, b) . Solución al sistema A~x = ~b
. Paso a SL triangular superior
2: for k ← 1 ÷ n do
3: for i ← k + 1 ÷ n do
4: mik ← aakkik

5: for j ← k + 1 ÷ n do
6: aij ← aij − mik akj
7: end for
8: end for
9: end for
. Resolución del SL triangular superior
10: xn ← abnn
n

11: for k ← n − 1 P ÷ 1 do
bk − nj=k+1 akj xj
12: xk ←
akk
13: end for
14: return x . x és la solució al sistema A~x = ~b
15: end procedure

(k)
Nota. En la implementación del método de Gauss en un computador, el aij que se obtiene
(k)
sustituye los aij originales, de la misma forma que los bi sustituyen los bi . Los mik , llamados
multiplicadores, se pueden guardar en el sitio de los aik , que ya no serán usados nuevamente
una vez se ha calculado el mik .
Algebra Lineal Numérica 17

Numero de operaciones

(i) Paso a SL triangular superior

• Obtener A(n−1) :
n(n−1)(2n−1)
- Multiplicaciones: 6
- Sumas: n(n−1)(2n−1)
6
n(n−1)
- Divisiones: 2
• Obtener b(n−1)
n(n−1)
- Multiplicaciones: 2
n(n−1)
- Sumas: 2

(ii) Resolución del SL triangular superior:

• n divisiones
2
• n 2−n sumas
n2 −n
• 2 multiplicaciones
4n3 + 9n2 − 7n 2n3
∴ Total: '
6 3

Tiempo de ejecución (máquina de 1 GFlops)

n Número de operaciones Tiempo

5 115 1.2 × 10−7 segons
15 2570 2.6 × 10−6 segons
20 11 325 1.1 × 10−5 segons
100 681 550 6.8 × 10−4 segons

2.4. Estrategias de pivote.

En la eliminación gaussiana se ha supuesto que los pivotes eran no nulos pero podrı́a no ser ası́.
Por ejemplo, si nos encontramos con el siguiente sistema de ecuaciones:

x1 + x2 + x3 = 1 

x1 + x2 + 2x3 = 2


x1 + 2x2 + 2x3 = 1 
mediante el proceso de escalonamiento por el método de Gauss obtendrı́amos

x1 + x2 + x3 = 1 

x3 = 1


x2 + x3 = 0 
18 Rodriguez Muñoz, Adrián

(1)
En este ejmplo a22 = 0, de manera que es necesario intercambiar las filas 2 y 3.

x1 + x2 + x3 = 1  

x2 + x3 = 0


x3 = 1 

(1)
Ası́, a22 6= 0 y se puede continuar.
Este procedimiento se sigue en general. Si en el paso i-ésimo (A(i−1) a A(i) )
 
a11 a12 . . . . . . . . . . . . . a1n
 

 a22 . . . . . . . . . . . . . a2n 
 .. .. 

 . . . . . . . . . . . . 
(i)
A = akk . . . . . . . ak,n 
 
 

 ak+1,k+1 · · · ak+1,n 
 .. .. .. 

 . . . 

an,k+1 · · · an,n

(i−1) (i−1)
Si aii = 0 entonces existe algún aji 6= 0 (ya que det A 6= 0) y entonces se intercambian las
ecuaciones i-ésima y j-ésima y se continua el proceso con el pivote no nulo.
Por otro lado, si el pivote es no nulo pero muy pequeño se pueden generar errores numéricos.
Por ejemplo, en el sistema: ! ! !
−10−5 1 x1 1
=
2 1 x2 0
que tiene como solución x1 = −0,4999975... y x2 = 0,999995, si suponemos que se trabaja con
−2 (1)
un ordenador de tres dı́gitos, entonces m21 = = −0,2 · 106 y a2 = 1 − (−0,2 · 106) =
10−5
(1)
0,000001 · 106 + 0, 2 · 106 = 0, 2 · 106 . Por lo tanto, b2 = 0 − (−0,2 · 106 ) = 0,2 · 106 . De manera
el sistema es ! ! !
−105 1 x1 1
=
0 0,2 · 106 x2 0,2 · 106
1 − x2
y, solucionándolo obtenemos que x2 = 1 y x1 = = 0.
−10−5
(0) (0)
El número a11 es muy pequeño, m21 es muy grande y no deja intervenir al elemento a22 en el
cálculo de a122 .
La estrategia es intentar evitar multiplicadores grandes. En el ejemplo se intercambian las dos
ecuaciones ! ! !
2 1 x1 0
=
−10−5 1 x2 1
Operando como antes se tienen las soluciones x2 = 1 y x1 = 0,5.
Algebra Lineal Numérica 19

Por tanto, para prevenir posibles errores numéricos se sigue la estrategia de pivotaje parcial. En
el paso k-ésimo se escoge r como el menor natural tal que
(k−1) (k−1)
ark = max aik , k≤i≤n

y se intercambian las ecuaciones k y r. Se plantea ahora otro ejemplo:

! ! !
1 10000 x1 10000
=
1 0,0001 x2 1
con las solución x1 = x2 = 0,999.... Imaginemos que trabajamos con aritmética de coma flotante
(0)
de tres decimales. Se toma, por pivotaje parcial, a11 = 1 como pivote y m21 = 1 y se obtiene
! ! !
1 10000 x1 10000
=
0 −10000 x2 −10000
por lo que x2 = 1 y x1 = 0. El problema está en que la matriz está “desequilibrada”.
Si se considera el mismo sistema lineal intercambiando al principio las ecuaciones 1 y 2, entonces
! ! !
1 0,0001 x1 1
=
1 10000 x2 10000
con soluciones x2 = x1 = 1 que es una buena solución.
(0)
Si primero se equilibra el sistema lineal de partida de manera que max aij = 1, es decir
! ! !
1 10000 x1 10000
=
1 0,0001 x2 1
,dividiendo la primera ecuación,
! ! !
10−4 1 x1 1
=
1 0,0001 x2 1
solucionamos el problema y podemos aplicar el pivotaje parcial con resultados satisfactorios.
Parece pues razonable equilibrar el sistema y después aplicar pivotaje parcial.
En la práctica no se equilibra el sistema explı́citamente. Sólo es necesario modificar la estrategia
de pivotaje parcial; este procedimiento recibe el nombre de pivotaje parcial escalonado. En
el paso k-ésimo se elige
(k−1)
a (k−1)
max ik donde si = max aij .
k≤i≤n s1 k≤i≤n

Ejemplo 2.4.1. Se tiene el sistema

    
2 3 4 x 3
   1  
 4 4 −3 x2  = 3
    
−2 3 −1 x3 1
20 Rodriguez Muñoz, Adrián

En el primer paso,

(0)
s1 = max a1j = max(2, 3, | − 1|) = 3
1≤i≤3
(0)
s2 = max a2j = max(4, 4, | − 3|) = 4
1≤i≤3
(0)
s3 = max a3j = max(| − 2|, 3, | − 1|) = 3
1≤i≤3

y se calcula
(0)  
ai1 2 4 2
max = max , , = 1.
1≤i≤3 si 3 4 3

Por tanto, se intercambian las ecuaciones 1o y 2o .

       
4 4 −3
x1 3 4 4 −3 x
  1
!
     
 2 3 4  x2  = 3 ∼ 0 1 1/2  x2  = 3
        3/2
−2 3 −1 x3 1 0 5 −5/2 x3

Ahora se calculan se nuevo

 
1
s1 = max 1, =1
2
 
−5
s2 = max 5, =5
2

y por tanto
 
1 5
max , = 1.
1 5
No se hace el intercambio de filas y se sigue. El sistema final es
    
4 4 −3 x1 3
    
0 1 1/2 x2  = 3/2
    
0 0 −5 x3 −5

y como consecuencia x3 = 1, x2 = 1, x1 = 1/2.

Propiedades del método de Gauss

Lema. Cada paso de la eliminación gausiana se puede expresar en forma matricial:

A(k) = G(k) A(k−1) , b(k) = G(k) b(k−1) K =1÷n−1

Algebra Lineal Numérica 21

Con  
1
 .
..

 
 
 
1
G(k)
 
=
 
 −mk−1,k 1 

 .. .. 

 . . 

−mn,k 1 · · · 1

Demostración. G(k) es la matriz elemental que resta la fila k multiplicada por −mi,k a las filas
i, i = k + 1 ÷ n.

Corol·lari (Descomposición LU ).
(n−1) (n−1)
A
| {z } = |G · G(1) G(0)} A(0)
· ·{z
TS UTI

Si notamos L̃ ≡ G(n−1) · · · G(1) G(0) , U ≡ A(n−1) , obtenemos A = A(0) como producto de una
matriz triangular inferior con unos en la diagonal con una matriz triangular superior (TS: Trian-
gular superior, UTI: Unitriangular inferior).

U = LA ⇐⇒ A = L̃−1 U = LU

Donde  −1  
1 1
 ..   .. 
−1
 −m
2,1 .  m
 2,1 . 
L = L̃ = . = .
  
 .. .. 
 .. .. 
 . 
  . 

−mn,1 . . −mn,n−1 1 mn,1 . . mn,n−1 1

Proposición 2.4.1. Los pasos de la eliminación de Gauss (sin intercambios o pivotaje) man-
tienen el valor del determinante.

Demostración. A(k) = G(k) A(k−1) =⇒ det A(k) = det G(k) det A(k−1) = 1 · det A(k−1)

Notación. Llamamos matrices principales de orden m de una matriz A a:

A[m]11 submatriz formada por las primeras m files y las primeras m columnas.

A[m]12 submatriz formada por las primeras m files y las ultimes n − m columnas.

A[m]21 submatriz formada por las últimas n − m files y las primeras m columnas.

A[m]22 submatriz formada por las últimas n − m files y las últimas n − m columnas.
Lema. Sea P = LA con P, A, L ∈ Mn×n (R), y L triangular inferior. Entonces P[m]11 =
L[m]11 A[m]11 , m = 1, 2, .., n
(k)
Corol·lari. En el algoritmo de Gauss, det A[m]11 = det A[m]11 ∀k = 1, .., n − 1,∀m = 1, .., n
22 Rodriguez Muñoz, Adrián

Demostración. Utilizando el lema anterior y una demostración análoga a la proposición de la

conservación del determinante de la matriz entera.

Proposición 2.4.2. Sea A ∈ Mn (R) regular. Entonces podemos aplicar el algoritmo de Gauss
(i−1)
sin pivotaje, o equivalentemente: aii 6= 0 ⇐⇒ det A[m]11 6= 0 ∀m = 1, .., n

Demostración. Tenemos que:

 k−1 
(k−1) (i−1) (k−1)
Y
det A[k]11 = det A[k]11 = aii akk
i=1

De donde se deduce el enunciado

Nota. Este resultado nos da una condición necesaria y suficiente para que podamos aplicar el
método de Gauss sin intercambios, pero es difı́cil comprobar “a priori” la regularidad de matrices
generales. No obstante, es útil para matrices de tipos especiales (simétricas definidas positivas y
estrictamente diagonales dominantes) que aparecen, por ejemplo, en la resolución numérica de
EDPs.

Método de Gauss para matrices simétricas

Si A es simétrica, podemos reescribir el algoritmo de Gauss con aproximadamente la mitad e

(k)
memoria y número de operaciones ya que el algoritmo mantiene la simetrı́a i.e. A[k]22 es simétrica
para k = 1, .., n − 1.
(k)
Proposición 2.4.3. A[k]22 es simétrica para k = 1, .., n − 1.

Demostración. Demostremos por inducción:

(0)
Caso base: Tenemos por hipótesis que A[0]22 es simétrica
(k−1) (k)
Paso inductivo: Supongamos que A[k−1]22 es simétrica. Vemos que A[k]22 también lo es: en
el paso k del algoritmo de Gauss hacemos:
(k−1) (k−1) (k−1)
(k) (k−1) aik (k−1) (k−1) a (k−1) (k−1) ajk (k−1) (k)
aij ← aij − a
(k−1) kj
= aji − ki
(k−1)
a jk = aji − a
(k−1) ki
→ aji
akk | {z } akk | {z } akk
HI | {z } HI
HI

Nota. Por esta razón sólo es necesario calcular y guardar los coeficientes de la parte superior
y de la diagonal de A(k) , k = 1, .., n − 1. La parte inferior tiene 0’s o bien coeficientes que
conocemos por simetrı́a. Ası́ pues pues el algoritmo se escribe:
Observación
" # 2.4.1. Si A es simétrica, no podemos asegurar que no será necesario el pivotaje.
0 1 (0)
Ej: (tenemos a11 = 0)
1 1
Algebra Lineal Numérica 23

Método de Gauss para matrices simétricas definidas positivas

Definición 2.4.1. Una matriz simétrica A ∈ Mn (R) es definida positva si: xt Ax > 0 ∀x 6= 0

Teorema 2.4.1 (Criterio de Sylvester). Sea A una matriz simétrica. Entonces A  0 ⇐⇒

det A[k]11 > 0 para k = 1, ..., n

Corol·lari. Por la proposición anterior y el criterio de Sylvester, tenemos que podemos aplicar
Gauss sin pivotaje a las matrices simétricas definidas positivas.

Lema. Sea A simétrica definida positiva. Entonces:

(i) A[k]11 es simétrica definida positiva ⇐⇒ det A[k]11 > 0 para k = 1, ..., n

(ii) aii > 0 ∀ i = 1, .., n

(i−1)
(iii) aii > 0 ∀ i = 1, ..., n
(k)
(iv) Las sub-matrices A[k]22 k = 0, ..., n − 1 encontradas aplicando el método de Gauss son
también simétricas definidas positivas.

Demostración.

(i) Sylvester

(ii) (e(i) )t A(e(i) ) = aii > 0 (definición de definida positiva)

Qi (k−1) (k)
(iii) Inducción y utilizando la igualtad: k=1 akk = det A[k]11 = det A[k]11 > 0

(k)
(iv) Hemos visto anteriormente que A[k]22 es simétrica ∀ k = 0, ..., n−1. Por tanto solo hace falta
demostrar que se satisface el criterio de Sylvester. Podemos expresar el menor principal
(k)
que queramos de A[k]22 k = 0, ..., n − 1, como el cociente de un menor de A[∗]11 > 0 y
productos de pivotes de la diagonal > 0.

Método de Gauss para matrices estrictamente diagonalmente dominantes por filas

Definición 2.4.2. Una matriz A ∈ Mn (R) es estrictamente diagonalmente dominante por filas
(edd) si X
aii > |aij |
j∈[n]\i

Para todo i = 1 ÷ n

Lema. Si A es edd, entonces akk > 0 para todo k = 1 ÷ n

P
Demostración. |akk | > | j∈[n]\k akj | ≥ 0
24 Rodriguez Muñoz, Adrián

Proposición 2.4.4. A edd =⇒ A no singular.

Demostración. Consideramos Ax con x 6= 0 y queremos ver Ax 6= 0. Sea |xk | ≥ |xi | i = 1, .., n

P P
el elemento màximo de x. Tenemos que |xk | > 0, y que j∈[n]\k |akj xj | ≤ j∈[n]\k |akj ||xj | ≤
P P
|akk |||xk |. Por tanto tenemos que (Ax)k = j∈[n] akj xj = akk xk + j∈[n]\k akj xj tiene el mismo
signo que akk xk 6= 0, y por tanto (Ax)k 6= 0

Proposición 2.4.5. Todas les submatrices principales de una matriz edd son edd.
Corol·lari.

=⇒ Todos los menores principales son no nulos

=⇒ Podemos hacer Gauss sin pivotaje (intercambios)

2.5. Métodos de descomposición

2.5.1. Descomposición LU

Dada A no singular, queremos factorizar A = LU , donde L es una matriz unitriangular inferior

y U triangular superior.

{zL} det U = det U = u11 · · · unn

Nota. det A = |det
=1

Con esta descomposición podemos resolver sistemas de la forma Ax = b ⇐⇒ LU x = b

resolviendo los dos sistemas triangulares Ly = b,U x = y, con coste computacional de orden n2
cada uno.
Nota. La utilidad de la descomposición LU viene de la necesidad de resolver muchos sistemas
con la misma matriz sin saber necesariamente los vectores coeficientes de los sistemas. Si utili-
zamos Gauss sobre todos por separado estarı́amos haciendo los mismos cálculos repetidamente,
lo cual es ineficiente. En vez de esto, hacemos la factorización LU una sola vez, y resolvemos
cada sistema con coste cuadrático.
Proposición 2.5.1. Si A = LU , con L unitriangular inferior y U triangular superior, entonces
l, U son únicas.

Demostración. Si tenemos A = L1 U1 = L2 U2 =⇒ L−1 −1

2 L1 = U2 U1 . Al saber que la inversa de
matriz unitriangular inferior es unitriangular inferior, y la inversa de una triangular superior es
triangular superior, tenemos que L−1 −1
2 L1 = U2 U1 = In =⇒ L1 = L2 ,U1 = U2 .

Nota. Si todos los menores de A son no nulos, entonces podemos aplicar Guass sin pivotaje y
obtener A = LU . Pero esto no pasa en general, puede ser necesario hacer pivotaje; en estos caso
como hacemos LU ?
Definición 2.5.1. Diremos que una matriz P ∈ Mn (R) es una matriz de permutación si
sus filas o columnas son una reordenación de la identidad.
Algebra Lineal Numérica 25

Ejemplo 2.5.1 (Intercambio de 

filas: 0
  A = P A).   
0 0 1 a a12 a13 a31 a32 a33
   11   
1 0 0 a21
  a22 a23  =
 
a11 a12 a13 

0 1 0 a31 a32 a33 a21 a22 a23

Es decir, si P ≡ (ai,∗ := eσ(i) ), P A ≡ (a0i,∗ := aσ(i),∗ )

Ejemplo 2.5.2 (Intercambio de columnas: 0

   A = AP
 ).  
a a12 a13 0 0 1 a12 a13 a11
 11    
 1 0 0 = a22 a23 a21 
a21 a22 a23     

a31 a32 a33 0 1 0 a32 a33 a31

Es decir, si P ≡ (a∗,i := eσ(i) ), AP ≡ (a0∗,i := a∗,σ(i) )

Nota. P −1 = P T

Teorema 2.5.1 (Descomposición LU con matrices de permutación). Sea A no singular cual-

quiera. Entonces ∃ P matriz de permutación,L unitriangular inferior y U triangular superior
tal que P A = LU .

Demostración. Aplicamos el método de Gauss con intercambios si hacen falta, y obtenemos la

siguiente igualdad.

U = A(n−1) = G(n−1) Pn−1,ln−1 G(n−2) · · · P2,l2 G(1) P1,l1 A(0)

Vemos que el producto de una matriz de permutación que intercambia solo dos filas por una
matriz unitriangular inferior atómica (que solo tiene entradas no cero en una sola columna excep-
tuando la diagonal), se puede expresar como el producto de otra matriz unitriangular atómica
por la misma matriz de permutación en orden inverso, si los ı́ndices de las filas intercambiadas
son superiores al de la columna no nula de la matriz atómica. Es decir:

Pr,lr Uk = Ũk Pr,lr k < r, lr

Por tanto podemos ir aplicando esta propiedad sucesivamente a la igualdad anterior (del método
de Gauss) para obtener:

U = A(n−1) = G(n−1) Pn−1,ln−1 G(n−2) · P2,l2 G(1) P1,l1 A(0) ⇐⇒

U = A(n−1) = G(n−1) G̃(n−2) P · G˜(1) P P A(0)
n−1,ln−1 2,l2 1,l1 ⇐⇒
U = A(n−1) = |G(n−1) G̃(n−2) ˜ (n−3) ... P · · · P1,l1 A(0) ⇐⇒
{z G̃ } | n−1,ln−1
{z }
L̃ P
−1
U = L̃P A = L PA ⇐⇒
∴LU = P A
26 Rodriguez Muñoz, Adrián

Ejemplo 2.5.3 (Aplicaciones).

Cálculo del determinante: det A = (−1)s det L det U = (−1)s u11 ...unn On s es el número
de intercambios

Cálculo de la inversa: Equivalente a resolver LU x = P e(k) ⇐⇒ Lyk = P e(k) ,U xk = yk

k = 1, .., n

Cálculo de P

Partimos de A, y calculemos P1 , ..., Pn−1 tal que Pn−1 ...P1 A = LU

| {z }
P

Ejemplo 2.5.4. Aplicamos Gauss con pivotaje parcial escalonado al sistema siguiente:
   
1 2 2 4 1
   
1 2 1 2
 , b = 0
 
A= 3 2 0
 6
 7
 
4 2 1 1 2

Calcular P A = LU
Pivote 1: max(|1/4|, |1/2|, |3/6|, |4/4|) = 1 ∴ f 1 ↔ f 4 l1 = 4 p1 = (4, 2, 3, 1)
Eliminación 1:    
4 2 1 1 4 2 1 1
   
1 2 1 2
 7→ 1/4 3/2 3/4 7/4 

 
3 2 0 6 3/4 1/2 −3/4 21/4
   
1 2 2 4 1/4 3/2 7/4 15/4

Pivote 2: max(|(3/2)/(7/4)|, |(1/2)/(21/4)|, |(3/2)/(15/4)|) = 6/7 ∴ f 2 ↔ f 2 l2 = 2 p2 =

(4, 2, 3, 1)
Eliminación 2:
   
4 2 1 1 4 2 1 1
   
1/4 3/2 3/4 7/4 
 7→ 1/4 3/2 3/4 7/4
 

3/4 1/2 −3/4 21/4 3/4 1 1 2 
   
1/4 3/2 7/4 15/4 1/4 1/3 −1 2

Pivote 3: max(|(−1)/(14/3)|, |1/2|) ∴ f 3 ↔ f 4 l3 = 4 p3 = (4, 2, 1, 3)

Eliminación 3:
   
4 2 1 1 4 2 1 1
   
1/4 3/2 3/4 7/4
 7→ 1/4 3/2 3/4 7/4 
 

3/4 1/3 1 2  1/4 1 1 2 
 
 
1/4 1 −1 2 3/4 1/3 −1 20/3
Algebra Lineal Numérica 27

Resultados:
 
0 0 0 1
 
0 1 0 0
P=
  ↔ p = (4, 2, 1, 3)
1 0 0 0


0 0 1 0
 
1 0 0 0
 
1/4 1 0 0
L=
 
1/4 1 1 0
3/4 1/3 −1 1
 
4 2 1 1
 
0 3/2 3/4 7/4 
U=
 
0 0 1 2 
0 0 0 20/3

Cálculo del determinante: det A = (−1)2 · 4 · 3/2 · 1 · 20/3 = 40

h i
Resolución del sistema: Ax = b ⇐⇒ |{z}
P A x = P b = b̂ =⇒ Ly = b̂ = 2 0 1 7 ↔
h i LU h i
y = 2 −1/2 1 4/3 =⇒ U x = y ↔ x = 1 −1 −1 1

Cálculo de la inversa: LU xk = P e(k) ↔

 
1/20 −2/5 1/20 3/10
 
(−1)
−21/20 6/5 −1/40 −3/20
A = 
 7/10
 −3/5 −3/10 1/5 
3/20 −1/5 3/20 −1/10

2.5.2. Esquemas compactos

Objetivo: Encontrar la descomposición LU , suponiendo que no hay intercambios, de manera

directa. La filosofı́a de este tipo de algoritmos es que los métodos basados en Gauss realizan
muchas operaciones repetidamente sobre los coeficientes de la matriz, la cual cosa puede dar
lugar a errores propagados originados por errores de cálculo indeseados.

Descomposición LU general

Queremos A = LU , con L triangular inferior, U triangular superior:

28 Rodriguez Muñoz, Adrián

    
m11 u11 · · · u1n a11 · · · a1n
 . ..  .. ..   . .. .. 
 ..  .
. . = .
.  . . 

  
mn1 · · · mnn unn an1 · · · ann

El algoritmo general de los esquemas compactos trata de obtener una expresión analı́tica para
el producto de LU , y resolver de manera recurrente y exacta el sistema lineal que se obtiene
equivalente al de A. Observamos que en la ecuación general escrita arriba tenemos un total de
n2 + n variables y solo n2 ecuaciones (una para cada coeficiente de A). Esto nos permite fijar n
términos de cualquiera de las dos matrices, aunque normalmente seran los de la diagonal como
en los métodos particulares que se explican a continuación,

Método de Doolittle

Para el método de Doolittle fijamos los coeficientes de la diagonal de la matriz L a 1; obtenemos

la misma LU que con el método de Gauss:
    
1 u11 · · · u1n a11 · · · a1n
 . ..  .. ..   . .. .. 
 ..  .
.  . = .
.  . . 

  
mn1 · · · 1 unn an1 · · · ann

De la cual obtenemos:
LU Compacto - Método de Doolitle
1: procedure Doolittle(A) . Descomposición LU (sin pivotaje) de Doolittle
2: for k ← 1...n do
3: for j ← k...n doP
4: ukj ← akj − k−1p=1 mkp upj
5: end for
6: for i ← k + 1...n
 do P 
7: mik ← u1kk aik − k−1p=1 m u
ip pk
8: end for
9: end for
10: return L = (mij ), U = (uij ) . Tenemos A = LU ; L unitriangular inferior, U triangular
superior
11: end procedure

Método de Craut

Para el método de Craut fijamos los coeficientes de la diagonal de la matriz U a 1:

    
m11 1 · · · u1n a11 · · · a1n
 . .
..   ... ..  =  ... .. .. 
   
 .. . . . 
    
mn1 · · · mnn 1 an1 · · · ann
Algebra Lineal Numérica 29

De la cual obtenemos
LU Compacto - Método de Craut
1: procedure Craut(A) . Descomposición LU (sin pivotaje) de Craut
2: for k ← 1...n do
3: for i ← k...n doP
4: mik ← aik − k−1p=1 mip upk
5: end for
6: for j ← k + 1...n
 do P 
7: ukj ← mkk akj − k−1
1
p=1 m kp upj
8: end for
9: end for
10: end procedure
11: return L = (mij ), U = (uij ) . Tenemos A = LU ; L triangular inferior, U unitriangular
superior

Factoritzación d’una matriz simétrica

En el caso que A sea una matriz simétrica queremos obtener la siguiente factorización: A =
LDLT :
   T  
1 d11 1 a11 · · · a1n
 . .  . .
.. . .  =  ... .. .. 
  
 .. ..   .   .. . . 
     
ln1 · · · 1 dnn ln1 · · · 1 a1n · · · ann

De la cual obtenemos:

Factoritzación LDLT compacto para matrices simétricas

1: procedure LDLT (A) . Descomposición LDLT para matrices simetricas
2: for k ← 1...n doP
3: dkk ← akk − k−1 2
r=1 lkr drr
4: for i ← k + 1...n
 doP 
5: lik ← dkk1
aik − k−1 r=1 lir drr l kr
6: end for
7: end for
8: end procedure
9: return L = (lij ), D = (dij ) . Tenemos A = LDLT ; L unitriangular inferior, D diagonal

Nota. Hemos de tener dkk 6= 0 ∀k. Entonces LDLT x = b → Ly = b,LT x = D−1 y

Método de Cholesky - Factoritzación de una matriz simétrica definida positiva

Tomamos L = LD1/2 , y por lo tanto queremos A = LLT

30 Rodriguez Muñoz, Adrián

  T  
l11 l11 a11 · · · an1
 . .  . .
 =  ... .. .. 
 
 .. ..   .. .. . . 
    
ln1 · · · lnn ln1 · · · lnn an1 · · · ann

De la cual obtenemos:

LLT Compacto - Factoritzación de Cholesky

1: procedure Cholesky(A) . Descomposició LLT de Cholesky per a matrius simètriques
definides positives
2: for j ← 1...n
q do P
j−1 2
3: ljj ← ajj − k=1 ljk
4: for i ← j + 1, ..., n do 
1 Pj−1
5: lij ← ljj aij − k=1 lik lkj
6: end for
7: end for
8: end procedure
9: return L = (lij ) . Tenim A = LLT ; L triangular inferior

Nota.
Pj−1 2
Las cantidades ajj − k=1 ljk son todas positives.

Podemos guardar A, L con aproximadamente la mitad de memoria. Los lij se pueden ir

guardando en los sitios de los aij a medida que se vayan calculando.
√
El numero de operaciones es: 13 n3 + 23 n2 − 13 n + (coste de )n

2.6. Sistemas mal condicionados y análisis del error

Cuando resolvemos un sistema lineal numéricamente, la solución obtenida es en general una

aproximación de la solución exacta ya que intervienen los errores de redondeo en los cálculos
intermedios.

Definición 2.6.1. Diremos que un sistema lineal está mal condicionado si pequeños cambios en
la matriz de coeficientes y/o el término independiente producen grandes cambios en la solución.

Ejemplo 2.6.1. Consideramos el siguiente sistema lineal, y lo resolvemos por Gauss con pivotaje
parcial escalonado y aritmética de tres decimales:
(
0,832x1 + 0,448x2 = 1
0,784x1 + 0,421x2 = 0
Algebra Lineal Numérica 31

Con las soluciones obtenidas x1 = −506, x2 = 942. La solución exacta es x1 = −439, x2 = 817.
¿Cuál es el problema? Las rectas son casi paralelas; un pequeño cambio en la pendiente de una
de las dos envı́a la intersección muy lejos de su posición inicial. Si tomamos el siguiente sistema
lineal perturbado:
(
0,832x1 + 0,448x2 = 1
0,784x1 + (0,421 + ε)x2 = 0

Entonces si ε crece de 0 a 0,0012, las rectas se vuelven paralelas.

Querrı́amos tener una medida del mal o buen condicionamiento de un sistema. El tamaño del
determinante parece a priori un buen indicador: en el ejemplo anterior: det A = 0,00096. Pero
esta métrica puede ser engañosa, consideramos el sistema:

10−10 x1 + 0x2 = 0
−10
0x1 + 10 x2 = 0

con det A = 10−20 ' 0. Este sistema esta muy bien condicionado al ser las dos rectas perpendi-
culares.

Hemos visto que el determinante puede no ser un buen indicador, pero si lo será si la matriz
esta bien escalada. En efecto, en dimensión 2 un buen indicador del mal condicionamiento es el
ángulo que forman las dos rectas que definen el sistema lineal. Lo que haremos en la práctica
es tomar un indicador que es equivalente al ángulo: el área del paralelogramo de lado 1 definido
por las dos rectas. Es decir, tomamos el punto de intersección como el origen y calculamos el
área del paralelogramo formado por los vectores unitarios en las direcciones de la rectas. Vemos
que el área V viene dada por la siguiente expresión:

det A
q q
V = α1 = kA1,∗ k = a11 + a12 , α2 = k2, ∗k = a221 + a222
2 2
α1 α2

Ya que es conocido que el volumen n dimensional de n vectores de Rn (en el caso n = 2, el

volumen 2-dimensional es el área) viene dado por el determinante de la matriz que los tiene
como filas. Esta medida se extiende a n ecuaciones de manera natural:
 
a11 a1n
α1 ··· α1
| det A|  . q
.. .. 
V = = . . αi = Ai,∗ = a2i1 + · · · + a2in
α1 · · · αn  . . 

an1 ann
αn ··· αn

Tenemos que 0 ≤ V ≤ 1 (ya que todos los vectores fila son unitarios); si V ∼ 0 =⇒ A mal
condicionada, y V ∼ 1 =⇒ A bien condicionada.
Seguidamente veamos otra manifestación del mal condicionamiento: consideremos el sistema
lineal general Ax = b, y sea x la solución aproximada que obtenemos al resolverlo (sea cual sea
el método que usemos).
32 Rodriguez Muñoz, Adrián

Definición 2.6.2. Definimos el vector residual r como: r ≡ b − Ax = A(x − x); x − x = A−1 r

Si x ≈ x =⇒ r ≈ 0, x es una buena aproximación. La implicación en la otra dirección: r ≈ 0

=⇒ x buena aproximación, no es cierta en general.
Ejemplo 2.6.2. Consideremos el sistema lineal:
! !
1,2969 0,8648 0,8642
A= ,b =
0,2161 0,1441 0,1440

Obtenemos la solución x = (0,9911, −0,4870), y por lo tanto r = b − Ax = (−10−8 , −10−8 ). La

solución exacta es x = (2, −2). Esto resulta de la casi-singularidad de A =⇒ su inversa tiene
entradas muy grandes (ya que hemos de dividir por el determinante) =⇒ x − x = A−1 r puede
ser grande aunque r sea pequeño. En este caso tenemos:
!
−1 0.1441 × 108 −0.8648 × 108
A =
−0.2161 × 108 1.2969 × 108

2.7. Normas vectoriales y normas matriciales

Definición 2.7.1. Una norma en E R e.v. es una aplicación:

k∆k : E → R+
x 7→ kxk

que satisface:

(i) kxk = 0 ⇐⇒ x = 0E
(ii) αx = |α|kxk ∀α ∈ R,∀x ∈ E
(iii) kx + yk ≤ kxk + kyk ∀x, y ∈ E (desigualdad triangular)

En E = Rn , les normas tı́picas son:

kxk1 = ni=1 |xi |

P

Pn 2

1
kxk2 = i=1 xi
2
(norma euclidiana)
kxk∞ = maxi=1÷n {|xi 1}
1
kxkp = ( ni=1 |xi |p ) p p ≥ 1 (norma p de Hölder)
P

Definición 2.7.2. Una norma matricial es una norma en el espacio vectorial de matices
E = Mn (R) que, aparte de la definición usual, también cumple:

kABk ≤ kAk · kBk

Esto recibe el nombre de propiedad sub-multiplicativa.

Algebra Lineal Numérica 33

Definición 2.7.3. Diremos que una norma matricial es consistente con una norma vectorial
(que denotemos igual) si:

kAxk ≤ kAkkxk ∀A ∈ Mn (R), ∀x ∈ Rn

Definición 2.7.4. Dada una norma vectorial k∆k, se puede definir una norma matricial
subordinada consistente con ella a partir de:
 
kAxk
kAk := sup tq x 6= 0 = sup {kAxk tq kxk = 1}
kxk
kAxk
entonces tenemos ∀x 6= 0, kxk ≤ kAk ⇐⇒ kAxk ≤ kAkkxk

Ejemplo 2.7.1. Las normas matriciales subordinadas para las normas vectoriales habituales
son:
Pn
kAk1 = max ( i=1 |aij |) (máxima suma de columnas)
1≤j≤n
P 
n
kAk∞ = j=1 |aij | (máxima suma de filas)
1≤i≤n
p q
kAk2 = ρ(AT A) = max |ρi |, on ρ1 , ..., ρn son los diferentes VAP’s de AT A
1≤i≤n

Definición 2.7.5. El radio espectral d’una matriz B: ρ(B) := max |λi |, on {λi } son los diferentes
VAP’s de B.

Ejemplo 2.7.2.

1
D = diag(d11 , ..., dnn ) diagonal, antonces: kDk2 = max |dii |, D−1 2
=
i=1÷n min |dii |
i=1÷n

∀k∗k matricial subordinada: kIn k = max kIkxk

n xk
= max kxk
kxk = 1

Proposición 2.7.1. Resultados sobre normas matriciales subordinadas:

ρ(A) ≤ kAk, ∀k∗k matricial subordinada

k(I − A)uk ≥ (1 − kAk)kuk, ∀u ∈ Rn , A ∈ Mn (R)

Demostración.
|λ|kxk kλxk
ρ(A) = |λ| = kxk = kxk ≤ kAk

k(I − A)uk = ku − Auk ≥ kuk − kAuk ≥ kuk − kAkkuk = (1 − kAk)kuk

34 Rodriguez Muñoz, Adrián

2.8. Propagación de los errores en los datos de un SL

Supongamos que queremos resolver Ax = b (det A 6= 0,b 6= 0) pero que tenemos un error en los
datos, y por lo tanto tenemos A + (δA),b + (δb) en vez de A, b (A = A + (δA),b = b(δb). Entonces
al resolver el sistema lineal no obtnemeos la solcuión exacta x, sino una aproximación x + (δx)
(x = x + (δx)) tq: (A + (δA))(x + (δx)) = (b + (δb)). Con tal de obtener una cota para le error
cometido en nuestra aproximación e x hacemos los siguiente:

(A + (δA))(x + (δx)) = (b + (δb)) =⇒

Ax=b
Ax + A(δx) + (δA)(x + (δx)) = b + (δb) =⇒
(δx) = A−1 ((δb) − (δA)(x + δx)) =⇒
−1 −1 −1
k(δx)k = A ((δb) − (δA)(x + δx)) ≤ A kδbk + A k(δA)kkx + (δx)k =⇒
k(δx)k A−1 k(δb)k A−1 kδAkkx + (δx)k
≤ +
kxk kxk kxk

1 kAk
Sabemos que kbk = kAxk ≤ kAkkxk, y por lo tanto kxk ≤ kbk

k(δx)k A−1 k(δb)kkAk A−1 kAkk(δA)kkx + (δx)k (δb) k(δA)k kx + (δx)k

≤ + = kAk A−1 ( + )
kxk kbk kAkkxk kbk kAk kxk

Definición 2.8.1. Definimos el numero de condición de una matriz A como µ(A) = kAk A−1

En función del número de condición obtenemos:

 
k(δx)k k(δb)k kδ(A)k kx + (δx)k
∴ ≤ µ(A) +
kxk kbk kAk kxk

Interpretación: Si pensamos kx + (δx) ≈ kxkk, entonces µ(A) representa el factor máxima de

amplificación de los errores relativos de A y b. Dicho de otra manera, µ(A) es un indicador de si
el sistema lineal esta bien o mal condicionado i.e. pequeños cambios en A i/o b pueden producir
grandes cambios en la solución.

Proposición 2.8.1.

µ(I) = 1

µ(A) ≥ 1

Demostración.

µ(I) = kIn kkIn k = 1 · 1 = 1

µ(A) = kAk A−1 ≥ AA−1 = 1

Algebra Lineal Numérica 35

Ejemplo 2.8.1. Consideramos el sistema lineal Ax = b (que vimos que era muy mal condicio-
nado)
! !
1,2969 0,8648 0,8642
A= ,b =
0,2161 0,1441 0,1440
con
!
−1 0.1441 × 108 −0.8648 × 108
A =
−0.2161 × 108 1.2969 × 108

Entonces calculamos:

µ∞ (A) = kAk∞ A−1 ∞

= 2,1617 · 1.513 × 108 ≈ 3.3 × 108

Es decir, amplifica el error relativo como mı́nimo por ∼ 108 =⇒ SL muy mal condicionado.

2.9. Sistemas Lineales Sobre-determinados

Consideramos el sistema Ax = b, donde A ∈ Mm×n (R), m ≥ n, rank A = n. Queremos cal-
cular arg minx kb − Axk. Vimos en álgebra lineal que esto es equivalente a resolver el sistema
compatible determinado At Ax = At b
Ejemplo 2.9.1. Consideramos el problema de calcular (aproximar) el peso atómico del hidrógeno
i el oxı́geno con los pesos atómicos de los seis óxidos de nitrógeno: NO : 30,006,N2 O : 40,013,NO2 :
46,006,N2 O3 : 76,012,N2 O5 : 180,010,N2 O4 : 92,011. Expresamos este problema con el sistema
sobre-determinado:    
1 1 30,006
   
2 1  40,013 
  ! 
1 2 x  46,006 
   
   
2 3 y  76,012 
   
2 5 180,010
   

2 4 92,011
Resolvemos: AT Ax = AT b ! ! !
18 29 x 716,014
29 56 y 1302,161
∴ x = 14,0069 y = 15,99929

2.10. Aproximació discreta per mı́nims quadrats

2.10.1. Preliminars

Definición 2.10.1. Direm que un conjunt de funcions φ0 (x), ..., φn (x) es linealment independent
(l.i.) en un interval I ⊂ R si és satisfà: a0 φ0 + · · · + an φn ≡ 0 en I ⇐⇒ a0 = · · · = an = 0.
36 Rodriguez Muñoz, Adrián

xi x0 x1 ··· xN
yi y0 y1 ··· yN

Lema. Sigui {φj }j=0÷n un conjunt de funcions definides en un interval I ⊂ R,i x0 , .., xm ∈ I
una xarxa de punts. Aleshores son equivalents:

El conjunt {φj }j=0÷n és l.i. i la matriu Φ ≡ (φj (xi )) 1 ≤ i ≤ m, 0 ≤ j ≤ n te rang màxim.

El conjunt de vectors (φi (x0 ), ..., φi (xm )) ∈ Rm+1 0 ≤ i ≤ n és l.i.

Definición 2.10.2. Donat f, g definides en un interval I, i una xarxa de punts x0 , ..., xm ∈ I;

definim el producte escalar de f i g en x0 , ..., xm com:
m
X
(f, g) = f (xi )g(xi )
i=0

Atenció: Tot i que l’anomenem producte escalar, no es compleix la propietat de no degeneració

habitual. Es a dir: (f, f ) = 0 f ≡ 0 en I, només que les imatges de f son zero a tots els punts
de la xarxa. Tanmateix, aquest ”producte escalar”és bilineal, simètric i semidefinit positiu.

Definición 2.10.3. Anàlogament definim la norma d’una funció definida en el interval I respecte
una una xarxa de punts x0 , ..., xm ∈ I com:
m
p X 1
kf k = (f, f ) = ( f (xi )2 ) 2
i=0

Atenció: El mateix comentari d’abans; kf k f ≡ 0 en I, només que les seves imatges en la xarxa
son totes 0. Per tant, aquesta aplicació rep el nom de semi-norma. Tanmateix, la semi-norma
si que compleix les altres propietats de les normes sobre espais vectorials com ara:

kαf k = |α|kf k

kf k ≥ 0

kf + gk ≤ kf k + kgk

Teorema 2.10.1 (Teorema de Pythagoras per a funcions). Si (f, g) = 0, aleshores kf + gk2 =

kf k2 + kgk2 .

2.10.2. Recta de regressió

Suposem que tenim una taula de dades de la forma:

Notació: x ≡ (x0 , ..., xN ), y ≡ (y0 , ..., yN )
Ens agradaria aproximar la taula amb una funció; tenim dues opcions:

1. Interpolació: Obtenim una funció exacte

Algebra Lineal Numérica 37

(+) Error 0

(−) Car de calcular

(−) Si hi ha xi coincidents, aleshores no el podem fer servir.

2. Aproximació: Busquem minimitzar la norma de la diferencia entre la taula i un tipus de

funció; en el nostre cas, la emprarem k∗k2 i una recta.

(+) Més fàcil/barat de calcular

(−) Hi ha error
qP
N
Per tant, volem M ∗, B∗
∈ R tq M ∗, B∗
= arg min M, Bky − (M x + b)k2 = 2
i=0 (yi − (M xi + B)) .
En termes de sistemes sobre determinats, voldrı́em la millor solució al sistema:
   
1 x0 ! y0
.
 .. .. 
 B  .. 
 . 
M
=
  ⇐⇒ Aa = y
. 
1 xN yN

Resultats d’àlgebra lineal ens diuen que el vector a que minimitza kAx − bk és la solució al
sistema compatible determinat:
! ! !
1, 1 1, x B 1, y
= ⇐⇒ At Aa = At y
x, 1 x, x M x, y

Que anomenem sistema d’equacions normals (SEN). Més en general, si volem aproximar la
taula per una combinació lineal d’un conjunt funcions l.i. φ0 , ..., φm , f ≡ a0 φ0 + · · · + am φm ;
començarı́em amb el sistema sobre-determinat:
    
φ0 (x0 ) · · · φm (x0 ) a0 y0
 .. .. ..  .   . 
. . 
. . .   .  =  .  ⇐⇒ Aa = y
    

φm (xN ) · · · φm (xN ) am yN

I passarı́em al sistema d’equacions normals associat:

    
(φ0 , φ0 ) · · · (φ0 , φm ) a0 (φ0 , y)
 .. .. ..  .   . 
  ..  =  ..  ⇐⇒ At Aa = At y

 . . .    
(φm , φ0 ) · · · (φm , φm ) am (φm , y)

On els productes escalars estan definits sobre la xarxa de punts x0 , ..., xN .

Ejemplo 2.10.1. Volem calcular la recta de regressió sobre la següent taula amb el SEN.
38 Rodriguez Muñoz, Adrián

x 1 2 3 4 5 6 7 8 9 10 11
y 0 0.6 1.77 1.92 3.31 3.52 4.59 5.31 5.79 7.06 7.17

Calculem directament els productes escalars involucrats:

(, ) φ0 = 1 φ1 = x
φ0 = 1 11 66
φ1 = x 66 506
y 41.04 328.05

De manera que el SEN és:

! ! !
11 66 a0 41,04
=
66 506 a1 382,05

∴ a∗0 = −0,7314, a∗1 = 0,7437, Recta: y ≡ −0,7314 + 0,7437x

Comentari sobre les equacions normals: Tenim dues opcions per resoldre aquest sistema:

Mètode 1:

(•) Calcular M = At A, b = At y
(•) Fer Cholesky sobre M = Lt L
(•) Resoldre els sistemes triangulars Lu = b,Lt a = u

Mètode 2:

(•) Fer descomposició A = QR, amb Q ortonormal (per l’esquerra) i R triangular superior
no singular.
(•) Resoldre el sistema triangular Ra = Qt y (Que s’obté de: At Ax = At y ⇐⇒
Rt Qt QRa = Rt Qt y ⇐⇒ Rt Ra = Rt Qt y ⇐⇒ Ra = Qt y)

El segon mètode serà en general més ràpid i per tant el que farem servir. Altres aspectes positius
de la descomposició QR es que el nombre de condició de R es el mateix que el de A; i per tant
el mal condicionament del sistema es manté constant al fer QR (el qual no es cert amb altres
mètodes com ara LU o Cholesky; i per tant, probablement empitjora)

2.11. Descomposició QR
Sigui A ∈ Mm×n (R), m ≥ n, i de rang màxim n. Volem calcular la descomposició A = QR
amb Q ∈ Mm×n (R) ortogonal per l’esquerra (Qt Q = I) i R ∈ Mm×n (R) triangular superior no
singular. Emprarem el mètode de Gram-Schmidt modificat:
Algebra Lineal Numérica 39

Mètode de Gram-Schmidt modificat

1: procedure QRDecompose(A) . Descomposició QR de la matriu A de dimensions m × n
2: Q←A
3: for i ← 1..n do
4: R(i, i) ← kQ(:, i)k
5: Q(:, i) ← Q(:, i)/R(i, i)
6: for j ← i + 1..n do
7: R(i, j) ← A(:, i), A(:, j)
8: Q(:, j) ← Q(:, j) − R(i, j) · Q(:, i)
9: end for
10: end for
11: end procedure
12: return Q, R . Tenim A = QR; Q esquerra ortogonal, R triangular superior

Ejemplo 2.11.1. Calculem la descomposició QR de la matriu:

 
2 1
 
A =  1 0


−2 4

Normalitzem la primera columna:

(1)
3r11 = a1 =3
1 (1)
q1 = a = (2/3, 1/3, −2/3)t
r11 1

(1)
Ortogonalitzem a2 :

(1)
3r12 = q1 , a2 = −2
(2) (1)
a2 = a2 − r12 q1 = (7/3, 2/3, 8/3)t

Normalitzem la segona columna:

(2) √
3r22 = a2 = 13
(2)  t
a2 7 2 8
q2 = = √ , √ , √
r22 3 13 3 13 3 13

I hem acabat amb:  

2 √7
 13
!
3 13 
3 −2
Q=  √2  R = √
3 3 13  0 13
8 √8
3 3 13
40 Rodriguez Muñoz, Adrián

2.12. Descomposició SVD

2.12.1. Teoremes de la SVD

Teorema 2.12.1 (Teorema Espectral Real). Sigui M ∈ Mn (R) simètrica (M t = M ), aleshores:

(i) Tots els valors propis de M son reals

(ii) Veps de vaps diferents son ortogonals
(iii) Existeix una base ortonormal de VEP’s de M
Teorema 2.12.2 (Teorema de la SVD: Descomposició reduı̈da). Sigui A ∈ Mm×n (R) amb
m ≥ n. Llavors existeix una descomposició de la matriu A de la forma: A = U ΣV T , tq:

U ∈ Mm×n (R) és esquerra ortogonal (U T U = In ). Els vectors columna de U ≡ (u1 , .., un )
s’anomenen vectors singulars per la esquerra.
V ∈ Mn×n (R) és ortogonal (V T V = V V T = In ). Els vectors columna de V ≡ (v1 , .., vn )
s’anomenen vectors singulars per la dreta.
Σ ∈ Mn×n (R) les entrades de Σ ≡ diag(σ1 , ..., σn ) s’anomenen valors singulars de la matriu
A. Endemés σ1 ≥ · · · ≥ σn ≥ 0; σ1 rep el nom del valor singular principal (el més gran).

Demostración. Considerem S = AT A ∈ Mn (R), simètrica semi-definida positiva. Per el teorema

espectral, existeix una base ortonormal de VEP’s, amb VAP’s reals i no-negatius.

Ordenem els VAP’s de major a menor: λ1 ≥ · · · ≥ λn ≥ 0, i els VEP’s similarment

(Avi = λi a la nova ordenació).
Sigui r tq λ1 ≥ · · · ≥ λr > λr+1 = · · · = λn = 0
√
Definim ara σi := + λi i = 1 ÷ r; clarament σ1 ≥ · · · ≥ σr > 0.
1 m
Introduı̈m els vectors {ui := σi Avi }1≤i≤r ⊂ R . Aquests vectors formen conjunt ortonormal:
σ 2
σ
ui , uj = σi1σj Avi , Avj = σi1σj vit At Avj = σi σj j vit vj = σji vi , vj = δij . A continuació, ampliem els vec-
tors u1 , .., ur a un conjunt ortonormal de tamany n de Rm , que anomenem u1 , .., ur , ur+1 , .., un .
Veiem que es compleix:

Avj = σj uj , j = 1 ÷ r
Avj = 0, j = r + 1 ÷ n (ja que son VEP’s de VAP 0 de At A; son el nucli de A)

Això ho escrivim com:

   
a11 · · · a1n   u11 · · · un1  
1 1
 v1 · · · vn   ... σ1
 . ..   .. 
 .. . . 
  .. ..    .. 
. =

.
 .
 . ..   ... .. 
 
 .. . .  
  v1 · · · vn   σn
n n 1 n
an1 · · · ann um · · · um
Algebra Lineal Numérica 41

AV = U Σ ⇐⇒ A = U ΣV −1 = U ΣV t

Teorema 2.12.3 (Teorema de la SVD: Descomposició completa). Sigui A ∈ Mm×n (R) amb
m ≥ n. Llavors existeix una descomposició de la matriu A de la forma: A = Û Σ̂V T , tq:

U ∈ Mm×m (R) és ortogonal (U T U = U U T = In ).

V ∈ Mn×n (R) és ortogonal (V T V = V V T = In ).

Σ ∈ Mm×n (R) de la forma: " #

diag(σ1 , ..., σn )
Σ=
0(m−n)×n

Demostración. Procedim exactament igual que a la descomposició reduı̈da, excepte que ampliem
el conjunt {u1 , .., ur } a una base ortonormal de Rm . Aleshores tenim:
    
a11 · · · a1n   u11 · · · · · · um 1 σ1
 . .  v11 · · · vn1  . .  .. 
 .. ..    .. ..   .
  .. .. 

=
   
 . ..   . .  .

 ..
 
. ..   
σn 
 .  v1 · · · vn  . . 
 
n n ..
an1 · · · ann u1m · · · · · · um m 0 . 0

Corol·lari (Rang de la matriu A). Sigui A = Û Σ̂V T la descomposició SVD completa de A, i

r ∈ N tq σi > 0 si i ≤ r, i σi = 0 si i > r. Aleshores rank A = r.

Demostración. Per la SVD completa tenim: A ∼ Σ̂ (son semblants; en bases adients son iguals).
Per tant rank A = rank Σ̂ = r.
" #
4 3
Ejemplo 2.12.1. Calculem la SVD de A =
8 6
" #
80 60
• S= AT A =
60 45

• VAP’s de S: λ1 = 125,λ2 = 0

• VEP’s de S:
1 = [(4, 3)], v = (4/5, 3/5)
→ E125 1

→ E01 = [(3, −4)], v2 = (3/5, −4/5)

" #
4/5 3/5
• V ≡
3/5 −4/5
42 Rodriguez Muñoz, Adrián

" √ #
5 5 0
• Σ≡ ,r=1
0 0

• Per tot σi 6= 0 fem ui := σ1i Avi

√ √
→ u1 = σ11 Av1 = (1/ 5, 2/ 5)

• Ampliem {u1 } a un conjunt ortonormal de Rm=2 de tamany n = 2

√ √
→ u2 = (2/ 5, −1/ 5)
" √ √ #
1/ 5 2/ 5
• U≡ √ √
2/ 5 −1/ 5
 
1 2
 
Ejemplo 2.12.2. Calculem la SVD de la matriu A =  2 2 

2 4
!
9 8
• S = AT A =
8 9

• VAP’s de S: λ1 = 17,λ2 = 1

• VEP’s de S:
√ √
1 = [(1, 1)], v = (1/ 2, 1/ 2)
→ E17 1
√ √
→ E11 = [(1, −1)], v2 = (1/ 2, −1/ 2)
√ √ !
1/ 2 1/ 2
• V ≡ √ √
1/ 2 −1/ 2
√ !
17 0
• Σ≡
0 1

• Per tot σi 6= 0 fem ui := σ1i Avi

√ √ √
→ u1 = σ11 Av1 = (3/ 34, 4/ 34, 3/ 34)
√ √
→ u2 = σ12 Av2 = (−1/ 2, 0, 1/ 2)
 √ √ 
3/ 34 −1/ 2
 √ 
• U ≡ 4/ 34
 0 
√ √ 
3/ 34 1/ 2

• Per obtenir descomposició completa, ampliem u1 , u2 a una base ortonormal de Rm=3 :

√ √ √
→ u3 = (2/ 17, −3/ 17, 2/ 17)
Algebra Lineal Numérica 43

 √ √ √ 
3/ 34 −1/ 2 2/ 17
 √ √ 
• Û ≡ 
4/√34 0 −3/ 17
√ √ 
3/ 34 1/ 2 2/ 17
√ 
17 0
 
• Σ̂ = 
 0 1
0 0

2.12.2. Comentaris/propietats de la SVD

SVD de una matriu amb mes files que columnes: Si A ∈ Mn×m , amb n ¡m; aleshores
prenem B = AT ∈∈ Mn×m amb m > n i calculem la seva SVD: B = Ũ Σ̃Ṽ T , i obtenim
A = B T = Ṽ Σ̃T Ũ = U ΣV T . Atenció: En aquest cas U ∈ Mn×n ,Σ ∈ Mn×n i V ∈ Mm×n . El que
hem fet es equivalent a considerar la matriu AAT i fer un procés anàlog al prèviament vist.
VAP’s i VEP’s de AT A: Els VAP’s de la matriu semidefinida positiva AT A son σ12 , .., σn2 ; i
una base ortonormal de VEP’s corresponents és v1 , .., vn .

Demostración. Per construcció en la demostració del teorema de la SVD.

VAP’s i VEP’s de AAT : Els VAP’s de la matriu semidefinida positiva AAT son σ12 , .., σn2 i
m − n zeros (suposem A ∈ Mm×n amb m > n). Una base ortonormal de VEP’s corresponents
és u1 , .., un .

Demostración.
AAT = Û Σ̂V T V Σ̂T Û T = Û Σ̂Σ̂T UˆT
" #
Σ2 0n×n
= Û Û T
0n×n 0(m−n)×(m−n)

On Σ = diag(σ1 , .., σn )

Unicitat dels valors singulars: Els valors singulars σ1 , .., σn son únics.

Demostración. Els valors propis de AT A son únics.

VEP’s i VAP’s de la matriu H: Sigui A ∈ Mn×n (R); prenem A = U ΣV T , U = (u1 , .., un ),

Σ = (σ1 , .., σn ), V = (v1 , .., vn ). Considerem la matriu per blocs:
" #
0 AT
H ≡=
A 0

Els VAP’s de H son {+σ1 , −σ1 , ..., +σn , −σn }, amb base de VEP’s ortonormals corresponents
{ √12 (v1 , u1 ), √12 (v1 , −u1 ), ..., √12 (vn , un ), √12 (vn , −un )}
44 Rodriguez Muñoz, Adrián

Demostración. " # " #" # " #

vi 0 AT vi AT ui
H = =
ui A 0 ui Avi

En cas que σi = 0, aleshores AT ui = Avi = 0; i per tant:

" # " #
vi vi
H =0
ui ui

En cas que σi 6= 0
σi2
" # " # " # " #
AT ui AT σ1i Avi σi vi
vi
= = = σi
Avi Avi σi ui ui
" # " #
vi vi
Anàlogament H . Al ser {vi },{ui } ortogonals, el nou conjunt també es
= σi
−ui ui
q
ortogonal, i √12 [vi | ± ui ] = √12 kvi k2 + kui k2 = 1.

Rang,nucli i imatge de A: Si σ1 ≥ · · · ≥ σr > σr + 1 = · · · σn = 0

rank A = r

A = {vr+1 , ..., vn }

Im A = {u1 , ..., ur }

Demostración. rank A = r vist a corol·lari del teorema de la SVD. {u1 , ..., ur } és un s.e.v de
dimensió r contingut dins Im A, també de dimensió r ∴ {u1 , ..., ur } = Im A. {vr+1 , ..., vn } és un
s.e.v. de dimensió n − r contingut dins A, també de dimensió n − r ∴ {vr+1 , ..., vn } = A.

Solució analı́tica del problema de mı́nims quadrats: Si A té rang màxim (= n), llavors:
arg minkAx − bk = V Σ−1 U T b
x

2
Demostración. arg minkAx − bk = arg min U ΣV T x − b . Al ser rank A = n, tenim que Σ és
x x
invertible. Prenem ara la matriu Û ≡ [U |Ũ ] (de la SVD completa) ortonormal m × m. Aleshores:
! 2 ! 2
T 2 T T
2 UT T ΣV T x − U T b
U ΣV x − b = Û (U ΣV x − b) = (U ΣV x − b) =
Ũ T −Ũ T b
2 2
= ΣV T x − U T b + Ũ T b

Per tant: arg min U ΣV T x − b = arg min ΣV T x − U T b . Al ser Σ,V T invertibles, existeix solu-
x x
ció exacte (amb norma de la diferencia 0): ∴ arg minkAx − bk = (ΣV T )−1 U T b = V Σ−1 U T b
x
Algebra Lineal Numérica 45

2.13. Mètodes Iteratius

2.13.1. Resolució numèrica de mètodes iteratius

Considerem el sistema Ax = b. Volem una successió {x(k) }k≥0 tq lim x(k) = x. Com ho fem?
La idea es trobar una matriu B i un vector c tq: Ax = b ⇐⇒ x = Bx + C, i prendre la
successió x(k+1) = Bx(k) + c. Una manera d’establir un mètode iteratiu és descomposant la
matriu A = P − (P − A), amb P una matriu no singular anomenada precondicionador de A.
Aleshores fem: Ax = b ⇐⇒ P x − (P − A)x = b ⇐⇒ x = P −1 (P − A)x + P −1 b ⇐⇒
x = (I − P −1 A)x + P −1 b; i procedim: x(k+1) = (I − P −1 A)x(k) + P −1 b per k ≥ 0, x(0) donat.
En termes del format anterior: B ≡ (I − P −1 A),c ≡ P −1 b.
Nota. Comentaris sobre la descomposició amb matriu P

1. Podem escriure x(k+1) = x(k) + P −1 (b − Ax(k) ). Això ho interpretem com una correcció de
P −1 (b − Ax(k) ) feta la pas k + 1 per millor la solució.

2. Cost computacional:

Cada iteració: 2n2 operacions

Després de m iteracions: 2mn2 operacions.

3. Prenem P tq P −1 és fàcil de calcular: El mètode de Jacobi empra P diagonal, i el mètode

de Gauss-Seidel empra P triangular.

Ara ens preguntem sota quines condicions el mètode iteratiu convergirà.

Definición 2.13.1. Direm que el mètode iteratiu donat per (B, c) és convergent si ∀x(0) , la
successió {x(k+1) = Bx(k) + c} és convergent.
Proposición 2.13.1. El mètode iteratiu {x(k+1) = Bx(k) + c} és convergent ⇐⇒ ρ(B) < 1.

Demostración. Sigui x(0) la aproximació incial i ek := x(k) −x. Es compleix que x(k) = Bx(k−1) +c
i que x = Bx + c. Restant les dues equacions obtenim: ek = Bek−1 = · · · = B k e0 .

=⇒ Suposem que el mètode és convergent. Aleshores ek → 0. Prenem ara x(0) tq e0 sigui VEP
de B. Beo = λe0 =⇒ ek = B k e0 = λk e0 → 0. ∴ |λ| < ∴ ρ(B) < 1

⇐= Suposem que ρ(B) < 1. Volem veure que lim ek = 0. Sigui J la matriu de Jordan associada
a B. Si ρ(B) < 1, lim J k = 0 =⇒ lim B k = lim SJ k S −1 = 0 =⇒ B k e0 → 0 ∴ ek → 0.

Nota. Comentaris sobre mètodes iteratius convergents

1. Observem que quan menor sigui ρ(B), la convergència de ek cap a 0 és més ràpida. Definim
la velocitat de convergència com − ln(ρ(B))

2. Si en alguna norma (consistent), kBk < 1, aleshores ρ(B) ≤ kBk < 1, i per tant el mètode
és convergent.
46 Rodriguez Muñoz, Adrián

2.13.2. Fites de l’error en mètodes iteratius

Lema. Si kBk < 1, llavors:

kBk
(a) x − x(k) ≤ 1−kBk x(k) − x(k−1)

(b) x − x(k) ≤ kBkk x − x(0)

x − x(k−1)
Demostración. (a) x − x(k) = B(x − x(k−1) ) ≤ kBk x − x(k−1) ≤ kBk =
1 − kBk
kBk
x(k) − x(k−1) . Ja que x − x(k−1) = x − x(k) + x(k) − x(k−1) = B(x − x(k−1) ) +
1 − kBk
x(k) − x(k−1) ⇐⇒ (I − B)(x − x(k−1) ) = x(k) − x(k−1) =⇒ (1 − kBk) x − x(k−1) ≤
x(k) − x(k−1)
(I − B)(x − x(k−1) ) = x(k) − x(k−1) ⇐⇒ x − x(k−1) ≤
(1 − kBk)

(b) x − x(k) = B k (x − x(0) ) ≤ B k x − x(0) ≤ kBkk x − x(0)

2.13.3. Mètode de Jacobi

Si ai i 6= 0 ∀i, prenem P ≡ D := diag(a11 , .., ann ); obtenim el mètode de Jacobi: x(k+1) = BJ x(k) + cJ k≥0 ;
 

on BJ = I − D−1 A,cJ = D−1 b. Més explı́citament:

 
0 − aa12 · · · − aa1n  
11 11 b1

 − a21 .. ..  ( a11
0 . .  0 i=j  .  nb
BJ =  a. 22 ≡ . 
 .  ≡ aii
cJ =  i
 
 .. .. .. .
.. 
 a
− aijii i 6= j
 . .  bn
ann
− aann
n1
··· ··· 0

h i
(k+1) 1 P (k)
En format de components: xi = aii bi − i6=j aij xj

Nota. Una manera més natural d’introduir el mètode de Jacobi es aı̈llar xi de l’equació i-èssima
de Ax = b:  
1  X
xi = bi − aij xj 
aii
i6=j

Proposición 2.13.2. Si A és e.d.d, aleshores el mètode de Jacobi aplicat a A és convergent.

P   P 
aij 1
Demostración. kBJ k∞ = max1≤i≤n |
j6=i aii | = max 1≤i≤n aii j6=i |aij | < 1 per definició
de e.d.d.
Algebra Lineal Numérica 47

2.13.4. Mètode de Gauss-Seidel

Una manera natural de intentar millorar el mètode de Jacobi és utilitzar a cada iteració les
millors aproximacions dels xi ja calculats i.e:
 
i−1 n
(k+1) 1  X X (k)
xi = bi − aij x(k+1) − aij xj 
aii
j=1 j=i+1

(k+1) (k+1) (k+1) (k) (k)

És a dir, al MGS xi és calcula a partir de x1 , .., xi−1 , xi+1 , .., xn . Per escriure el mètode
de Gauss-Seidel en termes de la matriu P , considerem les següents matrius:
( ( (
aij i=j aij i>j aij i<j
D≡ L≡ U≡
0 i 6= j 0 i≤j 0 i≥j

(les parts diagonals, estrictament inferior i estrictament superior de la matriu A, tq A = D +

L + D, i prenem P ≡ D + L. De la igualtat vista P (x(k+1) − x(k) ) = b − Ax(k) veiem que:

(D + L)(x(k+1) − x(k) ) = b − (D + L + U )x(k) ⇐⇒ (2.13.2.13.0.)

Dx(k+1) + Lx(k+1) − Dx(k) − Lx(k) = b − Dx(k) − Lx(k) − U x(k) ⇐⇒ Dx(k+1) + Lx(k+1) = b − U x(k) ⇐⇒ x
(2.13.2.13.0.)

que expressa per components dona la igualtat prèviament vista. En format x(k+1) = BGS x(k) +
cGS , tenim BGS = I − (L + D)−1 A, cGS = (D + L)−1 b

Nota. Convergència del mètode de Gauss-Seidel:

Si A és e.d.d. per files o simètrica definida positiva, aleshores el MGS és convergent.

No hi ha resultats generals que assegurin que Gauss-Seidel convergeixi més ràpidament

que Jacobi (però a la pràctica en general s’observa que és més ràpid.)

Ejemplo 2.13.1. Considerem el sistema:

   
10 1 1 x 12
   1  
 1 10 1  x2  12
   
1 1 10 x3 12

(a) Demostreu que el mètode de Jacobi i de Gauss-Seidel son convergents per a

aquest sistema: La matriu A és edd per files, per tant hi ha convergència.

(b) Analitzeu el error comés per els dos algorismes si paressin a les següents ite-
racions (inclòs Jacobi, per Gauss-Seidel es fa exactament igual):
48 Rodriguez Muñoz, Adrián

Mètode de Jacobi:
   
x(5) = 1, 00032 1, 00032 1, 00032 x(6) = 0, 999936 0, 999936 0, 999936

Emprem la següent fita del error:

kBJ k∞ 1/5
x − x(6) ≤ x(6) − x(5) = 0,00141091 = 0,000096/
∞ 1 − kBJ k∞ 1 − 1/5
Ejemplo 2.13.2. Considerem el sistema:
   
4 −1 −1 x1 7
   
 4 −8 1  x2  −21
   
−2 1 5 x3 15

(a) Demostreu que el mètode de Jacobi i de Gauss-Seidel son convergents per a

aquest sistema: La matriu A és edd per files, per tant hi ha convergència.

(b) Calculeu vJ ,vGS

ρ(BJ ) = 0,3347; vJ = − log(ρ(BJ )) = 0,476

ρ(BGS ) = 0,125; vGS = − log(ρ(BGS )) = 0,9031

Veiem que Gauss-Seidel té major velocitat de convergència.

2.13.5. Mètode de relaxació

(k)
Suposem aii 6= 0 ∀i = 1 ÷ n. Si sumem i restem xi en el mètode de GS, tenim:
 
i−1 n
(k+1) (k) 1 bi −
X (k+1)
X (k)
xi = xi + aij xj − aij xj 
aii
j=1 j=i

(k+1) (k)
Interpretem això com que cada iterat xi s’obté de l’anterior xi més una correcció. El mètode
de relaxació consisteix en multiplicar la correcció per un parametre w, anomenat paràmetre
de relaxació, de manera que s’acceleri la convergència, i.e.
" Pi−1 (k+1) Pn (k) #
(k+1) (k) bi − j=1 aij xj − j=i aij xj
xi = xi +w
aii

En termes d’una matriu P , el mètode de relaxació s’obté amb: P ≡ w1 D + L = w1 (D + wL) (on

D i L son la part diagonal i part inferior estricte de A). Com a problema de punt fix, el mètode
de relaxació s’expressa com:

x(k+1) = Bw X (k) + cw Bw := (D + wL)((1 − w)D − wU ), cw := w(D + wL)−1 b

Algebra Lineal Numérica 49

Al ser Bw = P −1 (P −A) = w(D +wL)−1 ( w1 D +L−L−D −U ) = (D +wL)((1−w)D −wU ),cw =

P −1 b = w(D + wL)−1 b.
Veiem que a partir d’aquesta descomposició matricial s’obté el mètode iteratiu per components
vist abans:

1
P (x(k+1) − x(k) ) = b − Ax(k) ⇐⇒ ( D + L)(x(k+1) − x(k) ) = b − Lx(k) − Dx(k) − U x(k) ⇐⇒
w
(k+1) (k) (k+1)
Dx = Dx + w[b − Lx − (D + U )x(k) ] ⇐⇒
(k+1) (k) −1 (k+1) (k)
x =x + wD [b − Lx − (D + u)x ]

Proposición 2.13.3. Quins w son raonables per emprar en el mètode de relaxació?

(i) |w − 1| ≤ ρ(B)

(ii) El mètode es convergent =⇒ 0 < w < 2

Demostración. (i) Bw = (D+wL)((1−w)D−wU ) =⇒ | det Bw | = | det (D + wL)−1 || det[(1 − w)D − wU ] =

 
1
| a11 ···ann
||a11 · · · ann (1 − w)n | = |1 − w|n = |λ1 | · · · |λn ∴ |1 − w| ≤ max |λi | = ρ(B)

(ii) El mètode es convergent =⇒ |w − 1| ≤ ρ(B) < 1 =⇒ 0 < w < 2

Nota. Comentaris sobre el mètode de relaxació:

1. Si w = 1, obtenim GS

2. Si A és simètrica definida positiva, aleshores ρ(Bw ) < 1 ∀w ∈ (0, 2); el mètode de relaxació
es convergent ∀m ∈ (0, 2)

3. No hi ha cap resultat que ens digui com triar el w òptim, encara que hi alguns tipus de
matrius concretes per als quals si es coneix.
Ejemplo 2.13.3. Considerem el sistema
    
2 5 1 x 12
   1  
5 2
 5 x2  = 12
   

5 5 2 x3 12

Justifiqueu si el mètode de Jacobi es convergent: Calculem ρ(BJ )

 
0 −5/2 −1/2
 
BJ = −5/2 0 −5/2
−5/2 −5/2 0

Obtenim que PBJ (λ) ≡ det(BJ − λI) = λ3 + 55 75

4 λ − 4 . Veiem que P (1) < 0,P (2) > 0, per tant
per el teorema del valor mig ∃ VAP entre 1 i 2 =⇒ ρ(BJ ) > 1. Per tant el mètode de Jacobi
no convergeix per a aquest sistema.
50 Rodriguez Muñoz, Adrián

Ejemplo 2.13.4. Calculeu per a quins β convergeix per al següent sistema el mètode de
Jacobi/Gauss-Seidel: !
−10 2
A=
β 5
Trobem les matrius de Jacobi i GS
! !
0 1/5 0 1/5
BJ = BGS =
−β/5 0 0 −β/25

Calculem els valors propis:

β p p
PBJ (λ) = λ2 + , β/25 λ2 = − β/25
λ1 =
25
β 1 1 p 1 p
PBGS (λ) = λ2 − λ − , λ1 = (β − β 2 + 500) λ2 = (β + β 2 + 500)
25 5 50 50

Jc max{kλi k} = |β|/25 < 1 ⇐⇒ |β| < 25

1
p
GS max{kλi k} = 50 (|β| + |β|2 + 500) < 1 ⇐⇒ |β| < 20
Capı́tulo 3

Càlcul numèric de valors i vectors

propis

3.1. Conceptes bàsics

La motivació al darrere del càlcul de valors i vectors propis (abreviat sovint-ment a VAP i VEP
respectivament) és senzillament la seva rellevància i utilitat. Definim seguidament uns conceptes
bàsics de VAP’s i VEP’s

Definición 3.1.1. Donada una matriu A ∈ Mn (R), un vector propi v de valor propi λ és
v ∈ Rn \{0} tq Av = λv.

Definición 3.1.2. Donada una matriu A i un vector v, el quocient de Rayleigh està definit
v T Av
com: T .
v v
Proposición 3.1.1. Si A és una matriu i v és un vector propi, aleshores el seu quocient de
T
Rayleigh es igual al seu vector propi: vvTAv
v
=λ

v T Av v T (λv) T
Demostración. vT v
= vT v
= λ vvT vv = λ

Nota. El fet que una matriu no sigui diagonalitzable dificulta clarament el càlcul de VAP’s i
VEP’s

Teorema 3.1.1 (Teorema de localització dels VAP’s i VEP’s de Gerschgorin). Sigui A ∈ Mn (R),
n
Pnpropis de A pertanyen tots al conjunt R = ∪i=1 Ri , on Ri
aleshores els valors := {z ∈ C :
kz − aii k ≤ ri = j=1,j6=i kaij k}

Demostración. Sigui λ VAP de A, aleshores ∃x 6= 0 VEP, i.e. Ax = λx. Això és equivalent a:
n
X n
X
λxi = aij xj i = 1, .., n ⇐⇒ (λ − aii )xi = aij xj i = 1, .., n
j=1 j=1,j6=i

51
52 Rodriguez Muñoz, Adrián

Sigui xk = max kxj k. Aleshores:

n
X
k(λ − aii )xk k = aij xj =⇒
j=1,j6=k
Pn
aij xj n n
j=1,j6=k X kxj k X
kλ − aii k ≤ ≤ |aij | ≤ kaij k
kxk k kxk k
j=1,j6=k j=1,j6=k

El qual implica λ ∈ Rk ⊂ R

Ejemplo 3.1.1. Considerem la matriu:

 
−8 −2 −4
 
A=
−1 −4 2 

2 2 −10

Trobem la regió de Gerschgorin dels VAP’s i VEP’s:

−8
R1 ≡ z − 16 ≤ | − 2/16| + |4/16| = 3/8
−4
R2 ≡ z − 16 ≤ | − 1/16| + |2/16| = 3/16
−10
R3 ≡ z − 16 ≤ |2/16| + |2/16| = 1/4

3.2. Mètode de la potència

El mètode de la potència serveix per a calcular el VAP dominant (més gran en mòdul) d’una
matriu.

Proposición 3.2.1. Sigui A una matriu amb valor propi dominant, i.e. kλ1 k > kλ2 k ≥ · · · ≥ λ1 .
Aleshores el vap dominant λ1 és un nombre real.

Demostración. Es conegut que λ VAP ⇐⇒ λ el conjugat complex de λ també és VAP. Al

haver-hi VAP dominant, el conjugat ha de ser ell mateix ja que si no hi hauria dos VAP’s
diferents amb el mateix módul (que hem suposat que era el més gran).

Proposición 3.2.2. Sigui A una matriu amb VAP dominant |λ1 | > |λ2 | ≥ · · · ≥ |λ1 , i base de
VEP’s x1 , .., xn associada. Aleshores la seqüencia: zk+1 = Azk , z0 = α1 x1 + ... + αn xn qualsevol
amb α1 6= 0 és tal que:

(zk+1 )i
lim (zk )i = λ1 , ∀i = 1, .., n

lim αzkk = α1 x1
1

zk x1
lim kzk k = ± kxk
Algebra Lineal Numérica 53

Demostración. Veiem que:

zk = Ak z0 = Ak (α1 x1 + ... + αn xn ) = λk1 α1 x1 + ... + λkn αn xn

 k  k !
λ2 λn
= λk1 α1 x1 + α2 x2 + ... + αn xn
λ1 λ1

Aleshores:
  k+1     k+1 
Pn λj
α1 x1 + j=2 αj xj α1 x1 + O( λλ21 )
(zk+1 )i λk+1
1
λ1
i i
lim = lim k   k  = lim λ1    k  = λ1
(zk )i λ1 Pn λ λ2
α1 x1 + j=2 αj λ1j xj α1 x1 + O( λ1 )
i i
k
De les equacions obteses també obtenim lim λzkk = lim α1 x1 + O( λ2
λ1 ) = α1 x1 , lim kzzkk k =
1
α 1 x1
kα1 x1 k = ± kxx11 k (l’ultima igualtat presa en el sentit que el quocient convergeix a x1
kx1 k o − kxx11 k ).

Nota. Comentaris:

  
(zk+1 )i k
λ2
1. Velocitat de convergència del mètode de la potencia: (zk )i = λ1 1 + O λ1

z˜k
2. Podria passar que zk creixes ràpidament. Aleshores prenem les successions: yk := kz˜k k ,zk+1
˜ =
Az˜k Ayk
Ayk , z˜0 = z0 inicial. En termes d’una sola successió: zk+1
˜ = kz˜k k , yk+1 = kAyk k . Tenim
que:
(a) kyk k = 1
(b) kzk+1
˜ k = kAyk k ≤ kAkkyk k = kAk
˜ )
(zk+1 (zk+1 )i
(c) lim (yk )i
i
= lim (zk )i = λ1
(d) lim yk = VEP unitari de VAP λ1
3. Cas particular A simètrica: Podem accelerar la convergència a partir del quocient de Ray-
T
leigh: xxTAx
x
. En efecte, per ser A simètrica, existeix una base ortonormal de vectors propis
zkT Azk
x1 , .., xn . A partir d’aquı́ definim σk+1 = zkT zk
; z0 inicial donat. Aleshores tenim que:

zkT zk+1
σk+1 =
zkT zk
h P λj k ih P λj k+1 i
λ2k+1
1 α 1 x T +
1 ( λ1 ) α j x T
j α1 x 1 + ( λ1 ) α j x j
= h ih i
λ λ
λ2k α1 xT1 + ( λ1 )k αj xTj α1 x1 + ( λ1j )k αj xj
j
P P
1
P λ  !
α12 + ( λ1j )2k+1 αj2 λ2 2k

= λ1 P λ ∼ λ1 1 + O
α12 + ( λ1j )2k αj2 λ1

Obtenint el doble de velocitat de convergència (el logaritme del factor de decreixement).

54 Rodriguez Muñoz, Adrián

Ejemplo 3.2.1. Exemple de no convergència del mètode de la potència si no hi ha VAP domi-

nant: considerem
" # !
0 1 1
A= , z (0) =
−1 0 1

(z ) i
Es pot veure que (zk+1k )i
= 1, −1, 1, −1, ... i.e. no convergeix. Això es deu al fet que λ1 = 1,λ2 =
−1, amb |λ1 | = |λ2 |, i.e. no hi ha VAP dominant.

(zk )i
Ejemplo 3.2.2. Matriu A simètrica de VAP dominant λ1 = 3, 28868. En 3 iteracions: (zk−1 )i =
T
zk−1 Azk−1
30, 25987, σk = T
zk−1 zk−1
m = 30, 288662 i.e. 3 i 7 dı́gits significatius respectivament (observem
2 ∗ 3 ' 7).

3.2.1. Variants del mètode de la potència:

Mètode de la potència inversa

Sabem de l’àlgebra lineal que λ VAP de A ⇐⇒ 1/λ VAP de A−1 . Per tant aplicant el mètode de
la potència a A−1 (amb VAPS de mòdul | λ1n | > | λn−1 1
| ≥ · · · > | λ11 |; suposant que A té un VAP
”dominatı̈.e. de mòdul estrictament més petit) obtindrem el (invers multiplicatiu) del VAP de
mòdul mı́nim de A. Nota: Per implementar aquesta estratègia, podem calcular A−1 per aplicar
el mètode de la potència usual, o bé resoldre el sistema lineal Azk+1 = zk a cada iteració fent la
descomposició LU (que seria més eficient i segur).

Mètode de la potència desplaçada

Recordem que si λ VAP de A de VEP v, aleshores λ − q és VAP de A − qI de VEP v, q ∈ R.

Aleshores podem aplicar el mètode de la potència a A − qI (de VAPs λ1 − q, .., λn − q) amb q’s
astuts per accelerar la convergència: voldrı́em λλ21 −q λ2
−q < λ1 ; clarament suposant que λ1 −q,λ2 −q
segueixen sent els VAP’s de mòdul més gran i segon més gran. Per exemple, suposem que tenim
eigA = 14, 13, 12, 11,5 i fem A0 := A − 12I. Aleshores tenim eigA0 = 2, 1, 0, −0,5, on ara la
convergència es O(|1/2|k ) en comptes de O(|13/14|k ) (bastant millor).

Mètode de la potència inversa desplaçada

Per tal de refinar un VAP aproximat q, prendrem A0 = (A − qI)−1 i li aplicarem el mètode de la

potència. Obtenim el VAP µ dominant de A0 ; si hem escollit q prou bé, aleshores q−q
1
serà molt
gran i precisament aquest µ trobat. Aleshores podem obtindre una millor aproximació de q amb
q̃ := q + µ1 . Exemple: matriu A amb VAP proper a 4; aplicant potència a (A − 4I)−1 obtenim
1
µ = 24,919816; aleshores q̃ = 4 + 24,919816 = 4, 040129.
Algebra Lineal Numérica 55

3.3. Mètodes per obtenir la totalitat dels VAPs

En aquesta secció tractarem mètodes basats en transformacions de semblança: volem obtenir a
partir de la matriu inicial una matriu més senzilla (desde el punt de vista del càlcul de VAPs)
amb els mateixos VAP’s que A.

3.3.1. Mètode de Jacobi

Suposarem A simètrica. Per el teorema espectral sabem que existeixen n VAPs relas i una base
ortonormal, representada per una matriu ortonormal Q, tq Q−1 AQ = D, amb: cols(Q) = VEPs,
D = diag(eig(A)).

Definición 3.3.1. Anomenem Q−1 AQ transformació de semblança de A, i direm que A ∼

Q−1 AQ semblants (o similars) si tenen els mateixos VAPs. Prendrem tı́picament Q matriu
de rotació (un tipus particular de matriu ortonormal). Nota: Per el que fa a la definició de
semblança Q podria ser qualsevol matriu no singular, però només una matriu ortonormal manté
la simetria.

Definición 3.3.2. Escriurem Rpq (φ) := Matriu de rotació d’angler φ en pla 0 + [ep , eq ].


cos(φ) i = p, j = p

− sin(φ) i = p, j = q



Rpq (φ)ij = sin(φ) i = q, j = p

cos(φ) i = q, j = q





1i=j altrament


On 1i=j = 1 si i = j i 0 altrament.

Observem que quan fem A := Rpq (φ)T ARpq (φ), només canviem la fila p, fila q, columna p i
columna q.

Proposición 3.3.1. Definim Ak = QTk Ak−1 Qk , k ≥ 1 A0 = A inicial simètrica i Qk ortonormal

per tot k. Aleshores es cert que:

Ak té els mateixos VAP’s que A, i els VEPs estan relacionats per: x = Q1 ...Qk xk , si x i
xk son VEPs de A i Ak del mateix VAP.

Les matrius Ak son simètriques per k ≥ 1.

Demostración.

Definim Q̃ := Q1 · · · Qk , aleshores tenim que Ak = Q̃T AQ̃ i x = Q̃xk . Veiem que Ax =

Q̃Ak Q̃T Q̃xk = Q̃Ak xk = Q̃λxk = λx

ATk = Q̃T AT Q̃ = Q̃T AQ̃ = Ak , al ser A = A0 simètrica per hipòtesi.

56 Rodriguez Muñoz, Adrián

El mètode de Jacobi es basa en triar la seqüencia adequada {Qk } de matrius de rotació ortonor-
mals tq Ak → D matriu diagonal. En el mètode de Jacobi busquem a cada iteració el element
(k) (k+1)
apq més gran en valor absolut de Ak fora de la diagonal i triem Qk := Rpq (φ), tal que apq = 0.
(k) (k+1) 0
Simplifiquem la notació referint-nos a aij per aij , aij per aij , cos(φ) per c, sin(φ) per s i
tan(φ) per t. Per trobar φ fem el següent càlcul:

Es pot veure que

a0pq = apq (c2 − s2 ) + (aqq − app )sc
Imposant la igualtat amb zero i dividint per c2 obtenim
 
2 2 aqq − app
app (1 − t ) + (aqq − app )t = 0 ⇐⇒ 1 − t + t=0
apq
aqq −app
p
Definim 2η := apq , d’on obtenim t = −η ± η 2 + 1.
p
Si η ≥ 0, prenem t = −η − η 2 + 1 < 0, c = √ 1 > 0, s = ct < 0
1+t2
p
Si η < 0, prenem t = −η + η 2 + 1 > 0, c = √ 1 > 0, s = ct > 0
1+t2

La elecció del valor de t pres està feta per tal de no sumar coses de signes diferents, el qual podria
(k+1) (k+1)
donar lloc a errors numèrics. Observem que a cada iterat, introduı̈m dos zeros en apq =qp =
0, però podem introduir components no nul·les ales posicions que abans eren zero. Malgrat això,
es pot demostrar que la successió A0 , A1 , ..., Ak tendeix a D. Al la pràctica, pararem a un cert k
i tindrem Ak = QTk · · · QT1 A0 Q1 · · · Qk ≡ Q̃T A0 Q̃, i prendrem els elements de la diagonal de Ak
i les columnes de Q̃ com una aproximació dels VAPs i els VEPs.
En el cas que A siguin una matriu qualsevol (no simètrica), haurı́em de fer servir altres mètodes
(que no s’expliquen al curs). Per exemple, si només volem els VAPs podrı́em transformar la
matriu inicial A en una de semblant de tipus Hessenberg (triangular superior + subdiagonal),
per a la qual és fàcil trobar el polinomi caracterı́stic i les seves arrels mitjançant un mètode
recurrent. Si volguéssim tots els VAPs i VEPs, farı́em servir un altre mètode diferent; un d’usual
s’anomena el mètode QR.