numpy 拟合系数多项式的线性组合

Question

我有要用多项式拟合的数据。我有 200,000 个数据点，所以我想要一个高效的算法。我想使用 numpy.polynomial 包，这样我就可以尝试不同的多项式族和次数。有什么方法可以将其表述为像 Ax=b 这样的方程组吗？有没有比 scipy.minimize 更好的方法来解决这个问题？

import numpy as np
from scipy.optimize import minimize as mini

x1 = np.random.random(2000)
x2 = np.random.random(2000)
y = 20 * np.sin(x1) + x2 - np.sin (30 * x1 - x2 / 10)
def fitness(x, degree=5):

    poly1 = np.polynomial.polynomial.polyval(x1, x[:degree])
    poly2 = np.polynomial.polynomial.polyval(x2, x[degree:])
    return np.sum((y - (poly1 + poly2)) ** 2 )

# It seems like I should be able to solve this as a system of equations
# x = np.linalg.solve(np.concatenate([x1, x2]), y)

# minimize the sum of the squared residuals to find the optimal polynomial coefficients
x = mini(fitness, np.ones(10))
print fitness(x.x)

Answer 1

你的直觉是对的。您可以将其求解为 Ax = b.

形式的方程组

但是：

系统超定义，想得到最小二乘解，所以需要用np.linalg.lstsq代替np.linalg.solve。
~~不能用polyval因为需要把自变量的系数和幂分开~~

这是构造方程组并求解的方法：

A = np.stack([x1**0, x1**1, x1**2, x1**3, x1**4, x2**0, x2**1, x2**2, x2**3, x2**4]).T             
xx = np.linalg.lstsq(A, y)[0]
print(fitness(xx))  # test the result with original fitness function

当然你可以对学位进行概括：

A = np.stack([x1**p for p in range(degree)] + [x2**p for p in range(degree)]).T

对于示例数据，最小二乘法的运行速度比 minimize 快得多（在我的笔记本电脑上为 800µs 与 35ms）。但是，A 可能会变得非常大，因此如果内存有问题，minimize 可能仍然是一个选项。

更新：

如果不了解多项式函数的内部结构，事情就会变得棘手，但可以将项和系数分开。这是从 polyval:

这样的函数构造系统矩阵 A 的一种有点丑陋的方法

def construct_A(valfunc, degree):
    columns1 = []
    columns2 = []
    for p in range(degree):
        c = np.zeros(degree)
        c[p] = 1
        columns1.append(valfunc(x1, c))
        columns2.append(valfunc(x2, c))
    return np.stack(columns1 + columns2).T

A = construct_A(np.polynomial.polynomial.polyval, 5)           
xx = np.linalg.lstsq(A, y)[0]
print(fitness(xx))  # test the result with original fitness function

numpy 拟合系数多项式的线性组合

numpy fit coefficients to linear combination of polynomials

python

math

numpy

linear-algebra

polynomials