Juegos previos

La solución numérica de los sistemas lineales de ecuaciones es un paso indispensable en muchas áreas de la matemática aplicada, la ingeniería y la industria de TI, ya sea trabajando con gráficos, calculando la aerodinámica de un avión u optimizando la logística. La "máquina" ahora de moda tampoco habría progresado mucho sin ella. Además, la solución de los sistemas lineales, por regla general, consume el mayor porcentaje de todos los costos informáticos. Está claro que este componente crítico requiere la optimización de la velocidad máxima.

A menudo trabajan con los llamados matrices dispersas: aquellas cuyos ceros son órdenes de magnitud mayores que otros valores. Esto, por ejemplo, es inevitable si se trata de ecuaciones diferenciales parciales u otros procesos en los que los elementos que surgen en las relaciones lineales definitorias están conectados solo con "vecinos". Aquí hay un posible ejemplo de una matriz dispersa para la ecuación de Poisson unidimensional conocida en física clásica $$$- \ phi ^ {''} = f$ en una cuadrícula uniforme (sí, aunque no haya tantos ceros en ella, ¡pero al rectificar la cuadrícula estarán sanos!):

$$

$$\ begin {pmatrix} 2 y -1 y 0 y 0 y 0 \\ -1 y 2 y -1 y 0 y 0 \\ 0 y -1 y 2 y -1 y 0 \\ 0 y 0 y -1 & 2 & -1 \\ 0 & 0 & 0 & -1 & 2 \ end {pmatrix}$$

Las matrices opuestas a ellas, aquellas en las que no prestan atención al número de ceros y tienen en cuenta todos los componentes sin excepción, se denominan densas.

Las matrices dispersas a menudo, por diversas razones, se presentan en un formato de columna comprimido: Columna dispersa comprimida (CSC). Esta descripción utiliza dos matrices enteras y un punto flotante. Deja que la matriz tenga todo $$$nnz_A$ distinto de cero y $$$N$ columnas Los elementos de la matriz se enumeran en columnas de izquierda a derecha, sin excepciones. Primer conjunto $$$i_A$ longitudes $$$nnz_A$ contiene números de fila de componentes matriciales distintos de cero. Segunda matriz $$$j_A$ longitudes $$$N + 1$ contiene ... no, no los números de columna, porque entonces la matriz se escribiría en formato de coordenadas (coordenadas) o ternario (triplete). Y la segunda matriz contiene los números de serie de los componentes de la matriz con los que comienzan las columnas, incluida una columna ficticia adicional al final. Finalmente, la tercera matriz. $$$v_A$ longitudes $$$nnz_A$ ya contiene los propios componentes, en el orden correspondiente a la primera matriz. Aquí, por ejemplo, bajo el supuesto de que la numeración de filas y columnas de matrices comienza desde cero, para una matriz particular

$$

$$A = \ begin {pmatrix} 0 y 1 y 0 y 4 y -1 \\ 7 y 0 y 2.1 y 0 y 3 \\ 0 y 0 y 0 y 10 y 0 \ end {pmatrix}$$

las matrices serán $$$i_A = \ {1, 0, 1, 0, 2, 0, 1 \}$ , $$$j_A = \ {0, 1, 2, 3, 5, 7 \}$ , $$$v_A = \ {7, 1, 2.1, 4, 10, -1, 3 \}$ .

Los métodos para resolver sistemas lineales se dividen en dos grandes clases: directa e iterativa. Las líneas rectas se basan en la posibilidad de representar la matriz como un producto de dos matrices más simples, para luego dividir la solución en dos etapas simples. Iterative utiliza las propiedades únicas de los espacios lineales y trabaja en el antiguo como el método mundial de aproximación sucesiva de incógnitas a la solución deseada, y en el proceso de convergencia en sí, las matrices se usan, por regla general, solo para multiplicarlas por vectores. Los métodos iterativos son mucho más baratos que los directos, pero funcionan lentamente en matrices mal acondicionadas. Cuando la fiabilidad del hormigón armado es importante, utilice líneas rectas. Aquí estoy aquí y quiero tocar un poco.

Digamos para una matriz cuadrada $$$A$ la descomposición del formulario está disponible para nosotros $$$A = LU$ donde $$$L$ y $$$U$ - respectivamente, las matrices triangular inferior y triangular superior, respectivamente. El primero significa que tiene uno ceros sobre la diagonal, el segundo significa que está debajo de la diagonal. Cómo obtuvimos exactamente esta descomposición: no nos interesa ahora. Aquí hay un ejemplo simple de descomposición:

$$

$$\ begin {pmatrix} 1 & -1 & -1 \\ 2 & - 1 & -0.5 \\ 4 & -2 & -1.5 \ end {pmatrix} = \ begin {pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & 2 & 1 \ end {pmatrix} \ begin {pmatrix} 1 & -1 & -1 \\ 0 & 1 & 1.5 \\ 0 & 0 & -0.5 \ end {pmatrix}$$

¿Cómo en este caso resolver el sistema de ecuaciones? $$$Ax = f$ por ejemplo, con el lado derecho $$$f = \ begin {pmatrix} 4 \\ 2 \\ 3 \ end {pmatrix}$ ? La primera etapa es un movimiento directo (resolución hacia adelante = sustitución hacia adelante). Denotamos $$$y: = Ux$ y trabajar con el sistema $$$Ly = f$ . Porque $$$L$ - triangular inferior, sucesivamente en el ciclo encontramos todos los componentes $$$y$ de arriba a abajo:

$$

$$\ begin {pmatrix} 1 y 0 y 0 \\ 2 y 1 y 0 \\ 4 y 2 y 1 \ end {pmatrix} \ begin {pmatrix} y_1 \\ y_2 \\ y_3 \ end {pmatrix} = \ begin {pmatrix} 4 \\ 2 \\ 3 \ end {pmatrix} \ implica y = \ begin {pmatrix} 4 \\ -6 \\ -1 \ end {pmatrix}$$

La idea central es que, al encontrar $$$i$ componentes vectoriales $$$y$ se multiplica por una columna con el mismo número de matriz $$$L$ que luego se resta del lado derecho. La matriz misma parece colapsar de izquierda a derecha, disminuyendo de tamaño a medida que se encuentran más y más componentes vectoriales $$$y$ . Este proceso se llama eliminación de columna.

La segunda etapa es un movimiento hacia atrás (resolución hacia atrás = sustitución hacia atrás). Tener un vector encontrado $$$y$ nosotros decidimos $$$Ux = y$ . Aquí ya pasamos por los componentes de abajo hacia arriba, pero la idea sigue siendo la misma: $$$i$ La columna th se multiplica por el componente recién encontrado $$$x_i$ y se mueve hacia la derecha, y la matriz colapsa de derecha a izquierda. Todo el proceso se ilustra para la matriz mencionada en el ejemplo en la imagen a continuación:

$$

$$\ small \ begin {pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & 2 & 1 \ end {pmatrix} \ begin {pmatrix} y_1 \\ y_2 \\ y_3 \ end {pmatrix} = \ begin {pmatrix} 4 \\ 2 \\ 3 \ end {pmatrix} \ implica \ begin {pmatrix} 1 & 0 \\ 2 & 1 \ end {pmatrix} \ begin {pmatrix} y_2 \\ y_3 \ end {pmatrix } = \ begin {pmatrix} -6 \\ -13 \ end {pmatrix} \ implica \ begin {pmatrix} 1 \ end {pmatrix} \ begin {pmatrix} y_3 \ end {pmatrix} = \ begin {pmatrix} -1 \ end {pmatrix}$$

$$

$$\ small \ begin {pmatrix} 1 & -1 & -1 \\ 0 & 1 & 1.5 \\ 0 & 0 & -0.5 \ end {pmatrix} \ begin {pmatrix} x_1 \\ x_2 \\ x_3 \ end { pmatrix} = \ begin {pmatrix} 4 \\ -6 \\ -1 \ end {pmatrix} \ Rightarrow \ begin {pmatrix} 1 & -1 \\ 0 & 1 \ end {pmatrix} \ begin {pmatrix} x_1 \ \ x_2 \ end {pmatrix} = \ begin {pmatrix} 6 \\ -9 \ end {pmatrix} \ implica \ begin {pmatrix} 1 \ end {pmatrix} \ begin {pmatrix} x_1 \ end {pmatrix} = \ begin {pmatrix} -3 \ end {pmatrix}$$

y nuestra decisión será $$$x = \ begin {pmatrix} -3 \\ - 9 \\ 2 \ end {pmatrix}$ .

Si la matriz es densa, es decir, está completamente representada en forma de algún tipo de matriz, unidimensional o bidimensional, y el acceso a un elemento específico tiene lugar con el tiempo $$$O (1)$ , entonces un procedimiento de solución similar con la descomposición existente no es difícil y es fácil de codificar, por lo que ni siquiera perderemos tiempo en ello. ¿Qué pasa si la matriz es escasa? Aquí, en principio, tampoco es difícil. Aquí está el código C ++ para el movimiento hacia adelante en el que la solución $$$x$ Está escrito para el lugar en el lado derecho, sin verificar la exactitud de los datos de entrada (las matrices CSC corresponden a la matriz triangular inferior):

Algoritmo 1:

void forward_solve (const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<double>& vA, std::vector<double>& x) { size_t j, p, n = x.size(); for (j = 0; j < n; j++) //    { x[j] /= vA[jA[j]]; for (p = jA[j]+1; p < jA[j+1]; p++) x[iA[p]] -= vA[p] * x[j] ; } } 

Toda discusión adicional se referirá solo al golpe hacia adelante para resolver el sistema triangular inferior $$$Lx = f$ .

Empate

Pero qué pasa si el lado derecho, es decir vector a la derecha del signo igual $$$Lx = f$ ¿Tiene en sí una gran cantidad de ceros? Entonces tiene sentido omitir los cálculos asociados con las posiciones cero. Los cambios en el código en este caso son mínimos:

Algoritmo 2:

 void fake_sparse_forward_solve (const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<double>& vA, std::vector<double>& x) { size_t j, p, n = x.size(); for (j = 0; j < n; j++) //    { if(x[j]) { x[j] /= vA[jA [j]]; for (p = jA[j]+1; p < jA[j+1]; p++) x[iA[p]] -= vA[p] * x[j]; } } } 

Lo único que agregamos es la if , cuyo propósito es reducir el número de operaciones aritméticas a su número real. Si, por ejemplo, todo el lado derecho consiste en ceros, entonces no necesita contar nada: la solución será el lado derecho.

Todo se ve muy bien y, por supuesto, funcionará, pero aquí, después de un largo juego previo, el problema finalmente es visible: el rendimiento asintóticamente bajo de este solucionador en el caso de sistemas grandes. Y todo por el hecho de tener un bucle for . Cual es el problema Incluso si la condición dentro del if resulta ser extremadamente rara, no hay forma de alejarse del ciclo, y esto crea la complejidad del algoritmo $$$O (N)$ donde $$$N$ Es el tamaño de la matriz. No importa cómo los compiladores modernos optimicen los bucles, esta complejidad se hará sentir con grandes $$$N$ . Especialmente ofensivo cuando todo el vector $$$f$ consiste completamente en ceros, porque, como dijimos, ¡no haga nada en este caso! ¡Ni una sola operación aritmética! Que demonios $$$O (N)$ acción?

Bueno, digamos Aun así, ¿por qué no puede soportar la ejecución inactiva, porque los cálculos reales con números reales, es decir, los que quedan bajo if todavía serán pocos? Pero el hecho es que este procedimiento de flujo directo con un lado derecho enrarecido a menudo se usa en ciclos externos y subyace a la descomposición de Cholesky $$$A = LL ^ T$ y descomposición LU de aspecto izquierdo. Sí, sí, una de esas expansiones, sin la capacidad de hacer que todos estos movimientos directos e inversos para resolver sistemas lineales, pierdan todo significado práctico.

Teorema Si la matriz es simétrica positiva definida (SPD), entonces se puede representar como $$$A = LL ^ T$ la única forma donde $$$L$ - matriz triangular inferior con elementos positivos en diagonal.

Para matrices SPD altamente dispersas, se utiliza la descomposición de Cholesky hacia arriba. Representación esquemática de la descomposición en forma de bloque de matriz

$$

$$\ begin {pmatrix} \ mathbf {L} _ {11} & \ mathbf {0} \\ \ mathbf {l} _ {12} ^ T & l_ {22} \ end {pmatrix} \ begin {pmatrix} \ mathbf {L} _ {11} ^ T & \ mathbf {l} _ {12} \\ \ mathbf {0} & l_ {22} \ end {pmatrix} = \ begin {pmatrix} \ mathbf {A} _ { 11} & \ mathbf {a} _ {12} \\ \ mathbf {a} _ {12} ^ T & a_ {22} \ end {pmatrix},$$

Todo el proceso de factorización se puede dividir lógicamente en solo tres pasos.

Algoritmo 3:

1. método superior de descomposición de Cholesky de menor dimensión $$$\ mathbf {L} _ {11} \ mathbf {L} _ {11} ^ T = \ mathbf {A} _ {11}$ (recurrencia!)
2. sistema triangular inferior con escaso lado derecho $$$\ mathbf {L} _ {11} \ mathbf {l} _ {12} = \ mathbf {a} _ {12}$ (¡ Aquí está! )
3. calculo $$$l_ {22} = \ sqrt {a_ {22} - \ mathbf {l} _ {12} ^ T \ mathbf {l} _ {12}}$ (operación trivial)

En la práctica, esto se implementa de tal manera que los pasos 3 y 2 se realizan en un ciclo grande y en este orden. Por lo tanto, la matriz se ejecuta diagonalmente de arriba a abajo, aumentando la matriz $$$L$ línea por línea en cada iteración del bucle.

Si en tal agloritmo la complejidad hacia adelante en el paso 2 será $$$O (N)$ donde $$$N$ Es el tamaño de la matriz triangular inferior. $$$\ mathbf {L} _ {11}$ en una iteración arbitraria de un ciclo grande, entonces la complejidad de la descomposición completa será al menos $$$O (N ^ 2)$ ! ¡Oh, no me gustaría eso!

Estado del arte

Muchos algoritmos se basan de alguna manera en imitar las acciones humanas para resolver problemas. Si le das a una persona un sistema lineal triangular inferior con un lado derecho, en el que solo 1-2 no son cero, entonces, por supuesto, primero ejecuta el vector del lado derecho con los ojos de arriba a abajo (ese maldito ciclo de complejidad $$$O (N)$ ! ) para encontrar estos nonzeros. Luego los usará solo, sin perder tiempo en los componentes cero, porque este último no afectará la solución: no tiene sentido dividir ceros en los componentes diagonales de la matriz, así como mover la columna multiplicada por cero a la derecha. Este es el algoritmo anterior 2. No hay milagros. Pero, ¿qué pasa si a una persona se le da de inmediato el número de componentes distintos de cero de otras fuentes? Por ejemplo, si el lado derecho es una columna de alguna otra matriz, como es el caso de la descomposición de Cholesky, entonces tenemos acceso instantáneo a sus componentes distintos de cero, si se solicita secuencialmente:

 //     j,     : for (size_t p = jA[j]; p < jA[j+1]; p++) //    vA[p]     iA[p]   j 

La complejidad de este acceso es $$$O (nnz_j)$ donde $$$nnz_j$ Es el número de componentes distintos de cero en la columna j. ¡Gracias a Dios por el formato CSC! Como puede ver, no solo se usa para ahorrar memoria.

En este caso, podemos captar la esencia de lo que sucede al resolver por el método de curso directo $$$Lx = f$ para matrices dispersas $$$L$ y lado derecho $$$f$ . Aguante la respiración: tomamos un componente distinto de cero $$$f_i$ en el lado derecho, encontramos la variable correspondiente $$$x_i$ dividiendo por $$$L_ {ii}$ y luego, multiplicando toda la columna i-ésima por esta variable encontrada, ¡introducimos no ceros adicionales en el lado derecho restando esta columna del lado derecho! Este proceso está bien descrito en el lenguaje de ... gráficos. Por otra parte, orientado y no cíclico.

Para una matriz triangular inferior que no tiene ceros en la diagonal, definimos un gráfico conectado. Suponemos que la numeración de filas y columnas comienza desde cero.
Definición Un gráfico de conectividad para una matriz triangular inferior $$$L$ el tamaño $$$N$ , que no tiene ceros en la diagonal, se denomina conjunto de nodos $$$V = \ {0, ..., N-1 \}$ y bordes orientados $$$E = \ {(j, i) | L_ {ij} \ ne 0 \}$ .

Aquí, por ejemplo, es cómo se ve el gráfico de conexión para una matriz triangular inferior.

en el que los números simplemente corresponden al número ordinal del elemento diagonal, y los puntos indican componentes distintos de cero debajo de la diagonal:



Definición Gráfico orientado a la accesibilidad $$$G$ en múltiples índices $$$W \ subconjunto V$ llamar tantos índices $$$R_G (W) \ subconjunto V$ que en cualquier índice $$$z \ en R_G (W)$ se puede obtener pasando el gráfico desde algún índice $$$w (z) \ en W$ .

Ejemplo: para un gráfico de una imagen $$$R_G (\ {0, 4 \}) = \ {0, 4, 5, 6 \}$ . Está claro que siempre tiene $$$W \ subconjunto R_G (W)$ .

Si representamos cada nodo del gráfico como el número de columna de la matriz que lo generó, entonces los nodos vecinos a los que se dirigen sus bordes corresponden a los números de fila de componentes distintos de cero en esta columna.

Dejar $$$nnz (x) = \ {j | x_j \ ne 0 \}$ denota el conjunto de índices correspondientes a posiciones distintas de cero en el vector x.

Teorema de Hilbert (no, no por cuyo nombre se nombran los espacios) $$$nnz (x)$ donde $$$x$ hay un vector solución de un sistema triangular inferior escaso $$$Lx = f$ con escaso lado derecho $$$f$ coincide con $$$R_G (nnz (f))$ .

Suma por nuestra cuenta: en el teorema no tenemos en cuenta la improbable posibilidad de que los números distintos de cero en el lado derecho, al destruir las columnas, se puedan reducir a cero, por ejemplo, 3 - 3 = 0. Este efecto se llama cancelación numérica. Tener en cuenta tales ceros que ocurren espontáneamente es una pérdida de tiempo, y se perciben como todos los demás números en posiciones distintas de cero.

Un método efectivo para realizar un movimiento directo con un lado derecho escaso dado, bajo el supuesto de que tenemos acceso directo a sus componentes distintos de cero por índices, es el siguiente.

1. Recorremos el gráfico utilizando el método de "búsqueda en profundidad", comenzando secuencialmente desde el índice $$$i \ en nnz (f)$ cada componente distinto de cero del lado derecho. Escribir nodos encontrados en una matriz $$$R_G (nnz (f))$ mientras hacemos esto en el orden en que volvemos a la gráfica! Por analogía con el ejército de invasores: ocupamos el país sin crueldad, pero cuando comenzamos a regresar, nosotros, enojados, destruimos todo a su paso.
   Vale la pena señalar que no es absolutamente necesario que la lista de índices $$$nnz (f)$ se ordenó de forma ascendente cuando se introdujo en la entrada del algoritmo de "búsqueda en profundidad". Puedes comenzar en cualquier orden en el set $$$nnz (f)$ . Secuencia diferente perteneciente al conjunto $$$nnz (f)$ Los índices no afectan la solución final, como veremos en el ejemplo.
   
   Este paso no requiere ningún conocimiento de una matriz real. $$$v_A$ ! ¡Bienvenido al mundo del análisis simbólico cuando trabaje con solucionadores dispersos directos!
   
2. Procedemos a la solución en sí, teniendo a nuestra disposición una matriz $$$R_G (nnz (f))$ del paso anterior Destruimos las columnas en el orden inverso del registro de esta matriz . Voila!

Ejemplo

Considere un ejemplo en el que ambos pasos se demuestran en detalle. Supongamos que tenemos una matriz de 12 x 12 de la siguiente forma:

El gráfico de conectividad correspondiente tiene la forma:

Deje que el distinto de cero en la parte derecha esté solo en las posiciones 3 y 5, es decir $$$nnz (f) = \ {3, 5 \}$ . Repasemos el gráfico, comenzando por estos índices en el orden escrito. Entonces el método de "búsqueda profunda" se verá de la siguiente manera. Primero visitaremos los índices en orden $$$3 \ a 8 \ a 11$ , sin olvidar marcar los índices como visitados. En la siguiente figura están sombreados en azul:

Al regresar, ingresamos índices en nuestra matriz de índices de componentes de soluciones diferentes de cero $$$nnz (x): = \ {\ color {blue} {11}, \ color {blue} 8, \ color {blue} 3 \}$ . Luego, intenta correr $$$5 \ to8 \ to ...$ , pero nos encontramos con un nodo 8 ya marcado, por lo que no tocamos esta ruta y no vamos a la sucursal $$$5 \ to9 \ to ...$ . El resultado de esta carrera será $$$5 \ to9 \ to10$ . No podemos visitar el Nodo 11, ya que ya está etiquetado. Entonces, regresa y complementa la matriz $$$nnz (x)$ nueva captura marcada en verde: $$$nnz (x): = \ {\ color {blue} {11}, \ color {blue} 8, \ color {blue} 3, \ color {green} {10}, \ color {green} 9, \ color {verde} 5 \}$ . Y aquí está la imagen:

Los nodos verdes fangosos 8 y 11 son los que queríamos visitar durante la segunda carrera, pero no pudimos, porque Ya visitó durante el primero.

Entonces la matriz $$$nnz (x) = R_G (nnz (f))$ formado Pasamos al segundo paso. Un paso simple: corremos a través de la matriz $$$nnz (x)$ en el orden inverso (de derecha a izquierda), encontrando los componentes correspondientes del vector solución $$$x$ división por los componentes diagonales de la matriz $$$L$ y moviendo las columnas a la derecha. Otros componentes $$$x$ como eran ceros, permanecieron así. Esquemáticamente, esto se muestra a continuación, donde la flecha inferior indica el orden en que se destruyen las columnas:

Nota: en este orden de destrucción de columnas, el número 3 se encontrará más tarde que los números 5, 9 y 10. Las columnas no se destruyen en orden ascendente, lo que sería un error para las densas, es decir matrices incompletas Pero para las matrices dispersas, esto es común. Como se puede ver en la estructura matricial distinta de cero en este ejemplo, el uso de las columnas 5, 9 y 10 a la columna 3 no distorsionará los componentes en la respuesta $$$x_5$ , $$$x_9$ , $$$x_ {10}$ en absoluto, tienen c $$$x_3$ simplemente "sin intersecciones". ¡Nuestro método tomó esto en cuenta! Del mismo modo, usar la columna 8 después de la columna 9 no estropeará el componente $$$x_9$ . Gran algoritmo, ¿no?

Si vamos alrededor del gráfico primero desde el índice 5, y luego desde el índice 3, entonces nuestra matriz será $$$nnz (x): = \ {\ color {blue} {11}, \ color {blue} 8, \ color {blue} {10}, \ color {blue} {9}, \ color {blue} 5, \ color {verde} 3 \}$ . ¡Destruir las columnas en el orden inverso de esta matriz dará exactamente la misma solución que en el primer caso!

¡La complejidad de todas nuestras operaciones se escala por el número de operaciones aritméticas reales y el número de componentes distintos de cero en el lado derecho, pero no por el tamaño de la matriz! Hemos logrado nuestro objetivo.

Critica

SIN EMBARGO! Un lector crítico notará que la suposición misma al principio de que tenemos "acceso directo a los componentes distintos de cero del lado derecho por índice" ya significa que una vez antes, ejecutamos el lado derecho de arriba a abajo para encontrar estos índices y organizarlos en la matriz $$$nnz (f)$ es decir, ya pasé $$$O (N)$ acción! Además, la ejecución del gráfico en sí requiere que primero asignemos la memoria de la longitud máxima posible (¡en otro lugar debe escribir los índices encontrados al buscar en profundidad!) Para no sufrir una nueva reasignación a medida que crece la matriz $$$nnz (x)$ , y esto también requiere $$$O (N)$ operaciones! ¿Por qué, entonces, dicen, todo el alboroto?

Pero, de hecho, para una solución única de un sistema triangular inferior escaso con un lado derecho escaso, originalmente definido como un vector denso, no tiene sentido perder el tiempo del desarrollador en todos los algoritmos mencionados anteriormente. Pueden ser incluso más lentos que el método de la frente presentado por el algoritmo 2 anterior. Pero, como ya se mencionó, este aparato es indispensable para la factorización de Cholesky, por lo que no debe arrojarme tomates. De hecho, antes de ejecutar el algoritmo 3, toda la memoria necesaria de longitud máxima se asigna inmediatamente y requiere $$$O (N)$ tiempo En otro ciclo de columna $$$A$ todos los datos solo se sobrescriben en una matriz de longitud fija $$$N$ , y solo en aquellas posiciones en las que esta reescritura es relevante, gracias al acceso directo a elementos distintos de cero. ¡Y es precisamente por esto que surge la eficiencia!

Implementación de C ++

El gráfico en sí como una estructura de datos en el código no es necesario para construir. Es suficiente usarlo implícitamente cuando se trabaja con la matriz directamente. El algoritmo tendrá en cuenta toda la conectividad requerida. La búsqueda de profundidad, por supuesto, se implementa convenientemente mediante recursión banal. Aquí, por ejemplo, parece una búsqueda de profundidad recursiva basada en un único índice de inicio:

 /* j -   iA, jA -    ,    CSC top -     nnz(x);     0 result -   N,    nnz(x)   0  top-1  marked -    N;      */ void DepthFirstSearch(size_t j, const std::vector<size_t>& iA, const std::vector<size_t>& jA, size_t& top, std::vector<size_t>& result, std::vector<int>& marked) { size_t p; marked[j] = 1; //   j   for (p = jA[j]; p < jA[j+1]; p++) //       j { if (!marked[iA[p]]) //  iA[p]   { DepthFirstSearch(iA[p], iA, jA, top, result, marked); //      iA[p] } } result[top++] = j; //  j   nnz(x) } 

Si pasa la variable top igual a cero en la primera llamada a DepthFirstSearch , luego de completar todo el procedimiento recursivo, la variable top será igual al número de índices encontrados en la matriz de result . Tarea: escriba otra función que tome una matriz de índices de componentes distintos de cero en el lado derecho y los pase secuencialmente a la función DepthFirstSearch . Sin esto, el algoritmo no está completo. Tenga en cuenta: una serie de números reales $$$vA$ no pasamos a la función en absoluto, porque no es necesario en el proceso de análisis simbólico.

A pesar de su simplicidad, es obvio un defecto en la implementación: para sistemas grandes, los desbordamientos de la pila están a la vuelta de la esquina. Bueno, entonces hay una opción basada en un bucle en lugar de recurrencia. Es más difícil de entender, pero ya es adecuado para todas las ocasiones:

 /* j -   N -   iA, jA -    ,    CSC top -     nnz(x) result -   N,     nnz(x)   top  ret-1 ,  ret -    marked -    N work_stack -     N */ size_t DepthFirstSearch(size_t j, size_t N, const std::vector<size_t>& iA, const std::vector<size_t>& jA, size_t top, std::vector<size_t>& result, std::vector<int>& marked, std::vector<size_t>& work_stack) { size_t p, head = N - 1; int done; result[N - 1] = j; //    while (head < N) { j = result[head]; //  j     if (!marked[j]) { marked[j] = 1; //   j   work_stack[head] = jA[j]; } done = 1; //    j      for (p = work_stack[head] + 1; p < jA[j+1]; p++) //    j { //   c   iA[p] if (marked[iA[p]]) continue; //  iA[p]  ,   work_stack[head] = p; //       j result[--head] = iA[p]; //       iA[p] done = 0; //   j    break; } if (done) //      j  { head++; //  j    result[top++] = j; //  j    } } return (top); } 

Y así es como se ve el generador de una estructura vectorial de solución distinta de cero $$$nnz (x)$ :

 /* iA, jA -   ,    CSC iF -        f nnzf -      f result -   N,    nnz(x)   0  ret-1 ,  ret -    marked -    N,    work_stack -     N */ size_t NnzX(const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<size_t>& iF, size_t nnzf, std::vector<size_t>& result, std::vector<int>& marked, std::vector<size_t>& work_stack) { size_t p, N, top; N = jA.size() - 1; top = 0; for (p = 0; p < nnzf; p++) if (!marked[iF[p]]) //        top = DepthFirstSearch(iF[p], N, iA, jA, vA, top, result, marked, work_stack); for (p = 0; p < top; p++) marked[result[p]] = 0; //   return (top); } 

Finalmente, combinando todo en uno, escribimos el solucionador triangular inferior para el caso del lado derecho enrarecido:

 /* iA, jA, vA -  ,    CSC iF -        f nnzf -      f vF -       f result -   N,    nnz(x)   0  ret-1 ,  ret -    marked -    N,    work_stack -     N x -    N,    ;     */ size_t lower_triangular_solve (const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<double>& vA, const std::vector<size_t>& iF, size_t nnzf, const std::vector<double>& vF, std::vector<size_t>& result, std::vector<int>& marked, std::vector<size_t>& work_stack, std::vector<double>& x) { size_t top, p, j; ptrdiff_t px; top = NnzX(iA, jA, iF, nnzf, result, marked, work_stack); for (p = 0; p < nnzf; p++) x[iF[p]] = vF[p]; //    for (px = top; px > -1; px--) //     { j = result[px]; // x(j)   x [j] /= vA[jA[j]]; //   x(j) for (p = jA[j]+1; p < jA[j+1]; p++) { x[iA[p]] -= vA[p]*x[j]; //  j-  } } return (top) ; } 

Vemos que nuestro ciclo se ejecuta solo a través de los índices de la matriz $$$nnz (x)$ , no todo el conjunto $$$0, 1, ..., N - 1$ . Hecho

Hay una implementación que no utiliza una matriz de etiquetas markedpara ahorrar memoria. En cambio, se utiliza un conjunto adicional de índices.$$$V_1$ no se cruza con el conjunto $$$V = \{0, ..., N - 1\}$en el que hay un mapeo uno a uno mediante una operación algebraica simple como un procedimiento de marcado de nodos. Sin embargo, en nuestra era de memoria barata, guardarla en una sola matriz de longitud$$$N$ Parece completamente redundante.

En conclusión

El proceso de resolver un sistema disperso de ecuaciones lineales por el método directo, como regla, se divide en tres etapas:

1. Análisis de personajes
2. Factorización numérica basada en los datos del análisis sivol
3. La solución de los sistemas triangulares obtenidos con un denso lado derecho

El segundo paso, la factorización numérica, es la parte que consume más recursos y devora la mayor parte (> 90%) del tiempo estimado. El propósito del primer paso es reducir el alto costo del segundo. Un ejemplo de análisis simbólico se presentó en esta publicación. Sin embargo, es el primer paso que requiere el mayor tiempo de desarrollo y los costos mentales máximos por parte del desarrollador. Un buen análisis simbólico requiere el conocimiento de la teoría de gráficos y árboles y la posesión de un "instinto algorítmico". El segundo paso es incomparablemente más fácil de implementar.
Bueno, el tercer paso, tanto en la implementación como en el tiempo estimado, es en la mayoría de los casos el más sencillo.

Una buena introducción a los métodos directos para sistemas dispersos se encuentra en Tim Davis , profesor del Departamento de Informática de la Universidad de Texas A&M , titulado "Métodos directos para sistemas lineales dispersos ". El algoritmo descrito anteriormente se describe allí. No estoy familiarizado con Davis, aunque yo solía trabajar en la misma universidad (aunque en una facultad diferente). Si no me equivoco, el propio Davis participó en el desarrollo de solucionadores (!). utilizado en Matlab por otra parte, se vio involucrado incluso para los generadores de energía foto en Google Street View (MCO) por cierto, en la misma facultad que aparece ninguno otro que el propio Bjarne Stroustrup - .. el creador de C ++

Tal vez algún día Escribiré algo más sobre matrices dispersas o métodos numéricos en general Bueno para todos!