使用联合边界/符号约束进行优化

Optimization with joint bounds / symbolic constraint

我正在尝试对向量 [a1,a2,a3] 执行最小二乘优化,具有以下约束,其中 k 是常数:

-k < a3/a2 < k

我将 post 我的代码的简化版本,希望这能让我清楚地知道我想要的是什么。

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, 0)


optimize_results = optimize.minimize(
    loss_function,
    a_init,  # Initial guess
    constraints=constraint)

由于初始猜测 a_init 的生成方式,它位于 feasibilityFalse 的区域。 (我们首先需要使用数值方法的原因是早期的封闭形式方法返回了一个不可行的解决方案。有可能提供一个非常差的可行猜测,如 (0,0,0),但这将是离真正的解决方案更远)。

由于constraint的梯度几乎处处为零,优化例程无法找到离开这个不可行(不可接受)区域的出路,并且没有成功终止。使用 SLSQP 它仅在 1 次迭代后停止,并显示消息 Singular matrix C in LSQ subproblem。使用 trust-constr 求解器,它达到了函数计算的最大次数,我相信它没有离开不可行区域,因为 constr_violation1.0.

据我了解,在 SciPy 中不可能在 a2a3 上提供 'joint bound'(对发明的术语表示歉意),意思是我被迫使用 NonlinearConstraint 方法。

正确的做法是什么? (一些搜索建议我可能想尝试带有符号约束的 mystic 包。但在我花时间学习这个新包之前,我想看看 Whosebug 是否有基于 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, variables=['a1','a2','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,不允许优化器 select 不可行的解决方案。并不是说它用布尔值拒绝不好的解决方案,而是它将“坏”候选者转换为“好”候选者,如下所示:

>>> c([2,10,4])
[2, 10, 4]
>>> c([2,1,4])
[2, 1, 1.6671099999999974]