1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
|
# The SDP example of section 8.7 (Exploiting structure).
from cvxopt import blas, lapack, solvers, matrix, normal
def mcsdp(w):
"""
Returns solution x, z to
(primal) minimize sum(x)
subject to w + diag(x) >= 0
(dual) maximize -tr(w*z)
subject to diag(z) = 1
z >= 0.
"""
n = w.size[0]
def Fs(x, y, alpha = 1.0, beta = 0.0, trans = 'N'):
"""
y := alpha*(-diag(x)) + beta*y.
"""
if trans=='N':
# x is a vector; y is a matrix.
blas.scal(beta, y)
blas.axpy(x, y, alpha = -alpha, incy = n+1)
else:
# x is a matrix; y is a vector.
blas.scal(beta, y)
blas.axpy(x, y, alpha = -alpha, incx = n+1)
def cngrnc(r, x, alpha = 1.0):
"""
Congruence transformation
x := alpha * r'*x*r.
r and x are square matrices.
"""
# Scale diagonal of x by 1/2.
x[::n+1] *= 0.5
# a := tril(x)*r
a = +r
tx = matrix(x, (n,n))
blas.trmm(tx, a, side='L')
# x := alpha*(a*r' + r*a')
blas.syr2k(r, a, tx, trans = 'T', alpha = alpha)
x[:] = tx[:]
def Fkkt(W):
rti = W['rti'][0]
# t = rti*rti' as a nonsymmetric matrix.
t = matrix(0.0, (n,n))
blas.gemm(rti, rti, t, transB = 'T')
# Cholesky factorization of tsq = t.*t.
tsq = t**2
lapack.potrf(tsq)
def f(x, y, z):
"""
Solve
-diag(z) = bx
-diag(x) - inv(rti*rti') * z * inv(rti*rti') = bs
On entry, x and z contain bx and bs.
On exit, they contain the solution, with z scaled
(inv(rti)'*z*inv(rti) is returned instead of z).
We first solve
((rti*rti') .* (rti*rti')) * x = bx - diag(t*bs*t)
and take z = -rti' * (diag(x) + bs) * rti.
"""
# tbst := t * zs * t = t * bs * t
tbst = matrix(z, (n,n))
cngrnc(t, tbst)
# x := x - diag(tbst) = bx - diag(rti*rti' * bs * rti*rti')
x -= tbst[::n+1]
# x := (t.*t)^{-1} * x = (t.*t)^{-1} * (bx - diag(t*bs*t))
lapack.potrs(tsq, x)
# z := z + diag(x) = bs + diag(x)
z[::n+1] += x
# z := -rti' * z * rti = -rti' * (diag(x) + bs) * rti
cngrnc(rti, z, alpha = -1.0)
return f
c = matrix(1.0, (n,1))
# Initial feasible x: x = 1.0 - min(lambda(w)).
lmbda = matrix(0.0, (n,1))
lapack.syevx(+w, lmbda, range='I', il=1, iu=1)
x0 = matrix(-lmbda[0]+1.0, (n,1))
s0 = +w
s0[::n+1] += x0
# Initial feasible z is identity.
z0 = matrix(0.0, (n,n))
z0[::n+1] = 1.0
dims = {'l': 0, 'q': [], 's': [n]}
sol = solvers.conelp(c, Fs, w[:], dims, kktsolver = Fkkt,
primalstart = {'x': x0, 's': s0[:]}, dualstart = {'z': z0[:]})
return sol['x'], matrix(sol['z'], (n,n))
n = 100
w = normal(n,n)
x, z = mcsdp(w)
|
