问题描述
我正在尝试对向量 [a1,a2,a3] 执行最小二乘优化,具有以下约束,其中 k
是一个常量:
-k < a3/a2 < k
我将发布我的代码的简化版本,希望这能让我清楚地了解我的目标。
from scipy import optimize
def loss_function(a_candidate):
return MyObject(a_candidate).mean_square_error()
def feasibility(a_candidate):
# Returns a boolean
k = 1.66711
a2 = a_candidate[1]
a3 = a_candidate[2]
return -k < a3/a2 < k
def feasibility_scipy(a_candidate):
# As I understand it,for SciPy the constraint has to be a number
boolean = feasibility(a_candidate)
if boolean:
return 0
else:
return 1
# Equality constraint forcing feasibility_scipy to be 0.
constraint = optimize.NonlinearConstraint(feasibility_scipy,0)
optimize_results = optimize.minimize(
loss_function,a_init,# Initial guess
constraints=constraint)
由于初始猜测 a_init
的生成方式,它位于 feasibility
为 False
的区域。 (首先我们需要使用数值方法的原因是早期的封闭形式方法返回了一个不可行的解决方案。可以提供一个非常差的可行猜测,如 (0,0),但这将是很多离真正的解决方案更远)。
由于 constraint
的梯度几乎处处都是零,因此优化例程无法找到离开这个不可行(不可接受)区域的方法,并且不会成功终止。使用 SLSQP
,它会在 1 次迭代后停止,并显示消息 Singular matrix C in LSQ subproblem
。使用trust-constr
求解器,达到了最大函数求值次数,我相信它没有离开不可行区域,因为constr_violation
是1.0
。
据我所知,在 SciPy 中不可能在 a2
和 a3
上提供“联合绑定”(对发明的术语表示歉意),这意味着我不得不使用 {{1 }} 方法。
这样做的正确方法是什么? (一些搜索表明我可能想尝试带有符号约束的 NonlinearConstraint
包。但在我投入时间学习这个新包之前,我想看看 StackOverflow 是否有基于 SciPy 的解决方案。或者如果您知道如何在 mystic
中执行此操作,一些示例代码会很有帮助。)
解决方法
没关系,我认为解决方案可能只是:
def a_ratio(a_candidate):
a2 = a_candidate[1]
a3 = a_candidate[2]
return a3/a2
feasibility_constraint = optimize.NonlinearConstraint(a_ratio,-k,k)
,
我不使用 -k < a3/a2 < k
而是使用这两个线性约束:
-k*a2 <= a3
a3 <= k*a2
,
我知道问题已经得到解答,但您想知道如何在 mystic
中解决它。我是 mystic
作者。像这样:
>>> import mystic as my
>>> import numpy as np
>>>
>>> truth = np.array([1,1,1])
>>>
>>>
>>> class MyObject(object):
... def __init__(self,candidate):
... self.candidate = candidate
... self.truth = truth
... def mean_square_error(self):
... return ((self.truth - self.candidate)**2).sum()
...
>>>
>>> def loss_function(a_candidate):
... return MyObject(a_candidate).mean_square_error()
...
>>>
然后是解决方案:
>>> equations = """
... -k < a3/a2
... a3/a2 < k
... """
>>>
>>> eqn = my.symbolic.simplify(equations,variables=['a1','a2','a3'],locals=dict(k=1.66711),all=False,target='a3')
>>> c = my.symbolic.generate_constraint(my.symbolic.generate_solvers(eqn,'a3']),join=my.constraints.or_)
>>> res = my.solvers.fmin(loss_function,x0=[2,3,4],constraints=c,disp=True,xtol=1e-8)
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 121
Function evaluations: 229
>>> res
array([1.,1.,1.])
mystic
的不同之处在于它构建了一个“运算符”c
,不允许优化器选择不可行的解决方案。并不是说它用布尔值拒绝糟糕的解决方案,而是将“坏”候选者转换为“好”候选者,如下所示:
>>> c([2,10,4])
[2,4]
>>> c([2,1.6671099999999974]