具有关节边界/符号约束的优化

Question

我正在尝试对向量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, 0)


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

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

由于constraint的梯度几乎处处为零，优化例程无法找到走出这个不可行(不可接受)区域的方法，因此它不会成功终止。使用SLSQP，优化例程仅在1次迭代后停止，消息为Singular matrix C in LSQ subproblem。使用trust-constr求解器，它达到函数求值的最大次数，我相信它没有离开不可行域，因为constr_violation为1.0。

据我所知，在SciPy中不可能在a2和a3上提供“联合边界”(为发明的术语道歉)，这意味着我被迫使用NonlinearConstraint方法。

这样做的正确方法是什么？(一些搜索表明，我可能想尝试使用带有符号约束的mystic包。但是在我花时间学习这个新的包之前，我想看看StackOverflow是否有一个基于SciPy的解决方案。或者，如果您知道如何在mystic中做到这一点，那么一些示例代码将非常有用。)

Answer 1

不要紧，我认为解决方案可能就是：

def a_ratio(a_candidate):
    a2 = a_candidate[1]
    a3 = a_candidate[2]
    return a3/a2

feasibility_constraint = optimize.NonlinearConstraint(a_ratio,-k,k)

Answer 2

我不使用-k < a3/a2 < k，而是使用以下两个线性约束：

     -k*a2 <= a3
     a3 <= k*a2 

Answer 3

我知道这个问题已经回答了，但是您想知道如何在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，不允许优化器选择不可行的解决方案。这并不是说它用布尔值拒绝糟糕的解决方案，而是它将“糟糕的”候选者转换为“好的”候选者，如下所示：

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