11

Given code that computes a function $f(x)$, automatic differentiation tools produce a code that can compute $f(x)$ and its derivatives at the same time. Solving a differential equation is an entirely different problem and AD doesn't solve differential equations (although AD tools are sometimes useful in connection with PDE constrained optimization.) AD ...


9

Source-to-source transformation is considered the gold standard in terms of performance. OO approaches seem to be almost as good, in that there are more OO packages out there, and performance is not mentioned as a significant drawback. If you find an OO library you like for the language you're working in, I'd use it first, and then figure out later if you ...


9

It's state-dependent control flow that's an issue. function f(u) z = 0.0 while z < 10.0 z += u u += z^2 end return u end What's the program for computing the derivative? Automatic differentiation would give you: function f(_u) u = (_u,1.0) # seed the input derivative for the jvp in direction of basis e1 z = (0.0,0.0) while z[1] &...


8

Julia has a whole ecosystem for generating sparsity patterns and doing sparse automatic differentiation in a way that mixes with scientific computing and machine learning (or scientific machine learning). Tools like SparseDiffTools.jl, ModelingToolkit.jl, and SparsityDetection.jl will do things like: Automatically find sparsity patterns from code Generate ...


5

The SINTEF Matlab Reservoir Simulation Toolbox includes a GPL-licensed AD library. The usage is mostly geared towards numerical applications in subsurface flow, but the library itself is usable for more general purposes. Here is a basic runthrough of your example as you would run it from the base directory of MRST: startup; % Load ad based module ...


5

So you said that you have a few thousand variables. AD will usually be slower than hand coded derivatives, and if you're using forward AD, then you essentially need to evaluate the cost function once for each design variable,s and you have thousands of them. Typically in these cases, people use the reverse mode of AD, which scales independent of the number ...


4

The paper you linked answers the question. Autodiff (or hand differentiation) can differentiate branched program statements. For example, limiters, entropy fixes branching in flux statements, and the like. It can be rather helpful for min max statements as well. You can see an example below: Function(Vn_bar, a_bar, ul, cl, ur, cr) lambda1 = abs(...


4

Is there any out-of-the-shelf method I can use in this case to solve the equation without fully calculating the Jacobian? Using any integrator for stiff ODEs where the implicit equation is solved via Newton-Krylov methods will only require Jacobian-vector products. So something like DASKR, IDA, etc. can be made to do this.


4

This article in SIAM Review describes how to implement operator-overloading automatic differentiation in MATLAB, and gives a good introduction to the technology.


3

This is what Griewank et al. call "Piecewise linearization in secant mode", see for instance https://opus4.kobv.de/opus4-zib/files/6164/newton_secant_approx_paper.pdf. The aim of that research was to capture the kinks of absolute value operations with the same precision a tangent or a secant captures the local behavior of a smooth function, with an ...


3

One approach is to use the Faà di Bruno formula: $$ \big((1+x^M)^{1/M}\big)^{(n)} = \sum \frac{n!}{m_1!1!^{m_1}\cdots m_n! n!^{m_n}} (1/M)^{\underline{m_1+\cdots+m_n}}(1+x^M)^{1/M-m_1-\cdots-m_n} \prod_{j=1}^{n} \big(M^{\underline{\phantom{,}j}} x^{M-j}\big)^{m_j}, \qquad M\notin\mathbb{Z} $$ Here $a^{\underline{b}}$ is the falling factorial, and the sum is ...


3

If you were using Julia, then SparsityDetection.jl solves (1) by using non-standard analysis with concolic fuzzing to directly determine the sparsity pattern from the function's AST, and SparseDiffTools.jl performs matrix coloring and colored automatic differentiation to do a sparse calculation of the Jacobian to solve (2). There are other systems that can ...


2

Are you sure it is supposed to be sign function ? According to this https://projects.coin-or.org/ADOL-C/browser/stable/2.1/ADOL-C/doc/adolc-manual.pdf?format=raw in section 1.8, condassign(a,b,c,d) is equal to a = (b > 0) ? c : d So the function posted is actually giving Heaviside function and not the sign function. It implements this function $$ f(x) =...


2

Potentially related and useful (I gave these resources to my students when teaching Intro to Computational Mathematics, kinda useful pedagogically too): "Automatic Source-to-Source Error Compensation of Floating-Point Programs" by Laurent Thévenoux, Philippe Langlois and Matthieu Martel : https://hal.archives-ouvertes.fr/hal-01158399/document ...


2

Assuming that the issue is that DAEPACK works in 32-bit mode but not 64-bit mode, here are some approaches to address this. Compile DAEPACK in 32-bit mode on a 64-bit OS See how to do this in another question. This may be an easy way to continue to work with it. If the issue is that the generated code is having issues in 64-bit mode, compile the generated ...


1

Define the iteration $$ x_{n+1} = f(x_n) $$ Then $$ e_n := \frac{\|x_{n+1} - x_n\|}{\| x_n - x_{n-1}\|} \le L $$ Compute $e_n$ and that should tell you something about L. If the iterations converge ($L < 1$), then $e_n$ should converge to $L$.


1

There's a good discussion of globalization strategies for Newton's method root finding in Numerical Methods for Unconstrained Optimization and Nonlinear Equations by J. E. Dennis, Jr. and Robert B. Schnabel. See section 6.5 in particular. Dennis and Schnabel advocate treating the problem of finding a root of $f(p)=0$ as a nonlinear least squares problem: $\...


1

I would also like to point at MatlabAutoDiff, which supports sparse Jacobians. Have tried it myself: it is possible to compute large Jacobians (tried with N=1e5) in a small amount of time.


1

So the size of the problem doesn't necessarily stop you from using the full Jacobian if you store it in a compressed/sparse storage format, you also should be able to know which residuals are dependent on which unknowns which should make it cheaper as well. That said, just use Jacobian free Newton-Krylov (JFNK) solvers. FGMRES is a really strong linear ...


1

An automatic differentiation program could be refactored to output a symbolic representation, instead of a numerical one. Therefore the 2 forms would be equivalent. Please see the non-accepted answer here: https://stackoverflow.com/a/55607008/104910


1

Starting with $b=Ax$ the forward differentiation gives $\dot b=\dot A x+A\dot x$, or, with $\dot y = \dot A x$: $$\dot x = -A^{-1} \dot y + A^{-1} \dot b$$ Let's write this in terms of the local transformation Jacobian for the tangent linear variables $\left( \dot A, \dot y, \dot b, \dot x\right)$: $$ \left[ \begin{matrix} \dot A \\ \dot y\\ \dot b\\ \...


1

The identity connecting forward and reverse mode for $x=\text{solve}(A,b)$ is $$\bar x^⊤\dot x=\bar A(\dot A)+\bar b^⊤\dot b.$$ Where $\bar A$ is a linear functional on matrices which may be realized as the $trace(\bar A^⊤\dot A)$ which is related to the (Frobenius) scalar product for matrices $\langle X,Y\rangle=trace(X^⊤Y)=trace(XY^⊤)$, which is kind-of ...


1

I need to calculate the Jacobian matrix of a subroutine F(U). Both F and U are of size $N(=O(10^5))$. Using Tapenade, I differentiated the routine in tangent mode. I cannot calculate the full Jacobian directly because of the large memory requirement. There are AD (automatic differentiation) software packages that can do this calculation; I know DAEPACK can, ...


1

I realize this is an old question, but when looking for this myself today I found ADiGator, which is open source, and seems to handle vectors. I haven't tested it yet myself, but it seems to be actively developed.


1

For gradient computation, you use the reverse mode of AD. This requires in both cases to build an operand stack, the OO version also needs to build up an operations stack, which has to be interpreted in the reverse traversal of the code. Source transformed code writes out the reverse-ordered operations as additional source code that is compiled. The overhead ...


1

Based on so many simulations that I've run using my manual divergence calculation. The answer is yes. I also compared manual calculation with the built-in one and no error or different was observed.


Only top voted, non community-wiki answers of a minimum length are eligible