SciPy 在约束函数中使用附加变量进行优化
SciPy optimize with additional variables in constraint function
我大多是 Python 和 Whosebug 的新手,所以请指出我是否正在做任何事情 wrong/need 以包含更多详细信息 e.t.c。我的问题围绕着最小化函数 en_sum
,它取决于变量数组 en
。此优化必须受到约束 nlc
:取决于变量 en
以及另一组变量 x2
的矩阵的最小特征值必须为正。
我想做的基本上是提供一些 en
和 x2
的猜测值,根据 nlc
优化 en_sum
并让它产生优化结果en
和 x2
值。
我写的最小代码如下:
import numpy as np
from scipy import optimize as spo
from scipy.optimize import NonlinearConstraint
def en_min(varia):
#en_min is a function that (aims to) return the optimised values of en and x2.
m = 2
n = 2
s = 2
en = varia[0:m]
x2 = varia[m:len(varia)]
# use_mat(a,b) defines the matrix which depends on en and x2 variables.
def use_mat(a, b):
mat = np.empty([m * s, n * s])
for i in range(0, m * s):
for j in range(0, m * s):
mat[i, j] = (b[j] * a[i - s]) ** 2 + a[j - s] ** 2 + (b[j] + b[i]) ** 2
return mat
# Defining the optimisation function.
def en_sum(q):
return sum(q[iii] ** 2 for iii in range(len(q)))
# This function provides the minimal eigenvalue of a given matrix.
def min_eigenvalue(a):
y = np.linalg.eigvals(a)
return np.amin(y)
# My attempt at setting the constraint (which is that the eigenvalues of the matrix
# defined through use_mat(a,b) are non-negative) is nlc
nlc = NonlinearConstraint(min_eigenvalue(use_mat(en, x2)), 0, np.inf)
# The guess values in opt_res2 are just for the en variables, as en_sum must
# only be a function of those (not x2). However, I do not know how to properly give guess
# values to x2 in the constraint, so I am quite sure this is implemented incorrectly.
opt_res2 = spo.minimize(en_sum(en), np.array([1, 4]), constraints=nlc)
return opt_res2
print(en_min([1, 3, 1.5, 0, 0, 4]))
当然,这行不通(我收到错误 TypeError: 'numpy.complex128' object is not callable
等),我想知道是否有人对我可能出错的地方有任何提示,以及在这种情况下问题是否更容易处理?我的主要困惑来自不知道如何正确处理约束以及如何为 x2
提供猜测值(这与约束字典中的 args 相关吗?)。非常感谢任何想法,谢谢。
有趣的问题。 IMO,最佳做法是将 objective 函数和所有约束编写为一个(矢量)优化变量 z 的函数。在您的函数中,您可以简单地设置 z = (en, x2) 并将变量 z 拆分回 en 和 x2。因此,您的 objective 函数可以写成如下形式:
def en_sum(z):
en, x2 = z[:m], z[m:]
return np.sum(en**2)
另请注意,en_sum(en)
和 min_eigenvalues(use_mat(en, x2))
不是函数,因此不可调用。此外,最小特征值必须为正(use_mat 正定返回的矩阵)的约束等同于所有特征值均为正。这部分有点混乱,因为找到一个可行的起点来满足这个约束并不是一件容易的事。然而,我们可以通过只考虑复杂特征值的实部来作弊。因此,即使当前点不可行,约束也会得到满足。在这里,我们基本上希望它在前几次迭代中仍然收敛到可行域(所有特征值都是实数和正数)。
在代码中:
import numpy as np
from scipy.optimize import minimize, NonlinearConstraint
def en_min(varia):
#en_min is a function that (aims to) return the optimised values of en and x2.
m = 2
n = 2
s = 2
en = varia[0:m]
x2 = varia[m:len(varia)]
# use_mat(a,b) defines the matrix which depends on en and x2 variables.
def use_mat(z):
a, b = z[:m], z[m:]
mat = np.empty([m * s, n * s])
for i in range(0, m * s):
for j in range(0, m * s):
mat[i, j] = (b[j]*a[i-s])** 2 + a[j-s]**2 + (b[j]+b[i])**2
return mat
# Defining the optimisation function.
def en_sum(z):
en, x2 = z[:m], z[m:]
return np.sum(en**2)
# This function provides the minimal eigenvalue of a given matrix.
def min_eigenvalue(a):
return np.real(np.linalg.eigvals(a))
# constraint
nlc = NonlinearConstraint(lambda z: min_eigenvalue(use_mat(z)), 1e-12, np.inf)
# Bounds on the variables
bounds = [(None, None)]*en.size + [(0, None)]*x2.size
# initial guess
x0 = np.hstack((en, 0.1*np.ones(x2.size)))
opt_res2 = minimize(en_sum, x0=x0, bounds=bounds, constraints=nlc)
return opt_res2
print(en_min(np.array([1, 3, 1.5, 0, 0, 4])))
或者,由于矩阵 A 是正定的当且仅当存在矩阵 V 使得 A = V'V,因此可以解决辅助问题以找到可行的初始猜测。
还值得一提的是,由于雅可比行列式(=我们的约束函数的导数)的逼近,上面的代码并不是很有效,并且对于大型优化问题会很慢:目前,each 雅可比行列式的计算需要 N*N 次计算 min_eigenvalue
,其中 N 是变量的总数。因此,我们计算 N*N 次每次迭代的特征值。总而言之,对于更大的问题,绝对值得提供精确的雅可比矩阵或使用算法微分
我大多是 Python 和 Whosebug 的新手,所以请指出我是否正在做任何事情 wrong/need 以包含更多详细信息 e.t.c。我的问题围绕着最小化函数 en_sum
,它取决于变量数组 en
。此优化必须受到约束 nlc
:取决于变量 en
以及另一组变量 x2
的矩阵的最小特征值必须为正。
我想做的基本上是提供一些 en
和 x2
的猜测值,根据 nlc
优化 en_sum
并让它产生优化结果en
和 x2
值。
我写的最小代码如下:
import numpy as np
from scipy import optimize as spo
from scipy.optimize import NonlinearConstraint
def en_min(varia):
#en_min is a function that (aims to) return the optimised values of en and x2.
m = 2
n = 2
s = 2
en = varia[0:m]
x2 = varia[m:len(varia)]
# use_mat(a,b) defines the matrix which depends on en and x2 variables.
def use_mat(a, b):
mat = np.empty([m * s, n * s])
for i in range(0, m * s):
for j in range(0, m * s):
mat[i, j] = (b[j] * a[i - s]) ** 2 + a[j - s] ** 2 + (b[j] + b[i]) ** 2
return mat
# Defining the optimisation function.
def en_sum(q):
return sum(q[iii] ** 2 for iii in range(len(q)))
# This function provides the minimal eigenvalue of a given matrix.
def min_eigenvalue(a):
y = np.linalg.eigvals(a)
return np.amin(y)
# My attempt at setting the constraint (which is that the eigenvalues of the matrix
# defined through use_mat(a,b) are non-negative) is nlc
nlc = NonlinearConstraint(min_eigenvalue(use_mat(en, x2)), 0, np.inf)
# The guess values in opt_res2 are just for the en variables, as en_sum must
# only be a function of those (not x2). However, I do not know how to properly give guess
# values to x2 in the constraint, so I am quite sure this is implemented incorrectly.
opt_res2 = spo.minimize(en_sum(en), np.array([1, 4]), constraints=nlc)
return opt_res2
print(en_min([1, 3, 1.5, 0, 0, 4]))
当然,这行不通(我收到错误 TypeError: 'numpy.complex128' object is not callable
等),我想知道是否有人对我可能出错的地方有任何提示,以及在这种情况下问题是否更容易处理?我的主要困惑来自不知道如何正确处理约束以及如何为 x2
提供猜测值(这与约束字典中的 args 相关吗?)。非常感谢任何想法,谢谢。
有趣的问题。 IMO,最佳做法是将 objective 函数和所有约束编写为一个(矢量)优化变量 z 的函数。在您的函数中,您可以简单地设置 z = (en, x2) 并将变量 z 拆分回 en 和 x2。因此,您的 objective 函数可以写成如下形式:
def en_sum(z):
en, x2 = z[:m], z[m:]
return np.sum(en**2)
另请注意,en_sum(en)
和 min_eigenvalues(use_mat(en, x2))
不是函数,因此不可调用。此外,最小特征值必须为正(use_mat 正定返回的矩阵)的约束等同于所有特征值均为正。这部分有点混乱,因为找到一个可行的起点来满足这个约束并不是一件容易的事。然而,我们可以通过只考虑复杂特征值的实部来作弊。因此,即使当前点不可行,约束也会得到满足。在这里,我们基本上希望它在前几次迭代中仍然收敛到可行域(所有特征值都是实数和正数)。
在代码中:
import numpy as np
from scipy.optimize import minimize, NonlinearConstraint
def en_min(varia):
#en_min is a function that (aims to) return the optimised values of en and x2.
m = 2
n = 2
s = 2
en = varia[0:m]
x2 = varia[m:len(varia)]
# use_mat(a,b) defines the matrix which depends on en and x2 variables.
def use_mat(z):
a, b = z[:m], z[m:]
mat = np.empty([m * s, n * s])
for i in range(0, m * s):
for j in range(0, m * s):
mat[i, j] = (b[j]*a[i-s])** 2 + a[j-s]**2 + (b[j]+b[i])**2
return mat
# Defining the optimisation function.
def en_sum(z):
en, x2 = z[:m], z[m:]
return np.sum(en**2)
# This function provides the minimal eigenvalue of a given matrix.
def min_eigenvalue(a):
return np.real(np.linalg.eigvals(a))
# constraint
nlc = NonlinearConstraint(lambda z: min_eigenvalue(use_mat(z)), 1e-12, np.inf)
# Bounds on the variables
bounds = [(None, None)]*en.size + [(0, None)]*x2.size
# initial guess
x0 = np.hstack((en, 0.1*np.ones(x2.size)))
opt_res2 = minimize(en_sum, x0=x0, bounds=bounds, constraints=nlc)
return opt_res2
print(en_min(np.array([1, 3, 1.5, 0, 0, 4])))
或者,由于矩阵 A 是正定的当且仅当存在矩阵 V 使得 A = V'V,因此可以解决辅助问题以找到可行的初始猜测。
还值得一提的是,由于雅可比行列式(=我们的约束函数的导数)的逼近,上面的代码并不是很有效,并且对于大型优化问题会很慢:目前,each 雅可比行列式的计算需要 N*N 次计算 min_eigenvalue
,其中 N 是变量的总数。因此,我们计算 N*N 次每次迭代的特征值。总而言之,对于更大的问题,绝对值得提供精确的雅可比矩阵或使用算法微分