线性递归向量化实现优化需求及Python与C向量化性能对比咨询

我需要对如下线性递归数学表达式进行向量化优化：

f(i) = f(i-1)/c + g(i)，其中i从1开始，f(0)、c为给定常数。

我当前用Python的列表推导式实现来提升速度，代码示例如下：

def function():
    txt= [0,2,0,2,0,2,0,2,2,2,0,2,2,2,0,2,0,2,2,2,0,2,0,2,0,2,0,2,2,2,0,2,0,2,0,2,0,2,2,2]
    indices_0=[]
    vl=0
    sb_l=10
    CONST=512
    # 初始计算部分
    [vl := vl+pow(txt[i],(i+1)) for i in range(sb_l)]
    if (vl==1876):
        indices_0=[0]
    # 核心递归滑动计算部分
    p=[i for i in range(1,len(txt)-sb_l+1) if (vl := (vl-txt[i-1])/2+ txt[i+sb_l-1]*CONST)==1876]
    print(indices_0+p)
function()

核心优化目标 #

核心需要向量化的部分是：

1. 初始累积计算：[vl := vl+pow(txt[i],(i+1)) for i in range(sb_l)]
2. 滑动递归计算：p=[i for i in range(1,len(txt)-sb_l+1) if (vl := (vl-txt[i-1])/2+ txt[i+sb_l-1]*CONST)==1876]
   其中递归式对应：f(i-1)= (vl-txt[i-1])，c=2，g(i)= txt[i+sb_l-1]*CONST。

额外问题 #

我当前用Python实现，请问采用C语言的向量化实现是否会大幅提升运行速度？是否存在比向量化更高效的实现方式？

一、Python环境下的向量化实现（基于NumPy） #

Python中最常用的向量化工具是NumPy，它能把循环操作转换为底层的C级数组运算，避免Python解释器的循环开销。

1. 初始累积计算的向量化 #

初始的pow(txt[i], i+1)可以用NumPy的广播机制实现，然后用cumsum完成累积求和：

import numpy as np

def vectorized_initial_calc(txt_np, sb_l):
    indices = np.arange(1, sb_l+1)  # i+1从1到sb_l
    pow_vals = np.power(txt_np[:sb_l], indices)
    vl = pow_vals.cumsum()[-1]  # 取最后一个累积和，等价于循环累加
    return vl

2. 滑动递归计算的向量化 #

你的递归式可以展开为滑动窗口的递推公式，我们可以把整个递推过程转换为数组运算：
观察递推式：
vl_new = (vl_old - txt[i-1])/2 + txt[i+sb_l-1] * CONST
我们可以先构造所有需要的txt[i-1]（对应txt[0 : len(txt)-sb_l]）和txt[i+sb_l-1]（对应txt[sb_l : ]），然后用累积运算模拟递推：

def vectorized_sliding_calc(initial_vl, txt_np, sb_l, CONST):
    # 构造递推的输入数组
    left_terms = txt_np[:len(txt_np)-sb_l]
    right_terms = txt_np[sb_l:] * CONST
    
    # 递推式可以改写为：vl[i] = vl[i-1]/2 + (right_terms[i-1] - left_terms[i-1]/2)
    # 这是一个线性递推，可以用NumPy的累积运算结合广播实现
    n_steps = len(left_terms)
    # 先计算每一步的增量项
    increments = right_terms - left_terms / 2
    # 构造递推的系数数组：(1/2)^k，k从0到n_steps-1
    coeffs = np.power(0.5, np.arange(n_steps, 0, -1))
    # 累积计算所有vl值
    vl_sequence = initial_vl * np.power(0.5, np.arange(1, n_steps+1)) + (increments * coeffs).cumsum()
    
    # 找到等于1876的索引（注意这里的索引对应原代码中的i，从1开始）
    matches = np.where(vl_sequence == 1876)[0] + 1  # 原代码中i从1开始
    return matches

把两部分结合起来，完整的向量化函数：

def vectorized_function():
    txt = [0,2,0,2,0,2,0,2,2,2,0,2,2,2,0,2,0,2,2,2,0,2,0,2,0,2,0,2,2,2,0,2,0,2,0,2,0,2,2,2]
    txt_np = np.array(txt, dtype=np.float64)
    sb_l = 10
    CONST = 512
    TARGET = 1876
    
    # 初始计算
    initial_vl = vectorized_initial_calc(txt_np, sb_l)
    indices_0 = [0] if initial_vl == TARGET else []
    
    # 滑动递推计算
    p = vectorized_sliding_calc(initial_vl, txt_np, sb_l, CONST).tolist()
    
    print(indices_0 + p)

二、比向量化更高效的Python实现：Numba JIT编译 #

如果你的递归逻辑比较复杂，NumPy的向量化可能需要额外的内存开销（比如构造系数数组），这时用Numba的JIT编译可以直接优化原生Python循环，性能接近C语言：

from numba import jit

@jit(nopython=True)  # 编译为机器码，无Python解释器开销
def numba_optimized_function():
    txt= [0,2,0,2,0,2,0,2,2,2,0,2,2,2,0,2,0,2,2,2,0,2,0,2,0,2,0,2,2,2,0,2,0,2,0,2,0,2,2,2]
    indices_0=[]
    vl=0.0
    sb_l=10
    CONST=512
    TARGET=1876
    
    # 初始计算
    for i in range(sb_l):
        vl += pow(txt[i], i+1)
    if vl == TARGET:
        indices_0=[0]
    
    # 滑动递推
    p = []
    max_i = len(txt)-sb_l
    for i in range(1, max_i+1):
        vl = (vl - txt[i-1])/2 + txt[i+sb_l-1]*CONST
        if vl == TARGET:
            p.append(i)
    
    print(indices_0 + p)

Numba会把这段代码编译成机器码，避免Python的循环开销，性能通常比NumPy向量化还要好，尤其是当数据量很大时。

三、C语言下的向量化实现及性能对比 #

1. C语言的向量化实现 #

在C语言中，你可以利用SIMD指令集（比如SSE、AVX）来实现向量化，或者用编译器的自动向优化（比如GCC的-O3 -mavx2编译选项）。以下是一个简化的C示例：

#include <stdio.h>
#include <stdint.h>
#include <math.h>

#define TARGET 1876
#define CONST 512

int main() {
    int txt[] = {0,2,0,2,0,2,0,2,2,2,0,2,2,2,0,2,0,2,2,2,0,2,0,2,0,2,0,2,2,2,0,2,0,2,0,2,0,2,2,2};
    int txt_len = sizeof(txt)/sizeof(int);
    int sb_l = 10;
    double vl = 0.0;
    int indices_0[1] = {-1};
    int p[100]; // 假设最多100个匹配项
    int p_count = 0;

    // 初始计算
    for(int i=0; i<sb_l; i++){
        vl += pow(txt[i], i+1);
    }
    if(vl == TARGET){
        indices_0[0] = 0;
    }

    // 滑动递推（编译器会自动向量化这个循环，加-O3选项）
    int max_i = txt_len - sb_l;
    for(int i=1; i<=max_i; i++){
        vl = (vl - txt[i-1])/2.0 + (double)txt[i+sb_l-1] * CONST;
        if(vl == TARGET){
            p[p_count++] = i;
        }
    }

    // 输出结果
    if(indices_0[0] != -1){
        printf("%d ", indices_0[0]);
    }
    for(int i=0; i<p_count; i++){
        printf("%d ", p[i]);
    }
    printf("\n");
    return 0;
}

编译时使用gcc -O3 -mavx2 your_code.c -o optimized_program，编译器会自动把循环转换为SIMD指令，实现向量化。

2. 性能对比 #

• Python原生列表推导式：最慢，因为每次循环都要经过Python解释器，还有:=赋值的额外开销。
• NumPy向量化：比原生快10-100倍，取决于数据量大小，底层是C实现的数组运算，但需要额外的内存存储中间数组。
• Numba JIT：性能接近C语言，通常比NumPy快，尤其是当循环逻辑复杂时，不需要额外内存。
• C语言向量化：比Python所有实现都快，通常是Numba的1.2-2倍左右，尤其是当数据量非常大时，SIMD指令的优势更明显。

所以，如果你的代码需要处理超大数据集，C语言的向量化实现会有大幅的速度提升；如果数据集中等，Numba JIT已经足够满足性能需求，同时保留Python的开发便捷性。

内容的提问来源于stack exchange，提问作者Michael