¿Cómo lidiar con la complejidad en el código numérico, por ejemplo, cuando se trata de matrices jacobianas grandes?

Estoy resolviendo un sistema no lineal de ecuaciones acopladas, y he calculado el jacobiano del sistema discretizado. El resultado es realmente complicado, a continuación se encuentran (¡solo!) Las primeras 3 columnas de una matriz de $3\times 9$ ,

(La complejidad surge, en parte, porque el esquema numérico requiere un ajuste exponencial para la estabilidad).

Tengo una pregunta bastante general sobre la implementación de códigos numéricos usando jacobianos.

Puedo seguir adelante e implementar esta matriz en código. Pero mi intuición me dice que espere unos días (¡quizás semanas!) De depuración tediosa debido a la gran complejidad y la inevitabilidad de la introducción de errores. ¿Cómo hacer frente a una complejidad como esta en el código numérico, parece inevitable? ¿Utiliza la generación automática de código a partir de paquetes simbólicos (luego modifica el código a mano)?

Primero planeo depurar al analítico jacobiano con una aproximación de diferencia finita, ¿debo estar al tanto de cualquier error? ¿Cómo manejas problemas similares en tu código?

Actualizar

Estoy codificando esto en Python y usé sympy para generar el jacobiano. ¿Quizás pueda usar la función de generación de código ?

Una palabra: modularidad .

Hay muchas expresiones repetidas en su jacobiano que podrían escribirse como su propia función. No hay razón para que escriba la misma operación más de una vez, y eso facilitará la depuración; si solo lo escribe una vez, solo hay un lugar para un error (en teoría).

El código modular también facilitará las pruebas; Puede escribir pruebas para cada componente de su Jacobian en lugar de intentar probar toda la matriz. Por ejemplo, si escribe su función am () de forma modular, puede escribir fácilmente pruebas de cordura para ella, verificar si la está diferenciando correctamente, etc.

Otra sugerencia sería buscar bibliotecas de diferenciación automática para ensamblar el jacobiano. No hay garantía de que estén libres de errores, pero probablemente habrá menos depuración / menos errores que escribir el suyo. Aquí hay algunos que quizás desee ver:

• Sacado (Sandia Labs)
• ADIC (Argonne)

Lo siento, acabo de ver que estás usando Python. ScientificPython tiene soporte para AD.

• ScientificPython

Permítanme analizar aquí con algunas palabras de precaución, precedidas de una historia. Hace mucho tiempo, trabajé con un compañero cuando recién comenzaba. Tenía que resolver un problema de optimización, con un objetivo bastante desordenado. Su solución fue generar los derivados analíticos para una optimización.

El problema que vi fue que estos derivados eran desagradables. Generados usando Macsyma, convertidos a código fortran, cada uno tenía docenas de declaraciones de continuación. De hecho, el compilador de Fortran se molestó por eso, ya que excedió el número máximo de declaraciones de continuación. Si bien encontramos una bandera que nos permitió solucionar ese problema, hubo otros problemas.

• En expresiones largas, como las que comúnmente generan los sistemas de CA, existe el riesgo de una cancelación sustractiva masiva. Calcule muchos números grandes, solo para descubrir que todos se cancelan entre sí para producir un número pequeño.

• A menudo, los derivados generados analíticamente son en realidad más costosos de evaluar que los derivados generados numéricamente utilizando diferencias finitas. Un gradiente para n variables puede tomar más de n veces el costo de evaluar su función objetivo. (Es posible que pueda ahorrar algo de tiempo porque muchos de los términos se pueden reutilizar a través de varias derivadas, pero eso también lo obligará a realizar una codificación manual cuidadosa, en lugar de usar expresiones generadas por computadora. Y cada vez que codifique matemática desagradable expresiones, la probabilidad de un error no es trivial. Asegúrese de verificar la precisión de estos derivados).

El punto de mi historia es que estas expresiones generadas por CA tienen sus propios problemas. Lo curioso es que mi colega estaba realmente orgulloso de la complejidad del problema, que claramente estaba resolviendo un problema realmente difícil porque el álgebra era muy desagradable. Lo que no creo que haya considerado es si ese álgebra en realidad estaba calculando lo correcto, si lo estaba haciendo con precisión y si lo estaba haciendo de manera eficiente.

Si hubiera sido la persona mayor en el momento en este proyecto, le habría leído el acto antidisturbios. Su orgullo le hizo usar una solución que probablemente era innecesariamente compleja, sin siquiera comprobar que un gradiente basado en diferencias finitas era adecuado. Apuesto a que hemos pasado tal vez una semana hombre para hacer que esta optimización funcione. Por lo menos, le habría aconsejado que probara cuidadosamente el gradiente producido. ¿Fue exacto? ¿Qué tan preciso fue, en comparación con los derivados de diferencias finitas? De hecho, hoy existen herramientas que también devolverán una estimación del error en su predicción derivada. Esto es ciertamente cierto para el código de diferenciación adaptativa, (derivación) que he escrito en MATLAB.

Prueba el código. Verificar los derivados.

Pero antes de hacer CUALQUIERA de esto, considere si otros esquemas de optimización mejores son una opción. Por ejemplo, si está haciendo un ajuste exponencial, entonces hay una muy buena posibilidad de que pueda usar un mínimo cuadrado no lineal particionado (a veces llamado mínimo cuadrado separable. Creo que ese fue el término utilizado por Seber y Wild en su libro). La idea es dividir el conjunto de parámetros en conjuntos intrínsecamente lineales e intrínsecamente no lineales. Use una optimización que funcione solo en los parámetros no lineales. Dado que esos parámetros son "conocidos", entonces los parámetros intrínsecamente lineales pueden estimarse usando mínimos cuadrados lineales simples. Este esquema reducirá el espacio de parámetros en la optimización. Hace que el problema sea más robusto, ya que no necesita encontrar valores iniciales para los parámetros lineales. Reduce la dimensionalidad de su espacio de búsqueda, por lo que el problema se ejecuta más rápidamente. De nuevo he suministradouna herramienta para este propósito , pero solo en MATLAB.

Si utiliza los derivados analíticos, codifíquelos para reutilizar los términos. Esto puede ser un gran ahorro de tiempo, y en realidad puede reducir los errores, ahorrando su propio tiempo. ¡Pero luego revisa esos números!

Hay varias estrategias a considerar:

1. Encuentre los derivados en forma simbólica usando un CAS, luego exporte el código para calcular los derivados.

2. Use una herramienta de diferenciación automática (AD) para producir código que calcule las derivadas del código para calcular las funciones.

3. Use aproximaciones de diferencias finitas para aproximar el jacobiano.

La diferenciación automática podría producir un código más eficiente para calcular todo el jacobiano y luego usar el cálculo simbólico para producir una fórmula para cada entrada en la matriz. Las diferencias finitas son una buena manera de verificar sus derivados.

Aquí hay un ejemplo de dónde hemos usado la diferenciación automática usando Sacado en un código: http://www.dealii.org/developer/doxygen/deal.II/step_33.html

Además de las excelentes sugerencias de BrianBorcher, otro enfoque posible para las funciones de valor real es utilizar la aproximación derivada de pasos complejos (consulte este artículo (pago) y este artículo ). En algunos casos, este enfoque produce derivadas numéricas más precisas a costa de cambiar los valores de las variables en su función de real a complejo. El segundo artículo enumera algunos casos en los que la aproximación de la función de paso complejo podría romperse.