使用 mystic 解决依赖于参数的优化问题
Using mystic to solve a parameter-dependent optimization problem
我有一个类型为
的非凸二次优化问题
x' * B * x,
其中 x 的所有条目都在 0 和 1 之间,并且所有条目的总和等于 1。
在scipy.optimize中我会尝试通过
解决这个优化问题
import numpy as np
from scipy.optimize import minimize, LinearConstraint
N = 2 # dimension 2 for this example
B = np.array([[2,-1],[-1,-1]]) # symmetric, but indefinite matrix
fnc = lambda x: x.T @ B @ x
res = minimize(fnc, x0 = np.ones((N,))/N, bounds = [(0,1) for m in range(N)], constraints = (LinearConstraint(np.ones((N,)),0.99, 1.01)))
所以我从初始猜测 [0.5, 0.5] 开始,我在每个维度上应用边界 (0,1) 并且等式约束由非常窄的双重不等式约束处理。
现在我想把它翻译成 mystic,因为 scipy 不能很好地处理高维非凸设置(我感兴趣)。
我没能找到的是如何以一种形式编写约束,以便我只需要提供具有可变维度的矩阵 B。到目前为止,我在 mystic 中找到的所有示例都是这样做的:
def objective(x):
x0,x1,x2,x3,x4,x5,x6,x7,x8,x9 = x
return x0**2 + x1**2 + x0*x1 - 14*x0 - 16*x1 + (x2-10)**2 + \
4*(x3-5)**2 + (x4-3)**2 + 2*(x5-1)**2 + 5*x6**2 + \
7*(x7-11)**2 + 2*(x8-10)**2 + (x9-7)**2 + 45.0
bounds = [(-10,10)]*10
from mystic.symbolic import generate_constraint, generate_solvers, simplify
from mystic.symbolic import generate_penalty, generate_conditions
equations = """
4.0*x0 + 5.0*x1 - 3.0*x6 + 9.0*x7 - 105.0 <= 0.0
10.0*x0 - 8.0*x1 - 17.0*x6 + 2.0*x7 <= 0.0
-8.0*x0 + 2.0*x1 + 5.0*x8 - 2.0*x9 - 12.0 <= 0.0
3.0*(x0-2)**2 + 4.0*(x1-3)**2 + 2.0*x2**2 - 7.0*x3 - 120.0 <= 0.0
5.0*x0**2 + 8.0*x1 + (x2-6)**2 - 2.0*x3 - 40.0 <= 0.0
0.5*(x0-8)**2 + 2.0*(x1-4)**2 + 3.0*x4**2 - x5 - 30.0 <= 0.0
x0**2 + 2.0*(x1-2)**2 - 2.0*x0*x1 + 14.0*x4 - 6.0*x5 <= 0.0
-3.0*x0 + 6.0*x1 + 12.0*(x8-8)**2 - 7.0*x9 <= 0.0
"""
cf = generate_constraint(generate_solvers(simplify(equations, target=['x5','x3'])))
pf = generate_penalty(generate_conditions(equations))
这非常冗长,需要手动插入所有约束和参数等作为我想避免的字符串:每次我需要 运行优化方法。我想要的(在一个完美的世界里)会是这样的
def objective(x):
return x @ B @ x # numpy syntax
equations = """
np.ones((1,N)) @ x == 1.0
"""
# constraint in a form which can handle variable dimension of x
这可能吗?
Mystic 默认使用列表,因此您必须在成本函数中转换为数组。还有很多其他方法可以在不使用符号字符串的情况下创建约束,在您的特定情况下,有一种开箱即用的方法。我会这样做:
>>> import mystic as my
>>> import numpy as np
>>> N = 2 # dimension 2 for this example
>>> B = np.array([[2,-1],[-1,-1]]) # symmetric, but indefinite matrix
>>> c = my.constraints.normalized()(lambda x:x)
>>> bounds = [(0,1)]*N
>>> mon = my.monitors.VerboseMonitor(10)
>>> fnc = lambda x: np.array(x).T @ B @ x
>>> res = my.solvers.diffev2(fnc, x0=bounds, npop=10, bounds=bounds, ftol=1e-4, gtol=100, full_output=1, itermon=mon, constraints=c)
Generation 0 has ChiSquare: -0.920151
Generation 10 has ChiSquare: -0.999667
Generation 20 has ChiSquare: -1.000000
Generation 30 has ChiSquare: -1.000000
Generation 40 has ChiSquare: -1.000000
Generation 50 has ChiSquare: -1.000000
Generation 60 has ChiSquare: -1.000000
Generation 70 has ChiSquare: -1.000000
Generation 80 has ChiSquare: -1.000000
Generation 90 has ChiSquare: -1.000000
Generation 100 has ChiSquare: -1.000000
Generation 110 has ChiSquare: -1.000000
STOP("ChangeOverGeneration with {'tolerance': 0.0001, 'generations': 100}")
Optimization terminated successfully.
Current function value: -1.000000
Iterations: 113
Function evaluations: 1140
>>> res[0]
array([1.07421473e-07, 9.99999993e-01])
>>> res[1]
-1.0000001999996087
>>> my.scripts.log_reader(mon)
我有一个类型为
的非凸二次优化问题x' * B * x,
其中 x 的所有条目都在 0 和 1 之间,并且所有条目的总和等于 1。
在scipy.optimize中我会尝试通过
解决这个优化问题import numpy as np
from scipy.optimize import minimize, LinearConstraint
N = 2 # dimension 2 for this example
B = np.array([[2,-1],[-1,-1]]) # symmetric, but indefinite matrix
fnc = lambda x: x.T @ B @ x
res = minimize(fnc, x0 = np.ones((N,))/N, bounds = [(0,1) for m in range(N)], constraints = (LinearConstraint(np.ones((N,)),0.99, 1.01)))
所以我从初始猜测 [0.5, 0.5] 开始,我在每个维度上应用边界 (0,1) 并且等式约束由非常窄的双重不等式约束处理。
现在我想把它翻译成 mystic,因为 scipy 不能很好地处理高维非凸设置(我感兴趣)。
我没能找到的是如何以一种形式编写约束,以便我只需要提供具有可变维度的矩阵 B。到目前为止,我在 mystic 中找到的所有示例都是这样做的:
def objective(x):
x0,x1,x2,x3,x4,x5,x6,x7,x8,x9 = x
return x0**2 + x1**2 + x0*x1 - 14*x0 - 16*x1 + (x2-10)**2 + \
4*(x3-5)**2 + (x4-3)**2 + 2*(x5-1)**2 + 5*x6**2 + \
7*(x7-11)**2 + 2*(x8-10)**2 + (x9-7)**2 + 45.0
bounds = [(-10,10)]*10
from mystic.symbolic import generate_constraint, generate_solvers, simplify
from mystic.symbolic import generate_penalty, generate_conditions
equations = """
4.0*x0 + 5.0*x1 - 3.0*x6 + 9.0*x7 - 105.0 <= 0.0
10.0*x0 - 8.0*x1 - 17.0*x6 + 2.0*x7 <= 0.0
-8.0*x0 + 2.0*x1 + 5.0*x8 - 2.0*x9 - 12.0 <= 0.0
3.0*(x0-2)**2 + 4.0*(x1-3)**2 + 2.0*x2**2 - 7.0*x3 - 120.0 <= 0.0
5.0*x0**2 + 8.0*x1 + (x2-6)**2 - 2.0*x3 - 40.0 <= 0.0
0.5*(x0-8)**2 + 2.0*(x1-4)**2 + 3.0*x4**2 - x5 - 30.0 <= 0.0
x0**2 + 2.0*(x1-2)**2 - 2.0*x0*x1 + 14.0*x4 - 6.0*x5 <= 0.0
-3.0*x0 + 6.0*x1 + 12.0*(x8-8)**2 - 7.0*x9 <= 0.0
"""
cf = generate_constraint(generate_solvers(simplify(equations, target=['x5','x3'])))
pf = generate_penalty(generate_conditions(equations))
这非常冗长,需要手动插入所有约束和参数等作为我想避免的字符串:每次我需要 运行优化方法。我想要的(在一个完美的世界里)会是这样的
def objective(x):
return x @ B @ x # numpy syntax
equations = """
np.ones((1,N)) @ x == 1.0
"""
# constraint in a form which can handle variable dimension of x
这可能吗?
Mystic 默认使用列表,因此您必须在成本函数中转换为数组。还有很多其他方法可以在不使用符号字符串的情况下创建约束,在您的特定情况下,有一种开箱即用的方法。我会这样做:
>>> import mystic as my
>>> import numpy as np
>>> N = 2 # dimension 2 for this example
>>> B = np.array([[2,-1],[-1,-1]]) # symmetric, but indefinite matrix
>>> c = my.constraints.normalized()(lambda x:x)
>>> bounds = [(0,1)]*N
>>> mon = my.monitors.VerboseMonitor(10)
>>> fnc = lambda x: np.array(x).T @ B @ x
>>> res = my.solvers.diffev2(fnc, x0=bounds, npop=10, bounds=bounds, ftol=1e-4, gtol=100, full_output=1, itermon=mon, constraints=c)
Generation 0 has ChiSquare: -0.920151
Generation 10 has ChiSquare: -0.999667
Generation 20 has ChiSquare: -1.000000
Generation 30 has ChiSquare: -1.000000
Generation 40 has ChiSquare: -1.000000
Generation 50 has ChiSquare: -1.000000
Generation 60 has ChiSquare: -1.000000
Generation 70 has ChiSquare: -1.000000
Generation 80 has ChiSquare: -1.000000
Generation 90 has ChiSquare: -1.000000
Generation 100 has ChiSquare: -1.000000
Generation 110 has ChiSquare: -1.000000
STOP("ChangeOverGeneration with {'tolerance': 0.0001, 'generations': 100}")
Optimization terminated successfully.
Current function value: -1.000000
Iterations: 113
Function evaluations: 1140
>>> res[0]
array([1.07421473e-07, 9.99999993e-01])
>>> res[1]
-1.0000001999996087
>>> my.scripts.log_reader(mon)