如何为“scipy.interpolate.make_lsq_spline”选择好的节点序列

e7arh2l6 于 2023-01-17 发布在其他

关注(0)|答案(1)|浏览(131)

我想创建一个B样条平滑使用scipy.interpolate.make_lsq_spline的2D数据序列。

x = [0., 0.37427465, 0.68290943, 0.83261929, 1. ]
y = [-1.0, 3.0, 4.0, 2.0, 1.0]

但是，我不知道如何选择正确的t，错误消息对我来说没有意义。

In [1]: import numpy as np

In [2]: from scipy.interpolate import make_lsq_spline

In [3]: x = [0., 0.37427465, 0.68290943, 0.83261929, 1. ]

In [4]: y = [-1.0, 3.0, 4.0, 2.0, 1.0]

In [5]: t = [0.,0.,0.,0.,0.25,0.5,0.75,1.,1.,1.,1 ]

In [6]: spl = make_lsq_spline(x, y, t)
---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
<ipython-input-6-4440a73d26f0> in <cell line: 1>()
----> 1 spl = make_lsq_spline(x, y, t)

/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/scipy/interpolate/_bsplines.py in make_lsq_spline(x, y, t, k, w, axis, check_finite)
   1513
   1514     # have observation matrix & rhs, can solve the LSQ problem
-> 1515     cho_decomp = cholesky_banded(ab, overwrite_ab=True, lower=lower,
   1516                                  check_finite=check_finite)
   1517     c = cho_solve_banded((cho_decomp, lower), rhs, overwrite_b=True,

/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/scipy/linalg/_decomp_cholesky.py in cholesky_banded(ab, overwrite_ab, lower, check_finite)
    280     c, info = pbtrf(ab, lower=lower, overwrite_ab=overwrite_ab)
    281     if info > 0:
--> 282         raise LinAlgError("%d-th leading minor not positive definite" % info)
    283     if info < 0:
    284         raise ValueError('illegal value in %d-th argument of internal pbtrf'

LinAlgError: 5-th leading minor not positive definite

是否有选择正确打结顺序t的指南？

scipy

来源：https://stackoverflow.com/questions/73592565/how-to-select-good-knot-sequences-for-scipy-interpolate-make-lsq-spline

1条答案

按热度按时间

a1o7rhls1#

我有一个类似的问题。由于你的例子，我想我可以告诉哪里出了问题。从线性代数的Angular 来看，你要求一个不能唯一解的问题的解，你提供11个节点t，这意味着有11-3-1 = 7个系数要确定，因为你试图用k=3次样条拟合（make_lsq_spline的默认值）。在5个点x上求值，方程组的左边由一个5 × 7的矩阵D给出。2在这个例子中D是满秩的，但这没有帮助。3 7 × 7的矩阵N = D。4 T@D只是半正定的。5两个特征值是0。它不能被反转，因此你的问题不能被唯一地解决。一个解决方案是去掉2个结点，比如0.25和0.75处的结点。当处理3次样条时，你应该保留边界处的四折结点，因为你很可能希望你的插值样条跳到那里。总之，节点必须以这样一种方式来选择，即插值问题是唯一可解的。2我也试着添加一些代码来说明我想说的。3希望能有所帮助。

import numpy as np
import scipy.interpolate as sciint
import matplotlib.pyplot as plt

x = [0., 0.37427465, 0.68290943, 0.83261929, 1.]
y = [-1.0, 3.0, 4.0, 2.0, 1.0]
t = [0.,0.,0.,0.,0.25,0.5,0.75,1.,1.,1.,1 ]

splines = []

for k in range(7):
    coeff    = np.zeros(7)
    coeff[k] = 1.
    splines.append(sciint.BSpline(t,coeff,3))

fig,ax = plt.subplots(3,3)
dom    = np.linspace(0.,1.,1000)

for count,axes in enumerate(ax.flat):
    axes.plot(dom,splines[count](dom))
    if count == len(splines)-1:
        break
    
data = []

for spline in splines:
    data.append(np.vstack(spline(x)))
    
D = np.hstack(tuple(data))    
N = D.T @ D

sing = np.linalg.eigvalsh(N)

print(sing)

t2 = [0.,0.,0.,0.,0.5,1.,1.,1.,1 ]

bspline = sciint.make_lsq_spline(x,y,t2)

ax[2,1].plot(x,y,'r')
ax[2,2].plot(dom,bspline(dom),'m')

赞(0）回复(0）举报 2023-01-17

我来回答

如何为“scipy.interpolate.make_lsq_spline”选择好的节点序列

1条答案

相关问题

热门标签

最新问答