寻找求解带唯一正解非线性方程的快速Python优化算法

问题目标 #

寻找一种快速Python算法，求解下述函数f(x)的正解：

def f(x):
    return (l / ((np.tile(r, (n, 1)).transpose() / D / (np.tile(x, (n, 1)) / D).sum(axis = 0)).sum(axis = 1))) - x

其中l、r、x为n维向量，D为n×n矩阵。已知：

• 方程存在正解
• 解在缩放因子下唯一（即若x是解，则k*x也是解，k>0）
• 需要支持n最大约4000的规模

已尝试方案 #

试过scipy.optimize下的多个函数，但都存在问题：

• fsolve：偶尔会得到含负元素的解，不符合正解要求
• minimize（带正约束）：只有当初始值非常接近真实解时才能找到全局最优，初始值稍有偏差就会陷入局部最优，无法得到准确解
• differential_evolution：能找到正确解，但n≥250时速度极慢，仅适合n很小的测试场景

合适算法推荐与加速技巧 #

1. 带边界约束的拟牛顿法（L-BFGS-B） #

使用scipy.optimize.minimize的L-BFGS-B方法，它天然支持变量上下界约束，直接设置bounds=[(1e-8, None)]*n（用极小正数代替0避免数值问题）。同时结合归一化约束减少自由度：因为解缩放后唯一，可固定x的几何均值为1，把n维问题降为n-1维，避免算法在无效的缩放方向浪费计算资源。

2. 信任域约束优化（trust-constr） #

scipy.optimize.minimize的trust-constr方法适合带约束的非线性优化，能处理等式/不等式约束，若D是稀疏矩阵还可利用稀疏性优化，对于n=4000的规模也能保持较好效率。

3. 自定义不动点迭代法 #

观察原函数f(x)=0可变形为不动点迭代格式：
$$x = l / \left( \sum_j \frac{r_j / D_{ji}}{\sum_k x_k / D_{kj}} \right)$$
构造正初始值（比如全1向量、l/r向量），迭代更新x后每次归一化保持几何均值为1，直到收敛。该方法实现简单、计算效率高，适合大规模n的场景。

关键加速细节 #

• 利用缩放唯一性降维：每次迭代后将x归一化（如几何均值为1），减少无效搜索方向
• 优化函数计算效率：用numpy广播代替原代码中的np.tile，减少内存占用与计算时间，优化后的函数如下：

def f(x):
    x_over_D = x / D  # shape (n,n)
    sum_xD = x_over_D.sum(axis=0)  # shape (n,)
    rD_over_sum = (r[:, np.newaxis] / D) / sum_xD  # shape (n,n)
    denom = rD_over_sum.sum(axis=1)  # shape (n,)
    return l / denom - x

• 选择合理初始值：用全1向量或l/r作为初始值，避免初始值偏离最优解过远，减少迭代次数

最小工作示例 #

import numpy as np
from scipy.optimize import minimize

np.random.seed(1)
n = 250  # 测试可改用n=5，大规模测试用n=4000

# 生成测试数据
r = np.random.rand(n)
D = 1 + np.random.rand(n, n)
x_true = np.random.rand(n)
# 归一化真实解，几何均值为1
x_true = x_true / np.prod(x_true) ** (1/n)
# 计算对应的l
x_over_D_true = x_true / D
sum_xD_true = x_over_D_true.sum(axis=0)
rD_over_sum_true = (r[:, np.newaxis] / D) / sum_xD_true
denom_true = rD_over_sum_true.sum(axis=1)
l = denom_true * x_true

# 优化后的目标函数（最小化平方和）
def opt(x):
    x_over_D = x / D
    sum_xD = x_over_D.sum(axis=0)
    rD_over_sum = (r[:, np.newaxis] / D) / sum_xD
    denom = rD_over_sum.sum(axis=1)
    return np.sum((l / denom - x)**2)

# 初始值选择与归一化
x0 = np.ones(n)
x0 = x0 / np.prod(x0)**(1/n)

# 使用L-BFGS-B方法求解
result = minimize(
    opt,
    x0=x0,
    method='L-BFGS-B',
    bounds=[(1e-8, None)]*n,
    tol=1e-10
)

# 归一化结果并验证
x_sol = result.x / np.prod(result.x)**(1/n)
print(f"L-BFGS-B平均误差：{abs(x_sol - x_true).mean():.6e}")
print(f"L-BFGS-B目标函数值：{result.fun:.6e}")

# 不动点迭代求解示例
def fixed_point_iter(x_init, max_iter=1000, tol=1e-10):
    x = x_init.copy()
    for i in range(max_iter):
        x_over_D = x / D
        sum_xD = x_over_D.sum(axis=0)
        rD_over_sum = (r[:, np.newaxis] / D) / sum_xD
        denom = rD_over_sum.sum(axis=1)
        x_new = l / denom
        # 归一化
        x_new = x_new / np.prod(x_new)**(1/n)
        # 检查收敛
        if np.max(abs(x_new - x)) < tol:
            print(f"不动点迭代收敛，迭代次数：{i+1}")
            return x_new
        x = x_new
    print("不动点迭代未收敛")
    return x

x_fp = fixed_point_iter(x0)
print(f"不动点迭代平均误差：{abs(x_fp - x_true).mean():.6e}")

内容的提问来源于stack exchange，提问作者eigenvector