使用NumPy、Numba和Python异步编程的高性能大数据分析与对比

介绍 

几个月前，一位客户问我：“目前大数据分析中较快的Python数据结构对象是什么？ ”我总被问到类似的问题。其中有一些问题很难解决，通常需要多花一些时间才能找到合适的优化解决方案。我一般会在周末和晚上做这些事，并以此为乐。

以前关于这个问题，我第一个简单的答案是Python List对象。我曾在许多数据科学项目中使用List对象，包括数据管道和提取-转换-加载（ETL）生产系统等。然后我就想到了以下问题：我可以使用List对象进行数百万或数十亿行数据的操作和分析吗？如果我将数据科学项目切分成许多小任务，然后使用最新的Python asyncio库异步运行它们，会怎么样呢？基于这些问题，我决定抽出一些时间，借助Python数据生态系统库寻找可用于大数据分析的一些实用解决方案。为了便于读者理解和快速验证结果，程序将计算由浮点数组成的一维NumPy数组的算术平均值、中值和样本标准差。为了对比程序的运行时间，我将使用以下库： 

1. NumPy  - NumPy是用于科学计算的基础Python包。
2. Numba  -  Numba提供了由Python直接编写的高性能函数来加速应用程序的能力。通过几个注释，面向数组和数学计算较多的Python代码就可以被实时编译为原生机器指令。而且Numba拥有类似于C、C++和FORTRAN的性能，无需切换语言或Python解释器。
3. asyncio  -  asyncio是Python异步编程库。

为什么要使用NumPy？ 

正如NumPy网站所说：“NumPy是用于科学计算的基础Python包。它提供了强大的N维数组对象和复杂的（广播）功能。”导入NumPy库之后，Python程序的性能更好、执行速度更快、更容易保证一致性并能方便地使用大量的数学运算和矩阵功能。也许正因为如此，我们不再需要使用Python List对象了？重要的是，许多Python数据生态系统库都基于NumPy之上，像Pandas、SciPy、Matplotlib等等。

用到的Python算法 

在此我将通过算术平均值、中值和样本标准差的简单计算，来展示不同Python程序的运行时间并进行对比。测试数据来自由64位浮点数构成的一维NumPy数组。我将实现以下三种Python算法并对它们进行分析：

1. NumPy数组
2. 结合asyncio异步库的NumPy数组
3. 结合Numba库的NumPy数组

NumPy数组程序 

我们来看看每个算法的代码。每个算法都有遵循面向对象编程（OOP）方法的类对象和主调用程序。类对象包含以下五种方法：

• calculate_number_observation() - 计算观测数
• calculate_arithmetic_mean() - 计算算术平均值
• calculate_median() - 计算中值
• calculate_sample_standard_deviation() - 计算样本标准差
• print_exception_message() - 如果出现异常，打印异常信息

列表1显示了单独使用NumPy数组的汇总统计类对象代码。

import sys
import traceback
import time
from math import sqrt

class SummaryStatistics(object):
    """
    使用标准过程计算观测数、算术平均值、中值和样本标准差 
    """
    def __init__(self):
        pass        

    def calculate_number_observation(self, one_dimensional_array):        
        """
        计算观测数 
        :参数 one_dimensional_array: numpy一维数组 
        :返回值 观测数 
        """

        number_observation = 0
        try:
            number_observation = one_dimensional_array.size   
        except Exception:
            self.print_exception_message()
        return number_observation 

    def calculate_arithmetic_mean(self, one_dimensional_array, number_observation):    

        """
        计算算术平均值 
        :参数 one_dimensional_array: numpy一维数组 
        :参数 number_observation: 观测数 
        :返回值 算术平均值 
        """
        
        arithmetic_mean = 0.0
        try:
            sum_result = 0.0
            for i in range(number_observation):       
                sum_result += one_dimensional_array[i]    
                arithmetic_mean = sum_result / number_observation
        except Exception:
            self.print_exception_message()
        return arithmetic_mean

    def calculate_median(self, one_dimensional_array, number_observation):      
        """
        计算中值 

        :参数 one_dimensional_array: numpy一维数组 
        :参数 number_observation: 观测数 
        :返回值 中值 

        """
        median = 0.0
        try:
            one_dimensional_array.sort()    
            half_position = number_observation // 2
            if not number_observation % 2:
                median = (one_dimensional_array[half_position - 1] + one_dimensional_array[half_position]) / 2.0
            else:
                median = one_dimensional_array[half_position]        
        except Exception:
            self.print_exception_message()
        return median

    def calculate_sample_standard_deviation(self, one_dimensional_array, number_observation, arithmetic_mean):    
        """
        计算样本标准差 
        :参数 one_dimensional_array: numpy一维数组 
        :参数 number_observation: 观测数 
        :参数 arithmetic_mean: 算术平均值 
        :返回值 样本标准差值 
        """
        sample_standard_deviation = 0.0
        try:
            sum_result = 0.0
            for i in range(number_observation):                   
                sum_result += pow((one_dimensional_array[i] - arithmetic_mean), 2)            
            sample_variance = sum_result / (number_observation - 1)            
            sample_standard_deviation = sqrt(sample_variance)        
        except Exception:
            self.print_exception_message()
        return sample_standard_deviation 

    def print_exception_message(self, message_orientation = "horizontal"):
        """
        打印完整的异常信息
        :参数 message_orientation: 水平或垂直 
        :返回值 空
        """
        try:
            exc_type, exc_value, exc_tb = sys.exc_info()            
            file_name, line_number, procedure_name, line_code = traceback.extract_tb(exc_tb)[-1]       
            time_stamp = " [Time Stamp]: " + str(time.strftime("%Y-%m-%d %I:%M:%S %p")) 
            file_name = " [File Name]: " + str(file_name)
            procedure_name = " [Procedure Name]: " + str(procedure_name)
            error_message = " [Error Message]: " + str(exc_value)        
            error_type = " [Error Type]: " + str(exc_type)                    
            line_number = " [Line Number]: " + str(line_number)                
            line_code = " [Line Code]: " + str(line_code) 

            if (message_orientation == "horizontal"):
                print( "An error occurred:{};{};{};{};{};{};{}".format(time_stamp, file_name, procedure_name, error_message, error_type, line_number, line_code))
            elif (message_orientation == "vertical"):
                print( "An error occurred:\n{}\n{}\n{}\n{}\n{}\n{}\n{}".format(time_stamp, file_name, procedure_name, error_message, error_type, line_number, line_code))
            else:
                pass                    
        except Exception:
            pass

列表1.  单独使用NumPy数组的汇总统计类对象代码

列表2显示了汇总统计的主程序。可以看到，程序创建了summary_statistics类对象，然后调用其中的方法。该程序导入NumPy库以生成一维数组，并使用time模块中的clock()方法来计算程序的运行时间。

import time
import numpy as np
 
from class_summary_statistics import SummaryStatistics
 
def main(one_dimensional_array):
    
#     创建汇总统计类对象 
    summary_statistics = SummaryStatistics()
    
#     计算观测数 
    number_observation = summary_statistics.calculate_number_observation(one_dimensional_array)
    print("Number of Observation: {} ".format(number_observation))
    
#     计算算术平均值 
    arithmetic_mean = summary_statistics.calculate_arithmetic_mean(one_dimensional_array, number_observation)
    print("Arithmetic Mean: {} ".format(arithmetic_mean))
    
#     计算中值 
    median = summary_statistics.calculate_median(one_dimensional_array, number_observation)
    print("Median: {} ".format(median))
    
#     计算样本标准差 
    sample_standard_deviation = summary_statistics.calculate_sample_standard_deviation(one_dimensional_array, number_observation, arithmetic_mean)
    print("Sample Standard Deviation: {} ".format(sample_standard_deviation))
 
if __name__ == '__main__':
    start_time = time.clock()  
    one_dimensional_array = np.arange(100000000, dtype=np.float64)        
    main(one_dimensional_array)
    end_time = time.clock()
    print("Program Runtime: {} seconds".format(round(end_time - start_time, 1)))

列表2.  使用NumPy数组的汇总统计类对象代码

当行数达到1百万时，汇总统计主程序将得到以下结果：

观测数：       1000000 
算术平均值：    499999.5 
中值：         499999.5 
样本标准差：    288675.27893349814 
程序运行时间：  1.3秒 

结合asyncio异步库的NumPy数组

列表3显示了结合Python asyncio异步库的汇总统计asyncio类对象代码。请注意，main()方法使用calculate_number_observation()作为第一个也是唯一任务来启动事件循环异步进程。

import sys
import time
import traceback
import asyncio
from math import sqrt

class SummaryStatisticsAsyncio(object):
    """
    使用asyncion库计算观测数、算术平均值、中值和样本标准差 
    """
    def __init__(self):
        pass
    
    async def calculate_number_observation(self, one_dimensional_array):    
        """
        计算观测数 
        :参数 one_dimensional_array: numpy一维数组 
        :返回值 空
        """
        try:
            print('start calculate_number_observation() procedure')   
            await asyncio.sleep(0)
            number_observation = one_dimensional_array.size
            print("Number of Observation: {} ".format(number_observation))    
            print("观测数：  {} ".format(number_observation))    
            await self.calcuate_arithmetic_mean(one_dimensional_array, number_observation)
            print("finished calculate_number_observation() procedure")   
        except Exception:
            self.print_exception_message()
            
    async def calcuate_arithmetic_mean(self, one_dimensional_array, number_observation):    
        """
        计算算术平均值 
        :参数 one_dimensional_array: numpy一维数组 
        :参数 number_observation: 观测数 
        :返回值 空
        """
        try:
            print('start calcuate_arithmetic_mean() procedure')   
            await self.calculate_median(one_dimensional_array, number_observation)
            sum_result = 0.0
            await asyncio.sleep(0)
            for i in range(number_observation):       
                sum_result += one_dimensional_array[i]    
            arithmetic_mean = sum_result / number_observation
            print("Arithmetic Mean: {} ".format(arithmetic_mean))    
            await self.calculate_sample_standard_deviation(one_dimensional_array, number_observation, arithmetic_mean)
            print("finished calcuate_arithmetic_mean() procedure")   
        except Exception:
            self.print_exception_message()
            
    async def calculate_median(self, one_dimensional_array, number_observation):      
        """
        计算中值 
        :参数 one_dimensional_array: numpy一维数组 
        :参数 number_observation: 观测数 
        :返回值 空
        """
        try:
            print('starting calculate_median()')   
            await asyncio.sleep(0)
            one_dimensional_array.sort()    
            half_position = number_observation // 2
            if not number_observation % 2:
                median = (one_dimensional_array[half_position - 1] + one_dimensional_array[half_position]) / 2.0
            else:
                median = one_dimensional_array[half_position]        
            print("Median: {} ".format(median))
            print("finished calculate_median() procedure")   
        except Exception:
            self.print_exception_message()

    async def calculate_sample_standard_deviation(self, one_dimensional_array, number_observation, arithmetic_mean):    
        """
        计算样本标准差 
        :参数 one_dimensional_array: numpy一维数组 
        :参数 number_observation:  观测数 
        :参数 arithmetic_mean: 算术平均值 
        :返回值 空
        """
        try:
            print('start calculate_sample_standard_deviation() procedure')   
            await asyncio.sleep(0)
            sum_result = 0.0
            for i in range(number_observation):                   
                sum_result += pow((one_dimensional_array[i] - arithmetic_mean), 2)           
            sample_variance = sum_result / (number_observation - 1)            
            sample_standard_deviation = sqrt(sample_variance)        
            print("Sample Standard Deviation: {} ".format(sample_standard_deviation))
            print("finished calculate_sample_standard_deviation() procedure")   
        except Exception:
            self.print_exception_message()
            
    def print_exception_message(self, message_orientation = "horizontal"):
        """
        打印完整异常消息 
        :参数 message_orientation: 水平或垂直 
        :返回值 空
        """
        try:
            exc_type, exc_value, exc_tb = sys.exc_info()            
            file_name, line_number, procedure_name, line_code = traceback.extract_tb(exc_tb)[-1]       
            time_stamp = " [Time Stamp]: " + str(time.strftime("%Y-%m-%d %I:%M:%S %p")) 
            file_name = " [File Name]: " + str(file_name)
            procedure_name = " [Procedure Name]: " + str(procedure_name)
            error_message = " [Error Message]: " + str(exc_value)        
            error_type = " [Error Type]: " + str(exc_type)                    
            line_number = " [Line Number]: " + str(line_number)                
            line_code = " [Line Code]: " + str(line_code) 
            if (message_orientation == "horizontal"):
                print( "An error occurred:{};{};{};{};{};{};{}".format(time_stamp, file_name, procedure_name, error_message, error_type, line_number, line_code))
            elif (message_orientation == "vertical"):
                print( "An error occurred:\n{}\n{}\n{}\n{}\n{}\n{}\n{}".format(time_stamp, file_name, procedure_name, error_message, error_type, line_number, line_code))
            else:
                pass                    
        except Exception:
            pass
    
    def main(self, one_dimensional_array):    
        """
        启动事件循环异步进程 
        :参数 one_dimensional_array: numpy一维数组 
        """
        try:
            ioloop = asyncio.get_event_loop()
            tasks = [ioloop.create_task(self.calculate_number_observation(one_dimensional_array))]
            wait_tasks = asyncio.wait(tasks)
            ioloop.run_until_complete(wait_tasks)
            ioloop.close()
        except Exception:

列表3.  结合Python异步库的汇总统计asyncio类对象代码

汇总统计asyncio主程序如列表4所示。可以看到main() 方法是唯一被调用的方法。

import time
import numpy as np
 
from class_summary_statistics_asyncio import SummaryStatisticsAsyncio
 
def main(one_dimensional_array):
    
#     新建汇总统计asyncio类对象
    summary_statistics_asyncio = SummaryStatisticsAsyncio()
 
#     调用main方法
    summary_statistics_asyncio.main(one_dimensional_array)
 
if __name__ == '__main__':
    start_time = time.clock()  
    one_dimensional_array = np.arange(1000000000, dtype=np.float64)        
    main(one_dimensional_array)
    end_time = time.clock()
    print("Program Runtime: {} seconds".format(round(end_time - start_time, 1)))

列表4.  结合Python异步库的汇总统计asyncio主程序代码

当行数达到1百万时，结合asyncio库的汇总统计主程序将得到以下结果。我加入了 开始/结束过程的打印，以展示异步过程在这种特殊情况下的工作原理。

start calculate_number_observation() procedure
观测数：  1000000000
start calcuate_arithmetic_mean() procedure
starting calculate_median()
中值：  499999.5 
finished calculate_median() procedure
算术平均值：  499999.5 
start calculate_sample_standard_deviation() procedure
样本标准差：  288675.27893349814 
finished calculate_sample_standard_deviation() procedure
finished calcuate_arithmetic_mean() procedure
finished calculate_number_observation() procedure
程序运行时间：  1504.4秒 

结合Numba库的NumPy数组

结合Numba库的汇总统计类对象代码如列表5所示。你可以访问Numba在GitHub上的目录，以了解更多关于这个Python的开源NumPy感知优化编译器的信息。值得一提的是Numba支持CUDA GPU编程。下面的代码中，调试代码已被删除，以便在编译模式下运行该程序。

import time
from numba import jit
import numpy as np
from math import sqrt
 
class SummaryStatisticsNumba(object):
    """
    结合numba库计算观测数、算术平均值、中值和样本标准差 
    """
    def __init__(self):
        pass
        
    @jit
    def calculate_number_observation(self, one_dimensional_array):    
        """
        计算观测数 
        :参数 one_dimensional_array: numpy一维数组 
        :返回值 观测数 
        """        
        number_observation = one_dimensional_array.size
        return number_observation
    
    @jit
    def calcuate_arithmetic_mean(self, one_dimensional_array, number_observation):    
        """
        计算算术平均值 
        :参数 one_dimensional_array: numpy一维数组 
        :参数 number_observation: 观测数 
        :返回值 算术平均值 
        """
        sum_result = 0.0
        for i in range(number_observation):       
            sum_result += one_dimensional_array[i]    
        arithmetic_mean = sum_result / number_observation
        return arithmetic_mean
    
    @jit
    def calculate_median(self, one_dimensional_array, number_observation):      
        """
        计算中值 
        :参数 one_dimensional_array: 指numpy一维数组 
        :参数 number_observation: 观测数 
        :返回值 中值 
        """
        one_dimensional_array.sort()    
        half_position = number_observation // 2
        if not number_observation % 2:
            median = (one_dimensional_array[half_position - 1] + one_dimensional_array[half_position]) / 2.0
        else:
            median = one_dimensional_array[half_position]        
        return median
    
    @jit
    def calculate_sample_standard_deviation(self, one_dimensional_array, number_observation, arithmetic_mean):    
        """
        计算样本标准差 
        :参数 one_dimensional_array: numpy一维数组 
        :参数 number_observation: 观测数 
        :参数 arithmetic_mean: 算术平均值 
        :返回值 样本标准差值 
        """
        sum_result = 0.0
        for i in range(number_observation):                   
            sum_result += pow((one_dimensional_array[i] - arithmetic_mean), 2)            
        sample_variance = sum_result / (number_observation - 1)            
        sample_standard_deviation = sqrt(sample_variance)        
        return sample_standard_deviation

列表5  结合Numba库的汇总统计类对象代码

汇总统计Numba主程序如列表6所示。 

import time
import numpy as np

from class_summary_statistics_numba import SummaryStatisticsNumba

def main(one_dimensional_array):

#    创建类汇总统计numba类对象 
    class_summary_statistics_numba = SummaryStatisticsNumba()    

#     计算观测数 
    number_observation = class_summary_statistics_numba.calculate_number_observation(one_dimensional_array)
    print("Number of Observation: {} ".format(number_observation))
 
#     计算算术平均值 
    arithmetic_mean = class_summary_statistics_numba.calcuate_arithmetic_mean(one_dimensional_array, number_observation)
    print("Arithmetic Mean: {} ".format(arithmetic_mean))

#     计算中值 
    median = class_summary_statistics_numba.calculate_median(one_dimensional_array, number_observation)
    print("Median: {} ".format(median))

#     计算样本标准差 
    sample_standard_deviation = class_summary_statistics_numba.calculate_sample_standard_deviation(one_dimensional_array, number_observation, arithmetic_mean)
    print("Sample Standard Deviation: {} ".format(sample_standard_deviation))
   
if __name__ == '__main__':
    start_time = time.clock()  
    one_dimensional_array = np.arange(1000000000, dtype=np.float64)        
    main(one_dimensional_array)
    end_time = time.clock()
    print("Program Runtime: {} seconds".format(round(end_time - start_time, 1))) 

列表6.  汇总统计Numba主程序

当行数达到十亿时，汇总统计Numba主程序将得到以下结果。

观测数：  1000000000 
算术平均值：  499999999.067109 
中值：  499999999.5 
样本标准差：  288675134.73899055 
程序运行时间：  40.2秒 

NumPy数组达到十亿行时，计算在40.2秒内就完成了，这真是一个激动人心的结果。我认为现在是时候在大数据项目中更多地使用Numba库和NumPy数组了。当然某些特殊场景可能还需要进一步研究和测试。

笔记本硬件参数

下面是我运行上面这些Python程序所使用的笔记本电脑硬件参数：

• Windows 10 64位操作系统
• 英特尔酷睿™i7-2670QM CPU @2.20 GHz
• 16 GB 内存

程序运行时间对比 

表1显示了数据行数为100万、1000万、1亿和10亿时不同程序的运行时间。

表1：  程序运行时间对比

结论 

1. 单独使用NumPy数组与结合asyncio异步库的NumPy数组之间没有明显差别。由于本次计算量不足以证明Python数据科学项目中asyncio异步库的性能，因此可能需要进行更多的研究来找到它适合的应用场景。
2. 与单独使用NumPy数组或结合asyncio异步库的NumPy数组相比，将NumPy数组与Numba库组合具有最佳的数据操作和分析性能。当Numpy数组中有十亿行数据时，执行时间竟然只需要40.2秒，令我印象深刻。我怀疑目前的R程序是否也可以达到这样的速度。如果不能，也许现在就是R程序的程序员学习Python及其数据生态系统库的时候了。除此之外，请务必采用持续集成软件开发和部署实践，使用面向对象编程方法论来为实际的生产环境编写Python程序。
3. 在具有不同类型数据源的数据管道和提取-转换-加载（ETL）系统项目中，使用结合Numba库的NumPy数组是目前大数据分析的最佳编程实践之一。你不再需要使用Python List对象。

最后，献上我在软件业务应用设计和开发领域摸爬滚打25年之后最喜欢的句子之一：

“一个差的程序员写代码只是为了能运行，而一个好的程序员写代码不只是为了运行程序，也为了以后可以更方便地维护和更新” -  Ernest Bonat博士