Problemas de autovalores e autovetores no Python

Depois de termos falado sobre sistemas lineares, este é o último texto da série sobre álgebra linear com Python, que trata da solução de problemas de autovalores.

Algumas definições

Em algumas direções, uma transformação linear se comporta como multiplicação por escalar. Estas direções e seus respectivos fatores multiplicativos são propriedade da transformação (isto é, da matriz) e são chamados de autovetores e autovalores, respectivamente.

Ou seja, se $A$ for uma matriz quadrada, um vetor não-nulo $v$ é autovetor de $A$ com autovalor $\sigma$ se

$$A v = \sigma v$$

Se rearranjarmos a equação, podemos observar que $v$ é solução do sistema de equações homogêneo

$$(A - \sigma I) v = 0$$

onde $I$ é a matriz identidade de mesma dimensão de $A$. Soluções não-nulas só vão existir se $A - \sigma I$ representar uma transformação singular, ou seja, se $p(\sigma) = \det(A - \sigma I) = 0$. O polinômio $p(\sigma)$ chamado de polinômio característico de $A$.

In [1]:
import numpy as np

Os métodos automáticos: numpy.linalg.eig*

No futuro pretendo escrever sobre obtenção de raízes de polinômios mas, na prática, obter autovalores dessa maneira é ineficiente. O pacote NumPy de álgebra linear para Python é capaz de encontrar autovalores e autovetores de matrizes de maneira muito mais eficiente.

Primeiro, vamos definir uma matriz e, depois, calcular apenas seus autovalores com numpy.linalg.eigvals:

In [2]:
A = np.array([[1., 0.], [0., 4.]])
A
Out[2]:
array([[1., 0.],
       [0., 4.]])
In [3]:
sigma = np.linalg.eigvals(A)
sigma
Out[3]:
array([1., 4.])

Como a matriz acima é diagonal, seus autovalores correspondem às entradas da diagonal. Vamos calcular agora para uma matriz não-diagonal:

In [4]:
B = np.array([[1., 2.], [3., 4.]])
B
Out[4]:
array([[1., 2.],
       [3., 4.]])
In [5]:
sigma = np.linalg.eigvals(B)
sigma
Out[5]:
array([-0.37228132,  5.37228132])

Para matrizes hermitianas (ou simétricas) existe uma função especializada mais eficiente, numpy.linalg.eigvalsh:

In [6]:
S = np.array([[2., 3.], [3., 4.]])
S
Out[6]:
array([[2., 3.],
       [3., 4.]])
In [7]:
np.linalg.eigvalsh(S)
Out[7]:
array([-0.16227766,  6.16227766])

A diferença entre as duas funções acima está (i) nos valores que acessam (a numpy.linalg.eigvalsh só precisa usar metade dos elementos da matriz, já que ela é simétrica) e, principalmente, (ii) na eficiência. Abaixo comparamos os tempos de execução para ambas as funções para uma matriz simétrica um pouco maior:

In [8]:
B = np.block([[np.eye(2), np.zeros([2, 2])],
              [np.zeros([2, 2]), S]])
B
Out[8]:
array([[1., 0., 0., 0.],
       [0., 1., 0., 0.],
       [0., 0., 2., 3.],
       [0., 0., 3., 4.]])
In [9]:
%%timeit
np.linalg.eigvalsh(B)
9.22 µs ± 373 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
In [10]:
%%timeit
np.linalg.eigvals(B)
27.3 µs ± 1.43 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Dá logo para perceber que, quando a matriz for simétrica ou hermitiana, vale a pena usar a função especializada: o tempo de execução pode chegar a ser duas vezes menor. Outra diferença é que, como matrizes simétricas ou hermitianas possuem sempre autovalores reais, funções especializadas como numpy.linalg.eigvalsh irão sempre garantir que o resultado seja não-complexo. (Se tivéssemos usado numpy.linalg.eigvals e tivéssemos certeza de que nossos autovalores são todos reais, poderíamos ter tomado apenas a parte real do vetor fazendo sigma.real.)

Mas e os autovetores? O conjunto de vetores com seus autovalores correspondentes pode ser obtido com a função numpy.linalg.eig (ou numpy.linalg.eigh, se a matriz for simétrica ou hermitiana):

In [11]:
# se é diagonal, é simétrica
sigma, U = np.linalg.eigh(A)
display(sigma, U)
array([1., 4.])
array([[1., 0.],
       [0., 1.]])

Ambas funções numpy.linalg.eig e numpy.linalg.eigh retornam dois objetos, um array unidimensional e um bidimensional. O primeiro, um vetor (chamado aqui de $\sigma$), contém os autovalores. O segundo, uma matriz (aqui chamado de $U$), contém os respectivos autovetores, um autovetor por coluna, na posição correspondente do seu autovalor em $\sigma$. Vale também observar que os autovetores são retornados normalizados (normas euclidianas iguais a um).

Como dito, a ordem deles "bate", de forma que o autovalor $\sigma_0$ (sigma[0]) está associado ao autovetor apresentado na zeroésima coluna de $U$ ($U_{:0}$ ou U[:, 0]):

In [12]:
U[:, 0]
Out[12]:
array([1., 0.])

(Observe que, quando extraímos a primeira coluna de $U$, obtemos um array unidimensional, o que é típico do NumPy.)

Agora fica fácil verificar que a condição de autovalor/autovetor é satisfeita para $\sigma_0$ e $U_{:0}$. Simplesmente realizamos ambas multiplicações $A v$ e $\sigma v$ da condição e verificamos se ambas resultam iguais. A verificação disso fica fácil com numpy.allclose:

In [13]:
np.allclose(A @ U[:, 0], sigma[0] * U[:, 0])
Out[13]:
True

Ambos resultados são portanto o mesmo, como devem ser. É claro que o mesmo ocorre com o segundo par $\sigma_1$ e $U_{:1}$:

In [14]:
np.allclose(A @ U[:, 1], sigma[1] * U[:, 1])
Out[14]:
True

Uma propriedade muito útil das matrizes simétricas (ou hermitianas) é que seus autovetores são em geral ortogonais (quando não degenerados). Vamos verificar isso abaixo para a matriz $S$, para a qual havíamos calculado apenas autovalores:

In [15]:
sigma, U = np.linalg.eigh(S)
display(sigma, U)
array([-0.16227766,  6.16227766])
array([[-0.81124219,  0.58471028],
       [ 0.58471028,  0.81124219]])
In [16]:
U[:, 0] @ U[:, 1]
Out[16]:
0.0

O produto interno entre autovetores associados a autovalores distintos é zero, o que caracteriza ortogonalidade. Já que os autovetores também estão normalizados, isso quer dizer que $U$ é unitária, ou seja,

$$U^T U = I$$
In [17]:
np.allclose(U.T @ U, np.eye(2))
Out[17]:
True

Por fim, em uma comparação análoga à anterior, observe o efeito no tempo de execução ao se utilizar a função numpy.linalg.eigh ao invés de numpy.linalg.eig para matrizes simétricas ou hermitianas:

In [18]:
%%timeit
sigma, U = np.linalg.eigh(S)
12 µs ± 904 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
In [19]:
%%timeit
sigma, U = np.linalg.eig(S)
34 µs ± 2.86 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Reordenando autovalores e autovetores

Os autovalores retornados de numpy.linalg.eigh ou numpy.linalg.eigvalsh são sempre retornados em ordem ascendente e repetidos de acordo com suas multiplicidades. No entanto, as funções numpy.linalg.eig e numpy.linalg.eigvals não garantem ordenamento dos autovalores, mas existe uma maneira de ordená-los e corrigir a ordem dos autovetores de acordo usando numpy.argsort.

Vamos gerar uma matriz aleatória $4 \times 4$ usando numpy.random.rand, calcular seus autovalores e autovetores, e reordená-los em ordem ascendente:

In [20]:
A = np.random.rand(4, 4)
sigma, U = np.linalg.eig(A)
display(sigma)

indices = np.argsort(sigma)
sigma = sigma[indices]
U = U[:, indices]
display(sigma)
array([ 1.91515211+0.j        ,  0.38355489+0.20195601j,
        0.38355489-0.20195601j, -0.2697303 +0.j        ])
array([-0.2697303 +0.j        ,  0.38355489-0.20195601j,
        0.38355489+0.20195601j,  1.91515211+0.j        ])

Diagonalização

Uma matriz quadrada é diagonalizável se ela é similar à uma matriz diagonal. Ou seja, $A$ é diagonalizável se existe uma matriz invertível $V$ tal que

$$\Sigma = V^{-1} A V$$

seja diagonal. Um resultado interessante da álgebra linear é que uma matriz quadrada $A$ de tamanho $n$ é diagonalizável se e somente se ela possui $n$ autovetores linearmente independentes. As colunas de $V$ terão então os autovetores de $A$ e $\Sigma$ terá os autovalores na diagonal.

Portanto, se tivermos todos os autovalores e autovetores de uma matriz, toda a informação para sabermos se ela é diagonalizável e toda a informação para a diagonalização já estará disponível. Podemos construir a matriz diagonal a partir dos autovalores com numpy.diag.

In [21]:
np.diag(sigma)
Out[21]:
array([[-0.2697303 +0.j        ,  0.        +0.j        ,
         0.        +0.j        ,  0.        +0.j        ],
       [ 0.        +0.j        ,  0.38355489-0.20195601j,
         0.        +0.j        ,  0.        +0.j        ],
       [ 0.        +0.j        ,  0.        +0.j        ,
         0.38355489+0.20195601j,  0.        +0.j        ],
       [ 0.        +0.j        ,  0.        +0.j        ,
         0.        +0.j        ,  1.91515211+0.j        ]])

Além disso, se os autovetores forem ortonormais, podemos nos valer do fato de que $U^{-1} = U^T$ (ou seja, como visto acima, $U$ será unitária):

In [22]:
np.allclose(np.linalg.inv(U), U.T)
Out[22]:
False
In [23]:
U.T @ np.diag(sigma) @ U
Out[23]:
array([[ 0.54099418+0.17689684j, -0.1086356 +0.03839203j,
        -0.22377567-0.09720118j, -0.01058118-0.07770512j],
       [-0.1086356 +0.03839203j,  0.08516411-0.24702832j,
         0.31191995-0.06718322j,  0.13681642+0.21232781j],
       [-0.22377567-0.09720118j,  0.31191995-0.06718322j,
         0.13608851+0.03815938j,  0.19298666-0.34693977j],
       [-0.01058118-0.07770512j,  0.13681642+0.21232781j,
         0.19298666-0.34693977j,  0.31502865-0.03682457j]])

Também podemos usar este truque para construírmos matrizes com autovalores e autovetores desejados, fazendo apenas a operação inversa ($C = V \Sigma V^{-1}$):

In [24]:
V = np.array([[1, 1, 0], [1, -1, 0], [1, 1, 1]])
V
Out[24]:
array([[ 1,  1,  0],
       [ 1, -1,  0],
       [ 1,  1,  1]])
In [25]:
sigma = np.array([3, 1, 2])
sigma
Out[25]:
array([3, 1, 2])
In [26]:
C = V @ np.diag(sigma) @ np.linalg.inv(V)
C
Out[26]:
array([[2., 1., 0.],
       [1., 2., 0.],
       [0., 1., 2.]])

Se calcularmos os autovalores e autovetores de C, veremos que obtivemos êxito em construir uma matriz com exatamente os autovalores e autovetores desejados:

In [27]:
sigma, V = np.linalg.eig(C)
display(sigma, V)
array([2., 3., 1.])
array([[ 0.        ,  0.57735027,  0.57735027],
       [ 0.        ,  0.57735027, -0.57735027],
       [ 1.        ,  0.57735027,  0.57735027]])

Observe acima que (i) a ordem dos autovalores não é necessáriamente a que usamos (embora possamos reordená-los como visto anteriormente) e que (ii) os autovetores diferem por terem sido retornados normalizados (embora preservem colinearidade com os que havíamos escolhido).

Potenciação de matrizes por diagonalização

Seja $A$ uma matriz quadrada. Calcular as potências de $A$ por multiplicação matricial

$$A^k = \overbrace{A A \cdots A}^{k \text{ multiplicações}}$$

é ineficiente. Já vimos em outro texto, por exemplo, que a função numpy.linalg.matrix_power não faz todas as multiplicações intermediárias e é mais rápida. Mas lá não havíamos dito como isso funciona. Aqui vamos ver uma maneira eficiente de potenciar matrizes, semelhante à implementada em numpy.linalg.matrix_power.

Se tivermos diagonalizado a matriz $A = V \Sigma V^{-1}$, podemos fazer

$$A^k = \left(V \Sigma V^{-1}\right)^k = \overbrace{V \Sigma V^{-1} V \Sigma V^{-1} \cdots V \Sigma V^{-1}}^{k \text{ multiplicações}} = V \Sigma^k V^{-1}$$

Ou seja, $V$ também diagonaliza a $k$-ésima potência de $A$, cujos autovalores são a $k$-ésima potência dos de $A$.

Vamos gerar uma matriz aleatória $4 \times 4$ usando numpy.random.rand, calcular $A^{50}$ das duas maneiras e medir o tempo de execução. Depois vamos comparar com o tempo de cálculo de numpy.linalg.matrix_power.

In [28]:
A = np.random.rand(4, 4)
A
Out[28]:
array([[0.63793558, 0.73715957, 0.89984734, 0.70444143],
       [0.94404201, 0.22125809, 0.60741571, 0.56164526],
       [0.35590991, 0.71898612, 0.60455103, 0.05662982],
       [0.64455821, 0.03599551, 0.0969152 , 0.52545749]])
In [29]:
sigma, V = np.linalg.eig(A)
display(sigma, V)
array([ 2.14018072,  0.44379691, -0.36299497, -0.23178048])
array([[-0.66012609,  0.00078128, -0.33430846, -0.09979794],
       [-0.54561856, -0.08857098,  0.77481098, -0.73386412],
       [-0.41955082,  0.6578778 , -0.46813645,  0.67106746],
       [-0.30085054, -0.74789797,  0.26220997,  0.03394541]])

(Observe que, para a matriz obtida, tanto os autovalores quanto os autovetores possuem entradas complexas.)

O código abaixo se vale do fato de que a potenciação de uma matriz diagonal é equivalente à potenciação individual das entradas da diagonal:

In [30]:
%%timeit
V @ np.diag(sigma**50) @ np.linalg.inv(V)
61.2 µs ± 5.43 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [31]:
%%timeit
A @ A @ A @ A @ A @ A @ A @ A @ A @ A \
@ A @ A @ A @ A @ A @ A @ A @ A @ A @ A \
@ A @ A @ A @ A @ A @ A @ A @ A @ A @ A \
@ A @ A @ A @ A @ A @ A @ A @ A @ A @ A \
@ A @ A @ A @ A @ A @ A @ A @ A @ A @ A
70.5 µs ± 1.26 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Fica claro que diagonalização computa potenciação muito mais rapidamente. Porém vale comentar que não estamos levando em consideração o tempo necessário para se computar $V$ e $\Sigma$. Se levarmos isso em conta, o método não é assim tão competitivo:

In [32]:
%%timeit
sigma, V = np.linalg.eig(A)
V @ np.diag(sigma**50) @ np.linalg.inv(V)
176 µs ± 24.3 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

No entanto, nenhum deles chega a ganhar da função numpy.linalg.matrix_power:

In [33]:
%%timeit
np.linalg.matrix_power(A, 50)
16 µs ± 255 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Talvez eu escreva um texto no futuro sobre os motivos pelos quais funções como a numpy.linalg.matrix_power são tão rápidas. Mas, por hora, completamos essa série sobre álgebra linear e Python. O próximo texto será o primeiro de uma série sobre gráficos e visualização de dados.

Deixe uma resposta