¿Cómo resuelve el operador de barra invertida de MATLAB

36

Estaba comparando algunos de mis códigos con los códigos MATLAB "de inventario". Estoy sorprendido por los resultados.

Ejecuté un código de muestra (matriz dispersa)

n = 5000;
a = diag(rand(n,1));
b = rand(n,1);
disp('For a\b');
tic;a\b;toc;
disp('For LU');
tic;LULU;toc;
disp('For Conj Grad');
tic;conjgrad(a,b,1e-8);toc;
disp('Inv(A)*B');
tic;inv(a)*b;toc;

Resultados:

    For a\b
    Elapsed time is 0.052838 seconds.

    For LU
    Elapsed time is 7.441331 seconds.

    For Conj Grad
    Elapsed time is 3.819182 seconds.

    Inv(A)*B
    Elapsed time is 38.511110 seconds.

Para matriz densa:

n = 2000;
a = rand(n,n);
b = rand(n,1);
disp('For a\b');
tic;a\b;toc;
disp('For LU');
tic;LULU;toc;
disp('For Conj Grad');
tic;conjgrad(a,b,1e-8);toc;
disp('For INV(A)*B');
tic;inv(a)*b;toc;

Resultados:

For a\b
Elapsed time is 0.575926 seconds.

For LU
Elapsed time is 0.654287 seconds.

For Conj Grad
Elapsed time is 9.875896 seconds.

Inv(A)*B
Elapsed time is 1.648074 seconds.

¿Cómo diablos es a \ b tan increíble?

linear-algebra performance matlab Encuesta
fuente
1
La barra invertida incorporada de MATLAB, en otras palabras, el solucionador directo para un sistema de ecuaciones lineales, utiliza el método Multifrontal para matriz dispersa, es por eso que A \ B es tan impresionante.
Shuhao Cao
1
Utiliza el código de Tim Davis disponible en cise.ufl.edu/research/sparse . También la genialidad desaparece cuando tienes un problema no trivial
stali
1
¿Qué es "LULU"? ¿Por qué crees que es una buena implementación de una factorización LU y su posterior resolución directa?
Jed Brown
3
@Nunoxic ¿Qué implementación? ¿Lo escribiste tú mismo? El álgebra lineal densa de alto rendimiento, aunque generalmente se entiende algorítmicamente, no es fácil de implementar de manera eficiente en el hardware moderno. Las mejores implementaciones de BLAS / Lapack deberían estar cerca del pico para una matriz de ese tamaño. Además, a partir de sus comentarios, tengo la impresión de que cree que LU y Gaussian Elimination son algoritmos diferentes.
Jed Brown
1
Llama a un código Fortran escrito con Intel MKL.
Investigación

Respuestas:

37

En Matlab, el comando '\' invoca un algoritmo que depende de la estructura de la matriz A e incluye controles (pequeña sobrecarga) en las propiedades de A.

  1. Si A es escaso y con bandas, emplee un solucionador con bandas.
  2. Si A es una matriz triangular superior o inferior, emplee un algoritmo de sustitución hacia atrás.
  3. Si A es simétrico y tiene elementos diagonales positivos reales, intente una factorización de Cholesky. Si A es escaso, emplee la reordenación primero para minimizar el relleno.
  4. Si ninguno de los criterios anteriores se cumple, realice una factorización triangular general utilizando la eliminación gaussiana con pivote parcial.
  5. Si A es escaso, entonces emplee la biblioteca UMFPACK.
  6. Si A no es cuadrado, emplee algoritmos basados ​​en la factorización QR para sistemas indeterminados.

Para reducir los gastos generales, es posible utilizar el comando linsolve en Matlab y seleccionar un solucionador adecuado entre estas opciones usted mismo.

Allan P. Engsig-Karup
fuente
Suponiendo que estoy tratando con una matriz densa no estructurada de 10000x10000 con todos los elementos distintos de cero (Alto nivel de densidad), ¿cuál sería mi mejor apuesta? Quiero aislar ese algoritmo 1 que funciona para matrices densas. ¿Es LU, QR o eliminación gaussiana?
Investigación
1
Suena como un Paso 4 donde se invoca Gaussian Elimination que corresponde al caso más general donde no se puede explotar ninguna estructura de A para aumentar el rendimiento. Entonces, básicamente se trata de una factorización LU y una posterior hacia adelante seguida de un paso de sustitución hacia atrás.
Allan P. Engsig-Karup
¡Gracias! Creo que eso me da una dirección para pensar. Actualmente, la Eliminación Gaussiana es lo mejor que tenemos para resolver problemas tan desestructurados, ¿es correcto?
Investigación
37

Si desea ver qué a\bhace para su matriz particular, puede establecer spparms('spumoni',1)y calcular exactamente qué algoritmo le impresionó. Por ejemplo:

spparms('spumoni',1);
A = delsq(numgrid('B',256));
b = rand(size(A,2),1);
mldivide(A,b);  % another way to write A\b

saldrá

sp\: bandwidth = 254+1+254.
sp\: is A diagonal? no.
sp\: is band density (0.01) > bandden (0.50) to try banded solver? no.
sp\: is A triangular? no.
sp\: is A morally triangular? no.
sp\: is A a candidate for Cholesky (symmetric, real positive diagonal)? yes.
sp\: is CHOLMOD's symbolic Cholesky factorization (with automatic reordering) successful? yes.
sp\: is CHOLMOD's numeric Cholesky factorization successful? yes.
sp\: is CHOLMOD's triangular solve successful? yes.

entonces puedo ver que "\" terminó usando "CHOLMOD" en este caso.

dranxo
fuente
3
+1 para la nueva configuración de MATLAB de la que nunca había oído hablar. Bien jugado, señor.
Geoff Oxberry
2
¡Hey gracias! Está en help mldivide.
dranxo
16

Para matrices dispersas, Matlab usa UMFPACK para la \operación " ", que, en su ejemplo, básicamente ejecuta los valores de a, los invierte y los multiplica con los valores de b. Sin embargo, para este ejemplo, debe usar b./diag(a), que es mucho más rápido.

Para sistemas densos, el operador de barra invertida es un poco más complicado. Una breve descripción de lo que se hace cuando se da aquí . De acuerdo con esa descripción, en su ejemplo, Matlab resolvería a\busando la sustitución hacia atrás. Para las matrices cuadradas generales, se utiliza una descomposición LU.

Pedro
fuente
Regd. Sparsity, el inv de una matriz diag sería simplemente el recíproco de los elementos diagonales, por lo que b./diag(a) funcionaría, pero a \ b también funciona de maravilla para matrices dispersas generales. ¿Por qué no es linsolve o LULU (Mi versión optimizada de LU) más rápido que a \ b para matrices densas en ese caso?
Investigación
@Nunoxic ¿Su LULU hace algo para detectar la diagonalidad o la triangularidad de la matriz densa? De lo contrario, tomará el mismo tiempo con cada matriz, independientemente de su contenido o estructura.
Pedro
Algo. Pero, con los indicadores OPT de linsolve, definí todo lo que hay que definir sobre la estructura. Sin embargo, a \ b es más rápido.
Investigación del
@Nunoxic, cada vez que llama a una función de usuario, genera gastos generales. Matlab hace todo para la barra invertida internamente, por ejemplo, la descomposición y la posterior aplicación del lado derecho, con muy poca sobrecarga, y por lo tanto será más rápido. Además, en sus pruebas, debe usar más de una llamada para obtener tiempos confiables, por ejemplo tic; for k=1:100, a\b; end; toc.
Pedro
5

Como regla general, si tiene una matriz escasa de complejidad razonable (es decir, no tiene que ser la plantilla de 5 puntos, pero de hecho puede ser una discretización de las ecuaciones de Stokes para las cuales el número de no ceros por fila es mucho más grande que 5), entonces un solucionador directo escaso como UMFPACK generalmente supera a un solucionador iterativo de Krylov si el problema no es mayor que alrededor de unas 100,000 incógnitas.

En otras palabras, para la mayoría de las matrices dispersas resultantes de discretizaciones 2D, un solucionador directo es la alternativa más rápida. Solo para los problemas 3D en los que superas rápidamente las 100.000 incógnitas, es imprescindible utilizar un solucionador iterativo.

Wolfgang Bangerth
fuente
3
No tengo claro cómo responde esto a la pregunta, pero también estoy en desacuerdo con la premisa. Es cierto que los solucionadores directos tienden a funcionar bien para problemas 2D de tamaño modesto, pero los solucionadores directos tienden a sufrir en 3D mucho antes de las 100k incógnitas (los separadores de vértices son mucho más grandes que en 2D). Además, afirmo que en la mayoría de los casos (p. Ej., Operadores elípticos), los solucionadores iterativos se pueden hacer más rápido que los solucionadores directos, incluso para problemas 2D de tamaño moderado, pero puede requerir un esfuerzo considerable (p. Ej., Usar los ingredientes correctos para hacer que la cuadrícula funcione) .
Jed Brown
1
Si sus problemas son razonablemente pequeños y no quiere pensar en diseñar solucionadores implícitos, o si sus problemas son muy difíciles (como Maxwell de alta frecuencia) y no desea dedicar su carrera a diseñar buenos solucionadores, entonces yo De acuerdo en que los solucionadores directos dispersos son una gran opción.
Jed Brown
En realidad, mi problema es diseñar un buen solucionador denso. No tengo una aplicación particular como tal (Por lo tanto, no hay estructura). Quería modificar algunos de mis códigos y hacerlos competitivos con lo que se usa actualmente. Tenía curiosidad acerca de a \ by cómo siempre supera a mi código.
Investigación
@JedBrown: Sí, tal vez con una cantidad significativa de esfuerzo uno puede obtener solucionadores iterativos para vencer a un solucionador directo incluso para pequeños problemas 2D. ¿Pero por qué hacerlo? Para problemas con <100k incógnitas, los solucionadores directos dispersos en 2d son bastante rápidos. Lo que quería decir es: para problemas tan pequeños, no pierdas el tiempo jugando y descubriendo la mejor combinación de parámetros, solo toma la caja negra.
Wolfgang Bangerth
Ya existe un orden de magnitud en la diferencia de tiempo de ejecución entre una escasa cuadrícula geométrica de Cholesky y (basada en matriz) para problemas 2D "fáciles" con dofs de 100k usando una plantilla de 5 puntos (~ 0.5 segundos en comparación con ~ 0.05 segundos). Si su plantilla utiliza segundos vecinos (por ejemplo, discretización de cuarto orden, Newton para algunas opciones de reología no lineal, estabilización, etc.), el tamaño de los separadores de vértices se duplica aproximadamente, por lo que el costo de la resolución directa aumenta aproximadamente 8x (el costo es más dependiente del problema para MG). Con muchos pasos de tiempo o ciclos de optimización / UQ, estas diferencias son significativas.
Jed Brown