# Inaccurate Derivatives from Adjoint method for ODE-constrained problem

I have this very simple ODE-contrained optimization problem:

1. $h(x,x',p,t) = x'-A(p)x-b(p) = 0$, the constraint
2. $g(x(0)) = x_0$, the initial condition with no parameters involved
3. $F = \int (X-X_{obs})^2 dt$, the objective equation

According to adjoint method, I need to

1. Integrate constraint equation: $$x'=A(p)x+b(p)$$

2. Integrate adjoint equation and reverse $\lambda$ in $t$: $$\lambda'= A(p)^T \lambda-(X-X_{obs})$$

3. Calculate $\frac{dF}{dp}$: $$\frac{dF}{dp} = \int \lambda^T \frac{\partial h}{\partial p} dt,$$ since $$\frac{\partial f}{\partial p} = \frac{\partial g}{\partial p}=0$$

But for parameters only show up in $\mathbf{b(p)}$ term, derivatives from adjoint method is inconsistent with derivatives estimated using $\frac{\partial F}{\partial p}$ directly while derivatives for other parameters seem OK.

I’m thinking that this inconsistency maybe due to the fact that parameters in $b(p)$ doesn’t affect the calculation of $\lambda$ directly, namely it doesn’t show up in the adjoint equation? But there is also the possibility that I did something wrong in coding.

Any body have any similar experience? Thanks!

• I edited your post to use LaTeX/mathjax. I hope this is OK with you, and that I didn't change any of the meaning. Aug 28, 2015 at 18:07
• Your formula for the derivative $dF/dp$ is incorrect. I'll try to post a more detailed derivation soon, but in the meanwhile you should revisit steps 2. and 3... Aug 28, 2015 at 18:20
• Maybe I missed something, but other than the missing $\lambda$ after $A(p)^T$ in the adjoint equation (which must be a typo, otherwise the code wouldn't run), I don't see any errors in the equations. Aug 28, 2015 at 19:00
• @NickAlger I think the $\partial h/\partial p$ part is wrong, equation no 3 should read $\displaystyle \frac{\partial F}{\partial p} = \int_{0}^T \lambda^T [A'(p) x + b'(p)] dt$, where $x$ is the forward trajectory. Also, there should be a negative sign in front of $A(p)^T$ in the adjoint equation. Aug 28, 2015 at 19:20
• @user2186862 Doesn't that end out being the same thing? $\frac{\partial h}{\partial p}|_{x,p} = - A'(p)x - b'(p)$. Aug 28, 2015 at 19:37

As discussed in the comments, the equations look fine to me, with the exception of a missing $\lambda$ in the adjoint equation. It should read, $$\lambda'= A(p)^T\lambda-(X-X_{obs})$$

I assume the missing $\lambda$ in the original post is a typo, otherwise you'd be adding a matrix to a vector so your code wouldn't even run.

Since you haven't posted your formula for $\frac{\partial h}{\partial p}$, there could be a mistake there. It should be, $$\frac{\partial h}{\partial p} = -A'(p)x - b'(p).$$

If you forgot the $b'(p)$ term here, it might cause the issues you describe.

Otherwise, if the error is large there is probably a bug in the code somewhere.

If the error is small, then it could be a optimize-then-discretize vs. discretize-then-optimize issue. These operations only commute if you are discretizing with a Galerkin method. You can check if this is the problem by using very find discretizations in space and time, and seeing if the gradient error goes away.

• I only meant I missed it in this typing but not in code, also thanks for your editing! I'm pretty sure I included b'(p) term in the derivative cause otherwise derivative for that part is zero. The difference is 2-5 times underestimated and this only occurs for these parameters only showing up in b(p) term. If it's not a theoretical problem, then there must be something wrong with my code. I will check more. Aug 28, 2015 at 20:16
• @NickAlger: three minor observations (1) $A(p)^T \lambda$ should be $-A'(p)^T \lambda$ in the adjoint equation; (2) the cost function should read $F = \frac{1}{2} \int (x - x_{obs}) dt$ (otherwise the forcing term adjoint equation should read $-2(x - x_{obs})$), and (3) the correct form of the gradient is $\frac{dF}{dp} = - \int \lambda^T \frac{dh}{dp}$. Sep 3, 2015 at 13:03

I run the model forward with artificial parameters and then used my 'adjoint method' code to retrieve parameters. I found that as long as I assigned a good enough initial parameter values, the 'adjoint method' code indeed gives derivatives as expected. So, the issue I encountered at first is due to, surprisingly, the inaccuracy of direct estimation of dF/dp with finite differences. This is of course because of the complexity/instability of my model equations, which is a very dangerous sign... Anyway, issue cleared. Thanks everybody!

• It is often useful to take a sequence of progressively smaller finite difference step sizes, and plot the error as a function of step size 'h' on a log-log plot. If things are working correctly, there should be 3 distinct regions: 1) for large h, a high error region due to nonlinearity, 2) for medium h, a region where the curve decreases linearly, and 3) for small h, a noisy region due to ill-conditioning/numerical precision. Basically it should look something like this: i.stack.imgur.com/LuJMP.png Aug 29, 2015 at 9:05