My DFP algorithm doesn't converge at all and generates a lot of curvature errors. I'll send you my implementation. Then, I send you the implementation that I want to follow: # Define DFP Optimization
def DFP(f, func_grad, x0, tol=1e-6, max_iter=10000, time_limit=100, stop="time", H0=None):
    """
    Performs the DFP (Davidon-Fletcher-Powell) optimization algorithm with Armijo line search.

    Parameters:
    - f: objective function.
    - func_grad: gradient function.
    - x0: initial point as a NumPy array.
    - tol: tolerance for the gradient norm.
    - max_iter: maximum number of iterations.
    - time_limit: time limit for optimization (in seconds).
    - stop: stopping criterion ("time" or "max_iter").
    - H0: initial Hessian approximation matrix. Defaults to identity matrix.

    Returns:
    - x_k: optimized variable.
    - k: number of iterations performed.
    - total_armijo: total Armijo iterations.
    - f_current: final function value.
    - f_list: list of function values per iteration.
    - error_ratios: list of error ratios per iteration.
    - total_time: total time taken for optimization.
    """
    k = 0
    x_k = copy.deepcopy(x0)
    f_current = f(x_k)
    f_list = [f_current]  # Initialize with the initial function value
    f_prev = f_current
    error_ratios = []
    
    n = len(x0)
    H_k = np.eye(n) if H0 is None else H0
    start_time = time.time()
    
    total_armijo = 0
    grad_k = func_grad(x_k)
    
    while True:
        # Check stopping conditions
        elapsed_time = time.time() - start_time
        if stop == "time" and elapsed_time >= time_limit:
            break
        if stop == "max_iter" and k >= max_iter:
            break
        if np.linalg.norm(grad_k) <= tol:
            break
        
        # Compute search direction
        d_k = -H_k @ grad_k
        
        # Perform Armijo line search
        step, armijo_iter = armijo_line_search(f, func_grad, x_k, d_k)
        total_armijo += armijo_iter
        
        # Update position
        x_next = x_k + step * d_k
        p_k = x_next - x_k
        grad_next = func_grad(x_next)
        q_k = grad_next - grad_k
        
        # Compute denominator to avoid division by zero
        p_q = np.dot(p_k, q_k)
        if p_q == 0:
            print(f"Violation of curvature condition at iteration {k}.")
            break
        
        # Update Hessian approximation using DFP formula
        H_k = H_k + np.outer(p_k, p_k) / p_q - (H_k @ np.outer(q_k, p_k) + np.outer(p_k, q_k) @ H_k) / np.dot(q_k, H_k @ p_k)
        
        # Update variables for next iteration
        x_k = x_next
        grad_k = grad_next
        k += 1
        
        # Compute new function value
        f_current = f(x_k)
        f_list.append(f_current)
        
        # Calculate error ratio
        if k > 1:
            if abs(f_prev) > 1e-12:
                error_ratio = abs(f_current - f_prev) / abs(f_prev)
            else:
                error_ratio = np.nan
            error_ratios.append(error_ratio)
        else:
            error_ratios.append(np.nan)
        
        f_prev = f_current
    
    total_time = time.time() - start_time
    return x_k, k, total_armijo, f_current, f_list, error_ratios, total_time