具有向量、标量和不等式约束的非线性优化

Question

我有一组方程式：Y=aA+bB 其中 Y 是知道浮点数的向量（只有这个是已知的！）； a、b 是未知标量（浮点数），A、B 是浮点数的未知向量。每个方程都有自己的 Y、a、b，而所有方程共享相同的未知向量 A 和 B。

我已经设置了这样的方程，所以我的问题是最小化函数： (Y-aA-bB)+(Y'-a'A-b'B)+.... 我也有许多类型的不等式约束：Ai>Aj（向量 A 的 Ai i-th 元素），Bi>= Bk，Bi>0，a>a'，...

是否有任何软件或库（最适合 python）可以处理这个问题？

Answer 1

我同意 sascha 的观点，即这是一个线性问题。由于我不太喜欢约束，实际上我更喜欢使它成为没有约束的非线性。我这样做是通过像这样设置向量A=(a1**2, a1**2+a2**2, a1**2+a2**2+a3**2, ...) 来确保它都是正的，A_i > A_j 是i>j。这使得错误有点问题，因为您现在必须考虑错误传播以获取 A1、A2 等，包括相关性，但我将在最后提出重要的一点。 “简单”的解决方案如下所示：

import numpy as np
from scipy.optimize import leastsq
from random import random
np.set_printoptions(linewidth=190)


def generate_random_vector(n, sortIt=True):
    out=np.fromiter( (random() for x in range(n) ),np.float)
    if sortIt:
        out.sort()
    return out


def residuals(parameters,dataVec,dataLength,vecDims):
    aParams=parameters[:dataLength]
    bParams=parameters[dataLength:2*dataLength]
    AParams=parameters[-2*vecDims:-vecDims]
    BParams=parameters[-vecDims:]
    YList=dataVec
    AVec=[a**2 for a in AParams]##assures A_i > 0
    BVec=[b**2 for b in BParams]
    AAVec=np.cumsum(AVec)##assures A_i>A_j for i>j
    BBVec=np.cumsum(BVec)
    dist=[ np.array(Y)-a*np.array(AAVec)-b*np.array(BBVec) for Y,a,b in zip(YList,aParams,bParams)  ]
    dist=np.ravel(dist)
    return dist


if __name__=="__main__":

    aList=generate_random_vector(20, sortIt=False)
    bList=generate_random_vector(20, sortIt=False)
    AVec=generate_random_vector(5)
    BVec=generate_random_vector(5)
    YList=[a*AVec+b*BVec for a,b in zip(aList,bList)]

    aGuess=20*[.2]
    bGuess=20*[.3]
    AGuess=5*[.4]
    BGuess=5*[.5]

    bestFitValues, covMX, infoDict, messages ,ier = leastsq(residuals, aGuess+bGuess+AGuess+BGuess ,args=(YList,20,5) ,full_output=True)


    print "a"
    print aList
    besta = bestFitValues[:20]
    print besta

    print "b"
    print bList
    bestb = bestFitValues[20:40]
    print bestb

    print "A"
    print AVec
    bestA = bestFitValues[-2*5:-5]
    realBestA = np.cumsum([x**2 for x in bestA])
    print realBestA

    print "B"
    print BVec
    bestB = bestFitValues[-5:]
    realBestB = np.cumsum([x**2 for x in bestB])
    print realBestB
    print covMX

关于错误和相关性的问题是问题的解决方案不是唯一的。如果Y = a A + b B 是一个解决方案并且我们，例如，旋转使得A = c E + s F 和B = -s E + c F 那么Y = (ac-bs) E + (as+bc) F =e E + f F 也是一个解决方案。因此，参数空间在“解决方案”处完全平坦，导致巨大的错误和世界末日的相关性。

Answer 2

一般说明

这是一个线性问题（至少在线性最小二乘意义上，继续阅读）！

它也没有完全指定，因为不清楚在您的情况下是否应该始终有一个可行的解决方案，或者您是否希望总体上最小化一些给定的损失。您的文字听起来像后者，但在这种情况下，必须选择损失（这在可能的算法方面有所不同）。让我们采用欧几里得范数（可能是这里最好的选择）！

暂时忽略约束，我们可以将此问题视为基本的线性矩阵方程的最小二乘解问题（欧几里得范数与平方欧几里德范数没有区别！） .

min || b - Ax ||^2

这里：

M = number of Y's
N = size of Y
b = (Y0,
     Y1,
     ...) -> shape: M*N (flattened: Y_x = (y_x_0, y_x_1).T)
A = ((a0, 0, 0, ..., b0, 0, 0, ...),
     (0, a0, 0, ..., 0, b0, 0, ...),
     (0, 0, a0, ..., 0, 0, b0, ...),
     ...
     (a1, 0, 0, ..., b1, 0, 0, ...)) -> shape: (M*N, N*2)
x = (A0, A1, A2, ... B0, B1, B2, ...) -> shape: N*2 (one for A, one for B)

你应该怎么做

• 如果不受约束：
  • 转换为标准格式并使用numpy的lstsq
• 如果受约束：
  • 要么使用定制的优化算法，要么：
    • 线性规划（如果最小化绝对差异/l1-norm）
      • 我懒得为scipy的linprog制定它
      • 没那么难，但使用 scipy 的 API，l1-norm 并非易事
      • 使用cvxpy (obj=cvxpy.norm(X, 1)) 更容易制定
    • 二次规划/二阶锥规划（如果最小化欧几里得范数/l2-范数）
      • 再次，懒得制定它； scipy 上还没有特殊的求解器
      • 可以使用 cvxpy (obj=cvxpy.norm(X, 2)) 轻松制定
    • 紧急情况：使用通用的约束非线性优化算法，如 SLSQP -> 参见代码

一些 hacky 代码（不是最好的方法！）

这段代码：

• 只是一个演示！
• 使用来自 scipy 的一般非线性优化算法
  • 因此：
    • 更容易制定
    • 不如 LP、QP、SOCP 快速和稳健
    • 但在保证凸优化问题收敛的情况下，将获得大致相同的结果
  • 在需要时使用自动微分
    • （作者懒得加渐变）
    • 如果性能很重要，这真的很糟糕
• 在 np.repeat 与广播方面真的很难看！

代码：

import numpy as np
from scipy.optimize import minimize
np.random.seed(1)

""" Fake-problem (usually the job of the question-author!) """
def get_partial(N=10):
    Y = np.random.uniform(size=N)
    a, b = np.random.uniform(size=2)
    return Y, a, b

""" Optimization """
def optimize(list_partials, N, M):
    """ General approach:
        This is a linear system of equations (with constraints)
        Basic (unconstrained) form: min || b - Ax ||^2
     """
    Y_all = np.vstack(map(lambda x: x[0], list_partials)).ravel()       # flat 1d
    a_all = np.hstack(map(lambda x: np.repeat(x[1], N), list_partials)) # repeat to be of same shape
    b_all = np.hstack(map(lambda x: np.repeat(x[2], N), list_partials)) # """

    def func(x):
        A = x[:N]
        B = x[N:]
        return np.linalg.norm(Y_all - a_all * np.repeat(A, M) - b_all * np.repeat(B, M))

    """ Example constraints: A >= B element-wise """
    cons = ({'type': 'ineq',
             'fun' : lambda x: x[:N] - x[N:]})

    res = minimize(func, np.zeros(N*2), constraints=cons, method='SLSQP', options={'disp': True})
    print(res)
    print(Y_all - a_all * np.repeat(res.x[:N], M) - b_all * np.repeat(res.x[N:], M))


""" Test """
M = 4
N = 3
list_partials = [get_partial(N) for i in range(M)]
optimize(list_partials, N, M)

输出：

Optimization terminated successfully.    (Exit mode 0)
            Current function value: 0.9019356096498999
            Iterations: 12
            Function evaluations: 96
            Gradient evaluations: 12
     fun: 0.9019356096498999
     jac: array([  1.03786588e-04,   4.84041870e-04,   2.08129734e-01,
         1.57609582e-04,   2.87599862e-04,  -2.07959406e-01])
 message: 'Optimization terminated successfully.'
    nfev: 96
     nit: 12
    njev: 12
  status: 0
 success: True
       x: array([ 1.82177105,  0.62803449,  0.63815278, -1.16960281,  0.03147683,
        0.63815278])
[  3.78873785e-02   3.41189867e-01  -3.79020251e-01  -2.79338679e-04
  -7.98836875e-02   7.94168282e-02  -1.33155595e-01   1.32869391e-01
  -3.73398306e-01   4.54460178e-01   2.01297470e-01   3.42682496e-01]

我没有检查结果！如果有错误，那就是实现错误，而不是概念错误（我的观点）！