BCM0505-15 - Processamento da Informação
Aula A07 – 20/05: Matrizes
- Sistemas lineares e eliminação gaussiana.
- Cálculo de auto-valores e auto-vetores.
- Google’s page rank (breve descrição).
Produtos: Vetores e Matrizes
Problema A7.1: Dada uma matriz $A\in\mathbb{R}^{m\times n}$ e um vetor $x\in\mathbb{R}^n$, determine o vetor $z=Ax$.
Temos que $z\in\mathbb{R}^m$ e que pode ser calculado com uma variação do algoritmo de multiplicação de matrizes.
def matrix_prod_vector (A: [[float]], x: [float]) -> [float]:
m, n = len (A), len (A[0])
z = [0] * m
for i in range (m):
for j in range (n):
z[i] += A[i][j] * x[j]
return z
Uma versão com list comprehensions:
def matrix_prod_vector (A: [[float]], x: [float]) -> [float]:
return [sum (aij*xj for aij, xj in zip (ai, x)) for ai in A]
Problema A7.2: Dada uma matriz $A\in\mathbb{R}^{m\times n}$ e um vetor $y\in\mathbb{R}^m$, determine o vetor $w=y^\top A$.
Temos que $w\in\mathbb{R}^n$. Como $(y^\top A)^\top = A^\top y$, podemos simplesmente escrever:
def vector_prod_matrix (y: [float], A: [[float]]) -> [float]:
return matrix_prod_vector (matrix_transpose (A), y)
Apesar de correta, a solução acima precisa criar a matriz transposta $A^\top$ antes da multiplicação por $y$ e isso produz um código ligeiramente mais lento que o necessário na prática. Modificando levemente a primeira resposta ao Problema A7.1, obtemos um código direto.
def vector_prod_matrix (y: [float], A: [[float]]) -> [float]:
m, n = len (A), len (A[0])
w = [0] * n
for j in range (n):
for i in range (m):
w[j] += y[i] * A[i][j]
return w
Sem realizarmos a transposição de $A$, no entanto, não há como evitar a indexação explícita de $A$ em uma versão via list comprehensions:
def vector_prod_matrix (y: [float], A: [[float]]) -> [float]:
return [sum (y[i]*A[i][j] for i in range (len (A))) for j in range (len (A[0]))]
Observe que se $x=\mathbb{1}_n$ e $y=\mathbb{1}_m$, vetores de dimensões apropriadas em que todas as entradas são iguais a $1$, temos que $z[i]$ é a soma dos elementos da linha $i$ de $A$ e que $w[j]$ é a soma dos elementos da coluna $j$ de $A$.
Matrizes Triangulares
Uma matriz quadrada $U\in\mathbb{R}^{n\times n}$ é triangular superior se $U[i][j] = 0$ quando $j < i$. Em sendo o caso,
$U$ tem a seguinte estrutura:
$$U = \left(\begin{array}{ccccc}
u_{0,0} & u_{0,1} & u_{0,2} & \cdots & u_{0,n-1} \\
0 & u_{1,1} & u_{1,2} & \cdots & u_{1,n-1} \\
0 & 0 & u_{2,2} & \cdots & u_{2,n-1} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
0 & 0 & 0 & \cdots & u_{n-1,n-1}
\end{array}\right).$$
Ou seja, todas as entradas abaixo da diagonal principal de $U$ são iguais a zero. Dizemos que $U$ é estritamente triangular superior
se além de triangular superior tivermos que $u_{0,0} = u_{1,1} = \cdots = u_{n-1,n-1} = 0$. Uma matriz quadrada
$L\in\mathbb{R}^{n\times n}$ é (estritamente) triangular inferior se $L^\top$ é (estritamente) triangular superior.
Problema A7.3: Dada uma matriz quadrada $A\in\mathbb{R}^{n\times n}$, determine se $A$ é triangular superior.
def is_triang_sup (A: [[float]]) -> bool:
for i in range (1, len (A)):
for j in range (i):
if A[i][j] != 0:
return False
return True
Sistemas de Equações Lineares
Considere o seguinte
Problema A7.4: Sejam $A\in\mathbb{R}^{m\times n}$ uma matriz e $b\in\mathbb{R}^m$ um vetor. Determine um vetor $x\in\mathbb{R}^n$ tal que $$Ax=b.$$
Como enunciado, $A$ e $b$ são dados e $x=[x_0,x_1,\ldots,x_{n-1}]$ é um vetor de variáveis que tomam valores reais. $Ax=b$ é um sistema de $m$ equações em $n$ variáveis (incógnitas); o grau de cada variável $x_j$, em qualquer equação $i$, é igual a $1$. Logo, $Ax=b$ é um sistema de equações lineares ou, simplesmente, um sistema linear.
Talvez vocês estejam mais acostumados com a seguinte notação,
$$
\left.\begin{array}{ccccccccc}
a_{0,0} x_0 &+& a_{0,1} x_1 &+& \cdots &+& a_{0,n-1} x_{n-1} & = & b_0\\
a_{1,0} x_0 &+& a_{1,1} x_1 &+& \cdots &+& a_{1,n-1} x_{n-1} & = & b_1\\
\vdots & & \vdots & & \ddots & & \vdots & & \vdots\\
a_{m-1,0} x_0 &+& a_{m-1,1} x_1 &+& \cdots &+& a_{m-1,n-1} x_{n-1} & = & b_{m-1}
\end{array}\right.
$$
em que as equações são escritas explicitamente. Essa forma produz um enfoque natural na organização do sistema por linhas e às vezes também é
escrita como
$$
\left.\begin{array}{ccc}
a_0^\top x & = & b_0 \\
a_1^\top x & = & b_1 \\
\vdots & \vdots & \vdots \\
a_{m-1}^\top x & = & b_{m-1} \\
\end{array}\right.
$$
em que $a_i = [a_{i,0},a_{i,1},\ldots,a_{i,n-1}]$.
Uma variação consiste em organizar o sistema por colunas, denominada notação vetorial:
$$ x_0 \left(\begin{array}{c} a_{0,0} \\ a_{1,0} \\ \vdots \\ a_{m-1,0} \end{array}\right) +
x_1 \left(\begin{array}{c} a_{0,1} \\ a_{1,1} \\ \vdots \\ a_{m-1,1} \end{array}\right) +
\cdots +
x_{n-1} \left(\begin{array}{c} a_{0,n-1} \\ a_{1,n-1} \\ \vdots \\ a_{m-1,n-1} \end{array}\right) =
\left(\begin{array}{c} b_0 \\ b_1 \\ \vdots \\ b_{m-1} \end{array}\right).
$$
Fazendo
$$
A = \left(\begin{array}{cccc}
a_{0,0} & a_{0,1} & \cdots & a_{0,n-1} \\
a_{1,0} & a_{1,1} & \cdots & a_{1,n-1} \\
\vdots & \vdots & \ddots & \vdots \\
a_{m-1,0} & a_{m-1,1} & \cdots & a_{m-1,n-1}
\end{array}\right),
\qquad
x = \left(\begin{array}{c}
x_0 \\ x_1 \\ \vdots \\ x_{n-1}
\end{array}\right)
\qquad\text{e}\qquad
b = \left(\begin{array}{c}
b_0 \\ b_1 \\ \vdots \\ b_{m-1}
\end{array}\right),
$$
temos a notação matricial $Ax=b$, a mais sucinta das três. É conveniente alternarmos entre as três notações, pois cada uma
ilumina características diferentes do sistema.
O sistema $Ax=b$ pode ter uma única solução, ter múltiplas soluções, ou não ter solução alguma. A classificação depende dos valores de $m$, $n$, $a_{i,j}$ e $b_i$. Não é nosso intuito desenvolver a (bela e simples) teoria dos sistemas lineares sobre $\mathbb{R}$; o curso de Álgebra Linear trata deste e de outros assuntos relacionados. Logo, vamos tratar somente de alguns casos particulares do Problema A7.4. O primeiro deles consiste do chamado sistema diagonal.
Uma matriz quadrada $D\in\mathbb{R}^{n\times n}$ é diagonal se $D[i][j]=0$ quando $i\neq j$.
Problema A7.5: Dados uma matriz diagonal $D\in\mathbb{R}^{n\times n}$ tal que $D[i][i]\neq 0$ e um vetor $b\in\mathbb{R}^n$, determine um vetor $x\in\mathbb{R}^n$ tal que $Dx=b$.
O Sistema A7.5 tem a forma $$ \left(\begin{array}{cccc} d_{0,0} & 0 & 0 & 0 \\ 0 & d_{1,1} & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & d_{n-1,n-1} \end{array}\right) \cdot \left(\begin{array}{c} x_0 \\ x_1 \\ \vdots \\ x_{n-1} \end{array}\right) = \left(\begin{array}{c} d_{0,0}x_0 \\ d_{1,1}x_1 \\ \vdots \\ d_{n-1,n-1}x_{n-1} \end{array}\right) = \left(\begin{array}{c} b_0 \\ b_1 \\ \vdots \\ b_{n-1} \end{array}\right). $$ Como $D[i][i] = d_{i,i}\neq 0$, segue que $$x_i = \frac{b_i}{d_{i,i}}$$ é a solução (única, no caso). Em forma matricial, $x=D^{-1}b$, em que $D^{-1}D=I_n$.
def linsys_diagonal (D: [[float]], b: [float]) -> [float]:
return [b[i]/D[i][i] for i in range (len (D))]
Algumas considerações sobre o problema/código acima:
- (1) se sabemos que a matriz é diagonal, ela pode ser representada em um vetor; poupamos $n^2-n$ entradas iguais a zero com isso;
- (2) se permitirmos zeros na diagonal principal de $D$, o problema fica mais interessante:
- caso $d_{i,i} = b_i = 0$, o valor de $x_i$ pode ser qualquer; caso o sistema seja consistente, ele tem então infinitas soluções;
- caso $d_{i,i} = 0$ e $b_i \neq 0$, o sistema é inconsistente e não tem solução.
Exercício 1: escreva um código em Python que leve em consideração os pontos acima. Dica: além do vetor $x$, devolva um
bool
para indicar se o sistema é consistente.
Problema A7.6: Dados uma matriz triangular superior $U\in\mathbb{R}^{n\times n}$ tal que $U[i][i]\neq 0$ e um vetor $b\in\mathbb{R}^n$, determine um vetor $x\in\mathbb{R}^n$ tal que $Ux=b.$
O Sistema A7.6 tem a forma
$$\begin{array}{ccccccccccc}
u_{0,0}x_0 &+& u_{0,1}x_1 &+& u_{0,2}x_2 &+& \cdots &+& u_{0,n-1}x_{n-1} &=& b_0\\
& & u_{1,1}x_1 &+& u_{1,2}x_2 &+& \cdots &+& u_{1,n-1}x_{n-1} &=& b_1\\
& & & & u_{2,2}x_2 &+& \cdots &+& u_{2,n-1}x_{n-1} &=& b_2\\
& & & & & & \ddots & & \vdots & & \vdots \\
& & & & & & & & u_{n-1,n-1}x_{n-1} &=& b_{n-1}
\end{array}$$
denominada escalonada. Da última equação, obtemos que $x_{n-1} = b_{n-1} / u_{n-1,n-1}$, já que $u_{i,i}\neq 0$. Substituindo o
valor de $x_{n-1}$ na penúltima equação, obtemos que
$$x_{n-2} = \frac{b_{n-2}-u_{n-2,n-1}x_{n-1}}{u_{n-2,n-2}} = \frac{b_{n-2}-u_{n-2,n-1}b_{n-1} / u_{n-1,n-1}}{u_{n-2,n-2}}.$$
Não é difícil perceber que este processo pode ser repetido, de baixo para cima, até a primeira equação, com $x_{n-1}, x_{n-2},\ldots, x_0$
sendo determinados nesta ordem. Esse processo é chamado de retro-substituição. O passo genérico é dado por:
$$x_{i} = \frac{1}{u_{i,i}}\left(b_i - \sum_{j=i+1}^{n-1} u_{i,j}x_j\right),$$
com $x_j$ sendo resolvido antes de $x_i$ (o que implica em $i=n-1,n-2,\ldots,0$, nesta ordem).
def linsys_triangular (U: [[float]], b: [float]) -> [float]:
n = len (U)
x = [0] * n
for i in reversed (range (n)):
x[i] = (b[i] - sum (U[i][j]*x[j] for j in range (i+1, n))) / U[i][i]
return x
Por ser triangular superior e não possuir zeros na diagonal principal, temos que $U$ é não singular: existe uma matriz $U^{-1}$, dita inversa à esquerda de $U$, tal que $U^{-1}U=I_n$. Com isso, $Ux=b$ se e somente se $$x = I_n x = U^{-1}Ux=U^{-1}b,$$ e o sistema sempre tem uma única solução. Considerações semelhantes às feitas no caso diagonal aplicam-se caso seja permitido $u_{i,i}=0$.
Eliminação Gaussiana
Claramente, nem todas as matrizes $A\in\mathbb{R}^{m\times n}$ do Problema A7.4 possuem as características das matrizes $U\in\mathbb{R}^{n\times n}$ do Problema A7.6. Um dos resultados fundamentais em Álgebra Linear estabelece que $A$ pode ser transformada a ponto de ser “parecida” com $U$, desde que satisfaça algumas condições estruturais.
Vamos descrever a transformação e identificar as condições em matrizes quadradas. Sejam $A\in\mathbb{R}^{n\times n}$ e $b\in\mathbb{R}^n$ e defina a matriz estendida $$\big[A\big|b\big]=\left(\begin{array}{cccc|c} a_{0,0} & a_{0,1} & \cdots & a_{0,n-1} & b_0 \\ a_{1,0} & a_{1,1} & \cdots & a_{1,n-1} & b_1 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ a_{n-1,0} & a_{n-1,1} & \cdots & a_{n-1,n-1} & b_{n-1} \end{array}\right)$$ que codifica o sistema $Ax=b$.
Considere as seguintes três operações envolvendo linhas de $[A|b]$:
- troca de posições entre duas linhas;
- multiplicação uma linha por um escalar $c\in\mathbb{R}$ não nulo;
- soma do múltiplo de uma linha a uma outra linha.
Utilizando essas três operações adequadamente e um certo número de vezes, é possível transformar $[A|b]$ em uma matriz $[A'|b']$, em que $A'$ é triangular superior. Vamos fornecer um exemplo: $$\big[A\big|b\big]= \left(\begin{array}{rrr|r} 2 & 1 & -1 & 8 \\ -3 & -1 & 2 & -11 \\ -2 & 1 & 2 & -3 \end{array}\right) \simeq \left(\begin{array}{rrr|r} 2 & 1 & -1 & 8 \\ 0 & \frac{1}{2} & \frac{1}{2} & 1 \\ 0 & 2 & 1 & 5 \end{array}\right) \simeq \left(\begin{array}{rrr|r} 2 & 1 & -1 & 8 \\ 0 & \frac{1}{2} & \frac{1}{2} & 1 \\ 0 & 0 & -1 & 1 \end{array}\right)= [A'|b'], $$ em que a segunda matriz é obtida da primeira somando-se a segunda linha com $3/2$ da primeira e a terceira linha com a primeira, e a terceira matriz é obtida da segunda somando-se a terceira linha com $-4$ vezes a segunda.
Não é difícil provar que $Ax=b$ e $A’x=b'$ tem o mesmo conjunto de soluções. Logo, podemos utilizar a
função linsys_triangular
para determinar uma solução para o sistema acima. O processo de escalonamento ilustrado é conhecido como
eliminação gaussiana. Uma versão simples da eliminação para $[A|b]$ é descrita como:
para cada coluna $k = 0, 1, \ldots, n-2$
$\quad$seja $\mathsf{pivot} = A[k][k]$
$\quad$para cada linha $i = k+1, k+2, \ldots, n-1$
$\quad\quad$seja $c = -A[i][k]/\mathsf{pivot}$
$\quad\quad$faça $A[i][k] = 0$
$\quad\quad$para cada coluna $j = k+1, k+2, \ldots, n-1$
$\quad\quad\quad$faça $A[i][j] = A[i][j] + c\cdot A[k][j]$
$\quad\quad$faça $b[i] = b[i] + c\cdot b[k]$
É um bom exercício simular este pseudo-código para o exemplo $[A|b]$ acima.
O pseudo-código fornecido funciona muito bem quando a matriz $A$ tem posto completo, isto é, o número de linhas linearmente independentes em $A$ é igual ao número de linhas em $A$. Caso $A$ não tenha posto completo, divisões por zero podem acontecer. Isto é ilustrado em outro exemplo: $$\big[C\big|d\big]= \left(\begin{array}{rrr|r} 1 & 3 & 1 & 9 \\ 1 & 1 & -1 & 1 \\ 3 & 11 & 5 & 35 \end{array}\right) \simeq \left(\begin{array}{rrr|r} 1 & 3 & 1 & 9 \\ 0 & -2 & -2 & -8 \\ 0 & 2 & 2 & 8 \end{array}\right) \simeq \left(\begin{array}{rrr|r} 1 & 3 & 1 & 9 \\ 0 & -2 & -2 & -8 \\ 0 & 0 & 0 & 0 \end{array}\right)= [C'|d']. $$ Note que $[C'|d']$ — e, portanto, $C'$ — tem uma linha nula, pois $$[3,11,5,35] - 3\cdot[1,3,1,9] + ([1,1,-1,1] - [1,3,1,9]) = 0,$$ o que mostra que $[3,11,5,35] = 4\cdot[1,3,1,9] - [1,1,-1,1]$. Isto é, a terceira linha é uma combinação linear das outras duas.
Temos algumas observações importantes e úteis:
-
uma linha nula durante o processo de eliminação não tem pivô (pivot); deve ser ignorada ou o processo deve ser interrompido, dependendo do propósito; insistir em pivotação em tal linha produzirá uma divisão por $0$;
-
se o propósito for ignorar as linhas nulas e continuar a eliminação, é conveniente escolher, para cada coluna, a entrada da linha cujo valor seja o maior em módulo (fato conhecido como pivotação parcial); isto aumente a estabilidade do método;
-
a eliminação pode ser utilizada para determinar o posto da matriz $C$, que é igual a $n$ menos o número de linhas nulas em $C'$; se o posto de $C$ for menor que $n$, $C$ é singular (isto é, não invertível); quando $C$ é singular, seu determinante é igual a $0$;
-
o sistema $Cx=d$ tem infinitas soluções já que a última linha de $[C'|d']$ é nula; de forma similar, $[C'|d']$ poderia ter uma linha nula em $C'$ e não nula em $d'$ — se trocássemos o $35$ por $36$ em $d$, por exemplo — o que implicaria que $Cx=d$ não tem solução;
Uma versão melhor de eliminação gaussiana é:
para cada coluna $k = 0, 1, \ldots, n-2$
$\quad$seja $\ell = \textrm{argmax} \big\{|A[i][k]|: i = k, k+1, \ldots, n-1\big\}$
$\quad$se $A[\ell][k] = 0$, ignore a coluna $k$ e continue com o próximo passo do laço
$\quad$troque as linhas $(A[k][\cdot],b[k])\leftrightarrow (A[\ell][\cdot],b[\ell])$
$\quad$seja $\mathsf{pivot} = A[k][k]$
$\quad$para cada linha $i = k+1, k+2, \ldots, n-1$
$\quad\quad$seja $c = -A[i][k]/\mathsf{pivot}$
$\quad\quad$faça $A[i][k] = 0$
$\quad\quad$para cada coluna $j = k+1, k+2, \ldots, n-1$
$\quad\quad\quad$faça $A[i][j] = A[i][j] + c\cdot A[k][j]$
$\quad\quad$faça $b[i] = b[i] + c\cdot b[k]$
A codificação em Python do pseudo-código acima é imediata. Vamos fazer as transformações na própria matriz $A$ e vetor $b$.
def gaussian (A: [[float]], b: [float]):
n = len (A)
for k in range (n-1):
l = max (range (k, n), key = lambda i: abs (A[i][k]))
if A[l][k] == 0: continue
A[l], A[k] = A[k], A[l]
b[l], b[k] = b[k], b[l]
pivot = A[k][k]
for i in range (k+1, n):
c = -A[i][k] / pivot
A[i][k] = 0
for j in range (k+1, n):
A[i][j] += c * A[k][j]
b[i] += c * b[k]
Não é difícil perceber que o código acima consome tempo $O(n^3)$ no pior caso.
Calculando Determinantes
Formalmente, o determinante de $A$ é uma função dada por $$ \mathsf{det}(A) := \sum_{\sigma\in\widehat{\mathfrak{S}_n}}\left(\textrm{sgn}(\sigma)\prod_{i=0}^{n-1} A[i][\sigma[i]]\right), \qquad\qquad (*) $$ em que $\widehat{\mathfrak{S}_n}$ é o conjunto de todas as permutações de $N:=\{0,1,\ldots,n-1\}$, $\sigma$ é uma permutação de $N$, e $\mathrm{sgn}(\sigma)$ é igual a $1$ se o número de inversões em $\sigma$ é par e igual a $-1$ em caso contrário.
Se $A$ é triangular superior, temos que $$\mathsf{det}(A) = \prod_{i=0}^{n-1} A[i][i],$$ pois todos os demais produtos em $(*)$ são iguais a zero.
Observe que o número de termos no lado direito de $(*)$ é $n!$ no pior caso. Logo, é vantajoso transformar $A$ em uma equivalente triangular superior $A'$ para calcular o determinante de $A$ em tempo $O(n^3)$. Para tal, observe que:
- multiplicar uma linha de $A$ por um escalar $c\in\mathbb{R}$ implica em multiplicar $\mathsf{det}(A)$ por $c$;
- trocar duas linhas de $A$ de lugar implica em multiplicar $\mathsf{det}(A)$ por $-1$;
- adicionar um múltiplo escalar de uma linha a outra linha de $A$ mantem o valor de $\mathsf{det}(A)$ inalterado.
Mantendo o controle dessas operações durante uma eliminação gaussiana permite calcular o determinante de $A$ à partir do de $A'$.
Exercício 2: Escreva uma função que recebe uma matriz quadrada $A\in\mathbb{R}^{n\times n}$ e calcula $\mathsf{det}(A)$ usando o método de eliminação gaussiana.
Eliminação Completa
Voltemos ao exemplo $[A|b]$ fornecido acima. Mais precisamente, foquemos na forma triangular superior $$ [A'|b'] = \left(\begin{array}{rrr|r} 2 & 1 & -1 & 8 \\ 0 & \frac{1}{2} & \frac{1}{2} & 1 \\ 0 & 0 & -1 & 1 \end{array}\right). $$ Não é difícil perceber que podemos realizar uma eliminação gaussiana reversa (de baixo para cima, de trás para frente) em $[A'|b']$ e obtermos $$ [A'|b'] = \left(\begin{array}{rrr|r} 2 & 1 & -1 & 8 \\ 0 & \frac{1}{2} & \frac{1}{2} & 1 \\ 0 & 0 & -1 & 1 \end{array}\right) \simeq \left(\begin{array}{rrr|r} 2 & 1 & 0 & 7 \\ 0 & \frac{1}{2} & 0 & \frac{3}{2} \\ 0 & 0 & -1 & 1 \end{array}\right) \simeq \left(\begin{array}{rrr|r} 2 & 0 & 0 & 4 \\ 0 & \frac{1}{2} & 0 & \frac{3}{2} \\ 0 & 0 & -1 & 1 \end{array}\right) = [A''|b'']. $$ Finalmente (ou durante a eliminação reversa), podemos simplificar as entradas de $[A''|b'']$ e escrever $$ [A^*|b^*] = \left(\begin{array}{rrr|r} 1 & 0 & 0 & 2 \\ 0 & 1 & 0 & 3 \\ 0 & 0 & 1 & -1 \end{array}\right). $$ Temos que $A^*=I_n$, a matriz identidade: claramente, $A$ tem posto completo!
Exercício 3: Escreva uma função que implemente a eliminação gaussiana completa.
Caso a matriz $A$ tenha posto completo, a eliminação gaussiana completa permite determinarmos sua inversa $A^{-1}$ em tempo $O(n^3)$. E o que é ainda melhor: podemos descobrir se $A$ tem posto completo durante a determinação de $A^{-1}$. Para tal, basta fazermos a eliminação completa em $[A|I_n]$ até:
- encontrarmos uma linha nula, o que implica que $A$ é singular, ou
- obtermos $[I_n|B]$; neste caso, temos que $A\cdot B = B\cdot A = I_n$ e assim, $B$ é a inversa de $A$.
Nota: erros numéricos podem ocorrer na determinação de $B$ de forma que $A\cdot B \neq I_n$; não vamos entrar nos detalhes de como lidar com e evitar isso neste curso.
Exercício 4: Escreva uma função que recebe uma matriz quadrada real $A$ e calcula sua inversa $A^{-1}$ ou determina que $A$ é singular.
Auto-valores e Auto-vetores
Seja $A\in\mathbb{R}^{n\times n}$ uma matriz quadrada e não nula de valores reais. Dizemos que um vetor não nulo $x\in\mathbb{R}^n$ é um auto-vetor ou um vetor característico de $A$ se existe um escalar $\lambda\in\mathbb{R}$ tal que $$Ax = \lambda x.$$ Em sendo o caso, $\lambda$ é um auto-valor associado a $x$. Em outras palavras, $x$ é um auto-vetor de $A$ se $x$ é um múltiplo escalar de $Ax$ e $\lambda$ é o fator de escala aplicado por $A$ em $x$; quando $\lambda<0$, a direção de $x$ é revertida.
Por exemplo, se
$$
A = \left(\begin{array}{cc} 2 & 1 \\ 1 & 2 \end{array}\right)
\qquad\text{temos que}\qquad
x = \left(\begin{array}{c} 1 \\ 1 \end{array}\right)
$$
é um auto-vetor de $A$, pois
$$
Ax = \left(\begin{array}{c} 3 \\ 3 \end{array}\right) =
3\cdot \left(\begin{array}{c} 1 \\ 1 \end{array}\right) = 3x,
$$
com $\lambda=3$. Já $z=(0~~1)^\top$ não é um auto-vetor de $A$, pois $Az = (1~~2)^\top$ e não existe $\lambda$ real
(nem complexo!) tal que $(1~~2)^\top = \lambda\cdot(0~~1)^\top$.
A equação $Ax = \lambda x$ pode ser re-escrita como $$ (A-\lambda I_n)x = \boldsymbol{0}, $$ em que $\boldsymbol{0}$ é o vetor com $n$ entradas iguais a zero. Temos que esta equação tem uma solução não nula se e somente se $$ \mathsf{det}(A-\lambda I_n) = 0, $$ isto é, se e somente se $\lambda$ for uma raiz do polinômio mônico $P(\lambda)$ de grau $n$ do lado esquerdo. Pelo teorema fundamental da álgebra, $P(\lambda)$ tem $n$ raízes $\lambda_0, \lambda_1, \ldots, \lambda_{n-1}\in\mathbb{C}$, isto é, complexas (algumas podem ser reais). Isso implica que $$ \mathsf{det}(A-\lambda I_n) = (\lambda_0-\lambda)\cdot (\lambda_1-\lambda)\cdots (\lambda_{n-1}-\lambda). \qquad(\triangle) $$ Logo, $$\mathsf{det}(A) = \prod_{i=0}^{n-1} \lambda_i.$$ Em outras palavras, o determinante de $A$ é igual ao produto de seus auto-valores. Isto sugere a existência de uma matriz diagonal $\Lambda$ equivalente a $A$: $\Lambda[i][i] = \lambda_i$ e $\Lambda[i][j] = 0$ para $i\neq j$. De fato, seja $Q$ uma matriz quadrada tal que a $j$-ésima coluna de $Q$ é um auto-vetor associado a $\lambda_j$. Temos que $$ AQ = Q\Lambda, $$ o que resulta em $$ A = Q\Lambda Q^{-1} \qquad\text{ou}\qquad Q^{-1}AQ = \Lambda. $$
Além disso, segue de $(\triangle)$ que dos $n$ auto-valores de $A$, alguns (todos) podem não ser reais e alguns podem ser iguais — falamos em multiplicidades de auto-valores, da mesma forma que falamos de multiplicidades de raízes.
Exemplos de matrizes com todos os $n$ auto-valores reais e imaginários são as matrizes simétricas ($A=A^\top$) e anti-simétricas ($A=-A^\top$), respectivamente.
Auto-vetor Dominante
Suponha que $A$ tem todos os auto-valores $\lambda_0, \lambda_1, \ldots, \lambda_{n-1}$ reais, ordenados de forma que $$|\lambda_0| \geq |\lambda_1| \geq \cdots \geq |\lambda_{n-1}|.$$
Um auto-vetor associado a $\lambda_0$ é dito dominante pois para $\Lambda$ como acima, temos que
$$
\Lambda^d =
\left(\begin{array}{cccc}
\lambda_0^d & 0 & \cdots & 0 \\
0 & \lambda_1^d & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & \lambda_{n-1}^d
\end{array}\right) =
\lambda_0^d
\left(\begin{array}{cccc}
1 & 0 & \cdots & 0 \\
0 & \Big(\frac{\lambda_1}{\lambda_0}\Big)^d & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & \Big(\frac{\lambda_{n-1}}{\lambda_0}\Big)^d
\end{array}\right)
\overset{d\to\infty}{\longrightarrow}
\lambda
\left(\begin{array}{cccc}
1 & 0 & \cdots & 0 \\
0 & 0 & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & 0
\end{array}\right),
$$
pois $|\lambda_i/\lambda_0|<1$ para $i\neq 0$, e que
$$
\lambda = \begin{cases}
0 & \text{se}~ -1 < \lambda_0 < +1,\\
1 & \text{se}~ \lambda_0 = +1,\\
\searrow & \text{se}~ \lambda_0 = -1,\\
\pm\infty & \text{se}~ \lambda_0 > \pm 1,
\end{cases}
$$
em que $d\geq 0$ é um inteiro e o símbolo $\searrow$ significa que o limite não existe:
fica oscilando entre $1$ e $-1$ à medida que $d$ vai de par para impar, e vice-versa.
Certamente, o caso interessante é quando $\lambda_0=1$:
$$
A^d Q = \Lambda^d Q = \lambda_0 x_0.
$$
Claramente, o produto de $A^d$ por $Q$ não é necessário; $Q$ pode ser substituída por $x_0$.
Caso o limite exista, normalizando $A^d x_0$, conseguimos $\lambda_0=1$.
Power Method
O método iterativo a seguir, conhecido como power method (método de potência), calcula uma aproximação para um auto-vetor
dominante de $A$. Seja $x_0[0\,..n-1]$ um vetor não nulo (uma estimativa) inicial e seja $m > 0$ o número de
passos a serem realizados. Para $A$, $x$ e $m$ como enunciados, faça:
$$
x_{k+1} = \frac{Ax_k}{\|Ax_k\|}
\qquad\text{enquanto}\qquad
k \leq m.
\qquad\qquad (\diamond)
$$
A função que calcula o produto de $A$ com $x_k$ foi feita no início da aula. Repetimos seu código abaixo por questões de conveniência.
def matrix_prod_vector (A: [[float]], x: [float]) -> [float]:
return [sum (aij*xj for aij, xj in zip (ai, x)) for ai in A]
A função que calcula a norma euclidiana $\|\cdot\|$ de um vetor foi feita antes do ECE, na aula 05 . Novamente, por questões de conveniência, exibimos seu código a seguir.
from math import sqrt
def euclid_norm (z: [float]) -> float:
return sqrt (sum (zi*zi for zi in z))
Problema A7.7: Escreva uma função que recebe $A$, $x_0\neq\vec{\boldsymbol{0}}$ e $m$ e devolve, se possível, uma aproximação para um auto-vetor dominante de $A$ de acordo com a Equação $(\diamond)$.
def power_method (A: [[float]], x: [float], m: int) -> (bool, [float]):
n = len (A)
for k in range (m):
x = matrix_prod_vector (A, x)
c = euclid_norm (x)
if c == 0: return (False, x)
for j in range (n):
x[j] /= c
return (True, x)
Google’s PageRank
No caso de $A$ ser estocástica e aperiódica, temos que $\lambda_0=1$. O auto-vetor dominante $x_0$, com $\|x_0\|=1$ é tal que $x_0[i]$ representa a probabilidade de, começando-se em um site aleatório da web, seguindo-se aleatoriamente os links das páginas, chagar-se à página $i$ no equilíbrio. Logo, os sites “mais populares” têm maiores valores em $x_0$.
Forneceremos mais detalhes após introduzirmos o conceito de grafos e digrafos na aula A09.