TESIS DOCTORAL Alta Precisión Relativa en Problemas de Álgebra Lineal Numérica en Matrices con Estructura

Autor:

Johan Armando Ceballos Cañón Director:

Juan Manuel Molera Molera

DEPARTAMENTO DE MATEMÁTICAS

Leganés, Julio de 2013

TESIS DOCTORAL

Alta Precisión Relativa en Problemas de Álgebra Lineal Numérica en Matrices con Estructura

Autor:

Johan Armando Ceballos Cañón

Director:

Juan Manuel Molera Molera

Firma del Tribunal Calificador: Firma Presidente:

José Javier Martínez Fernández de las Heras

Vocal:

Esther Dopazo González

Secretario:

Fernando de Terán Vergara

Calificación:

Leganés,

de

de

A mi hija Mariana

Este trabajo se desarroll´o en el Departamento de Matem´aticas de la Universidad Carlos III de Madrid Campus de Legan´es bajo la direcci´on del profesor Juan Manuel Molera Molera. Se cont´o con una beca de la UC3M para el periodo del M´aster y posteriormente de un contrato de Personal Investigador de la Comunidad de Madrid. Adicionalmente se recibi´o ayuda parcial de los siguientes proyectos de investigaci´on: MINISTERIO DE EDU´ Y CIENCIA DIRECCION ´ GENERAL DE CACION ´ 2006/03964/001, MINISTERIO DE INVESTIGACION ´ 2010/00014/001, COMUCIENCIA E INNOVACION NIDAD DE MADRID - UC3M 2011/00120/001 y MINISTERIO DE ECONOM´IA Y COMPETITIVIDAD 2013/00065/001.

Agradecimientos El principal y m´as merecido agradecimiento se lo debo a mi director de tesis, Juan Manuel Molera, por haberme introducido y guiado en el mundo de la investigaci´on lo cual hizo posible esta tesis. Tambi´en quiero expresar mi agradecimiento a mis compa˜ neros del grupo de Algebra Lineal Num´erica, en especial a Froil´an M. Dopico por haberme permitido formar parte del mismo, de los proyectos de investigaci´on, por sus consejos y charlas que fueron vitales para el trabajo; sin duda aprend´ı mucho gracias a ellos. Debo extender este agradecimiento al Departamento de Matem´aticas por el apoyo financiero, a su planta de profesores y en especial al personal administrativo: Natalia Delgado y Alberto Calvo por su apoyo y buena disposici´on. A mi gran amigo Javier P´erez que siempre estuvo ah´ı para brindarme una palabra de aliento en el momento preciso. Por u ´ltimo, quisiera expresar mi m´as profundo agradecimiento a mi hija Mariana y a mis padres por su constante cari˜ no, confianza y apoyo, a ellos quiero dedicar esta tesis.

´Indice general

Agradecimientos

III

Notaci´ on

IX

Resumen

XI

1. Introducci´ on 2. RRD y HRA: Preliminares y resultados previos

1 17

2.1. Notaci´on . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

2.2. Teor´ıa de perturbaciones de EVP y SVD . . . . . . . . . . . . . . . .

18

2.3. Teor´ıa de perturbaciones de SEL y LSP . . . . . . . . . . . . . . . . .

20

2.4. Descomposiciones que revelan el rango y teor´ıa de perturbaciones . .

23

2.4.1. Descomposici´on en valores singulares y problema sim´etrico de autovalores . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24

2.4.2. Sistemas de ecuaciones lineales v´ıa la RRD . . . . . . . . . . .

26

2.5. Algoritmos y errores . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

3. Teor´ıa de perturbaciones multiplicativas y soluciones precisas del LSP 35 3.1. Preliminares y conceptos b´asicos . . . . . . . . . . . . . . . . . . . . .

38 v

3.2. Perturbaciones multiplicativas para la pseudoinversa de Moore-Penrose 40 3.2.1. Cotas de perturbaciones multiplicativas para la pseudoinversa de Moore-Penrose . . . . . . . . . . . . . . . . . . . . . . . . .

44

3.3. Perturbaciones multiplicativas para problemas de m´ınimos cuadrados

48

3.3.1. ¿Por qu´e el factor kA† k2 kbk2 /kx0 k2 es peque˜ no? . . . . . . . .

52

3.3.2. El n´ umero de condici´on para perturbaciones multiplicativas del problema de m´ınimos cuadrados . . . . . . . . . . . . . . .

54

3.3.3. Cotas de perturbaciones multiplicativas para otras soluciones del problema de m´ınimos cuadrados . . . . . . . . . . . . . . .

55

3.4. Perturbaci´on del problema de m´ınimos cuadrados en los factores . . .

56

3.5. Algoritmo y an´alisis de errores . . . . . . . . . . . . . . . . . . . . . .

57

3.6. Experimentos Num´ericos . . . . . . . . . . . . . . . . . . . . . . . . .

63

3.6.1. Matrices de Cauchy . . . . . . . . . . . . . . . . . . . . . . . .

64

3.6.2. Matrices de Vandermonde . . . . . . . . . . . . . . . . . . . .

66

3.6.3. Matrices Graduadas . . . . . . . . . . . . . . . . . . . . . . .

69

3.6.4. Experimentos num´ericos controlando el residuo . . . . . . . .

75

3.6.5. Experimentos donde kA† k2 kbk2 / kx0 k2 no es peque˜ no . . . . .

76

4. C´ alculo de autovalores y autovectores precisos de matrices sim´ etricas graduadas 79 4.1. Matrices sim´etricas indefinidas y m´etodo de pivotaje diagonal . . . .

83

4.2. Teor´ıa de perturbaciones de matrices sim´etricas graduadas . . . . . .

86

4.3. Algoritmos y an´alisis de errores . . . . . . . . . . . . . . . . . . . . .

91

4.4. Experimentos y ejemplos num´erico . . . . . . . . . . . . . . . . . . .

99

4.4.1. Experimento num´erico . . . . . . . . . . . . . . . . . . . . . . 100 4.4.2. Ejemplo num´erico . . . . . . . . . . . . . . . . . . . . . . . . . 100 5. Conclusiones y trabajos futuros

105

Bibliograf´ıa

109

vi

´Indice de figuras

3.6.1.Error relativo progresivo kb x0 − x0 k2 /kx0 k2 frente κ2 (C). C matrices de Cauchy aleatorias de orden 100 × 50. . . . . . . . . . . . . . . . .

67

3.6.2.Error relativo progresivo kb x0 − x0 k2 /kx0 k2 frente n, para matrices de Cauchy de orden m × n con m = 50 y n = 10 : 2 : 40 . . . . . . . . . .

68

3.6.3.Error relativo progresivo kb x0 − x0 k2 /kx0 k2 frente κ2 (V ) para matrices de Vandermonde aleatorias V de tama˜ nos 100 × 10 : 5 : 60. . . . . . . .

70

3.6.4.Error relativo progresivo kb x0 − x0 k2 /kx0 k2 frente κ2 (B) para matrices graduadas aleatorias A = S1 BS2 de tama˜ no 100 × 40. . . . . . . . . .

73

4.4.1.Error relativo progresivo en los autovalores. Φ como en la ecuaci´on 4.4.3 frente τ para matrices sim´etricas graduadas aleatorias A = SBS de orden 100. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 4.4.2.Error relatico progresivo en los autovectores. Ψ como en la ecuaci´on 4.4.4 frente τ para matrices sim´etricas graduadas aleatorias A = SBS de orden 100. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

vii

Tabla de S´ımbolos Para comodidad del lector presentamos las siguientes tablas de notaciones, que se seguir´an a lo largo de toda la tesis.

S´ımbolo R, C Rn , C n Rm×n , Cm×n In AT rango(A) R(A) A−1 A† A≤B |A| θ(x, y) PX λi bi λ σi σ bi k · k2 k · kF

Conjunto de los n´ umeros reales, complejos Conjunto de vectores de n componentes en R, C Conjunto de matrices de m filas y n columnas con componentes en R, C Matriz identidad de orden n Traspuesta de A Rango de A Espacio columna de A Inversa de A Pseudoinversa de Moore-Penrose de A Desigualdades matriciales componente a componente: aij ≤ bij ∀i, j Valor absoluto de A: |A|ij = |aij | ∀i, j ´ Angulo agudo entre los vectores x e y Proyecci´on ortogonal sobre el subespacio X i-´esimo autovalor de A i-´esimo autovalor calculado de A i-´esimo valor singular de A i-´esimo valor singular calculado de A Norma vectorial 2 o Eucl´ıdea: kxk2 = (xT x)1/2 1/2 P P m n 2 |a | Norma matricial de Frobenius: kAkF = ij i=1 j=1

k · k2 O(·) u κ(A) relgapλ (A, λi )

Norma matricial 2 o norma espectral: kAk2 = σm´ax (A) Notaci´on O-grande de Landau Unidad de redondeo N´ umero de condici´on de A: κ(A) = kAkkA−1 k gap relativo en los autovalores: relgapλ (A, λi ) = m´ınj6=i

|λj −λi | |λi |

relgapσ (A, σi ) gap relativo en los valores singulares: relgapσ (A, σi ) = m´ınj6=i

|σj −σi | σi

ix

Tabla de Siglas

Siglas ALN HRA RRD EVP SEVP SVD LSE LSP GE GECP GEPP BP-LDLT flops

x

´ Algebra Lineal Num´erica Alta Precisi´on Relativa (High Relative Accuracy) Descomposici´on que revela el rango (Rank revealing decomposition) Problema de autovalores (Eigenvalue problem) Problema de autovalores de matrices sim´etricas (Symmetric eigenvalue problem) Descomposici´on en valores singulares (Singular value decomposition) Sistemas de ecuaciones lineales (Linear systems equations) Problema de m´ınimos cuadrados (Least squares problems) Eliminaci´on Gaussiana (Gaussian elimination) Eliminaci´on Gaussiana est´andar con pivote completo (GE complete pivoting) Eliminaci´on Gaussiana est´andar con pivote parcial (GE partial pivoting) Factorizaci´on sim´etrica por bloques de Bunch y Parlett con pivote completo Operaciones en coma flotante (Floating point operations)

Resumen

Esta tesis se enmarca dentro del campo de la Alta Precisi´ on Relativa (HRA) ´ en Algebra Lineal Num´erica (ALN). Sus l´ıneas maestras son dos. Por un lado, el ´ dise˜ no y an´alisis de algoritmos que permitan resolver problemas de Algebra Lineal con m´as precisi´on de la habitual para matrices con estructura. Y por otro el estudio de la teor´ıa espec´ıfica de perturbaciones necesaria para tratar los problemas que nos ocupan. En nuestra investigaci´on hemos tratado dos: La obtenci´on de soluciones precisas del problema de m´ınimos cuadrados para matrices con estructura (Cap´ıtulo 3). La obtenci´on de autovalores y autovectores precisos para matrices sim´etricas graduadas (Cap´ıtulo 4). ´ Cl´asicamente, el Algebra Lineal Num´erica (ALN) es una disciplina de los M´etodos Num´ericos que desarrolla algoritmos eficientes y estables para Resolver sistemas de ecuaciones lineales (LSE): Ax = b con A ∈ Cn×n y b ∈ Cn , Resolver problemas de m´ınimos cuadrados (LSP): m´ınx∈Cn kAx − bk2 , donde k · k2 es la norma eucl´ıdea usual, A ∈ Cm×n y b ∈ Cm , Calcular autovalores y autovectores de matrices (EVP): Ax = λx, A ∈ Cn×n , x ∈ Cn y λ ∈ C, y Calcular la Descomposici´on en Valores Singulares (SVD) de matrices: A = U ΣV ∗ con A ∈ Cm×n , U y V matrices unitarias y Σ ∈ Cm×n diagonal con entradas no negativas. xi

Resumen ´ Los cuatro problemas cl´asicos del Algebra Lineal Num´erica (ALN) siguen siendo objeto de intensa investigaci´on debido a cuatro causas principales relacionadas entre s´ı: (a) las aplicaciones actuales requieren c´alculos eficientes con matrices cada vez m´as grandes para las que los algoritmos existentes son muy lentos; (b) las arquitecturas de los ordenadores evolucionan continuamente, lo que obliga a crear o modificar algoritmos para buscar la m´axima eficiencia; (c) los requerimientos de precisi´on cada vez son mayores; (d) aparecen continuamente nuevas clases de matrices estructuradas para las cuales se debe intentar desarrollar algoritmos espec´ıficos m´as eficientes y/o precisos que los existentes. Hoy en d´ıa el ALN es una disciplina muy amplia que engloba muchos otros problemas motivados por las aplicaciones. Los problemas cl´asicos y modernos del ALN est´an completamente permeados de la idea clave de aprovechar la estructura concreta de las clases de matrices que aparecen en las aplicaciones en el desarrollo de algoritmos espec´ıficos para dichas matrices, que resuelvan num´ericamente los problemas con m´as rapidez y/o precisi´on que los algoritmos v´alidos para matrices generales y que reduzcan la memoria utilizada por el ordenador. Los mejores algoritmos actualmente existentes en el ALN son estables en sentido regresivo, pero esto, a veces, no es suficiente porque se pueden cometer errores inadmisiblemente grandes cuando se aplican a matrices con n´ umeros de condici´on muy grandes. El objetivo general de la investigaci´on en algoritmos de alta precisi´on dentro del ALN es aprovechar la estructura de ciertas clases de matrices para calcular magnitudes con un error mucho menor que el de los algoritmos tradicionales v´alidos para matrices generales y esencialmente con el mismo coste computacional que ellos, esto es, coste O(n3 ) operaciones en coma flotante (flops) en matrices n × n. Alcanzar alta precisi´on (HRA) s´olo es posible para tipos particulares de matrices mediante algoritmos espec´ıficos para las mismas que explotan la estructura del problema para acelerar la velocidad de los c´alculos, disminuir los requerimientos de memoria y mejorar la precisi´on de la soluci´on en comparaci´on con algoritmos est´andar. Nuestra investigaci´on se enmarca en un programa general para encontrar solu´ ciones precisas de los cuatro problemas cl´ asicos del Algebra Lineal Num´erica, aprovechando las propiedades de las Descomposiciones que Revelan el Rango (RRD) xii

Resumen usando un esquema b´asico de algoritmo de 2 pasos que puede ser usado para resolver diferentes problemas con HRA: Algoritmo (Algoritmo de 2 pasos para obtener HRA) Paso 1: Calcular una RRD precisa de A = XDY.

Paso 2: Aplicar algoritmos, espec´ıficos a cada problema tratado, a los factores X, D, Y para obtener la respuesta.

La ventaja de este esquema es su modularidad y que es adaptable a distintos problemas y tipos de matrices. En principio permite obtener HRA para cualquier tipo de matriz para la que se pueda obtener una RRD precisa. Este esquema de Algoritmo en 2 pasos, partiendo de una RRD precisa, se aplica ya a la SVD, a Sistemas de Ecuaciones Lineales (LSE), al Problema Sim´etrico de Autovalores (SEVP). Nosotros lo hemos ampliado al Problema de M´ınimos Cuadrados (LSP), Cap´ıtulo 3, y al Problema Sim´etrico de Autovalores (SEVP) para matrices graduadas, Cap´ıtulo 4. En el Cap´ıtulo 3 se ha desarrollado una nueva teor´ıa de perturbaciones multiplicativa para la pseudoinversa de Moore-Penrose (Secci´on 3.2) y para la soluci´on de m´ınima longitud del problema de M´ınimos Cuadrados (Secci´on 3.3), y se ha implementado un algoritmo para el LSP y se ha realizado su an´alisis de errores (Secci´on 3.5), y los correspondientes experimentos num´ericos (Secci´on 3.6). El resultado que dan los m´etodos cl´asicos para el error de la soluci´on de m´ınima longitud de un problema de M´ınimos Cuadrados es: kb x0 − x0 k2 ≤ u g(m, n) κ22 (A), kx0 k2

(0.0.1)

donde g(m, n) es una funci´on moderadamente creciente de m y n. Esta cota no garantiza ning´ un d´ıgito de precisi´on en la soluci´on calculada si la matriz A est´a mal condicionada. Desafortunadamente, muchos tipos de matrices estructuradas que aparecen en las aplicaciones est´an mal condicionadas y los algoritmos convencionales calcular´an soluciones del problema de m´ınimos cuadrados con errores relativos progresivos grandes. La nueva teor´ıa de perturbaciones multiplicativas, unida al an´alisis de errores del algoritmo que proponemos, muestran que el error progresivo cometido en la soluci´on de m´ınima longitud (v´alido para problemas de rango completo, rango xiii

Resumen deficiente, o incluso para sistemas lineales infraterminados m < n) viene dado por:   kb x0 − x0 k2 kA† k2 kbk2 ≤ c u py (m, n) κ2 (Y ) + px (m, n) κ2 (X) + O(u2 ) , (0.0.2) kx0 k2 kx0 k2 donde c es una constante entera peque˜ na, py (m, n) y px (m, n) son funciones de m y n moderadamente creciente. El u ´nico factor potencialmente grande es kA† k2 kbk2 /kx0 k2 , pero un an´alisis cuidadoso de este factor muestra que es peque˜ no para la mayor´ıa de vectores b. En el Cap´ıtulo 4 se muestra el trabajo realizado para encontrar soluciones precisas del Problema Sim´etrico de Autovalores (SEVP) para matrices graduadas. Una matriz sim´etrica A = AT ∈ Rn×n es llamada graduada si S −1 AS −1 ≡ B es una matriz bien condicionada para alguna matriz diagonal escalada S = diag[s1 , . . . , sn ]. En primer lugar se ha demostrado una nueva teor´ıa de perturbaciones estructurada para matrices de la forma A = SBS. En la Secci´on 4.2 se presenta el Teorema 4.2.3 e = S(B + δB)S. que nos da la sensibilidad del problema a perturbaciones de tipo A Se usa la t´ecnica de transformar perturbaciones aditivas en perturbaciones multiplicativas. Se encuentra que la sensibilidad del problema viene gobernada, entre otros, por un factor nuevo τD . La sensibilidad bajo perturbaciones de ese tipo depende del n´ umero de condici´on de los correspondientes factores LDLT de B y de los elementos de la matriz diagonal S de dos maneras: su orden, despu´es del pivotaje, y el tama˜ no relativo de los elementos consecutivos en las posiciones de los bloques 2 × 2 de la matriz D. Este efecto es completamente nuevo y no se ha tenido en cuenta en previos an´alisis. El algoritmo que se ha usado en este caso tambi´en ha constado de 2 pasos. Con el fin de calcular una RRD se ha usado la factorizaci´on sim´etrica por bloques P AP T = LDLT con estrategia de pivote completo de Bunch y Parlett (factorizaci´on BP-LDLT ). Para el segundo paso se ha usado el Algoritmo de Jacobi impl´ıcito [30]. El an´alisis de errores del algoritmo junto con la teor´ıa de perturbaciones multiplicativas muestra que los autovalores y los autovectores se calculan con errores que, a primer orden en la unidad de redondeo u, son   b − λi | |λ i b ≤ q(n) u τ ΞB + κ2 (L) + O(u2 ) |λi |    2 b sin θ(qi , qbi ) ≤ q(n) u τ ΞB + κ2 (L) 1 + + O(u2 ) relgapλb (A, λi ) b es el factor L de la donde ΞB es una cantidad peque˜ na si A est´a bien escalada, L factorizaci´on BP-LDLT de A y θ(qi , qbi ) es el a´ngulo agudo entre los autovectores xiv

Resumen exacto y calculado, respectivamente, qi y qbi . El factor τ controla el escalamiento y es definido como el m´aximo de tres factores,

τ := max{1, τL , τD },

sk τL := m´ax j η 0 . En el caso de valores singulares m´ ultiples, o grupos de valores singulares muy pr´ oximos, hay una cota similar para los senos de los ´angulos can´ onicos entre los subespacios correspondientes. Hemos modificado ligeramente el enunciado que se presenta en [20, Teorema 2.1] a˜ nadiendo el punto 1 de los resultados. Queremos destacar de forma expl´ıcita que bajo las hip´otesis del teorema, la perturbaci´on de los factores se puede escribir como una peque˜ na perturbaci´on multiplicativa de la matriz. Los otros dos resultados del teorema se obtienen inmediatamente sin m´as que aplicar el teorema 2.2.4. De esta manera destacamos el papel crucial de la teor´ıa multiplicativa de perturbaciones en establecer, en este caso, que valores y vectores singulares son poco sensibles a perturbaciones de la RRD. De una forma similar al problema de valores singulares anterior se puede tratar el problema sim´etrico de autovalores. Dopico y Koev demostraron en [28] tambi´en que a partir de una RRD sim´etrica y precisa de una matriz tambi´en sim´etrica es posible calcular sus autovalores y autovectores con alta precisi´on relativa. 25

Cap´ıtulo 2 e=X eD eX e T RRDs de las Teorema 2.4.4 [28, Teorema 2.1] Sean A = XDX T y A e ∈ Rn×n respectivamente. Sean λ1 ≥ · · · ≥ λn los matrices sim´etricas A ∈ Rn×n y A e1 ≥ · · · ≥ λ en los autovalores de A. e Sean q1 , . . . , qn y qe1 , . . . , qen autovalores de A y λ los correspondientes autovectores ortonormales. Supongamos que e − Xk2 kX ≤ β, kXk2 e ii − Dii | |D ≤ β para todo i, |Dii | donde 0 ≤ β < 1. Sea η = β(2 + β)κ2 (X) menor que 1; entonces ei | ≤ (2η + η 2 )|λi |, |λi − λ

1 ≤ i ≤ n,

 1+

2+η relgapλe (A, λi )

y η sin θ(qi , qei ) ≤ 1−η

 ,

1 ≤ i ≤ n,

donde θ(qi , qei ) es el ´angulo agudo entre qi y qei . Para el caso de autovalores m´ ultiples, o cl´ usters de autovalores muy cercanos en el sentido relativo, una cota similar se cumple para los senos de los ´angulos can´ onicos de los correspondientes subespacios invariantes.

2.4.2.

Sistemas de ecuaciones lineales v´ıa la RRD

Las t´ecnicas anteriores tambi´en se han extendido a los otros dos grandes pro´ blemas del Algebra Lineal Num´erica, los sistemas de ecuaciones lineales [32] y los problemas de m´ınimos cuadrados [11]. Este u ´ltimo problema es el objeto del Cap´ıtulo 3 de esta memoria, y all´ı nos ocuparemos de ´el. En [32] se demostr´o que si la matriz de un LSE Ax = b se puede factorizar de forma precisa entonces la soluci´on del LSE cambia poco al variar los factores de la RRD. Se presenta aqu´ı s´olo el resultado a primer orden (para el resultado completo ver [32, Teorema 3.2]) Teorema 2.4.5 [32, Teorema 3.2] Sea A ∈ Cn×n con RRD XDY y b ∈ Cn . Sea XDY x = b y (X + ∆X)(D + ∆D)(Y + ∆Y )e x = b + h, donde k∆Xk ≤ kXk, k∆Y k ≤ kY k |∆D| ≤ |D| y khk ≤ kbk y supongamos que κ(Y ) < 1 y que (2 + )κ(X) < 1. Entonces, a primer orden en ,     ke x − xk kX −1 bk kA−1 k kbk ≤  κ(Y ) + 1 + 2 kXk + O(2 ). (2.4.4) kxk kbk kxk 26

RRD y HRA: Preliminares y resultados previos O, en una forma m´as d´ebil, pero de mayor utilidad en la pr´ actica:   ke x − xk kA−1 k kbk ≤  κ(Y ) + (1 + 2 κ(X)) + O(2 ), kxk kxk ya que no requieren que se calcule kX −1 bk y no sobreestima significativamente la variaci´on de la soluci´on si κ(X) es peque˜ no. −1

Ya se present´o en la Secci´on 2.3 una breve explicaci´on de por qu´e el factor kA kxkk kbk es habitualmente peque˜ no. En la Secci´on 3.3.1 se presentar´a el factor equivalente para kA† k kbk el LSP kxk y se ver´a tambi´en que es usualmente peque˜ no.

2.5.

Algoritmos y errores

El estudio de la alta precisi´on relativa tiene su fin u ´ltimo en poder calcular magnitudes con m´as precisi´on que la que garantizan los algoritmos convencionales. En esta secci´on vamos a describir los algoritmos presentes en la literatura que garantizan HRA. Describiremos los algoritmos y los resultados que se obtienen de sus an´alisis de errores, es decir la precisi´on esperada cuando se ejecutan en un ordenador. Antes de empezar recordaremos que en los an´alisis de errores presentados en la Secciones 3.5 y 4.3 usaremos el modelo de errores convencional para aritm´etica en coma flotante [43, Secci´on 2.2] f l(a b) = (a b)(1 + δ), donde a y b son n´ umeros reales en coma flotante, ∈ {+, −, ×, /}, y |δ| ≤ u. Recordemos que este modelo tambi´en se cumple para n´ umeros complejos en coma flotante si u es reemplazada por una constante ligeramente m´as grande, ver [43, Secci´on 3.6]. Suponemos tambi´en que no ocurre overflow ni underflow. El campo de la HRA tiene una larga tradici´on (ve´anse las referencias citadas en el Cap´ıtulo 1), sin embargo en la forma en que se enfoca en esta memoria data de los primeros a˜ nos de la d´ecada de 1990. Entre las referencias destacamos la de Demmel y Veseli´c [25]. En ella se inicia el esquema, que ha resultado tan fruct´ıfero despu´es, de Factorizaci´on+Algoritmo espec´ıfico para lograr la HRA. En ese trabajo se mostraba que el algoritmo de Jacobi impl´ıcito aplicado a los factores de Cholesky de una matriz sim´etrica definida positiva pod´ıa alcanzar HRA, al contrario que el algoritmo QR; ´ esto tambi´en sirvi´o para traer a primera plana de los algoritmos del Algebra Lineal Num´erica los algoritmo de tipo Jacobi que estaban ligeramente olvidados. 27

Cap´ıtulo 2 Los algoritmos de tipo Jacobi est´an basados en las rotaciones del mismo nombre. i 

R(i, j, c, s) =

    i      j    

j 

1 ..

      ,      

. −s

c ... s

c ..

.

(2.5.1)

1 donde los cosenos c, y los senos s, se eligen en cada caso para que la caja 2 × 2 adecuada se diagonalice: " #T " #" # " # c −s aii aij c −s µ1 0 = . s c aji ajj s c 0 µ2 Para m´as detalles ver los siguientes libros cl´asicos [17, Secci´on 5.3.5], [39, Secci´on 8.5] y [61, Cap´ıtulo 9]. Quiz´as la forma m´as conocida de usar las rotaciones de Jacobi es para diagonalizar una matriz sim´etrica haciendo rotaciones simult´aneas, por la derecha y por la izquierda. Es lo que se conoce como el algoritmo de Jacobi “two-sided”. Ligeramente menos conocido es el algoritmo de Jacobi “one-sided”que sirve para calcular la SVD de una matriz (ver por ejemplo [39]). La esencia del algoritmo “one-sided” es aplicar impl´ıcitamente el algoritmo “two´ sided” a la matriz G = AT A. Esta no se calcula nunca expl´ıcitamente, seg´ un se van necesitando sus elementos se calculan con la expresi´on:

gij =

n X

aki akj

(2.5.2)

k=1

Algoritmo 2.5.1 (Algoritmo de Jacobi “one-sided”) Input: A ∈ Rn×n Output: Σ = diag(σ1 , . . . , σn ) y V ∈ Rn×n matrices de valores y vectores, singulares derechos respectivamente. V =I 28

RRD y HRA: Preliminares y resultados previos repetir para i, j calcular gii , gij , gjj de G = AT A como en la ecuaci´ on (2.5.2) calcular una rotaci´ on de Jacobi R(i, j, c, s) tal que: " # " # g g µ 0 ii ij 1 RT (i, j, c, s) R(i, j, c, s) = gij gjj 0 µ2 A = A R(i, j, c, s) V = V R(i, j, c, s) fin para |g | hasta convergencia √ ij

|gii gjj |

calcular

< tol

σi = kG(:, i)k para i = 1, . . . , n.

los vectores singulares izquierdos, ser´ an las columnas normalizadas de G. A partir del algoritmo “one-sided” 2.5.1 para calcular la SVD, se puede dise˜ nar el siguiente algoritmo [25] para calcular los autovalores y los autovectores de matrices sim´etricas definidas positivas: Algoritmo 2.5.2 [25, Algoritmo 4.4] Input:

A = AT ∈ Rn×n definida positiva.

Output: Autovalores λi y autovectores Q de A. Paso 1: Calcular una descomposici´on de Cholesky con pivote completo de A = L LT . Paso 2: Calcular la SVD de L usando el Algoritmo 2.5.1: LR1 R2 · · · Rp = U Σ, donde Σ = diag(σ1 , . . . , σn ). Paso 3: λi = σi2 para i = 1 : n y Q = U .

Obs´ervese que es un algoritmo de tipo 2 pasos como los usados en esta memoria y en otras situaciones (ver el Cap´ıtulo 1). Se demostr´o en [25, 55, 56] que si los autovalores y los autovectores de matrices sim´etricas definidas positivas son calculados por el algoritmo 2.5.2, entonces, los errores en los autovalores calculados est´an dados por el siguiente resultado. 29

Cap´ıtulo 2 Teorema 2.5.3 Sea A = SBS una matriz sim´etrica definida positiva, con S = 1/2 1/2 bi }n los autovadiag[a11 , . . . , ann ], bii = 1 y {λi }ni=1 los autovalores de A. Sean {λ i=1 lores calculados por el Algoritmo 2.5.2 en aritm´etica finita con precisi´ on u. Entonces # " bi − λi | |λ 1 1 ≤ O(u) +p |λi | λn (B) λn (B)

(2.5.3)

Si la matriz A est´a bien escalada, es decir, si κ2 (B) & 1, entonces todos los autovalores est´an calculados con varios d´ıgitos de precisi´on. El siguiente hito en la historia reciente de los algoritmos de alta precisi´on relativa, es el Algoritmo 3.1 de [20], que calcula la soluci´on precisa del problema de valores singulares de A, v´ıa la RRD de A. Algoritmo 2.5.4 (Soluci´on precisa del problema de valores singulares) Input: X ∈ Rm×n , D ∈ Rn×n e Y ∈ Rn×n tales que A = XDY es una RRD de A. Output: Σ = diag(σ1 , . . . , σn ) la matriz de valores singulares y, U ∈ Rm×n y V ∈ Rn×n matrices de vectores singulares izquierdos y derechos respectivamente. Paso 1: Calcular la factorizaci´on QR con pivote por columnas XD = QRP . Paso 2: Multiplicar las matrices R, P e Y para obtener W = RP Y Paso 3: Calcular la SVD de W = U ΣV T usando el Algoritmo 2.5.1. Paso 4: multiplicar las matrices Q y U para obtener U = QU La alta precisi´on relativa la garantizamos a partir de la teor´ıa multiplicativa de perturbaciones, Teorema 2.4.3, y el an´alisis de errores del Algoritmo: Teorema 2.5.5 [33, Teorema 2.1] El Algoritmo 2.5.4 produce errores multiplicativos regresivos cuando es ejecutado con unidad de redondeo u, es decir, Si A = XDY ∈ bΣ b Vb T es la SVD Rm×n es la RRD calculada en el paso 1 del Algoritmo 2.5.4 y U calculada por el algoritmo, entonces, existen matrices U 0 ∈ Rm×r , V 0 ∈ Rn×r , E ∈ Rm×m y F ∈ Rn×n tales que U 0 y V 0 tienen columnas ortonormales,

30

b k = O(u), kU 0 − U

kV 0 − Vb k = O(u)

kEk = O(uκ(X)),

kF k = O(uκ(R0 )κ(Y )),

RRD y HRA: Preliminares y resultados previos donde R0 es la matriz con mejor n´ umero de condici´ on de todos los escalamientos por filas de la matriz triangular R que aparece en el Algorithm 3.1 de [20] y b 0T . (I + E)A(I + F ) = U 0 ΣV Se demostr´o en [20], que κ(R0 ) es a lo sumo de orden O(n3/2 κ(X)), pero en la pr´actica se observa mediante diferentes experimentos num´ericos que κ(R0 ) se comporta como O(n). A continuaci´on presentamos otro hito significativo. El algoritmo de Jacobi impl´ıcito [30, Algoritmo 1] aplicado a matrices sim´etricas, definidas o indefinidas, produce HRA si se parte de una matriz de la que se ha calculado una RRD precisa. Sea A = XΩX T una factorizaci´on sim´etrica de A, donde X ∈ Rn×n y Ω = diag(ω1 , . . . , ωn ) ∈ Rn×n son matrices no singulares. La idea es aplicar el algoritmo sim´etrico (“two-sided”) de Jacobi de forma impl´ıcita, eligiendo en cada paso una rotaci´on de Jacobi R tal que la posici´on (i, j) de A = XΩX T sea cero. En cada paso del proceso s´olo se actualiza el factor X y la matriz Ω, que es la que porta el mal condicionamiento, permanece constante, es decir, pasamos de XΩX T a RT (XΩX T )R, manteniendo la matriz factorizada (X → RT X). Se demostr´o en [30, Teorema 6] que si XΩX T es una descomposici´on precisa que revele el rango, entonces el Algoritmo de Jacobi impl´ıcito, calcula soluciones precisas del problema de autovalores y autovectores. Algoritmo 2.5.6 (Algoritmo de Jacobi impl´ıcito [30, Algoritmo 1]) Input: X ∈ Rn×n una matriz bien condicionada y Ω = diag(ω1 , . . . , ωn ) ∈ Rn×n no singular. Output: Λ = diag(λ1 , . . . , λn ) y U ∈ Rn×n matrices de autovalores y autovectores, respectivamente. κ b(X) es el valor calculado de κ(X) U = In p p G = X diag( |ω1 |, . . . , |ωn |) J = diag(sign(ω1 ), . . . , sign(ωn )) repetir para i = 1 hasta n − 1 para j = i + 1 hasta n calcular aii , aij , ajj de A = GJGT usando P aij = nk=1 gik gjk sign(ωk ) 31

Cap´ıtulo 2 calcular " " # # " # a a c −s µ 0 ii ij 1 T = , c2 + s2 = 1 tal que T T T = aij ajj s c 0 µ2 T G = R(i, j, c, s) G U = U R(i, j, c, s)2 fin para fin para 

calcular λi = aii

√ |aij |

≤  m´ax{n, κ b(X)} para todo i < j y  Pn 2 k=1 gik ≤ 2b κ(X) para todo i |aii | para i = 1, . . . , n.

hasta convergencia

|aii ajj |

Las cotas de errores multiplicativas regresivas dadas en el Teorema 2.5.7 pueden combinarse f´acilmente con los Teoremas 2.2.2 y 2.2.3 para probar que el Algoritmo 2.5.6 calcula los autovalores y los autovectores de una RRD XDX T con alta precisi´on relativa: Teorema 2.5.7 [30, Teorema 6] Si NR es el n´ umero de rotaciones de Jacobi aplicadas en el Algoritmo 2.5.6 hasta que se satisfaga el criterio de parada, κ b(X)γn+1 < 18 , √ b1 , . . . , λ bn ), b = diag(λ y nκ(X)γ2 < 12 , entonces la matriz de autovalores calculada, Λ b , est´ y la matriz calculada de autovectores, U an cerca de las matrices de autovalores y autovectores de una peque˜ na perturbaci´on multiplicativa de XDX T . Es decir, existe una matriz ortogonal U ∈ Rn×n tal que b T = (I + E) XDX T (I + E)T , U ΛU

(2.5.4)

b − U kF = O(NR u). con kEkF = O(u (n2 κ b(X) + NR κ(X))) y kU Por u ´ltimo, extendiendo las ideas de usar la RRD como buena parametrizaci´on que refleja la sensibilidad de los problemas espectrales a los sistemas de ecuaciones lineales, Dopico y Molera [32] mostraron c´omo se pueden calcular soluciones con alta precisi´on de LSE, Ax = b, para cualquier matriz de la que se pueda obtener una RRD precisa mediante el siguiente algoritmo: Algoritmo 2.5.8 [32, Algoritmo 4.1] Input: A ∈ Cn×n , b ∈ Cn 2

Ver (2.5.1), R(i, j, c, s) es la notaci´ on usual para las rotaciones de Jacobi en los ´ındices (i, j) y en el ´angulo θ, es decir: R(i, j, c, s) = In excepto en las posiciones R(i, j, c, s)ii = R(i, j, c, s)jj = cos(θ) = c y R(i, j, c, s)ij = −R(i, j, c, s)ji = − sin(θ) = −s.

32

RRD y HRA: Preliminares y resultados previos Output: x , soluci´on de A x = b Paso 1: Calcular una RRD precisa de A = XD Y , con D = diag(d1 , d2 , . . . , dn ). Paso 2: Resolver los tres sistemas de ecuaciones, Xs=b −→ s, D w = s −→ w, Y x = w −→ x, donde X s = b e Y x = w se resuelven por cualquier m´etodo regresivo estable, tal como, eliminaci´ on Gaussiana con pivote parcial (GEPP) o factorizaci´on QR, y D w = s es resuelto como wi = si /di , i = 1 : n. El Teorema 2.4.5 expresa la sensitividad de la soluci´on del sistema Ax = b, donde ´ A = XDY , para perturbaciones de los factores de A. Esta depende de los n´ umeros de condici´on de X e Y , los cuales son peque˜ nos si A = XDY es una RRD, y tambi´en −1 de la expresi´on kA kkbk/kxk, la cual es moderada para la mayor´ıa de vectores b, seg´ un los comentarios realizados al final de la Secci´on 2.4. El an´alisis de errores del Algoritmo 2.5.8 demuestra que ´estos se pueden escribir como peque˜ nos errores regresivos multiplicativos, y as´ı se obtiene que: Teorema 2.5.9 Sea k · k una norma vectorial en Cn cuya norma matricial inducida satisface que kΛk = m´axi=1:n |λi | para todas las matrices diagonales Λ = b D, b e Yb los factores de A calculados en el Paso 1 del diag(λ1 , . . . , λn ). Sean X, Algoritmo 2.5.8 y supongamos que satisfacen b − Xk kX ≤ u p(n), kXk

kYb − Y k ≤ u p(n), kY k

y

b − D | ≤ u p(n) |D |, |D

i = 1 : n,

donde X, D e Y son los correspondientes factores exactos de A, p(n) es una funci´on poco creciente de n, y u es la unidad de redondeo. Supongamos tambi´en que los sistemas X s = b e Y x = w se resolvieron con un algoritmo regresivamente estable tal que cuando es aplicado a cualquier sistema lineal Bz = c, B ∈ Cn×n y c ∈ Cn , calcula una soluci´on zb que satisface (B + ∆B )b z = c,

con k∆B k ≤ u q(n) kB k,

√ donde q(n) es una funci´on poco creciente n tal que q(n) ≥ 4 2/(1 − 4u). Sea g(n) := p(n) + q(n) + u p(n)q(n). 33

Cap´ıtulo 2 1. Si x b es la soluci´on calculada de A x = b usando el Algoritmo 2.5.8, entonces (X + ∆X )(D + ∆D )(Y + ∆Y )b x = b, donde k∆X k ≤ u g(n) kX k, |∆D | ≤ u g(n) |D |, y k∆Y k ≤ u g(n) kY k. 2. Adem´ as, si x es la soluci´on exacta de Ax = b, ( u g(n) κ(Y ) ) < 1 y ( u g(n) (2 + u g(n)) κ(X) ) < 1, entonces u g(n) kb x − xk ≤ kxk 1 − u g(n)κ(Y )

 κ(Y ) +

 1 + (2 + u g(n)) κ(X) kA−1 k kbk . 1 − u g(n)(2 + u g(n))κ(X) kxk (2.5.5)

Observaci´ on 2.5.10 La ecuaci´on (2.5.5) se puede simplificar considerablemente si s´olo prestamos atenci´on a los t´erminos de primer orden  
kb x − xk kA−1 k kbk ≤ u g(n) κ(Y ) + (1 + 2 κ(X)) + O (u g(n))2 , kxk kxk −1
kA k kbk ≤ 4u g(n) m´ax{κ(X), κ(Y )} + O (u g(n))2 . kxk Estas cotas muestran que el error relativo en la soluci´ on es O(u) si κ(X), κ(Y ) y −1 kA k kbk/kxk son peque˜ nos. En los Cap´ıtulos 1 y 2 hemos presentado el estado actual del campo de la alta precisi´on relativa. Dentro de las ideas mostradas aqu´ı se ha encuadrado la investigaci´on que presentamos en esta memoria. La idea clave es la obtenci´on de HRA mediante algoritmos que usan RRDs precisas y que luego aprovechan la estructura espec´ıfica del problema. Este esquema ya se ha aplicado con ´exito a problemas de valores singulares, a problemas sim´etricos de autovalores y a sistemas de ecuaciones lineales. Nosotros lo hemos extendido a problemas de m´ınimos cuadrados (Cap´ıtulo 3) y al problema sim´etrico de autovalores para un tipo especial de matrices: graduadas sim´etricas. Los trabajos existentes en el campo al comienzo de mi investigaci´on han dado la pauta, pero en este memoria se presentan nuevos algoritmos, Algoritmo 3.5.1 y 4.3.1; nuevos an´alisis de errores, Teorema 3.5.2 y Corolario 4.3.5; y nuevos resultados de teor´ıa de perturbaciones de la pseudo inversa de Moore-Penrose, Teoremas 3.2.2, 3.2.3, 3.2.5, nuevos resultados de teor´ıa de perturbaciones de problemas de m´ınimos cuadrados 3.3.1 y nuevos resultados de teor´ıa de perturbaciones de problemas de autovalores de matrices sim´etricas graduadas Corolario 4.2.4.

34

3 Teor´ıa de perturbaciones multiplicativas y soluciones precisas del problema de m´ınimos cuadrados

Distintas clases de matrices con estructuras particulares aparecen frecuentemente en la teor´ıa y las aplicaciones [58, 59]. Como consecuencia, el dise˜ no y an´alisis de algoritmos especiales que involucran matrices estructuradas es un ´area cl´asica de ´ investigaci´on del Algebra Lineal Num´erica que atrae la atenci´on de muchos investigadores. Algoritmos especiales para resolver sistemas de ecuaciones estructurados o problemas estructurados de autovalores, aparecen en muchas referencias est´andar [17, 39, 43, 49, 73], pero es m´as raro encontrar algoritmos especiales para resolver problemas de m´ınimos cuadrados. En general, el objetivo, al considerar algoritmos especiales es explotar la estructura del problema para acelerar la velocidad de los c´alculos y disminuir los requerimientos de memoria y mejorar la precisi´on de la soluci´on en comparaci´on con algoritmos est´andar. Sobre este u ´ltimo objetivo, mencionemos que ´ desde los primeros d´ıas del Algebra Lineal Num´erica se han desarrollado algoritmos especiales para resolver sistemas de ecuaciones lineales estructurados de manera m´as precisa que con algoritmos est´andar (ver las referencias en [32, 43]). El desarrollo de algoritmos precisos para problemas de autovalores es mucho m´as reciente. Sus inicios datan del principio de los 90’s con el famoso art´ıculo [21], el cual recibi´o considerable atenci´on (ver por ejemplo, [3, 20, 25, 30, 33, 35, 36, 38, 48, 66, 75] entre otras muchas ´ referencias). El presente cap´ıtulo trata de una parte del “Algebra Lineal Num´erica Precisa” para la cual no existen muchas referencias disponibles en la literatura: algoritmos para resolver problemas estructurados de m´ınimos cuadrados m´ınx kb − Axk2 , donde A ∈ Rm×n y b ∈ Rm , con m´as precisi´on que la dada por algoritmos est´andar y m´as o menos el mismo coste computacional, es decir, O(n2 m) flops. Sobre este 35

Cap´ıtulo 3 tema s´olamente conocemos la referencia [53], la cual trata de una clase particular de problemas de m´ınimos cuadrados. El m´etodo est´andar para resolver problemas de m´ınimos cuadrados m´ınx kb − Axk2 con rango completo por columnas es utilizar la factorizaci´on QR calculada con el algoritmo de Householder [43, Cap´ıtulo 19 y 20]. Dicho m´etodo es regresivamente estable, es decir, la soluci´on calculada x b0 es la soluci´on exacta del problema de m´ınimos cuadrados m´ınx k(b + ∆b) − (A + ∆A)xk2 , donde k∆bk2 ≤ c u mn kbk2 , k∆Ak2 ≤ c u m n3/2 kAk2 , u es la unidad de redondeo del ordenador y c denota una peque˜ na constante entera [43, Teorema 20.3]. Resultados de errores regresivos an´alogos se cumplen para otros m´etodos de soluci´on del problema de m´ınimos cuadrados basado en descomposiciones ortogonales. Por ejemplo, la descomposici´on en valores singulares (SVD)1 . Este resultado fuerte de errores regresivos, junto con la teor´ıa de perturbaciones cl´asica “normwise” del problema de m´ınimos cuadrados [74, Teorema 5.1] (ver tambi´en [5, Teorema 1.4.6, p. 30]), implica la siguiente cota de error progresiva en la soluci´on calculada x b0 con respecto a la soluci´on exacta x0 del problema de m´ınimos cuadrados   kA† k2 kbk2 kb x0 − x0 k 2 3/2 2 kb − Ax0 k2 ≤ (c u m n ) κ2 (A) + + κ2 (A) , (3.0.1) kx0 k2 kx0 k2 kAk2 kx0 k2 donde A† es la pseudoinversa de Moore-Penrose de A, kAk2 denota la norma espectral de A y κ2 (A) = kAk2 kA† k2 es el n´ umero de condici´on espectral de A. La cota en (3.0.1) es mayor que u κ2 (A) y puede ser mucho mayor en ciertas condiciones, y as´ı, (3.0.1) no garantiza ning´ un d´ıgito de precisi´on en la soluci´on calculada si κ2 (A) & 1/u, es decir, si A est´a mal condicionada con respecto a la inversa de la unidad de redondeo. Desafortunadamente, muchos tipos de matrices estructuradas que aparecen en las aplicaciones est´an extremadamente mal condicionadas y los algoritmos est´andar para problemas de m´ınimos cuadrados pueden calcular soluciones con errores relativos muy grandes. Dos ejemplos famosos son las matrices de Vandermonde, las cuales aparecen en problemas de interpolaci´on polinomial y las matrices de Cauchy [43, Cap´ıtulos 22 y 28]. Nuestro objetivo en este trabajo es presentar un marco te´orico general para la soluci´on del problema de m´ınimos cuadrados y probar rigurosamente que dicho m´etodo permite calcular, para muchas clases de matrices estructuradas soluciones, con cotas de error mucho m´as peque˜ nas que la presentada en (3.0.1). El marco te´orico del que 1

Cabe se˜ nalar que el error regresivo en A cometido por la soluci´on del problema de m´ınimos cuadrados v´ıa la factorizaci´ on QR de Householder es “columnwise”, es decir, k∆A(:, j)k2 ≤ c u mn kA(: , j)k2 para j = 1 : n (notaci´ on de MATLAB) y, por lo tanto, este resultado es m´as fuerte que el mencionado anteriormente. Sin embargo, esta cota “columnwise” no se cumple para la soluci´on v´ıa SVD, ya que las transformaciones ortogonales son aplicadas a la matriz A por ambos lados.

36

Teor´ıa de perturbaciones multiplicativas y soluciones precisas del LSP hablamos se basa en el concepto de descomposici´ on que revela el rango (RRD), inicialmente introducida en [20] para calcular la SVD con alta precisi´on relativa. nuestro m´etodo consiste en calcular la soluci´on m´ınima de m´ınx kb − Axk2 en 2 pasos: 1. Primer paso. Calcular una RRD precisa de A = XDY , en el sentido de [20] (revisamos el significado de RRD precisa en la Definici´on 2.4.2). 2. Segundo paso. Este paso lo realizaremos en tres etapas: (1) calcular la soluci´on u ´nica s de m´ınx kb − Xxk2 v´ıa la factorizaci´on QR de Householder; (2) calcular la soluci´on w del sistema lineal Dw = s como wi = si /di , i = 1 : r; y (3) calcular la soluci´on m´ınima x0 del sistema lineal infraterminado Y x = w utilizando el m´etodo Q [43, Cap´ıtulo 21]. El vector x0 es la soluci´on de m´ınima longitud de m´ınx kb − Axk2 . Veremos que el procedimiento descrito anteriormente calcula soluciones precisas incluso para matrices A extremadamente mal condicionadas, ya que el sistema lineal mal condicionado Dw = s es resuelto de manera precisa, y adem´as que m´ınx kb − Xxk2 e Y x = w tambi´en son resueltos de manera precisa ya que X e Y est´an bien condicionadas. Probaremos en la Secci´on 3.5 que el error relativo para la soluci´on de m´ınima longitud x b0 de m´ınx kb − Axk2 calculada por el m´etodo propuesto es   kA† k2 kbk2 kb x0 − x0 k2 ≤ u f (m, n) κ2 (Y ) + κ2 (X) , (3.0.2) kx0 k2 kx0 k2 donde f (m, n) es una funci´on moderadamente creciente de m y n. N´otese primero que (3.0.2) es mejor que (3.0.1), ya que X e Y est´an bien condicionadas, y as´ı, el u ´nico factor potencialmente grande en (3.0.2) es kA† k2 kbk2 /kx0 k2 , el cual tambi´en aparece en (3.0.1). Pero el verdadero punto importante en la cota (3.0.2) es que si A es fija, entonces kA† k2 kbk2 /kx0 k2 es peque˜ no para la mayor´ıa de los vectores b, incluso para matrices A muy mal condicionadas. Este hecho es bien conocido si A es cuadrada y no singular, como lo indicamos en la Secci´on 2.3 (ver tambi´en [2, 14] y [32, Secci´on 3.2]) y, como explicaremos en la Secci´on 3.3.1, tambi´en se cumple para matrices generales en dos sentidos: para la mayor´ıa de vectores aleatorios b, y para la mayor´ıa de vectores b con un valor fijo del residuo relativo kAx0 − bk2 /kbk2 no cercano a 1. En este cap´ıtulo la frase “para la mayor´ıa de los vectores b” se puede entender en cualquiera de esos dos sentidos. La idea y los resultados discutidos anteriormente se parecen a los presentados en [32] para calcular soluciones precisas de sistemas de ecuaciones lineales estructurados 37

Cap´ıtulo 3 Ax = b con A no singular. Sin embargo, el an´alisis para problemas de m´ınimos cuadrados es mucho m´as complicado y requiere t´ecnicas muy diferentes para desarrollar una nueva teor´ıa de perturbaciones multiplicativas, necesaria para probar la cota de error presentada en (3.0.2). Adem´as, los resultados y algoritmos que presentamos son generales: son v´alidos para matrices A con y sin rango completo. Aunque nos centramos principalmente en problemas de m´ınimos cuadrados, estos resultados pueden ser aplicados para resolver sistemas de ecuaciones lineales infradeterminados. La mayor´ıa de los algoritmos citados en la introducci´on del Cap´ıtulo 2 determinan exactamente el rango de matrices de rango deficiente y para estas clases de matrices listadas anteriormente, el m´etodo introducido en este art´ıculo resuelve problemas de m´ınimos cuadrados con alta precisi´on relativa acotados como en (3.0.2). Esta cota de error es O(u f (m, n)) para la mayor´ıa de los vectores b independientemente del n´ umero de condici´on tradicional de las matrices, y as´ı garantiza soluciones precisas. Este cap´ıtulo est´a organizado como sigue. En la Secci´on 3.1 se introduce la notaci´on b´asica, los conceptos y los resultados que se utilizar´an en todo el cap´ıtulo. La Secci´on 3.2 estudia la variaci´on de la pseudoinversa de Moore-Penrose bajo perturbaciones multiplicativas y, basado en estos resultados, en la Secci´on 3.3 presentamos cotas para perturbaciones multiplicativas de problemas de m´ınimos cuadrados. Como consecuencia, en la Secci´on 3.4 presentamos cotas perturbativas para problemas de m´ınimos cuadrados cuya matriz de coeficientes est´a dada como una RRD con factores perturbados. En la Secci´on 3.5 presentamos un nuevo algoritmo para calcular de manera precisa el problema de m´ınimos cuadrados v´ıa la RRD y el correspondiente an´alisis de errores regresivo y progresivo. La precisi´on de este algoritmo se verifica en la Secci´on 3.6 por medio de exhaustivos experimentos num´ericos.

3.1.

Preliminares y conceptos b´ asicos

Dado que consideramos problemas de m´ınimos cuadrados, usaremos la norma m´as natural para estos problemas: la norma vectorial eucl´ıdea, es decir, dado 1/2 x = [x1 , . . . , xn ]T ∈ Rn , kxk2 := (|x1 |2 + · · · + |xn |2 ) , y para matrices A ∈ Rm×n la correspondiente norma matricial inducida kAk2 := m´axkxk2 =1 kAxk2 , llamada la norma espectral o norma 2 de A. En la Secci´on 3.2.1, usaremos tambi´en normas matriciales unitariamente invariantes [69, Cap´ıtulo II. Secci´on 3], que ser´an denotadas por k · k. El s´ımbolo In denotar´a la matriz identidad de orden n × n, pero usaremos simplemente I si el tama˜ no es claro en el contexto, y AT denota la transpuesta de A. Usaremos la notaci´on de MATLAB para submatrices: A(i : j, :) indica la subma38

Teor´ıa de perturbaciones multiplicativas y soluciones precisas del LSP triz de A consistente de las filas i hasta la j y A(:, k : l) indica la submatriz de A consistente de las columnas k hasta la l. Dada A ∈ Rm×n , con m ≥ n, sus valores singulares se denotar´an como σ1 (A) ≥ · · · ≥ σn (A) ≥ 0. El Lema 3.1.1 ser´a utilizado para derivar algunas cotas de perturbaciones. Lema 3.1.1 Sean B, C ∈ Rm×n , S ⊆ Rm y W ⊆ Rn subespacios vectoriales y sean PS ∈ Rm×m y PW ∈ Rn×n proyectores ortogonales sobre S y W, respectivamente. Entonces se cumple: q (a) kPS B + (I − PS )Ck2 ≤ kBk22 + kCk22 . (b) kBPW + C(I − PW )k2 ≤

q

kBk22 + kCk22 .

Demostraci´on. (a). Sea x ∈ Rn tal que kxk2 = 1. Dado que los vectores PS Bx y (I − PS )Cx son ortogonales, entonces k(PS B + (I − PS )C)xk22 = kPS Bxk22 + k(I − PS )Cxk22 ≤ kBxk22 + kCxk22 ≤ kBk22 + kCk22 y q kPS B + (I − PS )Ck2 = m´ax k(PS B + (I − PS )C)xk2 ≤ kBk22 + kCk22 . kxk2 =1

La parte (b) se sigue de aplicar (a) a la transpuesta y de que para cualquier matriz, kAk2 = kAT k2 . Dada una matriz G ∈ Rm×n con entradas gij , denotaremos por |G| a la matriz con entradas |gij |. Expresiones como |G| ≤ |B|, donde B ∈ Rm×n , significan que |gij | ≤ |bij | para 1 ≤ i ≤ m, 1 ≤ j ≤ n. La pseudoinversa de Moore-Penrose de A ∈ Rm×n juega un papel importante en este trabajo. Es bien sabido que la u ´nica matriz Z ∈ Rn×m tal que (i) AZA = A,

(ii) ZAZ = Z,

(iii) (AZ)T = AZ,

y

(iv) (ZA)T = ZA, (3.1.1)

o equivalentemente, tal que AZ = PA

y

ZA = PZ ,

(3.1.2)

donde PA y PZ son los proyectores ortogonales sobre el espacio columna de A y de Z, respectivamente. La equivalencia de las cuatro condiciones en (3.1.1) y las dos condiciones en (3.1.2) pueden ser encontradas en [10, Teorema 1.1.1]. Denotaremos por A† ∈ Rm×n la pseudoinversa de Moore-Penrose de A ∈ Rm×n . Recordemos 39

Cap´ıtulo 3 que si A ∈ Rn×n es no singular, entonces A† = A−1 y tambi´en que la SVD de A nos permite obtener una expresi´on para A† y probar muchas de sus propiedades [69, Cap´ıtulo 3]. R(A) y N (A) denotar´an respectivamente el espacio columna y el espacio nulo de A. Es f´acil ver que R(AT ) = R(A† ), entonces, por (3.1.2), PA = AA† y PAT = PA† = A† A son respectivamente los proyectores ortogonales sobre R(A) y R(AT ). Enunciamos en el Lema 3.1.2 algunas propiedades muy conocidas [10] de la pseudoinversa de Moore-Penrose que necesitaremos a lo largo de este cap´ıtulo.

Lema 3.1.2 (a) Si A tiene rango completo por filas, entonces A† = AT (AAT )−1 y AA† = I. (b) Si A tiene rango completo por columnas, entonces A† = (AT A)−1 AT y A† A = I. (c) Sean F ∈ Rm×r y G ∈ Rr×n . Si rango (F ) = rango (G) = r, entonces (F G)† = G† F † . La soluci´on de m´ınima longitud del problema de m´ınimos cuadrados m´ınx∈Rn kb− Axk2 es x0 = A† b y la soluci´on de m´ınima longitud de un sistema lineal indeterminado Ax = b es tambi´en x0 = A† b. Si A = XDY ∈ Rm×n es una RRD de A, entonces dos aplicaciones del Lema 3.1.2-(c) implican que A† = Y † D−1 X † y la soluci´on de m´ınima longitud del problema de m´ınimos cuadrados o de un sistema lineal indeterminado es x0 = Y † D−1 X † b.

3.2.

Perturbaciones multiplicativas para la pseudoinversa de Moore-Penrose

En esta secci´on y en la Secci´on 3.3, consideramos una perturbaci´on multiplicativa e = (I + E)A(I + F ), donde de una matriz general A ∈ Rm×n , es decir, una matriz A (I + E) ∈ Rm×m e (I + F ) ∈ Rn×n son matrices no singulares. El objetivo final es acotar en la Secci´on 3.3, la expresi´on ke x0 − x0 k2 /kx0 k2 , donde x0 y x e0 son las soluciones de m´ınima longitud de los problemas de m´ınimos cuadrados m´ınx∈Rn kAx − bk2 e − ebk2 , respectivamente. Este objetivo se alcanza a trav´es del teorema y m´ınx∈Rn kAx principal de esta secci´on, el Teorema 3.2.2, donde obtenemos dos expresiones para e† en terminos de A† , (I + E)−1 e (I + F )−1 . Utilizamos las expresiones de A e† para A 40

Teor´ıa de perturbaciones multiplicativas y soluciones precisas del LSP e† −A† k/kA† k en cualquier norma unitariadesarollar en la Secci´on 3.2.1 cotas para kA mente invariante y en particular, en la norma 2. Aunque estas cotas no son necesarias para nuestro resultado final, resaltamos que encontrar cotas perturbativas para la pseudoinversa de Moore-Penrose es un tema cl´asico en el An´alisis Matricial (ver [74] y [69, Cap´ıtulo 3, Secci´on 3]) que ha llamado la atenci´on de muchos investigadores. Mostramos que los resultados en la Secci´on 3.2.1 son mejores que los presentados en [9], los cuales tienen una naturaleza diferente y son obtenidos a trav´es de un m´etodo diferente. La teor´ıa multiplicativa de perturbaciones de matrices es muy importante en el c´alculo de autovalores y valores singulares con alta precisi´on [37, 45, 46, 51, 52] y tambi´en el c´alculo de soluciones precisas de sistemas de ecuaciones [32, Lema 3.1]. Hasta donde sabemos, no ha sido estudiada en el contexto de encontrar soluciones precisas del problema de m´ınimos cuadrados. El Lema 3.2.1 es un resultado t´ecnico que utilizaremos en la demostraci´on del Teorema 3.2.2. e = (I +E)A(I +F ) ∈ Rm×n , donde (I +E) ∈ Rm×m Lema 3.2.1 Sean A ∈ Rm×n y A e (I + F ) ∈ Rn×n son matrices no singulares. Entonces las siguientes igualdades se cumplen: (a) PA (I + E T )(I − PAe) = 0. (b) (I − PAeT )(I + F T )PAT = 0. e = R((I +E)A) entonces (I −P e)(I +E)A = 0. Demostraci´on. (a) Dado que R(A) A As´ı que, (I − PAe)(I + E)AA† = (I − PAe)(I + E)PA = 0, lo cual es equivalente a PA (I + E T )(I − PAe) = 0. eT = (I + F T )AT (I + E T ), luego conjugar y transponer (b) Aplicar la parte (a) a A la igualdad. Ahora, enunciamos el resultado principal de esta secci´on, el cual es v´alido para matrices de rango completo, matrices de rango deficiente y para perturbaciones de cualquier magnitud. e = (I + E)A(I + F ) ∈ Rm×n , donde (I + E) ∈ Teorema 3.2.2 Sean A ∈ Rm×n y A Rm×m e (I + F ) ∈ Rn×n son matrices no singulares. Entonces e† = P eT (I + F )−1 A† (I + E)−1 P e A A A

(3.2.1) 41

Cap´ıtulo 3 y     e† = I + (I − P eT )F T − P eT Fb A† I + E T (I − P e) − EP b e , A A A A A

(3.2.2)

b = (I + E)−1 E y Fb = (I + F )−1 F . donde E Demostraci´on. Definamos Z := PAeT (I + F )−1 A† (I + E)−1 PAe. Probaremos que Z e Recordemos que P eT = A e† A e satisface las condiciones (3.1.2) reemplazando A por A. A eA e† . Entonces y PAe = A e = A(I e + F )−1 A† (I + E)−1 P e = (I + E)AA† (I + E)−1 P e AZ A A † −1 e e† † † e eA e† = P e. = (I + E)AA (I + E) AA = (I + E)AA A(I + F )A = A A An´alogamente, e = P eT (I + F )−1 A† (I + E)−1 A e = P eT (I + F )−1 A† A(I + F ) ZA A A † e −1 † † e e e† A e = P eT . = A A(I + F ) A A(I + F ) = A (I + E)AA† A(I + F ) = A A

(3.2.3) eT ) ⊆ R(Z) y la definici´on de Z implica que La igualdad (3.2.3) implica que R(A eT ). As´ı, R(Z) = R(A eT ) y la ecuaci´on (3.2.3) implica Z A e = PZ . Por lo R(Z) ⊆ R(A eyZ=A e† , con lo que se demuestra tanto, las condiciones (3.1.2) se cumplen para A (3.2.1). A continuaci´on, utilizamos (3.2.1) para demsotrar (3.2.2). Primero, escribimos b e (I + F )−1 = I − (I + F )−1 F = I − Fb. (I + E)−1 = I − (I + E)−1 E = I − E Sustituyendo estas expresiones en (3.2.1), tenemos e† = P eT (I − Fb)A† (I − E)P b e = P eT (PAT − Fb)A† (PA − E)P b e. A A A A A

(3.2.4)

Por el Lema 3.2.1-(a) se sigue que PA (I + E T (I − PAe)) = PA PAe. An´alogamente, del Lema 3.2.1-(b) tenemos que, ((I − PAeT )F T + I)PAT = PAeT PAT . Finalmente, sustituyendo estas relaciones en (3.2.4), y usando que A† PA = A† y PAT A† = A† , obtenemos (3.2.2). Si m = n y A es no singular, entonces PAe = PAeT = In , (3.2.1) y (3.2.2) se e−1 = (I +F )−1 A−1 (I +E)−1 . Si A tiene rango completo por columnas, convierten en A e tiene rango completo por columnas, P eT = In , (3.2.1) se simplifica a A e† = entonces A  A    e† = I − Fb A† I + E T (I − P e) − EP b e . (I + F )−1 A† (I + E)−1 P e y (3.2.2) a A A

A

A

e tambi´en tiene rango Finalmente, si A tiene rango completo por filas, entonces A 42

Teor´ıa de perturbaciones multiplicativas y soluciones precisas del LSP −1 † −1 e† completo por filas,  PAe = Im , (3.2.1) se simplifica  a A = PAeT (I + F ) A (I + E) e† = I + (I − P eT )F T − P eT Fb A† I − E b . y (3.2.2) a A A

A

Detacamos que la expresi´on (3.2.2) garantiza, que bajo peque˜ nas perturbaciones multiplicativas de A, obtenemos peque˜ nas perturbaciones multiplicativas de A† . e Esto Las hip´otesis del Teorema 3.2.2 garantizan que rango (A) = rango (A). simplifica considerablemente el an´alisis del cambio de la pseudoinversa de Mooree = A + ∆A [69, 74]. Adem´as, el Penrose con respecto a perturbaciones aditivas A Teorema 3.2.3 implica que si la condici´on d´ebil m´ax{kEk2 , kF k2 } < 1 se cumple, e = (I +E)A(I +F ) es una perturbaci´ entonces A on aguda de A (ver la definici´on en [74, Definici´on 7.2] y tambi´en en [69, Cap´ıtulo III, Definici´on 3.2]). Es bien conocido que perturbaciones agudas introducen simplificaciones a´ un para perturbaciones aditivas e A = A + ∆A. Las cotas en el Teorema 3.2.3 ser´an utilizadas en la Secci´on 3.3. e = (I + E)A(I + F ) ∈ Rm×n , donde (I + E) ∈ Teorema 3.2.3 Sean A ∈ Rm×n y A Rm×m e (I + F ) ∈ Rn×n son matrices no singulares, adem´ as consideremos PN (A) y e PN (A) e los proyectores ortogonales sobre el espacio nulo de A y A, respectivamente. Entonces: (a) kPAe − PA k2 = kPAe(I − PA )k2 = kPA (I − PAe)k2 ≤ kEk2 . (b) kPAeT − PAT k2 = kPAeT (I − PAT )k2 = kPAT (I − PAeT )k2 ≤ kF k2 . (c) kPN (A) e − PN (A) k2 ≤ kF k2 . e tienen la misma dimensi´on. Demostraci´on. (a). Los subespacios R(A) y R(A) As´ı que, de [69, Cap´ıtulo I, Teorema 5.5], kPAe − PA k2 = kPAe(I − PA )k2 = kPA (I − PAe)k2 . Adem´as, por el Lema 3.2.1-(a), kPA (I − PAe)k2 = k − PA E T (I − PAe)k2 ≤ kE T k2 = kEk2 . eT = (I +F T )AT (I +E T ). La parte (b) se demuestra aplicando (a) a la expresi´on A Finalmente, la parte (c) se obtiene de (b), PN (A) = I − PAT y PN (A) e = I − PA eT . e† − A† que se obtiene direcEl Corolario 3.2.4 proporciona una expresi´on para A tamente de (3.2.2). Este corolario ser´a utilizado en la Secci´on 3.2.1 y en el Teorema 3.3.1. e = (I + E)A(I + F ) ∈ Rm×n , donde (I + Corolario 3.2.4 Sean A ∈ Rm×n y A b = (I + E)−1 E y E) ∈ Rm×m e (I + F ) ∈ Rn×n son matrices no singulares, E Fb = (I + F )−1 F . Entonces e† − A† = A† ΘE + ΘF A† + ΘF A† ΘE , A

(3.2.5) 43

Cap´ıtulo 3 donde b e ΘE = E T (I − PAe) − EP A

3.2.1.

y

ΘF = (I − PAeT )F T − PAeT Fb .

(3.2.6)

Cotas de perturbaciones multiplicativas para la pseudoinversa de Moore-Penrose

Como explicamos en el primer p´arrafo de la Secci´on 3.2, los resultados en esta secci´on no ser´an utilizados en ning´ un otro lugar del cap´ıtulo. El objetivo principal en e† − A† k/kA† k. Supongamos que A 6= 0, ya que esta secci´on es presentar cotas para kA en otro caso el problema es trivial. Como resultado principal tenemos el Teorema 3.2.5.

e = (I + E)A(I + F ) ∈ Rm×n , donde (I + E) ∈ Teorema 3.2.5 Sean A ∈ Rm×n y A b = (I + E)−1 E y Fb = Rm×m e (I + F ) ∈ Rn×n son matrices no singulares, y sean E (I + F )−1 F . Denotemos por k · k una norma invariante unitariamente normalizada y por k · k2 la norma espectral. Entonces:

(a)

   e† − A† k kA b bk+ kEk + kEk b bk . ≤ kEk+k Ek+kF k+k F kF k + k F e† k2 } m´ın{kA† k2 , kA

e† − A† k2 kA ≤ (b) kA† k2 e† − A† k2 kA ≤ e† k2 kA

r

 2 b 2 + kFbk2 + kEk b 2 kFbk2 kEk22 + kF k22 + kEk

y

q b 22 + kFbk22 + (kEk2 + kF k2 + kEk2 kF k2 )2 . kEk

e† − A† k/kA† k2 se obtiene directamente de Demostraci´on. (a). La cota para kA (3.2.5), teniendo s´olamente en cuenta que para cualquier matriz B y C, kBCk ≤ e† − A† k/kA e† k2 se kBk2 kCk y kBCk ≤ kBkkCk2 [69, p´agina 80]. La cota para kA e† − A† k/kA† k2 intercambiando los roles de A y A, e es decir, obtiene de la cota para kA e o bien A = (I + E)−1 A(I e + considerando A como una perturbaci´on multiplicativa A, b A(I e − Fb), esto equivale a intercambiar kEk y kF k en las cotas por F )−1 = (I − E) b y kFbk, respectivamente, y vice versa. kEk 44

Teor´ıa de perturbaciones multiplicativas y soluciones precisas del LSP (b). Observemos que de (3.2.6), el Lema 3.2.1-(b) y PAT = A† A, tenemos (I + ΘF )A† A = (I + (I − PAeT )F T − PAeT Fb)PAT = PAT + (I − PAeT )F T PAT − PAeT FbPAT = PAT − (I − P eT )PAT − P eT FbPAT A

A

= PAeT (I − Fb)PAT , y si multiplicamos por la derecha la expresi´on anterior por A† (I + ΘF )A† = PAeT (I − Fb)A† .

(3.2.7)

Similarmente, pero usando el Lema 3.2.1-(a), tenemos b e. A† (I + ΘE ) = A† (I − E)P A

(3.2.8)

Ahora, usamos (3.2.5), (3.2.7)-(3.2.8) y (3.2.6) para demostrar que e† − A† = (I + ΘF )A† ΘE + ΘF A† A = PAeT (I − Fb)A† ΘE + ΘF A† h i = PAeT (I − Fb)A† ΘE − FbA† + (I − PAeT )F T A† h i † † b = PAeT A ΘE − F A (I + ΘE ) + (I − PAeT )F T A† i h b e + (I − P eT )F T A† . b e − FbA† (I − E)P = PAeT A† E T (I − PAe) − A† EP A A A (3.2.9) Queda entonces aplicar el Lema 3.1.1 a la ecuaci´on (3.2.9) y obtenemos e† − A† k22 ≤ kA† E T (I − P e) − [A† E b + FbA† (I − E)] b P ek22 + kF T A† k22 kA A A † T 2 †b † 2 b b ≤ kA E k2 + kA E + F A (I − E)k2 + kF k22 kA† k22   b 2 + kFbk2 + kEk b 2 kFbk2 )2 , ≤ kA† k22 kEk22 + kF k22 + (kEk e† − A† k2 /kA† k2 en la parte (b). Con lo que, lo cual nos proporciona una cota para kA e† − A† k2 /kA e† k2 se obtiene intercambiando los roles de A y A e como la cota para kA lo hicimos en la demostraci´on de la parte (a). Observaci´ on 3.2.6 Resaltamos los siguientes puntos del Teorema 3.2.5. (a) Las cotas en el Teorema 3.2.5 mejoran significativamente las cotas cl´asicas para el cambio relativo de la pseudoinversa de Moore-Penrose a perturbaciones e = A + ∆A (ver [74, Teorema 4.1] y [69, Cap´ıtulo III, Corolario aditivas A 3.10]). El punto importante es que las cotas en el Teorema 3.2.5 no dependen de κ2 (A), mientras que las cotas cl´asicas si. 45

Cap´ıtulo 3 (b) La cota en la parte (a) del Teorema 3.2.5 tiene la ventaja de ser v´alida para cualquier norma invariante unitariamente normalizada, pero para la norma 2, la cota en la parte (b) es siempre mejor que la dada en la parte (a), dado que p x2 + y 2 ≤ x + y, para n´ umeros reales x e y tales que x ≥ 0, y ≥ 0 y r  2 b 2 + kFbk2 + kEk b 2 kFbk2 ≤ kEk22 + kF k22 + kEk

≤

q

b 2 + kFbk2 + kEk b 2 kFbk2 kEk22 + kF k22 + kEk

b 2 + kF k2 + kFbk2 + (kEk2 + kEk b 2 ) (kF k2 + kFbk2 ) . ≤ kEk2 + kEk e tiene rango completo por filas (c) Si A tiene rango completo por filas, entonces A y adem´as PAe = Im . As´ı que, la expresi´on para ΘE que aparece en (3.2.6) se b y todos los t´erminos que contienen kEk o´ kEk2 en las simplifica a ΘE = −E b y kEk b 2 ). cotas del Teorema 3.2.5 se anulan (pero es necesario mantener kEk e tiene rango completo por (d) Si A tiene rango completo por columnas, entonces A columnas y adem´as PAeT = In . As´ı que, la expresi´on para ΘF que aparece en (3.2.6) se simplifica a ΘF = −Fb y todos los t´erminos que contienen kF k o´ kF k2 en las cotas del Teorema 3.2.5 se anulan (pero es necesario mantener kFbk y kFbk2 ). (e) Si en el Teorema 3.2.5 restringimos el tama˜ no de las perturbaciones a m´ax{kEk2 , kF k2 } < 1, una condici´on que garantiza que (I + E) e (I + F ) sean no singulares, entonces las desigualdades est´andar de normas matriciales [39] implican b 2≤ kEk

kEk2 1 − kEk2

y kFbk2 ≤

kF k2 . 1 − kF k2

(3.2.10)

Estas desigualdades pueden ser usadas en la parte (b) del Teorema 3.2.5 para obtener cotas en t´erminos de kEk2 y kF k2 . (g) Finalmente, con la restricci´on adicional m´ax{kEk2 , kF k2 } < 1, el Teorema 3.2.5 puede completarse con kA† k2 kA† k2 e† k2 ≤ ≤ kA . (1 + kEk2 )(1 + kF k2 ) (1 − kEk2 )(1 − kF k2 ) El lado derecho de la desigualdad proviene de (3.2.1), lo cual implica que e† k2 ≤ k(I + F )−1 k2 kA† k2 k(I + E)−1 k2 . Para la desigualdad del lado izkA e es decir, quierdo: consideremos A como una perturbaci´on multiplicativa de A, 46

Teor´ıa de perturbaciones multiplicativas y soluciones precisas del LSP e + F )−1 , y apliquemos (3.2.1) con los roles de A y A e inA = (I + E)−1 A(I e† (I + E)PA . Esto implica que tercambiados para obtener A† = PAT (I + F )A e† k2 (1 + kEk2 ). kA† k2 ≤ (1 + kF k2 )kA Recientemente en [9, Secci´on 4] fueron presentadas cotas de perturbaciones multiplicativas para la pseudoinversa de Moore-Penrose. Las cotas presentadas en [9] e† como las del Teorema 3.2.2, son obtenidas no est´an basadas en expresiones para A siguiendo un m´etodo diferente. En el resto de esta secci´on comparamos las cotas presentadas en el Teorema 3.2.5 con las propuestas en [9]. Escribimos las cotas multiplicativas 2 dadas en [9, Teoremas 4.1 y 4.2] con la misma notaci´on del Teorema 3.2.5: e† − A† k kA b + kF k + kFbk , ≤ kEk + kEk e† k2 } m´ax{kA† k2 , kA r q e† − A† k2 3 kA b 22 + kF k22 + kFbk22 . kEk22 + kEk ≤ † † e 2 m´ax{kA k2 , kA k2 }

(3.2.11) (3.2.12)

e† k2 } dificulta comparar las cotas presenNotemos que la presencia del m´ax{kA† k2 , kA tadas en (3.2.11)-(3.2.12) con las cotas presentadas en el Teorema 3.2.5. Por ejemplo, e† k2 , entonces la cota en (3.2.11) es obviamente mejor que la del Teosi kA† k2 = kA e† k2 , entonces (3.2.11) no da informaci´on sobre rema 3.2.5-(a), pero si kA† k2  kA e† − A† k/kA† k2 , mientras que el Teorema 3.2.5-(a) si. Sin embargo, como discutikA remos a continuaci´on, las cotas en el Teorema 3.2.5 son mejores que las presentadas en (3.2.11)-(3.2.12) ambas a primer orden. Si consideramos perturbaciones muy peque˜ nas e ignoramos los te´rminos de segune† k2 } y m´ın{kA† k2 , kA e† k2 }, do orden, entonces podemos reemplazar m´ax{kA† k2 , kA simplemente por kA† k2 , esto permite hacer comparaciones f´acilmente. El Teorema 3.2.5-(a) y la ecuaci´on (3.2.11) proporcionan la misma cota a primer orden, es decir, e† − A† k/kA† k2 ≤ 2(kEk + kF k). Sin embargo, a primer orden, el lado derecho kA √ p de la ecuaci´on (3.2.12) es Ξ = 3 kEk22 + kF k22 y la cota en el Teorema 3.2.5-(b) p es Γb = kEk22 + kF k22 + (kEk2 + kF k2 )2 . Para comparar Γb y Ξ, utilizamos que (x + y)2 ≤ 2(x2 + y 2 ), para n´ umeros reales x e y tales que x ≥ 0 e y ≥ 0 . As´ı que q Γb ≤ 3 kEk22 + 3 kF k22 = Ξ, implica, que a primer orden, la cota en el Teorema 3.2.5-(b) es siempre mejor que la presentada en (3.2.12). 2

la cota dada en (3.2.12) es presentada en [9] para la clase de normas unitariamente invariantes llamadas Q-normas, las cuales incluyen, entre otras, la norma espectral y la de Frobenius.

47

Cap´ıtulo 3 e† k2 } Para perturbaciones suficientemente grandes, la presencia del m´ax{kA† k2 , kA hace que las ecuaciones (3.2.11) y (3.2.12) sean inaplicables en ciertas situaciones, ya e† − A† k que uno de los objetivos est´andar de la teor´ıa de perturbaciones es acotar kA e† y teniendo s´olamente algunas cotas sobre las normas de las perturbasin conocer A ciones E y F . Ilustramos esto con un ejemplo. Sean A = [1 0; 0 1; 0 0] ∈ R3×2 (hemos utilizado la notaci´on de MATLAB para las matrices), E = diag(−4/5, −4/5, −4/5) y F = diag(−4/5, −4/5). Un c´alculo sencillo muestra que e† − A† k2 kA = 24, kA† k2

e† − A† k2 kA = 0,96 e† k2 kA

b 2 = 4, kEk2 = 0,8, kF k2 = 0,8, kEk

kFbk2 = 4,

los cuales dan: 9,6 para la cota en (3.2.11); 32,64 para la cota en el Teorema 3.2.5-(a); 7,07 para la cota en (3.2.12); 24,03 para la primera cota del Teorema 3.2.5-(b); y 6,084 para la segunda cota del Teorema 3.2.5-(b). Por lo tanto, en este ejemplo, las e† − A† k2 /kA† k2 , mientras ecuaciones (3.2.11)-(3.2.12) fallan al dar una cota para kA que las cotas en el Teorema 3.2.5 dan estimaciones fuertes (en particular aquellas en el Teorema 3.2.5-(b)).

3.3.

Perturbaciones multiplicativas para problemas de m´ınimos cuadrados

En esta secci´on consideramos el problema de m´ınimos cuadrados m´ın kAx − bk2 ,

x∈Rn

A ∈ Rm×n ,

b ∈ Rm ,

(3.3.1)

y el problema de m´ınimos cuadrados perturbado e − ebk2 , m´ınn kAx

x∈R

e = (I + E)A(I + F ) ∈ Rm×n , eb = b + h ∈ Rm , A

(3.3.2)

donde (I + E) ∈ Rm×m e (I + F ) ∈ Rn×n son matrices no singulares. Nuestro objetivo es encontrar una cota superior para el cambio relativo ke x0 − x0 k2 /kx0 k2 , † †e e donde x0 = A b y x e0 = A b son respectivamente las soluciones de m´ınima longitud de (3.3.1) y (3.3.2). Tambi´en examinamos el cambio de los residuos asociados r = b−Ax0 ex0 . El Teorema 3.3.1 es el resultado principal en esta secci´on. y re = eb − Ae Teorema 3.3.1 Sean x0 y x e0 respectivamente las soluciones de m´ınima longitud de ex0 los correspondientes residuos. (3.3.1) y (3.3.2), y sean r = b − Ax0 y re = eb − Ae 48

Teor´ıa de perturbaciones multiplicativas y soluciones precisas del LSP q −1 −1 b b b 22 y αF := Sea E = (I + E) E y F = (I + F ) F , definamos αE := kEk22 + kEk q kF k22 + kFbk22 , y asumamos que khk2 ≤ kbk2 . Entonces las siguientes cotas se cumplen ke x0 − x0 k2 ≤ αF kx0 k2 + [αE (1 + αF ) (1 + ) +  (1 + αF ) ] kA† k2 kbk2 y (3.3.3) ke r − rk2 ≤ kbk2

q ( + kEk2 )2 + kEk22 .

(3.3.4)

Demostraci´on. Veamos primero (3.3.3). La demostraci´on est´a basada en el Corolario 3.2.4, que implica: e† (b + h) − A† b x e0 − x0 = A   e† − A† (b + h) + A† h = A
= A† ΘE + ΘF A† + ΘF A† ΘE (b + h) + A† h
= A† ΘE + ΘF A† ΘE (b + h) + ΘF x0 + ΘF A† h + A† h . Tomando norma 2 en ambos lados y teniendo en cuenta las propiedades de las normas, tenemos ke x0 −x0 k2 ≤ kΘF k2 kx0 k2 +[kΘE k2 (1 + kΘF k2 ) (1 + ) +  (1 + kΘF k2 )] kA† k2 kbk2 . (3.3.3) se sigue del Lema 3.1.1 que implica kΘE k2 ≤ αE y kΘF k2 ≤ αF . Ahora, demostremos (3.3.4). Primero, observemos que ex re − r = h − A e0 + Ax0 eA e† eb + Ax0 =h−A eA e† ) h + Ax0 − A eA e† b = (I − A eA e† ) h + Ax0 − A eA e† (r + Ax0 ) = (I − A eA e† ) (h + Ax0 ) − A eA e† r . = (I − A

(3.3.5)

eA e† , Observemos que los sumandos en (3.3.5) son vectores ortogonales, ya que PAe = A recordemos que Ax0 = PA b, r = (I − PA )b y el Teorema 3.2.3, para obtener (3.3.4) como a continuaci´on ke r − rk22 = k(I − PAe) (h + PA b)k22 + kPAe (I − PA )bk22
2 ≤ khk2 + k(I − PAe)PA bk2 + kPAe (I − PA )bk22 ≤ (kbk2 + kEk2 kbk2 )2 + kEk22 kbk22 . 49

Cap´ıtulo 3

Si A tiene rango completo por filas o columnas la cota (3.3.3) se simplifica en el sentido explicado en las partes (d) y (e) de la Observaci´on 3.2.6. Si A tiene rango completo por filas, entonces re = r = 0 y ke r − rk2 = 0. Las cotas en el Teorema 3.3.1 mejoran significativamente las cotas cl´asicas para el cambio relativo de las soluciones de m´ınima longitud y para los residuos del proe = A+∆A [74, Teorema blema de m´ınimos cuadrados bajo perturbaciones aditivas A 5.1]. Con el objetivo de comparar, enunciamos estas cotas cl´asicas de perturbacioe y nes como aparecen en [5, Teorema 1.4.6] para el caso rango (A) = rango (A) η := κ2 (A) k∆Ak2 /kAk2 < 1. Sea x e0 la soluci´on m´ınima del problema de m´ınimos cuadrados m´ınx∈Rn k(b + ∆b) − (A + ∆A)xk2 y x0 la soluci´on m´ınima de m´ınx∈Rn kb − Axk2 , y sean re := (b + ∆b) − (A + ∆A)e x0 y r := b − Ax0 . Entonces, si suponemos que x0 6= 0 y definimos A := k∆Ak2 /kAk2 y b := k∆bk2 /kbk2 , por el Teorema 1.4.6 de [5] tenemos   1 kA† k2 kbk2 krk2 ke x0 − x0 k2 2 ≤ 2 κ2 (A) A + b + κ2 (A) A y kx0 k2 1−η kx0 k2 kAk2 kx0 k2 (3.3.6)   kAk2 kx0 k2 krk2 ke r − rk2 ≤ A + b + κ2 (A) A . (3.3.7) kbk2 kbk2 kbk2 En la ecuaci´on (3.3.7), es conveniente recordar que (kAk2 kx0 k2 )/kbk2 ≤ κ2 (A) y krk2 ≤ kbk2 . Luego, observemos que la cota para ke r − rk2 /kbk2 en (3.3.7) incluye t´erminos que pueden ser muy grandes aunque A y b sean muy peque˜ nos. Esto sucede si κ2 (A) es grande y si krk2 6= 0 no es muy peque˜ no. En cambio, si kEk2 y  son muy peque˜ nos, entonces la cota para ke r − rk2 /kbk2 que se obtiene de (3.3.4) siempre es muy peque˜ na. Con respecto a las cotas para ke x0 − x0 k2 /kx0 k2 : la cota en (3.3.6) aumenta las perturbaciones en los datos al menos por un factor κ2 (A) y el aumento puede ser m´as grande bajo ciertas condiciones. Adicionalmente, (3.3.6) incluye el factor kA† k2 kbk2 /kx0 k2 , el cual es el u ´nico factor potencialmente grande en la cota que se obtiene de (3.3.3). Mostraremos en la Secci´on 3.3.1 que kA† k2 kbk2 /kx0 k2 es un n´ umero moderado excepto para elecciones muy particulares de b. Por lo tanto, (3.3.3) siempre mejora (3.3.6) y, si kEk2 , kF k2 , y  son muy peque˜ nos, entonces (3.3.3) produce cotas muy peque˜ nas para ke x0 − x0 k2 /kx0 k2 para casi todo b. Las cotas en el Teorema 3.3.1 no pueden ser aplicadas directamente debido a b y Fb. El Corolario 3.3.2 resuelve esta deficiencia restringiendo la la presencia de E magnitud de las perturbaciones y usando la ecuaci´on (3.2.10). El Corolario 3.3.2 se 50

Teor´ıa de perturbaciones multiplicativas y soluciones precisas del LSP demuestra directamente del Teorema 3.3.1 y es enunciado de manera conveniente para su uso en la Secci´on 3.4. Corolario 3.3.2 Con la misma notaci´ on e hip´ otesis del Teorema 3.3.1, y suponiendo que kEk2 ≤ µ < 1, kF k2 ≤ ν < 1, x0 = 6 0 y b 6= 0. Definamos s s 1 1 θµ := µ 1 + y θν := ν 1 + . (3.3.8) 2 (1 − µ) (1 − ν)2 Entonces las siguientes cotas se cumplen kA† k2 kbk2 ke x0 − x0 k2 ≤ θν + [θµ (1 + θν )(1 + ) + (1 + θν )] , kx0 k2 kx0 k2 p ke r − rk2 ≤ ( + µ)2 + µ2 . kbk2 La cota (3.3.9) se cumple a primer orden en , µ y ν  √  kA† k2 kbk2 ke x0 − x0 k2 √ ≤ 2ν +  + 2µ + h.o.t , kx0 k2 kx0 k2

(3.3.9) (3.3.10)

(3.3.11)

donde h.o.t significa t´erminos de orden superior (high order terms) en , µ y ν. Probaremos en la Secci´on 3.3.2 que, a primer orden, la cota de perturbaciones para ke x0 − x0 k2 /kx0 k2 que se obtiene del Teorema 3.3.1 es o´ptima. Es decir, puede ser alcanzada salvo una constante. En este contexto, es interesante observar que otro procedimiento para obtener cotas de perturbaciones multiplicativas para el problema de m´ınimos cuadrados es utilizarel Teorema 3.2.5, como describimos a continuaci´on: e† − A† (b+h)+A† h, as´ı que ke e† −A† k2 (kbk2 + sabemos que x e0 −x0 = A x0 −x0 k2 ≤ kA khk2 )+kA† k2 khk2 , esto puede ser combinado con el Teorema 3.2.5-(b) y khk2 ≤ kbk2 para obtener una cota para ke x0 − x0 k2 /kx0 k2 . La cota obtenida (la llamaremos ΓLP ) incluye en todos sus t´erminos el factor kA† k2 kbk2 /kx0 k2 . Esto contrasta con la cota obtenida de (3.3.3) (lo llamaremos ΓLS ) la cual tiene un t´ermino que es simplemente √ αF (ver tambi´en los t´erminos θν en (3.3.9) ´o 2 ν en (3.3.11)). As´ı, ΓLP es mucho mayor que ΓLS si m´ax{kEk2 , }  kF k2 y kA† k2 kbk2 /kx0 k2 es grande (esto se puede comprobar f´acilmente a primer orden). Adicionalmente, podemos demostrar a primer orden que siempre se cumple que ΓLS ≤ 2 ΓLP y si kA† k2 kbk2 /kx0 k2 ≥ 2 entonces ΓLS ≤ ΓLP . Entonces, podemos concluir que el m´etodo considerado en el Teorema 3.3.1 es mejor que el usado en el Teorema 3.2.5. Finalmente, observemos que todos los resultados en esta secci´on, igual que en la Secci´on 3.2, son v´alidos para cualquier valor de m y n, esto es, si m ≥ n o´ m < n. As´ı que, tambi´en son v´alidos para perturbaciones multiplicativas de soluciones de sistemas lineales infradeterminados. 51

Cap´ıtulo 3

3.3.1.

¿Por qu´ e el factor kA† k2 kbk2 /kx0 k2 es peque˜ no?

Esta secci´on est´a relacionada con [32, Secci´on 3.2] (ver tambi´en Secci´on 2.3), que considera el mismo problema para una matriz no singular A. Aunque el hecho que A ∈ Rm×n sea rectangular, obliga a realizar modificaciones no triviales, las conclusiones m´as importantes permanecen igual. El Teorema 3.3.1 y el Corolario 3.3.2 prueban que la sensibilidad de la soluci´on de m´ınima longitud, x0 = A† b del problema de m´ınimos cuadrados bajo perturbaciones multiplicativas est´a regida por kA† k2 kbk2 /kx0 k2 . Esta cantidad es justamente el n´ umero de condici´on del problema de m´ınimos cuadrados, cuando s´olamente el lado derecho b de m´ınx∈Rn kAx − bk2 es perturbado, como lo demostramos en la ecuaci´on (2.3.9). M´as precisamente, es f´acil demostrar que si x0 = A† b, entonces   ke x0 − x0 k 2 kA† k2 kbk2 † = l´ım sup : x e0 = A (b + h), khk2 ≤  kbk2 . →0 kx0 k2 kx0 k2 Por lo tanto, el Teorema 3.3.1 prueba b´asicamente que las perturbaciones multiplicativas tienen un efecto sobre la soluci´on de m´ınima longitud del problema de m´ınimos cuadrados similar a perturbar s´olamente el lado derecho b. Notemos que 1 ≤ kA† k2 kbk2 /kx0 k2 , pero, en general, kA† k2 kbk2 /kx0 k2  κ2 (A), contrastando con el caso cuando A es no singular3 [32, Secci´on 3.2]. Sin embargo, (3.3.6) muestra que   kA† k2 kbk2 krk2 2 κLS (A, b) := 2 κ2 (A) + + κ2 (A) kx0 k2 kAk2 kx0 k2 puede ser considerado como un n´ umero de condici´on para el problema de m´ınimos cuadrados bajo perturbaciones aditivas muy peque˜ nas de A y b (en efecto, se ha demostrado en [74, Secci´on 6] que la cota (3.3.6) es aproximadamente alcanzada a primer orden en las perturbaciones), y kA† k2 kbk2 /kx0 k2 ≤ κLS (A, b). Pero el primer punto clave en esta secci´on es mostrar que si A es fija, entonces kA† k2 kbk2 /kx0 k2 es un n´ umero moderado para la mayor´ıa de los vectores b, incluso si κ2 (A)  1 y as´ı κLS (A, b)  1, lo cual implica que kA† k2 kbk2 /kx0 k2  κLS (A, b) para la mayor´ıa de los problemas de m´ınimos cuadrados mal condicionados cuya matriz de coeficientes es A. Sin embargo, esto no es suficiente para nuestro prop´osito, ya que si rango(A) < m, entonces para la mayor´ıa de vectores b el ´angulo agudo θ(b, R(A)) entre b y el espacio columna de A no es peque˜ no, lo cual es equivalente a decir que el residuo Si A es no singular o, si Ax0 = b, entonces kA† k2 kbk2 /kx0 k2 ≤ κ2 (A). Sin embargo, consideremos A = [1 0; 0 1; 0 0] ∈ R3×2 y b = [η; 0; 1] ∈ R3×1 . En este caso x0 = [η; 0] ∈ R2×1 y p kA† k2 kbk2 /kx0 k2 = |η|2 + 1/|η| tiende a ∞ si η → 0 mientras κ2 (A) = 1. 3

52

Teor´ıa de perturbaciones multiplicativas y soluciones precisas del LSP relativo kAx0 −bk2 /kbk2 = sin θ(b, R(A)) no es peque˜ no. Pero muy frecuentemente en la pr´actica, los problemas de m´ınimos cuadrados tienen residuos relativos peque˜ nos, dado que los problemas corresponden a sistemas lineales inconsistentes Ax ≈ b que est´an cercanos de ser consistentes. Por lo tanto, el segundo punto importante en esta secci´on es que, si fijamos A para considerar todos los vectores b tales que Υ = θ(b, R(A)) < π/2 tambi´en fijo, podemos demostrar que para la mayor´ıa de estos vectores b el factor kA† k2 kbk2 /kx0 k2 es un n´ umero moderado mucho m´as peque˜ no que κLS (A, b) cuando A est´a muy mal condicionada. Para explicar las propiedades mencionadas anteriormente, supongamos que rango (A) = r y sea A = U ΣV T la SVD de A, donde U ∈ Rm×r y V ∈ Rn×r tienen columnas ortonormales, Σ = diag(σ1 , . . . , σr ) ∈ Rr×r , y σ1 ≥ · · · ≥ σr > 0. Observemos que kx0 k2 = kA† bk2 = kΣ−1 U T bk2 ≥ |uTr b|/σr y 1 kbk2 kbk2 kA† k2 kbk2 , = ≤ T = kx0 k2 σr kx0 k2 |ur b| cos θ(ur , b)

(3.3.12)

donde ur es la u ´ltima columna de U y θ(ur , b) es el a´ngulo agudo entre ur y b. Note que la cota sobre kA† k2 kbk2 /kx0 k2 en (3.3.12) puede ser grande s´olo si b es “casi” ortogonal a ur . Por ejemplo, si A es una matriz fija extremadamente mal condicionada (por ejemplo κ2 (A) = 101000 ) y b es considerado como un vector aleatorio cuya direcci´on est´a uniformemente distribuida en todo el espacio, entonces la probabilidad que 0 ≤ θ(ur , b) ≤ π/2 − 10−6 es aproximadamente 1 − 10−6 . Observemos que la condici´on 0 ≤ θ(ur , b) ≤ π/2 − 10−6 implica que kA† k2 kbk2 /kx0 k2 . 106 , lo cual es un n´ umero moderado comparado a 101000 . En particular, si los par´ametros perturbativos µ, ν y  en el Corolario 3.3.2 son 10−16 , entonces kA† k2 kbk2 /kx0 k2 . 106 proporciona una excelente cota para la variaci´on de la soluci´on de m´ınima longitud del problema de m´ınimos cuadrados. Incluso, es posible que kA† k2 kbk2 /kx0 k2 sea moderado aunque cos θ(ur , b) ≈ 0. Esto puede verse como una extensi´on del caso de matrices no singulares a matrices generales del resultado original presentado por Chan y Foulser (Teorema 2.3.2). En el argumento anterior, el vector b puede estar en cualquier lugar en el espacio. Luego, consideremos vectores b tales que Υ = θ(b, R(A)) < π/2 se mantiene constante. Describiraremos todos estos vectores a continuaci´on: sea y ∈ Rr cualquier vector y sea U⊥ ∈ Rm×(m−r) tal que [ U U⊥ ] ∈ Rm×m es unitaria. Entonces elegimos cualquier z ∈ Rm−r tal que kzk2 = kyk2 tan Υ, y definimos b = U y + U⊥ z. Sabemos que Υ = θ(b, R(A)), ya que R(U ) = R(A). Adem´as, de (3.3.12), es f´acil demostrar 53

Cap´ıtulo 3 que estos vectores b satisfacen p kA† k2 kbk2 kbk2 kbk2 1 + tan2 Υ 1 = ≤ T = = , (3.3.13) kx0 k2 σr kx0 k2 |ur b| cos θ(er , y) (cos Υ) · (cos θ(er , y) ) donde er es la r-´esima columna de Ir . La cota en (3.3.13) es una cantidad “geom´etrica” que no depende de κ2 (A), si suponemos que Υ no es muy cercano a π/2, y que es un n´ umero moderado para la mayor´ıa de los vectores y, es decir, para la mayor´ıa de los vectores b tales que Υ = θ(b, R(A)). Finalmente, discutimos una interesante relaci´on entre kA† k2 kbk2 /kx0 k2 con el t´ermino de κLS (A, b) que depende de κ22 (A). Note que este t´ermino puede acotarse como sigue Φ := κ22 (A)

kA† k2 krk2 kA† k2 kbk2 krk2 = κ2 (A) ≤ κ2 (A) . kAk2 kx0 k2 kx0 k2 kx0 k2

(3.3.14)

En relaci´on a nuestra discusi´on kA† k2 kbk2 /kx0 k2 es un n´ umero moderado para la mayor´ıa de los vectores b. Por lo tanto, Φ es acotado superiormente por un n´ umero moderado de veces de κ2 (A) para la mayor´ıa de los vectores b y como consecuencia, κ22 (A) s´olamente afecta la sensibilidad del problema de m´ınimos cuadrados en situaciones muy particulares. Adem´as, Φ puede escribirse como:   kA† k2 kbk2 krk2 κ2 (A) = Φ, kx0 k2 kbk2

(3.3.15)

lo cual implica que para residuos relativos suficientemente grandes (pensar, por ejemplo, en krk2 /kbk2 ≥ 10−3 ) y matrices A muy mal condicionadas, tenemos kA† k2 kbk2 /kx0 k2  Φ ≤ κLS (A, b), incluso si kA† k2 kbk2 /kx0 k2 es grande.

3.3.2.

El n´ umero de condici´ on para perturbaciones multiplicativas del problema de m´ınimos cuadrados

En esta secci´on demostramos que kA† k2 kbk2 /kx0 k2 es el n´ umero de condici´on bajo perturbaciones multiplicativas del problema de m´ınimos cuadrados salvo una constante. El lector debe notar que, por simplicidad, consideramos en nuestra definici´on de n´ umero de condici´on que la variaci´on relativa de b y las perturbaciones multiplicativas izquierdas y derechas tienen el mismo orden.

54

Teor´ıa de perturbaciones multiplicativas y soluciones precisas del LSP Teorema 3.3.3 Usando la misma notaci´ on e hip´ otesis del Corolario 3.3.2 con los par´ametros µ, ν y  iguales a η, definamos el n´ umero de condici´ on  ke x0 − x0 k 2 (M ) κLS (A, b) := l´ım sup : x e0 = [(I + E)A(I + F )]† (b + h), η→0 η kx0 k2  kEk2 ≤ η, kF k2 ≤ η, khk2 ≤ ηkbk2 . Entonces, 1 kA† k2 kbk2 ) (M ) √ κ(M (A, b) ≤ ≤ κLS (A, b) . LS kx0 k2 1+2 2

(3.3.16)

Demostraci´on. De la ecuaci´on (3.3.11) y de 1 ≤ kA† k2 kbk2 /kx0 k2 , tenemos  √  kA† k2 kbk2 √  kA† k2 kbk2 ke x0 − x0 k 2 √  ≤ 2+ 1 + 2 + O(η) ≤ 1 + 2 2 + O(η) , η kx0 k2 kx0 k2 kx0 k2 lo cual implica la desigualdad izquierda de (3.3.16). Para demostrar la desigualdad del lado derecho elegimos una perturbaci´on tal que E = 0, F = 0 y h = ηw, donde kwk2 = kbk2 y kA† wk2 = kA† k2 kwk2 . Para esta perturbaci´on ke x0 − x0 k2 /kx0 k2 = † † kA hk2 /kx0 k2 = η kA k2 kbk2 /kx0 k2 . Por lo tanto, la expresi´on “sup” que aparece en (M ) (M ) la definici´on de κLS (A, b) implica que kA† k2 kbk2 /kx0 k2 ≤ κLS (A, b).

3.3.3.

Cotas de perturbaciones multiplicativas para otras soluciones del problema de m´ınimos cuadrados

Cotas para el cambio en las soluciones de problemas de m´ınimos cuadrados se pueden obtener f´acilmente a partir del Teorema 3.3.1 y Teorema 3.2.3-(c) que son una peque˜ na modificaci´on de la ecuaci´on (3.3.3). Dado que el residuo de un problema de m´ınimos cuadrados es el mismo para todas sus soluciones, no es necesario considerar nuevas cotas perturbativas para el residuo. Teorema 3.3.4 Si y ∈ Rn es una soluci´ on del problema de m´ınimos cuadrados (3.3.1), entonces existe una soluci´ on ye ∈ Rn del problema de m´ınimos cuadrados (3.3.2) tal que ke y − yk2 ≤ (αF + kF k2 ) kyk2 + [αE (1 + αF ) (1 + ) +  (1 + αF ) ] kA† k2 kbk2 , donde αE , αF y  se definen como en el enunciado del Teorema 3.3.1. 55

Cap´ıtulo 3 Demostraci´on. Dado y, existe un vector z ∈ Rn tal que y = x0 + PN (A) z, donde x0 es la soluci´on de m´ınima longitud de (3.3.1). Recordemos tambi´en que kyk22 = kx0 k22 + kPN (A) zk22 y por lo tanto, kPN (A) zk2 ≤ kyk2 . Elegimos la siguiente soluci´on de (3.3.2), ye = x e0 + PN (A) e0 es la soluci´on m´ınima de (3.3.2). Por lo e PN (A) z, donde x tanto ke y − yk2 ≤ ke x0 − x0 k2 + k(PN (A) x0 − x0 k2 + kF k2 kyk2 , e − PN (A) )PN (A) zk2 ≤ ke donde hemos utilizado el Teorema 3.2.3-(c). Ahora, usamos (3.3.3) y kx0 k2 ≤ kyk2 para obtener el resultado. Observemos que el cambio relativo ke y − yk2 /kyk2 est´a acotado por m´ax{1, kA† k2 kbk2 /kyk2 }, que es menor o igual que kA† k2 kbk2 /kx0 k2 . Por lo tanto, la soluci´on de m´ınima longitud es la m´as sensible de las soluciones bajo perturbaciones multiplicativas.

3.4.

Perturbaci´ on del problema de m´ınimos cuadrados en los factores

Como explicamos en la Introducci´on, presentamos en la Secci´on 3.5 un algoritmo preciso para la soluci´on del problema de m´ınimos cuadrados m´ınx∈Rn kb − Axk2 utilizando una RRD precisa XDY de A, el Algoritmo 3.5.1. El an´alisis de errores del Algoritmo 3.5.1 lo presentamos en el Teorema 3.5.2 y muestra que la soluci´on calculada es la soluci´on exacta del problema de m´ınimos cuadrados correspondiente a una RRD con factores cercanos a (X + ∆X)(D + ∆D)(Y + ∆Y ), donde las perturbaciones son “normwise” para los factores bien condicionados X e Y , y “componentwise” para los elementos diagonales y potencialmente mal condicionados del factor D. Por lo tanto, es necesario desarrollar cotas perturbativas para la soluci´on del problema de m´ınimos cuadrados cuya matriz de coeficientes est´a expresada como una RRD bajo perturbaciones en los factores. Esto es presentado en el Teorema 3.4.1, y su demostraci´on se basa en escribir las perturbaciones en los factores como una perturbaci´on multiplicativa de toda la matriz. Recordemos que, por el Lema 3.1.2-(c), si A = XDY es una RRD de A, entonces A† = Y † D−1 X † . Por lo tanto, la soluci´on de m´ınima longitud del problema de m´ınimos cuadrados m´ınx∈Rn kb − XDY xk2 es x0 = Y † D−1 X † b. Teorema 3.4.1 Sean X ∈ Rm×r , D ∈ Rr×r e Y ∈ Rr×n tales que rango (X) = rango (Y ) = r, D es diagonal y no singular y b ∈ Rm . Sea x0 la soluci´ on de 56

Teor´ıa de perturbaciones multiplicativas y soluciones precisas del LSP m´ınima longitud de m´ınx∈Rn kb − XDY xk2 y x e0 la soluci´ on de m´ınima longitud de m´ınx∈Rn k(b + h) − (X + δX)(D + δD)(Y + δY ) xk2 , donde kδXk2 ≤ αkXk2 , kδY k2 ≤ βkY k2 , |δD| ≤ ρ|D|, y khk2 ≤ kbk2 . Sean r = b − XDY x0 y re = (b + h) − (X + δX)(D + δD)(Y + δY ) x e0 . Supongamos que µ := α κ2 (X) < 1 y

ν := [β + ρ(1 + β)]κ2 (Y ) < 1,

(3.4.1)

y definamos los par´ametros θµ y θν como en la ecuaci´ on (3.3.8). Entonces, la cota † † −1 † (3.3.9) se satisface con A reemplazada por Y D X , la cota (3.3.10) se satisface y a primer orden en α, β, ρ y  tenemos   kY † D−1 X † k kbk √ ke x0 − x0 k2 √ 2 2 ≤ 2 (β+ρ) κ2 (Y )+  + 2 α κ2 (X) +h.o.t. (3.4.2) kx0 k2 kx0 k2 e = (X + δX)(D + δD)(Y + δY ). Escribamos Demostraci´on. Sean A = XDY y A e como una perturbaci´on multiplicativa de A: A e = (I + δXX † ) XD (I + D−1 δD)Y (I + Y † δY ) A = (I + δXX † ) XDY (I + Y † D−1 δD Y ) (I + Y † δY ) =: (I + E)A(I + F ), donde E = δXX † y F = Y † δY + Y † D−1 δD Y + Y † D−1 δD δY . Luego, tomando en cuenta que kδD D−1 k2 ≤ ρ, tenemos kEk2 ≤ α κ2 (X) = µ < 1,

kF k2 ≤ [β + ρ(1 + β)] κ2 (Y ) = ν < 1,

y el Teorema 3.4.1 se obtiene inmediatamente del Corolario 3.3.2. Dado que los factores X e Y de una RRD est´an bien condicionados, de la ecuaci´on (3.4.2) tenemos que la sensitividad respecto a perturbaciones de los factores de la soluci´on de m´ınima longitud del problema de m´ınimos cuadrados m´ınx∈Rn kb−XDY xk2 est´a nuevamente controlada por kA† k2 kbk2 /kx0 k2 , donde A = XDY , que es un n´ umero moderado para la mayor´ıa de los vectores b (ver Secci´on 3.3.1). Note que el Teorema 3.4.1 es v´alido si m ≥ n o si m < n, es decir, para problemas de m´ınimos cuadrados o para sistemas de ecuaciones lineales infradeterminados, ya que el Corolario 3.3.2 se cumple en ambos casos.

3.5.

Algoritmo y an´ alisis de errores

En esta secci´on presentamos el Algoritmo 3.5.1 para resolver el problema de m´ınimos cuadrados m´ınx∈Rn kb − A xk2 y demostramos que calcula la soluci´on de 57

Cap´ıtulo 3 m´ınima longitud con error relativo acotado por O(u) kA† k2 kbk2 /kx0 k2 , lo cual es simplemente O(u) para la mayor´ıa de vectores b de acuerdo con lo discutido en la Secci´on 3.3.1. El primer paso del algoritmo calcula una RRD precisa de A = XDY ∈ Rm×n en el sentido de la Definici´on 2.4.2, lo cual es posible para muchas clases de matrices estructuradas como hemos dicho en la Introducci´on. Los siguientes pasos del Algoritmo 3.5.1 est´an basados en el hecho de que, acorde al Lema 3.1.2-(c), la soluci´on de m´ınima longitud es x0 = Y † (D−1 (X † b)) y en las siguientes observaciones: (1) s = X † b es la soluci´on u ´nica del problema de m´ınimos cuadrados de rango completo por columnas m´ınx∈Rr kb − X xk2 ; (2) w = D−1 (X † b) es la soluci´on u ´nica † −1 † del sistema lineal Dw = s; y (3) Y (D (X b)) es la soluci´on u ´nica del sistema lineal indeterminado con rango completo por filas Y x = w. Observemos que este m´etodo es v´alido si m ≥ n y si m < n. Por lo tanto, en el u ´ltimo caso y, si rango (A) = m, el m´etodo resuelve de manera precisa el sistema lineal infradeterminado Ax = b. La soluci´on m´ınima x0 del sistema infraterminado Y x = x2 es calculado v´ıa el QM´etodo descrito en [43, Secci´on 21.1] y que recordamos brevemente a continuaci´on. La idea es calcular a primer orden una factorizaci´on QR reducida de Y T = W RY ∈ Rn×r utilizando el algoritmo de Householder, donde W ∈ Rn×r satisface W T W = Ir y RY ∈ Rr×r es triangular superior y no singular. As´ı, Y = RYT W T ∈ Rr×n y el Lema 3.1.2-(c) implican que Y † = (W T )† (RYT )† = W RY−T , donde RY−T denota la inversa de RYT . Finalmente, x0 = W (RY−T x2 ) y RY−T x2 es calculado resolviendo el sistema triangular RYT x = x2 por sustituci´on progresiva. En la pr´actica, es importante notar que el factor W no se requiere expl´ıcitamente, s´olamente necesitamos la habilidad de multiplicar W veces un vector y esto puede hacerse multiplicando las permutaciones de Householder resultantes de la factorizaci´on QR de Y T por [(RY−T x2 )T , 0]T ∈ Rn . Estamos ahora en posici´on de enunciar el Algoritmo 3.5.1.

Algoritmo 3.5.1 (Soluci´on precisa del problema de m´ınimos cuadrados v´ıa la RRD) Input: A ∈ Rm×n , b ∈ Rm Output: x0 la soluci´on de m´ınima longitud de m´ınx∈Rn kb − A xk2

Paso 1: Calcular una RRD precisa de A = XDY en el sentido de la Definici´on 2.4.2, donde X ∈ Rm×r , D ∈ Rr×r es diagonal, Y ∈ Rr×n y rango (A) = rango (X) = rango (Y ) = rango (D) = r.

58

Teor´ıa de perturbaciones multiplicativas y soluciones precisas del LSP Paso 2: Calcular la soluci´on u ´nica s de m´ınx∈Rr kb − X xk2 usando la factorizaci´on QR de Householder de X. Paso 3: Calcular la soluci´on u ´nica w del sistema lineal diagonal D w = s como wi = si /dii , i = 1, . . . , r. Paso 4: Calcular la soluci´on de m´ınima longitud x0 de Y x = w usando el Q-m´etodo, es decir, v´ıa la factorizaci´on QR de Householder de Y T .

Antes de realizar el an´alisis de errores del Algoritmo 3.5.1, observemos que por el teorema de perturbaciones de Weyl [69], las diferencias entre los valores singub est´an acotadas como |σi (X) b − σi (X)| ≤ kX b − Xk2 ≤ lares ordenados de X y X b u p(m, n) kXk2 , para i = 1 : r. Por lo tanto, |σi (X)−σ i (X)|/σi (X) ≤ u p(m, n) κ2 (X), para i = 1 : r. An´alogamente podemos acotar las diferencias entre los valores singulares ordenados de Y e Yb . Como consecuencia, la condici´on (2.4.3) implica que b = r, rango (Y ) = rango (Yb ) = r y rango (X) = rango (X) κ2 (X) b ≤ 3 κ2 (X) y ≤ κ2 (X) 3

κ2 (Y ) ≤ κ2 (Yb ) ≤ 3 κ2 (Y ) . 3

(3.5.1)

b La ecuaci´on (3.5.1) nos permite usar indistintamente κ2 (X), κ2 (Y ) ´o κ2 (X), κ2 (Yb ) en las cotas de error obtenidas con el coste de modificar ligeramente algunas de las constantes involucradas. El coste computacional en el Paso 1 del Algoritmo 3.5.1 depende del tipo espec´ıfico de matrices y del algoritmo usado entre los mencionados en la introducci´on. De cualquier manera, todos estos algoritmos tienen un coste de O(mn2 ) flops si m ≥ n y O(m2 n) flops si m < n. El t´ermino principal del coste de los Pasos 2, 3 y 4 es 2r2 (m − r/3), r y 2r2 (n − r/3) flops, respectivamente. Dado que r ≤ m´ın{m, n}, el coste total del Algoritmo 3.5.1 es O(mn2 ) flops si m ≥ n y O(m2 n) flops si m < n. Los errores de redondeo regresivos cometidos por el Algoritmo 3.5.1 son analizados en el Teorema 3.5.2. Usaremos la siguiente notaci´on introducida en [43, Secciones 3.1 y 3.4] γn :=

nu 1 − nu

and γ en :=

cnu , 1 − cnu

(3.5.2)

donde c denota una peque˜ na constante entera cuyo valor exacto no es esencial en el an´alisis. Antes de enunciar y demostrar el Teorema 3.5.2, comentamos la necesidad 59

Cap´ıtulo 3 y el significado de las hip´otesis usadas en el teorema. Primero, asumimos que los b D b e Yb calculados en el Paso 1 en aritm´etica en coma flotante satisfactores X, b = r, facen (2.4.1), (2.4.2) y (2.4.3) , lo cual implica que rango (X) = rango (X) b = r, rango (Y ) = rango (Yb ) = r y (3.5.1). Por lo tanto, rango (D) = rango (D) b y podemos usar κ2 (X) y κ2 (Y ) en los errores de los Pasos 2 y 4 en vez de κ2 (X) κ2 (Yb ) al coste de no prestar atenci´on a los valores exactos de las constantes num´ericas √ er m´ax{m,n} < 1 garantiza en las cotas de errores. La hip´otesis m´ax{κ2 (X), κ2 (Y )} r γ b sobre X b en el Paso 2 preservan el rango completo, es que los errores regresivos ∆X b = rango (X b + ∆X) b = r y lo mismo para los errores regresivos de decir, rango (X) Yb en el Paso 4. Finalmente, la hip´otesis t´ecnica κ2 (Y ) n r2 γ en < 1 se necesita para aplicar [43, Teorema 21.4] en el an´alisis de errores del Paso 4. Presentaremos en el Teorema 3.5.2 dos enunciados para los errores regresivos del b D b e Yb de A y otro Algoritmo 3.5.1, uno con respecto a los factores calculados X, con respecto a los factores exactos, el cual es el resultado usado generalmente. La raz´on para presentar estos dos enunciados es que el primero da errores regresivos por b e Yb m´as fuertes que el segundo. Esto puede ser usado para filas y columnas en X obtener errores regresivos m´as fuertes para algunas clases particulares de matrices, tales como matrices de Cauchy. b ∈ Rm×r , D b ∈ Rr×r e Yb ∈ Rr×n los factores de A calculados Teorema 3.5.2 Sean X en el Paso 1 del Algoritmo 3.5.1 y supongamos que satisfacen las cotas de error dadas en las ecuaciones (2.4.1) y (2.4.2) con respecto a los factores exactos X, D e Y de A. Supongamos tambi´en que (2.4.3), √ m´ax{κ2 (X), κ2 (Y )} r γ er m´ax{m,n} < 1 y (3.5.3) κ2 (Y ) n r2 γ en < 1

(3.5.4)

se cumplen. Sea x b0 la soluci´on de m´ınima longitud calculada del problema de m´ınimos cuadrados m´ınx∈Rn kb−A xk2 usando el Algoritmo 3.5.1 en precisi´ on finita con unidad de redondeo u. Entonces los siguientes enunciados se cumplen: (a) x b0 es la soluci´on de m´ınima longitud de b + ∆X)( b D b + ∆D)( b Yb + ∆Yb ) xk2 , m´ın k(b + ∆b) − (X

x∈Rn

(3.5.5)

donde b j)k2 ≤ γ b j)k2 , k∆X(:, emr kX(:, b ≤γ b |∆D| e1 |D|, 60

k∆Yb (j, :)k2 ≤ γ enr kYb (j, :)k2 , para j = 1, . . . , r k∆bk2 ≤ γ emr kbk2 .

Teor´ıa de perturbaciones multiplicativas y soluciones precisas del LSP (b) x b0 es la soluci´on de m´ınima longitud de m´ın k(b + ∆b) − (X + ∆X)(D + ∆D)(Y + ∆Y ) xk2 ,

x∈Rn

(3.5.6)

donde √ √ rγ emr + r γ emr u p(m, n)) kXk2 , √ √ k∆Y k2 ≤ (u p(m, n) + r γ enr + r γ enr u p(m, n)) kY k2 , k∆Xk2 ≤ (u p(m, n) +

|∆D| ≤ (u p(m, n) + γ e1 + γ e1 u p(m, n)) |D|,

k∆bk2 ≤ γ emr kbk2 .

(c) Si x0 es la soluci´on de m´ınima longitud exacta de m´ınx∈Rn kb − A xk2 , entonces kb x0 −x0 k2 /kx0 k2 puede acotarse como en el Teorema 3.4.1 con α = (u p(m, n)+ √ √ √ √ rγ emr + r γ emr u p(m, n)), β = (u p(m, n) + r γ enr + r γ enr u p(m, n)), ρ = (u p(m, n) + γ e1 + γ e1 u p(m, n)) y  = γ emr . En particular, a primer orden en u, y si c es una constante entera peque˜ na, entonces   kb x0 − x0 k 2 kA† k2 kbk2 ≤ c u py (m, n) κ2 (Y ) + px (m, n) κ2 (X) + O(u2 ) , kx0 k2 kx0 k2 donde py (m, n) := (p(m, n) + nr3/2 ) y px (m, n) := (p(m, n) + mr3/2 ). Demostraci´on. Para demostrar la parte (a) escribamos los errores regresivos en los pasos 2, 3 y 4 del Algoritmo 3.5.1. 1. Los errores regresivos del Paso 2 est´an dados en [43, Teorema 20.3]: la soluci´on calculada en el Paso 2, x b1 , es la soluci´on exacta del problema de m´ınimos cuadrados b + ∆X) b xk2 , m´ın k(b + ∆b) − (X

x∈Rr

(3.5.7)

b j)k2 ≤ γ b j)k2 , para j = 1, . . . , r, y k∆bk2 ≤ γ donde k∆X(:, emr kX(:, emr kbk2 . √ b b b b Por lo tanto, k∆Xk2 ≤ k∆XkF ≤ γ emr kXkF ≤ re γmr kXk2 . Note que, como mencionamos anteriormente, las ecuaciones (2.4.1) y (2.4.3) implib = r, por lo tanto el teorema de perturcan rango (X) = rango (X) b + baciones de Weyl [69] para valores singulares y (3.5.3) implican |σr (X √ b − σr (X)|/σ b b b b b < 1, y finalmente, ∆X) re γmr κ2 (X) r (X) ≤ k∆Xk2 /σr (X) ≤ b = rango (X b + ∆X) b = r. Como consecuencia, x rango (X) b1 satisface b + ∆X) b † (b + ∆b), x b1 = (X

(3.5.8)

b + ∆X b ∈ Rm×r tal que rango (X b + ∆X) b = r. con X 61

Cap´ıtulo 3 2. Como consecuencia de [43, Lema 3.5], la soluci´on, x b2 , calculada en el Paso 3 cumple que b + ∆D) b x (D b2 = x b1

b ≤γ b con |∆D| e1 |D|,

(3.5.9)

b + ∆D b ∈ Rr×r diagonal y no singular, dado que (2.4.2) y (2.4.3) implican con D b =ryγ rango (D) = rango (D) e1 < 1 por (3.5.3). 3. Los errores regresivos del Paso 4 est´an dados en [43, Teorema 21.4]. Para aplicar [43, Teorema 21.4] necesitamos que rango (Yb ) = r, el cual se sigue de (2.4.1) y (2.4.3), y la hip´otesis k |Yb † | |Yb | k2 r n γn < 1, se cumple por la ecuaci´on (3.5.4), dado que k |Yb † | |Yb | k2 r n γn ≤ κ2 (Yb ) r2 n γn < 1. Con esta condici´on, la soluci´on de m´ınima longitud calculada en el Paso 4 x b0 , es la soluci´on m´ınima exacta del sistema lineal indeterminado (Yb + ∆Yb )x = x b2 , con k∆Yb (j, :)k2 ≤ γ enr kYb (j, :)k2 , para j = 1, . . . , r. Adem´as, podemos demostrar que rango (Yb ) = rango (Yb + ∆Yb ) = r utilizando un argumento similar b + ∆X. b Por lo tanto, x al usado para demostrar el resultado an´alogo para X b0 satisface x b0 = (Yb + ∆Yb )† x b2 ,

(3.5.10)

con Yb + ∆Yb ∈ Rr×n y rango (Yb + ∆Yb ) = r. De las ecuaciones (3.5.8), (3.5.9) y (3.5.10) tenemos que b + ∆D) b −1 (X b + ∆X) b † (b + ∆b) x b0 = (Yb + ∆Yb )† (D h i† b + ∆X) b (D b + ∆D) b (Yb + ∆Yb ) (b + ∆b) , = (X

(3.5.11) (3.5.12)

donde la segunda igualdad se obtiene del Lema 3.1.2-(c). Esto y las cotas que hemos b D b e Yb demuestran la parte (a) del Teorema 3.5.2. desarrollado para X, La demostraci´on del Teorema 3.5.2-(b) se obtiene f´acilmente de la parte (a). Las b = X + EX , D b = D + ED e Yb = ecuaciones (2.4.1) y (2.4.2) nos permiten escribir X Y +EY , con kEX k2 ≤ u p(m, n) kXk2 , |ED | ≤ u p(m, n)|D| y kEY k2 ≤ u p(m, n)kY k2 . Por lo tanto, podemos escribir b + ∆X b = X + EX + ∆ X b =: X + ∆X , X 62

(3.5.13)

Teor´ıa de perturbaciones multiplicativas y soluciones precisas del LSP donde b 2 k∆Xk2 ≤ kEX k2 + k∆Xk ≤ u p(m, n) kXk2 +

√

b 2 rγ emr kXk

√ ≤ u p(m, n) kXk2 + r γ emr (kXk2 + kEX k2 ) √ √ ≤ (u p(m, n) + r γ emr + r γ emr u p(m, n)) kXk2 .

(3.5.14)

An´alogamente, b + ∆D b =: D + ∆D , con |∆D| ≤ (u p(m, n) + γ D e1 + γ e1 u p(m, n)) |D| , √ √ enr + r γ enr u p(m, n)) Yb + ∆Yb =: Y + ∆Y , con k∆Y k2 ≤ (u p(m, n) + r γ kY k2 . Si estas u ´ltimas ecuaciones junto con (3.5.13) y (3.5.14) se reemplazan en (3.5.5), entonces obtenemos la ecuaci´on (3.5.6) y as´ı demostramos la parte (b). Finalmente, la parte (c) es una consecuencia inmediata de (b) y el Teorema 3.4.1. Observemos que, como en una RRD los factores X e Y est´an bien condicionados, el Teorema 3.5.2-(c) garantiza que el error progresivo en la soluci´on calculada por el Algoritmo 3.5.1 est´a acotado por O(u)kA† k2 kbk2 /kx0 k2 .

3.6.

Experimentos Num´ ericos

En esta secci´on mostraremos experimentos num´ericos realizados usando MATLAB— los cuales ilustran como son los errores cometidos por el Algoritmo 3.5.1 comparados con las predicciones te´oricas y con los errores cometidos por el m´etodo usual para resolver problemas de m´ınimos cuadrados, es decir, usando la factorizaci´on QR calculada con el algoritmo tradicional de Householder implementado en MATLAB— [43, Secci´on 20.2]. Para esto, usaremos tres clases importante de matrices rectangulares estructuradas que pueden tener n´ umeros de condici´on muy grandes: Matrices de Cauchy, de Vandermonde y Graduadas. Para matrices pertenecientes a estas clases, RRDs precisas en el sentido de la Definici´on 2.4.2 pueden ser calculadas usando los algoritmos dados en [18] y [42]. Presentaremos s´olamente experimentos para matrices A ∈ Rm×n con entradas reales, m ≥ n y tales que rango(A) = n, lo cual significa que consideramos s´olamente problemas de m´ınimos cuadrados con soluci´on u ´nica. Sabemos de (3.0.1) y (3.3.14) que si x b0 es la soluci´on u ´nica de m´ınx∈Rn kAx − bk2 calculada por el algoritmo QR de MATLAB— y x0 es la soluci´on exacta, entonces kb x0 − x0 k2 kA† k2 kbk2 ≤ c m n3/2 u κ2 (A) , kx0 k2 kx0 k2

(3.6.1) 63

Cap´ıtulo 3 la cual es una cota m´as grande que (3.0.1) pero confiable en la mayor´ıa de situaciones. En cambio, el Algoritmo 3.5.1 satisface (ver (3.0.2) y el Teorema 3.5.2-(c))   kb x0 − x0 k2 kA† k2 kbk2 ≤ u f (m, n) κ2 (Y ) + κ2 (X) . (3.6.2) kx0 k2 kx0 k2 En nuestros experimentos, hemos calculado el error relativo en la soluci´on para ambos algoritmos, y tambi´en las cantidades     kA† k2 kbk2 kA† k2 kbk2 ΘQR := u κ2 (A) , ΘRRD := u κ2 (Y ) + κ2 (X) , (3.6.3) kx0 k2 kx0 k2 con el fin de comprobar, la veracidad de las cotas (3.6.1) y (3.6.2). Observamos que el Algoritmo 3.5.1, hasta ahora es el m´as preciso de los dos: para un vector aleatorio b, alcanza errores relativos en norma de orden la unidad de redondeo, lo que significa que ΘRRD es O(u) casi siempre, incluso para matrices A mal condicionadas.

3.6.1.

Matrices de Cauchy

Las componentes de una matriz de Cauchy, C ∈ Rm×n , m ≥ n, son definidas en t´erminos de dos vectores z = [z1 , . . . , zm ]T ∈ Rm e y = [y1 , . . . , yn ]T ∈ Rn como cij =

1 , zi + yj

i = 1, . . . , m, j = 1, . . . , n.

(3.6.4)

Matrices de la forma G = S1 CS2 , donde C es una matriz de Cauchy y, S1 y S2 son matrices diagonales y no singulares, en [18] son llamadas matrices quasi-Cauchy, las cuales incluyen, como un caso particular, las matrices de Cauchy para S1 = Im y S2 = In . Las matrices quasi-Cauchy tienen rango completo por columnas si zi 6= zj para calquier i 6= j, yk 6= yl para cualquier k 6= l y zi 6= −yj para todo i y j. El Algoritmo 3 en [18] utiliza una versi´on estructurada de GECP para calcular una RRD precisa de cualquier matriz quasi-Cauchy cuadrada. Este algoritmo puede ser extendido f´acilmente para el caso de matrices rectangulares, y esta versi´on la usaremos en los experimentos de esta secci´on para calcular la RRD en el Paso 1 del Algoritmo 3.5.1. El coste total de este paso es 2mn2 −2n3 /3+O(n2 +mn) operaciones m´as mn2 /2 − n3 /6 + O(n2 + mn) comparaciones. Con el objetivo de hacer sencillas referencias, resumamos y mencionemos los dos algoritmos que son usados en esta secci´on para resolver m´ınx∈Rn kCx − bk2 , donde C es una matriz de Cauchy: 64

Teor´ıa de perturbaciones multiplicativas y soluciones precisas del LSP LS-QR: dados los vectores z e y, las componentes de C son calculadas como en (3.6.4) y el problema de m´ınimos cuadrados es resuelto usando la factorizaci´on QR de Householder implementada en MATLAB—. LS-RRD: el problema de m´ınimos cuadrados es resuelto usando el Algoritmo 3.5.1 y la RRD del Paso 1 es calculada con la versi´on rectangular del Algoritmo 3 en [18] discutido anteriormente. Las factorizaciones QR necesarias en los Pasos 2 y 4 del Algoritmo 3.5.1 son calculadas con la rutina qr de MATLAB—. Observemos que en este caso, en el sistema lineal Y x = x2 del Paso 4 la matriz Y es no singular y por lo tanto GEPP se puede utilizar para obtener su soluci´on. En nuestros experimentos, hemos generado matrices de Cauchy con vectores z e y aleatorios, tambi´en hemos generado vectores b y hemos calculado la soluci´on de m´ınx∈Rn kCx − bk2 usando los algoritmos LS-QR y LS-RRD. Para calcular los errores relativos kb x0 −x0 k2 /kx0 k2 , tomamos como soluci´on “exacta” x0 la calculada utilizando el comando svd de MATLAB— ejecutado en precisi´on aritm´etica variable. En cada experimento hemos fijado la precisi´on a 2 log10 |D1 /Dn | + 30 d´ıgitos decimales, donde D1 y Dn son respectivamente, las entradas diagonales m´as grande y m´as peque˜ na (en valor absoluto) de la matriz diagonal D en la RRD de C calculada en el Paso 1 del Algoritmo 3.5.1. La motivaci´on de tomar 2 log10 |D1 /Dn |+30 d´ıgitos decimales, viene del hecho de que |D1 /Dn | tiene una magnitud similar a κ2 (C), ya que X e Y est´an bien condicionadas, y esto, de acuerdo con (3.6.1) y a la discusi´on en la Secci´on 3.3.1, el error en el algoritmo tradicional para m´ınimos cuadrados es casi siempre mucho menor que u κ22 (C). Los vectores aleatorios z, y y b se han elegido de la distribuci´on uniforme en el intervalo [0, 1] (comando rand en MATLAB—), o bien, de la distribuci´on normal (comando randn en MATLAB—). En todos los experimentos hemos analizado las ocho posibilidades resultantes en la selecci´on de las distribuciones aleatorias para z, y y b. Dos tipos de experimentos se han realizado. En el primer grupo hemos fijado el tama˜ no de la matriz: m × n = 100 × 50, 50 × 30 ´o 25 × 10. Para cada tama˜ no hemos generado 50 × 8 conjuntos diferentes de vectores aleatorios z, y y b, por lo tanto, generamos un total de 400 problemas de m´ınimos cuadrados diferentes para cada tama˜ no. En la Figura 3.6.1 mostramos los resultados para el caso 100 × 50 cuando los vectores z y b son seleccionados de la distribuci´on normal y el vector y de la distribuci´on uniforme en [0, 1]. Hemos graficado en una escala log-log el error relativo kb x0 −x0 k2 /kx0 k2 de la soluci´on, contra el n´ umero de condici´on de las matrices 65

Cap´ıtulo 3 (calculado a partir de la SVD “exacta” en precisi´on aritm´etica variable usada para calcular la soluci´on “exacta” x0 ) para los algoritmos LS-QR y LS-RRD. Aparte de esto, tambi´en hemos representado las cantidades ΘQR y ΘRRD que aparecen en la ecuaci´on (3.6.3). Observamos que el error relativo en el algoritmo LS-RRD es de orden u veces una peque˜ na constante, como se predijo, mientras que el error para LS-QR es escalado casi lineal con respecto a κ2 (C) hasta saturarse. La dependencia lineal en κ2 (C) del error relativo en LS-QR es la predecida en la ecuaci´on (3.6.1) dado que kC † k2 kbk2 /kx0 k2 ha sido siempre moderada en estos experimentos. Tambi´en podemos observar que la cota ΘRRD es bastante buena y no sobreestima los errores actuales. Para otros tama˜ nos y otras formas de generar z, y y b, los resultados han sido similares. En nuestro segundo grupo de experimentos, hemos fijado el n´ umero de filas de la matriz y variado el n´ umero de las columnas. Hemos probado con matrices de tama˜ no m = 100, n = 10 : 10 : 90 (5 × 8 conjuntos de vectores aleatorios z, y y b para cada tama˜ no), m = 50, n = 10 : 2 : 40 (10 × 8 conjuntos de vectores aleatorios z, y y b para cada tama˜ no) y m = 25, n = 5 : 5 : 20 (20 × 8 conjuntos de vectores aleatorios z, y y b para cada tama˜ no). Esto hace un total de 2280 matrices. La Figura 3.6.2 muestra los resultados para m = 50, n = 10 : 2 : 40 para cuatro combinaciones diferentes de distribuciones aleatorias para z, y y b. Para cada tama˜ no graficamos el m´aximo error relativo de las 10 muestras. De nuevo, los errores relativos de la soluci´on del algoritmo LS-RRD son de orden u veces una constante peque˜ na, mientras que para el algoritmo LS-QR son muy grandes. Para otros tama˜ nos y otras formas de generar z, y y b se obtienen resultados similares. Para todos los experimentos con matrices de Cauchy, el rango del n´ umero de condici´on ha sido 100 . κ2 (C) . 10100 , el m´aximo valor del t´ermino kC † k2 kbk2 /kx0 k2 ha sido 1376, 8 ≤ κ2 (X) ≤ 72 y 13 ≤ κ2 (Y ) ≤ 58.

3.6.2.

Matrices de Vandermonde

Hemos desarrollado experimentos num´ericos del Algoritmo 3.5.1 similares a los presentados en la Secci´on 3.6.1 con matrices de Vandermonde. Las matrices de Vandermonde aparecen naturalmente en problemas de interpolaci´on polinomial. Dado 66

Teor´ıa de perturbaciones multiplicativas y soluciones precisas del LSP

[] LS−QR LS−RRD

10

10

+ .............. + +

ΘRRD ΘQR

",+

+

+ 5

10

k∆xk 2 0 kx 0 k 2 10

::: O

−5

10

o

+ O

+-I:t-+ o

*o

~©o

o

a bD

&60-0 -0JlJ-

O-Ql'i[l

-O

000- -O

-«!ID-O

~o

O O

o

−10

10

v

V '7 'lOO O O

−15

10

15

10

20

10

Vv

!ilJ o

25

10

κ2 (C) Figura 3.6.1: Error relativo progresivo kb x0 − x0 k2 /kx0 k2 frente κ2 (C). C matrices de Cauchy aleatorias de orden 100 × 50. Los vectores z y b son seleccionados de la distribuci´ on normal est´ andar y el vector y de la distribuci´ on uniforme en [0, 1].

un vector de puntos z = [z1 , . . . , zm ]T ∈ Rm tal que zi 6= zj si i 6= j y un conjunto de “valores de la funci´on” b = [b1 , . . . , bm ]T ∈ Rm , el problema de encontrar el polinomio Pn−1 (z), de grado menor o igual a n − 1, con n ≤ m, que mejor se ajusta a los datos z y b en el sentido de m´ınimos cuadrados es equivalente a resolver el problema de m´ınimos cuadrados m´ınc∈Rn kV c − bk2 donde V ∈ Rm×n es una matriz de Vandermonde, cuyas entradas est´an dadas por vij = zij−1 ,

i = 1, . . . , m, j = 1, . . . , n,

(3.6.5)

y la soluci´on buscada c ∈ Rn es el vector que contiene los coeficientes del polinomio Pn−1 (z) = c1 + c2 z + · · · + cn z n−1 . Un m´etodo que calcula una RRD precisa de cualquier matriz de Vandermonde fue presentado en [18, Secci´on 5] y est´a basado en el siguiente hecho, si F ∈ Rn×n es la transformada discreta de Fourier de orden n×n, entonces V F es una matriz quasi-Cauchy cuyos par´ametros pueden ser calculados de forma precisa en O(mn) operaciones, as´ı como las sumas y restas de cualquier par de estos par´ametros. Entonces, una RRD precisa de V F = XDY puede ser calculada con el Algoritmo 3 en [18] adaptado al caso de matrices rectangulares como en la Secci´on 3.6.1. Finalmente, V = X D (Y F T ) es una RRD precisa de V . El coste total 67

Cap´ıtulo 3 mode[z,y,b] = [1, 2,2] 0

mode[z,y,b] = [2, 1,2]

o

10

0

0-:-0--0

10

o O

O

o o o

−10

10

o

o

10

−10

10

o

D D D D 15

OD 20

LS−QR LS−RRD

LS−QR LS−RRD

D

D

CJ

D

D D 25

D D

30

D

35

40

10

15

mode[z,y,b] = [2, 2,1]

o

0

10

O D

25

30

35

40

mode[z,y,b] = [1, 1,1] 0

O O -O 0--0-:-0 O 0- -O- -0-: -0- -O

O

20

10

LS−QR LS−RRD

O D

LS−QR LS−RRD O

O O

−10

O O O

−10

10

10

O O O O

D D

10

D D

15

D

D

20

D D

D D

25

D D

30

O

D

O

D

35

D

40

10

tl

D D D D

15

O D

20

D D D D D D D

25

30

35

40

Figura 3.6.2: Error relativo progresivo kb x0 − x0 k2 /kx0 k2 frente n, para matrices de Cauchy de orden m × n con m = 50 y n = 10 : 2 : 40 (10 matrices por cada tama˜ no) para cuatro combinaciones diferentes de distribuciones aleatorias (mode=1 denota la distribuci´ on normal est´ andar y mode=2 denota la distribuci´ on uniforme en [0, 1]) para z, y y b.

es el mismo que para matrices de Cauchy: 2mn2 − 2n3 /3 + O(n2 + mn) operaciones m´as mn2 /2 − n3 /6 + O(n2 + mn) comparaciones. Este algoritmo ser´a usado para calcular la RRD en el Paso 1 del Algoritmo 3.5.1. Los dos algoritmos usados en esta secci´on para resolver m´ınx∈Rn kV x − bk2 son:

• LS-QR: dado el vector z, las componentes de V son calculadas como en (3.6.5) y el problema de m´ınimos cuadrados es resuelto usando la factorizaci´on QR de Householder implementada en MATLAB—.

• LS-RRD: el problema de m´ınimos cuadrados es resuelto usando el Algoritmo 3.5.1 y la RRD del Paso 1 es calculada con la versi´on rectangular discutida anteriormente del m´etodo presentado en [18, Secci´on 5]. En nuestros experimentos, generamos vectores aleatorios b y matrices de Vandermonde con vectores aleatorios z, y calculamos la soluci´on de m´ınx∈Rn kV x−bk2 usando los algoritmos LS-QR y LS-RRD. Para calcular el error relativo kb x0 − x0 k2 /kx0 k2 , seguimos el procedimiento presentado en la Secci´on 3.6.1 para obtener la soluci´on 68

Teor´ıa de perturbaciones multiplicativas y soluciones precisas del LSP “exacta” x0 . Los vectores aleatorios z y b han sido elegidos de la distribuci´on uniforme en [0, 1] o bien, de la distribuci´on normal. En todos los experimentos hemos probado las cuatro combinaciones resultantes de la elecci´on de las distribuciones aleatorias para z y b. Hemos probado con matrices m × n de tama˜ nos m = 50, n = 5 : 5 : 30; m = 100, n = 10 : 5 : 60; y m = 500, n = 100 : 50 : 250. Para cada tama˜ no hemos generado diferentes conjuntos de vectores aleatorios z y b (para m = 50, 100, 25 × 4 diferentes conjuntos y para m = 500, 10 × 4 diferentes conjuntos), generando un total de 1860 problemas de m´ınimos cuadrados diferentes. La Figura 3.6.3 muestra los resultados para los tama˜ nos m = 100 y n = 10 : 5 : 60, cuando los vectores z y b son seleccionados de la distribuci´on normal. Hemos graficado en una escala log-log el error relativo kb x0 −x0 k2 /kx0 k2 de la soluci´on frente el n´ umero de condici´on de las matrices (calculado a partir de la SVD “exacta” como en la Secci´on 3.6.1) para los algoritmos LS-QR y LS-RRD. Adem´as, hemos graficado las cantidades ΘQR y ΘRRD 4 que aparecen en la ecuaci´on (3.6.3). Los resultados son similares a los obtenidos en la Figura 3.6.1 para las matrices de Cauchy y referimos al lector a los comentarios hechos en la Secci´on 3.6.1. Para otros tama˜ nos y otras formas de generar z y b los resultados han sido similares, produciendo el algoritmo LS-RRD errores relativos que son siempre de orden u veces una peque˜ na constante y el algoritmo LS-QR errores relativos escalados linealmente con respecto a κ2 (V ) y que son muy grandes para κ2 (V ) muy grandes. Para todos nuestros experimentos con matrices de Vandermonde, el rango de los n´ umeros de condici´on ha sido 100 . κ2 (V ) . 1070 , el m´aximo valor del t´ermino kV † k2 kbk2 /kx0 k2 ha sido 1076, 4 ≤ κ2 (X) ≤ 65 y 3 ≤ κ2 (Y ) ≤ 87.

3.6.3.

Matrices Graduadas

Otra clase de matrices para la cual es posible calcular una RRD precisa bajo ciertas condiciones son las matrices graduadas: matrices de la forma A = S1 BS2 ∈ Rm×n , con S1 ∈ Rm×m y S2 ∈ Rn×n matrices diagonales y no singulares que pueden ser arbitrariamente mal condicionadas, B ∈ Rm×n una matriz bien condicionada, y rango(A) = n. Por lo tanto, la matriz A puede tener n´ umero de condici´on muy grande. Higham en [42] determina condiciones tales que si la factorizaci´on QR con 4

Para mantener la escala de la figura graficamos m´ın (ΘQR , 10) en vez de ΘQR , dado que valores de ΘQR mucho m´ as grandes que los que aparecen en la Figura 3.6.1 han aparecido en este experimento.

69

Cap´ıtulo 3 2

10

0

10

−2

10

LS−QR LS−RRD ΘRRD

−4

10

ΘQR

−6

k∆xk 2 10 kx 0 k 2

−8

10

−10

10

−12

10

−14

10

−16

10

0

10

10

10

20

10

30

10

40

10

50

10

60

10

κ2 (V ) Figura 3.6.3: Error relativo progresivo kb x0 − x0 k2 /kx0 k2 frente κ2 (V ) para matrices de Vandermonde aleatorias V de tama˜ nos 100 × 10 : 5 : 60. Los vectores aleatorios z y b son seleccionados de la distribuci´ on normal est´ andar.

pivote completo (pivote por columnas junto con pivote por filas u ordenamiento por filas, para m´as detalles ver [42] ) de una matriz graduada A es calculada y la SVD del factor R permutado es calculado por medio del Algoritmo 2.5.1, entonces la SVD de A se obtiene con alta precisi´on relativa. En esta secci´on mostramos que, bajo las mismas condiciones, la factorizaci´on QR con pivote completo puede usarse para resolver de manera precisa y eficiente el problema de m´ınimos cuadrados cuya matriz de coeficientes es graduada. Para este prop´osito, observemos que si PR APC = QR es la factorizaci´on QR reducida con pivote completo, donde PR y PC son matrices de permutaciones, Q ∈ Rm×n y R ∈ Rn×n , entonces una RRD A = XDY se puede obtener como se muestra a continuaci´on: A = PRT QRPCT = (PRT Q) D (D−1 R PCT ), donde D := diag(R), X := PRT Q e Y := D−1 R PCT . En este caso, el Paso 2 del Algoritmo 3.5.1 se reduce a x1 = QT PR b y los Pasos 3-4 pueden combinarse en uno s´olo sin afectar los errores de redondeo, dado que los errrores son componente a componente y D es diagonal. Combinar estos pasos se hace simplemente para resolver el sistema lineal (RPCT )x = x1 . Por lo tanto, D no es necesario y el Algoritmo 3.5.1 se simplifica al Algoritmo 3.6.1, el cual es casi siempre el algoritmo usual para resolver 70

Teor´ıa de perturbaciones multiplicativas y soluciones precisas del LSP problemas de m´ınimos cuadrados usando la factorizaci´on QR.

Algoritmo 3.6.1 (Soluci´on precisa para problemas de m´ınimos cuadrados v´ıa QR con pivote completo cuando la matriz de coeficientes es graduada) Input: A ∈ Rm×n matriz graduada, b ∈ Rm , m ≥ n = rango(A). Output: x0 soluci´on u ´nica de m´ınx∈Rn kb − A xk2 Paso 1: Calcular la descomposici´on QR reducida con pivote completo (ver [42]) de A, A = PRT QRPCT . Paso 2: Calcular x1 = QT PR b. Paso 3: Resolver el sistema lineal consistente R PCT x = x1 para obtener x0 .

El coste del Algoritmo 3.6.1 es esencialmente el mismo que el del m´etodo QR de Householder, es decir, 2mn2 − 23 n3 flops, ya que el coste del pivotaje es O(mn). Como es usual, no es necesario formar expl´ıcitamente la matriz Q. El Algoritmo 3.6.1 se puede aplicar tambi´en para matrices A que no sean graduadas, pero entonces la precisi´on no estar´ıa garantizada. Es importante notar que GECP tambi´en puede usarse para calcular bajo las mismas condiciones una RRD precisa de una matriz graduada [20, Secci´on 4] y el Algoritmo 3.5.1 puede usarse para resolver de manera precisa problemas de m´ınimos cuadrados cuando la matriz de coeficientes es graduada. Sin embargo, este procedimiento requiere calcular una factorizaci´on QR en el Paso 2 del Algoritmo 3.5.1 y por lo tanto es m´as costoso que el Algoritmo 3.6.1. Ahora, explicamos brevemente y de manera simplificada cuales son los errores cometidos por el Algoritmo 3.6.1. Estos errores son b´asicamente iguales a los del Algoritmo 3.5.1, si consideramos como RRD de A la siguiente factorizaci´on A = (PRT Q) D (D−1 R PCT ), y por lo tanto son determinados por los errores cometidos en el c´alculo de la factorizaci´on QR con pivote completo. Para hacer la notaci´on m´as simple, asumimos que la matriz A ha sido pre-pivoteada, es decir, que el Paso 1 del Algoritmo 3.6.1 produce PR = Im y PC = In . Observe que esto induce correspondientemente pre-permutaciones en S1 , B y S2 . Primero, se demostr´o en [42, Teorema 2.5] que si A = S1 BS2 ∈ Rm×n (m ≥ n) con S1 y S2 matrices diagonales arbitrarias 71

Cap´ıtulo 3 b calculado en el Paso 1 del Algoritmo 3.6.1, y no singulares, entonces el factor R satisface b S1 (B + ∆B)S2 = QR,

con k∆Bk2 = O(u)kBk2 ,

(3.6.6)

donde Q ∈ Rm×n es una matriz con columnas ortonormales. La notaci´on O-grande en (3.6.6) oculta algunos factores: factores de crecimiento y polinomios de grado menor en m y n, que pueden ser importantes en algunos casos. Para m´as detalles, ver [42]. Sin embargo, nuestros experimentos (como aquellos en [42]) muestran que en la pr´actica O(u) es una constante peque˜ na multiplicada por u. Si calcul´aramos el factor b tal que kQ − Qk b 2 = O(u) Q expl´ıcitamente, entonces obtendr´ıamos una matriz Q [43, ecuaci´on (19.13)], donde Q es la matriz que aparece en (3.6.6). Por lo tanto, en la siguiente discusi´on, no distinguiremos entre la Q exacta y la calculada. Ahora necesitamos obtener, del error regresivo en (3.6.6), errores progresivos sobre la RRD del mismo tipo que aparece en la Definici´on 2.4.2. Con este fin, recordemos que si B = LU es una factorizaci´on LU de B (sin pivote) con L ∈ Rm×n y U ∈ Rn×n , cuyos factores est´an bien condicionados y tales que kB † k2 ≈ kL† k2 kU −1 k2 , entonces en [20, Teorema 4.1] se demostr´o que S1 (B + ∆B)S2 puede escribirse como5 S1 (B + ∆B)S2 = (I + E)A(I + F ),

(3.6.7)

con m´ax(kEk2 , kF k2 ) = O(τ k∆Bk2 kB † k2 ) = O(u) τ κ2 (B),

(3.6.8)

donde el factor τ controla el gradiente (despu´es de las permutacion en el Paso 1 del Algoritmo 3.6.1) y est´a dado por τ = m´ax(1, τ1 , τ2 ),

|(S1 )kk | |(S2 )kk | , τ2 = m´ax . (3.6.9) 1≤j≤n 1≤j≤k≤n |(S2 ) | |(S ) | 1 jj jj j ≤k≤m

con τ1 = m´ax

Combinando (3.6.6) y (3.6.7), ignorando los t´erminos de segundo orden y definiendo b := diag(R), b obtenemos que A = (I −E)QD( b D b −1 R)(I b −F ). Por lo tanto, X = (I − D b e Y = (D b −1 R)(I b −F ) pueden ser consideradas como los factores de una RRD E)Q, D b = Q, D b e Yb = (D b −1 R). b Observemos exacta de A, cuyos factores calculados ser´ıan X que en este caso los factores diagonales exactos y calculados son los mismos. Por lo b − Xk2 /kXk b 2 , kYb − Y k2 /kYb k2 } = tanto, podemos usar (3.6.8) para obtener m´ax{kX O(u) τ κ2 (B), lo cual nos permite aplicar el Teorema 3.5.2-(c) reemplazando p(m, n) 5

Las actuales cotas son m´ as complicadas e incluyen κ2 (L) y κ2 (U ), y pueden ser mucho m´as grandes que κ2 (B), ya que el pivotaje se realiza sobre A y no sobre B (para m´as detalle ver [20, Secci´on 4] y [42]). En nuestra discusi´ on, pretendemos enfatizar los principales factores que controlan los errores en la pr´ actica y no presentar un an´alisis riguroso.

72

Teor´ıa de perturbaciones multiplicativas y soluciones precisas del LSP por τ κ2 (B) y obtener a primer orden el siguiente error progresivo en la soluci´on calculada por el Algoritmo 3.6.1   kb x0 − x0 k2 kA† k2 kbk2 −1 b b ≤ O(u) τ κ2 (B) κ2 (D R) + + O(u2 ) , (3.6.10) kx0 k2 kx0 k2 donde usamos el hecho que κ2 (Q) = 1. Un punto clave en la cota (3.6.10) es el par´ametro τ , el cual penaliza a κ2 (B). Puede ser muy grande ya que las matrices diagonales S1 y S2 pueden ser arbitrariamente mal condicionadas. Sin embargo, las permutaciones provenientes de QR con pivote completo casi siempre reordenan S1 y S2 de forma que τ sea de orden uno. 0

10

−2

10

−4

10

k∆xk 2 10−6 kx0 k 2

−8

10

m*u*κ2(B)

−10

10

LS−QR LS−RRD

−12

10

2

10

4

10

6

10

8

10

10

10

κ2 (B) Figura 3.6.4: Error relativo progresivo kb x0 −x0 k2 /kx0 k2 frente κ2 (B) para matrices graduadas aleatorias A = S1 BS2 de tama˜ no 100 × 40, con B, S1 y S2 generados con la opci´ on (b) que se explica en el texto.

Con el objetivo de comprobar la precisi´on del Algoritmo 3.6.1 para matrices graduadas A, hemos desarrollado diferentes experimentos num´ericos similares a los presentados en [42] y hemos utilizado los siguientes dos Algoritmos para resolver m´ınx∈Rn kAx − bk2 : 73

Cap´ıtulo 3 LS-QR: calcula la soluci´on usando la factorizaci´on QR de Householder con pivote por columnas implementada en MATLAB—. LS-RRD: usa el Algoritmo 3.6.1 y la factorizaci´on QR necesaria en el Paso 1 se obtiene usando ordenamiento por filas y pivotaje por columnas [42]. Observe que los algoritmos LS-QR en las Secciones 3.6.1 y 3.6.2 usaron la factorizaci´on QR de Householder sin pivote. Ahora, usaremos pivote por columnas para ilustrar que el ordenamiento por filas adicional en LS-RRD es fundamental para obtener soluciones precisas. Para calcular el error relativo kb x0 − x0 k2 /kx0 k2 , seguimos el procedimiento presentado en la Secci´on 3.6.1 con D = diag(R). En nuestros experimentos, hemos generado matrices aleatorias de la forma A = S1 BS2 , donde B es construida siempre usando el modo 3 en la rutina randsvd de “Test Matrix Toolbox” desarrollado por Higham [41], es decir, B es una matriz densa aleatoria con un n´ umero de condici´on dado y con sus valores singulares distribuidos geom´etricamante. Las matrices diagonales S1 y S2 son generadas tambi´en con randsvd usando una variedad de distribuciones para valores singulares, es decir, para sus entradas diagonales. Hemos llevado a cabo experimentos donde los tama˜ nos m×n de las matrices A y B han sido m = 50, n = 10 : 10 : 30 y m = 100, n = 20 : 20 : 60. Las matrices A se han construido de la siguiente manera: hemos seleccionado matrices B con n´ umeros de condici´on κ2 (B) = 10i , para i = 1 : 10. Hemos generado matrices diagonales S1 y S2 con entradas diagonales seleccionadas de uno de los tres pares de distribuciones: (a) las entradas de S1 y S2 tienen sus logar´ıtmos uniformemente distribuidos (modo 5 en randsvd), pero las entradas de S1 est´an ordenadas decrecientemente y las entradas de S2 est´an ordenadas crecientemente; (b) las entradas de S1 y S2 est´an geom´etricamente distribuidas (modo 3 en randsvd) pero las entradas de S1 est´an ordenadas crecientemente y las entradas de S2 est´an ordenadas decrecientemente; (c) entradas distribuidas geom´etricamente en orden decreciente para S1 y entradas con sus logar´ıtmos uniformemente distribuidos ordenadas crecientemente para S2 . Tomamos κ2 (S1 ) = κ2 (S2 ) = 10k con k = 2 : 2 : 16. Para cada tama˜ no, para cada opci´on (a), (b) o´ (c), y para cada tripleta (κ2 (B), κ2 (S1 ), κ2 (S2 )), fueron generadas diez matrices y diez vectores b, cuyas entradas siguen la distribuci´on normal. Esto hace, para cada tama˜ no y cada κ2 (B), un total de 10 × 8 × 3 = 240 matrices. En la Figura 3.6.4 mostramos los resultados para matrices S1 y S2 de tama˜ no 100 × 40 con opci´on (b) para las matrices. Hemos graficado en una escala log-log el m´aximo error relativo para todas las 80 matrices con un κ2 (B) en particular. Observamos que el error para el algoritmo LS-RRD se comporta como O(u) κ2 (B), el cual, 74

Teor´ıa de perturbaciones multiplicativas y soluciones precisas del LSP b −1 R) b y kA† k2 kbk2 /kx0 k2 han de acuerdo con (3.6.10) implica que los n´ umeros τ , κ2 (D sido de orden uno en todos estos experimentos. El error del algoritmo LS-QR pierde hasta seis d´ıgitos de precisi´on, comport´andose notablemente peor que LS-RRD. Hemos graficado una l´ınea punteada mostrando la cantidad m u κ2 (B). El comportamiento para todos los otros tama˜ nos de matrices y todos los otros modos de randsvd son similares. Para todos nuestros experimentos con matrices graduadas el rango de los n´ ume40 ros de condici´on ha sido 10 . κ2 (A) . 10 y el valor m´aximo del t´ermino b −1 R) b ≤ 18. kA† k2 kbk2 /kx0 k2 ha sido 108 y 2 ≤ κ2 (D

3.6.4.

Experimentos num´ ericos controlando el residuo

En los experimentos num´ericos realizados anteriormente el vector b ha sido seleccionado aleatoriamente. Por lo tanto, el residuo relativo ρr := kb − Ax0 k2 /kbk2 ha sido pocas veces muy peque˜ no, ya que si rango(A) = n < m, entonces es muy poco probable que θ(b, R(A)) sea muy peque˜ no. En esta secci´on, consideramos diferentes tipos de experimentos en los cuales generamos vectores aleatorios b con un valor fijo de ρr . Para este prop´osito, hemos desarrollado experimentos con matrices de Cauchy y Vandermonde de tama˜ no m × n de la siguiente manera. Primero obtenemos la RRD de la matriz A = XDY . Calculamos la SVD completa de la matriz X (lo cual se puede hacer de manera precisa utilizando el comando svd de MATLAB— dado que X est´a bien condicionada): X = UX ΣX VXT , donde UX ∈ Rm×m , ΣX ∈ Rm×n y VX ∈ Rn×n . Entonces, particionamos UX = [U1 U2 ], donde U1 ∈ Rm×n y U2 ∈ Rm×(m−n) , y generamos vectores aleatorios α ∈ Rn y β ∈ Rm−n tales que kαk2 = kβk2 = 1. Con esto definimos b0 := U1 α ∈ R(A),

∆b := t U2 β ∈ R(A)⊥

y b := b0 + ∆b ,

(3.6.11)

donde hemos usado que R(A) = R(X) y t ≥ 0 es un par´ametro. Observemos que de esta manera ρr =

kb − Ax0 k2 t =√ , kbk2 1 + t2

(3.6.12)

ya que la soluci´on x0 satisface Ax0 = b0 . Hemos usado matrices de Cauchy de tama˜ no m = 100 × n = 20 : 20 : 60 y matrices de Vandermonde de tama˜ no m = 50 × n = 5 : 5 : 25. Para cada tama˜ no hemos cambiado el valor de t para [−16:2:−2] obtener residuos relativos ρr = 10 . Para cada tama˜ no y para cada valor de ρr hemos generado 10 matrices y 10 vectores b, usando la distribuci´on normal para 75

Cap´ıtulo 3 los par´ametros de las matrices y para los vectores α y β. Finalmente, para cada valor de ρr , obtenemos el valor m´aximo de todos los errores relativos progresivos para todos los tama˜ nos y todos los experimentos aleatorios. Los resultados para las matrices de Vandermonde son mostrados en la Tabla 3.6.1. Resultados similares fueron obtenidos para las matrices de Cauchy. Puede observarse que el an´alisis en la Secci´on 3.3.1 se cumple independientemente del tama˜ no del residuo relativo, a´ un cuando este es muy peque˜ no. Resaltamos nuevamente que el verdadero punto importante en la cota (3.0.2) es que kA† k2 kbk2 /kx0 k2 es peque˜ no para casi todos los vectores b para cualquier tama˜ no fijo del residuo relativo no cercano a uno, independientemente del mal condicionamiento de la matriz A. log10 (kb − Ax0 k2 /kbk2 )

−16

−14

−12

−10

−8

−6

−4

−2

QR : log10 (k∆xk2 /kx0 k2 )

−2,7

−2,8

−2,1

−2,5

−3,8

−3,5

−3,3

−2,6

RRD : log10 (k∆xk2 /kx0 k2 )

−14,1

−14,0

−13,8

−13,9

−14,1

−14,0

−13,8

−13,8

Tabla 3.6.1: Experimentos controlando el residuo. El error relativo progresivo k∆xk2 /kx0 k2 se muestra para ambos algoritmos LS-QR y LS-RRD descritos en la Secci´ on 3.6.2 para diferentes valores del residuo relativo. Este experimento se hizo con matrices de Vandermonde de tama˜ nos m = 50 × n = 5 : 5 : 25. Todos los vectores aleatorios necesarios se eligieron de la distribuci´ on normal est´ andar.

3.6.5.

Experimentos donde kA† k2 kbk2 / kx0 k2 no es peque˜ no

En todos los experimentos aleatorios presentados en las Secciones 3.6.1, 3.6.2, 3.6.3 y 3.6.4 el factor kA† k2 kbk2 /kx0 k2 ha sido moderado y por lo tanto el Algoritmo 3.5.1 resuelve de manera precisa todos los problemas ejecutados de m´ınimos cuadrados. Por lo tanto, estos experimentos han confirmado la discusi´on en la Secci´on 3.3.1. Por supuesto, es posible preparar experimentos donde kA† k2 kbk2 /kx0 k2 no sea peque˜ no y el Algoritmo 3.5.1 no sea preciso, pero esto requiere seleccionar cuidadosamente los vectores b. Hemos procedido de la siguiente manera: Definimos los vectores b como b = u1 + b⊥ , donde u1 es el vector singular izquierdo de A correspondiente al valor singular m´as grande y b⊥ es cualquier vector aleatorio ortogonal a p R(A). Notemos que para vectores de este tipo kA† k2 kbk2 /kx0 k2 = κ2 (A) 1 + kb⊥ k22 y como consecuencia, los errores relativos progresivos en las soluciones cometidos por el Algoritmo 3.5.1 han sido grandes y proporcionales a la unidad de redondeo por el n´ umero de condici´on de las matrices. Sin embargo, sin importar este hecho, los errores del Algoritmo 3.5.1 han sido mucho m´as peque˜ nos que los cometidos por el algoritmo QR de Householder. La raz´on es que el error (3.0.1) para la QR de Householder incluye el t´ermino Φ definido en (3.3.14) y para los vectores b = u1 + b⊥ , 76

Teor´ıa de perturbaciones multiplicativas y soluciones precisas del LSP Φ = κ2 (A)2 kb⊥ k2 , el cual es muy grande si κ2 (A) es grande y kb⊥ k2 = krk2 no es muy peque˜ no. Recordemos en este contexto nuestra discusi´on de la ecuaci´on (3.3.15). Como ha sido resaltado en [32], note que para matrices mal condicionadas, los vectores b = u1 + b⊥ deben ser generados usando algoritmos altamente precisos para obtener u1 y b⊥ (ver [20] y Secci´on 3.6.4). Si el vector b es construido usando aritm´etica en coma flotante usual y el comando svd en MATLAB—, los errores de redondeo hacen imposible que b tenga exactamente la estructura requerida y el Algoritmo 3.5.1 calcule soluciones con alta precisi´on relativa.

77

4 C´ alculo de autovalores y autovectores precisos de matrices sim´ etricas graduadas

Cuando utilizamos algoritmos convencionales para calcular autovalores y autovectores de matrices sim´etricas reales mal condicionadas en aritm´etica en coma flotante, s´olamente los autovalores con valor absoluto grande son calculados con precisi´on relativa garantizada. Los autovalores muy peque˜ nos pueden calcularse sin precisi´on relativa, e incluso con el signo equivocado. Los autovectores est´an calculados con errores peque˜ nos con respecto al gap relativo de los autovalores. Esto significa que si u es la unidad de redondeo y, qi y qbi son los respectivos autovectores exactos y calculados correspondientes al autovalor λi , entonces el a´ngulo agudo entre estos vectores est´a acotado por sin θ(qi , qbi ) ≤ O(u)/gapi , donde gapi = m´ınj6=i |λi − λj |/ m´axk |λk |. Lo cual implica que si existe m´as de un autovalor muy peque˜ no, entonces los correspondientes autovectores est´an calculados con errores grandes, incluso si los autovalores muy peque˜ nos est´an bien separados en el sentido relatido. Ver [1, Secci´on 4.7] para un resumen sobre cotas de errores para problemas de autovalores sim´etricos. Esto puede ser muy diferente si consideramos matrices con cierta estructura. Una matriz sim´etrica A = AT ∈ Rn×n es llamada graduada (o escalada) si B := S −1 AS −1 es una matriz bien condicionada para alguna matriz diagonal S = diag[s1 , . . . , sn ], S ser´a el escalamiento de la matriz graduada A. En [54], Martin, Reinsch y Wilkinson demostraron que una matriz graduada tendr´a autovalores peque˜ nos insensibles a peque˜ nas perturbaciones relativas en los elementos de la matriz. En este cap´ıtulo estamos interesados en encontrar condiciones sobre las matrices S y B para obtener soluciones precisas del problema de autovalores cuando A es una matriz sim´etrica. Nuestro objetido es dise˜ nar un algoritmo para calcular autovalores y autovectores de matrices sim´etricas graduadas de orden n × n con alta precisi´on 79

Cap´ıtulo 4 relativa, respetando la simetr´ıa del problema, con coste O(n3 ), es decir, aproximadamente el mismo coste de los algoritmos convencionales para matrices sim´etricas. Por alta precisi´on relativa en un problema de autovalores, queremos decir que los bi y qbi , autovalores λi , los autovectores qi y las correspondientes partes calculadas λ satisfacen:

bi −λi | ≤ O(u)|λi | y |λ

sin θ(qi , qbi ) ≤

O(u) relgapλ (A, λi )

para i = 1, . . . , n. (4.0.1)

Estas condiciones garantizan que los nuevos algoritmos calculan todos los autovalores, incluyendo los m´as peque˜ nos con el signo y d´ıgitos principales correctos. Sin embargo, los autovectores correspondientes a autovalores muy peque˜ nos relativamente bien separados son calculados de manera precisa. En el caso de autovalores m´ ultiples, o un cl´ uster de autovalores muy cercanos en el sentido relativo, la cota dada anteriormente para sin θ(qi , qbi ) es muy grande o tiende a infinito. En este caso, entendemos por alta precisi´on relativa que los senos de los ´angulos can´onicos entre los subespacios invariantes perturbados y sin pertubar correspondientes al cl´ uster de autovalores acotados por O(u) sobre el gap relativo entre los autovalores en el cl´ uster y aquellos fuera de ´el [51]. Lo cual significa que el nuevo algoritmo calcular´a bases precisas de subespacios invariantes correspondientes a cl´ usters de autovalores bien separados en el sentido relativo del resto de los autovalores. Hasta ahora la posibilidad de calcular los autovalores y autovectores con alta precisi´on relativa de matrices sim´etricas bien escaladas ha estado restringida para algunas clases especiales de matrices. En [3] los autores demostraron que los autovalores y los autovectores de matrices sim´etricas diagonalmente dominantes bien escaladas pueden ser calculados de manera precisa. En [25] se estudiaron las matrices sim´etricas escaladas definidas positivas (ver tambi´en [55, 56]). Esto es un resultado importante ya que el problema de autovalores para matrices definidas positivas fue el primero en ser tratado con el objetivo de obtener alta precisi´on relativa utilizando una factorizaci´on de Cholesky y rotaciones de Jacobi. Existen resultados parciales para el caso indefinido. En [67, 70], la clase de matrices para las cuales es posible resolver el problema fue ampliado significativamente, aunque sus resultados se basaron en las propiedades del factor polar definido positivo de A, que no son f´aciles de obtener a partir de la matriz A. Por tanto, no existen resultados, f´aciles de comprobar, sobre las condiciones que S = diag[s1 , . . . , sn ] y B tienen que satisfacer para obtener soluciones precisas del problema de autovalores de matrices sim´etricas generales A = SBS. En este cap´ıtulo vamos a demostrar 80

C´alculo de autovalores y autovectores precisos de matrices sim´etricas graduadas que, contrario a lo convencional, y al caso definido positivo, no es suficiente que B est´e bien condicionada y que los elementos diagonales de la matriz escalada est´en ordenados decrecientemente, despu´es de la permutaci´on producida por el pivotaje. Demostramos en el Corolario 4.2.4 que bajo perturbaciones de la matriz B de tama˜ no kδBk2 ≤ kBk2 , el error relativo en los autovalores y en los autovectores, para i = 1, . . . , n, est´a dado por λ e − λ i i ≤  τ O (κ2 (B)) + O(2 ) λi | sin θ(qi , qei )| ≤  τ

(4.0.2)

O (κ2 (B)) + O(2 ) relgapλ (A, λi )

(4.0.3)

donde θ(qi , qei ) es el a´ngulo agudo entre qi y qei . El factor τ controla el escalamiento y es definido como el m´aximo de los siguientes tres t´erminos,

τ = max{1, τL , τD },

sk τL = m´ax j