使用Scipy的最小化（scipy.optimize.minimize）与大型等式约束矩阵

Question

我需要最小化一个有五个变量（x[0]到x[4]）的函数。

这个需要最小化的标量函数可以表示为 X'*H*X。目标函数看起来大概是这样的：

def objfun(x):
    H = 0.1*np.ones([5,5])
    f = np.dot(np.transpose(x),np.dot(H,x))[0][0]
    return f

这个函数会返回一个单一的标量值。

我的问题是，如何实现一个约束方程，形式如下：

A*X - b = 0

其中 A 和 b 在每次运行时可能会变化。举个随机的例子：

A = 
array([[ 1,  2,  3,  4,  5],
       [ 2,  1,  3,  4,  5],
       [-1,  2,  3,  0,  0],
       [ 0, -5,  6,  3,  2],
       [-3,  5,  6,  2,  8]])

B = 
array([[ 0],
       [ 2],
       [ 3],
       [-2],
       [-7]])

A 和 B 不能直接写死在约束函数里，因为每次运行时它们可能不同。变量没有上下限，优化方法也不需要特别指定。

编辑

我意识到，对于一个有5个变量的优化问题，如果有5个约束方程，那么通过解这些方程可以得到一个唯一的解。那么如果 A 被定义为：

A = 
array([[ 1,  2,  3,  4,  5],
       [ 2,  1,  3,  4,  5],
       [-1,  2,  3,  0,  0],
       [ 0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0]])


B = 
array([[ 0],
       [ 2],
       [ 3],
       [ 0],
       [ 0]])

那么我们就有一个有5个变量和3个线性约束的优化问题。

Answer 1

NLopt的文档提到了一种很不错的通用方法：
所有的方程 Ax = b 的解都可以表示成 xany + nullspace(A) z 的形式，
其中 xany 是其中一个解，而 dim(z) < dim(x) 。
所以我们可以在没有限制的情况下，最小化 f( xany + nullspace(A) z )，这里的 z 是可以自由变化的。

举个例子，在三维空间中，约束条件 x0 + x1 + x2 = 1 的零空间矩阵是：

[  1  0 ]  : [z0 z1] -> [z0, -z0 + z1, -z1]  -- sum 0
[ -1  1 ]
[  0 -1 ]

("在数值计算零空间时需要小心...")

Answer 2

你可以试试使用 scipy.optimize.fmin_cobyla 这个函数。我对具体的数值细节不是很了解，所以你应该用一些你知道答案的值来测试一下，看看它是否符合你的需求。同时，可以调整一下容差参数 rhoend 和 rhobeg，看看能否得到预期的结果。一个示例程序可能是这样的：

import numpy as np
import scipy.optimize

A = \
np.array([[ 1,  2,  3,  4,  5],
          [ 2,  1,  3,  4,  5],
          [-1,  2,  3,  0,  0],
          [ 0,  0,  0,  0,  0],
          [ 0,  0,  0,  0,  0]])

B = \
np.array([[0],
          [2],
          [3],
          [0],
          [0]])

def objfun(x):
    H = 0.1*np.ones([5,5])
    f = np.dot(np.transpose(x),np.dot(H,x))
    return f

def constr1(x):
    """ The constraint is satisfied when return value >=0 """
    sol = A*x
    if np.allclose(sol, B):
        return 0.01
    else:
        # Return the negative distance between the expected solution
        # and the actual solution for a somehow meaningful value
        return -np.linalg.norm(B-sol)

scipy.optimize.fmin_cobyla(objfun, [0.0, 0.0, 0.0, 0.0,0.0], [constr1])
#np.linalg.solve(A, b)

请注意，这个例子没有解决方案，最好用一个有解决方案的例子来试试。我不太确定约束函数是否定义得正确，建议你找一个适合你的情况的函数。为了获得更好的结果，建议你提供一个实际的初始猜测，而不是 [0.0, 0.0, 0.0, 0.0, 0.0]。

想了解更多细节，可以查看官方文档：http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fmin_cobyla.html#scipy.optimize.fmin_cobyla

编辑：根据你想要的解决方案类型，你可能可以构建一个更好的约束函数，允许一些值在预期解决方案的某个容差范围内，即使不是完全准确，并且返回一个值，越接近容差这个值就越高，而不是总是返回0.1等等……