Métodos iterativos para sistemas lineales de ecuaciones e inecuaciones basados en proyecciones sucesivas

   FacultaddeCiencias      Trabajo Fin de Grado Grado en Matemáticas Curso 2018-2019 Métodositerativosparasistemaslinealesdeecuacioneseinecuaciones basadosenproyeccionessucesivas   Autor:D.JuanGonzálezAntolín  Tutor:Dr.LuisMªAbiaLlera    ii Prefacio El trabajo de Fin de Grado que presentamos tiene como objetivo el an´alisis del algoritmo de Proyecciones Alternadas de Von Neumann y su aplicaci´on en distintos problemas del ´algebra lineal y la optimizaci´on. El algoritmo de proyecciones alternadas en su forma m´as general se formula en el ´ambito de un espacio de Hilbert, y requiere por tanto de algunos de los elementos b´asicos de la teor´ıa de estos espacios, que recopilamos sin prueba, en su mayor´ıa, en el cap´ıtulo 1. El n´ucleo del trabajo lo conforman los cap´ıtulos 2 y 3. En el cap´ıtulo 2 probamos el teorema de Von Neumann de las proyecciones alternadas entre dos subespacios cerrados y su extesnsion por Halperin, a la intersecci´on de un n´umero finito de subespacios (cerrados). Consideramos tambi´en la cuesti´on de la velocidad asint´otica de convergencia del m´etodo que requiere introducir nociones de ´angulos entre subespacios de un espacio de Hilbert. El cap´ıtulo 3 selecciona algunos algoritmos del ´algebra lineal num´erica que pueden reinterpretarse en t´erminos del m´etodo de proyecciones alternadas. Esta selecci´on no pretende ser exhaustiva, estando condicionada por la tipolog´ıa que se describe en algunos textos como algoritmos de acci´on por filas ([2, 3]). Los teoremas de convergencia de estos algoritmos est´an dispersos en diferentes publicaciones (cuando se apartan de las hip´otesis del teorema de las proyecciones alternadas) y hemos seleccionado algunos de ellos. En Valladolid, a 12 de julio de 2019 iii iv ´ Indice general 1. Preliminares: Espacios de Hilbert 1 1.1. Introducci´on............................ 1 1.2. Conceptos previos . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2.1. Productos internos . . . . . . . . . . . . . . . . . . . . 2 1.2.13. Ortogonalidad . . . . . . . . . . . . . . . . . . . . . . 4 1.3. Operadores de proyecci´on . . . . . . . . . . . . . . . . . . . . 5 2. Teoremas de Von Neumann y Halperin 9 2.1. Teorema de Von Neumann . . . . . . . . . . . . . . . . . . . . 9 2.2. Teorema de extensi´on de Halperin. . . . . . . . . . . . . . . . 11 2.3. Velocidad de convergencia. . . . . . . . . . . . . . . . . . . . . 14 2.3.1. ´ Angulo entre subespacios. . . . . . . . . . . . . . . . . 14 2.3.6. Velocidad de convergencia de MAP. . . . . . . . . . . 16 2.4. T´ecnicas de aceleraci´on. . . . . . . . . . . . . . . . . . . . . . 20 3. M´etodos de acci´on por filas 23 3.1. Introducci´on............................ 23 3.2. Algunos m´etodos de Acci´on por Filas. . . . . . . . . . . . . . 24 3.2.1. El M´etodo de Kaczmarz . . . . . . . . . . . . . . . . . 24 3.2.2. El algoritmo de Cimmino . . . . . . . . . . . . . . . . 25 3.2.3. El M´etodo de Relajaci´on de Agmon, Motzkin y Schoenberg(MAMS)....................... 26 3.2.4. El algoritmo de Hildreth . . . . . . . . . . . . . . . . . 27 3.2.5. El m´etodo de relajaci´on de Herman . . . . . . . . . . 28 3.3. El m´etodo de Kaczmarz y la inversa generalizada. . . . . . . 29 3.4. El m´etodo de reconstrucci´on algebraica (ART) . . . . . . . . 36 3.4.1. Formulaci´on matem´atica del problema . . . . . . . . . 38 3.4.2. T´ecnicas de reconstrucci´on algebraica (ART) . . . . . 40 3.5. Algunos experimentos num´ericos . . . . . . . . . . . . . . . . 41 3.6. Conclusiones ........................... 44 A. Programas de Matlab 49 v vi Capitulo 1 Preliminares: Espacios de Hilbert 1.1. Introducci´on El objetivo de esta trabajo es el estudio del m´etodo de proyecciones alternadas en subespacios y variedades afines y su relaci´on con diferentes algoritmos para el tratamiento de problemas como: 1. El problema de la resoluci´on de sistemas Ax =b, con Amatriz m×n real, y bno necesariamente en el rango de b(en cuyo caso consideraremos soluciones de m´ınima norma). 2. El problema de la factibilidad en sistemas de inecuaciones lineales. 3. Resoluci´on de problemas de optimizaci´on con restricciones lineales. 4. T´ecnicas de reconstrucci´on algebraica en el tratamiento de im´agenes. El m´etodo de proyecciones alternadas, al que nos referiremos gen´ericamente como MAP, en su versi´on m´as simple se debe a John Von Neumann que consider´o el problema de la proyecci´on de un punto en un espacio de Hilbert a la intersecci´on de dos subespacios cerrados del mismo. Posteriormente, este resultado fu´e generalizado por Halperin a la intersecci´on de un n´umero finito de subespacios cerrados, y por Dykstra a la intersecci´on de un n´umero finito de convexos cerrados no vac´ıos. La memoria s´olo considerar´a el m´etodo de proyecciones alternadas en subespacios cerrados en un espacio de Hilbert en las versiones de Von Neumann y de Halperin. Por ello en este primer cap´ıtulo reunimos resultados b´asicos de la teor´ıa de espacios de Hilbert [1, 7, 10], en su mayor´ıa sin prueba, por estar ´estos incluidos en los programas de asignaturas del Grado. B´asicamente nos limitaremos a probar s´olo que un operador T:H→H, Hespacio de Hilbert, es un operador de proyecci´on si y s´olo si es lineal, 1 acotado, autoadjunto e idempotente, y est´a caracterizado por el subespacio R(T), el rango de T. Dedicaremos el cap´ıtulo 2 a la prueba del teorema de Von Neumann y a su extesnsion por Halperin. En el contexto del teorema de Von Neumann, la velocidad asint´otica de convergencia est´a estrechamente relacionada con los ´angulos entre los subespacios sobre los que se alternan las proyecciones. Describimos esta vinculaci´on y algunas t´ecnicas generales de aceleraci´on de la convergencia. Finalmente el cap´ıtulo 3 selecciona algunos algoritmos que pueden reinterpretarse en t´erminos del m´etodo MAP en ´areas como la soluci´on de sistemas de ecuaciones lineales, soluci´on de sistemas de inecuaciones lineales, determinados o sobredeterminados, y algunos problemas de optimizaci´on con restricciones. Describimos tambi´en una aplicaci´on al tratamiento de im´agenes m´edicas, que se conoce gen´ericametne con el nombre de t´ecnicas de reconstrucci´on algebraica (ART). Un an´alisis experimental completo de la familia de algoritmos propuestos se escapa al alcance de este trabajo, sobre todo tras comprobar con algunos ejemplos que resultan poco competitivos con algoritmos directos (en el caso de sistemas lineales determinados, si el tama˜no del sistema es moderado), e incluso con otros m´etodos iterativos cl´asicos. No obstante tienen su inter´es cuando se considera la soluci´on de sistemas sobredeterminados, o indeterminados, sujetos a restricciones en las variables de la soluci´on. 1.2. Conceptos previos 1.2.1. Productos internos 1.2.2 Definici´on. Sea Eun espacio vectorial sobre K. Un producto interno sobre Ees una aplicaci´on h·,·i de E×Een Kque verifica las siguientes propiedades: 1. hx+y, zi=hx, zi+hy, zipara todos x, y, z ∈E. 2. hλx, yi=λhx, yipara todos x, y ∈E,λ∈K. 3. hx, yi=hy, xipara todos x, y ∈E. 4. hx, xi ≥ 0 para todo x∈Eyhx, xi= 0 si, y solo si, x= 0. Un espacio prehilbertiano o con producto interno es un par (E, h·,·i) donde Ees un espacio vectorial sobre Kyh·,·i es un producto interno sobre E. 1.2.3 Proposici´on. Sea Eun espacio vectorial y h·,·i un producto interno sobre E. Se verifican las siguientes propiedades: 1. hx, λyi=λhx, yipara todos x, y ∈E, λ ∈K 2 2. hx, y +zi=hx, yi+hx, zipara todos x, y, z ∈E. 3. hx, 0i=h0, yi= 0 para todos x, y ∈E. 4. Si hx, yi= 0 para todo y∈E, entonces x= 0. 1.2.4 Teorema. (Desigualdad de Cauchy-Schwarz) Sea Eun espacio prehilbertiano. Si x,y∈Ese verifica que |hx, yi|2≤ hx, xihy, yi La igualdad se da si, y s´olo si, los vectores xeyson linealmente dependientes. 1.2.5 Proposici´on. Sea Eun espacio prehilbertiano. La aplicaci´on k k de E en R, definida por kxk=hx, xi1/2, x ∈E. (1.1) es una norma sobre E. 1.2.6 Definici´on. Un espacio de Hilbert es un espacio prehilbertiano completo para la norma definida por el producto interno como en (1.1). 1.2.7 Lema. Sean Eun espacio prehilbertiano y x, y ∈E. Entonces: kx+yk1/2=kxk1/2+kyk1/2+ 2Re(hx, yi). 1.2.8 Proposici´on. Sea Eun espacio prehilbertiano y x, y ∈E. Se verifican las siguientes propiedades: 1. Identidad del paralelogramo: kx+yk+kx−yk= 2kxk2+ 2kyk2. 2. Identidades de polarizaci´on: a) Si Ees real, hx, yi=1 4(kx+yk2− kx−yk2). b) Si Ees complejo, Re(hx, yi) = 1 4(kx+yk2− kx−yk2) Im(hx, yi) = 1 4(kx+iyk2− kx−iyk2) 3 hTmx, Tnyi=hTm+n−δx, Tyi Donde δvale 1 ´o 0 dependiendo de si myntienen la misma paridad. Necesitamos probar que para todo x∈Hel l´ım n→∞ Tnxexiste, para ello veremos que la sucesi´on {Tnx}es de Cauchy: kTmx−Tnxk2=hTmx−Tnx, Tmx−Tnxi(2.1) =hTmx, Tmxi−hTmx, Tnxi−hTnx, Tmxi+hTnx, Tnxi (2.2) =hT2m−1x, xi+hT2n−1x, xi − 2hTm+n−δx, xi(2.3) =hT2m−1x, xi+hT2n−1x, xi − 2hT2k−1x, xi(2.4) en la ´ultima expresi´on hemos tenido en cuenta que como m+n−δes siempre impar lo podemos sustituir por 2k−1, para cierto k. Por otro lado se tiene que: hT2i−1x, xi=hTix, Tixi=kTixk2 y kTi+1xk2=hT2i+1x, xi Podemos escribir Ti+1 como PMTixoPNTixy teniedo en cuenta que la norma de una proyecci´on es siempre o cero o uno: kTi+1xk2⩽kTixk2 de lo que se sigue que: hT1x, xi⩾hT3x, xi⩾hT5x, xi⩾... ⩾0 Entonces existe el l´ım i→∞hT2i−1x, xiy de (2.4) se tiene que: l´ım m,n→∞ kTmx−Tnxk= 0 Finalmente podemos concluir que existe l´ım n→∞ Tnx. Lo denotaremos por x∗. El operador definido como Tx =x∗es lineal, continuo, univaluado, su dominio es H, y cumple que hTx, T yi=hTx, yi. Por lo tanto ser´a una proyecci´on de la forma PL. Veamos ahora que necesariamente L=M∩N. Si x∈M∩N, entonces PMx=PNx=x,Tnx=xyTx =x. Por lo tanto x∈Ly se tiene que (M∩N)⊂L. A continuaci´on, PMT2i=T2i+1 yPNT2i−1=T2i, haciendo tender ia infinito tenemos que: PMT=TyPNT=T. Para cada y∈Hsea Ty =x∈L. Entonces, PMx=PMTy =T y =x∈M, y PNx=PNT y = Ty =x∈N. Concluimos que L⊆(M∩N) y que L=M∩N. Finalmente intercambiando los papeles de PMyPNen el argumento previo tenemos que Σ2tiene limite T0=PM∩N. Necesariamente T=T0y la demostraci´on queda completada. 10 2.2. Teorema de extensi´on de Halperin. Consideramos ahora el caso en el que el n´umero de subespacios sea mayor que dos. Denotaremos por PMi(i= 1...r) el operador proyecci´on sobre el subespacio Mide un espacio de Hilbert. 2.2.1 Teorema. (Halperin [5]) Sean M1, M2, M3, ..., Mrsubespacios cerrados de un espacio de Hilbert H. Entonces: l´ım n→∞(PMrPMr−1...PM1)nx=P∩r 1Mix Realizaremos ahora un peque˜no inciso sobre notaci´on: T, Ti, P denotan como ya es usual operadores lineales y acotados. M+N={x+y:x∈M, y∈N}; donde MyNson subespacios cerrados. S0(T)≡ {x∈H:Tx = 0}; S1(T)≡ {x∈H:Tx =x}; TM ≡ {T x :x∈H}; K(T)≡supnkTnk; Tes un operador no expansivo si kTk ≤ 1. Decimos que Tes normal si T∗T=TT ∗.Todo operador autoadjunto es normal. Para la prueba del teorema de Halperin necesitaremos unos lemas previos. 2.2.2 Lema. Si Tes no expansivo entonces S1(T) = S1(T∗). Demostraci´on. Como Tes no expansivo se tiene que kTk=kT∗k,kT∗k ≤ 1 (T∗es no expansivo). kxk2=hx, xi=hTx, xi=hx, T ∗xi≤kxkkT∗xk ≤ kxk2; por lo tanto hx, T∗i=kxkkT∗xkykT∗xk=kxk. De donde se deduce que, kx−T∗xk2=kxk2− hx, T∗xi−hT∗x, xi+kT∗xk2= 0, lo que implica que T∗x=x. Intercambiando los papeles de Ty de T∗en el argumento anterior se sigue que S1(T) = S2(T∗), tal y como queriamos probar. 2.2.3 Lema. Si Tes no expansivo entonces R(I−T)es el complemento ortogonal de S1(T)yS1(T) + R(I−T) = H. Demostraci´on. Por el lema previo sabemos que S1(T) = S1(T∗). Por tanto (I−T)∗x= 0 es lo mismo que h(I−T)∗x, yi= 0 para todo y. 11 h(I−T)∗x, yi=hx, (I−T)yi= 0, Se sigue que (I−T)∗= 0 es equivalente a x⊥R(I−T) para todo x. Por lo tanto S1(T) = S1(T∗) = S0((I−T)∗) es el complemento ortogonal de R(I−T). De lo anterior tenemos que S1(T) + R(I−T) = H, tal y como quer´ıamos probar. 2.2.4 Lema. Tes no expansivo e idempotente si, y s´olo si, Tes una proyecci´on. Adem´as Tes la proyecci´on sobre S1(T). Demostraci´on. Si Tes una proyecci´on de las propiedades del operador proyecci´on sabemos que Tes no expansivo e idempotente. Supongamos ahora que Tes no expansivo e idempotente. Por el teorema (1.3.12) Tser´a una proyecci´on si se tiene que hTx, yi=hx, Tyi. Sabemos que H=S1(T) + R(I−T), por lo tanto para cada x∈Hse tiene que: x=xs+xrdonde xs∈S1(T) y xr∈R(I−T). En primer lugar tenemos que, T(xs) = xspor definici´on de S1(T). Por otro lado, xr= (I−T)ypara un cierto y∈H, T((I−T)y) = T(y−Ty) = Ty −Ty = 0. De donde se sigue que : Tx =T xs+Txr=xs. Finalmente, hTx, yi=hT(xs+xr), ys+yri=hxs, ys+yri=hxs, ysi+hxs, yri= hxs, ysi=hxr, ysi+hxs, ysi=hxr+xs, ysi=hx, Tyi luego hTx, yi=hx, Tyi,Tes autoadjunto y por lo tanto es una proyecci´on. Adem´as, T=PS1puesto que R(T) = S1. 2.2.5 Lema. Si para i= 1, ..., r se tiene que kTixk<kxkcuando Tix6=x. Sea T=T1...Trentonces kTxk<kxkcuando Tx 6=xyT x =xsi y solo si Tix=xpara todo i. Demostraci´on. Si Ti=xpara todo i, entonces Tx =T1...Trx=x. Por otro lado, si Tix6=xpara alg´un i, sea kel mayor ital que: kTxk=kT1...Tkxk ≤ kTkxk<kxk. Por lo tanto, Tx 6=x,kTxk ≤ kxk, y T x =xsi, y solo si Tix=xpara todo i. 2.2.6 Lema. Si para cada i= 1...r se tiene que: kx−Tixk2≤ki(kxk2− kTixk2) (2.5) Para alg´un ki∈(0,´ınf) y para todo x∈H. Adem´as, si T=T1...Trexiste alg´un k∈(0,∞)tal que para todo x∈H se tiene que: kx−Txk2≤k(kxk2− kT xk2) (2.6) 12 Demostraci´on. Tenemos que: kx−T1T2xk2≤(kx−T2xk+kT2x−T1T2xk)2 ≤[2 m´ax(kx−T2xk,kT2x−T2T1xk)]2 ≤4(kx−T2xk2+kT2x−T1T2xk2) ≤4 m´ax(k1, k2)(kxk2− kT2xk2+kT2xk2− kT1T2xk2) = 4 m´ax(k1, k2)(kxk2− kT1T2xk2), entonces T1T2verifican la propiedad (2.6). Por inducci´on concluimos que T=T1...Trcumple la misma propiedad, tal y como queriamos demostrar. Podemos observar que si se satisface (2.5) entonces: N X n=0 kTnx−Tn+1xk ≤ N X n=0 k(kTnxk2− kTn+1xk2) =k(kxk2− kTn+1xk2) ≤kkxk2 Por lo tanto kTnx−Tn+1xk → 0 cuando n→ ∞, de donde tenemos que (Tn−Tn+1)x→0 cuando n→ ∞ para todo x∈H. 2.2.7 Lema. Sea Tun operador que satisface kx−Txk2≤k(kxk2− kT xk2) para alg´un k∈(0,∞)si se tiene que Tx 6=xentonces kTxk<kxk. Demostraci´on. Si: Tx 6=x⇒x−T x 6= 0 ⇒ kx−Txk 6= 0 ⇒k(kxk2− kTxk2)>0 ⇒kkxk2> kkTxk2 ⇒ kxk2>kTxk2 ⇒ kxk>kTxk. Tal y como quer´ıamos demostrar. 2.2.8 Teorema. Sea Tun operador que satisface la propiedad (2.6) entonces la sucesi´on de operadores Tnconverge fuertemente a PS1. Demostraci´on. Se deduce de forma directa de los lemas (2.2.3), (2.2.4) y (2.2.7). 13 A continuaci´on expondremos un corolario que es consecuencia directa de los lemas (2.2.5), (2.2.6) y del teorema (2.2.8). 2.2.9 Corolario. Si Ticon i= 1...r cumple la propiedad (2.5) y T=T1...Tr, entonces cuando n−→ ∞: Tn−→ P≡P∩r 1S1(Ti) Adem´as P x =xsi y solo si Tix=xpara alg´un i. Para finalizar esta secci´on, y la demostraci´on del teorema de extensi´on de Halperin, veremos a continuaci´on que todas las proyecciones cumplen la propiedad (2.6) para k= 1 : Sea Tun operador proyecci´on, para todo x∈Hse sigue que: kx−Txk2=hx−T x, x −Txi =kxk−hTx, xi−hx, T xi+kT xk2 Como Tes autoadjunto e idempotente: hTx, xi=hT2x, xi=hT x, Txi=kTxk2 por lo tanto: kx−Txk2=kxk2− kT xk2 Es f´acil ver que el teorema de extensi´on de Halperin ha quedado ya demostrado por ser un caso particular del corolario anterior. 2.3. Velocidad de convergencia. La velocidad de convergencia de MAP depender´a de los ´angulos entre los distintos subespacios involucrados. Este concepto merece nuestra atenci´on. Veamos en primer lugar que si xey∈Hen ´angulo θentre xeyse define como el ´angulo cuyo coseno viene dado por cosθ =hx, yi kxkkyk. 2.3.1. ´ Angulo entre subespacios. La siguiente definici´on, introducida originariamente por Friedrichs en 1937, es la m´as aceptada en la literatura referente a MAP para trabajar con el ´angulo entre dos subespacios. 2.3.2 Definici´on. El ´angulo θ(M, N) entre los dos subespacios cerrados My Nde Hes el ´angulo en [0,π 2] cuyo coseno c(M, N) esta dado por: sup{|hx, yi| :x∈M∩(M∩N)⊥,kxk ≤ 1, y ∈N∩(M∩N)⊥,kyk ≤ 1} 14 Es frecuente ver que otros autores definen el ´angulo θ(M, N) sin considerar el factor (M∩N)⊥de la expresi´on anterior. 2.3.3 Definici´on. El ´angulo minimal θ0(M, N) entre MyNes el ´angulo en [0,π 2] cuyo coseno c0(M, N) esta dado por: sup{|hx, yi| :x∈M, kxk ≤ 1, y ∈N, kyk ≤ 1} 2.3.4 Nota. (a) Es claro que si M∩N={0}ambas definiciones coinciden, c0(M, N) = c(M, N). (b) A continuaci´on expondremos unas consecuencias que se deducen de forma inmediata de las definiciones: 1. 0 ≤c(M, N)≤c0(M, N)≤1. 2. c(M, N) = c(N, M) y c0(M, N) = c0(N, M). 3. c(M, N) = c0(M∩(M∩N)⊥, N ∩(M∩N)⊥). 4. |hx, yi| ≤ c0(M, N)kxkkykpara todo x∈M, y∈N(desigualdad de Schwarz refinada) En el siguiente lema incluimos algunas propiedades ´utiles: 2.3.5 Lema. 1. c(M, N) = c0(M, N ∩(M∩N)⊥) = c0(M∩(M∩N)⊥, N) 2. c0(M, N) = kPMPNk=kPMPNPNk1/2. 3. c(M, N) = kPMPN−PM∩Nk=kPMPNP(M∩N)⊥k. Demostraci´on. Demostraremos en primer lugar 1.: Si x∈M x=PM∩Nx+P(M∩N)⊥xcon P(M∩N)⊥x=x−PM∩Nx∈M. Sea y∈ N∩(M∩N)⊥. Entonces, hx, yi=hPM∩Nx+P(M∩N)⊥x, yi=hP(M∩N)⊥x, yi con P(M∩N)⊥x∈M∩(M∩N)⊥, luego: {|hx, yi| :x∈M, kxk ≤ 1, y ∈N∩(M∩N)⊥,kyk ≤ 1}= {|hx, yi| :x∈M∩(M∩N)⊥,kxk ≤ 1, y ∈N∩(M∩N)⊥,kyk ≤ 1} Pasaremos ahora a demostrar 3. En primer lugar debemos ver que PM∩N conmuta con PM. Basta ver que PMPM∩N=PM∩(M∩N)=PM∩N. Tenemos que PM∩Nx∈Mpor lo tanto PMPM∩Nx=PM(PM∩Nx) = PM∩Ny ambos operadores conmutan. c(M, N) = c(M∩(M∩N)⊥, N ∩(M∩N)⊥ =kPM∩(M∩N)⊥PN∩(M∩N)⊥k =kPMP(M∩N)⊥PNP(M∩N)⊥k =kPMPN(I−PM∩N)k =kPMPN−PM∩Nk. 15 N´otese que combinado el apartado primero del lema previo y el punto 4 de la nota es trivial obtener de la desigualdad de Schwarz refinada para c(M, N): |hx, yi| ≤ c0(M, N)kxkkykpara todo x∈M,y∈N, cuando al menos xoypertenecen a (M∩N)⊥. Cabe resaltar, aunque no lo demostraremos, que si empleamos la primera definici´on de ´angulo se tiene que el ´angulo entre dos subespcios es el mismo que entre sus complementos ortogonales. Esto no es cierto si empleamos la segunda definici´on, de ah´ı nuestra preferencia por la primera. 2.3.6. Velocidad de convergencia de MAP. Del teorema de extensi´on de Halperin concluimos que (Pr...P2P1)nxconverge a PMxpara cada xen H( donde M=∩r i=1MiyPi=PMi, i = 1, ..., r). Sin embargo la tasa de convergencia puede ser arbitrariamente lenta. De hecho, es posible encontrar ejemplos incluso para r= 2 que ilustran la lentitud de MAP. Sin embargo, se han desarrollado t´ecnicas para acelerar dicha convergencia . En la siguiente secci´on analizaremos algunas de estas t´ecnicas en profundidad. Ahora trataremos de analizar la tasa de convergencia de MAP en subespacios. En primer lugar notamos que para cada i= 1, ..., r,PiPM=PMy que PiPM⊥=PMi∩M⊥(se tiene que PiPM⊥=Pi(I−PM) = Pi−PiPM= Pi−PMPi=PM⊥Pi). De aqu´ı se puede deducir que para todo x∈H: k(Pr...P2P1)nx−PMxk≤k(Pr...P2P1)n−PMkkxk =k(Pr...P2P1PM⊥)nkkxk ≤ k(Pr...P2P1PM⊥knkxk. De lo anterior se observa que la tasa de convergencia depende de la norma del operador Pr...P2P1PM⊥. En particular, para el caso en que r= 2 aplicando el lema (2.3.5) se deduce que k(P2P1)n−PMk≤kP2P1PM⊥kn=c(M1, M2)n. Sin embargo para el caso de dos subespacios esta cota es mejorable. Para cada x∈Hy para cada entero n≥1, se tiene k(P2P1)nx−PMxk ≤ c(M, N)2n−1kxk, y veremos en el siguiente teorema que c(M, N)2n−1coincide con la norma del operador (P2P1)n−PMy por lo tanto la cota no se puede mejorar. 2.3.7 Teorema. k(P2P1)n−PMk=c(M, N)2n−1(n= 1,2, ...). 16 Demostraci´on. En primer lugar introduciremos la siguiente notaci´on Qi= PiPM⊥=PMi∩M⊥donde i= 1,2. Entonces kQ2Q1k=kP2PM⊥P1PM⊥k= kP2P1PM⊥k. Ahora se tiene que, [(Q2Q1)n]∗= [(Q2Q1)∗]n= (Q2Q1)n, entonces k(Q2Q1)nk2=k(Q2Q1)n[(Q2Q1)n]∗k =k(Q2Q1)n(Q2Q1)nk =k(Q2Q1Q2)2n−1k, y como Q2Q1Q2es autoadjunto y por lo tanto normal, se sigue que k(Q2Q1Q2)2n−1k=kQ2Q1Q2k2n−1. Adem´as, kQ2Q1Q2k=kQ2Q1Q1Q2k=k(Q2Q1)(Q2Q1)∗k=kQ2Q1k2. Por consiguiente, k(Q2Q1)nk2=kQ2Q1Q2k2n−1=kQ2Q1k2(2n−1), Finalmente, k(Q2Q1)nk=kQ2Q1k2n−1. y el resultado que quer´ıamos probar se deduce del lemma (2.3.5) La tasa de convergencia de MAP se puede especificar en t´erminos de los ´angulos entre los subespacios involucrados, por el contrario para el caso r≥2 no podemos presentar una expresi´on exacta para la norma del operador en t´erminos de los ´angulos. Sin embargo, proporcionar una cota superior si que nos es posible. 2.3.8 Teorema. Para cada i= 1, ..., r sea Miun subespacio cerrado de H. Entonces para cada x∈Hy para cada entero n≥1, se tiene que k(PMrPMr−1...PM1)nx−P∩r i=1Mixk=cn/2kx−P∩r i=1Mixk, donde c= 1 − r−1 Q i=1 sin2θi, θies el ´angulo entre Miy∩r j=i+1Mj. 17 Demostraci´on. Denotemos por Mla intersecci´on de los distintos subespacios ∩r i=1Mipor P=PMrPMr−1...PM1y por y=P∩r i=1Mix. Ser´a suficiente probar que kPnx−yk2≤cnkx−yk2. teniendo en cuenta que y∈MyPes la identidad en Mla desigualdad anterior puede escribirse com sigue kPn(x−y)k2≤cnkx−yk2. Llamando v=x−y(con v∈M⊥), se tiene que basta con probar kPvk2≤ckvk2. Lo demostraremos por inducci´on sobre r. Si r= 1 es claro que se verifica nuestro resultado. Sea M0=Mr∩Mr−1∩ ... ∩M2yP0=PMrPMr−1...P2. Para todo v∈M⊥se escribe v=w+v1 con w∈M1y con v1∈M⊥ 1, Pv=P0 w. A continuaci´on escribimos w=w0+w00 , con w0∈M0y con w00 ∈M0⊥, de tal forma que P0w=w0+P0w00 . Por otro lado tenemos que hP0w00 , w0i=hw00 , PM2PM3...PMrw0i=hw00 , w0i= 0, Vemos que P0w00 yw0son ortogonales, por el teorema de Pit´agoras se sigue que kP0wk2=kw0k2+kP0w00 k2. Finalmente por la hip´otesis de inducci´on se tiene que kP0w00 k2≤"1− r−1 Y i=2 sin2θi#kw00 k2. 18 De las ´ultimas expresiones se obtiene que kP0wk2≤ kw0k2+"1− r−1 Y i=2 sin2θi#kw00 k2(2.7) =kw0k2+"1− r−1 Y i=2 sin2θi#kwk2− kw0k2(2.8) ="1− r−1 Y i=2 sin2θi#kwk2+ r−1 Y i=2 sin2θikw0k2.(2.9) Por otra parte, se puede escribir w=v−v1con v∈M⊥yv1∈M⊥ 1, para todo aen M, hw, ai=hv−v1, ai=hv, ai−hv1, ai= 0, por lo tanto se tiene que w∈M1es ortogonal a M=M1∩M0. Adem´as, w0=w−w00 , con w⊥M=M1∩M0yw00 ∈M0⊥, para todo a∈M, hw0, ai=hw−w00 , ai=hw, ai−hw00 , ai= 0; luego, w0∈M0y es ortogonal a M=M0∩M1. Teniendo en cuenta que el ´angulo entre MyM0es como m´ınimo θ1, se sigue que kw0k2=hw0, w0i=hw−w00 , w0i=hw, w0i ≤ cosθ1kwkkw0k; entonces kw0k ≤ cosθ1kwk. Reemplazando esta ´ultima expresi´on en (2.9) y operando kP0wk2≤"1− r−1 Y i=2 sin2θi#kwk2+ r−1 Y i=2 sin2θi1−sin2θ1kwk2 ="1− r−1 Y i=1 sin2θi#kwk2. Para finalizar, como Pv =P0wykwk ≤ kvkse cumple por tanto que kPvk2≤ckvk2y la demostraci´on queda finalizada. En el caso en que la sucesi´on {xk} ⊂ Hconverja hacia un x∗∈H, estaremos interesados en incluir una serie de definiciones relacionadas con la velocidad de convergencia y que involucraran a los vectores de error ek= xk−x∗. 2.3.9 Definici´on. Diremos que la sucesi´on ekconverge hacia 0 con q−orden p si existe una constante c > 0 y un k0∈N, tal que 19 de Kaczmarz. En vez de proyectar los sucesivos iterantes de forma c´ıclica en los hiperplanos Hi,i= 1, . . . , m, en el algoritmo de Cimmino se proyecta x0simultaneamente en todos los subespacios Hi,i= 1, . . . , m, y se toma como nuevo iterante el promedio de estas proyecciones. Denotando con PHi al operador proyecci´on ortogonal sobre el hiperplano Hi, uno tiene xk+1 =1 m m X i=1 PHi(xk). Los m´etodos de Cimmino y Kaczmarz se aplican a la resoluci´on de Ax =b, produciendo una sucesi´on {xk}que converge a un punto de ∩m i=1Hi, intersecci´on de los hiperplanos Hi. El algoritmo de Kaczmarz con par´ametro de relajaci´on ωk= 1, despu´es de un ciclo de proyecciones sobre los Hi, i= 1, . . . , m, equivale a computar la imagen de xkpor el operador TK=PHmPHm−1· · · PH1, mientras que para el m´etodo de Cimmino, el operador involucrado es TC=1 m m X i=1 PHi. Aunque cada uno de los operadores de proyecci´on PHi,i= 1, . . . , m, es autoadjunto, no ocurre lo mismo con el operador TKque corresponde a un ciclo de iteraciones del m´etodo de Kaczmarz. Sin embargo TCsi es autoadjunto, una circunstancia que puede explotarse para producir esquemas de aceleraci´on de la convergencia. 3.2.3. El M´etodo de Relajaci´on de Agmon, Motzkin y Schoenberg (MAMS) En el m´etodo de relajaci´on de Agmon, Motzkin y Schoenberg se considera la soluci´on del sistema de inecuaciones Ax ≤b, con Auna matriz m×n, y b∈ Rm. Este problema equivale al problema de la factibilidad en programaci´on lineal. Al igual que con el m´etodo de Kaczmarz este problema se puede generalizar a espacios de Hilbert, puesto que se trata de encontrar un xen la intersecci´on de mhemiespacios que vendran dados por: Si={x∈H:hai, xi ≤ bi} para cada i∈M. El algoritmo toma la aproximaci´on inicial x0arbitraria y calcula: xk+1 =xk+δkaik, donde δk= min 0, ωk bik− haik, xki haik, aiki. 26 Como en el algoritmo de Kaczmarz, los factores de relajaci´on ωkvar´ıan en el rango 0 <  ≤ωk≤2− < 2 para todo k, con peque˜no y positivo. La sucesi´on de control puede ser la c´ıclica o, alternativamente, casi c´ıclica. La interpretaci´on geom´etrica es an´aloga a la del m´etodo de Kaczmarz: con factores de relajaci´on ωk= 1 para todo k, el algoritmo equivale a realizar de forma c´ıclica sobre los distintos hemiespacios del sistema proyecciones ortogonales de los sucesivos iterantes hasta alcanzar la convergencia. N´otese que la proyecci´on de zsobre el hemiespacio Sies la proyecci´on ortogonal sobre el hiperplano Hi, cuando zno est´a en el hemiespacio, o se reduce al propio zsi z∈Si. Es importante se˜nalar que en este tipo de m´etodos no se garantiza la convergenc´ıa al elemento del conjunto de soluciones factibles m´as pr´oximo al iterante inicial x0. Figura 3.2: Interpretaci´on geom´etrica del m´etodo MAMS con ω= 1 3.2.4. El algoritmo de Hildreth El algoritmo de Hildreth determina la soluci´on de norma eucl´ıdea m´ınima del conjunto de inecuaciones lineales Ax ≤b, con Auna matriz real m×n, yb∈Rm. El algoritmo parte de una aproximaci´on inicial z0∈Rn +, y define x0= −ATz0. Entonces en cada iteraci´on computa xk+1 =xk+δkaik, zk+1 =zk−δkeik, 27 con δk= m´ın (zk)ik, ωk bik− haik, xki haik, aiki, yeikel vector de la base can´onica, que tiene 1 en la coordenada ik-´esima. De nuevo, el factor de relajaci´on se toma en el rango 0 <  ≤ωk≤2− < 2. La sucesi´on de control puede ser c´ıclica o casi c´ıclica. La interpretaci´on geom´etrica es an´aloga a la del m´etodo de Agmon, Motzkin y Schoenberg, con la diferencia de que si xkest´a en el interior del hemiespacio Sikentonces xk+1 o es la proyecci´on ortogonal sobre Hiko se queda en el interior de Sik en el rayo ortogonal desde xkal hiperplano Hik. La sucesi´on de control es c´ıclica o casi c´ıclica. El m´etodo anterior se modifica para tratar la situaci´on frecuente en la que se tienen desigualdades lineales en ambos sentidos; es decir, las desigualdades a resolver aparecen como ci≤ hai, xi ≤ bi, i = 1, . . . , m. El algoritmo determina todav´ıa la soluci´on de norma eucl´ıdea m´ınima y procede del siguiente modo: Se pone inicialmente x0= 0 y z0= 0. Entonces en cada paso se computa xk+1 =xk+δkaik, zk+1 =zk−δkeik, con δk= med (zk)ik,bik− haik, xki haik, aiki,cik− haik, xki haik, aiki donde med(a, b, c) denota la mediana de los tres valores a, b, c. El m´etodo se ha probado convergente s´olo en el caso en que los par´ametros de relajaci´on ωkson igual a 1. 3.2.5. El m´etodo de relajaci´on de Herman El problema que se plantea es encontrar x∈Rntal que bi−αi≤ hai, xi ≤ bi+αi, i = 1, . . . , m. Este es un problema de factibilidad lineal. Goffin observ´o que para obtener convergencia r´apida una estrategia conveniente es usar la proyecci´on ortogonal cuando los iterantes est´an alejados de la banda definida por las restricciones de intervalo anteriores, y usar reflexiones cuando los iterantes est´an ya pr´oximos a dicha banda. El m´etodo propuesto por Herman procede de la siguiente forma. Se toma x0∈Rnarbitrario, y en cada paso se computa xk+1 =xk+δk aik kaikk2, 28 donde δk=           0 si |bi− hai, xki| ≤ αi, bi− hai, xkisi |bi− hai, xki| ≥ 2αi, 2(bi+αi− hai, xki) si bi+αi<|hai, xki| < bi+ 2αi, 2(−bi+αi+hai, xki) si bi−2αi<|hai, xki| < bi−αi, 3.3. El m´etodo de Kaczmarz y la inversa generalizada. En esta secci´on ampliamos la teor´ıa de convergencia del algoritmo de Kaczmarz para considerar el problema de la aproximaci´on de la pseusdoinversa de Penrose de una matriz A,m×n. La teor´ıa de convergencia que aporta el teorema de proyecciones alternadas de Von Neumann debe extenderse para cubrir esta situaci´on. Consideramos, pues, el sistema de ecuaciones lineales Ax =b donde Aes una matriz real de dimensiones m×n(posiblemente m>n) y xybson vectores columna reales ny m-dimensionales respectivamente. Dada la matriz Ay un vector bla soluci´on de dicho sistema no tiene por que existir necesariamente. El problema es decir cuando el sistema tiene soluci´on o no y tambi´en encontrar todos los vectores soluci´on si estos existen. A lo largo de esta secci´on veremos que el m´etodo de acci´on por fila de Kaczmarz converge para cualquier sistema de ecuaciones lineales con filas no nulas, incluso cuando el sistema es singular e inconsistente. Adem´as, el n´umero de operaciones requeridas en cada iteraci´on del m´etodo es peque˜no. Finalmente, el m´etodo de Kaczmarz nos proporcionar´a un algoritmo para calcular la inversa generalizada de la matriz A. Como ya es usual, ATdenotar´a la transpuesta de A. Debemos destacar que toda matriz puede ser identificada como un operador lineal acotado entre espacios de Hilbert. De esta manera, dada una matriz real Ade m filas y ncolumnas la notaci´on kAkhar´a referencia a la norma del operador lineal acotado asociado a la matriz A, es decir: kAk= sup kAxk kxk:x∈Rn. Tal y como hemos visto previamente el sistema lineal Ax =bpuede ser descrito mediante productos internos de la siguiente manera Hi=hx, aii=bipara i= 1, ..., m. donde aies la i-´esima columna de nuestra matriz AT. PHi=x−hx, aii − bi αi aipara i= 1, ..., m 29 Denota la proyecci´on sobre el hiperplano Hiyαi=hai, aii. Un ciclo completo de proyecciones sucesivas del m´etodo de Kaczmarz puede ser descrito mediante el sigiente operador F(b, x) = PH1◦PH2◦PH3... ◦PHm−1◦PHm =PH1(...(PHm−2(PHm−1(PHm(x)))...). De esta manera podemos generar nuestra ya conocida sucesi´on de iterantes mediante la relaci´on de recurrencia siguiente xi+1 =F(b, xi)i= 0,1,2, ... (3.1) Hemos visto previamente que si el sistema es compatible determinado la sucesi´on anterior converge a la soluci´on de dicho sistema. Ahora demostraremos que la sucesi´on de iterantes del m´etodo de Kaczmarz converge cualquiera que sea el iterante inicial x0, para cualquiera que sean Ayb. Cada proyecci´on PHipuede expresarse en t´erminos de la matriz Ade la siguiente manera PHi=Pix+bi αi aipara i= 1, ..., m, donde Pidenota el operador de proyecci´on sobre el hiperplano definido por hx, aii= 0 para i= 1, ..., m. A su vez cada Pipuede escribirse de manera matricial como sigue Pi=I−1 αi aiaT i=δkl −aikail αi. De esta manera vemos que cada PHiconstituye una transformaci´on af´ın. Sea Qi=P1P2...Pi(i= 1, ..., m), donde Q0=I. Sea Rla matriz real de nfilas y de mcolumnas cuyo i-´esimo vector columna viene dado por la expresi´on: 1 αi Qi−1ai. El producto matricial Rb adquiere entonces la siguiente forma Rb = m X i=1 bi αi Qi−1ai. Finalmente, nuestro ciclo completo de proyecciones puede escribirse: F(b, x) = Qx +Rb, donde Q=Qmy la matriz Rdepende ´unica y exclusivamente de la matriz del problema A. 30 3.3.1 Proposici´on. Q+RA =I Demostraci´on. La i-´esima columna de Ry la i-´esima fila de Ason respectivamente, 1 αi Qi−1aiyaT i. Se tiene entonces que. RA =1 α1 a1aT 1+1 α2 Q1a2aT 2+1 α3 Q2a3aT 3+... +1 αm Qm−1amaT m = (I−P1) + Q1(I−P2) + ... +Qm−1(I−Pm) =I−Qm, puesto que Qi−1Pi=Qipor definici´on. Nuestro espacio de Hilbert Rnpuede ser descrito mediante suma directa de los subespacios Ker(A) y Im(AT), es decir Rn= Ker(A)⊕Im(AT). donde Im(AT) es el subespacio generado por los vectores fila de A. 3.3.2 Lema. Ker(A) = m \ i−1 {x∈Rn; Pix=x} Im(AT) = *m [ i−1 {x∈Rn; Pix=0}+, 3.3.3 Lema. kQxk=kxksi, y solo si x∈Ker(A). Demostraci´on. Si xno pertenece a Ker(A), entonces existe un entero i0≤m tal que Pi0x6=x. Tomemos i0como el mayor de estos n´umeros. Entonces, kPi0Pi0+1...Pmxk=kPi0xk<kxk. Como cada proyecci´on Pitiene norma exactamente 1, se tiene que para i= 1,2, ..., m kQik=kP1P2...Pik≤kP1kkP2k...kPik= 1. Por lo tanto, kQmxk=kQi0−1kkPi0Pi0+1...Pmxk<kxk. Por otro lado, si x∈Ker(A), entonces Pix=xpara i= 1,2, ..., m; de donde se obtiene Qx =P1P2...Pmx=x. Concluy´endose as´ı la demostraci´on 31 3.3.4 Corolario. kQk ≤ 1.Si el rango de A<nentonces kQk= 1. 3.3.5 Corolario. Qx =xsi, y solo si x∈Ker(A). 3.3.6 Teorema. 1. Ker(A) e Im(AT) son subespacios invariantes por la aplicaci´on lineal Q, adem´as Q|Ker(A)=I, es decir Q=PKer(A)⊕e QyPKer(A)e Q=e QPKer(A)= 0, donde e Q=QPIm(AT). 2. ke Qk= sup {kQxk<1}cuando x∈Im(AT) y kxk= 1. Demostraci´on. De la demostraci´on del lema anterior tenemos Q|Ker(A)= IKer(A). De forma similar se tiene QT|Ker(A)=IKer(A), puesto que QT= PmPm−1...P1.Entonces para todo x∈Ker(A) y para todo y∈Im(AT) tenemos hx, Qyi=hQTx, yi=hx, yi= 0. As´ı Qy ∈(Ker(A))⊥= Im(AT). Teniendo en cuenta que PKer(A)+PIm(AT)= IyPKer(A)PIm(AT)=PIm(AT)PKer(A)= 0, tenemos demostrada la primera parte del teorema. Del corolario 3.3.4 se sique que ke Qk ≤ kQk ≤ 1. Si ke Qk= 1, entonces existe un vector no nulo x0∈Im(AT) tal que kQx0k= kx0k, esto es debido a que el operador kQxkes continuo en el conjunto compacto x∈Im(AT); kxk= 1. Del lema 3.3.3 se sigue que x0∈Ker(A), y por lo tanto x0= 0. Esto claramente contradice nuestra hip´otesis con lo que queda probado nuestro resultado. N´otese que si el rango de Aes n, entonces Q=e QykQk<1. Puesto que Q=PKer(A)⊕e Qse deduce el siquiente corolario 3.3.7 Corolario. l´ım i→∞ Qi=PKer(A). Dados una serie de vectores no nulos x, a1, a2, ..., am(posiblemente dependientes), es un problema frecuente el calcular la componente de xortogonal a cada uno de los vectores a1, a2, ..., am:PKer(A)x. El siguiente corolario nos proporciona un algoritmo basado en el m´etodo de Kaczmarz para solventar dicho problema. 32 3.3.8 Corolario. Sea b= 0, el algoritmo (3.1) xi+1 =F(0, xi) = Qxi, i = 0,1,2,3, ... (3.2) genera una secuencia de iterantes xila cual converge a PKer(A), donde x0es un vector incial. Aplicando el proceso iterativo (3.2) del corolario anterior a los nvectores iniciales que forman la base can´onica de Rn: e1= [1,0, ..., 0]T,e2= [0,1, ..., 0]T, ..., en= [0,0, ..., 1]T, podemos calcular la matriz del operador de proyecci´on PKer(A), esto se debe a que el vector PKer(A)eiconstituye la i-´esima columna de dicha matriz. 3.3.9 Teorema. l´ım i→∞ i P j=0 QiRexiste y adem´as se tiene que: I−Q0−1R= ∞ X j=0 QjR. Demostraci´on. Para todo x∈Ker(A), se tiene que para i= 1,2, ..., m hQi−1ai, xi=hai, QT i−1xi=hai, xi= 0, puesto que QT i=PiPi−1...P1. Entonces los vectores columna 1 αi Qi−1aide la matriz R estan en el subespacio Im(AT). Por consiguiente, se tiene que i X j=0 QjR= i X j=0 e QjR. Puesto que ke Qk<1, tenemos ∞ X j=0 e QjR=  ∞ X j=0 e Qj R= (I−e Q)−1R. Finalmente, (I−e Q)−1R= ∞ X j=0 QjR. 33 3.3.10 Teorema. La matriz de nfilas y mcolumnas G= (I−e Q)−1R es una inversa generalizada de la matriz Ade nuestro problema, es decir: AGA =A, GAG =G, GA =PKer(A), AG =P, donde Pes la proyecci´on sobre Im(A)a trav´es de Ker(R).GA =PKer(A) puede ser reemplazada por (GA)T=GA. Demostraci´on. Es claro que los vectores que constituyen las columnas de ATestan en el subespacio Im(AT), por lo tanto, (I−Q)AT= (I−e Q)AT. Por la proposici´on (3.3.1), se tiene que RA =I−Q. Entonces. RAAT= (I−e Q)AT, de donde (I−e Q)RAAT=AT. Esto quiere decir que AGA =Ay (GA)T=GA. Tenemos que RA(I−e Q)−1= (I−Q)(I−e Q)−1 = (I−e Q−PKer(A))(I−e Q)−1 =I−PKer(A)((I−e Q)−1 =I−PKer(A), ya que PKer(A)e Q= 0. Entonces tenemos RA(I−e Q)−1R= (I−PKer(A))R=R, puesto que los vectores columna de Restan en el subespacio Im(AT). Premultiplicando la ecuaci´on anterior por (I−e Q)−1, tenemos (I−e Q)−1RA(I−e Q)−1R= (I−e Q)−1R, 34 esto es, GAG =G. Vemos que AG es idempotente y Im(A)⊂Im(AG), ya que AGA =A. Es claro que Ker(AG) = Ker(A(I−e Q)−1R)⊃Ker(R). Como Im(AG)⊕Ker(AG) = Rm, entonces Im(A)∩Ker(R) = {0}. Pero Im(AT) = Im(R), porque el i-´esimo vector columna de la matriz Res de la forma 1 αi Qi−1ai=1 αi ai− i−1 X j=1 cjaj . Entonces dim(Im(A)) = dim(Im(RT)), por lo tanto, dim(Ker(AT)) = dim(Ker(R)). Finalmente deducimos que Im(A)⊕Ker(R) = Rn. lo que completa nuestra demostraci´on. Acto seguido probaremos la convergencia del algoritmo (3.1) 3.3.11 Corolario. Para toda matriz m×nreal Acon filas no nulas, y para cada vector m-dimensional b, el algoritmo proporcionado por el m´etodo de Kaczmarz (3.1) genera una sucesi´on convergente de iterantes {xi}tal que l´ım i→∞ xi=PKer(A)x0+Gb. donde x0∈Rnes un vector inicial arbitrario Demostraci´on. de la relaci´on de recurrencia dada por el algoritmo (3.1) y de la expresi´on (3.3), se deduce que xi=Qix0+  i−1 X j=0 QjR b. Es consecuencia directa del corolario (3.3.7) y del teorema previo que el primero y el segundo t´ermino de (3.3.11) converjan a PKer(A)x0y a Gb, respectivamente, cuando itiende a infinito. Es ya conocido que nuestro sistema de partida Ax =btiene soluci´on si, y s´olo, si AGb =b. En este caso, PKer(A)x0+Gb es soluci´on del sistema para el iterante inicial arbitrario x0, y Gb es la soluci´on de norma m´ınima. Entonces, el procedimiento para resolver el sistema Ax =bser´a el siguiente: 35 M´etodo de Kaczmarz tol 10−610−910−11 N. Ciclos 4 5 6 x10.9999999723706 0.9999999998321 0.9999999999999 x21.0000000260756 1.0000000001203 0.9999999999996 x31.0000000095306 1.0000000000498 0.9999999999999 x40.9999999970802 0.9999999999625 0.9999999999997 cuando se aplica el algoritmo de Kaczmarz (la norma infinito de dos iterantes internos consecutivos es menor que la tolerancia), seg´un diferentes tolerancias. El iterante inicial es siempre el vector de ceros. Con una tolerancia de 10−5se logran todas las cifras significativas de la soluci´on en 8 ciclos. Ejemplo 2. Consideramos ahora el comportamiento del m´etodo de Cimmino en relaci´on con el sistema del ejemplo anterior. M´etodo de Cimmino tol 10−610−910−11 N. Ciclos 14 20 24 x10.9999999851103 1.0000000000117 1.0000000000002 x20.9999997480712 0.9999999997751 0.9999999999976 x30.9999999873753 1.0000000000290 1.0000000000004 x40.9999996838586 0.9999999997327 0.9999999999973 Con una tolerancia de 10−5se logran todas las cifras significativas de la soluci´on en 33 ciclos. Obviamente, el m´etodo de Kaczmarz resulta ventajoso frente al de Cimmino. Ejemplo 3. Consideramos ahora el sistema lineal Ax =b, con matriz cuadrada 84 ×84 A=            6 1 8 6 1 8 6 1 ......... 8 6 1 8 6 1 8 6            y vector bde forma que la soluci´on sea la que tiene todas sus componentes igual a 1 (tomada de [9]). El algoritmo de Kaczmarz, con las mismas tolerancias de los ejemplos anteriores, da convergencia en los ciclos que se 42 M´etodo de Kaczmarz tol 10−610−910−11 N. Ciclos 38 69 91 x11.0000002706247 1.0000000003039 1.0000000000030 x51.0000005142534 1.0000000006798 0.9999999999996 x80 0.9646853769984 0.9646809944622 0.9646809896319 x84 0.7083318423493 0.7083333315746 0.7083333333154 indican (la tabla s´olo presenta valores num´ericos para las componentes x1, x5,x80 yx84. Un inter´es adicional de este ejemplo es que la soluci´on con eliminaci´on Gaussiana y sin pivotaje (ni refinamiento iterativo) produce soluciones num´ericas para x84 del orden de 1013. Sin embargo, estos m´etodos iterativos dan aproximaciones que aunque no tienen m´as que una cifra significativa correcta (para x84), son mejores que las proporcionadas por la eliminaci´on Gaussiana. Ejemplo 4. Consideramos el siguiente sistema lineal Ax =b, con matriz 6×4.         1.0 3.0 2.0−1.0 1.0 2.0−1.0 2.0 1.0−1.0 2.0 3.0 2.0 1.0 1.0 1.0 5.0 5.0 4.0 1.0 4.0−1.0 5.0 7.0         x=         5.0 0.0 5.0 5.0 15.0 15.0         El rango de dicha matriz es 3 y el sistema es consistente. La soluci´on computada por el m´etodo de Kaczmarz se muestra en la siguiente tabla. El vector empleado como iterante inical es x= [7,6,10,6]T. La soluci´on real es x= [1,1,1,1]T. M´etodo de Kaczmarz tol 10−610−910−11 N. Ciclos 38 60 74 x11.0000100433911 1.0000000091225 1.0000000001057 x20.9999981753989 0.9999999983426 0.9999999999807 x30.9999894855657 0.9999999904495 0.9999999998892 x41.0000015105722 1.0000000013720 1.0000000000159 Ejemplo 5. Encontrar la inversa generalizada Gde la matriz   1.0 0.0−1.0 1.0 0.0 1.0 1.0 0.0 1.0 0.0 1.0 1.0  43 Mostramos los resultados obtenidos mediante la aplicaci´on del m´etodo de Kaczmarz con distintas tolerancias. Para una tolerancia de 106:     0,249999958084479 −0,000000038852564 0,250000060293264 0,499999850715657 0,999999861624540 −0,499999785262356 −0,499999916168957 0,000000077705129 0,499999879413472 0,249999958084479 −0,000000038852564 0,250000060293264     Para tolerancia de 109:     0,249999999952388 −0,000000000044133 0,250000000040654 0,499999999830426 0,999999999842817 −0,499999999855209 −0,499999999904775 0,000000000088266 0,499999999918692 0,249999999952388 −0,000000000044133 0,250000000040654     Finalmente, para 1011:     0,249999999999564 −0,000000000000404 0,250000000000627 0,499999999998448 0,999999999998562 −0,499999999997768 −0,499999999999129 0,000000000000808 0,499999999998747 0,249999999999564 −0,000000000000404 0,250000000000627     Mientras, la soluci´on exacta es:     0.25 0.0 0.25 0.5 1.0−0.5 −0.5 0.0 0.5 0.25 0.0 0.25     3.6. Conclusiones Inicialmente se pens´o en la aplicaci´on de los m´etodos de acci´on por filas a grandes sistemas lineales como los que aparecen en la discretizaci´on de ecuaciones en derivadas parciales. Una ventaja potencial de los m´etodos de acci´on por fila es su eficiencia cuando se utilizan con matrices dispersas, ya que su implementaci´on s´olo requiere del c´alculo de productos internos de las filas de estas matrices por vectores y ´estos pueden calcularse econ´omicamente. Sin embargo, la memoria no refleja algunas experiencias negativas con el uso de estos m´etodos con sistemas como el que surge al discretizar el Laplaciano en el cuadrado unidad; los m´etodos adolecen de una convergencia exasperantemente lenta y no parecen a priori competitivos con ninguno de los m´etodos iterativos cl´asicos. 44 En las publicaciones sobre estos m´etodos que hemos revisado no hemos encontrado ninguna referencia a la elecci´on del par´ametro de relajaci´on ω que aparece en algunas de sus formulaciones. Los experiencias negativas que hemos tenido podr´ıan aliviarse de contar con una teor´ıa razonable para la elecci´on de dichos par´ametros. Las aplicaciones de estos m´etodos parecen concentrarse en matrices llenas de tama˜no moderado. Un ´area donde el m´etodo de Kaczmarz tiene bastante vigencia es el de la reconstrucci´on algebraica de im´agenes a partir de proyecciones de las mismas, sobre todo en relaci´on con im´agenes m´edicas. Hemos presentado el contexto de este tipo de aplicaciones, aunque el uso actual de estas t´ecnicas involucran datos con incertidumbre y requieren combinar estos m´etodos con t´ecnicas estad´ısticas. Proseguir esta l´ınea de desarrollo no entraba en los planes de este trabajo y hubiera requerido tiempo y espacio no disponibles. 45 46 Bibliograf´ıa [1] H. Br´ezis, An´alisis Funcional, Alianza Universidad Textos, Madrid, Espa˜na, 1984. (Cited on p. 1) [2] Y. Censor, Row-Action Methods for Huge and Sparse Systems and their Applications, SIAM Review 23(4), 1981, pp. 444–466. (Cited on pp. iii, 23) [3] R. Escalante y M. Raydan, Alternation Proyection Methods, SIAM, Philadephia, USA, 2011. (Cited on pp. iii, 23) [4] R. Gordon, R. Berman y G. T. Herman, Algebraic Reconstruction Techniques (ART) for Three-dimensional Electron Microscopy and X-ray Photography, J. Theor. Biol. 29, 1970, pp. 471. (Cited on p. 36) [5] I. Halperin, The Product of Projections Operators, Acta Sci. Math. (Szeged), 1962, 23:93-99. (Cited on p. 11) [6] G. T. Herman, ART: Mathematics and Applications. A Report on the Mathematical Foundation and on the Applicability to Real Data of the Algebraic Reconstruction Techniques, J. Theor. Biol. 42, 1973, pp. 1–32. (Cited on p. 36) [7] B.V. Limaye, Functional Analysis, Wiley, USA, 1981. (Cited on p. 1) [8] J. Von. Neumann, Functional Operators. Vol.II. The Geometry of orthogonal Spaces., Princeton University Press, Princeton, USA, 1950. (Cited on p. 9) [9] K. Tanabe, Projection method for solving a singular system of linear equations and its applications, Numer. Math. 17, 1971, pp. 203–214. (Cited on pp. 41, 42) [10] A. Vera L´opez y P. Alegr´ıa Ezquerra, Un curso de An´alisis Funcional, AVL, Murcia, Espa˜na, 1997. (Cited on p. 1) 47 48 Ap´endice A Programas de Matlab En este ap´endice incluimos los programas Matlab de los m´etodos de Kaczmarz, Cimmino, Hildreth y MAMS que hemos utilizado. Kaczmarz.m 1function [solsp, niter] = Kaczmarz(Asp, bsp, x0, tol, nmax, omega) 2% 3nrows = size(Asp,1); 4for i = 1 : nrows 5norm(i) = norm(Asp(i,:))ˆ2; 6end 7dxnorm = 1.0e20; 8niter = 0; 9while ((dxnorm > tol) && (niter < nmax)) 10 for i = 1 : nrows 11 if (norm(Asp(i,:),inf) == 0.0) 12 break 13 end 14 x1 = x0+omega*((bsp(i)-Asp(i,:)*x0)/norm(i))*( Asp(i,:)'); 15 dxnorm = norm(x1-x0, inf); 16 x0 = x1; 17 niter = niter+1; 18 end 19 end 20 solsp = x0; 21 end 49 Cimmino.m 1function [solsp, niter] = Cimmino(Asp, bsp, x0, tol, nmax) 2% 3nrows = size(Asp,1); 4norm = zeros(nrows,1); 5for i = 1 : nrows 6norm(i) = norm(Asp(i,:))ˆ2; 7end 8dxnorm = 1.0e20; 9niter = 0; 10 while ((dxnorm > tol) && (niter < nmax)) 11 dx = zeros(size(x0)); 12 for i = 1 : nrows 13 if (norm(Asp(i,:),inf) == 0.0) 14 break 15 end 16 dx = dx + ((bsp(i) - Asp(i,:) *x0)/norm(i)) * (Asp(i,:)'); 17 end 18 dx = dx/nrows; 19 x0 = x0 + dx; 20 dxnorm = norm(dx,inf); 21 niter = niter+1; 22 end 23 solsp = x0; 24 end Hildreth.m 1function [solsp, niter] = Hildreth(Asp, bsp, z0, tol, nmax, omega) 2% HILDRETH encuentra la solucion de norma minima de 3% A x <= b. z0 tiene componentes positivas; 4% x0 = -AˆT z0 5% 6nrows = size(Asp,1); 7for i = 1 : nrows 8norm(i) = norm(Asp(i,:))ˆ2; 9end 10 x0 = - Asp'*z0; 11 dxnorm = 1.0e20; 50 12 niter = 0; 13 while ((dxnorm > tol) && (niter < nmax)) 14 for i = 1 : nrows 15 if (norm(Asp(i,:),inf) == 0.0) 16 break 17 end 18 delta = omega*((bsp(i) - Asp(i,:) *x0)/norm(i )); 19 if (delta < z0(i)) 20 x1 = x0 + delta *(Asp(i,:)'); 21 z1(i) = z0(i) - delta; 22 else 23 x1 = x0 + z0(i) *(Asp(i,:)'); 24 z1(i) = z0(i) - delta; 25 end 26 dxnorm = norm(x1-x0, inf); 27 x0 = x1; 28 z0 = z1; 29 niter = niter+1; 30 end 31 end 32 solsp = x0; 33 end MAMS.m 1function [solsp, niter] = MAMS(Asp, bsp, x0, tol, nmax , omega) 2% 3nrows = size(Asp,1); 4for i = 1 : nrows 5norm(i) = norm(Asp(i,:))ˆ2; 6end 7dxnorm = 1.0e20; 8niter = 0; 9while ((dxnorm > tol) && (niter < nmax)) 10 for i = 1 : nrows 11 if (norm(Asp(i,:),inf) == 0.0) 12 break 13 end 14 delta = omega*((bsp(i) - Asp(i,:) *x0)/norm(i )); 15 if (delta < 0.0) 51