Python中带约束的对数似然函数参数优化求解

问题分析 #

你的问题是最大化带$\theta_n, \phi_m \in (0,1)$约束的对数似然函数，参数规模达3070维（3000+70）。Scipy默认优化器在这种规模下效率偏低，加上初始猜测全为1违反边界约束、未使用解析梯度等问题，导致收敛极慢。

优化方案 #

1. 优化现有Scipy实现 #

• 修正初始猜测：将初始值设为(0,1)区间内的随机值，避免一开始触达边界，例如np.random.uniform(0.1, 0.9, N+M)。
• 手动实现解析梯度：Scipy默认数值梯度对高维参数计算量极大，手动推导并实现梯度能大幅提速：
  负目标函数的梯度：
  • 对$\theta_n$的偏导：$- \sum_m \frac{x_{nm}\phi_m}{1 + x_{nm}\theta_n\phi_m}$
  • 对$\phi_m$的偏导：$- \sum_n \frac{x_{nm}\theta_n}{1 + x_{nm}\theta_n\phi_m}$
• 更换优化器：使用L-BFGS-B，它专门适配带边界约束的问题，无需额外定义重复的不等式约束（bounds已覆盖边界要求）。

修改后的Scipy代码：

import numpy as np
from scipy.optimize import minimize

def objective_function(params, X):
    N, M = X.shape
    theta = params[:N]
    phi = params[N:]
    term = X * theta.reshape(-1, 1) * phi.reshape(1, -1)
    return -np.sum(np.log(1 + term))

def gradient(params, X):
    N, M = X.shape
    theta = params[:N]
    phi = params[N:]
    term = X * theta.reshape(-1, 1) * phi.reshape(1, -1)
    denom = 1 + term
    # 计算theta的梯度
    grad_theta = -np.sum(X * phi.reshape(1, -1) / denom, axis=1)
    # 计算phi的梯度
    grad_phi = -np.sum(X * theta.reshape(-1, 1) / denom, axis=0)
    return np.concatenate([grad_theta, grad_phi])

N, M = A.shape
# 修正初始猜测，落在(0,1)区间内
initial_guess = np.random.uniform(0.2, 0.8, N + M)
bounds = [(0, 1)] * (N + M)

# 使用L-BFGS-B优化器，传入解析梯度
result = minimize(objective_function, initial_guess, args=(A,), 
                  bounds=bounds, jac=gradient, method='L-BFGS-B')

2. 更适合大规模场景的替代库 #

• PyTorch/TensorFlow：利用自动微分框架处理高维参数优化，Adam、SGD等优化器在大规模场景下比Scipy更高效。通过参数投影确保约束满足：
  PyTorch示例代码：

  import torch
  import torch.optim as optim
  
  N, M = A.shape
  X = torch.tensor(A, dtype=torch.float32)
  # 初始化参数，限制在(0.1, 0.9)区间内
  theta = torch.nn.Parameter(torch.rand(N, dtype=torch.float32) * 0.8 + 0.1)
  phi = torch.nn.Parameter(torch.rand(M, dtype=torch.float32) * 0.8 + 0.1)
  
  optimizer = optim.Adam([theta, phi], lr=1e-3)
  epochs = 1000
  
  for epoch in range(epochs):
      optimizer.zero_grad()
      term = X * theta.unsqueeze(1) * phi.unsqueeze(0)
      loss = -torch.sum(torch.log(1 + term))
      loss.backward()
      optimizer.step()
      # 投影到(0,1)区间，确保约束满足
      with torch.no_grad():
          theta.clamp_(0, 1)
          phi.clamp_(0, 1)
      if epoch % 100 == 0:
          print(f"Epoch {epoch}, Loss: {loss.item()}")
  

• CVXPY：若验证目标函数为凸函数，可使用CVXPY调用ECOS、OSQP等高效凸求解器，自动处理约束与优化逻辑。

额外优化技巧 #

• 利用矩阵A的稀疏性：仅对非0元素计算对数项，减少不必要的计算量，大幅降低求和耗时。

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