python - 在 numpy 中计算 L2 内积？

Question

我在考虑 L2 内积。

我对使用 numpy/scipy 执行这些计算特别感兴趣。我想出的最好办法是执行基于数组的积分，例如 numpy.trapz。

import numpy as np
n=100000.
h=1./n
X = np.linspace(-np.pi,np.pi,n)

def L2_inner_product(f,g):
    return np.trapz(f*g,dx=2*np.pi*h)/np.pi

print L2_inner_product(np.sin(X), np.sin(X))
print L2_inner_product(np.cos(2*X), np.cos(2*X))
print L2_inner_product(np.sin(X), np.cos(X))
print L2_inner_product(np.sin(X), np.cos(3*X))
print L2_inner_product(np.ones(n),np.ones(n))

0.99999
0.99999
-3.86525742539e-18
1.6565388966e-18
1.99998

明确地说，我对使用 Mathematica、Sage 或 Sympy 不感兴趣。我对 numpy/scipy 特别感兴趣，我在其中探索 numpy“数组空间”作为希尔伯特空间的有限子空间。在这些参数中，其他人是否实现了 L2 内积，可能使用 numpy.inner 或 numpy.linalg.norm？

Answer 1

就速度而言，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()