基于scipy.minimize实现贝尔曼方程的递归目标函数优化问题

问题背景 #

正在用Python研究最优执行问题，不熟悉动态规划，尝试用递归最小化实现，但当前代码依赖v和q的初始猜测值，希望让v和q完全由整体最小化过程生成。尝试的递归贝尔曼方程实现导致内核崩溃，需要可行方案。

变量定义：

• v：各期交易股票的数组
• q：对冲投资组合数组，q[0]已知，其余元素递归定义为q[i+1] = q[i] + v[i]*dt
• prices：股票价格数组，终值根据期末股价取值为N或0
• V：常数，J：最大时间周期，i：迭代时的时间步

当前代码问题分析 #

1. 依赖初始猜测：现有实现中v_tree和q_tree的初始化依赖vb、qB等初始猜测数组，没有让优化过程完全自主生成变量。
2. 递归贝尔曼实现的错误：
   • 递归函数内部嵌套循环，导致无限递归或重复计算，引发内核崩溃
   • 终端条件判断逻辑错误，循环和递归的流程混乱，无法正确终止递归
   • 每次优化都触发递归调用，计算量呈指数级增长，超出内存和计算资源

可行实现方案 #

方案1：反向递推（动态规划标准实现） #

放弃递归，采用反向递推的动态规划方式，从最后一期开始向前计算，存储每一期的最优值函数theta和最优交易速度v，避免重复计算和递归崩溃。

import numpy as np
from scipy.optimize import minimize

def compute_optimal_strategy(initial_q, prices, J, dt, V, gamma, eta, phi, rho, mu, sigma, epsilon, payoff_wcost, T, K, N):
    # 初始化存储数组，无需初始猜测
    v_opt = np.zeros(J)
    theta = np.zeros(J)
    q_opt = np.zeros(J+1)
    q_opt[0] = initial_q  # 已知初始仓位

    # 处理最后一期（终端条件）
    i = J-1
    def terminal_obj(v):
        current_q = q_opt[i]
        next_q = current_q + v * dt
        return np.exp(J-(i+1)*dt)*gamma * (
            current_q * prices[i] * (np.exp(dt)-1) +
            L(v/V, eta, phi) * V * dt -
            next_q * (mu*dt + sigma*np.sqrt(dt)*epsilon[i]) +
            payoff_wcost("physical", N, prices[i], K, next_q, T, (i+1)*dt, V)
        )
    
    # 优化终端期的v
    res = minimize(terminal_obj, x0=0, bounds=[(-rho*V, rho*V)], method="trust-constr")
    v_opt[i] = res.x[0]
    theta[i] = res.fun  # 存储终端期的最优值
    q_opt[i+1] = q_opt[i] + v_opt[i] * dt

    # 反向递推前面的期数
    for i in reversed(range(J-1)):
        def step_obj(v):
            current_q = q_opt[i]
            next_q = current_q + v * dt
            # 利用下一期已经计算好的theta[i+1]
            return np.exp(J-(i+1)*dt)*gamma * (
                current_q * prices[i] * (np.exp(dt)-1) +
                L(v/V, eta, phi) * V * dt -
                next_q * (mu*dt + sigma*np.sqrt(dt)*epsilon[i]) +
                theta[i+1]
            )
        
        # 优化当前期的v
        res = minimize(step_obj, x0=0, bounds=[(-rho*V, rho*V)], method="trust-constr")
        v_opt[i] = res.x[0]
        theta[i] = res.fun
        q_opt[i+1] = q_opt[i] + v_opt[i] * dt
    
    return v_opt, q_opt, theta

# 调用示例
initial_q = 0.5 * N
v_opt, q_opt, theta = compute_optimal_strategy(
    initial_q=initial_q,
    prices=prices,
    J=J,
    dt=dt,
    V=V,
    gamma=gamma,
    eta=eta,
    phi=phi,
    rho=rho,
    mu=mu,
    sigma=sigma,
    epsilon=epsilon,
    payoff_wcost=payoff_wcost,
    T=T,
    K=K,
    N=N
)

# 可视化仓位变化
import matplotlib.pyplot as plt
time_points = np.linspace(0, J*dt, J+1)
plt.plot(time_points, q_opt/N)
plt.xlabel("Time")
plt.ylabel("Normalized Position")
plt.show()

方案2：带记忆化的递归（避免重复计算） #

如果坚持用递归，必须加入记忆化缓存，存储已经计算过的(q, i)状态的最优值，避免重复递归计算导致的内核崩溃。

import numpy as np
from scipy.optimize import minimize
from functools import lru_cache

# 注意：q需要离散化才能用lru_cache，因为浮点数作为缓存键会有精度问题
def discretize_q(q, N, num_bins=1000):
    return round(q / N * num_bins)

def bellman_recursion(initial_q, prices, J, dt, V, gamma, eta, phi, rho, mu, sigma, epsilon, payoff_wcost, T, K, N):
    num_bins = 1000  # 离散化仓位的精度

    @lru_cache(maxsize=None)
    def value_function(discrete_q, i):
        q = discrete_q * N / num_bins
        if i == J:
            # 终端条件：返回payoff
            return payoff_wcost("physical", N, prices[-1], K, q, T, J*dt, V)
        
        def objective(v):
            next_q = q + v * dt
            discrete_next_q = discretize_q(next_q, N, num_bins)
            # 递归调用下一期的价值函数
            return np.exp(J-(i+1)*dt)*gamma * (
                q * prices[i] * (np.exp(dt)-1) +
                L(v/V, eta, phi) * V * dt -
                next_q * (mu*dt + sigma*np.sqrt(dt)*epsilon[i]) +
                value_function(discrete_next_q, i+1)
            )
        
        # 优化当前期的v
        res = minimize(objective, x0=0, bounds=[(-rho*V, rho*V)], method="trust-constr")
        return res.fun, res.x[0]
    
    # 从初始状态开始计算
    v_opt = []
    current_discrete_q = discretize_q(initial_q, N, num_bins)
    for i in range(J):
        val, v = value_function(current_discrete_q, i)
        v_opt.append(v)
        current_q = current_discrete_q * N / num_bins
        current_discrete_q = discretize_q(current_q + v*dt, N, num_bins)
    
    # 计算q的序列
    q_opt = np.zeros(J+1)
    q_opt[0] = initial_q
    for i in range(J):
        q_opt[i+1] = q_opt[i] + v_opt[i] * dt
    
    return v_opt, q_opt

# 调用示例
initial_q = 0.5 * N
v_opt, q_opt = bellman_recursion(
    initial_q=initial_q,
    prices=prices,
    J=J,
    dt=dt,
    V=V,
    gamma=gamma,
    eta=eta,
    phi=phi,
    rho=rho,
    mu=mu,
    sigma=sigma,
    epsilon=epsilon,
    payoff_wcost=payoff_wcost,
    T=T,
    K=K,
    N=N
)

# 可视化
time_points = np.linspace(0, J*dt, J+1)
plt.plot(time_points, q_opt/N)
plt.show()

关键改进点 #

1. 反向递推方案：
   • 从终端期开始计算，逐步向前推导，每一期的计算依赖下一期的已计算结果，避免递归的重复计算
   • 完全自主生成v_opt和q_opt，无需初始猜测数组
2. 记忆化递归方案：
   • 用lru_cache缓存已计算的状态，避免指数级计算量
   • 对仓位q进行离散化，解决浮点数作为缓存键的精度问题
3. 优化逻辑统一：
   • 每一期的目标函数逻辑一致，终端期直接用payoff，前期用下一期的价值函数，避免代码冗余

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