I guess you’re right theoretically but when I try it in practice I get vastly different runtimes. Maybe an implementation issue? Maybe I misunderstood something about R-op or L-op? This is what the comparison looks like (timing on CPU):
Pearlmutter:                       3.89 s ± 151  ms per loop (7 runs, 1 loop each)
Simple AD:                         1.21 s ± 64.5 ms per loop (7 runs, 1 loop each)
Finite difference approximation:   327 ms ± 11   ms per loop (7 runs, 1 loop each)

The first two return the same (correct) result and the last one returns a pretty good approximation. Here is more or less what the Pearlmutter implementation looks like:
def hvp(model, criterion, inputs, outputs, v, *args, **kwargs):

    # Zero the gradients
    model.zero_grad()

    # Forward-pass
    prediction = model(inputs)
    L = criterion(prediction, outputs)

    # Compute gradient
    y = torch.autograd.grad(L, model.parameters(), create_graph = True, retain_graph = True)
    # Apply R-operator on gradient
    vp = rop(parameters_to_vector(y), model.parameters(), v)

    return vp

def rop(y, x, v):

    # Adapted from: https://gist.github.com/apaszke/c7257ac04cb8debb82221764f6d117ad

    w = torch.zeros(y.size(), requires_grad = True)
    g = torch.autograd.grad(y, x, grad_outputs = w, create_graph = True, allow_unused = True)
    r = torch.autograd.grad(parameters_to_vector(g), w, grad_outputs = v, create_graph = False, allow_unused = True)
    return parameters_to_vector(r)