Diapositivas Tema 2 Sistemas de ecuaciones lineales

SAVE OFFLINE

Lecci´ on 2: Sistemas de ecuaciones lineales

´ GAUSSIANA. ELIMINACION

  a11x1 + a12x2 + · · · + a1nxn = b1   Sistema de n ecuaciones  a21x1 + a22x2 + · · · + a2nxn = b2 ... ...  con n inc´ ognitas:    an1x1 + an2x2 + · · · + annxn = bn

    

a11 a12 a21 a22 ... ... an1 an2

Ax = b







x1 . . . a1n    . . . a2n     x2  =     ...   ...  ...  xn . . . ann

≡

b1 b2 ... bn

M = [A|b]

    

A es invertible ⇒ soluci´ on ´ unica x = A−1b

A es invertible ⇒ soluci´ on ´ unica x = A−1b Nunca se calcula A−1 (m´ as costoso y sensible a errores de redondeo)

A es invertible ⇒ soluci´ on ´ unica x = A−1b Nunca se calcula A−1 (m´ as costoso y sensible a errores de redondeo) Ejemplo: Resolver 10x = 0.9

Proceso de Eliminaci´ on Gaussiana: restar m´ ultiplos de una ecuaci´ on a las siguientes    a11x1 + a12x2 + · · · + a1nxn = b1       a x + a x + ··· + a x  21 1 22 2 2n n = b2 ... ...        an1x1 + an2x2 + · · · + annxn = bn 

(si a11 6= 0)

 a11x1 + a12x2 + · · · + a1nxn = b1     (2) (2) (2)  a22 x2 + · · · + a2n xn = b2 ... ...      (2) (2) (2) an2 x2 + · · · + ann xn = bn

          

Proceso de Eliminaci´ on Gaussiana: restar m´ ultiplos de una ecuaci´ on a las siguientes    M = 



a11 a12 · · · a1n b1 a21 a22 · · · a2n b2   ...  ... ... ... ...  an1 an2 · · · ann bn



(si a11 6= 0)

    

a11 a12 · · · a1n 0 ... 0

b1



(2)  (2) (2) a22 · · · a2n b2  ... ... ...  ...   (2) (2) (2) an2 · · · ann bn

Para eliminar los elementos a21, a31, . . . , an1:

Para eliminar los elementos a21, a31, . . . , an1:

Calculamos los multiplicadores: l21 = a21/a11, l31 = a31/a11, . . . , ln1 = an1/a11

Para eliminar los elementos a21, a31, . . . , an1:

Calculamos los multiplicadores: l21 = a21/a11, l31 = a31/a11, . . . , ln1 = an1/a11

Realizamos operaciones elementales de fila: F2 − l21 · F1,

F3 − l31 · F1, · · ·

Fn − ln1 · F1

Proceso de Eliminaci´ on Gaussiana: restar m´ ultiplos de una ecuaci´ on a las siguientes    a11x1 + a12x2 + · · · + a1nxn = b1       a x + a x + ··· + a x  21 1 22 2 2n n = b2 ... ...        an1x1 + an2x2 + · · · + annxn = bn 

(si a11 6= 0)

 a11x1 + a12x2 + · · · + a1nxn = b1     (2) (2) (2)  a22 x2 + · · · + a2n xn = b2 ... ...      (2) (2) (2) an2 x2 + · · · + ann xn = bn

          

Proceso de Eliminaci´ on Gaussiana: despu´ es de n − 1 etapas    a11x1 + a12x2 + · · · + a1nxn = b1       a x + a x + ··· + a x  21 1 22 2 2n n = b2 ... ...        an1x1 + an2x2 + · · · + annxn = bn 

 a11x1 + a12x2 + · · · + a1nxn = b1     (2) (2)  a(2) x + · · · + a x = b n 2 22 2n 2 . ... .  .     (n) (n) ann xn = bn

llegamos a un sistema triangular

          

Ejemplo: Resolver:    10x − 7y

= 7 −3x + 2y + 6z = 4   5x − y + 5z = 6 Soluci´ on x = 0, y = −1, z = 1

Algoritmo de sustituci´ on regresiva: xn = ebn/unn

Algoritmo de sustituci´ on regresiva: xn = ebn/unn xn−1 = (ebn−1 − un−1nxn)/un−1,n−1

Algoritmo de sustituci´ on regresiva: xn = ebn/unn xn−1 = (ebn−1 − un−1nxn)/un−1,n−1 ... x1 = (eb1 − u1,2x2 − · · · − u1nxn)/u11

Algoritmo de sustituci´ on regresiva: En MatLab: xn = ebn/unn for j = n − 1 : −1 : 1 xj = (ebj − uj,j+1xj+1 − · · · − ujnxn)/ujj end

Interpretaci´ on matricial del m´ etodo de eliminaci´ on:

[ A | b ] = [A(1)|b(1)]

···

[A(2)|b(2)]

[A(n−1)|b(n−1)]

···

[A(n)|b(n)] = [ U | eb ]

A = A(1)

A(2)

   (1) A = 

···

a11 a12 a21 a22 ... ... an1 an2

A(n−1)



A(n) = U

. . . a1n . . . a2n   =A ...  ...  . . . ann

A = A(1)

A(2)

   (n) A = 

···

A(n−1)

u11 u12 . . . u1n 0 u22 . . . u2n ... ... ... ... 0 . . . 0 unn

A(n) = U

   =U 

Relaci´ on entre la matriz inicial A = A(1) y la matriz triangular U = A(n)

Relaci´ on entre la matriz inicial A = A(1) y la matriz triangular U = A(n)

A = LU (factorizaci´ on LU de la matriz A)



1

0 ...  l ...  21 1 L =  .. . . ... .  . ln1 . . . ln,n−1

0 ... 0 1

   , 

   U = 

u11 u12 . . . u1n 0 u22 . . . u2n ... ... ... ... 0 . . . 0 unn

    

Requisito para poder realizar el proceso de eliminaci´ on: Los pivotes deben ser distintos de cero

Requisito para poder realizar el proceso de eliminaci´ on: Los pivotes deben ser distintos de cero Si nos encontramos con alg´ un pivote u11, u22, . . . , unn que se anule ⇒ no es posible la factorizaci´ on A = LU

Si alg´ un pivote es nulo ⇒ hay que permutar filas en A (permutar ecuaciones en el sistema)

Si alg´ un pivote es nulo ⇒ hay que permutar filas en A (permutar ecuaciones en el sistema) Una matriz de permutaci´ on P es una matriz que se obtiene intercambiando filas en la matriz identidad

Una matriz de permutaci´ on P es una matriz que se obtiene intercambiando filas en la matriz identidad Propiedades

Toda matriz de permutaci´ on P cumple P −1 = P T

Si P es una matriz de permutaci´ on y A es una matriz arbitraria, entonces P A coincide con la matriz A con las filas ordenadas seg´ un indica la matriz P

Si A es una matriz cuadrada con det(A) 6= 0, existe una matriz de permutaci´ on P tal que P A = LU

Si A es una matriz cuadrada con det(A) 6= 0, existe una matriz de permutaci´ on P tal que P A = LU

Podemos elegir distintas matrices de permutaci´ on P (⇒ existen diversas factorizaciones P A = LU )

Una de ellas se obtiene con el comando lu de MatLab: [L,U,P]=lu(A)

Una vez que tenemos una factorizaci´ on P A = LU , permutando las ecuaciones en el sistema Ax = b, obtenemos: P Ax = P b ! LU x = P b

Una vez que tenemos la factorizaci´ on P A = LU , permutando las ecuaciones en el sistema Ax = b, obtenemos: P Ax = P b ! LU x = P b Llamando eb = U x, se tiene que Leb = P b

Teorema. Si A es invertible y b ∈ Rn, entonces:

Existe P matriz de permutaci´ on tal que P A = LU, donde L es triangular inferior con unos en la diagonal y U es triangular superior

La eliminaci´ on gaussiana transforma el sistema Ax = b en otro equivalente U x = eb, donde eb es soluci´ on del sistema Leb = P b

Resolver Ax = b

!

(1) Obtener la factorizaci´ on P A = LU (2) Resolver el sistema triangular Leb = P b (3) Resolver el sistema triangular U x = eb

Costo de resolver el sistema Ax = b: n3 −n

Obtener la factorizaci´ on P A = LU :

3



=O

n3



flops

3

Resolver U x = eb (sustituci´ on regresiva):

n(n+1) flops 2

Resolver Leb = P b (sustituci´ on progresiva):

n(n−1) flops 2

Total:

n3 +3n2 −n 3



=O

n3 3



flops

La eliminaci´ on gaussiana es inestable num´ ericamente si los pivotes son muy peque˜ nos Supongamos que el primer pivote a11 es pr´ oximo a cero:

l21 = a21/a11 es muy grande si

a11 ≈ 0

(lo mismo para l31, . . . , ln1)

La eliminaci´ on gaussiana es inestable num´ ericamente si los pivotes son muy peque˜ nos

F2 − l21F1 ≈ −l21F1 (lo mismo para F3, . . . , Fn)

Las filas de A(2) son casi proporcionales

Para evitarlo, se utiliza la estrategia de pivoteo parcial (comando lu de Matlab)

Para evitarlo, se utiliza la estrategia de pivoteo parcial (comando lu de Matlab) ? Primera etapa (en las restantes etapas, se procede an´ alogamente)    (1) A = 

a11 a12 a21 a22 ... ... an1 an2



. . . a1n . . . a2n   =A ...  ...  . . . ann

? Primera etapa (en las restantes etapas, se procede an´ alogamente)    (1) A = 

a11 a12 a21 a22 ... ... an1 an2



. . . a1n . . . a2n   =A ...  ...  . . . ann

Se determina la posici´ on i tal que |ai1| = m´ ax {|a11| , . . . , |an1|} Se intercambian las filas 1 ↔ i

Situando el mayor elemento en m´ odulo de la fila en la posici´ on pivote, evitamos que los multiplicadores sean grandes: |l21| = |a21/a11| ≤ 1 (lo mismo para l31, . . . , ln1)

Normas matriciales

x ∈ Rn

norma eucl´ıdea kxk :=

q

2 x2 1 + · · · + xn

Normas matriciales

x ∈ Rn

norma eucl´ıdea kxk :=

q

2 x2 1 + · · · + xn

(

A m×n

kAk := m´ ax

kAxk : x ∈ Rn, x 6= 0 kxk

)

Propiedades de norma

• kAk ≥ 0

• kαAk = |α| kAk

• kA + Bk ≤ kAk + kBk

Propiedades de norma

• kAk ≥ 0

• kαAk = |α| kAk

• kA + Bk ≤ kAk + kBk

(⇒ podemos medir distancias entre matrices: dist(A, B) = kA − Bk)

C´ alculo de kAk

(

A m×n

kAxk : x ∈ Rn, x 6= 0 kAk := m´ ax kxk

)

C´ alculo de kAk

(

kAxk : x ∈ Rn, x 6= 0 kAk := m´ ax kxk

A m×n

Teorema. kAk =

q

ρ(AT A)

)

C´ alculo de kAk

(

kAxk : x ∈ Rn, x 6= 0 kAk := m´ ax kxk

A m×n

Teorema. kAk =

)

q

ρ(AT A)

M ρ(M ) = max λ se denomina radio espectral de la matriz M

(λM son los autovalores de M )

C´ alculo de kAk

(

kAxk : x ∈ Rn, x 6= 0 kAk := m´ ax kxk

A m×n

Teorema. kAk =

)

q

ρ(AT A)

M = AT A es sim´ etrica y semidef. posit. ⇒ autovalores reales no negativos

Otras propiedades

kOm×nk = 0, kInk = 1

kABk ≤ kAk kBk kQAk = kAk para toda matriz Q ortogonal (Q−1 = QT ). En particular: kQk = 1 para toda matriz Q ortogonal

Condicionamiento del producto matriz–vector Problema: Dado x, calcular y = Ax (A invertible)

Condicionamiento del producto matriz–vector Problema: Dado x, calcular y = Ax (A invertible) x∗ ≈ x

y ∗ = Ax∗ ≈ y

Condicionamiento del producto matriz–vector Problema: Dado x, calcular y = Ax (A invertible) x∗ ≈ x Error vectorial: e = x − x∗

y ∗ = Ax∗ ≈ y

Condicionamiento del producto matriz–vector Problema: Dado x, calcular y = Ax (A invertible) x∗ ≈ x Error absoluto: kek = kx − x∗k

y ∗ = Ax∗ ≈ y

Condicionamiento del producto matriz–vector Problema: Dado x, calcular y = Ax (A invertible) x∗ ≈ x kx − x∗k Error relativo: Rel(x) = kxk

y ∗ = Ax∗ ≈ y

Condicionamiento del producto matriz–vector Problema: Dado x, calcular y = Ax (A invertible) x∗ ≈ x

y ∗ = Ax∗ ≈ y

Rel(y) ≤ cond(A) · Rel(x) N´ umero de condici´ on: cond(A) = kA−1k kAk

Propiedades del n´ umero de condici´ on: cond(A) = kA−1k kAk

cond(A) ≥ 1 para toda matriz cond(Q) = 1 para toda matriz Q ortogonal (cond(I) = 1)

Condicionamiento de un sistema de ecuaciones lineales Problema: Dados A, b, calcular x tal que Ax = b

Condicionamiento de un sistema de ecuaciones lineales Problema: Dados A, b, calcular x tal que Ax = b  b∗ ≈ b

x∗ tal que Ax∗ = b∗

Condicionamiento de un sistema de ecuaciones lineales Problema: Dados A, b, calcular x tal que Ax = b  b∗ ≈ b

x∗ tal que Ax∗ = b∗ Rel(x) ≤ cond(A) · Rel(b)

Condicionamiento de un sistema de ecuaciones lineales Problema: Dados A, b, calcular x tal que Ax = b  A∗ ≈ A

x∗ tal que A∗x∗ = b

Condicionamiento de un sistema de ecuaciones lineales Problema: Dados A, b, calcular x tal que Ax = b  A∗ ≈ A

x∗ tal que A∗x∗ = b

Error matricial: E = A − A∗

Condicionamiento de un sistema de ecuaciones lineales Problema: Dados A, b, calcular x tal que Ax = b  A∗ ≈ A

x∗ tal que A∗x∗ = b Rel(x) . cond(A) · Rel(A)

CASOS PARTICULARES RELEVANTES

Matrices sim´ etricas y definidas positivas A es sim´ etrica si AT = A A definida positiva si xT Ax > 0 para todo x 6= 0

CASOS PARTICULARES RELEVANTES

Matrices sim´ etricas y definidas positivas A es sim´ etrica si AT = A A definida positiva si xT Ax > 0 para todo x 6= 0 2 (Q(x) = xT Ax = a11x2 1 + · · · + ann xn + 2a12 x1 x2 + 2a13 x1 x3 + · · ·

es la forma cuadr´ atica asociada a la matriz A)

A es sim´ etrica y definida positiva si, y s´ olo si:

Todos sus autovalores son positivos Todos las submatrices principales A(1 : k, 1 : k) tienen determinante positivo Todos los pivotes de la eliminaci´ on gaussiana son positivos

Para matrices sim´ etricas definidas positivas:

La eliminaci´ on gaussiana sin pivotaje es num´ ericamente estable La factorizaci´ on A = LU se traduce en la factorizaci´ on de Cholesky

Factorizaci´ on de Cholesky:

A = RT R,

con R triangular superior y diagonal positiva

PROBLEMAS SOBREDETERMINADOS: ´ DE M´INIMOS CUADRADOS. SOLUCION

Am × n, b ∈ Rm con m ≥ n

Ax = b sistema sobredeterminado

PROBLEMAS SOBREDETERMINADOS: ´ DE M´INIMOS CUADRADOS. SOLUCION

Am × n, b ∈ Rm con m ≥ n

Ax = b sistema sobredeterminado

x ˜ ∈ Rn soluci´ on en el sentido de los m´ınimos cuadrados si kb − A˜ xk ≤ kb − Axk ,

para todo x ∈ Rn

Geom´ etricamente: mejor aproximaci´ on de b al subespacio Col(A)

b

Col(A)

Ax

Geom´ etricamente: mejor aproximaci´ on de b al subespacio Col(A)

b

Col(A)

Ax

x ˜ es soluci´ on en m´ınimos cuadrados de Ax = b ⇔ A˜ x = ProyCol(A) (b)

Geom´ etricamente: mejor aproximaci´ on de b al subespacio Col(A)

b

Col(A)

Ax

x ˜ es soluci´ on en m´ınimos cuadrados de Ax = b ⇔ (b − A˜ x) ⊥ Col(A)

Teorema. Sean A m × n y b ∈ Rm. Entonces: x ˜ es soluci´ on en el sentido de los m´ınimos cuadrados de Ax = b si, y s´ olo si x ˜ es una soluci´ on del sistema compatible AT Ax = AT b (ecuaciones normales de Gauss)

Propiedades: Si x ˜ es soluci´ on de m´ınimos cuadrados, tambi´ en lo son x ˜ + x, para cualquier x ∈ Nul(A) (es decir, Ax = 0) Si el rango de A es m´ aximo ⇒ la soluci´ on de m´ınimos cuadrados es ´ unica.

Resolver Ax = b en el sentido de los m´ınimos cuadrados equivale a resolver las ecuaciones normales de Gauss (siempre compatibles), pero...

Resolver Ax = b en el sentido de los m´ınimos cuadrados equivale a resolver las ecuaciones normales de Gauss (siempre compatibles), pero... las ecuaciones normales de Gauss suelen tener mal condicionamiento

Descomposici´ on QR de una matriz A matriz cuadrada con determinante distinto de cero El m´ etodo de eliminaci´ on de Gauss aplicado a una matriz A devuelve una matriz U triangular superior factorizaci´ on LU

Descomposici´ on QR de una matriz A matriz m × n con columnas linealmente independientes El m´ etodo de ortogonalizaci´ on de Gram-Schmidt aplicado a las columnas de A devuelve una matriz Q con columnas ortonormales factorizaci´ on QR

Q = [q1| . . . |qn] tiene columnas ortonormales si: qiT qj = 0, para i 6= j (condici´ on de ortogonalidad) qiT qi = 1, para i = j (columnas unitarias)

Q = [q1| . . . |qn] tiene columnas ortonormales si: QT Q = I

Q = [q1| . . . |qn] tiene columnas ortonormales si: QT Q = I ! QT = Q−1

Q = [q1| . . . |qn] tiene columnas ortonormales si: QT Q = I ! QT = Q−1 ! Q es ortogonal

Teorema (Descomposici´ on QR). Sea Am×n, de rango n ≤ m. Entonces: A = QR, donde Qm×m es ortogonal (Q−1 = QT ), y Rm×n tiene rango n, sus n primeras filas forman una matriz cuadrada triangular superior y sus m − n ´ ultimas filas son nulas

Las operaciones elementales pueden cambiar las soluciones en m´ınimos cuadrados

Teorema (Descomposici´ on QR). Sea Am×n, de rango n ≤ m. Entonces: A = QR, donde Qm×m es ortogonal (Q−1 = QT ), y Rm×n tiene rango n, sus n primeras filas forman una matriz cuadrada triangular superior y sus m − n ´ ultimas filas son nulas Adem´ as, la soluci´ on de m´ınimos cuadrados de Ax = b coincide con la de Rx = QT b

(QRx = b)