Quociente de Iteração de Rayleigh

Question

Implemente o quociente de iteração de Rayleigh para computar o eigenvalue e o correspondente eigenvector de uma matriz. Teste seu programa na matriz do exercício anterior (4.3) usando um vetor inicial aleatório. (Livro, página 211, questão 4.4)

REFERÊNCIA DA QUESTÃO: HEATH, M. T. "Scientific Computing: An Introductory Survey". 2ª edição, ed.: McGraw-Hill. Nova York, 2002. Disponível em: < http://www-e6.ijs.si/~roman/files/tmp/M.Heath-SComputing/scientific-computing-michael-t-heath.pdf >

Comments

Excelente resposta, Luiz !

Se me permite, gostaria de trazer 2 contribuições a afim de melhor entender o quociente de iteração de Rayleigh.

A primeira trata-se do código em python, especificamente sobre o processo de iteração.
Faltou indentar o valor atualizado de mu, isto é. atualizar o autovalor na iteração, mediante sua inserção no for loop.
Aproveitando a deixa, segue uma proposta de reescrita do seu código.

def find_eigva(A,ite): 
# A square matrix
    # ite number of iteration

    v = randn(A.shape[0]) #eigvector
    #Guess
    mu = 1. + 1j #eigvalue

    #Iteration guess
    for k in range (ite):
        B = A - mu*eye(A.shape[0]) 
        try:
            w = solve(B,v) # find eigenvalue by characteristic polynomial
        except:
            print ("Converge solution")
            break

        """Rayleigh quotient"""

        v = w / norm(w,2) #normalize candidate eigenvector 
        mu = dot(conj(v.T), dot(A,v)) # Rayleigh quotient


    error= np.max(A - mu*eye(3))


    return mu.real, error


if __name__ == '__main__':
        A = np.array([[2,0,0],[0,3,4],[0,4,9]])
        print(find_eigva(A,10))

A segunda trata-se da parte matemática do problema.
Como explorado em seu código, especificamente na função solve(), é interessante lembrar que os autovalores podem ser interpretados como a raiz do polinômio característico definido pela matriz A, isto é, \( \lambda \) é um escalar tal que,
\[ p(\lambda) - det (A -\lambda I) = 0\]

Além disso, vale entender a construção do algoritmo (iteração + coeficiente de Rayleigh)
Entendendo o algoritmo

1. O algoritmo parte de um guess inicial para o autovalor, especificamente um complexo qualquer (definido no código por mu= 1+ 1j).
   Note que, se mu é um autovalor exato, então existe um número \( \rho \) tal que
   \[ A\textbf{q}=\rho\textbf{q} \]
   ou seja, \( \rho \) é um autovalor de q.
2. Reordenando a igualdade acima, temos o erro residual da iteração, isto é, para r o erro residual da iteração, segue que
   \[ r=A\textbf{q} -\rho\textbf{q} \]

Nesse sentido, podemos buscar encontrar \( \rho \) tal que \( \lVert \mathbf{r} \rVert \) seja mínimo, dado que quando q é um autovetor, ao minimizar \( \lVert \mathbf{r} \rVert \), \( \rho \) será exatamente o autovalor associado.

1. Sabemos que, para x uma solução do problema de mínimos quadrados, em termos do sistema \( A\textbf{b}=\textbf{b} \),isto é, para x que minimiza \( \lVert \mathbf{Ax-b} \rVert \), então existe \( A^*\), não nulotal que,
   \( A^{*}(A\textbf{x}-\textbf{b})=0 \) (condição de ortogonalidade)

De forma equivalente,
\[ A^*A\textbf{x}=A\textbf{b} \]

Nesse sentido, para o sistema \( \textbf{q} \rho=Aq \), tal que \( \rho\) é uma solução, existe **q*** tal que,

\[ (\textbf{q}^{*} \textbf{q}) \rho=A\textbf{q}\]

Logo,
\[\rho =\frac{\textbf{q}^*A\textbf{q}}{\textbf{q}^*\textbf{q}}\]

O \( \rho\) é chamado de Quociente de Rayleigh de q com respeito a A.

Note que no código o vetor q normalizado [v=w/norm(w,2)], isto é, \( \lVert \mathbf{r} \rVert=1\).

Nesse sentido, podemos reduzir \( \rho \) a
\[ \rho=\textbf{q}^*Aq \]

1. Agora vem o pulo do gato, qual a razão da escolha do quociente de Rayleigh?

1º) Como dissemos anteriormente, se q é um autovetor de A com autovalor \( \lambda \) associado, então \( \lambda \) será o coeficiente de Rayleigh. De fato,
\[ \rho=\frac{\textbf{q}*Aq}{\textbf{q*q}}=\frac{\textbf{q}*\lambda q}{\textbf{q*q}}=\lambda\]

2º) (O quociente de Rayleigh é a única solução). Teorema. Sejam \( A \in C^{n x m} \)e \( q \in C^n \). O único número complexo que minimiza \( \lVert \mathbf{\textbf{Aq}-\rho q}\rVert\) é o Quociente de Rayleigh.

3º) (O quociente de Rayleigh aproxima-se suficientemente de um autovalor, se q se aproxima de um autovetor).Teorema. Sejam \( A \in C^{nxm} \) e v um autovetor de A com autovalor \( \lambda \) associado. Assuma \( \lVert \mathbf{v} \rVert=1 \). Sejam \( q \in C^n \) com \( \lVert \mathbf{q} \rVert=1 \) e \( \rho=q^*Aq \) o Quociente de Rayleigh correspondente a q. Então,
\[ |\lambda - \rho| \leq 2 \lVert \mathbf{A}\rVert\ \lVert \mathbf{v-q} \rVert \]

Referência: https://repositorio.ufmg.br/bitstream/1843/EABA-978HQW/1/monografia_hortensia_americo.pdf

---

Obrigado pela excelente contribuição, Arthur. Já fiz a alteração no código você apontou.

Answer 1

Antes de adentrar na questão, é válido relembrar os conceitos de eigenvector e eigenvalue. O eigenvector (ou vetor característico) é um vetor "\(u\)" que, ao ser aplicado em um transformação linear \(T\) (ou seja, ao multiplicar pela direita uma matriz \(T\): \(Tu\)) retorna o próprio vetor "u" escalonado por um escalar c (tal escalar é o eigenvalue). Uma questão de interesse é investigar quem são os eigenvectors e eigenvalues de uma determinada matriz. Uma forma de conduzir tal investigação é a partir do uso do quociente de Rayleigh.

O quociente de Rayleigh é usado no teorema do min-max e em outros algoritmos com o intuito de encontrar os valores dos eigenvalues (exatos ou aproximados) de uma determinada matriz. Ele associa uma quantidade escalar a um vetor "\(u\)" (complexo ou real). Este quociente pode tanto ser empregado para uma dada matriz complexa Hermitiana (uma matriz quadrada de entradas complexas que é igual a sua própria matriz transposta conjugada) como para uma matriz real (esta, neste caso, deverá ser simétrica). De todas as formas, é válido apontar que ambas as matrizes são diagonalizáveis apenas com valores reais.

O método iterativo do quociente de Rayleigh entrega uma sequência de soluções aproximadas que, no limite, convergem para a solução verdadeira. Este algoritmo de iteração converge cubicamente para a matriz sob análise (Hermitiana ou real simétrica) desde que o vetor inicial usado para fazer tal iteração seja suficientemente próximo ao eigenvector de dita matriz. Primeiramente escolhemos algum valor inicial para "\(mu\)" ("chute" do eigenvalue inicial da matriz hermitiana) e algum valor inicial para o eigenvector "\(v\)". Nosso objetivo é chegar à convergência de ambos os valores "\(mu\)" e "\(v\)".

REFERÊNCIA: PARLETT, B. N. "The Rayleigh Quotient Iteration and Some Generalizations for Nonnormal Matrices". Mathematics of Computation, v. 28, n. 127, p. 679–693, jul. 1974. Disponível em: < https://www.ams.org/journals/mcom/1974-28-127/S0025-5718-1974-0405823-3/S0025-5718-1974-0405823-3.pdf >

A implementação em Python:

# Importar pylab é o mesmo que importar "numpy" e "pyplot" juntos  
from pylab import *

# A matriz do exercício anterior (4.3) é: 
A = array([[6,2,1],[2,3,1],[1,1,1]])

# Com este comando apenas mostro a matriz A que defini acima
print ("\nA = ")
print (A)

# Vou agora encontrar os eigenvectors da matriz A e mostrá-los
lam,X = eig(A)
eigenvalues = lam
print ("\nEigenvalues: ",lam)

# A questão pede para eu escolher um vetor aleatório como vetor inicial do algoritmo de 
iteração. Tomemos 
# alguns aleatórios para ver o que ocorre:
vetor_inicial = randn(3)
v = vetor_inicial
print ("\nVetor inicial v: ",v)
print (" ")

#"Chute" para valor inicial do eigenvalue da matriz Hermitiana 
mu = 1. + 1j
print ("k = %s, mu = %21.16f + %21.16fj" %  (0,mu.real,mu.imag))

# A iteração de Rayleigh
# A função "eye(3)" retorna a matriz identidade de dimensão 3
# A função solve vai resolver o sistema linear Bw=v
for k in range(1,10):
       B = A - mu*eye(3)
       try:
             w = solve(B,v)
       except:
             print ("A matriz é singular e a solução convergiu")
             break

# Atualizamos o valor de mu    
      v = w / norm(w,2)
      mu = dot(conj(v.T), dot(A,v))
      print ("k = %s, mu = %21.16f + %21.16fj" %  (k,mu.real,mu.imag))

Ao definirmos o \(vetor\_inicial=randn(3)\) acima, não há convergência; porém, ao colocarmos \(vetor\_inicial=eigenvalues\), a convergência é imediata.