为什么这个(非常简单的)矢量化代码比Numpy慢几个数量级？

Question

我非常清楚地知道朱莉娅的理念有关矢量化之类的,我可以接受有我的代码运行速度比numpy的甚至十几倍速度较慢,但为什么这个代码运行很多很多比NumPy的慢？也许代码中有错误; 我无法想象这个问题与使用矢量化而不是循环有关.

我正在使用矢量化,因为我的要求不强; 此外,内存分配似乎并不那么令人难以置信(并且Numpy证明它非常快).代码也很容易阅读.

下面的代码是用Python编写的; 它计算复平面的一部分上的广义连续分数.连续分数由两个不同的函数定义; 我使用此代码在此博客上绘制图片:https://kettenreihen.wordpress.com/

Python代码如下; 它在大约2秒内计算640x640的值:

def K(a, b, C):
    s = C.shape
    P1 = np.ones(s, dtype=np.complex)
    P2 = np.zeros(s, dtype=np.complex)
    Q1 = np.zeros(s, dtype=np.complex)
    Q2 = np.ones(s, dtype=np.complex)
    for n in range(1, 65):
        A, B = a(C, n), b(C, n)
        P = A*P2 + B*P1
        Q = A*Q2 + B*Q1
        P1, P2 = P2, P
        Q1, Q2 = Q2, Q
    return P2/Q2

以下Julia代码也应该这样做,但计算相同的东西需要2到3分钟.

function K(a, b, C)
    s = size(C)
    P1 = ones(Complex, s)
    P2 = zeros(Complex, s)
    Q1 = zeros(Complex, s)
    Q2 = ones(Complex, s)
    for n = 1:64
        println(n)
        A, B = a(C, n), b(C, n)
        P = A.*P2 + B.*P1
        Q = A.*Q2 + B.*Q1
        P1, P2 = P2, P
        Q1, Q2 = Q2, Q
    end
    return P2./Q2
end

Answer 1

您正在使用抽象元素类型分配矩阵:Complex类型是所有特定Complex{T}类型的抽象超类型.你想要的是一些具体元素类型的数组,如Complex128 == Complex{Float64}或Complex64 == Complex{Float32}.据推测,在NumPy dtype=np.complex中指的是一种特定的复杂类型,可能相当于Complex128.

如果你想编写通用的代码并且可以用于不同类型的C矩阵,那么假设C是一个复杂的矩阵,你想要的是创建相同元素类型和形状的1和0的矩阵,你可以调用ones和zeros函数on C以获得具有正确元素类型和形状的矩阵:

function K(a, b, C)
    P1 = ones(C)
    P2 = zeros(C)
    Q1 = zeros(C)
    Q2 = ones(C)
    for n = 1:64
        println(n)
        A, B = a(C, n), b(C, n)
        P = A.*P2 + B.*P1
        Q = A.*Q2 + B.*Q1
        P1, P2 = P2, P
        Q1, Q2 = Q2, Q
    end
    return P2./Q2
end

希望这有助于提高性能.通过预先分配三个矩阵并在矩阵中旋转来执行操作,您可以获得更高的性能.然而,方便的语法支持还没有使它成为一个稳定版本,所以在Julia 0.5上它仍然有点冗长,但可以让你比矢量化版本提升性能.