How to use the cvxopt.blas function in cvxopt
To help you get started, we’ve selected a few cvxopt examples, based on popular ways it is used in public projects.
Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
cvxopt / cvxopt / tests / test_custom_kkt.py View on Github 
# S = A*D^-1*A' + I
#
# where D = 2*D1*D2*(D1+D2)^-1, D1 = d[:n]**-2, D2 = d[n:]**-2.
d1, d2 = W['di'][:n]**2, W['di'][n:]**2
# ds is square root of diagonal of D
ds = math.sqrt(2.0) * div( mul( W['di'][:n], W['di'][n:]),
sqrt(d1+d2) )
d3 = div(d2 - d1, d1 + d2)
# Asc = A*diag(d)^-1/2
Asc = A * spdiag(ds**-1)
# S = I + A * D^-1 * A'
blas.syrk(Asc, S)
S[::m+1] += 1.0
lapack.potrf(S)
def g(x, y, z):
x[:n] = 0.5 * ( x[:n] - mul(d3, x[n:]) +
mul(d1, z[:n] + mul(d3, z[:n])) - mul(d2, z[n:] -
mul(d3, z[n:])) )
x[:n] = div( x[:n], ds)
# Solve
#
# S * v = 0.5 * A * D^-1 * ( bx[:n] -
# (D2-D1)*(D1+D2)^-1 * bx[n:] +
# D1 * ( I + (D2-D1)*(D1+D2)^-1 ) * bzl[:n] -
# D2 * ( I - (D2-D1)*(D1+D2)^-1 ) * bzl[n:] )# maximize -h'*z - b'*w
# subject to +/- c + G'*z + A'*w >= 0
# z >= 0
#
# with variables z (8), w (1).
cc = matrix(0.0, (9,1))
cc[:8], cc[8] = h, b
GG = spmatrix([], [], [], (n+8, 9))
GG[:n,:8] = -G.T
GG[:n,8] = -A.T
GG[n::n+9] = -1.0
hh = matrix(0.0, (n+8,n))
hh[:n,:] = matrix([i>=j for i in range(n) for j in range(n)],
(n,n), 'd') # upper triangular matrix of ones
l = [-blas.dot(cc, solvers.lp(cc, GG, hh[:,k])['x']) for k in range(n)]
u = [blas.dot(cc, solvers.lp(cc, GG, -hh[:,k])['x']) for k in range(n)]
def f(x,y): return x+2*[y]
def stepl(x): return functools.reduce(f, x[1:], [x[0]])
def stepr(x): return functools.reduce(f, x[:-1], []) + [x[-1]]
try: import pylab
except ImportError: pass
else:
pylab.figure(1, facecolor='w')
pylab.plot(stepl(a), stepr(p), '-')
pylab.title('Maximum entropy distribution (fig. 7.2)')
pylab.xlabel('a')
pylab.ylabel('p = Prob(X = a)')
pylab.figure(2, facecolor='w')def Fi(x, y, alpha = 1.0, beta = 0.0, trans = 'N'):
if trans == 'N':
# y := alpha * [P, -I; -P, -I] * x + beta*y
blas.gemv(P, x, u)
y[:m] = alpha * ( u - x[n:]) + beta*y[:m]
y[m:] = alpha * (-u - x[n:]) + beta*y[m:]
else:
# y := alpha * [P', -P'; -I, -I] * x + beta*y
blas.copy(x[:m] - x[m:], u)
blas.gemv(P, u, y, alpha = alpha, beta = beta, trans = 'T')
y[n:] = -alpha * (x[:m] + x[m:]) + beta*y[n:]raise TypeError("'b' must have length %d" %A.size[0])
if b is None and customy:
raise ValueEror("use of non vector type for y requires b")
if xnewcopy is None: xnewcopy = matrix
if xdot is None: xdot = blas.dot
if xaxpy is None: xaxpy = blas.axpy
if xscal is None: xscal = blas.scal
def xcopy(x, y):
xscal(0.0, y)
xaxpy(x, y)
if ynewcopy is None: ynewcopy = matrix
if ydot is None: ydot = blas.dot
if yaxpy is None: yaxpy = blas.axpy
if yscal is None: yscal = blas.scal
def ycopy(x, y):
yscal(0.0, y)
yaxpy(x, y)
# The problem is solved by applying cpl() to the epigraph form
#
# minimize t
# subject to f0(x) - t <= 0
# f1(x) <= 0
# ...
# fmnl(x) <= 0
# G*x <= h
# A*x = b.
#
# The epigraph form variable is stored as a list [x, t].# If not F['singular']:
# x := L^{-1} * P * (x + GGs'*z)
# = L^{-1} * P * (x + GG'*W^{-1}*W^{-T}*bz)
#
# If F['singular']:
# x := L^{-1} * P * (x + GGs'*z + A'*y))
# = L^{-1} * P * (x + GG'*W^{-1}*W^{-T}*bz + A'*y)
if mnl: base.gemv(F['Dfs'], z, x, trans = 'T', beta = 1.0)
base.gemv(F['Gs'], z, x, offsetx = mnl, trans = 'T',
beta = 1.0)
if F['singular']:
base.gemv(A, y, x, trans = 'T', beta = 1.0)
if type(F['S']) is matrix:
blas.trsv(F['S'], x)
else:
cholmod.solve(F['Sf'], x, sys = 7)
cholmod.solve(F['Sf'], x, sys = 4)
# y := K^{-1} * (Asc*x - y)
# = K^{-1} * (A * S^{-1} * (bx + GG'*W^{-1}*W^{-T}*bz) - by)
# (if not F['singular'])
# = K^{-1} * (A * S^{-1} * (bx + GG'*W^{-1}*W^{-T}*bz +
# A'*by) - by)
# (if F['singular']).
base.gemv(Asct, x, y, trans = 'T', beta = -1.0)
if type(F['K']) is matrix:
lapack.potrs(F['K'], y)
else:
cholmod.solve(Kf, y)pylab.ylabel('xcor[i]')
pylab.xlabel('i')
# A = D'*D is an n by n tridiagonal matrix with -1.0 on the
# upper/lower diagonal and 1, 2, 2, ..., 2, 2, 1 on the diagonal.
Ad = matrix([1.0] + (n-2)*[2.0] + [1.0])
As = matrix(-1.0, (n-1,1))
nopts = 50
deltas = -10.0 + 20.0/(nopts-1) * matrix(list(range(nopts)))
cost1, cost2 = [], []
for delta in deltas:
xr = +corr
lapack.ptsv(1.0 + 10**delta * Ad, 10**delta *As, xr)
cost1 += [blas.nrm2(xr - corr)]
cost2 += [blas.nrm2(xr[1:] - xr[:-1])]
# Find solutions with ||xhat - xcorr || roughly equal to 8.0, 3.1, 1.0.
mv1, k1 = min(zip([abs(c - 8.0) for c in cost1], range(nopts)))
xr1 = +corr
lapack.ptsv(1.0 + 10**deltas[k1] * Ad, 10**deltas[k1] *As, xr1)
mv2, k2 = min(zip([abs(c - 3.1) for c in cost1], range(nopts)))
xr2 = +corr
lapack.ptsv(1.0 + 10**deltas[k2] * Ad, 10**deltas[k2] *As, xr2)
mv3, k3 = min(zip([abs(c - 1.0) for c in cost1], range(nopts)))
xr3 = +corr
lapack.ptsv(1.0 + 10**deltas[k3] * Ad, 10**deltas[k3] *As, xr3)
if pylab_installed:
pylab.figure(2, facecolor='w')
pylab.plot(cost1, cost2, [blas.nrm2(corr)], [0], 'bo',def f(x, y, z):
# Solve for x[:n]:
#
# A*x[:n] = bx[:n] + P' * ( ((D1-D2)*(D1+D2)^{-1})*bx[n:]
# + (2*D1*D2*(D1+D2)^{-1}) * (bz[:m] - bz[m:]) ).
blas.copy(( mul( div(d1-d2, d1+d2), x[n:]) +
mul( 2*D, z[:m]-z[m:] ) ), u)
blas.gemv(P, u, x, beta = 1.0, trans = 'T')
lapack.potrs(A, x)
# x[n:] := (D1+D2)^{-1} * (bx[n:] - D1*bz[:m] - D2*bz[m:]
# + (D1-D2)*P*x[:n])
base.gemv(P, x, u)
x[n:] = div( x[n:] - mul(d1, z[:m]) - mul(d2, z[m:]) +
mul(d1-d2, u), d1+d2 )
# z[:m] := d1[:m] .* ( P*x[:n] - x[n:] - bz[:m])
# z[m:] := d2[m:] .* (-P*x[:n] - x[n:] - bz[m:])
z[:m] = mul(di[:m], u-x[n:]-z[:m])
z[m:] = mul(di[m:], -u-x[n:]-z[m:])for k in range(m):
edge = X[[k,k+1],:] + 0.1 * matrix([1., 0., 0., -1.], (2,2)) * \
(X[2*[k],:] - X[2*[k+1],:])
pylab.plot(edge[:,0], edge[:,1], 'k')
# 1000 points on the unit circle
nopts = 1000
angles = matrix( [ a*2.0*pi/nopts for a in range(nopts) ], (1,nopts) )
circle = matrix(0.0, (2,nopts))
circle[0,:], circle[1,:] = cos(angles), sin(angles)
# ellipse = L^-T * circle + xc where Hac = L*L'
lapack.potrf(Hac)
ellipse = +circle
blas.trsm(Hac, ellipse, transA='T')
ellipse += xac[:, nopts*[0]]
pylab.fill(ellipse[0,:].T, ellipse[1,:].T, facecolor = '#F0F0F0')
pylab.plot([xac[0]], [xac[1]], 'ko')
pylab.title('Analytic center (fig 8.7)')
pylab.axis('equal')
pylab.axis('off')
pylab.show()# Hi = Fi' * (diag(yi) - yi*yi') * Fi
# = Fisc' * Fisc
# where
# Fisc = diag(yi)^1/2 * (I - 1*yi') * Fi
# = diag(yi)^1/2 * (Fi - 1*gradfi')
Fsc[:K[i], :] = F[start:stop, :]
for k in range(start,stop):
blas.axpy(Df, Fsc, n=n, alpha=-1.0, incx=mnl+1,
incy=Fsc.size[0], offsetx=i, offsety=k-start)
blas.scal(math.sqrt(y[k]), Fsc, inc=Fsc.size[0],
offset=k-start)
# H += z[i]*Hi = z[i] * Fisc' * Fisc
blas.syrk(Fsc, H, trans='T', k=stop-start, alpha=z[i],
beta=1.0)
if z is None: return f, Df
else: return f, Df, H
