Resolución de matrices muy grandes en "trozos"

Question

Digamos que tienes una matriz muy densa de 30000x30000 elementos. La matriz muy densa viene de la ecuación de radiosidad, que he discutido aquí .

Digamos que tienes Ax = B. Tienes B, y A tiene 30000x30000 elementos. Resolver: Tienes que encontrar x dado B.

Resolver una matriz como ésta en código C introduce un problema interesante: no se puede almacenar todo en la memoria a la vez.

Así que quiero saber cómo resolver este problema "a trozos": quiero descomponer matemáticamente la matriz de 30000x30000 en submatrices más pequeñas (¿tal vez en trozos de 30000x20?), y resolverlo de alguna manera.

¿Qué algoritmos existen para descomponer/resolver una matriz "en trozos" o pasos, que puedan resolverse independientemente de las restricciones de memoria? Puede ser una técnica iterativa.

Answer 1

Esto es demasiado largo para un comentario, así que voy a lanzar esta pieza de elección por Forman Acton en su Métodos numéricos que funcionan sobre el tema del tratamiento de las matrices grandes:

Cada vez que una persona pregunta con entusiasmo si mi ordenador puede resolver un conjunto de 300 ecuaciones en 300 incógnitas, debo reprimir rápidamente la tentación de responder: "Sí, pero ¿por qué molestarse?". Hay, en efecto, conjuntos legítimos de ecuaciones tan grandes. Surgen al sustituir una ecuación diferencial parcial de una ecuación diferencial parcial en un conjunto de puntos de la cuadrícula, y la persona que sabe lo suficiente para abordar este tipo de problema también suele saber qué tipo de ordenador necesita. Lo más probable es que nuestro amigo preguntón sufra un problema muy diferente: probablemente ha recogido un conjunto de datos experimentales datos experimentales y ahora está intentando ajustar un modelo de 300 parámetros a los mismos, ¡por mínimos cuadrados! ¡cuadrados! Cuanto antes se pueda sacar a este tipo de su oficina, antes podrá podrás volver al trabajo útil, pero estos tipos son persistentes. Han montado modelos de tres parámetros en máquinas de escritorio sin dificultades serias y ahora el ordenador electrónico les permite una visión más visiones grandiosas. Saltan de lo conocido a lo desconocido con una aterradora inocencia y la perenne confianza en que cada parámetro está totalmente justificado. No sirve de nada señalar que varios parámetros que es casi seguro que compiten para "explicar" las mismas variaciones en los datos y, por tanto, el sistema de ecuaciones será casi indeterminado. Sí sirve No sirve de nada señalar que todas las grandes matrices de mínimos cuadrados se esfuerzan por ser subconjuntos adecuados de la matriz de Hilbert, que es prácticamente indeterminada y no invertible, por lo que incluso si los 300 parámetros fueran independientes, las ecuaciones de ajuste seguirían siendo violentamente inestables. Todo esto, repito, no sirve de nada, y se termina por enfadado y echando al tipo de tu oficina...

...El director del centro de cómputo debe evitar el saqueo de valiosos tiempo de computación por parte de estos aspirantes a ajustadores de muchos parámetros. La tarea no es agradable, pero los usuarios legítimos de los ordenadores tienen derechos también. La alternativa compromete a todo el mundo a dos miserables semanas de chapotear en grandes cantidades de "Resultados" que son manifiestamente imposibles, sin una forma visible de encontrar el problema. El problema, por supuesto Por supuesto, surge de la búsqueda de una buena respuesta a un problema mal planteado, pero un director de informática rara vez sabe lo suficiente sobre el tema para ganar de esas discusiones con el proponente del problema, y el callejón sin salida finalmente tiene que romper con la violencia, que por lo tanto bien podría que, por lo tanto, podría usarse desde el principio.

Answer 2

La primera pregunta que debe hacerse a cualquier persona que le presente una matriz de gran tamaño: "¿la matriz es densa o dispersa?" En este último caso, no es necesario asignar espacio a las entradas nulas, por lo que hay que buscar técnicas especiales para almacenarlas (que, según tengo entendido, se basan en una buena dosis de teoría de grafos en el caso general, aunque las matrices de bandas se siguen manejando almacenando sus bandas en un formato adecuado).

Ahora bien, si incluso después de eso tienes la excusa perfecta para tener una matriz densa grande (que sigo creyendo que es bastante improbable), hay una forma de invertir y tomar el determinante de una matriz grande mediante la partición.

Digamos que tenemos

$\textbf{A}=\begin{pmatrix}\textbf{E}&\textbf{F}\\ \textbf{G}&\textbf{H}\end{pmatrix}$

donde $\textbf{E}$ y $\textbf{H}$ son matrices cuadradas de dimensiones $m\times m$ y $n\times n$ respectivamente, y $\textbf{F}$ y $\textbf{G}$ están dimensionados adecuadamente (por lo que la dimensión de $\textbf{A}$ es $(m+n)\times(m+n)$ ). La inversa puede calcularse entonces como

$\textbf{A}^{-1}=\begin{pmatrix}\textbf{E}^{-1}+\left(\textbf{E}^{-1}\textbf{F}\right)(\textbf{H}-\textbf{G}\textbf{E}^{-1}\textbf{F})^{-1}\left(\textbf{G}\textbf{E}^{-1}\right)&-\left(\textbf{E}^{-1}\textbf{F}\right)(\textbf{H}-\textbf{G}\textbf{E}^{-1}\textbf{F})^{-1}\\-(\textbf{H}-\textbf{G}\textbf{E}^{-1}\textbf{F})^{-1}\left(\textbf{G}\textbf{E}^{-1}\right)&(\textbf{H}-\textbf{G}\textbf{E}^{-1}\textbf{F})^{-1}\end{pmatrix}$

donde los paréntesis enfatizan los factores comunes que podría aprovechar para calcular una sola vez.

En cuanto al determinante, viene dado por

$\det\;\textbf{A}=\det\left(\textbf{H}-\textbf{G}\textbf{E}^{-1}\textbf{F}\right)\;\det\;\textbf{E}$

EDITAR:

No estoy del todo seguro de lo que parecía insatisfactorio en esta respuesta, así que aquí hay un poco más de exposición: como siempre, la "inversa" aquí es simplemente notación. Lo más seguro es que primero se realice la descomposición LU en $\textbf{E}$ y $\textbf{H}$ primero. También se dividiría el lado derecho $\textbf{b}$ en consecuencia:

$\textbf{b}=\begin{pmatrix}\textbf{c}\\ \textbf{d}\end{pmatrix}$

así que $\textbf{c}$ es un vector de longitud m y $\textbf{d}$ es un vector de longitud n.

Multiplicando formalmente la inversa partitiva con el lado derecho partitivo se obtiene

$\begin{pmatrix}\textbf{E}^{-1}+\left(\textbf{E}^{-1}\textbf{F}\right)(\textbf{H}-\textbf{G}\textbf{E}^{-1}\textbf{F})^{-1}\left(\textbf{G}\textbf{E}^{-1}\right)&-\left(\textbf{E}^{-1}\textbf{F}\right)(\textbf{H}-\textbf{G}\textbf{E}^{-1}\textbf{F})^{-1}\\-(\textbf{H}-\textbf{G}\textbf{E}^{-1}\textbf{F})^{-1}\left(\textbf{G}\textbf{E}^{-1}\right)&(\textbf{H}-\textbf{G}\textbf{E}^{-1}\textbf{F})^{-1}\end{pmatrix}\cdot\begin{pmatrix}\textbf{c}\\ \textbf{d}\end{pmatrix}$

que ampliado y simplificado es

$\begin{pmatrix}\textbf{E}^{-1}\textbf{c}+\left(\textbf{E}^{-1}\textbf{F}\right)(\textbf{H}-\textbf{G}\textbf{E}^{-1}\textbf{F})^{-1}\left(\textbf{G}\textbf{E}^{-1}\textbf{c}-\textbf{d}\right)\\-(\textbf{H}-\textbf{G}\textbf{E}^{-1}\textbf{F})^{-1}\left(\textbf{G}\textbf{E}^{-1}\textbf{c}-\textbf{d}\right)\end{pmatrix}$

Llegados a este punto, deberías ser capaz de averiguar cómo utilizarías una descomposición disponible de $\textbf{E}$ o $\textbf{H}$ y qué subexpresiones comunes se pueden calcular una sola vez.

Answer 3

Si lo único que te interesa es resolver $A x = B$ En este caso, sugeriría que se buscaran métodos iterativos y/o de bloque. Teniendo en cuenta su otra pregunta sobre la positividad de las soluciones, le sugiero que busque Aplicada y Computacional Álgebra lineal: Un primer curso por Charlie Byrne. El algoritmo MART descrito en el mismo se utiliza para encontrar $x$ soluciones a $A x = B$ m, y aborda muchos problemas similares de la tomografía y el procesamiento de señales. Si te gusta, consigue su 'Applied Iterative Methods'.

Answer 4

Buscar "Out of core LU" en google podría ayudarte.

Esencialmente, hay dos opciones simples:

1. Un método iterativo que sólo depende de calcular repetidamente A*v y quizás A'*v para un vector "v" dado. Dicho cálculo puede realizarse sin tener la matriz A completa en la memoria RAM

2. Un método directo, para esto el "out of core LU" podría ayudar.

En general, para las dos opciones anteriores, la búsqueda de algoritmos de memoria externa será útil.