Partial Derivatives of Multivariable Functions

import numpy as np

def partial_derivative(f, x: list, var_index: int,
                       epsilon: float = 1e-7) -> float:
    """
    Compute partial derivative df/dx_i at point x.
    f: function taking a list/array of inputs
    x: point at which to evaluate
    var_index: index of the variable to differentiate with respect to
    """
    x = list(x)
    x_plus = x.copy()
    x_minus = x.copy()
    x_plus[var_index] += epsilon
    x_minus[var_index] -= epsilon
    return (f(x_plus) - f(x_minus)) / (2 * epsilon)

def gradient(f, x: list, epsilon: float = 1e-7) -> list:
    """Compute the gradient vector [df/dx_0, df/dx_1, ..., df/dx_n]."""
    return [partial_derivative(f, x, i, epsilon) for i in range(len(x))]

def hessian(f, x: list, epsilon: float = 1e-5) -> list:
    """Compute the Hessian matrix of second partial derivatives."""
    n = len(x)
    H = [[0.0] * n for _ in range(n)]
    for i in range(n):
        for j in range(n):
            x_pp = list(x)
            x_pm = list(x)
            x_mp = list(x)
            x_mm = list(x)
            x_pp[i] += epsilon
            x_pp[j] += epsilon
            x_pm[i] += epsilon
            x_pm[j] -= epsilon
            x_mp[i] -= epsilon
            x_mp[j] += epsilon
            x_mm[i] -= epsilon
            x_mm[j] -= epsilon
            H[i][j] = (f(x_pp) - f(x_pm) - f(x_mp) + f(x_mm)) / (
                4 * epsilon ** 2)
    return H

def directional_derivative(f, x: list, direction: list,
                            epsilon: float = 1e-7) -> float:
    """Compute the directional derivative along a unit vector."""
    d = np.array(direction, dtype=float)
    d = d / np.linalg.norm(d)
    grad = np.array(gradient(f, x, epsilon))
    return float(np.dot(grad, d))