¿Cuál es la forma más rápida de calcular el valor propio más grande de una matriz general?

27

EDITAR: Estoy probando si algún valor propio tiene una magnitud de uno o más.

Necesito encontrar el valor propio absoluto más grande de una gran matriz dispersa, no simétrica.

He estado usando la eigen()función de R , que usa el algoritmo QR de EISPACK o LAPACK para encontrar todos los valores propios y luego uso abs()para obtener los valores absolutos. Sin embargo, necesito hacerlo más rápido.

También he intentado usar la interfaz ARPACK en el igraphpaquete R. Sin embargo, dio un error para una de mis matrices.

La implementación final debe ser accesible desde R.

Probablemente habrá múltiples valores propios de la misma magnitud.

¿Tienes alguna sugerencia?

EDITAR: Precisión solo necesita ser 1e-11. Una matriz "típica" ha sido hasta ahora . He podido hacer una factorización QR en él. Sin embargo, también es posible tener unos mucho más grandes. Actualmente estoy empezando a leer sobre el algoritmo Arnoldi. Entiendo que está relacionado con Lanczsos.386×386

EDIT2: si tengo varias matrices que estoy "probando" y sé que hay una gran submatriz que no varía. ¿Es posible ignorarlo / descartarlo?

poder
fuente
Vea mi respuesta aquí: scicomp.stackexchange.com/a/1679/979 . Este es un tema de investigación actual y los métodos actuales pueden funcionar mejor que Lanczos. El problema de calcular valores singulares es equivalente al problema de calcular valores propios.
dranxo
2
Matriz 400x400! = Grande. Además, ¿qué significa mayor si "Probablemente habrá múltiples valores propios de la misma magnitud"? En tierra numpy: linalg.eig (random.normal (size = (400,400))) toma aproximadamente medio segundo. ¿Esto es demasiado lento?
meawoppl
@meawoppl sí, medio segundo es demasiado lento. Esto se debe a que es parte de otro algoritmo que ejecuta este cálculo muchas veces.
poder
1
@power gotcah. ¿Tiene una aproximación al vector propio. es decir, ¿es probable que sea similar a la última solución, o puede hacer una conjetura acerca de su estructura?
meawoppl

Respuestas:

14

Depende mucho del tamaño de su matriz, en el caso a gran escala también de si es escasa y de la precisión que desea lograr.

Si su matriz es demasiado grande para permitir una factorización única y necesita una alta precisión, el algoritmo de Lanczsos es probablemente la forma más rápida. En el caso no simétrico, se necesita el algoritmo Arnoldi, que es numéricamente inestable, por lo que una implementación debe abordar esto (es algo incómodo de curar).

Si este no es el caso en su problema, brinde información más específica en su pregunta. Luego agregue un comentario a esta respuesta, y lo actualizaré.

Editar: [Esto era para la versión anterior de la pregunta, buscando el valor propio más grande.] Como su matriz es pequeña y aparentemente densa, haría la iteración de Arnoldi en B = (IA) ^ {- 1}, usando una inicial la factorización triangular permutada de IA tiene una multiplicación barata por B. (O calcule un inverso explícito, pero esto cuesta 3 veces más que la factorización). Desea probar si B tiene un valor propio negativo. Al trabajar con B en lugar de A, los valores propios negativos están mucho mejor separados, por lo que si hay uno, debe converger rápidamente.

Pero tengo curiosidad acerca de dónde viene su problema. Las matrices no simétricas generalmente tienen valores propios complejos, por lo que "" más grande "" ni siquiera está bien definido. Por lo tanto, debe saber más sobre su problema, lo que podría ayudarlo a sugerir cómo resolverlo aún más rápido y / o de manera más confiable.

Edit2: Es difícil obtener con Arnoldi un subconjunto particular de interés. Para obtener los valores propios absolutamente más grandes de manera confiable, debe hacer una iteración del subespacio utilizando la matriz original, con un tamaño del subespacio que coincida o supere el número de valores propios que se espera que sea de 1 o mayor en magnitud. En matrices pequeñas, esto será más lento que el algoritmo QR, pero en matrices grandes será mucho más rápido.

Arnold Neumaier
fuente
Necesito probar si el valor propio más grande es mayor que 1. La precisión solo debe ser de 1e-11. Una matriz "típica" ha sido hasta ahora 386 x 386. He podido hacer una factorización QR en ella. Sin embargo, también es posible tener unos mucho más grandes. Actualmente estoy empezando a leer sobre el algoritmo Arnoldi. Entiendo que está relacionado con Lanczsos.
poder
Esta información pertenece a su pregunta, así que edítela y también agregue más información (¿por qué los valores propios son reales o qué significa mayor?), Vea la edición de mi respuesta.
Arnold Neumaier
Lamento no haberme explicado claramente. Tampoco expliqué claramente que los valores propios son complejos. Estoy probando si algún valor propio tiene una magnitud de uno o más.
poder
1
Esto tiene más sentido, pero ahora mi receta con funciona bien solo si el valor propio deficiente es real> 1. Por otro lado, la nueva información probablemente implica que tienes pocas opciones, pero calcular todos los valores propios. - Actualice su pregunta para transmitir la información adicional. (IA)1
Arnold Neumaier el
1
ver edición 2 en mi respuesta
Arnold Neumaier
7

La iteración de potencia (o método de potencia), por ejemplo, lo que describe Dan, siempre debe converger, aunque a la velocidad .El |λnorte-1/ /λnorteEl |

Si está cerca de λ n , será lento, pero puede usar la extrapolación para evitarlo . Puede parecer complicado, pero en el documento se da una implementación en pseudocódigo.λnorte-1λnorte

Pedro
fuente
1
¿Qué pasa si | λ (n − 1) | = | λ (n) | ?
poder
@power, entonces la iteración de potencia normal no convergerá. No sé qué tan bien los métodos de extrapolación distinguirán entre los diferentes valores propios, tendrá que leer el documento para eso.
Pedro
2
@power: Todas las cosas reconsideradas, si , entonces la iteración de potencia seguirá convergiendo al valor propio correcto. El vector propio resultante, en el que no parece estar interesado de ninguna manera, será una combinación lineal de los vectores propios correspondientes a λ n y λ n - 1 . El |λnorte-1El |=El |λnorteEl |λnorteλnorte-1
Pedro
¿Tiene alguna referencia a un trabajo o libro académico que lo respalde? Además, ¿qué pasa si \ lambda_ {n} es complejo?
poder
55
Si hay varios valores propios diferentes de módulo máximo, la iteración de potencia converge solo en circunstancias excepcionales. Generalmente oscila de una manera algo impredecible.
Arnold Neumaier
5

Ha habido una buena investigación sobre esto recientemente. Los nuevos enfoques utilizan "algoritmos aleatorios" que solo requieren unas pocas lecturas de su matriz para obtener una buena precisión en los valores propios más grandes. Esto contrasta con las iteraciones de potencia que requieren varias multiplicaciones de matriz-vector para alcanzar una alta precisión.

Puede leer más sobre la nueva investigación aquí:

http://math.berkeley.edu/~strain/273.F10/martinsson.tygert.rokhlin.randomized.decomposition.pdf

http://arxiv.org/abs/0909.4061

Este código lo hará por usted:

http://cims.nyu.edu/~tygert/software.html

https://bitbucket.org/rcompton/pca_hgdp/raw/be45a1d9a7077b60219f7017af0130c7f43d7b52/pca.m

http://code.google.com/p/redsvd/

https://cwiki.apache.org/MAHOUT/stochastic-singular-value-decomposition.html

Si su idioma de elección no está allí, puede rodar su propia SVD aleatoria con bastante facilidad; solo requiere una multiplicación de vector de matriz seguida de una llamada a un SVD estándar.

dranxo
fuente
4

Aquí encontrará una introducción algorítmica al algoritmo Jacobi-Davidson, que calcula el valor propio máximo.

En este artículo se exploran los aspectos matemáticos. JD permite matrices generales (reales o complejas) y puede usarse para calcular rangos de valores propios.

Aquí puede encontrar varias implementaciones de biblioteca JDQR y JDQZ (incluida una interfaz C, a la que debería poder vincular desde R).

Respiro de muerte
fuente
No he podido encontrar ninguna literatura que explícitamente indique que el método Jacobi-Davidson funciona para una matriz real y general.
poder
A menos que cada artículo establezca explícitamente una restricción y el argumento de convergencia se base en la restricción que no importa.
Deathbreath
Aquí hay otra explicación de JD. Las matrices consideradas son completamente generales. No se explota ninguna estructura especial y los resultados específicos de las matrices hermitianas se comparan y contrastan, por ejemplo, la convergencia para las matrices generales es cuadrática, pero cúbica para las matrices hermitianas.
Deathbreath
gracias por esto. No encuentro ningún código C para una matriz general, así que tendré que escribir el mío. Los enlaces a los algoritmos parecen ser solo para matrices hermetianas.
poder
1
@power tampoco encontrará en la literatura un resultado que establezca que las implementaciones estándar de QR convergen para una matriz general real, eso es un problema abierto y, de hecho, no hace mucho tiempo se encontró un contraejemplo para el código QR en LAPACK.
Federico Poloni
2

En tu publicación original, dices:

"También intenté usar la interfaz ARPACK en el paquete igraph R. Sin embargo, me dio un error en una de mis matrices".

Me interesaría saber más sobre el error. Si puede hacer que esta matriz esté disponible públicamente en algún lugar, me interesaría probar ARPACK en ella.

Según lo que he leído anteriormente, esperaría que ARPACK hiciera un muy buen trabajo extrayendo los valores propios más grandes (o algunos de los más grandes) de una matriz dispersa. Para ser más específico, esperaría que los métodos de Arnoldi funcionen bien para este caso y, por supuesto, en eso se basa ARPACK.

La convergencia lenta del método de potencia cuando hay valores propios estrechamente espaciados en la región de interés se mencionó anteriormente. Arnoldi mejora esto al iterar con varios vectores en lugar del método de potencia.

Bill Greene
fuente
Veré si puedo encontrar mi trabajo desde entonces. Trabajé en esto hace un año.
poder
0

No es el mas rapido forma rápida, pero una forma razonablemente rápida es simplemente golpear un vector (inicialmente aleatorio) con la matriz repetidamente y luego normalizar cada pocos pasos. Eventualmente convergerá al vector propio más grande, y la ganancia en la norma para un solo paso es el valor propio asociado.

Esto funciona mejor cuando el valor propio más grande es sustancialmente más grande que cualquier otro valor propio. Si otro valor propio es cercano en magnitud al mayor, esto tardará un tiempo en converger, y puede ser difícil determinar si ha convergido.

Dan
fuente
1
Gracias Dan, sin embargo: en mis matrices, algunos de los otros valores propios tendrán una magnitud similar (si no la misma) que la más grande. ¿Su método es similar a la iteración de potencia y la iteración de cociente de Rayleigh? Batterson y Smillie (1990) escriben que para algunas matrices no simétricas, la iteración del cociente de Rayleigh no convergerá. Batterson, S., Smillie, J (1990) "Iteración del cociente de Rayleigh para matrices no simétricas", Matemáticas de la computación, vol 55, num 191, P 169-178
potencia
Si otros valores propios tienen la misma magnitud que el más grande ... ¿entonces esos valores no son también "el más grande" también?
ely
@EMS: seguirían siendo los "valores propios más grandes", pero la presencia de más de uno más grande aún eliminaría la convergencia.
Dan
Me pregunto a qué valor propio desea que converja. Cosas como el método de cociente / potencia de Rayleigh se entiende cuando existe un valor propio distintivo más grande. Su pregunta pide encontrar el valor propio más grande, pero luego parece que en realidad no está bien definido para su problema. Simplemente estoy confundido por el título de la publicación.
ely
-1

El paquete R RARPACK funciona para mí. Y parece ser muy rápido, ya que es solo una interfaz para ARPACK, el paquete estándar para álgebra lineal dispersa (lo que significa calcular algunos valores propios y vectores propios).

HoangDT
fuente
Bienvenido a SciComp! Como dice la pregunta, ARPACK no funciona para el OP, por lo que esta respuesta no es realmente útil.
Christian Clason
@HoangDT Esta pregunta es anterior al RARPACK
poder el