unesp UNIVERSIDADE ESTADUAL PAULISTA Instituto de Biociências, Letras e Ciências Exatas DEPARTAMENTO DE CIÊNCIAS DE COMPUTAÇÃO E ESTATÍSTICA SOBRE UM MÉTODO ASSEMELHADO AO DE FRANCIS PARA A DETERMINAÇÃO DE AUTOVALORES DE MATRIZES Danilo Elias de Oliveira Dissertação de Mestrado Pós-Graduação em Matemática Aplicada Rua Cristovão Colombo, 2265 15054-000 - São José do Rio Preto - SP - Brasil Telefone: (017) 3221-2201 Fax: (017) 3221-2203 Sobre um método assemelhado ao de Francis para a determinação de autovalores de matrizes Danilo Elias de Oliveira Dissertação apresentada ao Instituto de Biociências, Letras e Ciências Exatas da Universidade Estadual Paulista “Júlio de Mesquita Filho”, Câmpus de São José do Rio Preto, São Paulo, para a obtenção do título de Mestre em Matemática Aplicada. Orientadora: Profa. Dra. Eliana Xavier Linhares de Andrade São José do Rio Preto Fevereiro de 2006 Aos meus pais, Gilberto e Júlia, e ao meu irmão, Vilmar, ofereço. À minha mocinha, Cristiane, dedico. Agradecimentos Primeiramente, devo agradecer a Deus por ter me dado a vida e oportunidade de poder conviver ao lado de pessoas maravilhosas; ao Espírito Santo por me iluminar e abençoar e à Nossa Senhora que, através da Legião de Maria, me proporcionou uma verdadeira experiência de vida. Aos meus pais e irmão que sempre me apoiaram e me incentivaram a cursar o mestrado. Além disso, por todo amor, carinho e atenção que tiveram comigo. À minha namorada Cristiane, querida mocinha, por todo carinho, atenção, paciência e compreensão mesmo nas horas em que mais precisei. Um agradecimento especial à Profa. Dra. Eliana Xavier Linhares de Andrade pela dedi- cação, paciência e amizade que sempre teve comigo. Ao meu amigo Fabiano por me incentivar a fazer o mestrado. Aos meus amigos e irmãos do Philadelpho, Ronaldo e Nilo, que sempre me apoiaram e incentivaram. A todos os professores e funcionários que de alguma forma contribuíram para a realização deste trabalho. Em especial aos funcionários do DCCE, Getúlio, Luiza, Olga e Marcos. A todos os meus amigos de Pós-Graduação, em especial, à turma do futebol. Às minhas queridas Mayla e Uly. À CAPES - Coordenação de Aperfeiçoamento de Pessoal de Nível Superior, pelo auxílio financeiro. i “ O livro do universo está escrito em linguagem matemática. ” Galileu Galilei ii Resumo O principal objetivo deste trabalho é apresentar, discutir as qualidades e desempenho e provar a convergência de um método iterativo para a solução numérica do problema de autova- lores de uma matriz, que chamamos de “Método Assemelhado ao de Francis (MAF)”. O método em questão distingue-se do QR de Francis pela maneira, mais simples e rápida, de se obter as matrizes ortogonais Qk, k = 1, 2, . . .. Apresentamos, também, uma comparação entre o MAF e os algoritmos QR de Francis e LR de Rutishauser. Palavras-chave: autovalor; algoritmo QR; algoritmo LR; Decomposição de Cholesky. iii Abstract The main purpose of this work is to presente, to discuss the qualities and performance and to prove the convergence of an iterative method for the numerical solution of the eigenvalue problem, that we have called the ”Método Assemelhado ao de Francis (MAF)´´. This method differs from the QR method of Francis by providing a simpler and faster technique of getting the unitary matrices Qk, k = 1, 2, . . . . We present, also, a comparison analises between the MAF and the QR of Francis and LR of Rutishauser algorithms. Keywords: eigenvalue; QR algorithm; LR algorithm; Cholesky decomposition. iv Sumário 1 Introdução 1 2 Pré-Requisitos 5 2.1 Espaço vetorial e produto interno . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 Normas de vetores e matrizes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.3 Tipos especiais de matrizes quadradas . . . . . . . . . . . . . . . . . . . . . . . . 8 2.4 Autovalores e autovetores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.5 Reflexões de Householder . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.6 Rotações de Givens . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.7 Eliminação de Gauss e Decomposição LU . . . . . . . . . . . . . . . . . . . . . . 23 2.8 Decomposição QR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.9 Decomposição de Cholesky . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3 Algoritmo LR de Rutishauser 30 3.1 Algoritmo LR de Rutishauser . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.1.1 Prova da convergência de As . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.1.2 Matrizes hermitianas definidas-positivas . . . . . . . . . . . . . . . . . . . 35 3.1.3 Autovalores complexos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.1.4 Algoritmo LR modificado . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.1.5 Matrizes na forma de Hessenberg . . . . . . . . . . . . . . . . . . . . . . . 42 3.1.6 Aceleração da convergência . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.1.7 Deflação de uma matriz . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.1.8 Shift na origem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 v Sumário 4 Algoritmo QR de Francis 50 4.1 Algoritmo QR de Francis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.1.1 Convergência do algoritmo QR . . . . . . . . . . . . . . . . . . . . . . . . 51 4.1.2 Prova formal da convergência . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.1.3 Autovalores desordenados . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 4.1.4 Autovalores de mesmo módulo . . . . . . . . . . . . . . . . . . . . . . . . 55 4.1.5 Decomposição de As . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.1.6 Procedimento prático . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.1.7 Autovalores múltiplos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.1.8 Matrizes simétricas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.1.9 Relação entre os algoritmos QR e LR . . . . . . . . . . . . . . . . . . . . . 61 4.1.10 Convergência do algoritmo LR-Cholesky . . . . . . . . . . . . . . . . . . . 62 4.2 Método assemelhado ao de Francis . . . . . . . . . . . . . . . . . . . . . . . . . . 64 5 Comparação dos métodos e conclusões 66 5.1 Análise comparativa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.1.1 Gráficos comparativos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 5.1.2 Tabelas comparativas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 5.2 Conclusões . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 Referências Bibliográficas 73 vi Capítulo 1 Introdução Com o desenvolvimento da teoria quântica nas décadas de 1920 e 1930, o estudo de matrizes tornou-se de grande importância. No final da década de 1940, alguns matemáticos se pergun- tavam como o computador digital poderia ser empregado para resolver o problema de autovalores de matrizes. A determinação dos autovalores de uma matriz tem merecido grande atenção, devido a sua grande aplicação aos diversos ramos das Ciências Aplicadas, como, por exemplo: • a teoria das vibrações, quer sejam mecânicas ou elétricas, dos tipos macroscópica ou mi- croscópica; • as vibrações de pontes ou outra estrutura sólida; • as vibrações das asas de um avião; • a análise do insumo-produto, introduzida por Leontief, com sua matriz, que liga indústrias individuais a todo o trabalho da Economia; • a oscilação de partículas atômicas e moleculares nas ondas mecânicas de Schroedinger; • a teoria dos operadores lineares diferenciais e integrais, etc. . . A obtenção dos autovalores λ de uma matriz quadrada de ordem n envolve cálculos com- plexos por serem os zeros do polinômio característico det(A − λI) = 0, que é um polinômio de grau n em λ. A busca de soluções para este problema acarretou o desenvolvimento de uma série de algoritmos, os mais diversos, cada vez mais poderosos e mais eficientes. O primeiro a surgir foi o mais óbvio, composto de apenas dois passos. Primeiro, calcula-se os coeficientes do polinômio característico e, então, as raízes desse polinômio. 1 Introdução Embora existam várias maneiras de se realizar o primeiro passo, ainda não existe um algoritmo eficaz para o segundo, ou seja, para matrizes de ordem maior do que 10. Isso porque, em geral, se os coeficientes dos polinômios sofrem pequenas alterações, suas raízes são alteradas. Por outro lado, os autovalores não são alterados se os coeficientes da matriz sofrem pequenas alterações. Dessa forma, surgiu uma outra alternativa para se calcular os autovalores de uma matriz. Utilizando-se transformações de similaridade, podemos obter uma matriz com vários elementos nulos e que possua os mesmos autovalores da matriz original. Uma matriz diagonal seria perfeita. Porém, se conseguirmos uma matriz triangular já é o suficiente, pois seus autovalores são os elementos da diagonal principal. Examinaremos, neste trabalho, o problema de obtenção dos autovalores de uma matriz quadrada A, de ordem n, através de transformações de similaridade, que possibilitem a criação de uma seqüência {Ak} de matrizes semelhantes, que tendam a uma matriz triangular superior ou diagonal (se A for simétrica), A∞, cujos autovalores são os elementos diagonais. Como a seqüência {Ak} é constituída de matrizes semelhantes, temos que A e A∞ são semelhantes e, portanto, possuem os mesmos autovalores, com as mesmas multiplicidades. Assim, a determi- nação dos autovalores de A resume-se à busca dos elementos da diagonal principal da matriz A∞. Deve-se a Jacobi, em 1846, o início desta linha de pesquisa. A seguir, Givens, em 1954, usando matrizes de rotação de Jacobi, estabeleceu um método para tridiagonalização de matrizes. Em 1958, Rutishauser publicou um método, atualmente conhecido como algoritmo LR [24], que, utilizando decomposições triangulares, cria uma seqüência de matrizes semelhantes que, sob certas condições, convergem para uma matriz triangular superior ou diagonal, cujos elementos da diagonal principal fornecem os autovalores. Entretanto, o algoritmo LR tem sua aplicação restrita, pois depende do fato de a matriz dada, bem como cada uma das demais matrizes da seqüência, ser decomponível em LU, onde L é uma matriz triangular inferior com todos os elementos da diagonal principal iguais a 1 e U é triangular superior. O Teorema de Schur (Teorema 2.8) assegura a existência dessa matriz U que, por ser unitária (U∗ = U−1), nos permite garantir a semelhança entre A e T . Entretanto, a afirmação da existência dessa matriz não soluciona nosso problema. 2 Introdução Em 1961, Francis [9], utilizando, também, matrizes unitárias, a partir da idéia original de Rutishauser, desenvolveu o algoritmo QR que, atualmente, é considerado como a mais poderosa técnica para a solução do problema de autovalor. Seu método modifica o algoritmo LR quanto à decomposição de cada matriz da seqüência. Ao invés da decomposição triangular (LU), utiliza a decomposição unitária-triangular (QR), onde Q é uma matriz unitária e R triangular superior. Para calcular a matriz Q, Francis utilizou as matrizes de rotação de Jacobi. Dessa maneira, Francis eliminou a restrição da aplicabilidade do processo LR, tornando, sob certas condições, possível a obtenção de seqüência de matrizes similares, como a idealizada por Rutishuaser. Até os dias de hoje, o algoritmo mais utilizado é o algoritmo QR de passo duplo, também desenvolvido por Francis [10]. Porém, em uma tentativa de deixá-lo mais rápido, surgem as vari- ações do algoritmo QR. Por exemplo, o algoritmo QR paralelo, o algoritmo SDC, os algoritmos γ-híbridos, 0 ≤ γ ≤ 1, e os algoritmos SR e HR. No algoritmo QR paralelo, apresentado por Kaufman em [15], n 2 transformações de sim- ilaridade podem ser realizadas simultaneamente. Comparando este com o algoritmo original, o algoritmo paralelo necessita do triplo de iterações do algoritmo original, porém, como as opera- ções são realizadas simultaneamente, o tempo computacional é menor. Os algoritmos SR [6] e HR [4] são membros da família dos algoritmos GR [27] para calcular os autovalores e espaços invariantes de matrizes. Os algoritmos SR e HR são de grande importân- cia, pois preservam certas características da matriz. O algoritmo SR preserva a propriedade das matrizes hamiltonianas, enquanto que o algoritmo HR mantém a propriedade de matrizes tridi- agonais sign-simétricas, isto é, matrizes simétricas a menos dos sinais dos elementos. Em um trabalho mais recente, Benner e Faßbender [3] conseguiram demonstrar que os dois algoritmos são equivalentes se tomados os shifts adequados. Em um trabalho conjunto com Demmel e Gu, Bai [2] desenvolveu o algoritmo Spectral Divide and Conquer (SDC), que fora baseado na teoria de Bulgakov, Godunov e Malyshev. Com armazenagem e custo aritmético razoáveis, este algoritmo pode ser aplicado ao problema de autovalor e ao problema generalizado de autovalor. Haag e Watkins [13] descreveram os algoritmos γ-híbridos que combinam elementos dos algoritmos QR e LR. Além disso, também desenvolveram um algoritmo semelhante ao LR de Rutishauser (ALR). O algoritmo 1-híbrido, com γ = 1, foi o que apresentou melhor desempenho, chegando a ser 15% mais rápido que o QR. 3 Introdução O objetivo deste trabalho é estudar um método equivalente ao QR de Francis e, portanto, de aplicação mais ampla. Este método difere-se do de Francis pela maneira mais simples e mais eficiente de se obter as matrizes Qk, k = 1, 2, . . ., das decomposições QR. Foi denominado "Método Assemelhado ao de Francis"(MAF) e apresentado por Andrade e Bracciali em [1]. O MAF é mais simples que o de Francis pois utiliza a decomposição de Cholesky para determinar as decomposições QkRk das matrizes Ak, para k = 1, 2, . . .. Para isso, em cada estágio k, a matriz Bk = At kAk é decomposta no produto Rt kRk e, assim, (Rt k) −1At kAk = Rk. Construímos, então, a matriz Qt k = (Rt k) −1At k. As matrizes Qk e Rk obtidas desta forma são as matrizes da decomposição QR da matriz Ak. Para desenvolver o estudo a que nos propusemos, organizamos esta dissertação da seguinte forma: No Capítulo 2, apresentamos alguns conceitos e resultados da Álgebra Linear no que diz respeito ao problema de autovalores de matrizes. Podemos dizer que este capítulo servirá de base para os demais. Nos Capítulo 3, fazemos um estudo sobre o algoritmo LR de Rutishauser. Mostramos sua convergência, eficácia e apresentamos o algoritmo LR modificado, uma alternativa para as matrizes que não possuem decomposição LU. No Capítulo 4, apresentamos o algoritmo QR de Francis e o MAF, Método Assemelhado ao de Francis, para obtenção dos autovalores de uma matriz. Mostramos suas convergências, eficácia dos métodos e a relação entre os algoritmo QR e LR. No capítulo 5, apresentamos os resultados e conclusões obtidas dos testes realizados, com- parando os algoritmos QR, LR e MAF. Além disso, apresentamos, também, uma das matrizes utilizadas na comparação e as matrizes obtidas com os algoritmos programados. 4 Capítulo 2 Pré-Requisitos Neste capítulo, apresentamos algumas definições de matrizes quadradas que desempenham pa- pel importante na Álgebra Linear, tais como as matrizes diagonais, triangulares, ortogonais, definidas-positivas e na forma de Hessenberg. Além disso, consideramos, também, as normas de vetores e matrizes bem como as definições de autovalor e autovetor. Nas seções 2.5 e 2.6, apresentamos as reflexões de Householder e as rotações de Givens, que podem ser encontrados em Hager [14]. Apresentamos, também, as decomposições LU, QR e de Cholesky. Todas as decomposições podem ser encontrados nos livros Lipschutz [19], Rutishauser [25] e Noble [22]. 2.1 Espaço vetorial e produto interno Definição 2.1. Dizemos que um conjunto V 6= ∅ é um espaço vetorial sobre um conjunto de escalares K se, para quaisquer u, v, w ∈ V e α, β ∈ K, (I) Existe uma operação adição V × V em V , (u, v) 7→ u + v, com as seguintes propriedades: • u + v = v + u (comutativa); • (u + v) + w = u + (v + w) (associativa); • Existe um vetor em V , chamado 0 ou vetor zero, tal que u + 0 = u (elemento neutro); • Para cada vetor u existe um vetor em V , denotado por −u, tal que u + (−u) = 0 (oposto). (II) Está definida uma multiplicação de K × V em V , o que significa que a cada par (α, u) de K×V está associado um único elemento de V (que se indica por αu) e, para essa multiplicação, tem-se o seguinte: 5 Espaços vetoriais e produto interno • α(u + v) = αu + αv; • (α + β)u = αu + βu; • (αβ)u = α(βu); • Existe um elemento de K, denotado por 1 ou escalar unitário, tal que 1u = u. Usa-se freqüentemente a notação Kn para denotar o conjunto de todas as n-uplas de elementos de K. Assumiremos que Kn é um espaço vetorial sobre K, onde a adição de vetores e a multiplicação por escalar se definem como (a1, a2, . . . , an) + (b1, b2, . . . , bn) = (a1 + b1, a2 + b2, . . . , an + bn), k(a1, a2, . . . , an) = (ka1, ka2, . . . , kan), onde k, ai, bi ∈ K, i = 1, . . . , n. Além disso, o vetor zero em Kn é a n-upla (0, 0, . . . , 0) e o negativo de um vetor é: −(a1, a2, . . . , an) = (−a1,−a2, . . . ,−an). Definição 2.2. Seja W um subconjunto de um espaço vetorial V sobre um corpo K. W é um subespaço de V se W é ele próprio um espaço vetorial sobre K em relação à adição de vetores e à multiplicação por escalar. Definição 2.3 (Produto Interno). Seja V um espaço vetorial real. Suponha que a cada par de vetores u, v ∈ V esteja associado um número real, denotado por < u, v >. Esta função é chamada produto interno (real) em V se satisfaz às seguintes propriedades: 1. < au1 + bu2, v >= a < u1, v > +b < u2, v >; 2. < u, v >=< v, u >; 3. < u, u >≥ 0 e < u, u >= 0 ⇔ u ≡ 0. Um espaço vetorial V onde está definido um produto interno é chamado espaço com produto interno (real). Como exemplo de um espaço com produto interno temos o Rn, com < u, v >= u1v1 + u2v2 + . . . + unvn, 6 Normas de vetores e matrizes onde u = (u1, u2, . . . , un) e v = (v1, v2, . . . , vn), ui, vi ∈ R, i = 1, . . . , n, n > 1. 2.2 Normas de vetores e matrizes Definição 2.4. Seja V um espaço vetorial sobre R e x ∈ V . Chama-se norma de x, ||x||, qualquer função definida em V , com valores em R, satisfazendo às seguintes condições: i) ||x|| ≥ 0 e ||x|| = 0 ⇔ x ≡ 0; ii) ||λx|| = |λ| ||x||, para todo escalar λ ∈ R; iii) ||x + y|| ≤ ||x||+ ||y||, para quaisquer x, y ∈ V (desigualdade triangular). Como exemplos de normas no Rn, temos a) ||x||2 = √√√√ n∑ i=1 x2 i , b) ||x||p = p √√√√ n∑ i=1 |xi|p, 0 < p < ∞ e c) ||x||∞ = max 1≤i≤n |xi|. De maneira similar, definimos norma de matrizes. Definição 2.5. Seja Rm×n o espaço vetorial das matrizes reais com m linhas e n colunas. Chama-se norma de uma matriz A ∈ Rm×n, ||A||, qualquer função definida em Rm×n, com valores em R, satisfazendo às seguintes condições: i) ||A|| ≥ 0 e ||A|| = 0 ⇔ A ≡ 0; ii) ||κA|| = |κ|||A||, para todo escalar κ ∈ R; iii) ||A + B|| ≤ ||A||+ ||B||, para quaisquer A,B ∈ Rm×n (desigualdade triangular). Vejamos alguns exemplos de normas de matrizes: a) ||A||E = √√√√ m,n∑ i,j=1 a2 ij , b) ||A||1 = max 1≤j≤n m∑ i=1 |aij | e c) ||A||∞ = max 1≤i≤m n∑ j=1 |aij |. Definição 2.6. Dada uma norma de vetor no Rn, podemos definir uma norma de matriz, que será subordinada a ela, por ||A|| = sup ||x||=1 ||Ax||. As normas de matrizes assim definidas podem ser interpretadas como sendo o comprimento do maior vetor do conjunto imagem Ax da esfera unitária {x, tal que||x|| = 1} pela transfor- mação x → Ax. Uma norma de matriz subordinada à uma norma de vetor possui a seguinte propriedade: 7 Tipos especiais de matrizes quadradas ||AB|| ≤ ||A|| ||B||. Em nosso estudo, vamos considerar apenas as normas subordinadas de matrizes e que satisfazem à propriedade abaixo. Definição 2.7. Dada uma norma no Rn e uma norma de matriz, dizemos que elas são consis- tentes se, para qualquer x ∈ Rn, ||Ax|| ≤ ||A|| ||x||. Note que existe um vetor x0 ∈ Rn tal que ||Ax0|| = ||A|| ||x0||. 2.3 Tipos especiais de matrizes quadradas Definição 2.8 (Matriz Diagonal). Uma matriz quadrada D = (dij) é diagonal se os elementos fora da diagonal principal são todos nulos, ou seja, dij = 0 para i 6= j. Denota-se freqüentemente tais matrizes por D = diag(d11, d22, . . . , dnn). Como exemplo de matrizes diagonais temos A = diag(3,−7, 2) e B = diag(4, 5), que também podemos escrever da forma   3 0 0 0 −7 0 0 0 2   e   4 −5   . A matriz I = diag(1, 1, . . . , 1) é chamada matriz identidade. Observamos que a soma, o produto por escalar e o produto de matrizes diagonais são também diagonais. Observe que a operação de multiplicação entre duas matrizes diagonais n×n quaisquer é comutativa. Definição 2.9 (Matriz Triangular Superior). Uma matriz quadrada A = (aij) é chamada triangular superior se todos os elementos abaixo da diagonal principal são nulos, isto é, se aij = 0 para i > j. Vejamos dois exemplos de matrizes triangulares superiores:   1 0 6 6 3 −9   e   11 −4 1   . 8 Tipos especiais de matrizes quadradas Podemos, agora, apresentar o seguinte teorema para matrizes triangulares superiores: Teorema 2.1. Sejam A e B matrizes triangulares superiores. Então, (i) C = A + B é triangular superior, com cii = aii + bii, i = 1, . . . , n. (ii) C = kA é triangular superior, com cii = kaii, i = 1, . . . , n. (iii) C = AB é triangular superior, com cii = aiibii, i = 1, . . . , n. (iv) Para qualquer polinômio f(x), a matriz B = f(A) é triangular superior e sua diagonal é dada por bii = f(aii) (v) A é inversível se, e somente se, aii 6= 0, i = 1, . . . , n. Demonstração: Veja Lipschutz [19], página 135. Analogamente, uma matriz triangular inferior é uma matriz quadrada cujos elementos acima da diagonal principal são todos nulos, isto é, aij = 0 para i < j. Para tais matrizes, vale teorema análogo ao Teorema 2.1. Note que os zeros podem ser omitidos das matrizes diagonais e, também, das matrizes triangulares. Definição 2.10 (Matriz em Banda). Se A é uma matriz quadrada, tal que aij = 0 para |i − j| > m, onde 2m + 1 < n, então dizemos que A é uma matriz em banda, com largura da banda igual a m. Definição 2.11 (Transposta de uma matriz). A transposta de uma matriz A, denotada por At, é a matriz obtida escrevendo-se as linhas de A, em ordem, como colunas, isto é,   a11 a12 . . . a1n a21 a22 . . . a2n ... ... ... ... am1 am2 . . . amn   t =   a11 a21 . . . am1 a12 a22 . . . am2 ... ... ... ... a1n a2n . . . amn   . Note que a transposta de um vetor coluna é um vetor linha e vice-versa. A operação de tranposição em matrizes satisfaz às seguintes propriedades: Teorema 2.2. 1. (A + B)t = At + Bt; 2. (At)t = A; 3. (kA)t = k(A)t, ∀ k; 9 Tipos especiais de matrizes quadradas 4. (AB)t = BtAt. Demonstração: Veja Lipschutz [19], página 111. Definição 2.12 (Matriz de Hessenberg). Dizemos que A é uma matriz de Hessenberg se aij = 0 para i > j + 1. A matriz abaixo é um exemplo de matriz de Hessenberg,   11 2 −2 1 1 4 3 −4 0 10 7 8 0 0 −1 7   . Definição 2.13 (Matriz Simétrica). Uma matriz quadrada A = (aij) é simétrica se At = A. Equivalentemente, A é simétrica se os elementos simétricos com relação à diagonal principal são iguais, isto é, aij = aji, para todo j e i. A matriz A é uma matriz simétrica, onde A =   2 −3 5 −3 6 7 5 7 −8   . Definição 2.14 (Matriz Anti-Simétrica). Uma matriz quadrada A é anti-simétrica se At = −A. Equivalentemente, A é anti-simétrica se aij = −aji, para todo j e i. Obviamente, os elementos diagonais de uma matriz anti-simétrica devem ser zero, pois bii = −bii ⇒ bii = 0. Por exemplo, B =   0 3 −6 −3 0 1 6 −1 0   . Se A e B são matrizes simétricas, então A + B e kA são simétricas. Porém, AB não é necessariamente simétrica. Teorema 2.3. Seja A uma matriz quadrada. Então, • A + At é simétrica; 10 Tipos especiais de matrizes quadradas • A−At é anti-simétrica e • A = B + C para alguma matriz simétrica B e alguma matriz anti-simétrica C. Demonstração: Veja Lipschutz [19], página 136. Definição 2.15. Uma matriz real A é ortogonal se AAt = AtA = I. Note que A−1 = At. Como exemplo de uma matriz ortogonal, temos A = 1 9   1 8 −4 4 −4 7 8 1 4   . Teorema 2.4. Seja A uma matriz real. Então, as seguintes afirmações são equivalentes: (i) A é ortogonal; (ii) As linhas de A formam um conjunto ortonormal e (iii) As colunas de A formam um conjunto ortonormal. Demonstração: Veja Lipschutz [19], página 137. Observação 2.1. Usaremos a notação A ∈ Rn×n(ou A ∈ Cn×n) quando A for uma matriz real (ou complexa) de ordem n. Para n = 2, temos o seguinte resultado, também demonstrado em Lipschutz [19], página 137 Teorema 2.5. Toda matriz ortogonal 2× 2 tem a forma   cosθ senθ −senθ cosθ   para algum θ real. Seja A = (aij) ∈ Cn×n. Definimos a matriz Ā = (bij) = (aij), como sendo a matriz conjugada de A. As operações de transposição e conjugação comutam para qualquer matriz complexa A, isto é, (A)t = (At). Usa-se a notação especial AH , A∗, para a transposta conjugada de A. Note que se A é real, então AH = At. Vejamos, agora, algumas matrizes complexas especiais. Definição 2.16. Seja A ∈ Cn×n 11 Autovalores e autovetores • Se AH = A, então A é uma matriz hermitiana; • Se AH = −A, então A é uma matriz anti-hermitiana; • Se AH = A−1, então A é uma matriz unitária. Definição 2.17 (Matriz Definida-Positiva). Uma matriz quadrada A, de ordem n, é definida- positiva se xtAx > 0 para qualquer vetor x 6≡ 0. 2.4 Autovalores e autovetores Sejam K = R, ou K = C, e A ∈ Kn×n, isto é, A =   a11 a12 · · · a1n a21 a22 · · · a2n ... ... . . . ... an1 an2 · · · ann   , com aij ∈ K. Definição 2.18. A matriz A− tI =   a11 − t a12 · · · a1n a21 a22 − t · · · a2n ... ... . . . ... an1 an2 · · · ann − t   , onde I é a matriz identidade de ordem n e t é uma incógnita, é chamada matriz característica de A. O determinante p(t) = det(A− tI) é um polinômio em t chamado polinômio característico de A. Dizemos, também, que p(t) = 0 é a equação característica de A. Desenvolvendo p(t) = det(A− tI), obtemos p(t) = (−1)ntn + (−1)n−1bn−1t n−1 + . . . + b0, onde bn−1 é o traço da matriz A e b0 = det(A). 12 Autovalores e autovetores Diz-se que um escalar λ ∈ K é um autovalor de A, se existe um vetor não-nulo v ∈ Kn para o qual Av = λv. Todo vetor que satisfaz a essa relação é, então, chamado um autovetor de A com relação ao autovalor λ. Observe que todo vetor múltiplo de v por um escalar também é um autovetor com relação a λ, pois A(kv) = k(Av) = k(λv) = λ(kv), onde k é um escalar. O produto Av pode variar o módulo e o sentido do vetor, mas nunca sua direção. O conjunto Eλ de todos os autovetores associados a λ é um subespaço de Kn, geralmente conhecido como autoespaço de λ. O conjunto de todas as raízes é chamado espectro de A e vamos denotar por λ(A). De acordo com a definição acima, os autovalores de At são os valores para os quais det(At − λI) = 0. Como o determinante de uma matriz é igual ao determinante de sua transposta, os autovalores de At são os mesmos de A. Os autovetores, em geral, são diferentes. Assim, seja λi o autovalor de At com relação ao autovetor yi, ou seja, Atyi = λiyi. (2.1) Agora, reescrevendo (2.1) da seguinte forma, temos yt iA = λiy t i . O vetor yi é chamado de autovetor à esquerda de A com relação ao autovalor λi. Teorema 2.6. Seja A ∈ Kn×n. Então, as afirmações seguintes são equivalentes: i) Um escalar λ ∈ K é um autovalor de A; ii) A matriz M = λI −A é singular; iii) O escalar λ é raiz do polinômio característica p(t) de A. 13 Autovalores e autovetores Demonstração: Veja Lipschutz [19], página 405. Utilizando o teorema acima, se p(λ) é o polinômio característico de uma matriz A de ordem n, então podemos decompô-lo da seguinte maneira: p(λ) = (λ− λ1)n1(λ− λ2)n2 . . . (λ− λl)nl , onde l ≤ n, n1 +n2 + . . .+nl = n e λi, i = 1, . . . , l, são os autovalores de A. A dimensão de Eλi é chamada de multiplicidade geométrica de λi e ni é chamado de multiplicidade algébrica de λi. Agora, seja a matriz A + σI, onde σ é um escalar qualquer e λ é um autovalor de A com relação ao autovetor v. Assim, (A + σI)v = Av + σv = λv + σv = (λ + σ)v. Logo, (λ + σ) é um autovalor de A + σI com relação ao autovetor v. Definição 2.19. Uma matriz B é similar, ou semelhante, a uma matriz A se existe uma matriz não-singular P tal que B = P−1AP . Uma matriz A é diagonalizável se existe uma matriz não-singular P tal que D = P−1AP é uma matriz diagonal, isto é, se A é similar a uma matriz diagonal D. O teorema seguinte caracteriza tais matrizes. Teorema 2.7. Uma matriz quadrada A de ordem n é similar a uma matriz diagonal D se, e somente se, A tem n autovetores linearmente independentes. Neste caso, os elementos diago- nais de D são os autovalores de A e D = P−1AP , onde P é a matriz cujas colunas são os correspondentes autovetores. Demonstração: Veja Lipschutz [19], página 406. Examinaremos, neste trabalho, o problema de obtenção dos autovalores de uma matriz quadrada A, de ordem n, através de transformações de similaridade que possibilitem a criação de uma seqüência de matrizes semelhantes, que tendam a uma matriz triangular superior A∞. Agora, os resultados abaixo nos garantem que, dado uma matriz A, podemos obter uma matriz triangular superior com os mesmos autovalores de A. Teorema 2.8 (Teorema de Schur). Seja A ∈ Cn×n. Então, existem uma matriz unitária 14 Autovalores e autovetores U ∈ Cn×n e uma matriz triangular T ∈ Cn×n tais que T = U∗AU . T é da forma T =   t11 t12 · · · t1n 0 t22 · · · t2n ... ... . . . ... 0 0 · · · tnn   , onde os elementos diagonais tii são autovalores de A. Teorema 2.9 (Teorema de Schur - para matrizes reais). Seja A ∈ Rn×n. Então, existem uma matriz unitária U ∈ Cn×n e uma matriz bloco triangular T ∈ Rn×n tais que T = U∗AU . T é da forma T =   T11 T12 · · · T1n 0 T22 · · · T2n ... ... . . . ... 0 0 · · · Tnn   , onde os blocos diagonais Tii são 1×1 ou 2×2. Cada bloco 1×1 é um autovalor real de A e cada bloco 2 × 2 tem como autovalores um par de autovalores complexos conjugados de A. A matriz U é chamada de matriz dos autovetores e U t a matriz dos autovetores à esquerda. A demonstração dos dois Teoremas acima podem ser encontradas em Golub [12], páginas 192 e 193. Teorema 2.10 (Teorema Espectral). Seja A ∈ Cn×n hermitiana. Então, existem uma matriz unitária U ∈ Cn×n e uma matriz diagonal D ∈ Rn×n tais que D = U∗AU . As colunas de U são os autovetores de A e os elementos da diagonal de D são os autovalores de A. Demonstração: Seja λ um autovalor de A com relação ao autovetor v. Assim, temos que Av = λv ⇒ (Av)H = (λv)H ⇒ vHAH = λ̄vH ⇒ vHA = λ̄vH ⇒ ⇒ vHAv = λ̄vHv ⇒ vHλv = λ̄vHv ⇒ λ||v||22 = λ̄||v||22 ⇒ λ = λ̄ Logo, todos os autovalores de A são reais. Agora, como A é hermitiana, então de acordo com o Teorema de Schur 2.8, T = UHAU = UHAHU = (UHAU)H = TH . Portanto, T só pode ser uma matriz diagonal. 15 Autovalores e autovetores Teorema 2.11 (Decomposição de Jordan). Seja A ∈ Cn×n. Então, existe uma matriz não-singular X ∈ Cn×n tal que X−1AX = JA = diag(Jm1(λ1), . . . , Jml (λl)), onde Jmi(λi) =   λi 1 . . . . . . λi 1 λi   é quadrada e de ordem mi e m1 + m2 + . . . + ml = n. Demonstração: Veja Golub [12], página 197. A matriz JA é chamada forma canônica de Jordan de A e as matrizes Jmi são chamadas de blocos de Jordan. O número e a dimensão dos blocos de Jordan associados a cada autovalor são únicos, porém a ordem em que eles aparecem na diagonal principal pode se alterar. Sejam A uma matriz com seis autovalores, sendo λ1 = . . . = λ5 6= λ6, e C sua forma canônica de Jordan, onde C =   J2(λ1) J2(λ1) J1(λ1) J1(λ6)   . Agora, em geral, a matriz (C − λI) é a soma, direta, das matrizes da forma Jr(λi − λ). Os determinantes das sub-matrizes da forma canônica (C −λI) são chamados divisores elementares da matriz A. Assim, os divisores elementares de qualquer matriz similar à C são (λ1 − λ)2, (λ1 − λ)2, (λ1 − λ) e (λ6 − λ). Podemos observar que o produto dos divisores elementares de uma matriz é igual ao seu polinômio característico. Quando a forma canônica de Jordan é diagonal, então seus divisores elementares são (λ1 − λ), (λ2 − λ), . . . , (λn − λ), ou seja, são lineares. Uma matriz com autovalores distintos sempre possui divisores elementares lineares, mas, como vimos, se uma matriz possui dois ou mais autovalores iguais, então não podemos afirmar se ela possui, ou não, divisores elementares lineares. 16 Autovalores e autovetores Se uma matriz possui um ou mais divisores elementares não-lineares, então a forma canônica de Jordan desta matriz possui pelo menos um bloco com dimensão maior que 1 e, além disso, m autovetores linearmente independente, onde m < n. Tais matrizes são chamadas defectivas. Os Teoremas de Gershgorin, das Matrizes Correspondentes e de Sylvester são de funda- mental importância para os Capítulos 3 e 4, por isso são apresentados a seguir. Teorema 2.12 (Teorema de Gerschgorin). Sejam A ∈ Cn×n, com A = (aij), e p′i = n∑ j=1 |aij |, onde p′i representa o somatório de j = 1, . . . , n, com i 6= j. Então, todos os autovalores de A estão em pelo menos um dos discos {z : |z − aij | ≤ p′i}, i = 1, . . . n, no plano complexo. Além disso, um conjunto de m discos que não possuem pontos em comum com os n−m discos restantes, possui exatamente m autovalores de A. Demonstração: Veja Peter [23], página 371 Teorema 2.13 (Teorema das Matrizes Correspondentes). Sejam A,B e C matrizes tais que C = AB, onde C é uma matriz quadrada de ordem n. Assim, o determinante de C é a soma dos produtos de todos menores de ordem n de A com o correspondente menor de mesma ordem de B. Demonstração: Veja Gantmacher [11], página 9. Teorema 2.14 (Teorema de Sylvester). Sejam A e B matrizes quadradas de mesma ordem. O determinante do produto AB é dado pelo somatório ∑ |A′B′|, onde A′ e B′ são iguais às matrizes A e B, respectivamente, porém com k colunas trocadas entre si. A troca das k colunas é feita fixando-se as k colunas de A que serão trocadas; dentre essas k colunas, a troca é feita sendo a primeira com a primeira de B, a segunda com a segunda de B, etc. Essas trocas são realizadas até que as k colunas de A tenham sido trocadas com todos os grupos possíveis de k colunas de B. Demonstração: Veja Muir [21], página 120. 17 Reflexões de Householder 2.5 Reflexões de Householder Seja w ∈ Rn tal que ||w||2 = 1. Uma reflexão de Householder é uma matriz da forma H = I − 2wwt. Por exemplo, se wt = 1 5(3, 4)t, então H = 1 25   7 −24 −24 −7   . As matrizes de Householder são simétricas e ortogonais. De fato: • Ht = (I − 2wwt)t = I − 2wwt = H. • HH = (I − 2wwt)(I − 2wwt) = I − 4wwt + 4wwtwwt. Como ||w||2 = wtw = 1, então HH = I e, portanto, H = H−1. Multiplicar uma matriz de Householder por um vetor x é equivalente a refletir o vetor x através do plano perpendicular a w. Vejamos: Hx = (I − 2wwt)x = x− 2wwtx = x− 2(wtx)w. Mas, (wtx)w é a projeção de x em w. Figura 2.1: Hx é a reflexão de x através do plano perpendicular a w. Devido a esta propriedade de reflexão, as matrizes de Householder são chamadas reflexões elementares. Vejamos outra propriedade das matrizes de Householder. Sejam x e y vetores não nulos tais que ||x||2 = ||y||2 = 1. Seja 18 Reflexões de Householder w = 1 ||x− y||2 (x− y). (2.2) Então, (I − 2wwt)x = y. Utilizando (2.2), podemos construir uma matriz de Householder que anula algumas compo- nentes de um vetor x, deixando todas as demais intactas, exceto uma. Assim, dado k ∈ Z, k ≥ 1, vamos construir a matriz de Householder tal que (Hx)i =    xi, i = 1, . . . , k − 1, 0, i = k + 1, . . . , n. Considerando (2.2), escolhemos y da seguinte forma: y1 = x1, y2 = x2, . . . , yk−1 = xk−1, yk = ± √ x2 k + x2 k+1 + ... + x2 n e yk+1 = · · · = yn = 0. Assim, ||x||2 = ||y||2 e, se escolhermos w como em (2.2), teremos Hx = y. Observe que existem duas possibilidades de sinais para o k-ésimo elemento de y. Se xk e yk possuem sinais opostos, então o erro cometido será menor. Por esta razão, podemos considerar w = 1√ 2s(s + |xk|)   0 ... 0 xk + sinal(xk)s xk+1 ... xn   , (2.3) onde s = √ x2 k + x2 k+1 + ... + x2 n. Exemplo 2.1. Considere o vetor x = (7, 4, 3)t. Vamos construir a matriz H que anula a terceira coordenada e mantém a primeira intacta. Temos k = 2 e, então, s = 5. Logo, w = 1√ 90   0 9 3   e H = −1 5   5 0 0 0 −4 −3 0 −3 4   . 19 Reflexões de Householder Portanto, obtemos Hx = y = (7,−5, 0)t. Podemos, também, utilizar as reflexões de Householder para obter duas matrizes similares, sendo que a primeira é dada e a segunda é de Hessenberg. Isto é, se A é uma matriz quadrada, então existe uma matriz H tal que H−1AH é uma matriz de Hessenberg, onde a matriz H pode ser obtida através das reflexões de Householder. Vejamos um exemplo com a matriz A =   450 75 −525 150 75 253 380 −79 150 5 325 −215 150 −604 160 322   . Seja H1 a matriz que anula o terceiro e quarto elementos da primeira coluna. Assim, utilizando (2.3), temos H1 =   1 0 0 0 0 −1 3 −2 3 −2 3 0 −2 3 2 3 −1 3 0 −2 3 −1 3 2 3   . Fazendo o produto H1AH1, obtemos A1 =   450 225 −450 225 −225 225 −495 −90 0 180 −18 −126 0 135 99 693   . Agora, resta apenas anular o quarto elemento da segunda coluna. Para isso, temos H2 =   1 0 0 0 0 1 0 0 0 0 −0.8 −0.6 0 0 −0.6 0.8   , que, ao fazermos o produto H2A1H2, obtemos A2 =   450 225 225 450 −225 225 −450 225 0 −225 225 −225 0 0 −450 450   . 20 Rotações de Givens Podemos ver que A2 = H2H1AH1H2 = H−1AH, onde H = H1H2. O algoritmo para reduzir uma matriz quadrada A, com dimensão n, à forma de Hessen- berg ocorre de modo análogo. No primeiro passo, construímos a matriz H1 que anula todos os elementos da primeira coluna, exceto os dois primeiros, e calculamos A1 = H1AH1. No k-ésimo passo, construímos a matriz Hk que anula todos os elementos da k-ésima coluna, exceto os k +1 primeiros, e calculamos Ak = HkAk−1Hk. No passo n− 2, obtemos a matriz An−2 na forma de Hessenberg, onde An−2 = Hn−2 . . .H1AH1 . . . Hn−2. Finalmente, H = H1 . . .Hn−2. 2.6 Rotações de Givens Uma rotação de Givens é uma matriz da forma Gij =   1 . . . 1 c 0 . . . 0 −s 0 1 0 ... . . . ... 0 1 0 s 0 . . . 0 c 1 . . . 1   , onde a linha com −s é a linha j, a linha com s é a linha i e c2 + s2 = 1. Note que a intersecção das linhas i e j com as colunas i e j de Gij forma   c −s s c   . 21 Rotações de Givens Esta matriz é ortogonal, pois   c −s s c     c s −s c   =   c2 + s2 0 0 c2 + s2   =   1 0 0 1   . Podemos concluir, assim, que as matrizes Gij são ortogonais. Essas matrizes são chamadas de rotações pois, se x e y satisfazem   y1 y2   =   cosφ −senφ senφ cosφ     x1 x2   , então y é o vetor x girado do ângulo φ. Figura 2.2: Rotação de Givens Assim, multiplicar um vetor por Gij é equivalente a girar duas componentes do vetor do ângulo φ = arctan ( s c ) , deixando as demais componentes intactas. Seja x um vetor com duas componentes e definimos c = x1√ x2 1 + x2 2 e s = − x2√ x2 1 + x2 2 . (2.4) Observe que c2 + s2 = 1 e   c −s s c     x1 x2   =   √ x2 1 + x2 2 0   . (2.5) Em um trabalho recente, Matioli, juntamente com Biloti e Jin-Yun Yuan, [20] apresentou uma nova transformação, a rotação de Givens generalizada. Essa nova transformação reduz uma matriz quadrada qualquer à forma tri-diagonal sem utilizar as matrizes de Hessenberg. 22 Eliminação de Gauss e Decomposição LU Além disso, é mais rápida, computacionalmente falando, do que as que utilizam as matrizes de Hessenberg. Observação 2.2. As rotações de Givens também podem ser utilizadas para construir uma matriz de Hessenberg similar à uma matriz quadrada qualquer. 2.7 Eliminação de Gauss e Decomposição LU Métodos para solução de sistemas de equações lineares são divididos principalmente em dois grupos: 1) Métodos Exatos: são aqueles que forneceriam a solução exata, não fossem os erros de arredondamento, com um número finito de operações e 2) Métodos Iterativos: são aqueles que permitem obter as raízes de um sistema, com uma dada precisão, através de um processo infinito convergente. A Eliminação de Gauss é um método exato e, ainda hoje, é o método mais utilizado para solução de sistemas de equações lineares. Vejamos um exemplo de sua aplicação no seguinte sistema com 3 equações e 3 variáveis:    2x1 + 4x2 + 2x3 = 4 4x1 + 7x2 + 7x3 = 13 −2x1 − 7x2 + 5x3 = 7 . (2.6) Primeiro, eliminamos a variável x1 da segunda e terceira equações subtraindo duas vezes a primeira equação da segunda e subtraindo -1 vezes a primeira equação da terceira. Obtemos, então um novo sistema linear    2x1 + 4x2 + 2x3 = 4 − x2 + 3x3 = 5 − 3x2 + 7x3 = 11 . (2.7) Em seguida, eliminamos a variável x2 da terceira equação subtraindo três vezes a segunda equação da terceira. Chegamos, assim, a uma estrutura triangular:    2x1 + 4x2 + 2x3 = 4 − x2 + 3x3 = 5 − 2x3 = −4 . (2.8) 23 Eliminação de Gauss e Decomposição LU Por último, basta calcular os valores de x3, x2 e x1. Obtemos, então, x3 = −4 −2 = 2, x2 = 5− 3x3 −1 = 5− 3× 2 −1 = 1, x1 = 4− 4x2 − 2x3 2 = 4− 4× 1− 2× 2 2 = −2. Agora, utilizando notação matricial, podemos reescrever o sistema (2.6) da seguinte forma:   2 4 2 4 7 7 −2 −7 5     x1 x2 x3   =   4 13 7   , onde a matriz à esquerda é chamada matriz dos coeficientes do sistema, que será denotada por A. Sejam as matrizes U e L dadas por: U =   2 4 2 0 −1 3 0 0 −2   e L =   1 0 0 2 1 0 −1 3 1   . A matriz U contém os coeficientes do sistema eliminado (2.8) e a matriz L representa os passos para a eliminação. Durante a eliminação, subtraímos 2 vezes a primeira equação da segunda; −1 vez a primeira equação da terceira; 3 vezes a segunda equação da terceira. Os fatores de eliminação 2, −1 e 3 são chamados de pivot. A matriz A é igual ao produto da matriz L por U :   1 0 0 2 1 0 −1 3 1     2 4 2 0 −1 3 0 0 −2   =   2 4 2 4 7 7 −2 −7 5   . Utilizando esse exemplo, podemos concluir que: A Eliminação de Gauss de uma matriz A é equivalente a decompô-la no produto de uma matriz triangular inferior, com os elementos da diagonal principal iguais a 1, e uma matriz triangular superior. Embora essa afirmação seja verdadeira, apenas um grupo de matrizes possuem decom- posição LU. Vejamos o seguinte resultado que caracteriza tais matrizes. 24 Eliminação de Gauss e Decomposição LU Teorema 2.15. Sejam A uma matriz real quadrada de ordem n e Ak a matriz menor principal constituída das k primeiras linhas e k primeiras colunas de A. Suponha que det(Ak) 6= 0 para k = 1, 2, . . . , n − 1. Então, existe uma única matriz triangular inferior L, com lii = 1 para i = 1, . . . , n, e uma única matriz triangular superior U tal que A = LU . Além disso, det(A) = u11u22 . . . unn. Demonstração: Para demonstrar esse teorema, usaremos indução sobre n. Se n = 1, temos que A = (a11) = 1(u11) = LU unicamente, e det(A) = u11. Suponhamos que o teorema seja verdadeiro para n = k − 1. Para n = k, particionamos A da seguinte forma: A =   Ak−1 x y akk   , onde y é um vetor linha com k − 1 componentes e x é um vetor coluna com, também, k − 1 componentes. Escrevendo L =   Lk−1 0 m 1   e U =   Uk−1 p 0 ukk   , onde m e p são vetores com dimensões apropriadas, temos que LU =   Lk−1Uk−1 Lk−1p mUk−1 mp + ukk   . Agora, pela hipótese de indução, Lk−1 e Uk−1 são unicamente determinados e Lk−1Uk−1 = Ak−1. Além disso, Lk−1 e Uk−1 são não-singulares, pois Ak−1 também é não-singular. Assim, LU = A é equivalente a Lk−1p = x, mUk−1 = y e mp + ukk = akk, ou seja, p = L−1 k−1x, m = U−1 k−1y e ukk = akk −mp. Então, p, m e ukk são determinados univocamente nesta ordem e, assim, L e U também. Finalmente, det(A) = det(L)det(U) = 1det(Uk−1)ukk = u11u22 . . . unn. Para matrizes que não possuem decomposição LU, podemos fazer uma modificação no algoritmo de forma que seja possível decompô-las no produto LU. Veremos essa alteração na seção 3.1.4. 25 Decomposição QR 2.8 Decomposição QR De maneira similar à decomposição LU, podemos decompor uma matriz A no produto QR, onde Q é uma matriz ortogonal e R é triangular superior. A decomposição QR de uma matriz A pode ser calculada utilizando-se uma sequência de matrizes de Householder para reduzir A a uma matriz triangular superior. Esta matriz triangular é a matriz R e o produto das matrizes de Householder é a matriz Q. Vejamos o exemplo a seguir. Exemplo 2.2. Considere a matriz A =   10 9 18 20 −15 −15 20 −12 51  . A primeira coluna de A é (10, 20, 20)t. Assim, usando (2.3) com k = 1, temos H1 = −1 3   1 2 2 2 −2 1 2 1 2   . Logo, A1 = H1A = −1 3   1 2 2 2 −2 1 2 1 2     10 9 18 20 −15 −15 20 −12 51   =   −30 15 −30 0 −12 −39 0 9 27   . Agora, queremos anular o terceiro elemento da segunda coluna de A1. Novamente usando (2.3), mas agora com k = 2, obtemos H2 =   1 0 0 0 −0.8 −0.6 2 −0.6 0.8   . Assim, A2 = H2A1 =   1 0 0 0 −0.8 −0.6 0 −0.6 0.8     −30 15 −30 0 −12 −39 0 9 27   =   −30 15 −30 0 15 15 0 0 45   . Observe que a primeira coluna de A2 é igual à primeira de A1. Além disso, A2 é triangular superior e, então, R = A2 = H2A1 = H2H1A. Portanto, A = H1H2R = QR, onde Q é ortogonal. A decomposição de uma matriz A qualquer, com dimensão n, ocorre da mesma forma. Primeiro pré-multiplicamos A pela matriz de Householder H1 que anula todos os elementos da 26 Decomposição QR primeira coluna, exceto o primeiro, e obtemos A1 = H1A. Em seguida, calculamos A2 = H2A1 pré-multiplicando A1 por H2, que anula todos os elementos, exceto os dois primeiros, da segunda coluna de A1. Após (n− 1) passos, teremos R = An−1 = Hn−1An−2 = · · · = Hn−1Hn−2...H1A. Fazendo Q = H1...Hn−2Hn−1, obtemos A = QR. A decomposição QR de uma matriz também pode ser feita através das rotações de Givens. Utilizando (2.4), podemos construir uma rotação de Givens que anula um elemento de uma matriz dada. Para decompormos A no produto QR, multiplicamos A por uma seqüência de rotações de Givens que anula os elementos abaixo da diagonal principal, obtendo, assim, a matriz R triangular superior. Finalmente, Q é o produto das rotações. Vejamos um exemplo. Exemplo 2.3. Seja A =   90 −153 114 120 −79 −223 200 −40 395  . Primeiramente, vamos anular o elemento a21 = 120. Utilizando (2.4) com x1 = 90 e x2 = 120, obtemos G21 =   90 150 120 150 0 −120 150 90 150 0 0 0 1   =   3 5 4 5 0 −4 5 3 5 0 0 0 1   , que anula a21. Assim, A21 = G21A =   150 −155 −110 0 75 −225 200 −40 395   = 5   30 −31 −22 0 15 −45 40 −8 79   . Novamente utilizando (2.4), mas agora com x1 = 30 e x2 = 40, temos G31 =   0, 6 0 0, 8 0 1 0 −0, 8 0 0, 6   , que anula a31. De fato, A31 = G31A21 = 25   10 −5 −10 0 3 −9 0 4 13   . 27 Decomposição QR Só nos resta agora anular o elemento a32. Para isso, fazemos G32 =   1 0 0 0 0, 6 0, 8 0 −0, 8 0, 6   e R = G32A31 = 125   2 −1 2 0 1 1 0 0 3   . Obtivemos então, R = G32A31 = G32G31A21 = G32G31G21A e, portanto, A = QR, onde Q = (G32 G31 G21)t. Para uma matriz qualquer, anulamos todos os elementos da primeira coluna, exceto o primeiro, com a seqüência de rotações G21, G31, . . . , Gn1. Em seguida, anulamos todos os ele- mentos da segunda coluna, exceto os dois primeiros, com a sequência de rotaçõesG32, G42, . . . Gn2. O último passo é anular o elemento an,n−1 com a rotação Gn,n−1. Finalmente, Q é o produto das matrizes transpostas dessas rotações, nesta mesma ordem. Podemos, agora, demonstrar o seguinte teorema. Teorema 2.16. Seja A ∈ Cn×n. Então, existe uma matriz unitária Q tal que A = QR, onde R é uma matriz triangular superior com os elementos da diagonal principal reais e não negativos. Além disso, Q é única se A é não-singular. Demonstração: A matriz Q pode ser determinada através das reflexões de Householder ou das rotações de Givens, como vimos. Agora, seja A uma matriz tal que A = Q1R1 = Q2R2. (2.9) Se A é não-singular, então R1 e R2 também são. Logo, R1R −1 2 = Q∗ 1Q2 ⇒ (R1R −1 2 )−1 = (R1R −1 2 )∗. Assim, R1R −1 2 é uma matriz diagonal. Considerando a diagonal principal dessa matriz, seus elementos são dados por r (2) ii r (1) ii = r̄ (1) ii r̄ (2) ii . Mas, como r (1) ii e r (2) ii são reais e positivos, r (1) ii = r (2) ii . Logo, R1R −1 2 = I ⇒ R1 = R2. Substituindo em (2.9), temos A = Q1R1 = Q2R1 ⇒ Q1 = Q2. 28 Decomposição de Cholesky 2.9 Decomposição de Cholesky A decomposição de Cholesky consiste em decompor uma matriz A simétrica e definida-positiva no produto RtR, onde R é uma matriz triangular superior com os elementos na diagonal principal positivos. Vamos demonstrar esse resultado utilizando indução sobre a dimensão de A. No caso em que k = 1, onde k é a dimensão de A, temos A = (a) = (r)(r) = RtR ⇒ r2 = a. Mas, como a é positivo, podemos escolher r de forma que r > 0. Agora, vamos supor válido para as matrizes de ordem k, para 1 ≤ k ≤ n − 1. Assim, se A = An, onde n é a dimensão de A, é definida-positiva, podemos escrever An =   An−1 b bt ann   , onde An−1 é a matriz menor principal de ordem (n− 1) formada pelas primeiras n− 1 linhas e colunas, que, por hipótese, é simétrica e definida-positiva. Logo, existe uma matriz triangular superior Rn−1 tal que, Rt n−1Rn−1 = An−1. É óbvio que Rn−1 é não-singular, pois [det(Rn−1)]2 = det(An−1). Assim, existe um vetor c tal que Rt n−1c = b. Agora, para qualquer valor de x, temos que   Rt n−1 0 ct x     Rn−1 c 0 x   =   An−1 b bt ctc + x2   . (2.10) Portanto, se escolhermos x de forma que ctc + x2 = ann, essa é uma decomposição triangular para An. Para mostrar que x é real, basta calcular o determinate em (2.10). Assim, [det(Rn−1)]2x2 = det(An) > 0, pois An é definida-positiva. Logo, podemos escolher x > 0. 29 Capítulo 3 Algoritmo LR de Rutishauser Neste capítulo, nosso objetivo é fazer um estudo sobre o Algoritmo LR para o cálculo de auto- valores de matrizes. Apresentado por Rutishauser no artigo [24], e por isso leva seu nome, esse algoritmo usa transformações não unitárias para reduzir uma matriz à forma triangular superior. Além do trabalho [24], utilizamos, também, o texto de Wilkinson [30] como referência para o estudo que apresentaremos a seguir. 3.1 Algoritmo LR de Rutishauser Em 1958, Rutishauser descreveu um processo iterativo para encontrar os autovalores de uma matriz A. O método consiste em formar uma seqüência de matrizes As da seguinte forma: A1 = A, As = LsRs, s = 1, 2, . . . , (3.1) As+1 = RsLs, onde Ls é uma matriz triangular inferior com todos os elementos da diagonal principal iguais a 1 e Rs é uma matriz triangular superior. Este método de decomposição é o conhecido método LU, como vimos. De (3.1), podemos escrever Rs = L−1 s As. Substituindo na expressão para As+1, obtemos As+1 = L−1 s AsLs, s = 1, 2, . . . Logo, as matrizes As+1 e As, s = 1, 2, . . . , são similares. Como conseqüência, podemos escrever As+1 = L−1 s L−1 s−1 . . . L−1 1 A1L1 . . . Ls−1Ls, s = 1, 2, . . . , (3.2) 30 Algoritmo LR de Rutishauser o que demonstra que todas as matrizes As, s = 2, 3, . . . , são similares a A1. Rutishauser mostrou que, sob certas condições, Ls → I e Rs → A∞ =   λ1 × × × 0 λ2 × × . . . 0 0 0 λn   , quando s →∞. (3.3) Para demonstrarmos este resultado, vamos reescrever (3.2) da seguinte forma: L1 ... Ls−2Ls−1As = A1L1 ... Ls−2Ls−1, s = 2, 3, . . . . (3.4) Agora, as matrizes Ts e Vs, definidas por Ts = L1L2 ... Ls e Vs = RsRs−1 ... R1, (3.5) são, respectivamente, triangular inferior com 1 na diagonal principal e triangular superior. Como TsVs = L1L2 ... Ls−1LsRsRs−1 ... R1 = L1L2 ... Ls−1AsRs−1 ... R1, de (3.4), temos TsVs = A1L1L2 ... Ls−1Rs−1 ... R1, = A1Ts−1Vs−1 = A2 1Ts−2Vs−2, s = 3, 4, . . . . Continuando dessa forma, obtemos TsVs = As 1, s = 1, 2, . . . , (3.6) ou seja, TsVs é a decomposição LU de As 1. 3.1.1 Prova da convergência de As A equação (3.6) é o resultado fundamental para esta análise. Vamos utilizá-la para demonstrar que, se os autovalores de A1 satisfazem à relação |λ1| > |λ2| > ... > |λn|, 31 Algoritmo LR de Rutishauser então, em geral, o resultado (3.3) é verdadeiro. Antes da prova formal, consideremos o caso mais simples de uma matriz de ordem três. Assim, definimos a matriz X dos autovetores à esquerda por X =   x1 y1 z1 x2 y2 z2 x3 y3 z3   e sua inversa Y por X−1 = Y =   a1 b1 c1 a2 b2 c2 a3 b3 c3   . Então, temos F = As 1 = X   λs 1 λs 2 λs 3  Y =   λs 1x1a1 + λs 2y1a2 + λs 3z1a3 ... λs 1x1b1 + λs 2y1b2 + λs 3z1b3 ... λs 1x1c1 + λs 2y1c2 + λs 3z1c3 λs 1x2a1 + λs 2y2a2 + λs 3z2a3 ... λs 1x2b1 + λs 2y2b2 + λs 3z2b3 ... λs 1x2c1 + λs 2y2c2 + λs 3z2c3 λs 1x3a1 + λs 2y3a2 + λs 3z3a3 ... λs 1x3b1 + λs 2y3b2 + λs 3z3b3 ... λs 1x3c1 + λs 2y3c2 + λs 3z3c3   . Agora, TsVs é a decomposição LU de As 1 e, então, os elementos da primeira coluna de Ts são dados por t (s) 11 = 1, t (s) 21 = λs 1x2a1 + λs 2y2a2 + λs 3z2a3 λs 1x1a1 + λs 2y1a2 + λs 3z1a3 , t (s) 31 = λs 1x3a1 + λs 2y3a2 + λs 3z3a3 λs 1x1a1 + λs 2y1a2 + λs 3z1a3 . Podemos ver que, se x1a1 6= 0, t (s) 21 = x2 x1 + O ( λ2 λ1 )s , t (s) 31 = x3 x1 + O ( λ2 λ1 )s e os elementos de Ts tendem, em geral, aos correspondentes elementos obtidos na decomposição LU de X. Para os elementos da segunda coluna de Ts, temos t (s) 22 = 1, t (s) 32 = f11f32 − f12f31 f11f22 − f12f21 = x1y3 − x2y1 x1y2 − x2y1 + O ( λ2 λ1 )s , (3.7) 32 Algoritmo LR de Rutishauser se (a1b2 − a2b1)(x1y2 − x2y1) 6= 0. (3.8) De (3.7), podemos ver que o valor limite de t (s) 32 é igual ao correspondente elemento obtido na decomposição LU de X. Assim, podemos concluir que se X = TU , x1a1 6= 0 e (a1b2−a2b1)(x1y2− x2y1) 6= 0, então Ts → T . Desta forma, com este simples caso mostramos que a matriz Ts tende à matriz L obtida da decomposição LU da matriz X dos autovetores, desde que os menores principais de X e Y sejam diferentes de zero. Mostraremos, agora, que este resultado vale, em geral, para qualquer matriz com autovalores distintos. Assim, escrevendo TsVs = As = B, temos que os elementos t (s) ji são dados pelas equações t (s) ji = det(Bji) det(Bii) , j = 1, . . . , n− 1, i = j, . . . , n, onde n é a dimensão da matriz A. Aqui, Bii denota a matriz menor principal formada pelas primeiras i linhas e i colunas de B e Bji é a matriz obtida de Bii substituindo-se a i-ésima linha pela j-ésima linha de B. Da relação B = As = X   λ1 . . . λn  Y, obtemos Bji =   x11 x12 . . . x1n x21 x22 . . . x2n . . . . . . . . . . . . xi−1,1 xi−1,2 xi−1,n xj1 xj2 . . . xjn   .   λs 1y11 λs 1y12 . . . λs 1y1i λs 2y21 λs 2y22 . . . λs 2y2i . . . . . . . . . . . . λs nyn1 λs nyn2 . . . λs nyni   . (3.9) Do Teorema das Matrizes Correspondentes (Teorema 2.13), det(Bji) é igual à soma dos produtos dos menores principais de ordem i das duas matrizes à direita em (3.9). Assim, podemos escrever t (s) ji = ∑ x (j) p1p2...pi y (i) p1p2...pi(λp1λp2 ...λpi) s ∑ x (i) p1p2...pi y (i) p1p2...pi(λp1λp2 ...λpi)s , (3.10) 33 Algoritmo LR de Rutishauser onde x (j) p1p2...pi é o menor de ordem i consistindo das linhas 1, 2, ..., i − 1,j e colunas p1, p2, ..., pi de X, e y (i) p1p2...pi é o menor de ordem i consistindo das linhas p1, p2, ..., pi e colunas 1, 2, ..., i de Y . Os termos dominantes no numerador e denominador de (3.10) são aqueles associados à (λp1λp2 ...λpi) s se os correspondentes coeficientes são diferentes de zero. O termo apropriado no denominador é det(Xii)det(Yii)(λp1λp2 ...λpi) s, onde Xii e Yii são as matrizes menores principais formadas pelas primeiras i linhas e primeiras i colunas de X e Y , respectivamente. Agora, se det(Xii)det(Yii) 6= 0, então t (s) ji → det(Xji)det(Yii) det(Xii)det(Yii) = det(Xji) det(Xii) , (3.11) o que implica que Ts → T , onde X = TU. De (3.4) e (3.5), temos As = T−1 s−1A1Ts−1 → T−1A1T = UX−1AXU−1 = U   λ1 . . . λn  U−1, (3.12) o que nos mostra que As converge para uma matriz triangular superior com λi na diagonal principal. Utilizando a relação Ls = T−1 s−1Ts e a equação (3.10), podemos mostrar que l (s) ji = O ( λi λj )s , quando s → ∞. Assim, da relação As = LsRs e do fato de Rs convergir para um limite, concluímos que, para i > j, a (s) ij = O ( λi λj )s , quando s →∞. Para obter esses resultados, consideramos as seguintes hipóteses: (i) todos os autovalores possuem módulos distintos; (ii) a decomposição LU de As existe para todo s. Existem matrizes em que isto não é possível, como, por exemplo, a matriz A =   0 1 −3 4   . Esta matriz possui autovalores iguais a 1 e 3, mas não possui decomposição LU. Note que a matriz (A + I) possuiu decomposição LU e, assim, poderíamos trabalhar com ela normalmente; 34 Algoritmo LR de Rutishauser (iii) os menores principais de X e Y são diferentes de zero. Existem algumas matrizes onde isto não ocorre. Os exemplos mais triviais são as matrizes da forma   A1 0 0 A2   , que claramente mantém esta estrutura em cada estágio s. As matrizes A1 e A2 são tratadas independentemente e a matriz limite possui os autovalores de A1 no canto esquerdo superior e os autovalores de A2 no canto direito inferior, independentemente de suas dimensões. Consideremos a seguinte matriz 2× 2:   a ε1 ε2 b   , onde |b| > |a|. Se os εi são iguais a zero, a matriz permanece inalterada pelo algoritmo LR e a matriz limite não possuiu os autovalores a e b na ordem correta. A matriz dos autovetores é   0 1 1 0   e o seu menor principal de ordem 1 é igual a zero. Agora, se os εi são diferentes de zero, mas pequenos, então os autovalores estão próximos de a e b. A matriz limite, agora, possui os autovalores na ordem correta, mas sua convergência pode ser lenta, especialmente se a e b não estão próximos. Neste exemplo, pudemos perceber que as duas metades da matriz podem ser “separadas” quando ε1 = ε2 = 0 e, se os εi estão próximos de zero, então as duas metades da matriz possuem o que chamaremos de "ligação fraca". Como exemplo mais complexo, podemos ter separações iniciais, isto é, quando, para algum i, tivermos det(Xii)det(Yii) = 0. Mas, erros de arredondamentos no processo LR podem causar "ligação fraca". 3.1.2 Matrizes hermitianas definidas-positivas Quando A1 é uma matriz hermitiana definida-positiva, podemos remover as restrições (i), (ii) e (iii). Temos, agora, A1 = X   λ1 . . . λn  XH , onde X é uma matriz unitária e os λi são reais e positivos. Assim, podemos escrever a equação (3.10) como t (s) ji = ∑ x (j) p1p2...pi x̄ (i) p1p2...pi(λp1λp2 ...λpi) s ∑ |x(i) p1p2...pi |2(λp1λp2 ...λpi)s . 35 Algoritmo LR de Rutishauser Considere, agora, a classe dos agrupamentos p1p2...pi para os quais x (i) p1p2...pi 6= 0. Seja q1q2...qi um membro desta classe, para o qual λq1λq2 ...λqi assume o valor máximo. O termo dominante no denominador é aquele em que (λq1λq2 ...λqi) s. A medida em que o numerador é calculado, sabemos, pela definição dos qi, que não pode haver termo de ordem maior do que (λq1λq2 ...λqi) s e, assim, este termo deve ter o valor zero como coeficiente. Logo, t (s) ji tende a zero. Nossa intenção, aqui, é de não fazermos restrições quanto aos menores principais de X e, então, não podemos garantir que o limite Ts é a matriz obtida da decomposição LR de X. Por essa razão, não agimos como em (3.12), mas sim da seguinte forma. Seja T∞ = lim s→∞Ts. Então, Ls = T−1 s−1Ts → T−1∞ T∞ = I. Desta forma, temos As = T−1 s−1A1Ts−1 → T−1 ∞ A1T∞, (3.13) o que implica que As converge para um limite, digamos A∞. Agora, Rs = L−1 s As → IT−1 ∞ A1T∞ e, assim, Rs converge para o mesmo limite que As. Como Rs é uma matriz triangular para todo s, seu limite é necessariamente uma matriz do mesmo tipo. Esta matriz deve ter os autovalores de A1 em alguma ordem na diagonal principal, pois, de (3.13), ela é similar a A1. Note que a demonstração não é afetada pela presença de autovalores de mesmo módulo e nem pela existência de algum menor principal de X que seja igual a zero. Porém, os autovalores podem não estar em ordem decrescente de multiplicidade na diagonal de A. Se calcularmos com precisão exata, obteremos uma matriz limite com os autovalores desordenados na diagonal principal. Novamente, há o risco de erros de arredondamento por parte do computador, o que pode causar convergência lenta para uma matriz triangular com os autovalores ordenados. 3.1.3 Autovalores complexos Quando A1 é uma matriz real e possui um ou mais pares de autovalores complexos, é óbvio que As não pode convergir para uma matriz triangular superior tendo esses autovalores na diagonal principal, pois todas as matrizes As são reais. Primeiramente, vamos considerar o caso em que a matriz possui apenas um par de autovalores complexos, digamos λm−1 e λm−1, e os outros autovalores satisfazem |λi| > |λj | se i < j. Vamos supor, também, que os menores principais de 36 Algoritmo LR de Rutishauser X e Y sejam todos diferentes de zero. De (3.10), temos que t (s) ji → det(Xji) det(Xii) , i 6= m− 1, exatamente como no caso anterior. Para i = m − 1, no entanto, existem termos da forma (λ1λ2...λm−2 λm−1)s e (λ1λ2...λm−2λm−1)s no numerador e no denominador. Como desejamos estudar o comportamento assintótico desses elementos, é conveniente introduzir uma notação mais simplificada para os menores de ordem (m− 1) de X e Y . Assim, sejam zj = ∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣ x11 x12 · · · x1,m−1 x21 x22 · · · x2,m−1 ... ... ... ... xm−2,1 xm−2,2 · · · xm−2,m−1 xj1 xj2 · · · xj,m−1 ∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣ e wm−1 = ∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣ y11 y12 · · · y1,m−1 y21 y22 · · · y2,m−1 ... ... ... ... ym−2,1 ym−2,2 · · · ym−2,m−1 ym−1,1 ym−1,2 · · · ym−1,m−1 ∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣ . Lembrando que as colunas (m − 1) e m de X são complexas conjugadas, assim como as linhas (m− 1) e m de Y , então zj = ∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣ x11 x12 · · · x1,m−2 x1,m x21 x22 · · · x2,m−2 x2,m ... ... ... ... ... xm−2,1 xm−2,2 · · · xm−2,m−2 xm−2,m xj1 xj2 · · · xj,m−2 xj,m ∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣ e wm−1 = ∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣ y11 y12 · · · y1,m−1 y21 y22 · · · y2,m−1 ... ... ... ... ym−2,1 ym−2,2 · · · ym−2,m−1 ym1 ym2 · · · ym,m−1 ∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣ . Da equação (3.10), obtemos t (s) j,m−1 ∼ zjwm−1(λ1λ2...λm−1)s + zjwm−1(λ1λ2 . . . λm−1)s zm−1wm−1(λ1λ2...λm−1)s + zm−1wm−1(λ1λ2...λm−1)s , (3.14) que não converge quando s →∞. Aqui, o símbolo ∼ significa lim s→∞ t (s) j,m−1 zjwm−1(λ1λ2...λm−1)s+zjwm−1(λ1λ2...λm−1)s zm−1wm−1(λ1λ2...λm−1)s+zm−1wm−1(λ1λ2...λm−1)s = 1. Podemos expressar (3.14) como t (s) j,m−1 ∼ aszj + aszj aszm−1 + aszm−1 e, assim, temos que t (s+1) j,m−1 − t (s) j,m−1 ∼ (zm−1zj − zm−1zj)(as+1as − as+1as) (as+1zm−1 + as+1zm−1)(aszm−1 + aszm−1) = Ps(zm−1zj − zm−1zj). 37 Algoritmo LR de Rutishauser A diferença entre as colunas (m − 1) de Ts+1 e Ts é, portanto, um múltiplo do vetor com componentes zm−1zj−zm−1zj . Aplicando o Teorema de Sylvester (Teorema 2.14) nas expressões de zj , zj , zm−1 e zm−1, obtemos zm−1zj − zm−1zj = (x1,2,...,m−2)(x (j) 1,2,...,m). O primeiro fator é independente de j e, então, lim s→∞ t (s) jm = x (j) 1,2,...,m x (m) 1,2,...,m . Logo, podemos escrever t (s+1) j,m−1 − t (s) j,m−1 ∼ Qst (s) jm. Demonstramos, assim, que todas as colunas de Ts convergem para um limite, exceto a coluna (m − 1). Por último, demonstramos que as colunas (m − 1) das matrizes Ts+1 e Ts diferem por múltiplos da coluna m de Ts. Como Ts+1 = TsLs+1, concluímos que Ls+1 ∼   1 . . . 1 cs 1 . . . 1   , onde o elemento cs está na posição (m,m− 1). Utilizando a relação T−1 s+1 = L−1 s+1T −1 s , temos que todas as linhas de T−1 s tendem a um limite, exceto a linha m. Vimos que a linha m de T−1 s+1 difere da linha m de T−1 s por um múltiplo, cs, da linha (m− 1). Assim, As+1 = T−1 s A1Ts tende para um limite, excetuando-se a linha m e a coluna (m− 1). Agora, Ls+1 é obtida da decomposição LR de As+1 e, do limite de Ls+1, temos que todos os elementos abaixo da diagonal principal de As+1 tendem a zero, exceto o da posição (m,m − 1), ou algum dos elementos da diagonal de As+1 tende a infinito. Agora, de (3.14), t (s) j,m−1 pode ser escrito da forma t (s) j,m−1 ∼ M cos(sθ + αj + βm−1) cos(sθ + αm−1 + βm−1) . A não ser que tenhamos θ = dπ e , onde d e e são inteiros, a expressão à direita assumirá um valor muito grande. No entanto, será da ordem de M para infinitos valores de s. Assim, podemos excluir a possibilidade dos elementos de As+1 tenderem a infinito. 38 Algoritmo LR de Rutishauser Podemos estender nossa análise para o seguinte resultado. Seja A uma matriz real, cujos autovalores estão ordenados da seguinte forma: |λ1| ≥ |λ2| ≥ . . . ≥ |λn|. Denotemos um autovalor real por λr e um par de autovalores complexos por λc−1 e λc, ou seja λc−1 = λc. Se: 1. λi = λi+1, exceto para os pares conjugados; 2. det(Xii) det(Yii) 6= 0, i = 1, . . . , m; 3. em cada estágio a matriz Ar pode ser decomposta no produto LU, então, • a (s) ij → 0 se i > j, exceto para a (s) c,c−1; • a (s) rr → λr; • todos os elementos a (s) ij , j > i, tendem para um limite, exceto aqueles que estão na linha c e coluna c− 1; • a submatriz   a (s) c−1,c−1 a (s) c−1,c a (s) c,c−1 a (s) c,c   não converge, mas seus autovalores convergem para os correspondentes λc−1 e λc; • os elementos t (s) ij tendem ao limite de (3.11), exceto aqueles da coluna (c− 1). Podemos, agora, utilizar o algoritmo LR como uma técnica prática. À primeira vista pode não parecer tão útil, pois (a) Existem matrizes que não possuem decomposição LR, no sentido de que seus problemas de autovalores não são bem condicionados. Sem uma modificação no algoritmo, não podemos calcular seus autovalores através do algoritmo LR. Além disso, existe uma grande classe de matrizes para as quais a decomposição LR é numericamente instável. Podemos ter instabilidade numérica em qualquer passo do processo iterativo, o que nos levaria a uma grande perda de precisão. (b) O tempo computacional é muito alto. Existem 2 3 n3 multiplicações em cada passo de cada iteração, sendo que metade é no processo de triangularização e a outra metade na pós- multiplicação. 39 Algoritmo LR de Rutishauser (c) Os elementos abaixo da diagonal principal convergem para zero a uma taxa de λr+1 λr e será muito devagar se os autovalores forem próximos entre si. Uma forma de contornar esses problemas no algoritmo LR é através da modificação que definiremos a seguir. 3.1.4 Algoritmo LR modificado Seja A uma matriz qualquer. Então, existem matrizes Ir,r′ e Nr tais que Nn−1In−1,(n−1)′ . . . N2I2,2′N1I1,1′A = R, (3.15) onde todos os elementos de Nr, em valor absoluto, são limitados por 1. As matrizes Ir,r′ são a matriz identidade com as linhas r e r′ trocadas, sendo que r′ é definido por |ar−1 r,r′ | = max |ar−1 i,r |, i = r, ..., n. Se r′ não for único, toma-se o de menor valor. Podemos completar a transformação de similaridade de A se pós-multiplicarmos R por I1,1′N −1 1 I2,2′N −1 2 . . . In−1,(n−1)′N −1 n−1 e, assim, obter Nn−1In−1,(n−1)′ . . . N2I2,2′N1I1,1′AI1,1′N −1 1 I2,2′N −1 2 . . . In−1,(n−1)′N −1 n−1 = = RI1,1′N −1 1 I2,2′N −1 2 . . . In−1,(n−1)′N −1 n−1. (3.16) No caso particular em que não é necessário a troca de linhas, ou seja, Ir,r′ = I, r = 1, . . . , n− 1, a matriz produto que precede A em (3.15) é a matriz L−1, para a qual L−1A = R. Neste caso, a matriz à direita, na equação (3.16), é a matriz RL. Assim, definimos a modificação no algoritmo LR. Em cada passo s, reduzimos a matriz As a uma matriz triangular superior Rs usando eliminação de Gauss com troca de linhas. Em seguida, pós-multiplicamos Rs pelas inversas dos fatores usados na redução. Obtemos, então, a matriz As+1. O esforço computacional é o mesmo que o do algoritmo LR original. A fim de evitar uma análise complexa, demonstraremos a relação entre A1 e A2 no caso em que n = 4. A análise pode ser estendida para o caso geral sem grandes dificuldades. Temos, então, N3I3,3′N2I2,2′N1I1,1′A1 = R1 40 Algoritmo LR de Rutishauser e A2 = {N3I3,3′N2I2,2′N1I1,1′}A1{I1,1′N −1 1 I2,2′N −1 2 I3,3′N −1 3 } = {N3I3,3 ′N2(I3,3 ′ I3,3 ′ )I2,2 ′N1(I2,2 ′ I3,3 ′ I3,3 ′ I2,2 ′ )I1,1 ′A1} (3.17) ×{I1,1′ (I2,2′ I3,3′ I3,3′ I2,2′ )N −1 1 I2,2′ (I3,3′ I3,3′ )N −1 2 I3,3′N −1 3 }, onde os termos entre parênteses são iguais a I. Lembramos que I3,3′ I2,2′N1I2,2′ I3,3′ = Ñ1 e I3,3′N2I3,3′ = Ñ2, onde Ñ1 e Ñ2 são da mesma forma que N1 e N2, mas com os elementos da sub-diagonal reorde- nados. Reagrupando os termos em (3.17), obtemos A2 = L̃−1 1 A1L̃1, L̃1 = Ñ−1 1 Ñ−1 2 Ñ−1 3 e Ã1 = I3,3′ I2,2′ I1,1′A1I1,1′ I2,2′ I3,3′ . Podemos ver que L̃1 é uma matriz triangular inferior com 1 na diagonal principal e Ã1 é A1 com suas linhas e colunas permutadas. Se escrevermos I3,3′ I2,2′ I1,1′ = P1, onde P1 é uma matriz de permutação, então Ã1 = P1A1P t 1, P1A1 = L̃1R1 e A2 = R1P t 1L̃1 = L̃−1 1 (P1A1P t 1)L̃1. A prova da convergência do algoritmo LR não se aplica ao algoritmo modificado. No entanto, se houver convergência para uma matriz triangular superior e se nenhum dos autovalores é zero, é óbvio que, a partir de um certo ponto, as trocas de linhas não são mais necessárias, pois os elementos abaixo da diagonal convergem para zero. Infelizmente, existem matrizes para as quais o algoritmo LR converge e o algoritmo modi- ficado não converge, como é o caso da matriz A1 =   1 3 2 0   , cujos autovalores são λ1 = 3, e λ2 = −2. Temos A2 =   1 2 3 0   e A3 =   1 3 2 0   = A1, com uma troca de linhas nas reduções de A1 e A2. Como A3 = A1, o processo é cíclico e, assim, não converge. Por outro lado, no algoritmo original, a matriz converge com os autovalores na ordem correta. 41 Algoritmo LR de Rutishauser 3.1.5 Matrizes na forma de Hessenberg Com as observações feitas, pode parecer que o algoritmo LR não é uma boa escolha se A1 é uma matriz com poucos zeros. Se a matriz A1 possui a característica de que os elementos nulos são mantidos ao longo do algoritmo, então o número de operações pode ser consideravelmente reduzido. Existem dois tipos de matrizes com essas características: as de Hessenberg superior, que iremos nos referir apenas como matrizes de Hessenberg, e as matrizes em banda. Primeiramente, demonstraremos que uma matriz de Hessenberg é invariante com respeito ao algoritmo modificado. No r-ésimo estágio da redução, a escolha do pivot está entre as linhas r e (r+1) e a matriz Nr é igual à matriz identidade exceto na posição (r, r+1), onde o elemento é diferente de zero. Vamos mostrar que a matriz A2 = R1I ′ 12N −1 21 I ′ 23N −1 32 . . . I ′ n−1,2N −1 n,n−1 é uma matriz de Hessenberg. Aqui, I ′ r,r+1 denota uma matriz na forma Ir,r+1 ou a matriz identidade, de acordo com a necessidade de troca de linhas ou não no r-ésimo estágio da redução triangular. A prova será feita por indução. Vamos supor que a matriz R1I ′ 12N −1 21 I ′ 23N −1 32 . . . I ′ r−1,rN −1 r,r−1 é uma matriz de Hessenberg nas primeiras (r − 1) colunas e triangular superior no restante da matriz. A seguir, mostramos um exemplo com n = 6 e r = 4:   × × × × × × × × × × × × × × × × × × × × × × × ×   . O próximo passo é pós-multiplicá-la por I ′ r,r+1. Se esta matriz for igual a Ir,r+1, devemos trocar as colunas r e (r + 1), caso contrário a matriz não se altera. Logo, a matriz resultante é de uma 42 Algoritmo LR de Rutishauser das formas (a) ou (b) abaixo: (a)   × × × × × × × × × × × × × × × × × × × × × × 0 × ×   , (b)   × × × × × × × × × × × × × × × × × × × × × × × ×   . A pós-multiplicação por N−1 r,r+1, resulta na adição de um múltiplo da coluna (r+1) à coluna r e, assim, a matriz resultante é de uma das formas (c) ou (d): (c)   × × × × × × × × × × × × × × × × × × × × × × 0 × ×   , (d)   × × × × × × × × × × × × × × × × × × × × × × × × ×   . Em qualquer um desses casos, a matriz é da forma de Hessenberg nas primeiras r colunas e triangular nas demais, o elemento igual a zero em (c) não é importante. Desta forma, a indução está provada. Para reduzir a matriz à forma triangular, são necessárias 1 2 n2 multiplicações e outras 1 2 n2 multiplicações na pós-multiplicação. Temos, assim, n2 multiplicações em cada iteração do algoritmo LR para matrizes de Hessenberg, enquanto que, para uma matriz qualquer, teríamos 2 3 n3. 3.1.6 Aceleração da convergência Vimos que a redução da matriz A à forma de Hessenberg, antes da aplicação do algoritmo, reduz consideravelmente o volume computacional. Porém, o método ainda possui taxa de convergência muita alta Vimos, também, que, em uma matriz qualquer, os elementos aij abaixo da diagonal prin- cipal convergem para zero com ( λi λj )s (raramente a convergência é rápida). Para as matrizes de Hessenberg, os únicos elementos diferentes de zero abaixo da diagonal são aqueles nas posições (r, r − 1), r = 2, . . . , n. Considere, agora, a matriz (A− σI), que possui autovalores (λi − σ) e, 43 Algoritmo LR de Rutishauser pelo algoritmo LR, a (s) n,n−1 tende para zero com ( λn−σ λn−1−ρ )s . Se σ é uma boa aproximação para λn, o elemento a (s) n,n−1 converge para zero rapidamente. Então, seria melhor trabalhar com a matriz (As − σI) do que com a matriz A. Observe que se σ é exatamente igual a λn, então, se o processo modificado fosse calculado com precisão, o elemento (n, n − 1) se tornaria zero em apenas uma iteração. Considerando o processo de triangularização, nenhum dos pivots pode ser igual a zero, exceto o último, pois, em cada estágio da redução, o pivot pode ser ar,r+1 ou um outro elemento e, neste caso, estamos supondo que todos os elementos ar,r+1 são diferentes de zero. Assim, como det(A− λnI) = 0, o último pivot tem que ser igual a zero e a matriz R é da forma   × × × × × × × × × × × × × × 0   , onde a última linha é inteiramente nula. Como a pós-multiplicação apenas combina duas colunas consecutivas, ao fim da iteração devemos obter uma matriz da forma   × × × × × × × × × × × × × × × × × 0 0   . É com o intuito de melhorar a convergência do algoritmo LR que apresentamos a mo- dificação a seguir. No s-ésimo estágio, ao invés de decompormos a matriz A no produto LR, decompomos a matriz (As − ksI), onde ks é um valor escolhido a priori. Desta forma, produzimos uma seqüência de matrizes da forma As − ksI = LsRs e RsLs + ksI = As+1. Logo, As+1 = RsLs + ksI = L−1 s (As − ksI)Ls + ksI = L−1 s AsLs (3.18) e, assim, as matrizes As são todas similares a A1. De fato, de (3.18), temos que As+1 = L−1 s AsLs = L−1 s L−1 s−1As−1Ls−1Ls = L−1 s ...L−1 2 L−1 1 A1L1L2...Ls, 44 Algoritmo LR de Rutishauser ou L1L2...LsAs+1 = A1L1L2...Ls. Esta forma modificada do algoritmo LR é conhecida como LR com shifts na origem. Em pro- gramas computacionais, utilizamos a seguinte alteração deste algoritmo: As − zsI = LsRs e RsLs = As+1. Assim, As+1 = RsLs = L−1 s (As − zsI)Ls = L−1 s AsLs − zsI = L−1 s (L−1 s−1As−1Ls−1 − zs−1I)Ls − zsI = L−1 s L−1 s−1As−1Ls−1Ls − (zs−1 + zs)I. Continuando com este raciocínio, obtemos As+1 = L−1 s ...L−1 2 L−1 1 A1L1L2...Ls − (z1 + z2 + ... + zs)I = L−1 s ...L−1 2 L−1 1 [A1 − (z1 + z2 + ... + zs)]L1L2...Ls, o que implica que os autovalores de As+1 diferem dos de A1 por ∑s i=1 zi. Voltando ao algoritmo com shifts na origem, temos L1L2...Ls−1(LsRs)Rs−1...R2R1 = L1L2...Ls−1(As − ksI)Rs−1...R2R1 = (As − ksI)L1L2...Ls−1Rs−1...R2R1 = (A1 − ksI)(A1 − ks−1I)L1L2...Ls−2Rs−2...R2R1 = (A1 − ksI)(A1 − ks−1I)...(A1 − k1I). (3.19) Logo, fazendo Ts = L1L2...Ls e Us = Rs...R2R1, TsUs é a decomposição LR de s∏ i=1 (A1 − kiI), sendo que a ordem dos fatores não importa. Claramente podemos combinar o uso de shifts na origem com o uso de trocas de linhas e, assim, cada elemento da seqüência de matrizes As é similar a A1. Porém, não podemos escrever como em (3.19). Temos, agora, um outro problema, o de encontrar a seqüência dos ks que produzirá con- vergência mais rápida. Se todos os autovalores possuem módulos distintos, então o elemento a (s) n,n−1 converge para zero enquanto que a (s) nn converge para λn, podendo, os autovalores, ser 45 Algoritmo LR de Rutishauser complexos ou reais. Assim, tomamos ks = a (s) nn se a (s) n,n−1 converge para zero ou a (s) nn demonstra sinais de convergência. De fato, podemos ver que, quando a (s) n,n−1 é da ordem de ε, então, se escolhermos ks = a (s) nn, temos que a (s+1) n,n−1 = O(ε2). Neste caso, com n = 6, temos As − a(s) nnI =   × × × × × × × × × × × × × × × × × × × × × × × × ε 0   . Reduzindo esta matriz à forma LR utilizando eliminação de Gauss com trocas de linhas, obtemos, como último estágio, a seguinte matriz:   × × × × × × × × × × × × × × × × × × a b ε 0   . O elemento a na linha (n−1) não está próximo de zero, a não ser que a (s) nn possua alguma relação especial com o menor principal de ordem (n − 1) de As. Por exemplo, se a (s) nn é um autovalor deste menor principal, então não é necessária a troca de linha neste último estágio e, assim, mn,n−1 = ε a , o que nos fornecerá a matriz   × × × × × × × × × × × × × × × × × × a b −bε a   . Ao pós-multiplicarmos esta matriz por todos os fatores, exceto I ′n−1,nN−1 n,n−1, a matriz 46 Algoritmo LR de Rutishauser obtida tomará uma das duas formas abaixo: (a)   × × × × × × × × × × × × × × × × × × × × × a 0 b −bε a   (b)   × × × × × × × × × × × × × × × × × × × × × × a b −bε a   . Se I ′n−2,n−1 = In−2,n−1, teremos a matriz em (a). Agora, se I ′n−2,n−1 = I, teremos a matriz em (b). Ao pós-multiplicarmos esta matriz pelo fator I ′n−1,nN−1 n,n−1, obtemos (c)   × × × × × × × × × × × × × × × × × × × × × × bε a b bε2 a2 −bε a   (d)   × × × × × × × × × × × × × × × × × × × × × × −a + bε a b −bε2 a2 −bε a   Se I ′n−1,n = In−1,n, teremos a matriz em (c). Caso contrário, teremos a matriz em (d). Ao completar o estágio do algoritmo LR com shift na origem, obtemos a(s+1) nn = a(s) nn − bε a , a (s+1) n,n−1 − bε2 a2 e, assim, temos que a (s) nn está convergindo e a (s+1) n,n−1 é da ordem de ε2, ou seja, está convergindo para zero. Observe que as trocas de linhas não surtem efeito na convergência. 3.1.7 Deflação de uma matriz Em geral, quando o elemento a (s) n,n−1 se torna pequeno, ele converge para zero rapidamente. Mas, como nem sempre dispomos de precisão exata, dificilmente teremos a (s) n,n−1 = 0 . Assim, ao atingir um valor suficientemente pequeno, o trataremos como zero e diremos que a (s) nn é um autovalor. Os autovalores restantes são os mesmos da matriz formada pelas primeiras (n − 1) linhas e colunas. Esta matriz também está na forma de Hessenberg e, assim, podemos continuar com o mesmo método mas, agora, com uma matriz de ordem (n−1). Como esperamos que todos os elementos abaixo da diagonal principal sejam iguais a zero, o elemento a (s) n−1,n−2 pode estar bem próximo de zero e, neste caso, podemos usá-lo como nosso próximo ks. 47 Algoritmo LR de Rutishauser Continuando com este processo, podemos obter, um a um, todos os autovalores trabalhando sempre com matrizes cada vez menores. Este método é conhecido como deflação. O último estágio de convergência de cada autovalor será, geralmente, quadrático e, além disso, imaginamos que, ao encontrar um autovalor, já tenhamos um bom ponto de partida para obtermos o autovalor seguinte. Este método de deflação é muito estável. De fato, os autovalores serão afetados apenas se forem “sensíveis” ao elementos logo abaixo da diagonal principal. Como supomos que, em geral, a matriz de Hessenberg original A1 foi obtida da redução de uma matriz A, se os autovalores forem “sensíveis” a estes elementos, então já devemos ter perdido a precisão antes mesmo de começarmos o algoritmo LR. O único perigo adicional é de que a condição sobre As possa se perder progressivamente. 3.1.8 Shift na origem A motivação para utilizarmos shift na origem é quando temos uma matriz real com autovalores complexos. Neste caso, teremos, como matriz limite, uma matriz bloco triangular, como vimos no Teorema de Schur. Se o autovalor de menor módulo é o de um par complexo, digamos λ± iµ, então o algoritmo com shift na origem deverá tender à uma matriz da forma   × × × × × × × × × × × × × × ε × × × ×   , onde os autovalores da submatriz 2×2 no canto inferior direito convergem para λ±iµ. Por outro lado, se o autovalor de menor módulo possui módulo igual a λ, deveremos obter uma matriz na forma   × × × × × × × × × × × × × × × × × ε ×   . Vamos considerar os autovalores da matriz 2×2 nos dois casos. No primeiro caso, os autovalores vão tender para λ ± iµ como já vimos e, no segundo, o menor dos dois autovalores vai tender 48 Algoritmo LR de Rutishauser para λ. Vamos analisar este caso. Para matrizes reais com autovalores também reais, utilizaremos a seguinte escolha para os ks. Calculamos os autovalores desta matriz 2 × 2 e, se os autovalores são complexos e iguais a γs ± iδs, tomamos ks = γs. Agora, se os autovalores são reais e iguais a γs e δs, então tomamos ks = fs, onde fs = min(|γs − a (s) nn|, |δs − a (s) nn|). Se a matriz original é complexa, então, em geral, teremos dois autovalores complexos αs e βs e, novamente, tomaremos ks como fs. Voltando ao caso em que o menor autovalor em módulo é o de um par complexo, esperamos que os autovalores da matriz 2×2 sejam aproximações para este par. Suponha que obtemos como autovalores da matriz 2 × 2, o par γ ± iδ. Assim, se nos dois próximos estágios do algoritmo LR utilizarmos, como shifts, γ − iδ e, depois, γ + iδ, certamente estaremos trabalhando com aritmética complexa, mas, se dispomos de precisão exata, a matriz que obteremos em seguida será real. Vamos denotar essas três matrizes por A1, A2 e A3. Logo, A1 − (γ − iδ)I = L1R1, R1L1 + (γ − iδ)I = A2, A2 − (γ + iδ)I = L2R2, R2L2 + (γ + iδ)I = A3. Utilizando o resultado (3.19), temos L1L2R2R1 = [A1 − (γ − iδ)I][A1 − (γ + iδ)I] (3.20) e, obviamente, a matriz à direita é real. Temos, também, que esta matriz é não singular, pois γ ± iδ não são os autovalores exatos de A1 e, assim, sua decomposição triangular, se existir, é única e as matrizes L e R são reais. Da equação (3.20), concluímos que L3 = L1L2 e R3 = R2R1. Logo, L3 é uma matriz real. Utilizando, ainda, o fato de que A3 = L−1 3 A1L3 = (L1L2)−1A1L1L2, temos que A3 é uma matriz real. Ao executarmos esses dois estágios do algoritmo LR no computador, deveremos obter A3 com valores imaginários. O processo é numericamente estável, no sentido de que os autovalores de A3 estão bem próximos dos autovalores de A1, mas não temos garantia de que a matriz final obtida esteja próxima da matriz esperada com cálculos exatos. 49 Capítulo 4 Algoritmo QR de Francis Fazemos, neste capítulo, um estudo sobre o algoritmo QR. Desenvolvido por Francis em [9], o algoritmo QR utiliza transformações unitárias para reduzir uma matriz à forma triangular superior. Assim como no capítulo anterior, o livro de Wilkinson [30] foi, também, uma referência para o estudo apresentado neste capítulo. Além disso, apresentamos o Método Assemelhado ao de Francis (MAF) e fazemos, também, um estudo sobre a eficiência deste novo método comparado com os algoritmos QR de Francis e LR de Rutishauser. O MAF distingue-se do método de Francis pela maneira mais simples e rápida de se obter as matrizes Qk e, do de Rutishauser, por se aplicar a classe mais geral de matrizes. 4.1 Algoritmo QR de Francis Uma outra maneira de se calcular os autovalores de uma matriz A é o algoritmo QR, descrito por Francis [9]. Ao invés de utilizar decomposição triangular, Francis utilizou a decomposição QR. O algoritmo é definido da seguinte forma: A1 = A, As = QsRs e As+1 = RsQs, s = 1, 2, . . . . Em cada passo do algoritmo, usamos uma transformação de similaridade com matrizes unitárias. Se As é real, então Qs e Rs são reais. Uma vantagem deste algoritmo, em relação ao algoritmo LR original, é que se As possui um menor principal igual a zero, o algoritmo continua funcionando. As iterações sucessivas satisfazem relações similares àquelas do algoritmo LR. Vejamos: As+1 = QH s AsQs = QH s QH s−1As−1Qs−1Qs = QH s ...QH 2 QH 1 A1Q1Q2...Qs 50 Algoritmo QR de Francis e, assim, Q1Q2...QsAs+1 = A1Q1Q2...Qs. Portanto, todas as matrizes As são similares a A1 e, se escrevermos Ps = Q1Q2...Qs e Us = RsRs−1...R1, então PsUs = Q1Q2...(QsRs)Rs−1...R1 = Q1...Qs−1AsRs−1...R1 = A1Q1...Qs−1Rs−1...R1 = A1Ps−1Us−1. Continuando este raciocínio, obtemos PsUs = As 1 e, então, demonstramos que Ps e Us são as matrizes da correspondente decomposição QR de As 1. Como exigimos que os elementos da diagonal principal de R sejam positivos, esta decomposição é única. 4.1.1 Convergência do algoritmo QR Em geral, as matrizes As convergem para uma matriz triangular superior sob menos condições do que as necessárias no algoritmo LR. Mas, antes, vejamos uma diferença entre os dois algoritmos. Se aplicarmos o algoritmo LR a uma matriz A1 triangular superior, temos L1 = I, R1 = A1 e, assim, As = A1 para todo s. Agora, se aplicarmos o algoritmo QR a esta mesma matriz, obteremos A1 = D(D−1A1) e A2 = D−1A1D, onde ajj = |ajj |eiθj e D = diag[eiθi ]. Então, embora A2 tenha a mesma diagonal principal que A1, os elementos acima dela estão multiplicados por fatores complexos de módulo 1. Obviamente, por esta razão, não devemos esperar que As convirja para um limite, a não ser que todos os λj sejam reais e positivos. Os fatores complexos, de módulo 1, são de pouca importância e, assim, podemos dizer que As converge essencialmente se existe alguma matriz diagonal unitária D tal que As+1 = D−1AsD,. 4.1.2 Prova formal da convergência Seja A uma matriz quadrada de dimensão n cujos autovalores satisfazem |λ1| > |λ2| > ... > |λn|. 51 Algoritmo QR de Francis Pelo algoritmo QR, temos A1 = A. Como A1 possui divisores lineares (autovalores diferentes), podemos escrever As 1 = XDsX−1 = X   λs 1 . . . λs n  Y. Sejam as matrizes Q, R, L e U de forma que X = QR e Y = LU, onde R e U são triangulares superiores, L triangular inferior com 1 na diagonal principal e Q é unitária. Note que: 1. as quatro matrizes não dependem de s; 2. a não singularidade de R segue de X; 3. a decomposição QR sempre existe, mas a decomposição LU de Y existe apenas se todos os menores principais formados pelas k primeiras linhas e k primeiras colunas, k = 1, ...n− 1, são diferentes de zero. Temos, então, As 1 = QRDsLU = QR(DsLD−s)DsU. (4.1) Observe que DsLD−s é uma matriz triangular inferior com 1 na diagonal principal: DsLD−s =   λs 1 λs 2 . . . λs n     1 l21 1 ... ... . . . ln1 ln2 · · · 1     1 λs 1 1 λs 2 . . . 1 λs n   =   λs 1 0 · · · 0 λs 2l21 λs 2 · · · 0 ... ... ... ... λs nln1 λs nln2 · · · λs n     1 λs 1 1 λs 2 . . . 1 λs n   =   1 0 · · · 0( λ2 λ1 )s l21 1 · · · 0 ... ... . . . ...( λn λ1 )s ln1 ( λn λ2 )s ln2 · · · 1   . 52 Algoritmo QR de Francis Assim, podemos escrever DsLD−s = I + Es, onde Es → 0 quando s →∞. Da equação (4.1), temos As 1 = QR(I + Es)DsU = Q(I + REsR −1)RDsU = Q(I + Fs)RDsU, (4.2) onde Fs é uma matriz estritamente triangular e, além disso, Fs = REsR −1. Mostremos que Fs → 0 quando s →∞. Como o elemento (i, j) de Fs satisfaz (Fs)i,j = n∑ k=1 (REs)ik(R−1)kj = n∑ k=1 n∑ l=1 Ril(Es)lk(R−1)kj , fazendo s →∞, temos que (Fs)ij → 0, pois (Es)lk → 0 para l, k = 1, ..., n. Agora, podemos escrever I + Fs = Q̃sR̃s. Como Fs → 0 quando s →∞, temos que Q̃sR̃s → I. Além disso, det(I + Fs) = 1 e, então, R̃s é inversível e det(R̃s) = 1 para todo s. De (4.2), obtemos As 1 = (QQ̃s)(R̃sRDsU). (4.3) A matriz do primeiro parênteses é unitário, enquanto que, a do segundo, é triangular superior. Como As 1 é não-singular, sua decomposição é única e, desta forma, temos que Ps e Q são possivelmente diferentes por uma pós-multiplicação por uma matriz diagonal unitária. Assim, podemos dizer que Ps → Q. Se insistirmos, ainda, que todas as matrizes Rs devem ter elementos positivos na diagonal principal, podemos encontrar a decomposição unitária diagonal de (4.3). Sejam D = |D|D1 e U = D2(D−1 2 U), onde D1 e D2 são matrizes unitárias diagonais, D−1 2 U possui os elementos da diagonal principal positivos (R̃s e R já possuem elementos positivos na diagonal principal) e |D| é a matriz D com os elementos tomados em módulo. Assim, de (4.3), temos As 1 = QQ̃sD2D s 1[(D2D s 1) −1R̃sR(D2D s 1)|D|s(D−1 2 U)]. A matriz entre colchetes é triangular superior com os elementos da diagonal principal positivos e, como Ps ∼ QD2D s 1, temos que Qs → D1, pois Ps = Q1Q2...Qs−1Qs = Ps−1Qs e, portanto, Qs = P ∗ s−1Ps. Logo, Qs ∼ (Ds−1 1 )∗D∗ 2Q ∗QD2D s 1 = (Ds−1 1 )∗Ds−1 1 D1 = D1. 53 Algoritmo QR de Francis 4.1.3 Autovalores desordenados Demonstramos até agora que, se os menores principais de Y são diferentes de zero, além da convergência, temos, também, que os λi da diagonal principal estão ordenados como gostaríamos. É óbvio que se um menor principal de X for igual a zero, ainda assim Ps converge para um limite. Apresentaremos, agora, uma demonstração para o caso em que existem menores principais de Y , formados pelas primeiras i linhas e i colunas, iguais a zero. Embora Y não possua decomposição LR quando um de seus menores principais é igual a zero, sempre existe uma matriz de permutação P tal que PY possui decomposição LR. Uma forma de se determinar P é utilizando a técnica de pivoteamento. Na eliminação de Gauss, em cada estágio, tomamos como pivot o maior elemento em módulo da coluna em questão. Ao invés disso, no estágio r, escolhemos como pivot o primeiro elemento diferente de zero dentre os elementos das posições (r, r), (r + 1, r),. . . , (n, r). Denotemos esta posição por (r′, r). Em seguida, rearranjamos as linhas na ordem r′, r, r + 1, ..., r′− 1, r′ + 1, ..., n. Finalmente, temos PY = LR, onde L possui zeros abaixo da diagonal principal nas posições relacionadas à matriz de permu- tação. Como Y é não-singular, sempre vai existir um pivot em cada estágio. Assim, podemos escrever As 1 = XDsY = XDsP tLU = XP t(PDsP t)LU. A matriz PDsP t é uma matriz diagonal com os elementos λs i em uma ordem diferente, enquanto XP t é obtida de X através da permutação de suas colunas. Se escrevermos, ainda, XP t = QR e PDsP t = Ds 3, então As 1 = QRDs 3LU = QR(Ds 3LD−s 3 )Ds 3U. Como os elementos λs i não estão na ordem natural em D3, é natural pensarmos que a matriz Ds 3LD−s 3 não convergirá para a matriz identidade. Porém, se denotarmos os elementos da diago- nal de D3 por λi1 , λi2 ,. . . ,λin , então o elemento (p, q), p > q, de Ds 3LD−s 3 é ( λip λiq )s lpq e o pivotea- mento em Y garante que lpq = 0 quando ip < iq. Então, podemos escrever Ds 3LD−s 3 = I + Es, onde Es → 0 quando s →∞. A matriz Ps converge essencialmente para a matriz Q, obtida da decomposição de XP t. 54 Algoritmo QR de Francis 4.1.4 Autovalores de mesmo módulo Vamos considerar, agora, o caso em que A1 possui alguns autovalores de mesmo módulo, mas com todos os divisores elementares lineares. Além disso, vamos supor que todos os menores principais de Y sejam diferentes de zero, pois já analisamos seu efeito. Temos, então, As 1 = XDsLU. Suponhamos que |λr| = |λr+1| = ... = |λt| e todos os demais autovalores são diferentes em módulo. O elemento (i, j) de DsLD−s, abaixo da diagonal principal, é ( λi λj )s lij e, então, converge para zero a não ser que t ≥ i > j ≥ r. (4.4) Neste caso, o módulo do elemento é igual a lij . Se todos os autovalores de mesmo módulo são iguais, então podemos escrever DsLD−s = L̂ + Es, onde Es → 0 quando s →∞, L̂ é uma matriz triangular inferior fixa com 1 na diagonal principal e igual à matriz I exceto nos elementos (i, j) que satisfazem à relação (4.4), onde são iguais a lij . Escrevendo XL̂ = QR, temos As 1 = QR(I + L̂−1Es)DsU = Q(I + RL̂−1EsR −1)RDsU = Q(I + F̂s)RDsU, onde F̂s = RL̂−1EsR −1 → 0 quando s →∞. Assim, As 1 = (QQ̂s)(R̂sRDsU), onde Q̂sR̂s é a decomposição QR de (I + F̂s). Então, a menos de uma pós-multiplicação por uma matriz diagonal unitária, Ps converge para Q, a matriz obtida da decomposição de XL̂. Note que as colunas de XL̂ formam um conjunto de autovetores de A1 linearmente independentes, que difere de X apenas nas colunas r, r + 1, . . . t, que são substituídas por uma combinação linear dessas colunas. Portanto, autovalores com multiplicidade algébrica maior que 1 não evitam a convergência do algoritmo. Se os autovalores de mesmo módulo não são iguais, então a matriz L̂ é da mesma forma que antes, mas, agora, os elementos diferentes de zero abaixo da diagonal principal não são fixos. Se λj = |λj |eiθj , temos l̂jk = ljke isθj−θ. 55 Algoritmo QR de Francis A matriz XL̂ é fixa, exceto nas colunas de r a t. Para cada valor de s essas colunas consistem de uma combinação linear das correspondentes colunas de X. Então, a matriz Q da decomposição QR de XL̂ possui todas as colunas fixas exceto as colunas de r a t, que consistem de uma combinação das colunas de r a t da matriz Qx, onde X = QxRx é a decomposição QR de X. Assim, exceto pelas colunas de r a t, Ps é essencialmente convergente. O caso mais importante é o de uma matriz real com alguns autovalores reais e outros complexos. Assim, como na subseção (3.1.3), podemos mostrar que As converge essencialmente para uma matriz triangular superior, exceto para um único elemento não nulo abaixo da diag- onal principal, associado a um par de autovalores complexos. Cada um desse elementos está associado a uma matriz 2×2 centrada na diagonal e com autovalores que convergem para o par de autovalores complexos ao qual esta sub-matriz está associada. Quando temos r autovalores complexos de mesmo módulo, então, no limite, As contém uma submatriz de ordem r centrada na diagonal principal, com autovalores que convergem para esses r valores. Como exemplo, vejamos a matriz abaixo: A1 =   0 0 0 1 1 0 0 0 0 1 0 0 0 0 1 0   . Aplicando o algoritmo QR, obtemos Q1 = A1, R1 = I e A2 = R1Q1 = A1, que possui, como autovalores, λj = ei jπ 2 r, r = 0, 1, 2, 3. O processo se torna invariante e, então, temos uma submatriz 4 × 4 centrada na diagonal principal. Tais blocos podem ser difíceis de serem detectados por um processo automático, principalmente se não sabemos sua dimensão. 4.1.5 Decomposição de As Seja A1 uma matriz na forma de Hessenberg com elementos da sub-diagonal diferentes de zero. Sabemos que esses elementos tendem a zero. Assim, em alguma iteração, podemos ter As com um ou mais elementos da sub-diagonal suficientemente pequenos em módulo, a ponto de podermos 56 Algoritmo QR de Francis substituí-los por zero. Dessa forma, podemos separar a matriz em duas As =   × × × × × × × × × × ε × × × × × × × ×   =   B C 0 D   , ou seja, se ε é suficientemente pequeno, então podemos encontrar os autovalores de D e B. Note que a matriz C não altera os autovalores. Assim, em cada iteração, podemos analisar os elementos da sub-diagonal e, se encontrarmos algum elemento da ordem de ε na posição (r + 1, r), podemos continuar iterando com a matriz de ordem (n − r) do canto direito inferior até encontrarmos todos os autovalores desta matriz. Em seguida, fazemos o mesmo para a matriz do canto esquerdo superior. Se r = n− 1, então o elemento na posição (n, n) é um autovalor de A1. Em seguida, podemos eliminar a última linha e a última coluna de As e trabalhar com a matriz de ordem (n− 1). Em 1961, Francis propôs uma extensão deste algoritmo. Seja A uma matriz de Hessenberg com os elementos das posições (r + 1, r) e (r + 2, r + 1) suficientemente pequenos em módulo, ε1 e ε2. Vejamos um exemplo com n = 6 e r = 2. Tomemos A =   × × × × × × × × × × × × ε1 × × × × ε2 × × × × × × × ×   =   X Y E W   , onde E é uma partição nula de A exceto pelo elemento ε1. Vamos mostrar que, se µ é um autovalor de W , então também é um autovalor da matriz A′, onde A′ difere de A apenas pelo elemento ε1ε2 (ar+1,r+1 − µ) na posição (r + 2, r). Como (W−µI) é singular, existe um vetor x 6= 0 tal que xt(W−µI) = 0. As duas primeiras componentes de x satisfazem (ar+1,r+1−µ)x1 +ε2x2 = 0 e, assim, ε1x1 + ( ε1ε2 ar+1,r+1 − µ ) x2 = 0. Desta forma, temos que as últimas (n− r) linhas de (A′ − µI) são linearmente dependentes, ou seja, µ é um autovalor de A′. Se ε1ε2 ar+1,r+1 − µ está próximo de zero, então µ é, com certeza, um autovalor de A. Logo, se (ar+1,r+1 − µ) não está próximo de zero, podemos separar W do resto da matriz. 57 Algoritmo QR de Francis Podemos ver este resultado na matriz a seguir A1 =   X Y E W   =   4 3 2 1 5 6 3 1 4 2 1 7 10−5 6 1 5 3 10−5 2 5 7 1 2 3 4 1   Essa matriz possui elementos 10−5 nas posições (3, 2) e (4, 3). Os autovalores de W são µ1 = −0.4641031621, µ2 = −0.9999985714, µ3 = 6.4641231611 e µ4 = 5.9999785724. Como apenas µ4 está próximo de a33 podemos esperar que a matriz A1 tenha os autovalores λ1, λ2 e λ3 diferindo de µ1, µ2 e µ3, respectivamente, por valores da ordem de 10−10. Já o autovalor λ4 difere de µ4 por um valor da ordem de 10−5. Os autovalores λ5 e λ6 de A1 não alteram os demais autovalores. 4.1.6 Procedimento prático Podemos tirar proveito deste fenômeno agindo da seguinte forma. Ao final de cada iteração, determinamos o shift ks a ser utilizado na próxima iteração. Em seguida, te