This came up in a conversation with @seantalts . When solving an algebraic equation, we find y^* such that f(y^*, \theta) = 0 where \theta are our model parameters. We compute the jacobian matrix of y with respect to \theta using the implicit function theorem: J = - [J^y(y^*, \theta)]^{-1} J^\thet…

That’s the way to do it – build up the full Jacobian and then partition it into the y and \theta blocks. It’s not so much about the number of sweeps here but rather the shared computations. That said there is the potential for some addition improvements. If you don’t compute everything at once th…

[image] charlesm93: Any objections to this idea? Nope. And Jacobian calcs should get faster and less memory intensive when we get forward-mode working. [image] charlesm93: The cost should be roughly halved Just make sure to test that.

There’s a potential for piggyback if Newton’s iteration is used to solve for y^{\ast}, since J^y is required anyway. When that’s the case, you don’t want to always compute J^y and J^{\theta} at the same time because J^y gets updated much more frequently. So maybe design accordingly based on the scen…

@yizhang Newton solver is still on the list – I’ve just had a busy semester – , and Powell’s method also uses derivatives, so your comment applies there too. In the solver, J^y is computed at each iteration except at the solution, which is what we need for the sensitivities. The computation of J^y …

[image] charlesm93: @betanalpha Calculating JyvJ^y v would work. Here vv is the adjoint of y∗y^* with respect to the log density, and is a member of the vari class. We would need a version of Jacobian() function that lets you specify an initial cotangent vector. Is there a way of doing this with…

If only J^y y^i is needed for the Powell updates then you should definitely only compute that! Ok. I think the devil will be in the details here. If memory serves well, Eigen’s function for Powell’s method accepts J^y as an argument, but not a method to compute J^y y. But I need to break the code…

There are a couple examples for how the Jacobian-adjoint stuff @betanalpha mentioned pieces together over here: Adj_jac_apply , Adj_jac_apply Lemme know if you end up trying it and it’s not magically working. The templating makes the compiler errors pretty annoying to dig through.

Alright, let’s give it a try. [image] Issue: Improving differentiation for algebraic solver opened by charlesm93 on 2019-06-04 Description There are two ways to improve the way in which we compute the Jacobian matrix of an algebraic solver: compute the Jacobian... …

For the INLA optimization problem, I think we would gain much from computing Jv with the Jacobian-vector product formalism. Indeed we have J^y in analytical form. J^\theta on the other hand is computed in n sweeps, where n is the dimension of y. In practice y is high-dimensional so this computation …

Better computation of the Jacobian matrix for algebraic solver

Developers

charlesm93 June 4, 2019, 12:28pm 9

Alright, let’s give it a try.

Scalable Bayesian multilevel modeling

Topic		Replies	Views
Improved differentiation of the algebraic solver Developers	1	588	May 26, 2020
Algebraic Solver Differentiation Speedup Developers performance	23	1648	September 22, 2021
Algebra solver details? Developers features	9	2234	August 14, 2017
Adj_jac_apply Developers	37	2621	April 11, 2020
Precomputed Jacobians or vector-Jacobian products? Algorithms performance	7	1207	January 31, 2022

Better computation of the Jacobian matrix for algebraic solver

Related topics