Como resolver 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 AA for uma matriz quadrada, um vetor não-nulo vv é autovetor de AA com autovalor σ\sigma se

Av=σv A v = \sigma v

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

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

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

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:

A = np.array([[1., 0.],
              [0., 4.]])
A
array([[1., 0.],
       [0., 4.]])
sigma = np.linalg.eigvals(A)
sigma
array([1., 4.])

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

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

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

S = np.array([[2., 3.],
              [3., 4.]])
S
array([[2., 3.],
       [3., 4.]])
np.linalg.eigvalsh(S)
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:

B = np.block([[np.eye(2), np.zeros([2, 2])],
              [np.zeros([2, 2]), S]])
B
array([[1., 0., 0., 0.],
       [0., 1., 0., 0.],
       [0., 0., 2., 3.],
       [0., 0., 3., 4.]])
%%timeit
np.linalg.eigvalsh(B)
6.7 µs ± 42.7 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
%%timeit
np.linalg.eigvals(B)
21.3 µs ± 225 ns 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):

# 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 UU), 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 σ0\sigma_0 (sigma[0]) está associado ao autovetor apresentado na zeroésima coluna de UU (U:0U_{:0} ou U[:, 0]):

U[:, 0]
array([1., 0.])

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

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

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

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

np.allclose(A @ U[:, 1], sigma[1] * U[:, 1])
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 SS, para a qual havíamos calculado apenas autovalores:

sigma, U = np.linalg.eigh(S)
display(sigma, U)
array([-0.16227766,  6.16227766])
array([[-0.81124219,  0.58471028],
       [ 0.58471028,  0.81124219]])
U[:, 0] @ U[:, 1]
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 UU é unitária, ou seja,

UTU=I U^T U = I
np.allclose(U.T @ U, np.eye(2))
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:

%%timeit
sigma, U = np.linalg.eigh(S)
8.9 µs ± 37.9 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
%%timeit
sigma, U = np.linalg.eig(S)
25.4 µs ± 63.6 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
2020-9-25. Compare com os tempos de execução da linguagem Julia para os mesmos problemas. Julia é pelo menos 2×2\times mais rápida em todos os casos.

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×44 \times 4 usando numpy.random.rand, calcular seus autovalores e autovetores, e reordená-los em ordem ascendente:

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.6625206 +0.j        ,  0.56711713+0.j        ,
       -0.25752503+0.22642322j, -0.25752503-0.22642322j])
array([-0.25752503-0.22642322j, -0.25752503+0.22642322j,
        0.56711713+0.j        ,  1.6625206 +0.j        ])

Diagonalização

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

Σ=V1AV \Sigma = V^{-1} A V

seja diagonal. Um resultado interessante da álgebra linear é que uma matriz quadrada AA de tamanho nn é diagonalizável se e somente se ela possui nn autovetores linearmente independentes. As colunas de VV terão então os autovetores de AA 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.

A = np.random.rand(4, 4)
A = A + A.T  # uma matriz simétrica
sigma, U = np.linalg.eigh(A)
Sigma = np.diag(sigma)
Sigma
array([[-0.85209919,  0.        ,  0.        ,  0.        ],
       [ 0.        ,  0.20935565,  0.        ,  0.        ],
       [ 0.        ,  0.        ,  1.03791616,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  3.98305659]])

Além disso, se os autovetores forem ortonormais (que de fato são para matrizes simétricas), podemos nos valer do fato de que U1=UTU^{-1} = U^T (ou seja, como visto acima, UU será unitária):

np.allclose(np.linalg.inv(U), U.T)
True
np.allclose(U @ Sigma @ U.T, A)
True

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

V = np.array([[1,  1, 0],
              [1, -1, 0],
              [1,  1, 1]])
V
array([[ 1,  1,  0],
       [ 1, -1,  0],
       [ 1,  1,  1]])
sigma = np.array([3, 1, 2])
sigma
array([3, 1, 2])
C = V @ np.diag(sigma) @ np.linalg.inv(V)
C
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:

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 AA uma matriz quadrada. Calcular as potências de AA por multiplicação matricial

Ak=AAAk multiplicaço˜es 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ΣV1A = V \Sigma V^{-1}, podemos fazer

Ak=(VΣV1)k=VΣV1VΣV1VΣV1k multiplicaço˜es=VΣkV1 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, VV também diagonaliza a kk-ésima potência de AA, cujos autovalores são a kk-ésima potência dos de AA.

Vamos gerar uma matriz aleatória 4×44 \times 4 usando numpy.random.rand, calcular A50A^{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.

A = np.random.rand(4, 4)
A
array([[0.86284278, 0.49848683, 0.77379656, 0.81607461],
       [0.88038112, 0.26791158, 0.39988671, 0.25028634],
       [0.05926472, 0.08051984, 0.04953885, 0.35667252],
       [0.37542371, 0.4042028 , 0.34657704, 0.74787015]])
sigma, V = np.linalg.eig(A)
display(sigma, V)
array([ 1.87574383+0.j        ,  0.09405826+0.20042815j,
        0.09405826-0.20042815j, -0.13569699+0.j        ])
array([[-0.72069797+0.j        ,  0.09149586+0.22918361j,
         0.09149586-0.22918361j,  0.43807158+0.j        ],
       [-0.50001752+0.j        ,  0.8093501 +0.j        ,
         0.8093501 -0.j        , -0.72470162+0.j        ],
       [-0.13541237+0.j        , -0.38053213+0.05450594j,
        -0.38053213-0.05450594j, -0.42942832+0.j        ],
       [-0.46069557+0.j        , -0.27604261-0.24511359j,
        -0.27604261+0.24511359j,  0.31383463+0.j        ]])

(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:

%%timeit
V @ np.diag(sigma**50) @ np.linalg.inv(V)
25.2 µs ± 137 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%%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
77.1 µs ± 186 ns 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 VV e Σ\Sigma. Se levarmos isso em conta, o método não é assim tão competitivo:

%%timeit
sigma, V = np.linalg.eig(A)
V @ np.diag(sigma**50) @ np.linalg.inv(V)
64.9 µs ± 330 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

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

%%timeit
np.linalg.matrix_power(A, 50)
15.5 µs ± 115 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
2020-9-25. Compare com os tempos de execução da linguagem Julia para os mesmos problemas. Julia é pelo menos 4×4\times mais rápida em todos os casos.

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.