基于scipy.optimize与cvxpy的凸函数优化问题求解

一、解决CVXPY的凹性错误 #

你的核心问题是目标函数的凸性判断偏差：直接优化Y时，目标函数T*log_det(Y) + trace(Y^{-1}@S)是凹函数+凸函数的组合，并非严格凸函数（当T较大时，凹项主导整体性质），而CVXPY要求最小化的目标必须是凸函数，因此触发报错。

正确的CVXPY实现：切换优化变量为Z=Y⁻¹ #

令Z=Y⁻¹（Z为对称正定矩阵），原目标函数可转换为严格凸函数：
原函数（忽略常数项）：
$$f(Y) = T\log\det(Y) + \text{trace}(Y^{-1}S)$$
代入Y=Z⁻¹，得：
$$f(Z) = -T\log\det(Z) + \text{trace}(ZS)$$
其中：

• -log_det(Z)是凸函数（CVXPY中log_det为凹函数，取负后变凸）
• trace(ZS)是线性函数（天然满足凸性要求）
  因此整体为凸函数，符合CVXPY的求解条件。

代码实现：

import cvxpy as cp
import numpy as np

def optimize_Y(T, term_data):
    m = term_data[0].shape[0]
    S = sum(term_data)  # 计算所有term_data的和S（正定矩阵）
    # 定义优化变量Z，约束为对称正定
    Z = cp.Variable((m, m), symmetric=True)
    constraints = [Z >> 0]
    # 构造凸目标函数
    objective = cp.Minimize(-T * cp.log_det(Z) + cp.trace(Z @ S))
    prob = cp.Problem(objective, constraints)
    # 求解可指定高效求解器，如SCS、ECOS
    prob.solve(solver=cp.SCS, verbose=True)
    # 转换回Y=Z⁻¹
    Y_opt = np.linalg.inv(Z.value)
    return Y_opt

二、提升Scipy.optimize在高维场景的性能 #

Scipy在高维下失效的主要原因是数值梯度误差大、参数化不合理，可通过以下两点优化：

1. 用Cholesky分解参数化Y，保证正定 #

将Y表示为Y = L@L.T，其中L是下三角矩阵（对角元为正），既保证Y正定，又减少参数冗余（参数个数为m(m+1)/2，与对称矩阵一致）。

2. 提供精确梯度，避免数值梯度 #

数值梯度在高维下计算慢、误差大，手动推导并实现精确梯度可大幅提升优化效率。

代码实现：

import numpy as np
from scipy.optimize import minimize

def unpack_lower_triangular(vec, m):
    """将向量转换为下三角矩阵，对角元用指数保证为正"""
    L = np.zeros((m, m))
    idx = 0
    for i in range(m):
        for j in range(i+1):
            if i == j:
                L[i, j] = np.exp(vec[idx])  # 指数确保对角元为正
            else:
                L[i, j] = vec[idx]
            idx += 1
    return L

def pack_lower_triangular(L):
    """将下三角矩阵转换为向量，对角元取对数"""
    m = L.shape[0]
    vec = []
    for i in range(m):
        for j in range(i+1):
            if i == j:
                vec.append(np.log(L[i, j]))
            else:
                vec.append(L[i, j])
    return np.array(vec)

def objective(vec, T, term_data):
    """目标函数计算"""
    m = int((np.sqrt(8*len(vec)+1)-1)/2)
    L = unpack_lower_triangular(vec, m)
    Y = L @ L.T
    S = sum(term_data)
    log_det = 2 * np.sum(np.log(np.diag(L)))  # log_det(Y)=2*sum(log(L的对角元))
    trace_term = np.trace(np.linalg.inv(Y) @ S)
    return T * log_det + trace_term

def gradient(vec, T, term_data):
    """精确梯度计算"""
    m = int((np.sqrt(8*len(vec)+1)-1)/2)
    L = unpack_lower_triangular(vec, m)
    Y = L @ L.T
    S = sum(term_data)
    Z = np.linalg.inv(Y)  # Z=Y⁻¹
    G = T * Z - Z @ S @ Z  # 目标函数对Y的梯度

    # 链式法则计算对L的梯度
    grad_L = G @ L
    grad_vec = []
    for i in range(m):
        for j in range(i+1):
            if i == j:
                # 对角元是指数变换，梯度需乘L[i,i]
                grad_vec.append(grad_L[i, i] * L[i, i])
            else:
                grad_vec.append(grad_L[i, j])
    return np.array(grad_vec)

# 使用示例
m = 3
T = 100
# 生成模拟的term_data（正定矩阵）
term_data = [np.random.randn(m,m) @ np.random.randn(m,m).T for _ in range(T)]
# 初始化参数（单位矩阵的Cholesky分解）
init_L = np.linalg.cholesky(np.eye(m))
init_vec = pack_lower_triangular(init_L)

# 调用L-BFGS-B优化器，传入精确梯度
result = minimize(
    objective, init_vec, args=(T, term_data),
    jac=gradient, method='L-BFGS-B',
    options={'maxiter': 1000, 'disp': True}
)
# 转换回最优Y
optimal_L = unpack_lower_triangular(result.x, m)
optimal_Y = optimal_L @ optimal_L.T

三、其他可行优化方法 #

1. 闭式解（无约束场景） #

如果你的问题是无约束的（即Y可以是任意对称正定矩阵），那么目标函数存在闭式解：
$$Y_{\text{opt}} = \frac{1}{T}\sum_{t=1}^T \text{term_data}[t]$$
这是多元正态协方差矩阵最大似然估计的经典结果，直接计算即可，无需优化器。

2. 流形优化（带结构约束场景） #

如果Y需要满足特殊结构约束（如低秩、对角、稀疏等），可以使用Manopt（黎曼流形优化库），它专门针对正定矩阵等流形上的优化问题，比通用优化器更高效。

3. 随机梯度下降（大数据集场景） #

当T极大时，可采用随机梯度下降（SGD）或Adam优化器，批量计算S的近似值，减少每次迭代的计算量，适合大规模数据集。

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