Primjeri korištenja PICOSa
Najprije rješavamo problem koji dolazi iz teorije kvantnih informacija.
Definiramo pojam vjernosti dvaju pozitivno semidefinitnih operatora $P$, $Q$ kao $$ F(P,Q) = \max_U \left\lvert \mathop{\mathrm{tr}}(P^{1/2}UQ^{1/2})\right\rvert,$$ gdje $U$ ide po svim unitarnim matricama. $F(P,Q)$ se može zapisati kao rješenje sljedećeg problema $$ \begin{gather*} \max_{Z\in\mathbb{C}^{n\times n}} \frac{1}{2} (Z+ Z^*), \\ \begin{pmatrix} P & Z \\ Z^* & Q \end{pmatrix}\ge 0. \end{gather*} $$
import picos as pic
import cvxopt as cvx
#Kreiramo P,Q pozitivne semidefinitne matrice
P = cvx.matrix([ [1-1j , 2+2j , 1 ],
[3j , -2j , -1-1j],
[1+2j, -0.5+1j, 1.5 ]
])
P = P * P.H
Q = cvx.matrix([ [-1-2j , 2j , 1.5 ],
[1+2j ,-2j , 2.-3j ],
[1+2j ,-1+1j , 1+4j ]
])
Q = Q * Q.H
n=P.size[0]
#Ubacimo P,Q u PICOS
P = pic.new_param('P',P)
Q = pic.new_param('Q',Q)
F = pic.Problem()
Z = F.add_variable('Z',(n,n),'complex')
F.set_objective('max','I'|0.5*(Z+Z.H)) #ili ('I' | Z.real) # | označava skalarni produkt, I je jedinična matrica
F.add_constraint(((P & Z) // (Z.H & Q))>>0 )
print(F)
F.solve(verbose = 0, solver = 'cvxopt')
print('vjernost: F(P,Q) = {0:.4f}'.format(F.obj_value()))
print('optimalna matrica Z:')
print(Z)
--------------------- optimization problem (Complex SDP): 18 variables, 0 affine constraints, 21 vars in 1 SD cones Z_IM : (3, 3), continuous Z_RE : (3, 3), continuous maximize trace(0.5·(Z + Zᴴ)) such that [P, Z; Zᴴ, Q] ≽ 0 --------------------- vjernost: F(P,Q) = 37.4742 optimalna matrica Z: [ 1.51e+01+j2.21e+00 -7.17e+00-j1.22e+00 2.52e+00+j6.87e-01] [-4.88e+00+j4.06e+00 1.00e+01-j1.57e-01 8.33e+00+j1.13e+01] [-4.33e-01+j2.98e-01 3.84e+00-j3.28e+00 1.24e+01-j2.05e+00]
Provjera:
import numpy as np
PP = np.matrix(P.value)
QQ = np.matrix(Q.value)
S,U = np.linalg.eig(PP)
sqP = U * np.diag([s**0.5 for s in S]) * U.H #P^{1/2}
S,U = np.linalg.eig(QQ)
sqQ = U * np.diag([s**0.5 for s in S]) * U.H #Q^{1/2}
fidelity = sum(np.linalg.svd(sqP * sqQ)[1]) #tr(P^{1/2} * Q^{1/2})
print('vjernost izračunata direktno: F(P,Q) = {0:.4f}'.format(fidelity))
vjernost izračunata direktno: F(P,Q) = 37.4742
U sljedećem primjeru nas zanima sljedeći minimizacijski problem $$ \begin{gather*} \min x_0^T X x_0,\\ X \ge 0,\\ A^T X + X A + R \le 0 . \end{gather*} $$
A = cvx.matrix([[-1,0,0],[0,-2,0],[0,0,-3]])
R = cvx.matrix( [ [1, 2, 1],
[3, -2, -1],
[1, 2, 3]])
R = R * R.T
x0 = cvx.matrix([[0,1,-1]])
A = pic.new_param('A',A)
R = pic.new_param('R',R)
n = 3
J = pic.Problem()
X = J.add_variable('X',(n,n),vtype='symmetric')
J.add_constraint(X>>0)
J.add_constraint(A.T*X + X*A + R<<0)
<3×3 LMI Constraint: Aᵀ·X + X·A + R ≼ 0>
J.set_objective('min', X*x0|x0)
print(J)
--------------------- optimization problem (SDP): 6 variables, 0 affine constraints, 12 vars in 2 SD cones X : (3, 3), symmetric minimize ⟨X·[3×1], [3×1]⟩ such that X ≽ 0 Aᵀ·X + X·A + R ≼ 0 ---------------------
J.solve()
=================================== PICOS 1.2.0.post13 =================================== Building a CVXOPT problem instance. Solving the problem via CVXOPT. ----------------------------------- Python Convex Optimization Solver via internal CONELP solver ----------------------------------- pcost dcost gap pres dres k/t 0: 7.6116e-01 3.4761e+01 3e+01 3e-01 5e+00 1e+00 1: 1.1162e+00 3.2958e+00 9e-01 2e-02 3e-01 2e-01 2: 8.3494e-01 1.0025e+00 6e-02 2e-03 2e-02 1e-02 3: 8.3440e-01 8.3942e-01 2e-03 5e-05 7e-04 2e-04 4: 8.3334e-01 8.3340e-01 2e-05 5e-07 7e-06 2e-06 5: 8.3333e-01 8.3333e-01 2e-07 5e-09 7e-08 2e-08 6: 8.3333e-01 8.3333e-01 2e-09 6e-11 9e-10 2e-10 Optimal solution found. ------------[ CVXOPT ]------------- Solution is optimal after 1.9e-02s. =============[ PICOS ]=============
{'cvxopt_sol': {'x': <6x1 matrix, tc='d'>, 'y': <0x1 matrix, tc='d'>, 's': <18x1 matrix, tc='d'>, 'z': <18x1 matrix, tc='d'>, 'status': 'optimal', 'gap': 2.4636188208613878e-09, 'relative gap': 2.956342555506986e-09, 'primal objective': 0.8333333349432771, 'dual objective': 0.8333333416563087, 'primal infeasibility': 6.420409562789395e-11, 'dual infeasibility': 9.034576311960793e-10, 'primal slack': 1.1930745639503172e-09, 'dual slack': 1.0096377193079867e-11, 'residual as primal infeasibility certificate': None, 'residual as dual infeasibility certificate': None, 'iterations': 6}, 'status': 'optimal', 'time': 0.018549680709838867, 'primals': {'X': [5.696959895488621, -0.9428109538225613, 3.000000028747104, 0.3535596486324303, 2.8284271645586343, 1.8333333624994705]}, 'duals': [<3x3 matrix, tc='d'>, <3x3 matrix, tc='d'>], 'obj': 0.833333338299793}
Provjera:
import numpy as np
import scipy.linalg as la
AA = np.matrix(A.value)
RR = np.matrix(R.value)
x00 = np.array(x0)
ly = la.solve_continuous_lyapunov(AA,-RR)
x00.T @ ly @ x00
array([[0.83333333]])