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.

Avatar
Aritanan Gruber
Assistant Professor

“See, if y’all haven’t the same feeling for this, I really don’t give a damn. If you ain’t feeling it, then dammit this ain’t for you!"
(desconheço a autoria; agradeço a indicação)