scipy fsolve和root哪个更快？

Question

我正在尝试找到一个函数的根。过去我使用过fsolve，但随着我的数据集变得越来越大，它似乎变得不太一致 (--> n = 187)。现在我正在寻找替代方案，并发现了scipy.root。我不明白两者之间的区别，以及在我的情况下哪一个更好。

我正试图解决以下3N个耦合方程，并尝试找到向量x y和z：

我的代码如下，其中inrec、outrec和rec是预先确定的值列表：

import scipy as sp
import numpy as np
from scipy.optimize import fsolve 
import math

def f(w, n, onrec, inrec, rec):
    F = [0]*3*n 
    for i in range(n): 
        F[i] = -onrec[i] #k_i>
        F[n+i] = -inrec[i] #k_i<
        F[(2*n)+i] = -rec[i] #k_i <>
        for j in range(n):
            if i == j:
                None
            else: #below the three functions are stated. w[i] = x_i, w[n+i] = y_i, w[2*n + i] = z_i
                F[i] += (w[i]*w[n+j])/(1+w[i]*w[n+j]+w[j]*w[n+i]+w[2*n+i]*w[2*n+j])
                F[n+i] += (w[j]*w[n+i])/(1+w[i]*w[n+j]+w[j]*w[n+i]+w[2*n+i]*w[2*n+j])
                F[2*n+i] += (w[(2*n)+i]*w[(2*n)+j])/(1+w[i]*w[n+j]+w[j]*w[n+i]+w[2*n+i]*w[2*n+j])
    return(F)
    
u = [1]*3*n
s = fsolve(f, u, args=(n, onrec, inrec, rec))

- Thomas Mc Donald

---

我不会尝试回答，因为我没有数据，但我猜问题在于您正在使用两个嵌套循环，函数的评估是O(N^2)，但您可以通过向量化大大降低常数因子。 - Bob

Answer 1

正如@Bob所建议的那样，第一步必须是将您的内部函数向量化。之后，回到您的主要问题：这不是正确的问题，因为

1. fsolve只是hybr算法的包装器，而该算法已作为root的一个选项提供；以及
2. 在性能之前担心正确性。

几乎可以肯定的是，优化器放弃了您的问题，结果是无效的。我唯一能够说服它收敛的情况是n=4和Levenberg-Marquardt算法。如果（四年后）您仍然需要解决此问题，我建议将问题带到像Mathematics StackExchange这样的不同社区。但与此同时，这里有一个具有一个收敛解的向量化示例：

import numpy as np
from numpy.random import default_rng
from scipy.optimize import root


def lhs(xyz: np.ndarray) -> np.ndarray:
    x, y, z = xyz[..., np.newaxis]
    n = len(x)

    coeff = (1 - np.eye(n)) / (1 + x*y.T + x.T*y + z*z.T)
    numerator = np.stack((
        x * y.T,
        x.T * y,
        z.T * z,
    ))
    result = (numerator * coeff).sum(axis=2)
    return result


def to_root(w: np.ndarray, recs: np.ndarray) -> np.ndarray:
    xyz = w.reshape((3, -1))
    rhs = lhs(xyz) - recs
    return rhs.ravel()


def test_data(n: int = 4) -> tuple[np.ndarray, np.ndarray]:
    rand = default_rng(seed=0)
    secret_solution = rand.uniform(-1, 1, (3, n))
    recs = lhs(secret_solution)
    return secret_solution, recs


def method_search() -> None:
    secret_solution, recs = test_data()

    for method in ('hybr', 'lm', 'broyden1', 'broyden2', 'anderson',
                   'linearmixing', 'diagbroyden', 'excitingmixing',
                   'krylov', 'df-sane'):
        try:
            result = root(
                to_root, x0=np.ones_like(recs),
                args=(recs,), method=method,
                options={
                    'maxiter': 5_000,
                    'maxfev': 5_000,
                },
            )
        except Exception:
            continue

        print(method, result.message,
              f'nfev={getattr(result, "nfev", None)} nit={getattr(result, "nit", None)}')
        print('Estimated RHS:')
        print(lhs(result.x.reshape((3, -1))))
        print('Estimated error:')
        print(to_root(result.x, recs))
        print()


def successful_example() -> None:
    n = 4
    print('n=', n)
    secret_solution, recs = test_data(n=n)

    result = root(
        to_root, x0=np.ones_like(recs),
        args=(recs,), method='lm',
    )
    print(result.message)
    print('function evaluations:', result.nfev)
    error = to_root(result.x, recs)
    print('Error:', error.dot(error))
    print()


if __name__ == '__main__':
    successful_example()

n= 4
The relative error between two consecutive iterates is at most 0.000000
function evaluations: 1221
Error: 8.721381160163159e-30

- Reinderien