Como lidar com a complexidade no código numérico, por exemplo, ao lidar com grandes matrizes jacobianas?

Estou resolvendo um sistema não linear de equações acopladas e calculei o jacobiano do sistema discretizado. O resultado é realmente complicado, abaixo estão (apenas!) As 3 primeiras colunas de uma matriz ,$3\times 9$

(A complexidade surge, em parte, porque o esquema numérico requer ajuste exponencial para estabilidade.)

Eu tenho uma pergunta bastante geral sobre a implementação de códigos numéricos usando jacobianos.

Eu posso ir em frente e implementar essa matriz no código. Mas minha intuição está me dizendo para esperar alguns dias (talvez semanas!) De depuração tediosa devido à pura complexidade e à inevitabilidade de introduzir erros. Como lidar com uma complexidade como essa no código numérico, parece inevitável ?! Você usa a geração automática de código a partir de pacotes simbólicos (depois ajusta o código manualmente)?

Primeiro, planejo depurar o jacobiano analítico com uma aproximação de diferença finita, devo estar ciente de alguma armadilha? Como você lida com problemas semelhantes no seu código?

Atualizar

Estou codificando isso em Python e usei o sympy para gerar o Jacobiano. Talvez eu possa usar o recurso de geração de código ?

Uma palavra: Modularidade .

Há muitas expressões repetidas no seu jacobiano que podem ser escritas como sua própria função. Não há motivo para você escrever a mesma operação mais de uma vez, e isso facilitará a depuração; se você escrever apenas uma vez, haverá apenas um lugar para um erro (em teoria).

O código modular também facilitará o teste; Você pode escrever testes para cada componente do seu jacobiano, em vez de tentar testar toda a matriz. Por exemplo, se você escrever sua função am () de forma modular, poderá escrever facilmente testes de sanidade para isso, verificar se a diferencia adequadamente, etc.

Outra sugestão seria procurar bibliotecas de diferenciação automática para montar os jacobianos. Não há garantia de que eles estejam livres de erros, mas provavelmente haverá menos erros de depuração / menos do que escrever seus próprios. Aqui estão alguns que você pode querer dar uma olhada:

• Sacado (Sandia Labs)
• ADIC (Argonne)

Desculpe, acabei de ver que você está usando python. ScientificPython tem suporte para AD.

• ScientificPython

Deixe-me pesar aqui com algumas palavras de cautela, precedidas de uma história. Há muito tempo, eu trabalhava com um colega quando estava começando. Ele tinha um problema de otimização a resolver, com um objetivo bastante confuso. Sua solução foi gerar os derivados analíticos para uma otimização.

O problema que eu vi foi que esses derivados eram desagradáveis. Gerados usando Macsyma, convertidos em código fortran, cada um deles tinha dezenas de instruções de continuação. De fato, o compilador Fortran ficou chateado com isso, pois excedeu o número máximo de instruções de continuação. Embora tenhamos encontrado uma bandeira que nos permitiu contornar esse problema, houve outros problemas.

• Em expressões longas, como geralmente são geradas pelos sistemas de CA, existe o risco de cancelamento subtrativo maciço. Calcule muitos números grandes, apenas para descobrir que todos se cancelam para produzir um número pequeno.

• Geralmente, os derivativos gerados analiticamente são mais caros para avaliar do que os derivados numericamente gerados usando diferenças finitas. Um gradiente para n variáveis ​​pode levar mais que n vezes o custo de avaliar sua função objetivo. (Você pode economizar algum tempo porque muitos dos termos podem ser reutilizados nos vários derivados, mas isso também o forçará a fazer uma codificação manual cuidadosa, em vez de usar expressões geradas por computador. E sempre que codificar matemática desagradável expressões, a probabilidade de um erro não é trivial. Verifique se essas derivações estão corretas.)

O ponto da minha história é que essas expressões geradas pela CA têm problemas próprios. O engraçado é que meu colega estava realmente orgulhoso da complexidade do problema, que estava resolvendo claramente um problema realmente difícil porque a álgebra era muito desagradável. O que eu não acho que ele considerou foi se essa álgebra estava realmente computando a coisa correta, estava fazendo com precisão e com eficiência.

Se eu fosse a pessoa mais velha na época nesse projeto, teria lido o ato de revolta. Seu orgulho o levou a usar uma solução que provavelmente era desnecessariamente complexa, sem ao menos verificar se um gradiente finito baseado em diferenças era adequado. Aposto que gastamos talvez uma semana de trabalho para executar essa otimização. No mínimo, eu o aconselharia a testar cuidadosamente o gradiente produzido. Foi preciso? Qual foi a precisão, em comparação com os derivados de diferenças finitas? De fato, existem ferramentas hoje em dia que também retornarão uma estimativa do erro em sua previsão derivada. Isso certamente é verdade para o código de diferenciação adaptável (derivest) que escrevi no MATLAB.

Teste o código. Verifique os derivados.

Mas antes de fazer QUALQUER disso, considere se outros esquemas de otimização melhores são uma opção. Por exemplo, se você estiver ajustando exponencialmente, há uma chance muito boa de poder usar mínimos quadrados não-lineares particionados (às vezes chamados de mínimos quadrados separáveis. Acho que esse foi o termo usado por Seber e Wild em seu livro.) é dividir o conjunto de parâmetros em conjuntos intrinsecamente lineares e intrinsecamente não lineares. Use uma otimização que funcione apenas nos parâmetros não lineares. Dado que esses parâmetros são "conhecidos", os parâmetros intrinsecamente lineares podem ser estimados usando mínimos quadrados lineares simples. Esse esquema reduzirá o espaço de parâmetros na otimização. Isso torna o problema mais robusto, pois você não precisa encontrar valores iniciais para os parâmetros lineares. Reduz a dimensionalidade do seu espaço de pesquisa, tornando o problema mais rápido. Mais uma vez eu forneciuma ferramenta para esse fim , mas apenas no MATLAB.

Se você usar os derivados analíticos, codifique-os para reutilizar termos. Isso pode economizar muito tempo e reduzir os bugs, economizando seu tempo. Mas então verifique esses números!

Existem várias estratégias a serem consideradas:

1. Encontre as derivadas na forma simbólica usando um CAS e exporte o código para calcular as derivadas.

2. Use uma ferramenta de diferenciação automática (AD) para produzir código que calcula as derivadas do código para calcular as funções.

3. Use aproximações de diferenças finitas para aproximar o jacobiano.

A diferenciação automática poderia produzir um código mais eficiente para computar todo o jacobiano e depois usar a computação simbólica para produzir uma fórmula para cada entrada na matriz. Diferenças finitas são uma boa maneira de verificar suas derivadas.

Aqui está um exemplo de onde usamos diferenciação automática usando o Sacado em um código: http://www.dealii.org/developer/doxygen/deal.II/step_33.html

Além das excelentes sugestões de BrianBorcher, outra abordagem possível para funções com valor real é usar a aproximação de derivada de etapas complexas (consulte este artigo (paywalled) e este artigo ). Em alguns casos, essa abordagem gera derivadas numéricas mais precisas ao custo de alterar os valores das variáveis ​​em sua função de real para complexo. O segundo artigo lista alguns casos em que a aproximação da função de etapa complexa pode ser interrompida.