使用CVXPY实现最大方差展开(MVU)时遇不可行问题的技术求助

我尝试复现Weinberger和Saul（2004）论文中的最大方差展开(MVU)算法，用它学习数据分布的低维表示。去掉秩约束后的MVU通用问题如下：

$$\max_G \text{tr}(G)$$
$$\text{s.t. } (e_k - e_l)^T (X - G) (e_k - e_l) = 0, \forall (k,l) \in \text{邻接对}$$
$$G \succeq 0$$
其中，$G$是低维嵌入的内积矩阵，$X$是原始数据的内积矩阵，$e_k$是仅第k位为1的指示向量；原论文的零中心化约束被我暂时忽略。

我基于400张76×101像素的连续旋转茶壶图像数据集做实验，步骤如下：

1. 导入数据与库 #

import scipy.io
import cvxpy as cvx
import numpy as np
from sklearn.neighbors import NearestNeighbors

data = scipy.io.loadmat('teapots.mat')
data = data["Input"][0][0][0]

# 转换为(400 x 23028)的样本-特征矩阵
data = data.T.astype(np.float64)

2. 构建k近邻邻接矩阵 #

n_points = data.shape[0]
# 生成4近邻的邻接矩阵
nn = NearestNeighbors(n_neighbors=4).fit(data)
nn = nn.kneighbors_graph(data).todense()
nn = np.array(nn)

3. CVXPY求解MVU嵌入 #

# 原始数据的内积矩阵
X = cvx.Constant(data.dot(data.T))
# 低维嵌入的内积矩阵，PSD约束保证对称且半正定
G = cvx.Variable((n_points, n_points), PSD=True)
G.value = np.zeros((n_points, n_points))

# 目标：最大化嵌入的方差（等价于trace(G)）
objective = cvx.Maximize(cvx.trace(G))
constraints = []

# 添加邻接对的距离保持约束
for i in range(n_points):
    for j in range(n_points):
        if nn[i, j] == 1:
            constraints.append(
                (X[i, i] - 2 * X[i, j] + X[j, j]) -
                (G[i, i] - 2 * G[i, j] + G[j, j]) == 0
            )

problem = cvx.Problem(objective, constraints)
problem.solve(verbose=True, max_iters=10000000)
    
print(problem.status)

运行后CVXPY返回问题状态为infeasible；移除约束则返回unbounded（符合最大化凸函数的预期），且已验证邻接图是连通的（k=4近邻满足），求解错误原因。

问题排查与修正方案 #

1. 重复约束导致的数值负担 #

kneighbors_graph生成的是无向邻接矩阵，每个邻接对(i,j)和(j,i)会被重复添加约束。虽然理论上不影响可行性，但会大幅增加求解器计算量，甚至因数值精度问题触发不可行判定。

修正代码：只处理上三角的邻接关系，避免重复约束：

# 替换原约束循环
for i in range(n_points):
    for j in range(i+1, n_points):  # 仅处理i<j的邻接对
        if nn[i, j] == 1:
            constraints.append(
                (X[i, i] - 2 * X[i, j] + X[j, j]) -
                (G[i, i] - 2 * G[i, j] + G[j, j]) == 0
            )

2. 高维数据内积的数值精度问题 #

原始图像数据维度达23028维，直接计算data.dot(data.T)会产生极大数值，易引发求解器的数值溢出或精度误差，导致误判不可行。

修正代码：先对数据做标准化处理：

# 数据标准化：零均值+单位方差
data = (data - data.mean(axis=0)) / data.std(axis=0)
# 再计算内积矩阵
X = cvx.Constant(data.dot(data.T))

3. 缺失零中心化约束的隐含影响 #

原论文的零中心化约束$(1^T G 1)/n = 0$是保证嵌入合理性的关键，缺失该约束会让求解器的自由度过多，可能引发数值层面的不可行。

修正代码：添加零中心化约束：

# 添加零中心化约束，1为全1向量
ones = np.ones((n_points, 1))
constraints.append(cvx.trace(G @ ones @ ones.T) == 0)

4. 初始化与求解器参数优化 #

初始值设为全零矩阵，求解器可能难以找到可行方向；过大的迭代次数也会导致不必要的数值波动。

修正代码：用PCA低秩近似初始化G，并调整迭代次数：

from sklearn.decomposition import PCA
# 用PCA的2维嵌入初始化G
pca = PCA(n_components=2)
embedding = pca.fit_transform(data)
G_init = embedding.dot(embedding.T)

G = cvx.Variable((n_points, n_points), PSD=True)
G.value = G_init

# 求解时设置合理的迭代次数
problem.solve(verbose=True, max_iters=10000)

5. 约束形式的向量化优化 #

双层循环的约束构建效率低，且易引入数值误差，改用向量化方式提升稳定性：

# 生成上三角邻接掩码，避免重复约束
mask = (nn == 1) & (np.triu(np.ones(nn.shape), k=1) == 1)
i, j = np.where(mask)

# 向量化构建约束
for u, v in zip(i, j):
    constraints.append(
        cvx.norm(G[u,:] - G[v,:], 2)**2 == cvx.norm(X[u,:] - X[v,:], 2)**2
    )

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