SymPy中derive_by_array结果变量替换与自定义多元函数求导问题

Hey there! Let’s work through your problem with building a custom function to calculate gradients and Hessians— I’ve tangled with this exact SymPy quirk before, so let’s break this down step by step.

First: Solving the derive_by_array Substitution Headache #

The main gotcha with derive_by_array is that it returns SymPy Array objects, not plain lists. To substitute numerical values correctly, you need to pass a dictionary that maps each symbol to its coordinate value to the .subs() method. Trying to substitute one variable at a time (like .subs(x, 1).subs(y, 2)) works but gets messy fast for multiple variables.

Full Robust Custom Function Implementation #

Here’s a tested, reusable function that handles any multivariable function, calculates both gradient and Hessian, and properly substitutes your coordinate values:

import sympy as sp

def compute_gradient_hessian(func, symbols, coords):
    # Quick sanity check: make sure symbols and coordinates match in count
    if len(symbols) != len(coords):
        raise ValueError("Number of symbols must match the number of coordinate values")
    
    # Create a clean substitution dictionary: {symbol: value, ...}
    subs_map = dict(zip(symbols, coords))
    
    # Calculate gradient (first-order partial derivatives)
    gradient = sp.derive_by_array(func, symbols)
    
    # Calculate Hessian (second-order partials: gradient of the gradient)
    hessian = sp.derive_by_array(gradient, symbols)
    
    # Substitute the coordinate values into both results
    grad_evaluated = gradient.subs(subs_map)
    hess_evaluated = hessian.subs(subs_map)
    
    # Optional: convert to numerical values (use sp.N() or .evalf())
    grad_numeric = sp.N(grad_evaluated)
    hess_numeric = sp.N(hess_evaluated)
    
    return grad_evaluated, hess_evaluated, grad_numeric, hess_numeric

Let’s Test This with an Example #

Let’s use a sample multivariable function to prove this works. Say we have f(x,y,z) = x² + 3yz + sin(xz) and want to evaluate at (1, 2, 3):

# Define our symbols and function
x, y, z = sp.symbols('x y z')
f = x**2 + 3*y*z + sp.sin(x*z)

# Coordinates to evaluate at (order matches the symbols list!)
eval_coords = (1, 2, 3)

# Run our custom function
grad_sym, hess_sym, grad_num, hess_num = compute_gradient_hessian(f, [x,y,z], eval_coords)

# Check the results
print("Symbolic Gradient:")
display(grad_sym)
print("\nNumerical Gradient:")
display(grad_num)
print("\nSymbolic Hessian:")
display(hess_sym)
print("\nNumerical Hessian:")
display(hess_num)

Common Fixes for Your Original Code #

If your initial custom function wasn’t working, it’s almost certainly one of these issues:

• Mismatched order: The order of symbols in derive_by_array determines the order of derivatives— your coordinate values must follow the exact same order.
• Bad substitution mapping: You probably weren’t using a dictionary to map symbols to values. The dict(zip(symbols, coords)) trick ensures every symbol gets paired with the right coordinate.
• Treating arrays like lists: SymPy Arrays have their own .subs() method— if you tried to loop through a list and substitute, you’d run into errors.
• Sanity checks missing: Forgetting to verify that the number of symbols matches the number of coordinates leads to silent errors or wrong results.

This implementation fixes all those gaps and works for any number of variables— perfect for your calculus practice!

内容的提问来源于stack exchange，提问作者Josmoor98