def multi_layer_second_directional_derivative(G, z, c, x, G_z, epsilon, w=None, Q=None):
    """
    Same as second_directional_derivative, but assumes G returns multiple outputs in a list
    """
    if w is None:
        _, G_to_x = G(z + x, c, return_bn=True, Q=Q)
        _, G_from_x = G(z - x, c, return_bn=True, Q=Q)
    else:
        _, G_to_x = G(z, c, w=w+x, return_bn=True, Q=Q)
        _, G_from_x = G(z, c, w=w-x, return_bn=True, Q=Q)

    eps_sqr = epsilon ** 2
    sdd = [(G2x - 2 * G_z_base + Gfx) / eps_sqr for G2x, G_z_base, Gfx in zip(G_to_x, G_z, G_from_x)]
    return sdd