Hessian calculation fails using tf.jacobian

Python: 3.8
TFQ: 0.4.0
Hi,

I am trying to get some Hessians from the one of my parameterized quantum circuits. This trainstep works as intended:

def build_train_step(circuit: cirq.Circuit, symbols: List,
                     paulisum: List[cirq.PauliSum], learning_rate: float) -> \
        Tuple[Any, tf.Variable, tfq.layers.Expectation]:
    model_params = tf.Variable(tf.random.uniform([1, len(symbols)]) * 2,
                               constraint=lambda x: tf.clip_by_value(x, 0, 4))
    expectation_layer = tfq.layers.Expectation()
    optimizer = tf.keras.optimizers.Adam(learning_rate=learning_rate)

    @tf.function
    def train_step():
        with tf.GradientTape() as tape:
            expectation_batch = expectation_layer(circuit,
                                                  symbol_names=symbols,
                                                  symbol_values=model_params,
                                                  operators=paulisum)
            energy = tf.reduce_sum(expectation_batch)
        gradients = tape.gradient(energy, model_params)
        optimizer.apply_gradients(zip([gradients], [model_params]))
        return energy

    return train_step, model_params, expectation_layer

Following this example in the TensorFlow 2 docs I was hoping I could get the Hessian with the following code:

def build_train_step_hessians(circuit: cirq.Circuit, symbols: List,
                     paulisum: List[cirq.PauliSum], learning_rate: float) -> \
        Tuple[Any, tf.Variable, tfq.layers.Expectation]:
    model_params = tf.Variable(tf.random.uniform([1, len(symbols)]) * 2,
                               constraint=lambda x: tf.clip_by_value(x, 0, 4))
    expectation_layer = tfq.layers.Expectation()
    optimizer = tf.keras.optimizers.Adam(learning_rate=learning_rate)

    @tf.function
    def train_step():
        with tf.GradientTape() as t2:
            with tf.GradientTape() as t1:
                expectation_batch = expectation_layer(circuit,
                                                      symbol_names=symbols,
                                                      symbol_values=model_params,
                                                      operators=paulisum)
                energy = tf.reduce_sum(expectation_batch)
            gradients = t1.gradient(energy, model_params)
        hess = t2.jacobian(gradients, model_params)
        optimizer.apply_gradients(zip([gradients], [model_params]))

        return energy, hess

    return train_step, model_params, expectation_layer

But this throws the error:

...
    LookupError: No gradient defined for operation 'TfqAdjointGradient' (op type: TfqAdjointGradient)

From which I conclude that calculating gradients of gradients is not supported yet ( I tried the other differentiators as well). Am I out of luck here? Or is there a hack I can use to get the Hessians from the circuit? Thanks! If you need an example where I use this train step I can throw one together.

P.S.
I am in the process of rewriting all my research code to TFQ and so far everything has worked like a charm. No more super slow graph building times and worrying about how to extract stuff the graph with my own TF1 simulator. And the adjoint differentiator in TFQ is amazing as well; I ran a VQE optimization with like 500 parameters the other day without any issues. Great stuff!

Start with the nested GradientTape example and the tfq.layers.Expectation call in the issue, then reproduce the LookupError for TfqAdjointGradient under TensorFlow's advanced autodiff pattern. Trace the differentiator entry point involved and define done as obtaining the Hessian without the missing-gradient failure, with coverage for the reported parameterized-circuit case.