什么是scipy fsolve vs root？

Question

我试图找到函数的根。我过去使用过fsolve，但是随着我的数据集越来越大，它似乎变得更加不一致（-> n = 187）。现在我正在寻找替代品并找到了scipy.root。我不明白两者之间的区别是什么，在我的情况下哪个更好。

我正在尝试求解以下3N耦合方程，并试图找到向量xy和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))

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()

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()