Valores propios de un lápiz de matriz lineal (Ax = λ Bx) mediante el método de iteración del subespacio

Question

Estoy tratando de encontrar los valores propios más a la derecha de la problema generalizado de valores propios ( $Ax = \lambda Bx$ ) utilizando el método de iteración del subespacio. Esta formulación surge del análisis de estabilidad del flujo donde la matriz ' B ' es singular y tanto A como B son escaso matrices del orden de 50,000 . Básicamente estoy tratando de comprobar si la matriz lápiz tiene valores propios con un positivo parte real (lo que significa que el flujo es inestable a pequeñas perturbaciones), y cuando lo hace, necesito saber el parte imaginaria de estos valores propios ya que da la frecuencia de las oscilaciones.

La idea es utilizar iteraciones de potencia (con ortonormalización) para aproximar un subconjunto de los vectores propios ( T ). Utilizando estos vectores como base para un subespacio invariante las matrices del subespacio reducido (A' = T T A T y B' = T T B T). A continuación, este problema de valores propios de dimensión reducida se resuelve directamente utilizando el algoritmo QZ.

El algoritmo es así -
Paso 1: Generar un conjunto aleatorio de $m$ vectores ortonormales ( $R^{(0)}$ ), donde m es mayor que $p$ el número de valores propios buscados y mucho menos que $n$ el orden de las matrices
Paso 2: Calcular iterativamente $R^{(i+1)} = (A^{-1} B) R^{(i)}$ resolviendo $AR^{(i+1)} = BR^{(i)}$ y normalizar $R$ después de cada iteración
Paso 3: Cuando $\kappa(R)$ aumenta por encima de un umbral, realice una descomposición QR con $T = R^{new} = Q$
Paso 4: Calcular las matrices del subespacio $A' = T^TA T$ y $B' = T^TB T$
Paso 5: Encontrar los valores y vectores propios de $A'y = \lambda B'y$ utilizando el algoritmo QZ
Paso 6: Recalcular $R$ utilizando $R = Ty$ y volver al paso 2.

Este proceso devolverá los valores propios más pequeños de la matriz lápiz. Para obtener los valores propios más grandes, precondiciono la matriz utilizando transformación cambio-inversión .

Todo el algoritmo parece ser perfecto en teoría, pero una reflexión sobre este proceso me hizo pensar que el subespacio abarcado por los vectores propios de la matriz lápiz no será realmente invariante en A y B de forma independiente. Considero que ésta es la razón por la que los valores propios no convergen cuando pruebo un problema de ejemplo utilizando mi código. Sin embargo, cuando elijo B para ser el identidad matriz, obtengo los valores propios correctos.

¿Es correcto mi razonamiento, o hay alguna otra fuente de error en este algoritmo?
Al darme cuenta de este problema, he modificado el algoritmo con $C = A^{-1}B$ calculado de antemano y sustituido en todos los lugares (Hace $C$ singular pero sí proporciona los valores propios). Sin embargo, todavía tengo curiosidad por saber si he identificado el problema correctamente, y por qué en absoluto funciona el algoritmo (se menciona en la literatura), si es defectuoso.

Gracias de antemano.

Nota : Soy estudiante de ingeniería y he estudiado sobre subespacios invariantes sólo para el proyecto en el que estoy trabajando.

Answer 1

i) $\det(A-\lambda B)=0$ se puede reescribir $\det(\dfrac{1}{\lambda}I-A^{-1}B)=0$ . Se calcula el mayor valor propio de $A^{-1}B$ que es el menor valor propio generalizado del lápiz $A-\lambda B$ .

ii) Usted calcula $A^{-1}$ (¿cuál es el tiempo de cálculo?).

Cuando una matriz es grande, suele estar mal condicionada; lo que es $CN(A)$ el número de condición de $A$ ? Si $CN(A)\approx 10^k$ entonces debe considerar $k$ dígitos significativos suplementarios a lo largo de los cálculos. A continuación, calcule $CN(A)$ (existen métodos de aproximación de $CN$ que tienen complejidad $O(n^2)$ ) antes de calcular $A^{-1}$ .

iii) Al calcular $A^{-1}Bx$ no calcule $(A^{-1}B)x$ , pero calcula $A^{-1}(Bx)$ , o mejor, resolver en $y$ la ecuación $Ay=Bx$ (en este caso, no es necesario $A^{-1}$ ).

EDITAR. Respuesta al OP.

i) Tal vez suponga que $A,B$ son simétricos reales (??); eso no implica que los valores propios de $A^{-1}B$ ¡son reales! -excepto si $A>0$ -

ii) Si he entendido bien, se obtiene (después de algunas iteraciones) una base ortonormal de un subespacio límite $E$ que es invariable para $A^{-1}B$ . En primer lugar, $E$ no es invariable por $A,B$ . En segundo lugar, ¿cómo se elige la dimensión de $E$ ? En tercer lugar, no entiendo por qué se obtiene ese límite (especialmente si hay valores propios muy cercanos); pero admitámoslo.

Debe escribir la matriz de $A^{-1}B$ en $E$ -utilizando su base ortonormal- y diagonalizar esta matriz. En particular, debes calcular explícitamente $A^{-1}B$ y no trabajas sobre $A,B$ .

iii) Normalmente, sólo se debe trabajar en $\mathbb{R}$ -las partes imaginarias son parásitas; elimínalas-

Sobre los dígitos -si se trabaja con los números reales-

complex*16 funciona con $16.8=128$ bits, es decir, con $38$ dígitos significantes. Debido a $CN(A)$ , sigue siendo $38-19=19$ dígitos del significante, eso es suficiente -excepto si se usan números complejos-