为什么Fortran中的一元Horner比NumPy的更快，而二元Horner是n

Question

我想用Python执行多项式演算。numpy中的polynomial包对我来说不够快。因此，我决定重写Fortran中的几个函数，并使用f2py创建可轻松导入Python的共享库。目前，我正在对照numpy对应的单变量和二元多项式求值例程。在

在单变量例程中，我使用Horner's method，就像numpy.polynomial.polynomial.polyval一样。我已经观察到Fortran例程比numpy对应程序快的因素随着多项式阶数的增加而增加。在

在双变量程序中，我使用霍纳方法两次。首先在y中，然后在x中。不幸的是，我注意到对于增加多项式阶，numpy对应的函数会赶上并最终超过我的Fortran例程。由于numpy.polynomial.polynomial.polyval2d使用了与我类似的方法，我认为第二个观察结果很奇怪。在

我希望这个结果源于我对Fortran和f2py缺乏经验。有人知道为什么一元程序总是显得优越，而二元程序只对低阶多项式优越？在

编辑 以下是我最新更新的代码、基准脚本和性能图：

多项式

subroutine polyval(p, x, pval, nx)

    implicit none

    real(8), dimension(nx), intent(in) :: p
    real(8), intent(in) :: x
    real(8), intent(out) :: pval
    integer, intent(in) :: nx
    integer :: i

    pval = 0.0d0
    do i = nx, 1, -1
        pval = pval*x + p(i)
    end do

end subroutine polyval

subroutine polyval2(p, x, y, pval, nx, ny)

    implicit none

    real(8), dimension(nx, ny), intent(in) :: p
    real(8), intent(in) :: x, y
    real(8), intent(out) :: pval
    integer, intent(in) :: nx, ny
    real(8) :: tmp
    integer :: i, j

    pval = 0.0d0
    do j = ny, 1, -1
        tmp = 0.0d0
        do i = nx, 1, -1
            tmp = tmp*x + p(i, j)
        end do
        pval = pval*y + tmp
    end do

end subroutine polyval2

subroutine polyval3(p, x, y, z, pval, nx, ny, nz)

    implicit none

    real(8), dimension(nx, ny, nz), intent(in) :: p
    real(8), intent(in) :: x, y, z
    real(8), intent(out) :: pval
    integer, intent(in) :: nx, ny, nz
    real(8) :: tmp, tmp2
    integer :: i, j, k

    pval = 0.0d0
    do k = nz, 1, -1
        tmp2 = 0.0d0
        do j = ny, 1, -1
            tmp = 0.0d0
            do i = nx, 1, -1
                tmp = tmp*x + p(i, j, k)
            end do
            tmp2 = tmp2*y + tmp
        end do
        pval = pval*z + tmp2
    end do

end subroutine polyval3

基准.py（使用此脚本生成绘图）

^{pr2}$

结果

编辑对steabert提案的修正

subroutine polyval(p, x, pval, nx)

    implicit none

    real*8, dimension(nx), intent(in) :: p
    real*8, intent(in) :: x
    real*8, intent(out) :: pval
    integer, intent(in) :: nx

    integer, parameter :: simd = 8
    real*8 :: tmp(simd), xpower(simd), maxpower
    integer :: i, j, k

    xpower(1) = x
    do i = 2, simd
        xpower(i) = xpower(i-1)*x
    end do
    maxpower = xpower(simd)

    tmp = 0.0d0
    do i = nx+1, simd+2, -simd
        do j = 1, simd
            tmp(j) = tmp(j)*maxpower + p(i-j)*xpower(simd-j+1)
        end do
    end do

    k = mod(nx-1, simd)
    if (k == 0) then
        pval = sum(tmp) + p(1)
    else
        pval = sum(tmp) + p(k+1)
        do i = k, 1, -1
            pval = pval*x + p(i)
        end do
    end if

end subroutine polyval

编辑测试代码，以验证上面的代码对x>；1的结果不佳

import polynomial as P
import numpy.polynomial.polynomial as PP

import numpy as np

for n in xrange(2,100):
    poly1n = np.random.rand(n)
    poly1f = np.asfortranarray(poly1n)

    x = 2

    print np.linalg.norm(P.polyval(poly1f, x) - PP.polyval(x, poly1n)), '\n'

Answer 1

按照其他建议，在计时器之前使用p=np.asfortranarray(p)确实可以使性能与我测试时的numpy不相上下。我将二元工作台的范围扩展到n_bi = np.array([2**i for i in xrange(1, 15)])，这样p矩阵将大于我的L3缓存大小。在

为了进一步优化这一点，我认为自动编译器选项不会有多大帮助，因为内部循环有一个依赖关系。只有当您手动展开它时，ifort才会向量化最里面的循环。对于gfortran，需要-O3和{}。对于受主内存带宽限制的矩阵大小，这将比numpy提高1到3倍的性能优势。在

Update：在将此方法应用于单变量代码并使用f2py --opt='-O3 -ffast-math' -c -m polynomial polynomial.f90进行编译之后，我得到以下源代码和结果基准.py公司名称：

subroutine polyval(p, x, pval, nx)

implicit none

real*8, dimension(nx), intent(in) :: p
real*8, intent(in) :: x
real*8, intent(out) :: pval
integer, intent(in) :: nx

integer, parameter :: simd = 8
real*8 :: tmp(simd), vecx(simd), xfactor
integer :: i, j, k

! precompute factors
do i = 1, simd
    vecx(i)=x**(i-1)
end do
xfactor = x**simd

tmp = 0.0d0
do i = 1, nx, simd
    do k = 1, simd
        tmp(k) = tmp(k)*xfactor + p(nx-(i+k-1)+1)*vecx(simd-k+1)
    end do
end do
pval = sum(tmp)


end subroutine polyval

subroutine polyval2(p, x, y, pval, nx, ny)

implicit none

real*8, dimension(nx, ny), intent(in) :: p
real*8, intent(in) :: x, y
real*8, intent(out) :: pval
integer, intent(in) :: nx, ny

integer, parameter :: simd = 8
real*8 :: tmp(simd), vecx(simd), xfactor
integer :: i, j, k

! precompute factors
do i = 1, simd
    vecx(i)=x**(i-1)
end do
xfactor = x**simd

! horner
pval=0.0d0
do i = 1, ny
    tmp = 0.0d0
    do j = 1, nx, simd
        ! inner vectorizable loop
        do k = 1, simd
            tmp(k) = tmp(k)*xfactor + p(nx-(j+k-1)+1,ny-i+1)*vecx(simd-k+1)
        end do
    end do
    pval = pval*y + sum(tmp)
end do

end subroutine polyval2

更新2：如前所述，此代码不正确，至少当大小不能被simd整除时。它只是展示了手动帮助编译器的概念，所以不要像这样使用它。如果大小不是2的幂次方，则一个小的余数循环必须处理悬空的索引。这样做并不难，以下是单变量情况下的正确步骤，应直接将其扩展到二元：

^{pr2}$

另外，对于非常小的尺寸也要小心，因为时间太短，无法获得准确的性能曲线。而且，关于numpy的相对时间可能是欺骗的，因为numpy的绝对时间可能非常糟糕。以下是最大案件的时间安排：

对于nx=220的单变量，numpy的时间为1.21s，而自定义fortran版本的时间为1.69e-3s。对于nxny=220的二元变量，numpy的时间是8e-3s，定制版本的时间是1.68e-3s。当nxny总大小相同时，单变量和双变量的时间是相同的这一事实非常重要，因为它支持这样一个事实，即代码的执行速度接近内存带宽限制。在

更新3：使用新的python脚本来实现更小的大小，simd=4我获得了以下性能：

更新4：对于正确性，结果在双精度精度精度内是相同的，如果您为单变量示例运行以下python代码，可以看到：

import polynomial as P
import numpy.polynomial.polynomial as PP

import numpy as np

for n in xrange(2,100):
    poly1n = np.random.rand(n)
    poly1f = np.asfortranarray(poly1n)

    x = 2

    print "%18.14e" % P.polyval(poly1f, x)
    print "%18.14e" % PP.polyval(x, poly1n)
    print (P.polyval(poly1f, x) - PP.polyval(x, poly1n))/PP.polyval(x,poly1n), '\n'

Answer 2

我猜，你的tmp数组太大了，以至于它需要二级、三级甚至主存访问，而不是缓存。最好把这些循环分开，一次只处理其中的大块（露天开采）。在

Answer 3

在二元情况下，p是一个二维数组。这意味着C与fortran数组的顺序不同。默认情况下，numpy函数提供C排序，显然fortran例程使用fortran排序。在

f2py足够聪明来处理这个问题，并且可以在C和fortran格式的数组之间自动转换。但是，这会导致一些开销，这是性能降低的可能原因之一。您可以通过在计时例程之外使用numpy.asfortranarray手动将p转换为fortran类型来检查原因。当然，要使这一点有意义，在实际用例中，您需要确保输入数组按fortran顺序排列。在

f2py有一个选项-DF2PY_REPORT_ON_ARRAY_COPY，它可以在任何时候警告您数组被复制。在

如果这不是原因，那么您需要考虑更深入的细节，例如您使用的是哪个fortran编译器，以及它应用的优化类型。可能会减慢速度的例子包括在堆上而不是堆栈上分配数组（通过昂贵的调用malloc），尽管我预计这样的影响对于较大的数组来说会变得不那么重要。在

最后，您应该考虑这样一种可能性：对于二元拟合，对于大的N，numpy例程基本上已经处于最佳效率。在这种情况下，numpy例程可能大部分时间都在运行经过优化的C例程，相比之下，python代码的开销变得微不足道。在这种情况下，您不会期望fortran代码显示任何显著的加速。在