问题描述
我正在尝试解决类似于下面描述的玩具示例的优化问题。正如评论中所指出的,我目前使用 scipy 的实现速度太慢,似乎没有收敛。我如何获得一个体面的解决方案?您可以使用 scipy、mystic 或任何您认为合适的库。
请注意,我不需要全局最小值,搜索可以在 loss(X) <= 1
后立即停止。现实世界的目标主要是用 sql 编写的,因此速度非常慢,所以我还希望优化在 loss
被评估约 200 次时终止。 (这不是硬性要求,您也可以在优化运行 5 分钟后终止。)
虽然这个问题类似于 Minimizing non-convex function with linear constraint and bound in mystic,但它绝对不是重复的。这两个问题甚至不是针对同一个目标。
import numpy as np
from scipy.optimize import differential_evolution,LinearConstraint
# Inhabitants in a rural city are voting for the name of a newborn. The city houses 40 dwarves
# and 10 humans in total,and there are 100 candidate names to Vote for.
dwarf_population,human_population = 40,10
population = dwarf_population + human_population
candidate_count = 100
# Each inhabitant has different number of Votes in their hand.
scores_per_citizen = np.random.randint(15,20,population)
# Let's say the original result looks like this. (Yes,one can cast a fraction of their Votes)
alpha = np.abs(np.random.normal(size=candidate_count))
scores = np.diag(scores_per_citizen).dot(np.random.dirichlet(alpha,population))
assert np.allclose(scores.sum(1),scores_per_citizen)
# Here is how the Votes are counted.
def count(scores_: np.ndarray) -> np.ndarray:
# Dwarves have a weird Tradition: the top 10 popular names chosen by dwarves will get all their Votes.
# (I guess this is what makes our objective non-convex.)
scores_by_dwarves = scores_[0:40,:]
score_per_candidate_by_dwarves_raw = scores_by_dwarves.sum(1)
top_10_popular_name_indices = np.argsort(-score_per_candidate_by_dwarves_raw)[:10]
score_per_candidate_by_dwarves = np.zeros(candidate_count)
score_per_candidate_by_dwarves[top_10_popular_name_indices] = score_per_candidate_by_dwarves_raw[
top_10_popular_name_indices]
score_per_candidate_by_dwarves = scores_by_dwarves.sum() \
* score_per_candidate_by_dwarves \
/ score_per_candidate_by_dwarves.sum()
assert np.allclose(score_per_candidate_by_dwarves.sum(),score_per_candidate_by_dwarves_raw.sum())
# Humans,on the other hand,just adds everyone's Votes together.
scores_by_humans = scores_[40:,:]
scores_per_candidate_by_humans = scores_by_humans.sum(0)
# The final result is the sum of the scores by dwarves and humans.
return score_per_candidate_by_dwarves + scores_per_candidate_by_humans
# So this is the legitimate result.
scores_per_candidate = count(scores)
# Now,you want to cheat a bit and make your proposal (assuming it's the first one) the most popular one.
target_scores_per_candidate = scores_per_candidate.copy()
argmax = scores_per_candidate.argmax()
target_scores_per_candidate[argmax] = scores_per_candidate[0]
target_scores_per_candidate[0] = scores_per_candidate[argmax]
assert np.allclose(scores_per_candidate.sum(),target_scores_per_candidate.sum())
# However,you cannot just manipulate the result,otherwise the auditors will find out!
# Instead,the raw scores must be manipulated such that your submission ranks top among others.
# You decide to solve for a multiplier to the original scores.
init_coef = np.ones_like(scores).reshape(-1)
# In other words,your goal is to find the minimum of the following objective.
def loss(coef_: np.ndarray) -> float:
scores_ = scores * coef_.reshape(scores.shape)
scores_per_candidate_ = count(scores_)
return ((scores_per_candidate_ - scores_per_candidate) ** 2).sum()
# This is a helper constant matrix. Ignore it for Now.
A = np.concatenate([np.tile(np.concatenate([np.repeat(1,candidate_count),np.repeat(0,population * candidate_count)]),population - 1),np.repeat(1,candidate_count)])
A = A.reshape((population,population * candidate_count))
# Meanwhile,some constraints must hold.
def constraints(coef_: np.ndarray):
# The total Votes of each citizen must not change.
coef_reshaped = coef_.reshape((population,candidate_count))
assert np.allclose((coef_reshaped * scores).sum(1),scores_per_citizen)
# Another way to express the same thing with matrices.
assert np.allclose(A.dot(np.diag(scores.reshape(-1))).dot(coef_),scores_per_citizen)
# Additionally,all scores must be positive.
assert np.all(coef_reshaped * scores >= 0)
# Let's express the constraint with a LinearConstraint.
score_match_quota = LinearConstraint(A.dot(np.diag(scores.reshape(-1))),scores_per_citizen,scores_per_citizen)
# Run optimization (Spoiler: this is extremely slow,and doesn't seem to converge)
rv = differential_evolution(loss,bounds=[(0,1000)] * init_coef.size,# the 1000 here is superficial and can be replaced by other large numbers
init=np.vstack((init_coef,init_coef,init_coef)),constraints=score_match_quota)
# Sanity check
constraints(rv.x)
解决方法
答案与您引用的问题几乎相同……但是,我必须进一步考虑一些限制来证明这一点。让我重写你的代码——只是使用一些更短的名字。
>>> import mystic as my
>>> import numpy as np
>>>
>>> # Inhabitants in a rural city are voting for the name of a newborn.
... # The city houses 40 dwarves and 10 humans in total,and there are
... # 100 candidate names to vote for.
... dwarves,humans = 40,10
>>> citizens = dwarves + humans
>>> names = 100
>>> # Each inhabitant has different number of votes in their hand.
... votes_per_citizen = np.random.randint(15,20,citizens)
>>>
>>> # Let's say the original result looks like this.
... # (Yes,one can cast a fraction of their votes)
... alpha = np.abs(np.random.normal(size=names))
>>> votes = np.diag(votes_per_citizen).dot(np.random.dirichlet(alpha,citizens))
>>> # NOTE: assert np.allclose(votes.sum(1),votes_per_citizen)
...
>>> # Here is how the votes are counted.
... def count(votes): #NOTE: votes.shape is (citizens,names)
... # Dwarves have a weird tradition: the top 10 popular names chosen
... # by dwarves will get all their votes.
... # (I guess this is what makes our objective non-convex.)
... dwarf_votes = votes[:dwarves]
... dwarf_votes_per_name_ = dwarf_votes.sum(1)
... top_10_idx = np.argsort(-dwarf_votes_per_name_)[:10]
... dwarf_votes_per_name = np.zeros(names)
... dwarf_votes_per_name[top_10_idx] = dwarf_votes_per_name_[top_10_idx]
... dwarf_votes_per_name = \
... dwarf_votes.sum() * dwarf_votes_per_name / dwarf_votes_per_name.sum()
... #NOTE: assert np.allclose(dwarf_votes_per_name.sum(),dwarf_votes_per_name_.sum())
... # Humans,on the other hand,just add everyone's votes together.
... human_votes = votes[dwarves:]
... human_votes_per_name = human_votes.sum(0)
... # The final result is the sum of the scores by dwarves and humans.
... return dwarf_votes_per_name + human_votes_per_name #NOTE: shape = (names,)
...
>>> # So this is the legitimate result.
... votes_per_name = count(votes)
>>>
>>> # Now,you want to cheat a bit and make your proposal
... # (assuming it's the first one) the most popular one.
... votes_per_name_ = votes_per_name.copy()
>>> argmax = votes_per_name.argmax()
>>> votes_per_name_[argmax] = votes_per_name[0]
>>> votes_per_name_[0] = votes_per_name[argmax]
>>> #NOTE: assert np.allclose(votes_per_name.sum(),votes_per_name_.sum())
...
>>> # However,you cannot just manipulate the result,otherwise the auditors
... # will find out! Instead,the raw scores must be manipulated such that your
... # submission ranks top among others. You decide to solve for a multiplier
... # to the original scores.
... coef = np.ones_like(votes).reshape(-1)
>>>
>>> # In other words,your goal is to find the minimum of the following objective.
... def loss(coef): #NOTE: coef.shape is (citizens*votes,)
... votes_ = votes * coef.reshape(votes.shape)
... votes_per_name_ = count(votes_)
... return ((votes_per_name_ - votes_per_name)**2).sum()
...
>>>
>>> # This is a helper constant matrix. Ignore it for now.
... A = np.concatenate([np.tile(np.concatenate([np.repeat(1,names),np.repeat(0,citizens * names)]),citizens - 1),np.repeat(1,names)]).reshape((citizens,citizens * names))
>>> A_ = A.dot(np.diag(votes.reshape(-1)))
>>>
到目前为止,代码与您的问题相同。 现在,让我们使用一些 mystic 的工具。
首先,让我们把它放在您引用的问题的上下文中。
>>> # Build constraints
... cons = my.symbolic.linear_symbolic(A_,votes_per_citizen)
>>> cons = my.symbolic.solve(cons) #NOTE: this may take a while...
>>> cons = my.symbolic.generate_constraint(my.symbolic.generate_solvers(cons))
>>>
>>> def cons_(x): #NOTE: preserve x
... return cons(x.copy())
...
>>>
现在,就像在引用的问题中一样,我们可以使用像 fmin_powell
和 constraints=cons_
这样的简单求解器来解决它。但是,正如我将要演示的那样,我发现神秘主义者在这里可能会有些挣扎。
>>> bounds = [(0,1000)] * coef.size
>>> x0 = np.random.randint(0,1000,coef.shape).astype(float).tolist()
>>>
>>> stepmon = my.monitors.VerboseMonitor(1)
>>> rv = my.solvers.fmin_powell(loss,x0=x0,bounds=bounds,itermon=stepmon,... constraints=cons_,... disp=1,maxfun=20000,gtol=None,ftol=500)
Generation 0 has ChiSquare: inf
Generation 1 has ChiSquare: inf
Generation 2 has ChiSquare: inf
Generation 3 has ChiSquare: inf
Generation 4 has ChiSquare: inf
Generation 5 has ChiSquare: inf
Generation 6 has ChiSquare: inf
Generation 7 has ChiSquare: inf
...
我在这里取消了优化……当您看到 inf
时,这意味着 mystic 正在努力解决约束。
问题是应该有额外的约束来帮助它,即 A_.dot(cons_(x0))
应该是整数。此外,最好能更好地控制底片的存在。
>>> np.round(A_.dot(cons_(x0)),10) - votes_per_citizen
array([0.,0.,0.])
>>> cons_(x0)[:5]
[-574537.1464429945,98.0,326.0,114.0,694.0]
事实上,每个第 citizens
元素都是负数。
所以,我尝试改用惩罚。
>>> def penalty1(x):
... return np.abs((np.array(x).reshape(citizens,names) * votes).sum(1) - votes_per_citizen).sum()
...
>>> def penalty2(x):
... return np.abs(A_.dot(x) - votes_per_citizen).sum()
...
>>> def penalty3(x):
... return np.abs((np.array(x).reshape(citizens,names) * votes).min(1).sum())
...
>>> #FIXME: utilize that A_.dot(x) should be integer-valued
...
>>> @my.penalty.quadratic_equality(penalty2)
... @my.penalty.quadratic_equality(penalty3)
... def penalty(x):
... return 0.0
...
>>>
>>> bounds = [(0,coef.shape).astype(float).tolist()
>>>
>>> stepmon = my.monitors.VerboseMonitor(1)
>>>
>>> rv = my.solvers.fmin_powell(loss,... penalty=penalty,# constraints=cons_,ftol=500)
Generation 0 has ChiSquare: 4316466725165.325684
Generation 1 has ChiSquare: 97299808.307906
Generation 2 has ChiSquare: 1125.438322
Generation 3 has ChiSquare: 1116.735393
Warning: Maximum number of function evaluations has been exceeded.
STOP("EvaluationLimits with {'evaluations': 20000,'generations': 500000}")
>>>
这似乎工作得很好(注意我使用了 20000 evals 而不是 200,当损失是 500 而不是 1 时,ftol 停止)。 rv
非常接近,而且优化速度相当快。
>>> penalty(rv)
4.979032041874226
>>> np.round(A_.dot(rv),2) - votes_per_citizen
array([0.,0.22,0. ])
>>>
>>> penalty(x0)
4313023244843.466
>>> np.round(A_.dot(x0),2) - votes_per_citizen
array([ 6996.61,5872.2,6398.82,7593.65,9137.81,10347.84,9204.44,6160.77,9626.64,7572.4,10771.24,8673.7,10212.18,7581.45,5097.13,7631.2,8274.92,9684.09,9504.27,9067.73,7332.77,10214.02,8255.38,9853.74,6613.19])
此处,A_.dot(rv)
不那么准确(舍入到 2 位而不是 10 位)...并且将再次受益于生成 A_.dot(rv)
整数的约束。
我将把它留给以后的例子。