Python:CVXPY 解算器错误

Question

目的:我试图在 python 中使用 cvxpy 来最大化 dual_func，但是我得到下面的 SolverError，我相信它可能是 Variable 具有不同的维度，但似乎无法弄清楚这个问题。我试过使用其他求解器，如 ECOS，但无济于事。我是凸优化的新手，如果有人能帮助我指出代码中错误的正确方向，那就太好了。

问题:代码在 cvxpy 版本 0.4 中有效，但在最新的 cvxpy 版本中无效，给我错误:

SolverError: Either candidate conic solvers (['CVXOPT']) do not support the
cones output by the problem (SOC, ExpCone, PSD), or there are not enough 
constraints in the problem.

dual_func是凸的，是DCP，所以我仍然不确定为什么不能解决。

附加问题:是否可以在 cvxpy 中使用严格的不等式约束？文档指出这是不可能的，因为它在现实世界中没有意义，但是，我真的不明白为什么会这样。

from cvxpy import *
alpha_k = Variable()
mu_k_1 = Variable(shape=(int(len(h_undl_list))))
mu_k_2 = Variable(shape=(int(len(h_undl_list))))
covar_k_1 = Variable(shape=(int(len(h_undl_list)), int(len(h_undl_list))), PSD=True)
covar_k_2 = Variable(shape=(int(len(h_undl_list)), int(len(h_undl_list))), PSD=True)

# fitting the first 1000 returns with normal GMM
from sklearn import mixture
clf = mixture.GaussianMixture(n_components=2, max_iter=1000, tol=1e-7).fit(ret_df_period)
clf = mixture.GaussianMixture(n_components=2, max_iter=500, tol=1e-7, weights_init=[0.2, 0.8]).fit(ret_df_period)

print('weights')
print(clf.weights_)

# sorting the values to the normal regime or the stressed regime
if clf.weights_[0] > clf.weights_[1]:
    alpha_k.value = clf.weights_[0]
    mu_k_1.value = clf.means_[0]
    mu_k_2.value = clf.means_[1]
    covar_k_1.value = clf.covariances_[0]
    covar_k_2.value = clf.covariances_[1]
else:
    alpha_k.value = clf.weights_[1]
    mu_k_1.value = clf.means_[1]
    mu_k_2.value = clf.means_[0]
    covar_k_1.value = clf.covariances_[1]
    covar_k_2.value = clf.covariances_[0]

iter_t = 3
N = ret_df_period.shape[0] #time step, e.g. return at t
K = 2

for t_iter in range(iter_t):
    q_k = [alpha_k.value, 1 - alpha_k.value]
    mu_sk_list = [np.array(mu_k_1.value).flatten(), np.array(mu_k_2.value).flatten()]
    covar_sk_list = [covar_k_1.value, covar_k_2.value]
    q_posterior_list = []
    q_posterior_map = {}
    
    
    # E-Step
    import math
    # marginal probability, lambda_k * multivariate_normal.pdf, returns a matrix with K rows, N columns
    mar_prob_list = np.zeros((K,N))
    for k_idx_temp in range(0, K):
        for j_temp in range(0, N):
            mar_prob_list[k_idx_temp, j_temp] = (q_k[k_idx_temp] * multivariate_normal.pdf(ret_df_period.iloc[j_temp].values, mean = mu_sk_list[k_idx_temp], cov = covar_sk_list[k_idx_temp]))
    
    # lambda_k_n, q_n, q_posterior_distribution, returns a matrix with K rows, N columns
    for k_idx in range(0, K):
        for j in range(0, N):
            joint_prob_k = q_k[k_idx] * multivariate_normal.pdf(ret_df_period.iloc[j], mean = mu_sk_list[k_idx], cov = covar_sk_list[k_idx])
            q_posterior = joint_prob_k / np.sum(mar_prob_list[:, j])
            q_posterior_list.append(q_posterior)
        
        q_posterior_map[k_idx] = q_posterior_list
        q_posterior_list = []
    
    
    #####a priori after E-Step, i.e. M-Step
    # Empirical probability, returns the two lambda_sk values
    alpha_sk_list = []
    for k_idx in range(0, K):
        alpha_sk = 0
        for j in range(0, N):
            alpha_sk += q_posterior_map[k_idx][j]
        alpha_sk_list.append(alpha_sk / N)
    alpha_sk_list[1] = 1 - alpha_sk_list[0]
        
    
    # Mean, returns two sets of means
    mu_sk_list = []
    weighted_xs_map = {}
    for k_idx in range(0, K):
        weighted_xs_map[k_idx] = np.empty([int(len(h_undl_list)), N])
    
    for k_idx in range(0, K):
        mu_sk = 0
        for j in range(0, N):
            weighted_x = (q_posterior_map[k_idx][j]) * ret_df_period.iloc[j].values
            mu_sk += weighted_x
            weighted_xs_map[k_idx][:, j] = weighted_x / alpha_sk_list[k_idx] # should divide twice of alpha_sk here? Come back and double check later!!!
        mu_sk_list.append(mu_sk / N / alpha_sk_list[k_idx])
    
    # Cov, returns two cov
    covar_sk_list = []
    for k_idx in range(0, K):
        sum_covar = 0
        for j in range(0, N):
            w0_temp = (ret_df_period.iloc[j] - mu_sk_list[k_idx]).values
            w0_covar = np.outer(w0_temp, w0_temp)
            weighted_covar = (q_posterior_map[k_idx][j]) * w0_covar
            sum_covar += weighted_covar
        sum_covar = sum_covar / (alpha_sk_list[k_idx] * N)
        covar_sk_list.append(sum_covar)


    #####
    #information params
    n_sk_list = []
    for k_idx in range(0, K - 1):
        n_sk = np.log(alpha_sk_list[k_idx] / (1 - (np.sum(alpha_sk_list[: -1]))))
        n_sk_list.append(n_sk)
    
    m_sk_list = []
    for k_idx in range(0, K):
        m_sk = np.dot(np.linalg.inv(covar_sk_list[k_idx]), mu_sk_list[k_idx])
        m_sk_list.append(m_sk)
    
    S_sk_list = []
    for k_idx in range(0, K):
        S_sk = np.linalg.inv(covar_sk_list[k_idx])
        S_sk_list.append(S_sk)
    
    
    alpha_sk_1 = Parameter(nonneg=True)
    alpha_sk_2 = Parameter(nonneg=True)
    n_sk_pr = Parameter()
    dim = Parameter(nonneg=True)
    m_sk_1 = Parameter(shape=(int(len(h_undl_list)),))
    m_sk_2 = Parameter(shape=(int(len(h_undl_list)),))
    alpha_sk_1.value = alpha_sk_list[0]
    alpha_sk_2.value = alpha_sk_list[1]
    n_sk_pr.value = n_sk_list[0]
    m_sk_1.value = m_sk_list[0]
    m_sk_2.value = m_sk_list[1]
    dim.value = int(len(h_undl_list))
    covar_k_tilde = Variable(shape=(int(len(h_undl_list)), int(len(h_undl_list))), PSD=True)
    
    
    dual_func = entr(alpha_k) + entr(1 - alpha_k) + \
    alpha_sk_1*0.5*log_det(covar_k_1) + alpha_sk_1*dim*0.5*log(2*np.e*np.pi) + alpha_sk_2*0.5*log_det(covar_k_1 + covar_k_tilde) + alpha_sk_2*dim*0.5*log(2*np.e*np.pi) + alpha_k*n_sk_pr + \
    alpha_sk_1*(mu_k_1.T)*m_sk_1 - alpha_sk_1*0.5*trace(covar_k_1*S_sk_list[0]) - alpha_sk_1*0.5*quad_form(mu_k_1, S_sk_list[0]) + \
    alpha_sk_2*(mu_k_1.T)*m_sk_2 - alpha_sk_2*0.5*trace((covar_k_1 + covar_k_tilde)*S_sk_list[1]) - alpha_sk_2*0.5*quad_form(mu_k_2, S_sk_list[1])


    prob = Problem(Maximize(dual_func), [alpha_k == 0.96, mu_k_1 >= 0, mu_k_2 - mu_k_1 <= 0, mu_k_2 >= - mu_k_1, mu_k_2 <= 0])
    
    print(prob.solve(solver=CVXOPT, verbose=True, kktsolver='robust', max_iters=1000))
    print(t_iter)
    covar_k_2.value = covar_k_1.value + covar_k_tilde.value
    print(is_pos_def(covar_k_1.value))
    print(is_pos_def(covar_k_tilde.value))
    print(is_pos_def(covar_k_2.value))
    
    covar_k_1_temp = np.array(covar_k_1.value)
    covar_k_2_temp = np.array(covar_k_2.value)
    covar_k_tilde_temp = np.array(covar_k_tilde.value)
    mu_k_1_temp = np.array(mu_k_1.value).flatten()
    mu_k_2_temp = np.array(mu_k_2.value).flatten()
    D_1 = np.diag(np.sqrt(np.diag(covar_k_1_temp)))
    correl_k_1 = np.dot(np.dot(np.linalg.inv(D_1), covar_k_1_temp), np.linalg.inv(D_1))
    D2 = np.diag(np.sqrt(np.diag(covar_k_2_temp)))
    correl_k_2 = np.dot(np.dot(np.linalg.inv(D_2), covar_k_2_temp), np.linalg.inv(D_2))

错误:

---------------------------------------------------------------------------
SolverError                               Traceback (most recent call last)
<ipython-input-28-1e1e9c7b0d73> in <module>
    109     prob = Problem(Maximize(dual_func), [alpha_k == 0.96, mu_k_1 >= 0, mu_k_2 - mu_k_1 <= 0, mu_k_2 >= - mu_k_1, mu_k_2 <= 0])
    110 
--> 111     print(prob.solve(solver=CVXOPT, verbose=True, kktsolver='robust', max_iters=1000))
    112     print(t_iter)
    113     covar_k_2.value = covar_k_1.value + covar_k_tilde.value

~\Anaconda3\lib\site-packages\cvxpy\problems\problem.py in solve(self, *args, **kwargs)
    287         else:
    288             solve_func = Problem._solve
--> 289         return solve_func(self, *args, **kwargs)
    290 
    291     @classmethod

~\Anaconda3\lib\site-packages\cvxpy\problems\problem.py in _solve(self, solver, warm_start, verbose, parallel, gp, qcp, **kwargs)
    565                     solver, warm_start, verbose, **kwargs)
    566 
--> 567         self._construct_chains(solver=solver, gp=gp)
    568         data, solving_inverse_data = self._solving_chain.apply(
    569             self._intermediate_problem)

~\Anaconda3\lib\site-packages\cvxpy\problems\problem.py in _construct_chains(self, solver, gp)
    508 
    509             except Exception as e:
--> 510                 raise e
    511 
    512     def _solve(self,

~\Anaconda3\lib\site-packages\cvxpy\problems\problem.py in _construct_chains(self, solver, gp)
    503                 self._solving_chain = \
    504                     construct_solving_chain(self._intermediate_problem,
--> 505                                             candidate_solvers)
    506 
    507                 self._cached_chain_key = chain_key

~\Anaconda3\lib\site-packages\cvxpy\reductions\solvers\solving_chain.py in construct_solving_chain(problem, candidates)
     93                       "enough constraints in the problem." % (
     94                           candidates['conic_solvers'],
---> 95                           ", ".join([cone.__name__ for cone in cones])))
     96 
     97 

SolverError: Either candidate conic solvers (['CVXOPT']) do not support the cones output by the problem (SOC, ExpCone, PSD), or there are not enough constraints in the problem.

代码引用自:https://www.maths.ox.ac.uk/system/files/attachments/TT18_dissertation_Vu_0.pdf

Answer 1

您的目标包含 PSD 矩阵、log_det 和 entr，因此您需要一个支持半定锥和指数锥的求解器。据我所知，只有 MOSEK 9 和 SCS 可以在 CVXPY 中为您完成这两项工作。

严格的不等式没有实际意义，因为浮点运算的精度和未达到的问题。考虑最小化 x 服从 x>0。具有严格不等式的问题往往在极限中有最佳点，这意味着您无论如何都在合并非严格闭包。