Python中旋转平移双曲线点集拟合的代码排查与实现

云上老周

2026-8-28

问题描述

我需要对xy平面内的一组数据点拟合通用形式的、经旋转和平移变换的双曲线，以此反推二次曲线一般方程的系数。我尝试了已有的直接最小二乘拟合双曲线/椭圆的方法，但代码始终无法正常运行：对已知属于双曲线的点集做拟合时，输出结果和真实系数偏差极大。需要定位代码错误，并提供可行的替代方案。

原始复现代码

import numpy as np
from sympy import plot_implicit, Eq
from sympy.abc import x, y

def fit_hyperbola(x, y):   
    D1 = np.vstack([x**2, x*y, y**2]).T
    D2 = np.vstack([x, y, np.ones(len(x))]).T
    S1 = D1.T @ D1
    S2 = D1.T @ D2
    S3 = D2.T @ D2

    # define the constraint matrix and its inverse
    C = np.array(((0, 0, -2), (0, 1, 0), (-2, 0, 0)), dtype=float)
    Ci = np.linalg.inv(C)

    # Setup and solve the generalized eigenvector problem
    T = np.linalg.inv(S3) @ S2.T

    S = Ci@(S1 - S2@T)

    eigval, eigvec = np.linalg.eig(S)

    # evaluate and sort resulting constraint values
    cond = eigvec[1]**2 - 4*eigvec[0]*eigvec[2]

    # [condVals index] = sort(cond)
    idx = np.argsort(cond)
    condVals = cond[idx]

    possibleHs = condVals[1:] + condVals[0]
    minDiffAt = np.argmin(abs(possibleHs))
    # minDiffVal = possibleHs[minDiffAt]
    alpha1 = eigvec[:, idx[minDiffAt + 1]]
    alpha2 = T@alpha1
    return np.concatenate((alpha1, alpha2)).ravel()

if __name__ == '__main__':

    # known hyperbola coefficients
    coeffs = [1., 6., -2., 3., 0., 0.]

    # hyperbola points
    x_ = [1.56011303e+00, 1.38439984e+00, 1.22595618e+00, 1.08313085e+00,
           9.54435408e-01, 8.38528681e-01, 7.34202759e-01, 6.40370424e-01,
           5.56053814e-01, 4.80374235e-01, 4.12543002e-01, 3.51853222e-01,
           2.97672424e-01, 2.49435970e-01, 2.06641170e-01, 1.68842044e-01,
           1.35644673e-01, 1.06703097e-01, 8.17157025e-02, 6.04220884e-02,
           4.26003457e-02, 2.80647476e-02, 1.66638132e-02, 8.27872926e-03,
           2.82211172e-03, 2.37095181e-04, 4.96740239e-04, 3.60375275e-03,
           9.59051203e-03, 1.85194083e-02, 3.04834928e-02, 4.56074477e-02,
           6.40488853e-02, 8.59999904e-02, 1.11689524e-01, 1.41385205e-01,
           1.75396504e-01, 2.14077865e-01, 2.57832401e-01, 3.07116093e-01,
           3.62442545e-01, 4.24388335e-01, 4.93599021e-01, 5.70795874e-01,
           6.56783391e-01, 7.52457678e-01, 8.58815793e-01, 9.76966133e-01,
           1.10813998e+00, 1.25370436e+00]

    y_ =  [-0.66541515, -0.6339625 , -0.60485332, -0.57778425, -0.5524732 ,
           -0.52865638, -0.50608561, -0.48452564, -0.46375182, -0.44354763,
           -0.42370253, -0.4040097 , -0.38426392, -0.3642594 , -0.34378769,
           -0.32263542, -0.30058217, -0.27739811, -0.25284163, -0.22665682,
           -0.19857079, -0.16829086, -0.13550147, -0.0998609 , -0.06099773,
           -0.01850695,  0.02805425,  0.07917109,  0.13537629,  0.19725559,
            0.26545384,  0.34068177,  0.42372336,  0.51544401,  0.61679957,
            0.72884632,  0.85275192,  0.98980766,  1.14144182,  1.30923466,
            1.49493479,  1.70047747,  1.92800474,  2.17988774,  2.45875143,
            2.76750196,  3.10935692,  3.48787892,  3.90701266,  4.3711261 ]


    plot_implicit (Eq(coeffs[0]*x**2 + coeffs[1]*x*y + coeffs[2]*y**2 + coeffs[3]*x + coeffs[4]*y, -coeffs[5]))

    coeffs_fit = fit_hyperbola(x_, y_)

    plot_implicit (Eq(coeffs_fit[0]*x**2 + coeffs_fit[1]*x*y + coeffs_fit[2]*y**2 + coeffs_fit[3]*x + coeffs_fit[4]*y, -coeffs_fit[5]))

代码错误定位

代码共有3处核心错误：

约束矩阵符号完全写反：当前使用的约束矩阵对应椭圆拟合的判别式约束b²-4ac < 0，双曲线拟合需要固定二次项判别式b²-4ac = 1（避免零解的正约束），对应约束矩阵应为np.array(((0, 0, 2), (0, -1, 0), (2, 0, 0)), dtype=float)。
特征向量筛选逻辑错误：原代码的排序、加和找最小值的逻辑没有理论依据，广义特征值分解后，3个特征向量里仅有1个满足双曲线判别式b²-4ac>0，直接筛选该特征向量即可，不需要额外排序计算。
数值稳定性缺陷：直接对矩阵求逆np.linalg.inv(S3)在点集坐标量级差异较大时会引入显著数值误差，应使用最小二乘求解代替显式求逆。

修正后代码

import numpy as np
from sympy import plot_implicit, Eq
from sympy.abc import x, y

def fit_hyperbola(x, y):   
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)
    D1 = np.vstack([x**2, x*y, y**2]).T
    D2 = np.vstack([x, y, np.ones(len(x))]).T
    S1 = D1.T @ D1
    S2 = D1.T @ D2
    S3 = D2.T @ D2

    # 双曲线约束：b²-4ac = 1 对应的约束矩阵
    C = np.array(((0, 0, 2), (0, -1, 0), (2, 0, 0)), dtype=float)
    Ci = np.linalg.inv(C)

    # 用最小二乘代替显式求逆提升稳定性
    T = np.linalg.lstsq(S3, S2.T, rcond=None)[0]
    S = Ci @ (S1 - S2 @ T)

    eigval, eigvec = np.linalg.eig(S)
    # 筛选满足双曲线判别式为正的特征向量
    for i in range(3):
        a, b, c = eigvec[:, i]
        disc = b**2 - 4*a*c
        if disc > 1e-8: # 排除数值误差引入的近零值
            alpha1 = eigvec[:, i]
            break
    alpha2 = T @ alpha1
    coeff = np.concatenate((alpha1, alpha2)).ravel()
    # 归一化系数方便和真实值对比
    coeff = coeff / coeff[0]
    return coeff

if __name__ == '__main__':
    coeffs = [1., 6., -2., 3., 0., 0.]
    x_ = [1.56011303e+00, 1.38439984e+00, 1.22595618e+00, 1.08313085e+00,
           9.54435408e-01, 8.38528681e-01, 7.34202759e-01, 6.40370424e-01,
           5.56053814e-01, 4.80374235e-01, 4.12543002e-01, 3.51853222e-01,
           2.97672424e-01, 2.49435970e-01, 2.06641170e-01, 1.68842044e-01,
           1.35644673e-01, 1.06703097e-01, 8.17157025e-02, 6.04220884e-02,
           4.26003457e-02, 2.80647476e-02, 1.66638132e-02, 8.27872926e-03,
           2.82211172e-03, 2.37095181e-04, 4.96740239e-04, 3.60375275e-03,
           9.59051203e-03, 1.85194083e-02, 3.04834928e-02, 4.56074477e-02,
           6.40488853e-02, 8.59999904e-02, 1.11689524e-01, 1.41385205e-01,
           1.75396504e-01, 2.14077865e-01, 2.57832401e-01, 3.07116093e-01,
           3.62442545e-01, 4.24388335e-01, 4.93599021e-01, 5.70795874e-01,
           6.56783391e-01, 7.52457678e-01, 8.58815793e-01, 9.76966133e-01,
           1.10813998e+00, 1.25370436e+00]
    y_ =  [-0.66541515, -0.6339625 , -0.60485332, -0.57778425, -0.5524732 ,
           -0.52865638, -0.50608561, -0.48452564, -0.46375182, -0.44354763,
           -0.42370253, -0.4040097 , -0.38426392, -0.3642594 , -0.34378769,
           -0.32263542, -0.30058217, -0.27739811, -0.25284163, -0.22665682,
           -0.19857079, -0.16829086, -0.13550147, -0.0998609 , -0.06099773,
           -0.01850695,  0.02805425,  0.07917109,  0.13537629,  0.19725559,
            0.26545384,  0.34068177,  0.42372336,  0.51544401,  0.61679957,