numpy 数组的第一个最小或鞍点（微积分导数）是多少？
What is the first minimum OR saddle point (calculus derivative) of a numpy array?
我有以下 numpy 数组，如上所示。
类似
的函数
print(arr.argsort()[:3])
将return三个最低值的三个最低值：
[69 66 70]
我如何 return first index 其中第一个 minimum 或第一个 saddle point （在微积分意义上) 以先到数组为准？
在这种情况下，index 2 和 3 处的两个数字 0.62026396 0.60566623 是第一个 saddle point（这不是真正的 saddle point，因为斜率不t变平了，但它显然打破了那里的硬下坡。也许添加一个"flattens"意味着什么的阈值）。由于函数在第一个 saddle point 之前不会上升，因此第一个 mimimum 出现在 saddle point 之后，这是我感兴趣的索引。
[1.04814804 0.90445908 0.62026396 0.60566623 0.32295758 0.26658469
0.19059289 0.10281547 0.08582772 0.05091265 0.03391474 0.03844931
0.03315003 0.02838656 0.03420759 0.03567401 0.038203 0.03530763
0.04394316 0.03876966 0.04156067 0.03937291 0.03966426 0.04438747
0.03690863 0.0363976 0.03171374 0.03644719 0.02989291 0.03166156
0.0323875 0.03406287 0.03691943 0.02829374 0.0368121 0.02971704
0.03427005 0.02873735 0.02843848 0.02101889 0.02114978 0.02128403
0.0185619 0.01749904 0.01441699 0.02118773 0.02091855 0.02431763
0.02472427 0.03186318 0.03205664 0.03135686 0.02838413 0.03206674
0.02638371 0.02048122 0.01502128 0.0162665 0.01331485 0.01569286
0.00901017 0.01343558 0.00908635 0.00990869 0.01041151 0.01063606
0.00822482 0.01312368 0.0115005 0.00620334 0.0084177 0.01058152
0.01198732 0.01451455 0.01605602 0.01823713 0.01685975 0.03161889
0.0216687 0.03052391 0.02220871 0.02420951 0.01651778 0.02066987
0.01999613 0.02532265 0.02589186 0.02748692 0.02191687 0.02612152
0.02309497 0.02744753 0.02619196 0.02281516 0.0254296 0.02732746
0.02567608 0.0199178 0.01831929 0.01776025]
您可以使用 np.gradient 或 np.diff 来评估差异（第一个计算中心差异，第二个只是 x[1:] - x[:-1]），然后使用 np.sign 获取渐变符号，另一个 np.diff 查看符号变化的位置。然后过滤正号变化（对应最小值）：
np.where(np.diff(np.sign(np.gradient(x))) > 0)[0][0]+2 #add 2 as each time you call np.gradient or np.diff you are substracting 1 in size, the first [0] is to get the positions, the second [0] is to get the "first" element
>> 8
x[np.where(np.diff(np.sign(np.gradient(x))) > 0)[0][0]+2]
>> 0.03420759
环顾四周并根据给出的两个建议（到目前为止），我这样做了：
import scipy
from scipy import interpolate
x = np.arange(0, 100)
spl = scipy.interpolate.splrep(x,arr,k=3) # no smoothing, 3rd order spline
ddy = scipy.interpolate.splev(x,spl,der=2) # use those knots to get second derivative
print(ddy)
asign = np.sign(ddy)
signchange = ((np.roll(asign, 1) - asign) != 0).astype(int)
print(signchange)
这给了我 second derivative，然后我可以对其进行分析，例如，查看符号变化发生的位置：
[-0.894053 -0.14050616 0.61304067 -0.69407217 0.55458251 -0.16624336
-0.0073225 0.12481963 -0.067218 0.03648846 0.02876712 -0.02236204
0.00167794 0.01886512 -0.0136314 0.00953279 -0.01812436 0.03041855
-0.03436446 0.02418512 -0.01458896 0.00429809 0.01227133 -0.02679232
0.02168571 -0.0181437 0.02585209 -0.02876075 0.0214645 -0.00715966
0.0009179 0.00918466 -0.03056938 0.04419937 -0.0433638 0.03557532
-0.02904901 0.02010647 -0.0199739 0.0170648 -0.00298236 -0.00511529
0.00630525 -0.01015011 0.02218007 -0.01945341 0.01339405 -0.01211326
0.01710444 -0.01591092 0.00486652 -0.00891456 0.01715403 -0.01976949
0.00573004 -0.00446743 0.01479495 -0.01448144 0.01794968 -0.02533936
0.02904355 -0.02418628 0.01505374 -0.00499926 0.00302616 -0.00877499
0.01625907 -0.01240068 -0.00578862 0.01351128 -0.00318733 -0.0010652
0.0029 -0.0038062 0.0064102 -0.01799678 0.04422601 -0.0620881
0.05587037 -0.04856099 0.03535114 -0.03094757 0.03028399 -0.01912546
0.01726283 -0.01392421 0.00989012 -0.01948119 0.02504401 -0.02204667
0.0197554 -0.01270022 -0.00260326 0.01038581 -0.00299247 -0.00271539
-0.00744152 0.00784016 0.00103947 -0.00576122]
[0 0 1 1 1 1 0 1 1 1 0 1 1 0 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 0 1 1 1 1 1
1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 0 1 0 1]
这就是我检测局部 maxima/minima、拐点和鞍点的方式。
首先定义以下函数
import numpy as np
def n_derivative(arr, degree=1):
"""Compute the n-th derivative."""
result = arr.copy()
for i in range(degree):
result = np.gradient(result)
return result
def sign_change(arr):
"""Detect sign changes."""
sign = np.sign(arr)
result = ((np.roll(sign, 1) - sign) != 0).astype(bool)
result[0] = False
return result
def zeroes(arr, threshold=1e-8):
"""Find zeroes of an array."""
return sign_change(arr) | (abs(arr) < threshold)
我们现在可以使用 derivative test
关键点的一阶导数将为零。
def critical_points(arr):
return zeroes(n_derivative(arr, 1))
如果临界点的二阶导数不为零，则该点要么是最大值，要么是最小值：
def maxima_minima(arr):
return zeroes(n_derivative(arr, 1)) & ~zeroes(n_derivative(arr, 2))
def maxima(arr):
return zeroes(n_derivative(arr, 1)) & (n_derivative(arr, 2) < 0)
def minima(arr):
return zeroes(n_derivative(arr, 1)) & (n_derivative(arr, 2) > 0)
如果二阶导数为零但三阶导数不为零，则该点为拐点：
def inflections(arr):
return zeroes(n_derivative(arr, 2)) & ~zeroes(n_derivative(arr, 3))
如果临界点的二阶导数为零，但三阶导数不为零，则这是一个鞍点：
def inflections(arr):
return zeroes(n_derivative(arr, 1)) & zeroes(n_derivative(arr, 2)) & ~zeroes(n_derivative(arr, 3))
请注意，此方法在数值上不稳定，因为一方面在某个任意阈值定义上检测到零，另一方面不同的采样可能导致函数/数组不可微分。
因此，根据这个定义，你所期望的实际上不是鞍点。
为了更好地逼近连续函数，可以对大部分过采样（根据代码中的 K）函数使用三次插值，例如：
import scipy as sp
import scipy.interpolate
data = [
1.04814804, 0.90445908, 0.62026396, 0.60566623, 0.32295758, 0.26658469, 0.19059289,
0.10281547, 0.08582772, 0.05091265, 0.03391474, 0.03844931, 0.03315003, 0.02838656,
0.03420759, 0.03567401, 0.038203, 0.03530763, 0.04394316, 0.03876966, 0.04156067,
0.03937291, 0.03966426, 0.04438747, 0.03690863, 0.0363976, 0.03171374, 0.03644719,
0.02989291, 0.03166156, 0.0323875, 0.03406287, 0.03691943, 0.02829374, 0.0368121,
0.02971704, 0.03427005, 0.02873735, 0.02843848, 0.02101889, 0.02114978, 0.02128403,
0.0185619, 0.01749904, 0.01441699, 0.02118773, 0.02091855, 0.02431763, 0.02472427,
0.03186318, 0.03205664, 0.03135686, 0.02838413, 0.03206674, 0.02638371, 0.02048122,
0.01502128, 0.0162665, 0.01331485, 0.01569286, 0.00901017, 0.01343558, 0.00908635,
0.00990869, 0.01041151, 0.01063606, 0.00822482, 0.01312368, 0.0115005, 0.00620334,
0.0084177, 0.01058152, 0.01198732, 0.01451455, 0.01605602, 0.01823713, 0.01685975,
0.03161889, 0.0216687, 0.03052391, 0.02220871, 0.02420951, 0.01651778, 0.02066987,
0.01999613, 0.02532265, 0.02589186, 0.02748692, 0.02191687, 0.02612152, 0.02309497,
0.02744753, 0.02619196, 0.02281516, 0.0254296, 0.02732746, 0.02567608, 0.0199178,
0.01831929, 0.01776025]
samples = np.arange(len(data))
f = sp.interpolate.interp1d(samples, data, 'cubic')
K = 10
N = len(data) * K
x = np.linspace(min(samples), max(samples), N)
y = f(x)
然后，所有这些定义都可以通过以下方式进行视觉测试：
import matplotlib.pyplot as plt
plt.figure()
plt.plot(samples, data, label='data')
plt.plot(x, y, label='f')
plt.plot(x, n_derivative(y, 1), label='d1f')
plt.plot(x, n_derivative(y, 2), label='d2f')
plt.plot(x, n_derivative(y, 3), label='d3f')
plt.legend()
for w in np.where(inflections(y))[0]:
plt.axvline(x=x[w])
plt.show()
但即使在这种情况下，该点也不是鞍点。
我有以下 numpy 数组，如上所示。
类似
的函数print(arr.argsort()[:3])
将return三个最低值的三个最低值：
[69 66 70]
我如何 return first index 其中第一个 minimum 或第一个 saddle point （在微积分意义上) 以先到数组为准？
在这种情况下，index 2 和 3 处的两个数字 0.62026396 0.60566623 是第一个 saddle point（这不是真正的 saddle point，因为斜率不t变平了，但它显然打破了那里的硬下坡。也许添加一个"flattens"意味着什么的阈值）。由于函数在第一个 saddle point 之前不会上升，因此第一个 mimimum 出现在 saddle point 之后，这是我感兴趣的索引。
[1.04814804 0.90445908 0.62026396 0.60566623 0.32295758 0.26658469
0.19059289 0.10281547 0.08582772 0.05091265 0.03391474 0.03844931
0.03315003 0.02838656 0.03420759 0.03567401 0.038203 0.03530763
0.04394316 0.03876966 0.04156067 0.03937291 0.03966426 0.04438747
0.03690863 0.0363976 0.03171374 0.03644719 0.02989291 0.03166156
0.0323875 0.03406287 0.03691943 0.02829374 0.0368121 0.02971704
0.03427005 0.02873735 0.02843848 0.02101889 0.02114978 0.02128403
0.0185619 0.01749904 0.01441699 0.02118773 0.02091855 0.02431763
0.02472427 0.03186318 0.03205664 0.03135686 0.02838413 0.03206674
0.02638371 0.02048122 0.01502128 0.0162665 0.01331485 0.01569286
0.00901017 0.01343558 0.00908635 0.00990869 0.01041151 0.01063606
0.00822482 0.01312368 0.0115005 0.00620334 0.0084177 0.01058152
0.01198732 0.01451455 0.01605602 0.01823713 0.01685975 0.03161889
0.0216687 0.03052391 0.02220871 0.02420951 0.01651778 0.02066987
0.01999613 0.02532265 0.02589186 0.02748692 0.02191687 0.02612152
0.02309497 0.02744753 0.02619196 0.02281516 0.0254296 0.02732746
0.02567608 0.0199178 0.01831929 0.01776025]
您可以使用 np.gradient 或 np.diff 来评估差异（第一个计算中心差异，第二个只是 x[1:] - x[:-1]），然后使用 np.sign 获取渐变符号，另一个 np.diff 查看符号变化的位置。然后过滤正号变化（对应最小值）：
np.where(np.diff(np.sign(np.gradient(x))) > 0)[0][0]+2 #add 2 as each time you call np.gradient or np.diff you are substracting 1 in size, the first [0] is to get the positions, the second [0] is to get the "first" element
>> 8
x[np.where(np.diff(np.sign(np.gradient(x))) > 0)[0][0]+2]
>> 0.03420759
环顾四周并根据给出的两个建议（到目前为止），我这样做了：
import scipy
from scipy import interpolate
x = np.arange(0, 100)
spl = scipy.interpolate.splrep(x,arr,k=3) # no smoothing, 3rd order spline
ddy = scipy.interpolate.splev(x,spl,der=2) # use those knots to get second derivative
print(ddy)
asign = np.sign(ddy)
signchange = ((np.roll(asign, 1) - asign) != 0).astype(int)
print(signchange)
这给了我 second derivative，然后我可以对其进行分析，例如，查看符号变化发生的位置：
[-0.894053 -0.14050616 0.61304067 -0.69407217 0.55458251 -0.16624336
-0.0073225 0.12481963 -0.067218 0.03648846 0.02876712 -0.02236204
0.00167794 0.01886512 -0.0136314 0.00953279 -0.01812436 0.03041855
-0.03436446 0.02418512 -0.01458896 0.00429809 0.01227133 -0.02679232
0.02168571 -0.0181437 0.02585209 -0.02876075 0.0214645 -0.00715966
0.0009179 0.00918466 -0.03056938 0.04419937 -0.0433638 0.03557532
-0.02904901 0.02010647 -0.0199739 0.0170648 -0.00298236 -0.00511529
0.00630525 -0.01015011 0.02218007 -0.01945341 0.01339405 -0.01211326
0.01710444 -0.01591092 0.00486652 -0.00891456 0.01715403 -0.01976949
0.00573004 -0.00446743 0.01479495 -0.01448144 0.01794968 -0.02533936
0.02904355 -0.02418628 0.01505374 -0.00499926 0.00302616 -0.00877499
0.01625907 -0.01240068 -0.00578862 0.01351128 -0.00318733 -0.0010652
0.0029 -0.0038062 0.0064102 -0.01799678 0.04422601 -0.0620881
0.05587037 -0.04856099 0.03535114 -0.03094757 0.03028399 -0.01912546
0.01726283 -0.01392421 0.00989012 -0.01948119 0.02504401 -0.02204667
0.0197554 -0.01270022 -0.00260326 0.01038581 -0.00299247 -0.00271539
-0.00744152 0.00784016 0.00103947 -0.00576122]
[0 0 1 1 1 1 0 1 1 1 0 1 1 0 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 0 1 1 1 1 1
1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 0 1 0 1]
这就是我检测局部 maxima/minima、拐点和鞍点的方式。
首先定义以下函数
import numpy as np
def n_derivative(arr, degree=1):
"""Compute the n-th derivative."""
result = arr.copy()
for i in range(degree):
result = np.gradient(result)
return result
def sign_change(arr):
"""Detect sign changes."""
sign = np.sign(arr)
result = ((np.roll(sign, 1) - sign) != 0).astype(bool)
result[0] = False
return result
def zeroes(arr, threshold=1e-8):
"""Find zeroes of an array."""
return sign_change(arr) | (abs(arr) < threshold)
我们现在可以使用 derivative test
关键点的一阶导数将为零。
def critical_points(arr):
return zeroes(n_derivative(arr, 1))
如果临界点的二阶导数不为零，则该点要么是最大值，要么是最小值：
def maxima_minima(arr):
return zeroes(n_derivative(arr, 1)) & ~zeroes(n_derivative(arr, 2))
def maxima(arr):
return zeroes(n_derivative(arr, 1)) & (n_derivative(arr, 2) < 0)
def minima(arr):
return zeroes(n_derivative(arr, 1)) & (n_derivative(arr, 2) > 0)
如果二阶导数为零但三阶导数不为零，则该点为拐点：
def inflections(arr):
return zeroes(n_derivative(arr, 2)) & ~zeroes(n_derivative(arr, 3))
如果临界点的二阶导数为零，但三阶导数不为零，则这是一个鞍点：
def inflections(arr):
return zeroes(n_derivative(arr, 1)) & zeroes(n_derivative(arr, 2)) & ~zeroes(n_derivative(arr, 3))
请注意，此方法在数值上不稳定，因为一方面在某个任意阈值定义上检测到零，另一方面不同的采样可能导致函数/数组不可微分。 因此，根据这个定义，你所期望的实际上不是鞍点。
为了更好地逼近连续函数，可以对大部分过采样（根据代码中的 K）函数使用三次插值，例如：
import scipy as sp
import scipy.interpolate
data = [
1.04814804, 0.90445908, 0.62026396, 0.60566623, 0.32295758, 0.26658469, 0.19059289,
0.10281547, 0.08582772, 0.05091265, 0.03391474, 0.03844931, 0.03315003, 0.02838656,
0.03420759, 0.03567401, 0.038203, 0.03530763, 0.04394316, 0.03876966, 0.04156067,
0.03937291, 0.03966426, 0.04438747, 0.03690863, 0.0363976, 0.03171374, 0.03644719,
0.02989291, 0.03166156, 0.0323875, 0.03406287, 0.03691943, 0.02829374, 0.0368121,
0.02971704, 0.03427005, 0.02873735, 0.02843848, 0.02101889, 0.02114978, 0.02128403,
0.0185619, 0.01749904, 0.01441699, 0.02118773, 0.02091855, 0.02431763, 0.02472427,
0.03186318, 0.03205664, 0.03135686, 0.02838413, 0.03206674, 0.02638371, 0.02048122,
0.01502128, 0.0162665, 0.01331485, 0.01569286, 0.00901017, 0.01343558, 0.00908635,
0.00990869, 0.01041151, 0.01063606, 0.00822482, 0.01312368, 0.0115005, 0.00620334,
0.0084177, 0.01058152, 0.01198732, 0.01451455, 0.01605602, 0.01823713, 0.01685975,
0.03161889, 0.0216687, 0.03052391, 0.02220871, 0.02420951, 0.01651778, 0.02066987,
0.01999613, 0.02532265, 0.02589186, 0.02748692, 0.02191687, 0.02612152, 0.02309497,
0.02744753, 0.02619196, 0.02281516, 0.0254296, 0.02732746, 0.02567608, 0.0199178,
0.01831929, 0.01776025]
samples = np.arange(len(data))
f = sp.interpolate.interp1d(samples, data, 'cubic')
K = 10
N = len(data) * K
x = np.linspace(min(samples), max(samples), N)
y = f(x)
然后，所有这些定义都可以通过以下方式进行视觉测试：
import matplotlib.pyplot as plt
plt.figure()
plt.plot(samples, data, label='data')
plt.plot(x, y, label='f')
plt.plot(x, n_derivative(y, 1), label='d1f')
plt.plot(x, n_derivative(y, 2), label='d2f')
plt.plot(x, n_derivative(y, 3), label='d3f')
plt.legend()
for w in np.where(inflections(y))[0]:
plt.axvline(x=x[w])
plt.show()
但即使在这种情况下，该点也不是鞍点。