bifurcationkit
Suggested ehancements for BifFunction for Hessian and parameter gradient caching
Currently, the BifFunction stores the second and third differentials as functions.
For some large-scale problems that arise in PDE, it is actually way more efficient to actually precompute and cache a sparse/matrix-free Hessian/dF2 for a given (u, p) similarly as it is currently done for the Jacobian.
I propose adding a hessian field in the same manner as the Jacobian, which would supplement dF2.
bf.d2F(u, p, x1, x1) = bf.hessian(u, p)(x1, x2)
The same interface could be used for the third-order tensor even if its applications are far more rare.
It also seems that the computation of the vector field gradient relative to the continuation parameter can only be done internally using finite difference, but in most cases, an analytic computation is possible.
A gradient method in BifFunction could be also added. Because p is a named tuple, there are several possibilities to either:
-
Compute the gradient to all the parameter
bf.grad(u, p)which returns a named tuple or dict containing the gradient associated with each keys ofp. -
Or specify explicitly which parameter is differentiated with a symbol
bf.grad(u, p, p_name::Symbol)that only returns a single vector.
Alternatively, the gradient computation function could also be directly passed as a parameter for the PALC algorithm (but that doesn't work really well for codim2 continuation).
It also seems that the computation of the vector field gradient relative to the continuation parameter can only be done internally using finite difference, but in most cases, an analytic computation is possible.
Yes, I know. I have been thinking about this for long. The problem is a problem of API, for making it simple for the user to specify. By default, we could use AD for the parameter differential. But there are cases when you want to use AD for the jacobian and finite differences for the parameters. How would you put this in the constructor BifurcationProblem? Add a keyword use_ad_for_param = true?
Alternatively, the gradient computation function could also be directly passed as a parameter for the PALC algorithm (but that doesn't work really well for codim2 continuation).
For codim2, we would need another differential but that's OK because the BifurcationProblem is then transformed into FoldProblemMinimallyAugmented and we can dispach the parameter differential for the second parameter axis
A gradient method in BifFunction could be also added. Because p is a named tuple, there are several possibilities to either:
What is the usecase for bf.grad(u, p, p_name::Symbol) ? (I would prefer an @optic to specify which parameter)
Currently, the BifFunction stores the second and third differentials as functions.
For some large-scale problems that arise in PDE, it is actually way more efficient to actually precompute and cache a sparse/matrix-free Hessian/dF2 for a given (u, p) similarly as it is currently done for the Jacobian.
You mean for Fold continuation with :minaug linearsolver here? If so, are you sure it is slower? It is only called once at for a given (u,p)
It also seems that the computation of the vector field gradient relative to the continuation parameter can only be done internally using finite difference, but in most cases, an analytic computation is possible.
Yes, I know. I have been thinking about this for long. The problem is a problem of API, for making it simple for the user to specify. By default, we could use AD for the parameter differential. But there are cases when you want to use AD for the jacobian and finite differences for the parameters. How would you put this in the constructor
BifurcationProblem? Add a keyworduse_ad_for_param = true?
Alternatively, the continuation parameter in the PALC could be set to a Dual Number (enabling it with keyword like params_grad = :FD or params_grad = :Dual) for evaluating the tangent at the same time as the residual.
The user would need to implement a dual number dispatch for the most generic case but if it only common functions are used, it should work without needing any changes.
Currently, the BifFunction stores the second and third differentials as functions.
For some large-scale problems that arise in PDE, it is actually way more efficient to actually precompute and cache a sparse/matrix-free Hessian/dF2 for a given (u, p) similarly as it is currently done for the Jacobian.
You mean for Fold continuation with :minaug linearsolver here? If so, are you sure it is slower? It is only called once at for a given (u,p)
Yes, you are right. I was under the impressions that multiple call were required for the procedure. If that's the case, then caching the Hessian is not very important for performance.
Hi, we can close for the hessian but for the parameter gradient, an interface would still be welcome.
Does BifurcationKit use ForwadDiff to compute it internally or only finite difference ?
You are right. Let's leave it open for parameter gradiant.
Right now, it is FD everywhere but it is just a matter of adding a function dpF(::BifurcationProblem to the interface.