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
|
import numpy as np
import numpy.linalg as la
from scipy.interpolate import PchipInterpolator
def decr_fun(a, b):
return np.sqrt(np.square(a).sum() + np.square(b).sum())
def KLfit(P, q, b):
"""
solves for x, min KL(x, q)
s.t. Px = b
x>=0
sum(x) = 1
"""
alpha = 0.4
beta = 0.8
if(len(P.shape)==1):
P = P[None]
P = np.vstack([P, np.ones(P.shape[1])])
b = np.hstack([b, 1])
#init x and nu
x = np.ones(P.shape[1])/P.shape[1]
nu = np.ones(P.shape[0])
decr = np.Inf
eps = 1e-12
niter = 1
while decr > eps and niter < 500:
rdual = 1 + np.log(x) - np.log(q) + P.T.dot(nu)
rprimal = P.dot(x) - b
S = -P.dot((x*P).T)
Dnu_pd = la.solve(S, (P*x).dot(rdual) - rprimal)
Dx_pd = -x*(P.T.dot(Dnu_pd)+rdual)
#backtracking search
phi = 1
newx = x + phi * Dx_pd
newnu = nu + phi * Dnu_pd
while sum(newx<0) > 0:
phi *= beta
newx = x + phi * Dx_pd
newnu = nu + phi * Dnu_pd
newrdual = 1 + np.log(newx) - np.log(q) + P.T.dot(newnu)
newrprimal = P.dot(newx) - b
while decr_fun(newrdual, newrprimal) > (1 - alpha * phi) * decr_fun(rdual, rprimal):
phi *= beta
newx = x + phi * Dx_pd
newnu = nu + phi * Dnu_pd
newrdual = 1 + np.log(newx) - np.log(q) + P.T.dot(newnu)
newrprimal = P.dot(newx) - b
x = newx
nu = newnu
decr = decr_fun(newrdual, newrprimal)
niter += 1
return {"obj": sum(x*(np.log(x) - np.log(q))), "weight":x, "status": decr}
def interpweights(w, v1, v2):
""" Given L=(w, v1), computes neww such that newL=(neww, v2) = L in distribution"""
cumw = np.cumsum(w)
neww = PchipInterpolator(v1, cumw)(v2, der = 1)
interpweights = new/np.sum(neww)
return interpweights
def interpvalues(w, v, neww):
""" Given a distribution D=(w,v), compute new weights
such that Dnew=(neww, newv) equals D in distribution """
cumw = np.cumsum(w)
cdf = PchipInterpolator(v, cumw)
eps = 1e-4
neww = np.zeros(neww.size)
cumneww = np.cumsum(neww)
mid = 0
for i in range(neww.size):
iter = 0
hi = cdf(1)
lo = mid
if hi < cumneww[i]:
neww[i] = hi
continue
if cdf(lo) > cumneww[i]:
neww[i] = lo
continue
mid = (lo+hi)/2
iter = 0
while( abs(cdf(mid)-cumneww[i]) > eps ):
if cdf(mid) > cumneww[i]:
hi = mid
else:
lo = mid
mid = (lo + hi)/2
newv[i] = mid
return newv
if __name__=="__main__":
#write some small test
P = np.array([[5, 4, 3, 2, 1],[3, 3, 4, 3, 3]])
w = np.ones(5)/5
b = np.array([2.5, 3.3])
test = KLfit(P, w, b)
|
