就速度而言，numpy.inner可能是固定的最佳选择n。numpy.trapz不过应该收敛得更快。无论哪种方式，如果您担心速度，您还应该考虑到功能本身的评估也需要一些时间。
下面是我使用不同的内部产品实现运行的一些简单基准。
计时
下图显示了仅计算积分的运行时间，即不计算函数评估。虽然numpy.trapz是一个常数因子，速度较慢，但numpy.inner与直接调用 BLAS 一样快。正如 Ophion 指出的那样，在numpy.inner内部调用 BLAS 可能需要一些输入检查开销。
查看评估函数本身所需的时间也很有趣，这当然必须用于计算内积。下面的图显示了标准超越函数numpy.sin和numpy.sqrt的评估numpy.exp。对于产品的评估和求和，缩放当然是相同的，并且所需的总时间是可比的
错误
最后，还应该考虑不同方法的准确性，这才是真正有趣的地方。下面是计算 的不同实现的收敛图
。在这里我们可以看到，它numpy.trapz实际上比其他两种实现的扩展性要好得多，在我用完内存之前，它们甚至没有达到机器精度。
结论
考虑到 的不良收敛特性numpy.inner，我会选择numpy.trapz。但即便如此，也需要大量的集成节点才能获得令人满意的精度。由于您的集成域是固定的，您甚至可以尝试使用更高阶的正交。
代码
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sls
from scipy.linalg.blas import ddot
import timeit
## Define inner product.
def l2_inner_blas( f, g, dx ):
return ddot( f, g )*dx / np.pi
def l2_inner( f, g, dx ):
return np.inner( f, g )*dx / np.pi
def l2_inner_trapz( f, g, dx ):
return np.trapz(f*g,dx=dx) / np.pi
sin1 = lambda x: np.sin( x )
sin2 = lambda x: np.sin( 2.0 * x)
## Timing setups.
setup1 = "import numpy as np; from __main__ import l2_inner,"
setup1 += "l2_inner_trapz, l2_inner_blas, sin1, sin2;"
setup1 += "n=%d; x=np.linspace(-np.pi,np.pi,n); dx=2.0*np.pi/(n-1);"
setup1 += "f=sin1(x); g=sin2(x);"
def time( n ):
setupstr = setup1 % n
time1 = timeit.timeit( 'l2_inner( f, g, dx)', setupstr, number=10 )
time2 = timeit.timeit( 'l2_inner_blas( f, g, dx)', setupstr, number=10 )
time3 = timeit.timeit( 'l2_inner_trapz( f, g, dx)', setupstr, number=10 )
return (time1, time2, time3)
setup2 = "import numpy as np; x = np.linspace(-np.pi,np.pi,%d);"
def time_eval( n ):
setupstr = setup2 % n
time_sin = timeit.timeit( 'np.sin(x)', setupstr, number=10 )
time_sqrt = timeit.timeit( 'np.sqrt(x)', setupstr, number=10 )
time_exp = timeit.timeit( 'np.exp(x)', setupstr, number=10 )
return (time_sin, time_sqrt, time_exp)
## Perform timing for vector product.
times = np.zeros( (7,3) )
for i in range(7):
times[i,:] = time( 10**(i+1) )
x = 10**np.arange(1,8,1)
f, ax = plt.subplots()
ax.set( xscale='log', yscale='log', title='Inner vs. BLAS vs. trapz', \
ylabel='time [s]', xlabel='n')
ax.plot( x, times[:,0], label='numpy.inner' )
ax.plot( x, times[:,1], label='scipy.linalg.blas.ddot')
ax.plot( x, times[:,2], label='numpy.trapz')
plt.legend()
## Perform timing for function evaluation.
times_eval = np.zeros( (7,3) )
for i in range(7):
times_eval[i,:] = time_eval( 10**(i+1) )
x = 10**np.arange(1,8,1)
f, ax = plt.subplots()
ax.set( xscale='log', yscale='log', title='sin vs. sqrt vs. exp', \
ylabel='time [s]', xlabel='n')
ax.plot( x, times_eval[:,0], label='numpy.sin' )
ax.plot( x, times_eval[:,1], label='numpy.sqrt')
ax.plot( x, times_eval[:,2], label='numpy.exp' )
plt.legend()
## Test convergence.
def error( n ):
x = np.linspace( -np.pi, np.pi, n )
dx = 2.0 * np.pi / (n-1)
f = np.exp( x )
l2 = 0.5/np.pi*(np.exp(2.0*np.pi) - np.exp(-2.0*np.pi))
err1 = np.abs( (l2 - l2_inner( f, f, dx )) / l2)
err2 = np.abs( (l2 - l2_inner_blas( f, f, dx )) / l2)
err3 = np.abs( (l2 - l2_inner_trapz( f, f, dx )) / l2)
return (err1, err2, err3)
acc = np.zeros( (7,3) )
for i in range(7):
acc[i,:] = error( 10**(i+1) )
x = 10**np.arange(1,8,1)
f, ax = plt.subplots()
ax.plot( x, acc[:,0], label='numpy.inner' )
ax.plot( x, acc[:,1], label='scipy.linalg.blas.ddot')
ax.plot( x, acc[:,2], label='numpy.trapz')
ax.set( xscale='log', yscale='log', title=r'$\langle \exp(x), \exp(x) \rangle$', \
ylabel='Relative Error', xlabel='n')
plt.legend()
