Python 中二维多项式的 `polyfit` 等价物

Question

我想找到

 中 

a

系数的最小二乘解

z = (a0 + a1*x + a2*y + a3*x**2 + a4*x**2*y + a5*x**2*y**2 + a6*y**2 +
     a7*x*y**2 + a8*x*y)

给定长度为 20 的数组

x

、

y

和

z

。基本上我正在寻找

numpy.polyfit

的等价物，但是对于 2D 多项式。

这个问题类似，但解决方案是通过MATLAB提供的。

Answer 1

以下示例展示了如何使用

numpy.linalg.lstsq

来完成此任务：

import numpy as np

x = np.linspace(0, 1, 20)
y = np.linspace(0, 1, 20)
X, Y = np.meshgrid(x, y, copy=False)
Z = X**2 + Y**2 + np.random.rand(*X.shape)*0.01

X = X.flatten()
Y = Y.flatten()

A = np.array([X*0+1, X, Y, X**2, X**2*Y, X**2*Y**2, Y**2, X*Y**2, X*Y]).T
B = Z.flatten()

coeff, r, rank, s = np.linalg.lstsq(A, B)

调整系数

coeff

为：

array([ 0.00423365,  0.00224748,  0.00193344,  0.9982576 , -0.00594063,
        0.00834339,  0.99803901, -0.00536561,  0.00286598])

请注意，

coeff[3]

和

coeff[6]

分别对应于

X**2

和

Y**2

，并且它们接近

1.

，因为示例数据是使用

Z = X**2 + Y**2 + small_random_component

创建的。

Answer 2

根据@Saullo 和@Francisco 的答案，我做了一个我发现有用的函数：

def polyfit2d(x, y, z, kx=3, ky=3, order=None):
    '''
    Two dimensional polynomial fitting by least squares.
    Fits the functional form f(x,y) = z.

    Notes
    -----
    Resultant fit can be plotted with:
    np.polynomial.polynomial.polygrid2d(x, y, soln.reshape((kx+1, ky+1)))

    Parameters
    ----------
    x, y: array-like, 1d
        x and y coordinates.
    z: np.ndarray, 2d
        Surface to fit.
    kx, ky: int, default is 3
        Polynomial order in x and y, respectively.
    order: int or None, default is None
        If None, all coefficients up to maxiumum kx, ky, ie. up to and including x^kx*y^ky, are considered.
        If int, coefficients up to a maximum of kx+ky <= order are considered.

    Returns
    -------
    Return paramters from np.linalg.lstsq.

    soln: np.ndarray
        Array of polynomial coefficients.
    residuals: np.ndarray
    rank: int
    s: np.ndarray

    '''

    # grid coords
    x, y = np.meshgrid(x, y)
    # coefficient array, up to x^kx, y^ky
    coeffs = np.ones((kx+1, ky+1))

    # solve array
    a = np.zeros((coeffs.size, x.size))

    # for each coefficient produce array x^i, y^j
    for index, (j, i) in enumerate(np.ndindex(coeffs.shape)):
        # do not include powers greater than order
        if order is not None and i + j > order:
            arr = np.zeros_like(x)
        else:
            arr = coeffs[i, j] * x**i * y**j
        a[index] = arr.ravel()

    # do leastsq fitting and return leastsq result
    return np.linalg.lstsq(a.T, np.ravel(z), rcond=None)

最终的拟合结果可以通过以下方式可视化：

fitted_surf = np.polynomial.polynomial.polyval2d(x, y, soln.reshape((kx+1,ky+1)))
plt.matshow(fitted_surf)

Answer 3

Saullo Castro 的出色回答。只需添加代码以使用 a 系数的最小二乘解来重建函数，

def poly2Dreco(X, Y, c):
    return (c[0] + X*c[1] + Y*c[2] + X**2*c[3] + X**2*Y*c[4] + X**2*Y**2*c[5] + 
           Y**2*c[6] + X*Y**2*c[7] + X*Y*c[8])

Answer 4

您也可以使用 scikit-learn 来实现此目的。

import numpy as np
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression

x = np.linspace(0, 1, 20)
y = np.linspace(0, 1, 20)
X, Y = np.meshgrid(x, y, copy=False)
X = X.flatten()
Y = Y.flatten()

# Generate noisy data
np.random.seed(0)
Z = X**2 + Y**2 + np.random.randn(*X.shape)*0.01

# Process 2D inputs
poly = PolynomialFeatures(degree=2)
input_pts = np.stack([X, Y]).T
assert(input_pts.shape == (400, 2))
in_features = poly.fit_transform(input_pts)

# Linear regression
model = LinearRegression(fit_intercept=False)
model.fit(in_features, Z)

# Display coefficients
print(dict(zip(poly.get_feature_names_out(), model.coef_.round(4))))

# Check fit
print(f"R-squared: {model.score(poly.transform(input_pts), Z):.3f}")

# Make predictions
Z_predicted = model.predict(poly.transform(input_pts))

出：

{'1': 0.0012, 'x0': 0.003, 'x1': -0.0074, 'x0^2': 0.9974, 'x0 x1': 0.0047, 'x1^2': 1.0014}
R-squared: 1.000

注释

我在定义线性回归模型时使用了

fit_intercept=False

参数，因为默认情况下多项式特征包含偏差项

'1'

。或者，您可以使用

禁用它

poly = PolynomialFeatures(degree=2, include_bias=False)

然后使用带有截距系数的常规

LinearRegression

模型：

model.intercept_  # 0.0011741137800872492

Answer 5

请注意，如果

kx != ky

代码将会失败，因为

j

和

i

索引在循环中反转。

您从

(j,i)

得到

enumerate(np.ndindex(coeffs.shape))

，但随后您将

coeffs

中的元素寻址为

coeffs[i,j]

。由于系数矩阵的形状是由您要求使用的最大多项式阶数给出的，因此如果

kx != ky

并且您将超过其维度之一，则该矩阵将是矩形。