Decomposicoes de matrizes utilizando conceitos de AutoVetores e Auto Valores

Felipe Schimith Batista1

1Instituto de Matematica e Estatısitica - Universidade Estadual do Rio de Janeiro

felipeschimith@gmail.com

Resumo. O trabalho visa contextualizar os teoremas e metodos da AlgebraLinear aplicados em solucoes computacionais. Este trabalho apresenta oreferencial teorico utilizado como base para resolver computacionalmentedecomposicoes que fazem pate do estudo de Auto Vetores e Auto Valores. Apre-sentaremos a ”Decomposicao de Cholesky”, a ”Decomposicao de A em QR”ea ”Decomposicao de Valor Singular”. A implementacao foi feita na linguagemde programacao Java e tambem sera documentada e comentada.

1. IntroducaoNa Matematica um sistema de equacoes lineares e um conjunto finito de equacoes linearesaplicadas num mesmo conjunto, igualmente finito, de variaveis.[ANTON 2006]

No estudo da Algebra Linear nos aprofundamos na analise de sistemas deequacoes lineares. Nela se utiliza alguns conceitos e estruturas fundamentais como veto-res, matrizes, transformacoes lineares e sistemas de equacoes lineares.

Com o desenvolvimento da teoria quantica nas decadas de 1920 e 1930, o estudode matrizes tornou-se de grande importancia. No final da decada de 1940, algunsmatematicos se perguntavam como o computador digital poderia ser empregado pararesolver o problema de autovalores de matrizes. A determinacao dos autovalores de umamatriz tem merecido grande atencao, devido a sua grande aplicacao aos diversos ramosdas Ciencias Aplicadas, como, por exemplo:[Oliveira 2006]

• a teoria das vibracoes, quer sejam mecanicas ou eletricas, dos tipos macroscopicaou microscopica;• as vibracoes de pontes ou outra estrutura solida;• as vibracoes das asas de um aviao;• a oscilacao de partıculas atomicas e moleculares nas ondas mecanicas;• a teoria dos operadores lineares diferenciais e integrais, etc...

A obtencao dos autovalores λ de uma matriz quadrada de ordem n envolvecalculos complexos por serem os zeros do polinomio caracterıstico det(A −λI) = 0, quee um polinomio de grau n em λ. A busca de solucoes para este problema acarretou odesenvolvimento de uma serie de algoritmos, os mais diversos, cada vez mais poderosose mais eficientes. O primeiro a surgir foi o mais obvio, composto de apenas dois passos.Primeiro, calcula-se os coeficientes do polinomio caracterıstico e, entao, as raızes dessepolinomio.[Oliveira 2006]

Neste trabalho apesentaremos um referencial teorico e uma forma computacio-nal de implementar a ”Decomposicao de Cholesky”, a ”Decomposicao de A em QR”e a”Decomposicao de Valor Singular”.

2. Teoria de Autovetores e Autovalores

Este capıtulo tem como objetivo apresentar as propriedades, teoremas e metodos daAlgebra Linear que foram usados como base para implementar solucoes computacionais.

2.1. Autovetores e Autovalores

Antes de nos aprofundarmos nas decomposicoes de matrizes, e necessario entender comoobter os autovetores e autovalores de uma matriz.

Uma transformacao especial T : V→W.

(I) T(v) = ∆v

Onde, ∆ e o autovalor (escalar) e v e autovetor. Como toda transformacao linearpode ser escrita pela multiplicacao de uma matriz por um vetor entao:

(II) T(v) = Av

Igualando (I) e (II), tem-se:

Av = ∆v ou Av – ∆v = 0 que resulta no sistema homogeneo:

(III) (A – ∆I) v = 0

Onde A e n x n, v = 0 e sempre solucao (trivial).

Os vetores v = 0 para os quais existe um ∆ que resolve a equacao (III) sao chama-dos de autovetores da matriz A e os valores de ∆, que conjuntamente com v resolvem aequacao sao chamados de autovalores da matriz A associados aos respectivos autovetores.Para que a equacao (III) tenha solucao alem da trivial e necessario que o determinante damatriz dos coeficientes seja zero, ou seja,

det(A – ∆I) = 0

ou

det

α11 α12 α13

α21 α22 α23

α31 α32 α33

−∆

1 0 00 1 00 0 1

= 0

o que resulta em um polinomio de grau n em ∆, conhecido como polinomio ca-racterıstico. As raızes do polinomio caracterıstico sao os autovalores da matriz A. Parase encontrar os autovetores basta substituir o valor do autovalor na equacao original eencontrar o autovetor. O autovalor sera, entao, associado ao autovetor encontrado. Naverdade, o autovetor encontrado forma uma base para o espaco de solucao da equacao(III), dado o respectivo autovalor. Logo, qualquer multiplo do autovetor tambem e umautovetor.[Stanford 2002]

2.2. Decomposicao de Cholesky

Se a matriz A for simetrica e inversıvel, vamos escrever esta decomposicao na forma A =LU de modo que

LU = A = AT = UTLT .

Como L e U sao inversıveis,

U(LT )−1 = L−1 UT = D

e diagonal pois U(LT )−1 e triangular superior e L−1UT e triangular inferior. As-sim, U = DLT e obtemos a decomposicao A = LDLT onde L e triangular inferior cujoselementos diagonais sao iguais a 1 e D = L−1UT e diagonal.

Como os elementos da diagonal principal de L sao iguais a 1, D = diag(U) ondediag(U) e uma matriz diagonal, cujos elementos da diagonal principal sao iguais aos ele-mentos da diagonal de U.

Uma matriz simetrica A e positiva definida se os elementos diagonais de D nadecomposicao A = LDLT forem todos maiores do que zero. Neste caso, podemos cal-cular a matriz

√D e definir R = L

√D, para assim obter a decomposicao A = RRT

denominada de decomposicao de Cholesky da matriz A.

2.2.1. Ordem de Complexidade

Analisando a Ordem de complexidade do programa desenvolvido para obter adecomposicao de Cholesky, obtemos a Ordem (n3). Comparando com as decomposicoesde A=QR e SVD, esse metodo e mais rapido. A pesar de a complexidade ter sido demesma Ordem do A=QR e SVD, e possıvel executa-lo em uma so operacao, aumentandoseu desempenho.

2.3. Decomposicao de A em QRToda matriz A ε <mxn com colunas linearmente independentes pode ser unicamentefatorada na forma A = Q · R com as colunas da matriz Q ε <mxn formando umabase ortonormal para o range(A) e a matriz triangular superior R ε <mxn tendo os ele-mentos da diagonal positivos. Esta fatoracao pode ser gerada pelo procedimento deortogonalizacao de Gram-Schmidt desenvolvido no artigo anterior, pois as colunas damatriz Q = (q1, q2, q3...qn) sao resultados da aplicacao do processo de Gram-Schmidt nascolunas da matriz A = (a1, a2, a3...an), e a matriz R e dada por:

v1 qT1 · a2 qT1 · a3 ... qT1 · an0 v2 qT2 · a3 ... qT2 · an0 0 v3 ... qT3 · an... ... ... ... ...0 0 0 ... vn

onde v1 = ‖ a1 ‖ e vk =‖ ak −

∑k−1i=1

qTi · ak · qi ‖ para k = 2, 3, ..., n.

Se a matriz A ε <mxn e regular, entao QT = Q−1, pois as colunas da matriz Q saoortonormais entre si.

2.3.1. Ordem de Complexidade

Analisando a ordem de complexidade do programa desenvolvido para obter adecomposicao de A=QR, obtemos a Ordem (n3). A pesar dele ser mais lento que a

decomposicao de Cholesky e praticamente de mesmo desempenho que a SVD, as pro-priedades numericas dessa decomposicao torna o programa mais importante matemati-camente, pois ele pode ser aplicado em qualquer matriz mxn, podendo facilitar outroscalculos torando-o mais benefico que a decomposicao de Cholesky.

2.4. Decomposicao de Valor Singular

A Decomposicao de Valor Singular (SVD) e uma tecnica de fatoracao de matrizes queconsiste em representar qualquer matriz Amn na forma apresentada abaixo, onde U e V T

sao matrizes ortogonais, ou seja, UUT = I e V V T = I e∑

e uma matriz diagonal quecontem os valores singulares σj em ordem nao crescente σ1 σ2... ≥ σj ≥ 0 para 1 ≤ j ≤min(m,n):

A = U∑V T

Esta funcao garante a existencia e unicidade da SVD para qualquer matriz Amn,sendo que ela pode ser completa ou reduzida.

Qualquer matriz Amxn possui uma decomposicao em valores singulares (SVD) daforma A = U

∑VT . Alem disso, tem-se que os valores singulares da mesma sao sempre

unicos. Porem, os vetores singulares a direita e os vetores singulares a esquerda somenteserao unicos se a matriz A for quadrada e possuir autovalores distintos.

2.4.1. SVD Completa

A forma completa da SVD de uma matriz Amxn (m ≥ n) e composta pelas matrizes or-togonais Umxm e V T

nxn e pela matriz∑

mxn. Se a matriz Amxn tiver posto completo (r= n) existira um conjunto de n vetores ortonormais de U ε <m e um conjunto de n ve-tores ortonormais de V T ε <n. No entanto, se n < m, o conjunto de n vetores ε <m

nao sera suficiente para formar uma base de <m e Umxm nao sera uma matriz quadradaortogonal. Portanto, para obter a forma completa da SVD, sera necessario a adicao demn colunas ortonormais em U para transforma-la em uma matriz ortogonal (Umxm) e mnlinhas nulas em

∑para possibilitar a multiplicacao Umxm

∑mxn. Por outro lado, se a

matriz Amxn possuir posto incompleto (r < n), sera necessario adicionar m−r colunasortonormais em U e n−r linhas ortonormais em V T para completar as bases dos espacos<m e <n respectivamente e transformar U e V em matrizes quadradas ortogonais. Alemdisso, sera necessario adicionar m−r linhas nulas e n−r colunas nulas em

∑para pos-

sibilitar a representacao Amxn = Umxm

∑mxnV

Tnxn . Em outras palavras, qualquer matriz

pode ser diagonalizada (∑

), desde que escolhamos um sistema de coordenadas ortogonalapropriado para <m (colunas de U) e <n (linhas de V ). [Rodrigues 2011]

Figura 1. SVD completa de uma matriz Amxn(m ≥ n).

2.4.2. SVD Reduzida

A forma reduzida da SVD de uma matriz Amxn (m ≥ n) de posto r consiste na remocaodas colunas e linhas extras colocadas com intuito de transformar Umxr e Vrxn em matrizesortogonais Umxr e Vrxn. Essa remocao tem por objetivo fornecer apenas a informacaofornecida pelo conjunto de vetores definidos pelo posto da matriz.[Rodrigues 2011]

Figura 2. SVD reduzida de uma matriz Amxn(m ≥ n).

2.4.3. Ordem de Complexidade

Analisando a ordem de complexidade do programa desenvolvido para obter a SVD, ob-temos a Ordem (n3). Esse programa e mais lento que a decomposicao de Cholesky epraticamente de mesmo desempenho que a A=QR. Matematicamente ela tambem se des-taca pois e de grande aplicacao em processamento de sinais e estatıstica.

3. ProgramasEste capıtulo apresenta os programas desenvolvidos utilizando como base os conheci-mentos adquiridos na algebra linear. Algumas funcoes secundarias como obtencao dasmatrizes a partir de arquivos textos e funcao de impressao de matrizes foram omitidas dosprogramas. Focarmos a apresentacao nas partes essenciais de processamento de matrizes.

3.1. Decomposicao de CholeskyEsta primeira funcao consiste em decompor A em LU. Ele percorre todos os elementosda matriz, fazendo o escalonamento da matriz U e guardando os divisores utilizado noescalonamento em L, sendo que a diagonal e preenchida com 1 e os valores da diagonalsuperior de L e preenchido por 0.

public static void TransformaTraingularLU(Double[][]matriz, int numLin,

int numCol, Double[][] matrizL) {/** Decomposi o de A em L e U

** L

* |1 0 0|

* |R 1 0|

* |R R 1|

*

* U

* |R R R|

* |0 R R|

* |0 0 R|

** */for (int k = 0; k < numLin; k++) {//Ordem n

for (int i = k + 1; i < numLin; i++) {//Ordem n*nmatrizL[i][k] = (double) matriz[i][k].floatValue()

/ matriz[k][k].floatValue();matriz[i][k] = 0.0;for (int j = k + 1; j < numLin; j++) {//Ordem n*n*n

matriz[i][j] = (double)matriz[i][j].floatValue()

-(matrizL[i][k].floatValue()

* matriz[k][j].floatValue());

}}

}}

A segunda funcao consiste em decompor U em L’ D. Ele cria uma matriz D todade 0, e depois preenche ela com os elementos da diagonal, sendo que depois divide oselementos das colunas da diagonal de L’ pelo valor da diagonal. Apos isso ele preenche adiagonal de o L’ com 1.

public static void DecompoeUemDU(Double[][] matriz, intnumLin,

int numCol, Double[][] matrizx) {

/** U*matrizy = b (matrizx)

* Calcula U[m,n]/D[x,y] = L’ sendo que x=y

** L’

* |1 R R|

* |0 1 R|

* |0 0 1|

** D

* |R 0 0|

* |0 R 0|

* |0 0 R|

**

* */

for (int k = 0; k < numLin; k++) {//Ordem nfor (int i = 0; i < numCol; i++) {//Ordem n*n

matrizx[i][k] = 0.0;}

}for (int p = numLin - 1; p >= 0; p--) {//Ordem n

matrizx[p][p]=(double)matriz[p][p].floatValue();for (int j = p; j < numCol; j++) {//Ordem n*n

matriz[p][j]=(double)matriz[p][j].floatValue()/matrizx[p][p].floatValue();

}}}

A terceira funcao percorre todos os elementos da diagonal de D e obtem a raiz quadradadeles, formando a Raiz de D.

// Retorna Raiz de Dpublic static void DecomposicaoD(Double[][] matriz, int

numLin, int numCol,Double[][] D) {

int N = numLin;

for (int i = 0; i < N; i++) {//Ordem nfor (int j = 0; j < N; j++) {//Ordem n*n

if (i==j){//Aplica a raiz dos elementos da diagonal

principalD[i][i] =

Math.sqrt(matriz[i][i].floatValue());

}else {D[i][j] =0.0;

}}

}}

A quarta funcao faz a multipicacao matricial de L * Raiz(D), obtendo R.

public static void produtoMatriz(Double[][] matriz, intnumLin, int numCol,

Double[][] matriz2, int numLin2, int numCol2,Double[][] matrizResultado) {

int i, j, k;Double soma = 0.0;

for (i = 0; i < numLin; ++i)//Ordem nfor (j = 0; j < numCol2; ++j) {//Ordem n*n

for (k = 0; k < numCol; ++k) {//Ordem n*n*nsoma += (double)

(matriz[i][k].floatValue() *matriz2[k][j].floatValue());

}matrizResultado[i][j] = soma;soma = 0.0;

}}

3.2. Decomposicao de A=QR

Este algoritmo consiste em decompor A em Q*R, onde a k-esima coluna da matriz A, eexpressa como combinacao linear das k primeiras colunas de Q (com coeficientes R(1,k),..., R(k,k)). As k primeiras colunas de Q formam uma base ortogonal para o subespacogerado pelas k primeiras colunas de A. Se a coluna k de A e uma combinacao linear das pprimeiras colunas de A, entao, as entradas de R(p+1,k), ..., R(k,k) sao zeros. Neste caso,R e trapezoidal superior. Se A tem posto rk, as linhas R(rk+1,:), R(rk+2,:), ... sao zeros.

//Decomposicao QR

public A_QR (Matrix A) {// Inicializa.QR = A.getArrayCopy();m = A.getRowDimension();n = A.getColumnDimension();Rdiag = new double[n];

// Loop Principal.for (int k = 0; k < n; k++) {//Ordem n

// Computa a segunda normal do k esima coluna semsobrecarga.

double nrm = 0;for (int i = k; i < m; i++) {//Ordem n*n

//Utilizei o a biblioteca Jama.util.Maths paraobter a Normal.

nrm = Maths.hypot(nrm,QR[i][k]);}//Caso a normal for diferente de zeroif (nrm != 0.0) {

// Para o k esimo vetor.if (QR[k][k] < 0) {

nrm = -nrm;}for (int i = k; i < m; i++) {//Ordem n*n

QR[i][k] /= nrm;}QR[k][k] += 1.0;

// Aplica a transforma o nas colunas restantes.for (int j = k+1; j < n; j++) {//Ordem n*n

double s = 0.0;for (int i = k; i < m; i++) {//Ordem n*n*n

s += QR[i][k]*QR[i][j];}s = -s/QR[k][k];for (int i = k; i < m; i++) {//Ordem n*n*n

QR[i][j] += s*QR[i][k];}

}}//Gurada o valor do elemento da diagoalRdiag[k] = -nrm;

}}

O programa abaixo monta matriz triangular superior R de mesma dimensao queA. O primeiro For percorre as linhas ”i”e o segundo For percorre as colunas ”j”. Dentrodo Segundo For, existem tres opcoes de entrada, caso ”i¡j”entao obtem o valor calculadode QR, caso ”i=j”obtem o valor da Diagonal e caso o valor de ”i¿j”zera os elementosabaixo da diagonal.

public Matrix getR () {Matrix X = new Matrix(n,n);double[][] R = X.getArray();for (int i = 0; i < n; i++) {//Ordem n

for (int j = 0; j < n; j++) { //Ordem n*nif (i < j) {

R[i][j] = QR[i][j];} else if (i == j) {

R[i][j] = Rdiag[i];} else {

R[i][j] = 0.0;}

}}return X;

}

O programa produz computacionalmente uma matriz ortogonal Q que e obtidocom a permutacao (de colunas) E, uma matriz triangular superior R com elementos nadiagonal decrescentes e uma matriz ortogonal (ou unitaria) Q tais que A*E = Q*R. Se rke o posto de A, as rk primeiras entradas ao longo da diagonal de R. R(1,1), R(2,2), ...,R(rk,rk) sao todas diferentes de zero. Isso produz uma ”economia de tamanho”: Se Aie m-por-n com m ¿ n, entao, apenas as n primeiras colunas de Q sao computadas tantoquanto as n priemiras linhas de R.

Como disse anteriormente o primeiro For percorre inversamente as colunas de”n-1”ate ”0”, o segundo for e para zerar os elementos de Q que serao computados poste-riormente. O Terceiro For faz o somatorio das multiplicacoes de ”QR*Q”que e guardadona varavel s e depois adicionado ao valor de ”Q[i][j].

public Matrix getQ () {Matrix X = new Matrix(m,n);double[][] Q = X.getArray();for (int k = n-1; k >= 0; k--) {//Ordem n

for (int i = 0; i < m; i++) {//Ordem n*nQ[i][k] = 0.0;

}Q[k][k] = 1.0;for (int j = k; j < n; j++) {//Ordem n*n

if (QR[k][k] != 0) {double s = 0.0;for (int i = k; i < m; i++) {//Ordem n*n*n

s += QR[i][k]*Q[i][j];}s = -s/QR[k][k];for (int i = k; i < m; i++) {//Ordem n*n*n

Q[i][j] += s*QR[i][k];}

}}

}return X;

}

3.3. Decomposicao SVD

Este algoritmo consiste em decompor A em USV T , onde Umm e V Tnn sao matrizes orto-

gonais, ou seja, UUT = Im e V V T = In e∑

mn e uma matriz diagonal que contem osvalores singulares σj em ordem nao crescente σ1 σ2... ≥ σj ≥ 0 para 1 ≤ j ≤ min(m,n)

public DecomposicaoValorSingular(Matrix Arg) {

// Inicializa.double[][] A = Arg.getArrayCopy();m = Arg.getRowDimension();n = Arg.getColumnDimension();

int nu = Math.min(m, n);s = new double[Math.min(m + 1, n)];U = new double[m][nu];V = new double[n][n];double[] e = new double[n];double[] trabalha = new double[m];boolean wantu = true;boolean wantv = true;

// Reduz A em uma forma bidiagonal, armazendando oselementos da

// diagonal em s e os elementos da diagonal superior em e.

int nct = Math.min(m - 1, n);int nrt = Math.max(0, Math.min(n - 2, m));for (int k = 0; k < Math.max(nct, nrt); k++) {//Ordem n

if (k < nct) {

// Computa a transformacao para a k-esimacoluna

// e coloca a k-esima diagonal em s[k].// Computa 2-normal da k-esima coluna sem

overflow.s[k] = 0;for (int i = k; i < m; i++) {//Ordem n*n

s[k] = hypot(s[k], A[i][k]);}if (s[k] != 0.0) {

if (A[k][k] < 0.0) {s[k] = -s[k];

}for (int i = k; i < m; i++) {

A[i][k] /= s[k];}A[k][k] += 1.0;

}s[k] = -s[k];

}for (int j = k + 1; j < n; j++) {//Ordem n*n

if ((k < nct) & (s[k] != 0.0)) {

// Aplica a transformacao

double t = 0;for (int i = k; i < m; i++)

{//Ordem n*n*n

t += A[i][k] * A[i][j];}t = -t / A[k][k];for (int i = k; i < m; i++)

{//Ordem n*n*nA[i][j] += t * A[i][k];

}}

// Coloca a k-esima lina de A em e para osubsequente

// calculo da trasnformacao da linha.e[j] = A[k][j];

}if (wantu & (k < nct)) {

// Coloca a transformacao em U para asubsequente multiplicacao

// retroativa.for (int i = k; i < m; i++) {

U[i][k] = A[i][k];}

}if (k < nrt) {

// Calcula a k-esima transforma o dalinha e coloca a

//k-esima diagonal superior em e[k].// Calcula segunda normal sem overflow.

e[k] = 0;for (int i = k + 1; i < n; i++) {//Ordem

n*ne[k] = hypot(e[k], e[i]);

}if (e[k] != 0.0) {

if (e[k + 1] < 0.0) {e[k] = -e[k];

}for (int i = k + 1; i < n; i++) {

e[i] /= e[k];}e[k + 1] += 1.0;

}e[k] = -e[k];if ((k + 1 < m) & (e[k] != 0.0)) {

// Aplica a transformacao.

for (int i = k + 1; i < m; i++){//Ordem n*n

trabalha[i] = 0.0;}for (int j = k + 1; j < n; j++)

{//Ordem n*nfor (int i = k + 1; i < m;

i++) {trabalha[i] +=

e[j] * A[i][j];}

}for (int j = k + 1; j < n; j++)

{//Ordem n*ndouble t = -e[j] / e[k +

1];for (int i = k + 1; i < m;

i++) {A[i][j] += t *

trabalha[i];}

}}if (wantv) {

// Coloca a transformacao em Vpara subsequente

//multiplicacao retroativa.

for (int i = k + 1; i < n; i++){//Ordem n*n

V[i][k] = e[i];}

}}

}

4. ConclusaoPodemos concluir que utilizando os conceitos de Autovetores e Autovalores podemosmanipular e decompor de maneiras diferentes as matrizes, o que possibilita obter ganho deprocessamento em calculos mais complexos. Devemos lembrar que essas decomposicoessao possıveis em casos especıficos onde as matrizes a serem trabalhadas apresentem aspropriedades como simetria e ou tambem simetrica definida positiva.

Com esse trabalho serviu para complementar o aprendizado da algebra linear

computacional no curso de Mestrado em ciencias computacionais, esperamos que essesprogramas sejam reutilizados e que inspirem a elaboracao de teses, nao apenas no escopoda algebra linear, mas tambem em outras situacoes onde seja possıvel aplicar conceitosmatematicos na solucao de problemas reais.

ReferenciasANTON, H. & BUSBY, R. (2006). Algebra Linear Contemporanea. BOOKMAN COM-

PANHIA EDITORA, 1th edition.

Oliveira, D. E. d. (2006). Sobre um metodo aseemelhado ao de francis para adeterminacao de autocalores de matrizes.

Rodrigues, C. F. (2011). Analise comparativa entre os metodos de composicaco em va-lores singulares e analise de componentes principais envolvendo matrizes esparsas degrande porte.

Stanford, A. (2002). Elemento de economia matematica i - conceitos fundamentais daAlgebra linear.