计算科学 - 关于 CVXPY 结果的问题 - 吾爱随笔录

我要优化功能

min_{X \in S_{+}^{n}} tr (C^{T} X) + tr (X^{- 1}),

$\min_{X \in \mathbb{S}^{n}_{+}} \mbox{tr} \left( C^T X \right) + \mbox{tr} \left( X^{-1} \right),$

其中我优化了等效问题

min tr (C^{T} X) + tr (Z)

$\min \mbox{tr}\left(C^T X\right) + \mbox{tr}(Z)$

s.t. \begin{matrix} (\begin{matrix} Z & I \\ I & X \end{matrix}) \end{matrix} ⪰ 0

$\text{s.t.} \begin{matrix} \begin{pmatrix} Z & I\\ I & X \\ \end{pmatrix} \end{matrix} \succeq 0$

使其遵循 cvxpy 的 DCP 形式。此外，请注意该函数具有梯度。因此，对于某些对称正定矩阵，我将定义为。由于问题是凸的，因此我们在。通过这种方式，我可以检查优化的性能。然而，当我将代码放在 CVXPY 中时，算法表明最佳点是另一个点。下面是代码和输出。 $C - X^{-2}$ $C$ $X_0^{-2}$ $X_0$ $0$ $X = X_0$

import cvxpy as cp
import numpy as np
import random
def symm(A):
  #symmetric operator
  return 0.5 *(A + A.T)

np.random.seed(1)
n = 50
CN = 10
data_choice = 'sparse'
solver_choice = 'CG'

D = 1000 * np.diag(np.logspace(-np.log(CN), 0, n))
Q, R = np.linalg.qr(np.random.normal(size = (n,n)))
A = Q @ D @ Q
A = (A.T + A)/2
C= np.linalg.inv(A) @ np.linalg.inv(A)


X = cp.Variable((n, n) , symmetric=True)
Z = cp.Variable((n, n) , symmetric=True)
V = cp.bmat([[Z, np.eye(n)], [np.eye(n), X]])

constraints = [V >> 0]
obj = cp.Minimize(cp.trace(C.T @ X) + cp.trace(Z))
prob = cp.Problem(obj,constraints)

history = prob.solve(verbose = True)

===============================================================================
                                     CVXPY                                     
                                    v1.1.18                                    
===============================================================================
(CVXPY) Feb 14 12:43:40 AM: Your problem has 5000 variables, 1 constraints, and 0 parameters.
(CVXPY) Feb 14 12:43:40 AM: It is compliant with the following grammars: DCP, DQCP
(CVXPY) Feb 14 12:43:40 AM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)
(CVXPY) Feb 14 12:43:40 AM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.
-------------------------------------------------------------------------------
                                  Compilation                                  
-------------------------------------------------------------------------------
(CVXPY) Feb 14 12:43:40 AM: Compiling problem (target solver=SCS).
(CVXPY) Feb 14 12:43:40 AM: Reduction chain: Dcp2Cone -> CvxAttr2Constr -> ConeMatrixStuffing -> SCS
(CVXPY) Feb 14 12:43:40 AM: Applying reduction Dcp2Cone
(CVXPY) Feb 14 12:43:40 AM: Applying reduction CvxAttr2Constr
(CVXPY) Feb 14 12:43:40 AM: Applying reduction ConeMatrixStuffing
(CVXPY) Feb 14 12:43:40 AM: Applying reduction SCS
(CVXPY) Feb 14 12:43:40 AM: Finished problem compilation (took 3.301e-02 seconds).
-------------------------------------------------------------------------------
                                Numerical solver                               
-------------------------------------------------------------------------------
(CVXPY) Feb 14 12:43:40 AM: Invoking solver SCS  to obtain a solution.
------------------------------------------------------------------
           SCS v3.1.0 - Splitting Conic Solver
    (c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem:  variables n: 2550, constraints m: 5050
cones:    s: psd vars: 5050, ssize: 1
settings: eps_abs: 1.0e-05, eps_rel: 1.0e-05, eps_infeas: 1.0e-07
      alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
      max_iters: 100000, normalize: 1, warm_start: 0
      acceleration_lookback: 10, acceleration_interval: 10
lin-sys:  sparse-direct
      nnz(A): 2550, nnz(P): 0
------------------------------------------------------------------
 iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)
------------------------------------------------------------------
     0| 1.46e+01  1.00e+00  7.55e+02 -3.51e+02  1.00e-01  3.16e-03 
   250| 1.83e-03  1.19e-03  3.89e-02  3.03e+00  1.00e-01  3.89e-01 
   500| 3.93e-03  1.43e-04  3.54e-03  3.21e+00  4.62e-03  7.70e-01 
   750| 2.25e-03  4.97e-05  2.02e-03  3.00e+00  4.62e-03  1.15e+00 
  1000| 2.02e-03  3.04e-05  1.82e-03  2.99e+00  4.62e-03  1.53e+00 
  1250| 8.64e-04  1.27e-05  7.78e-04  2.90e+00  4.62e-03  1.90e+00 
  1500| 8.21e-04  1.03e-05  7.39e-04  2.90e+00  4.62e-03  2.27e+00 
  1750| 7.89e-04  9.81e-06  7.11e-04  2.90e+00  4.62e-03  2.63e+00 
  2000| 7.63e-04  9.40e-06  6.87e-04  2.90e+00  4.62e-03  3.00e+00 
  2250| 7.42e-04  9.04e-06  6.68e-04  2.89e+00  4.62e-03  3.37e+00 
  2500| 2.22e-04  1.02e-05  2.00e-04  2.86e+00  4.62e-03  3.73e+00 
  2750| 2.44e-04  4.48e-06  2.20e-04  2.86e+00  4.62e-03  4.10e+00 
  3000| 2.34e-04  3.80e-06  2.11e-04  2.86e+00  4.62e-03  4.50e+00 
  3250| 2.30e-04  3.30e-06  2.07e-04  2.86e+00  4.62e-03  4.89e+00 
  3500| 5.20e-03  1.44e-05  1.36e-03  2.84e+00  4.62e-03  5.28e+00 
  3750| 2.24e-04  2.94e-06  2.02e-04  2.86e+00  4.62e-03  5.65e+00 
  4000| 2.21e-04  2.91e-06  1.99e-04  2.86e+00  4.62e-03  6.03e+00 
  4250| 2.19e-04  2.88e-06  1.97e-04  2.86e+00  4.62e-03  6.41e+00 
  4500| 2.16e-04  2.86e-06  1.95e-04  2.86e+00  4.62e-03  6.79e+00 
  4750| 2.14e-04  2.84e-06  1.93e-04  2.86e+00  4.62e-03  7.16e+00 
  5000| 2.12e-04  2.82e-06  1.91e-04  2.86e+00  4.62e-03  7.53e+00 
  5250| 2.10e-04  2.80e-06  1.89e-04  2.86e+00  4.62e-03  7.90e+00 
  5500| 7.40e-05  2.96e-06  6.67e-05  2.85e+00  4.62e-03  8.28e+00 
  5750| 5.86e-05  1.93e-06  5.28e-05  2.85e+00  4.62e-03  8.66e+00 
  6000| 5.93e-05  1.60e-06  5.34e-05  2.85e+00  4.62e-03  9.02e+00 
  6250| 6.10e-05  1.38e-06  5.50e-05  2.85e+00  4.62e-03  9.39e+00 
  6500| 1.43e-04  2.24e-06  7.92e-05  2.85e+00  4.62e-03  9.76e+00 
  6750| 6.32e-05  1.18e-06  5.69e-05  2.85e+00  4.62e-03  1.01e+01 
  7000| 6.37e-05  1.11e-06  5.73e-05  2.85e+00  4.62e-03  1.05e+01 
  7250| 6.38e-05  1.11e-06  5.75e-05  2.85e+00  4.62e-03  1.09e+01 
  7500| 6.38e-05  1.11e-06  5.75e-05  2.85e+00  4.62e-03  1.12e+01 
  7750| 6.37e-05  1.10e-06  5.73e-05  2.85e+00  4.62e-03  1.16e+01 
  8000| 6.34e-05  1.10e-06  5.71e-05  2.85e+00  4.62e-03  1.20e+01 
  8250| 6.31e-05  1.09e-06  5.68e-05  2.85e+00  4.62e-03  1.24e+01 
  8500| 6.27e-05  1.09e-06  5.64e-05  2.85e+00  4.62e-03  1.27e+01 
  8750| 6.22e-05  1.08e-06  5.60e-05  2.85e+00  4.62e-03  1.31e+01 
  9000| 6.18e-05  1.07e-06  5.57e-05  2.85e+00  4.62e-03  1.35e+01 
  9250| 4.37e-03  3.10e-06  7.69e-04  2.84e+00  4.62e-03  1.39e+01 
  9500| 6.09e-05  1.06e-06  5.48e-05  2.85e+00  4.62e-03  1.42e+01 
  9750| 6.04e-05  1.06e-06  5.44e-05  2.85e+00  4.62e-03  1.46e+01 
 10000| 6.00e-05  1.05e-06  5.40e-05  2.85e+00  4.62e-03  1.50e+01 
 10250| 5.95e-05  1.04e-06  5.36e-05  2.85e+00  4.62e-03  1.54e+01 
 10500| 5.91e-05  1.04e-06  5.32e-05  2.85e+00  4.62e-03  1.57e+01 
 10750| 5.87e-05  1.03e-06  5.28e-05  2.85e+00  4.62e-03  1.61e+01 
 11000| 5.82e-05  1.02e-06  5.24e-05  2.85e+00  4.62e-03  1.65e+01 
 11250| 5.78e-05  1.02e-06  5.21e-05  2.85e+00  4.62e-03  1.69e+01 
 11500| 5.74e-05  1.01e-06  5.17e-05  2.85e+00  4.62e-03  1.72e+01 
 11750| 5.70e-05  1.01e-06  5.14e-05  2.85e+00  4.62e-03  1.76e+01 
 12000| 3.48e-03  1.43e-06  6.26e-04  2.84e+00  4.62e-03  1.79e+01 
 12250| 5.63e-05  9.97e-07  5.07e-05  2.85e+00  4.62e-03  1.83e+01 
 12500| 5.60e-05  9.91e-07  5.04e-05  2.85e+00  4.62e-03  1.87e+01 
 12750| 5.56e-05  9.86e-07  5.01e-05  2.85e+00  4.62e-03  1.90e+01 
 13000| 5.53e-05  9.81e-07  4.98e-05  2.85e+00  4.62e-03  1.94e+01 
 13250| 5.50e-05  9.76e-07  4.95e-05  2.85e+00  4.62e-03  1.98e+01 
 13500| 5.47e-05  9.71e-07  4.92e-05  2.85e+00  4.62e-03  2.02e+01 
 13750| 5.44e-05  9.67e-07  4.90e-05  2.85e+00  4.62e-03  2.06e+01 
 13900| 1.01e-05  1.20e-06  9.11e-06  2.84e+00  4.62e-03  2.08e+01 
------------------------------------------------------------------
status:  solved
timings: total: 2.08e+01s = setup: 1.19e-02s + solve: 2.08e+01s
     lin-sys: 9.26e-01s, cones: 1.94e+01s, accel: 1.21e-01s
------------------------------------------------------------------
objective = 2.844050
------------------------------------------------------------------
-------------------------------------------------------------------------------
                                    Summary                                    
-------------------------------------------------------------------------------
(CVXPY) Feb 14 12:44:00 AM: Problem status: optimal
(CVXPY) Feb 14 12:44:00 AM: Optimal value: 2.844e+00
(CVXPY) Feb 14 12:44:00 AM: Compilation took 3.301e-02 seconds
(CVXPY) Feb 14 12:44:00 AM: Solver (including time spent in interface) took 2.081e+01 seconds

我想知道是否有人可以指出问题的表述或代码本身是否有任何错误。我从 cvxpy 得到的最优值显然不是函数的最优值。

np.trace(C.T @ X.value) + np.trace(np.linalg.inv(X.value)) 
2.8435372028588484

np.trace(C.T @ A) + np.trace(np.linalg.inv(A)) 
-0.32275885803468984

所以我想知道我对问题的表述是否有任何错误。