科学计算语言 Julia 在遗传算法中的性能优势实测

本文对比了Python与Julia在科学计算中的性能表现，以遗传算法解决背包问题为例，测试了两种语言在相同功能下的运行效率。通过在Orange Pi Zero 3上运行10,000次循环，结果显示Julia的平均单次运行时间仅为Python的33.8%（0.0842秒 vs 0.2492秒），总时间节省约66%。Julia凭借其高效的矩阵运算和即时编译优势，显著降低了计算成本，尤其适合需要高频率迭代的科学计算场景。文章提供了完整的代码实现与性能分析图表，为物理或计算科学领域的研究者提供了语言选型的参考依据。

引言

作为一名基本上只关心CS for Physics 或 AI for Physics的物理牲，高中的时候对在零基础的情况下对计算机由莫名的兴趣， 现在想起来依然不会觉得只是高中太压抑了云云，折腾了一整天弄一个高门槛的Typecho博客就说明我对CS/AI还有孩童的本能仰慕，亦或是机械崇拜。总之，对于一些基础的计算机知识，或者是新兴的产品，我还是有一些兴趣去了解的。

Julia

一群拥有各种语言丰富编程经验的Matlab高级用户，对现有的科学计算编程工具感到不满——这些软件对自己专长的领域表现得非常棒，但在其它领域却非常糟糕。他们想要的是一个开源的软件，它要像C语言一般快速而又拥有如同Ruby的动态性；要具有Lisp般真正的同像性而又有Matlab般熟悉的数学记号；要像Python般通用、像R般在统计分析上得心应手、像Perl般自然地处理字符串、像Matlab般具有强大的线性代数运算能力、像shell般胶水语言的能力，易于学习而又不让真正的黑客感到无聊；还有，它应该是交互式的，同时又是编译型的。

——WiKi for Julia

早在高中时期我就了解了Julia这门语言，甚至早于Python。对于物理专业，尤其是实验方向的学生，对于代码速度基本上没有太多的要求，MATLAB的可视化页面对我这种懒得Coding的人很友好。2024年我选修了《进化优化算法》课程，遗憾的是一些学习笔记本来存储在我的服务器上，但是数据丢失了。第一份大作业是一份遗传算法的题，解背包问题。我没有靠AI硬写了一份出来，在后来我逐渐理解，遗传算法终归是带有随机性的算法，且需要大量的计算成本，优化计算成本并提高准确率是算法优化的基本方向，同时也是我们作业评分的最重要标准。当时在《电动力学》课上测试代码，每次都要花比较长的时间来检验准确率，例如我刚刚测试了一下循环10000次需要311s。这个问题在后来的蚁群优化算法中更加严重，为了跑两张收敛速度对比图，会花上很长时间。

很大的问题是我选择了Python语言，出了名的慢，不过当时没有功夫折腾，我很需要一个高速的代码，又要好写，又要好用。很久之后我突然想起了Julia。本段开头的引言里说的算是比较明确了，这是一个很诱人的项目，实际上目前已经有高校用作科学计算。但是小众有小众的理由，有用户声称机器学习给出了错误的结果。具体的原因与解决方案我并没有继续了解，毕竟我现在还是一个Julia小白^^。

总之，Julia值得一试。我翻出来了我的Orange Pi Zero3，尽管我的数据库丢失就是它干的，但是我依然决定信任。我打算在上面做一些基础的测试。

不同代码运行速度的直观对比

折腾了好半天，终于在Pi上搭建好了最基础的Python和Julia，我将我的遗传算法代码传到Pi上，并尝试跑10000次。利用DS写了一个相同功能的jl代码，我试着分别输出它们的运行速度比较图。

测试结果

Python

测试的输出与结果如下。图一给出了每次循环所用的单次运行时间，图二给出的则是累计运行时间。实际上感觉图一没什么大用，由图二可以看出累计运行时间与单次运行时间的线性关系。

测试次数：10000, 错误次数: 0, 错误率: 0.0%
总运行时间：2492.22 秒
平均每次运行时间：0.2492 秒

Julia

测试的输出与结果如下。图一给出的数据出现了间隔性的激增，单次运行时间会在一定次数后突然增大，可能与垃圾回收（Garbage Collection, GC）有关:内存使用量达到阈值后触发GC，导致单次运行时间突然增加，之后内存被释放，时间恢复正常。图二给出了相同的线性关系，总体上程序高速运行。
测试次数：10000, 错误次数: 0, 错误率: 0.0%
总运行时间：842.06 秒
平均每次运行时间：0.0842 秒

测试时间对比

将单次运行时间对比绘出，可以发现Julia的运行速度几乎是Python的3倍。Julia具有强大的线性代数运算能力，在矩阵运算时对Python几乎形成降维打击，这也是我在《进化优化算法》课程上一直想要解决的问题：Python算矩阵实在太慢了。

结语

本文基于已有的遗传算法Python代码，将其重构为相同功能的Julia代码并在同一台机器上运行10000次循环，发现Julia的运行时间是Python的33.8%。这个结果发挥了Julia作为科学计算语言的长处。如果你也有关于减少运算成本的需求，想要一类简便开源的科学计算语言，这篇文章可以给你带来直观的感受，希望对你有帮助。

附录

测试中使用的代码。

import numpy as np
from collections import Counter
import time
import matplotlib.pyplot as plt
import os

# 创建监控目录
if not os.path.exists('monitor'):
    os.makedirs('monitor')

# 定义原始数据
c_j = np.array([91, 72, 90, 46, 55, 8, 35, 75, 61, 15, 77, 40, 63, 75, 29, 75, 17, 78, 40, 44])
w_j = np.array([84, 83, 43, 4,  44, 6, 82, 92, 25, 83, 56, 18, 58, 14, 48, 70, 96, 32, 68, 92])
w_limit = 879
p_c = 0.75  # 交叉率
p_m = 0.008  # 变异率
T = 400  # 迭代次数
N = 100  # 种群规模
N_2 = len(c_j)  # 维度

def initialize_population(N, N_2):
    """采用随机二进制编码初始化种群"""
    return np.random.randint(2, size=(N, N_2))

def calculate_sums(x):
    """计算个体总价值c_sum和个体重量w_sum"""
    c_sum = np.dot(c_j, x.T)
    w_sum = np.dot(w_j, x.T)
    return c_sum, w_sum

def evaluate(c_sum, w_sum):
    """计算适应度函数和轮盘赌归一化概率"""
    p = 1 - np.floor(w_sum / w_limit)  # 罚函数法：超重个体为p=0，反之为1
    eval = p * c_sum  # 适应度函数
    P = eval / np.sum(eval)  # 归一化概率
    cum_P = np.cumsum(P)  # 累计概率，用于轮盘赌
    return eval, cum_P

def select_parents(cum_P):
    """选择交叉操作2个不同的父代"""
    select_index = np.searchsorted(cum_P, np.random.rand())
    while True:
        select_index_2 = np.searchsorted(cum_P, np.random.rand())
        if select_index_2 != select_index:
            break
    return select_index, select_index_2

def crossover(parent1, parent2):
    """单点交叉操作"""
    r = np.random.randint(1, N_2 + 1)  # 确定交叉点
    offspring1 = np.concatenate((parent1[:r], parent2[r:]), axis=0)
    offspring2 = np.concatenate((parent2[:r], parent1[r:]), axis=0)
    return offspring1, offspring2

def mutate(x, mutation_rate):
    """变异操作"""
    mutation_mask = np.random.rand(*x.shape) < mutation_rate
    x[mutation_mask] = 1 - x[mutation_mask]

def perform_evolution(x):
    """执行进化过程"""
    x2 = x.copy()  # 复制种群作为子代

    c_sum, w_sum = calculate_sums(x2)
    eval, cum_P = evaluate(c_sum, w_sum)

    select_index, select_index_2 = select_parents(cum_P)
    parent1 = x2[select_index]
    parent2 = x2[select_index_2]

    # 交叉操作
    if np.random.rand() > p_c:
        offspring1, offspring2 = crossover(parent1, parent2)
        x2[select_index] = offspring1
        x2[select_index_2] = offspring2

    # 变异操作
    mutate(x2, p_m)

    # 计算新的c_sum,w_sum,eval
    c_sum2, w_sum2 = calculate_sums(x2)
    eval2, _ = evaluate(c_sum2, w_sum2)

    # 合并并排序选择
    eval0 = np.concatenate((eval, eval2), axis=0)
    x0 = np.concatenate((x, x2), axis=0)
    best_eval_indices = np.argsort(eval0)[-N:]  # 根据适应度排序选择前N个

    return x0[best_eval_indices]

def run_simulation(runs, target_value1, target_value2):
    """重复运行程序，测试程序稳定性并记录时间"""
    not_target_count = 0
    x = initialize_population(N, N_2)
    
    # 存储时间数据
    single_run_times = []
    cumulative_times = []
    total_start_time = time.time()
    
    for run_num in range(1, runs + 1):
        run_start_time = time.time()
        
        for _ in range(T):
            x = perform_evolution(x)
            
        # 找到对应于最大价值，容量及其对应的个体
        w_sums = np.dot(w_j, x.T)
        c_sums = np.dot(c_j, x.T)
        
        counts_w = Counter(w_sums)
        most_common_wsum, most_common_count_w = counts_w.most_common(1)[0]

        counts_c = Counter(c_sums)
        most_common_csum, most_common_count_c = counts_c.most_common(1)[0]

        best_individual = x[np.argmax(c_sums)]
        
        # 检查是否不等于目标值
        if most_common_wsum != target_value1 and most_common_csum != target_value2:
            not_target_count += 1
            
        # 记录时间
        run_time = time.time() - run_start_time
        single_run_times.append(run_time)
        cumulative_times.append(time.time() - total_start_time)
        
        # 每100次打印进度
        if run_num % 100 == 0:
            print(f"已完成 {run_num}/{runs} 次运行")
    
    # 绘制图表
    plt.figure(figsize=(12, 6))
    
    # 图表1: 单次运行时间
    plt.subplot(1, 2, 1)
    plt.plot(range(1, runs + 1), single_run_times)
    plt.xlabel('运行次数')
    plt.ylabel('单次运行时间 (秒)')
    plt.title('单次运行时间变化')
    
    # 图表2: 累计运行时间
    plt.subplot(1, 2, 2)
    plt.plot(range(1, runs + 1), cumulative_times)
    plt.xlabel('运行次数')
    plt.ylabel('累计运行时间 (秒)')
    plt.title('累计运行时间变化')
    
    plt.tight_layout()
    plt.savefig('monitor/python_performance.png')
    plt.close()
    
    return not_target_count

if __name__ == "__main__":
    target_value1 = 879  
    target_value2 = 1025
    runs = 10000
    
    print("开始运行遗传算法...")
    start_time = time.time()
    
    result = run_simulation(runs, target_value1, target_value2)
    rate = 100 * result / runs
    
    elapsed = time.time() - start_time
    print(f"测试次数：{runs}, 错误次数: {result}, 错误率: {rate}%")
    print(f"总运行时间：{elapsed:.2f} 秒")
    print(f"平均每次运行时间：{elapsed/runs:.4f} 秒")

using Random
using Plots
using Statistics
using Dates

# 创建监控目录
if !isdir("monitor")
    mkdir("monitor")
end

# 定义原始数据
const c_j = [91, 72, 90, 46, 55, 8, 35, 75, 61, 15, 77, 40, 63, 75, 29, 75, 17, 78, 40, 44]
const w_j = [84, 83, 43, 4,  44, 6, 82, 92, 25, 83, 56, 18, 58, 14, 48, 70, 96, 32, 68, 92]
const w_limit = 879
const p_c = 0.75  # 交叉率
const p_m = 0.008  # 变异率
const T = 400  # 迭代次数
const N = 100  # 种群规模
const N_2 = length(c_j)  # 维度

function initialize_population(pop_size::Int, dim::Int)
    """采用随机二进制编码初始化种群"""
    return rand(0:1, pop_size, dim)
end

function calculate_sums(x::Matrix{Int})
    """计算个体总价值c_sum和个体重量w_sum"""
    c_sum = x * c_j
    w_sum = x * w_j
    return c_sum, w_sum
end

function evaluate(c_sum::Vector{Int}, w_sum::Vector{Int})
    """计算适应度函数和轮盘赌归一化概率"""
    p = 1 .- floor.(Int, w_sum ./ w_limit)  # 罚函数法：超重个体为p=0，反之为1
    eval = p .* c_sum  # 适应度函数
    P = eval ./ sum(eval)  # 归一化概率
    cum_P = cumsum(P)  # 累计概率，用于轮盘赌
    return eval, cum_P
end

function select_parents(cum_P::Vector{Float64})
    """选择交叉操作2个不同的父代"""
    select_index = searchsortedfirst(cum_P, rand())
    select_index_2 = select_index
    while select_index_2 == select_index
        select_index_2 = searchsortedfirst(cum_P, rand())
    end
    return select_index, select_index_2
end

function crossover(parent1::Vector{Int}, parent2::Vector{Int})
    """单点交叉操作"""
    r = rand(1:N_2)  # 确定交叉点
    offspring1 = vcat(parent1[1:r], parent2[r+1:end])
    offspring2 = vcat(parent2[1:r], parent1[r+1:end])
    return offspring1, offspring2
end

function mutate!(x::Matrix{Int}, mutation_rate::Float64)
    """变异操作"""
    mutation_mask = rand(size(x)...) .< mutation_rate
    x[mutation_mask] .= 1 .- x[mutation_mask]
end

function perform_evolution(x::Matrix{Int})
    """执行进化过程"""
    x2 = copy(x)  # 复制种群作为子代

    c_sum, w_sum = calculate_sums(x2)
    eval, cum_P = evaluate(c_sum, w_sum)

    select_index, select_index_2 = select_parents(cum_P)
    parent1 = x2[select_index, :]
    parent2 = x2[select_index_2, :]

    # 交叉操作
    if rand() > p_c
        offspring1, offspring2 = crossover(parent1, parent2)
        x2[select_index, :] = offspring1
        x2[select_index_2, :] = offspring2
    end

    # 变异操作
    mutate!(x2, p_m)

    # 计算新的c_sum,w_sum,eval
    c_sum2, w_sum2 = calculate_sums(x2)
    eval2, _ = evaluate(c_sum2, w_sum2)

    # 合并并排序选择
    eval0 = vcat(eval, eval2)
    x0 = vcat(x, x2)
    best_eval_indices = sortperm(eval0, rev=true)[1:N]  # 根据适应度排序选择前N个

    return x0[best_eval_indices, :]
end

function count_most_common(arr)
    """统计最常见元素及其出现次数"""
    counts = Dict{Int,Int}()
    for val in arr
        counts[val] = get(counts, val, 0) + 1
    end
    sorted_counts = sort(collect(counts), by=x->x[2], rev=true)
    most_common = first(sorted_counts)
    return most_common[1], most_common[2]
end

function run_simulation(runs::Int, target_value1::Int, target_value2::Int)
    """重复运行程序，测试程序稳定性并记录时间"""
    not_target_count = 0
    x = initialize_population(N, N_2)
    
    # 存储时间数据
    single_run_times = Float64[]
    cumulative_times = Float64[]
    total_start_time = time()
    
    for run_num in 1:runs
        run_start_time = time()
        
        for _ in 1:T
            x = perform_evolution(x)
        end
        
        # 计算所有个体的重量和价值
        w_sums = x * w_j
        c_sums = x * c_j
        
        # 找到出现频率最高的重量和价值
        most_common_wsum, most_common_count_w = count_most_common(w_sums)
        most_common_csum, most_common_count_c = count_most_common(c_sums)
        
        best_individual = x[argmax(c_sums), :]
        
        # 检查是否不等于目标值
        if most_common_wsum != target_value1 && most_common_csum != target_value2
            not_target_count += 1
        end
        
        # 记录时间
        push!(single_run_times, time() - run_start_time)
        push!(cumulative_times, time() - total_start_time)
        
        # 每100次打印进度
        if run_num % 100 == 0
            println("已完成 $run_num/$runs 次运行")
        end
    end
    
    # 绘制图表
    p1 = plot(1:runs, single_run_times, xlabel="Times", ylabel="Time(s)", 
              title="single", legend=false)
    
    p2 = plot(1:runs, cumulative_times, xlabel="Times", ylabel="cum_Time (s)", 
              title="cum", legend=false)
    
    plot(p1, p2, layout=(1,2), size=(1200,600))
    savefig("monitor/julia_performance.png")
    
    return not_target_count
end

function main()
    target_value1 = 879
    target_value2 = 1025
    runs = 10000
    
    println("开始运行遗传算法...")
    start_time = time()
    
    result = run_simulation(runs, target_value1, target_value2)
    rate = 100 * result / runs
    
    elapsed = time() - start_time
    println("测试次数：$runs, 错误次数: $result, 错误率: $(round(rate, digits=2))%")
    println("总运行时间：$(round(elapsed, digits=2)) 秒")
    println("平均每次运行时间：$(round(elapsed/runs, digits=4)) 秒")
end

main()