计算科学 - CVX能否解决凸问题？ - 吾爱随笔录

CVX能否解决凸问题？

计算科学优化 matlab 约束优化凸优化简历

2021-12-22 02:38:18

我想解决凸优化如下：

\begin{aligned} min_{X_{1}, X_{2}} & - \frac{1}{N} \sum_{i = 1}^{N} \log det (I + H_{i}^{H} X_{2} H_{i}) - \log [1 + h^{H} (X_{1} + X_{2}) h] \\ + tr (A X_{1}) + tr (B X_{2}) \\ s.t & tr (X_{1} + X_{2}) \leq 1 \\ X_{1} \in S_{+}^{n}, X_{2} \in S_{+}^{n} \end{aligned}

$\begin{align} \underset{X_1,X_2}{\min} &\ -\frac{1}{N}\sum_{i=1}^N\log\det\left(I+H_i^HX_2H_i\right)-\log\left[1+h^H(X_1+X_2)h\right]\\ &\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+ \text{tr}(AX_1)+\text{tr}(BX_2)\\ \text{s.t} &\ \text{tr}(X_1+X_2) \leq 1 \\ &\ X_1 \in S^n_+,X_2 \in S^n_+ \end{align}$ 与，，，给出。和都是厄密矩阵。下标代表共轭转置。

H_{i}

$H_i$

h

$h$

A

$A$

B

$B$

A

$A$

B

$B$

H

$H$

我已经使用 CVX 来解决这个问题，但性能似乎并不好。因为在两个CVX码中使用了Sylvester的行列式恒等式，我们得到了完全不同的结果，这是不合理的。

如果 CVX 无法有效解决问题，SDPT3 或 SeDuMi 等其他工具箱能否解决问题？谢谢你。

matlab代码如下所示：

Sample = 10;
H = randn(4,2,Sample)+1i*randn(4,2,Sample);
h = randn(4,1)+1i*randn(4,1);
A = randn(4,4)+1i*randn(4,4); A = A*A';
B = randn(4,4)+1i*randn(4,4); B = B*B';

cvx_begin quiet
    variable X1(4,4) semidefinite hermitian
    variable X2(4,4) semidefinite hermitian
    t = -log(1+quad_form(h,X1+X2))+trace(A*X1)+trace(B*X2);
    for i = 1:Sample
        t =  t - log_det(eye(2) + H(:,:,i)'*X2 *H(:,:,i))/Sample; % Primal constraint
    end
    minimize(t)
    subject to
        trace(X1+X2)<=1;
cvx_end
X1
X2
t

cvx_begin quiet
    variable X1(4,4) semidefinite hermitian
    variable X2(4,4) semidefinite hermitian
    t = -log(1+quad_form(h,X1+X2))+trace(A*X1)+trace(B*X2);
    for i = 1:Sample
        t =  t - log_det(eye(4) + X2 *H(:,:,i)*H(:,:,i)')/Sample; % Sylvester's determinant identity is used
    end
    minimize(t)
    subject to
        trace(X1+X2)<=1;
cvx_end
X1
X2
t
```

1个回答

您应该使用建模语言，以便您的代码独立于底层求解器。

cvxpy是一个不错的选择。

当我在 cvxpy 中重写你的模型时：

#!/usr/bin/env python3
import cvxpy as cp
import numpy as np
from numpy.random import normal as randn

Sample = 10
H = randn(size=(4,2,Sample))+1j*randn(size=(4,2,Sample))
h = randn(size=(4,1))+1j*randn(size=(4,1))
A = randn(size=(4,4))+1j*randn(size=(4,4))
A = A*A.T;
B = randn(size=(4,4))+1j*randn(size=(4,4))
B = B*B.T;

X1 = cp.Variable((4,4), symmetric=True)
X2 = cp.Variable((4,4), symmetric=True)

t = -cp.log(cp.real(1+cp.quad_form(h,X1+X2)))+cp.trace(A*X1)+cp.trace(B*X2)
for i in range(Sample):
    t -= cp.log_det(np.eye(2) + H[:,:,i].T*X2*H[:,:,i])/Sample # Primal constraint

constraints = [cp.trace(X1+X2)<=1, X1>>0, X2>>0]

problem1 = cp.Problem(cp.Minimize(cp.real(t)), constraints)
optval1 = problem1.solve()
print("Optimval value", optval1)
print("X1.value", X1.value)
print("X2.value", X2.value)

t = -cp.log(cp.real(1+cp.quad_form(h,X1+X2)))+cp.trace(A*X1)+cp.trace(B*X2);
for i in range(Sample):
    t -= cp.log_det(np.eye(4) + X2*H[:,:,i]*H[:,:,i].T)/Sample #Sylvester's determinant identity is used

problem2 = cp.Problem(cp.Minimize(cp.real(t)), constraints)
optval2 = problem2.solve()
print("Optimval value", optval2)
print("X1.value", X1.value)
print("X2.value", X2.value)

我发现解决这两个问题大约需要 2 秒，并且两者同意小数点后几位（-2.1221305764766405 与 -2.122080730649998）。

如果您需要更好的性能，您可以让 cvxpy使用许多求解器中的任何一种。

如果您需要更好的性能，请在JuMP中重写模型，它提供了许多求解器的灵活性以及 Julia 的性能（在第一次运行代码后，您的模型将很快得到转换，让 JIT 有机会工作）。

其它你可能感兴趣的问题

上一篇检查 SOR 解决方案是否已针对 2d 矩阵收敛的最佳方法下一篇在 Python 中实现 Gelfand 的光谱半径公式 - 缺乏收敛性