Would it be possible to get some code like below working? The first example with the function f is meant to show that this definition of hessian_vector_product can work. The second example shows that this fails with g, which uses AlgebraicMultigrid. If Hessian-vector products could be computed efficiently this way, it would be really useful.
using Test
import AlgebraicMultigrid
import ForwardDiff
import LinearAlgebra
import SparseArrays
import Zygote
hessian_vector_product(f, x, v) = ForwardDiff.jacobian(s->Zygote.gradient(f, x + s[1] * v)[1], [0.0])[:]
n = 4
A = randn(n, n)
hessian = A + A'
f(x) = LinearAlgebra.dot(x, A * x)
x = randn(n)
v = randn(n)
hvp1 = hessian_vector_product(f, x, v)
hvp2 = hessian * v
@test hvp1 ≈ hvp2#the hessian_vector_product plausibly works!
function g(x)
k = x[1:n + 1]
B = SparseArrays.spdiagm(0=>k[1:end - 1] + k[2:end], -1=>-k[2:end - 1], 1=>-k[2:end - 1])
ml = AlgebraicMultigrid.ruge_stuben(B)
return sum(AlgebraicMultigrid.solve(ml, x[N + 2:end]))
end
x = randn(2 * n + 1)
v = randn(2 * n + 1)
hessian_vector_product(g, x, v)#seems to fail during the coarse solve in AlgebraicMultigrid
Would it be possible to get some code like below working? The first example with the function
fis meant to show that this definition ofhessian_vector_productcan work. The second example shows that this fails withg, which uses AlgebraicMultigrid. If Hessian-vector products could be computed efficiently this way, it would be really useful.