Numpy与SymPy混合编程中的类型转换陷阱及解决方案

TypeError: loop of ufunc does not support argument 0 of type Float which has no callable sqrt method

import sympy as sp
x = sp.symbols('x')
expr = x + 0.5
val = expr.subs(x, 1) # val 是 sympy.Float 对象
print(type(val)) # <class 'sympy.core.numbers.Float'>

import sympy as sp
import numpy as np

def grad(f):
    X = f.free_symbols
    Y = [f.diff(xi) for xi in X]
    return [x_k for x_k in X], Y

def descente_pas_opti(f, X0, eps = 1e-6):
    Xk = X0
    fonction = sp.sympify(f)
    X, gradform = grad(fonction)
    r=sp.symbols('r')
    dform= np.array([-df_k for df_k in gradform])

    while True:
        # 关键修改：在创建dk数组时，显式指定dtype为np.float32
        dk = np.array(
            [df_k.subs(
                [(X[k],Xk[k]) for k in range(len(X))])
                    for df_k in dform]
            , dtype=np.float32) # 或 np.float64，取决于所需的精度

        # ... 后续计算 ...
        # 计算最优步长rho
        # 注意：这里也需要确保传递给np.dot的参数是NumPy兼容的类型
        # grad_at_Xk_plus_r_dk = [df_k.subs([(X[k], Xk[k] + r*dk[k]) for k in range(len(X))]) for df_k in gradform]
        # dot_product_expr = np.dot(grad_at_Xk_plus_r_dk, dk)
        # rho = sp.solve(dot_product_expr, r)[0]

        # 为了避免类似的类型问题，确保np.dot的输入也是SymPy表达式列表
        # 如果dk已经被转换为np.float32，那么rho的计算逻辑可能需要调整
        # 这里假设sp.solve能够处理SymPy表达式与NumPy数组的混合运算，但更稳妥的做法是保持一致性
        # 在SymPy求解前，将dk转换为SymPy的数值或保持其符号形式

        # 修正rho的计算逻辑，确保点积操作是在SymPy的上下文进行的，以避免类型冲突
        # grad_at_Xk_plus_r_dk 仍是SymPy表达式列表
        grad_at_Xk_plus_r_dk = [df_k.subs(
                                    [(X[k], Xk[k] + r*dk[k]) for k in range(len(X))] )
                                        for df_k in gradform]

        # 将dk转换为SymPy的数值列表，以便与grad_at_Xk_plus_r_dk进行点积
        # 或者，如果dk已经是np.float32，需要确保点积的结果是SymPy表达式
        # 更安全的做法是，在计算rho时，dk应该仍然是SymPy表达式形式，或者将其元素转换为SymPy数值

        # 考虑到dk现在是np.float32数组，这里需要将dk的元素转换为SymPy的Float
        # 以便与SymPy表达式进行点积，并由sp.solve处理
        dk_for_sympy = [sp.Float(val) for val in dk] # 将np.float32转换为sympy.Float

        dot_product_expr = sum(g * d for g, d in zip(grad_at_Xk_plus_r_dk, dk_for_sympy))
        rho = sp.solve(dot_product_expr, r)[0]

        # 更新Xk
        Xk = [Xk[0]+rho*dk[0], Xk[1]+rho*dk[1]]

        # 检查收敛条件
        if (np.linalg.norm(dk) < eps): break

    return Xk

# 示例参数
# descente_pas_opti('5*x**2 + 0.5*y**2 -3*(x + y)', [-2,-7])