差分进化没有找到全局最小值

Question

我正在生成一个 Erdos-Renyi 网络并计算其度数分布。我知道它接近Poisson distribution，所以我想将我的经验分布拟合到理论分布。为了确定理论分布的参数，我使用来自 SciPy 的 differential evolution 最小二乘法。这是我的代码：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import networkx as nx
from scipy.stats import poisson, expon
from scipy.optimize import curve_fit, differential_evolution
from numpy.polynomial import Polynomial
from scipy.integrate import quad
from functools import partial

def get_deg_distribution(G):
    deg_sequence = sorted([deg for _, deg in G.degree()], reverse=True)
    deg_distribution = pd.DataFrame(np.column_stack(np.unique(deg_sequence, return_counts=True)),
                                    columns=['Degree', 'Count'])
    deg_distribution['Probability'] = deg_distribution['Count'] / deg_distribution['Count'].sum()

    deg_span = np.arange(deg_distribution['Degree'].min(), deg_distribution['Degree'].max() + 1)
    missing_degs = list(set(deg_span) - set(deg_distribution['Degree']))
    add_df = pd.DataFrame(0.0, index=missing_degs, columns=['Count', 'Probability']).reset_index()
    add_df.columns = deg_distribution.columns
    deg_distribution = pd.concat([deg_distribution, add_df]).sort_values('Degree').reset_index(drop=True)

    return deg_distribution

G = nx.erdos_renyi_graph(2500, 0.025, seed=42)
dist = get_deg_distribution(G)

xdata = list(dist['Degree'])
ydata = list(dist['Probability'])

def func(k, mu, loc):
    return poisson.pmf(k, mu, loc)

def param_func(params, ydata):
    return np.linalg.norm(partial(func, xdata)(*params) - ydata)**2

bounds = [(0, 100), (-5, 5)]

params_found = differential_evolution(param_func, bounds=bounds, args=(ydata,))['x']

这段代码有几个谜团。首先，不同的差分进化运行会产生不同的参数集。这本身不是问题，因为确实可能存在许多全局最小值。但是，这些参数集中的 func(xdata, *params) 总是产生零，就好像该方法正在最小化 func 本身而不是其平方误差 参数函数。第三，所有这些参数集的 param_func 值正好是 0.03708384。这是一个问题 - 这不是全球最低限度！只需采用 (65, 0) 参数集，它会给你合理的分布近似值，如果你检查图，误差实际上是 0.002 左右。我之前在其他函数上运行过管道本身，试试这样的：

def func(x, a, b):
    return a*x**2 + b

xdata = np.linspace(0, 10, 11)
ydata = func(xdata, 1, 5)

def param_func(params, ydata):
    return np.linalg.norm(partial(func, xdata)(*params) - ydata)**2

bounds = [(-20, 20), (-20, 20)]

params_found = differential_evolution(param_func, bounds=bounds, args=(ydata,))['x']

它工作得很好，但在这个特定的情况下它完全失败了，我不知道为什么。

Answer 1

问题是您的函数非常扁平，只需尝试使用各种参数调用param_func，例如添加打印以查看差异进化遇到的情况。

params= [41.59913589  4.89105046] param_func= 0.03708384
params= [63.91254532 -0.92131526] param_func= 0.03708384
params= [68.37125624  0.89639623] param_func= 0.03708384
params= [45.14573252  4.65577815] param_func= 0.03708384
params= [22.67673743 -0.18044346] param_func= 0.03708384
params= [97.72482339 -2.55722856] param_func= 0.03708384
params= [52.78006606  2.19661847] param_func= 0.03708384
params= [93.83834497  1.06309298] param_func= 0.03708384
params= [77.03260669  3.93668517] param_func= 0.03708384
params= [71.58743947 -3.15837317] param_func= 0.03708384
params= [92.41338486 -1.26579432] param_func= 0.03708384
params= [12.5600005  -1.45394569] param_func= 0.03708384
params= [85.19070035 -0.4268145 ] param_func= 0.03708384
params= [59.43444477  1.35851078] param_func= 0.03708384
params= [60.46778441 -3.82879569] param_func= 0.03708384
params= [23.51634608 -3.59756272] param_func= 0.03708384
params= [38.7543976   3.12218973] param_func= 0.03708384
params= [28.67896624  2.45580586] param_func= 0.03708384
params= [55.4596277   4.19064059] param_func= 0.03708384
params= [76.25562633 -1.78696273] param_func= 0.03708384
params= [17.45482587  0.40538534] param_func= 0.03708384
params= [15.07905094 -2.72161755] param_func= 0.03708384
params= [1.95968808 2.81329829] param_func= 0.03708384
params= [81.5581296  -2.17837356] param_func= 0.03708384
params= [47.88669722 -4.446492  ] param_func= 0.03708384
params= [7.08142657 0.03370213] param_func= 0.03708384
params= [88.3735352   1.88745077] param_func= 0.03708384
params= [30.65315351 -4.81303647] param_func= 0.03708384
params= [6.32994012 3.63894112] param_func= 0.03708384
params= [36.36195255 -4.16279484] param_func= 0.03708384
params= [41.59913589  4.03415129] param_func= 0.03708384
params= [23.21663858  3.1233522 ] param_func= 0.03708384
params= [68.50852913  0.17748603] param_func= 0.03708384
params= [34.12052892 -1.33749366] param_func= 0.03708384
params= [17.043581    4.36269403] param_func= 0.03708384
params= [97.72482339  2.18998272] param_func= 0.03708384
params= [43.40457859  1.88110669] param_func= 0.03708384
params= [10.69430062  1.06309298] param_func= 0.03708384
params= [31.40159659 -2.21271748] param_func= 0.03708384
params= [32.22739184 -1.70684696] param_func= 0.03708384
params= [75.13605335  1.20226297] param_func= 0.03708384
params= [32.57892338  2.33051894] param_func= 0.03708384
params= [49.927827    0.71320109] param_func= 0.03708384
params= [25.80448259  4.92162336] param_func= 0.03708384
params= [6.29401507 4.19227187] param_func= 0.03708384
params= [54.33720709  4.17796811] param_func= 0.03708384
params= [54.79444146 -2.00878436] param_func= 0.03708384
params= [82.92651292  4.03998162] param_func= 0.03708384
params= [55.4596277   2.05457226] param_func= 0.03708384
params= [0.50126545 2.90680459] param_func= 0.03708384
params= [49.63974499 -1.81133026] param_func= 0.03708384
params= [18.43226804  3.58469356] param_func= 0.03708384
params= [55.990353   1.5163623] param_func= 0.03708384
params= [40.18236935 -1.11509064] param_func= 0.03708384
params= [58.62511485  1.60337965] param_func= 0.03708384
params= [7.08142657 0.91677437] param_func= 0.03708384
params= [50.96674056  2.33974317] param_func= 0.03708384
params= [71.71883967  2.55729519] param_func= 0.03708384
params= [88.58496221  2.4634768 ] param_func= 0.03708384
params= [46.03618557  3.4831179 ] param_func= 0.03708384
params= [41.59913589  4.03415129] param_func= 0.03708384
params= [41.5991359   4.03415129] param_func= 0.03708384
params= [41.59913589  4.0341513 ] param_func= 0.03708384

所以从任何“合理”优化的角度来看，您无需做任何事情，您的功能是不变的。您可以调整表示、停止条件等，但我会从重新考虑您的问题/数据陈述开始。

特别是

partial(func, xdata)(*params)

对于几乎每组参数来说都是一个 0 的向量，因此它可能不会做你想做的事情。

让我们添加一些绘图

param1 = np.linspace(bounds[0][0], bounds[0][1], 100)
param2 = np.linspace(bounds[1][0], bounds[1][1], 100)

xs = []
ys = []
vs = []
for p1 in param1:
    for p2 in param2:
        xs.append(p1)
        ys.append(p2)
        vs.append(param_func([p1, p2], ydata))

from matplotlib import pyplot as plt
plt.xlabel("param 1")
plt.ylabel("param 2")
plt.scatter(xs, ys, c=vs, linewidth=0, alpha=0.5)
plt.colorbar()
plt.show()

如您所见，除了边缘之外，该函数基本上是完全平坦的。如果我们将第二个参数限制在边缘，进化会按预期找到最小值。