如何优化Python中傅里叶级数的计算耗时？

针对你当前的场景，除了已尝试的矩阵运算、并行循环，还有以下几种实用的优化方向：

1. 利用三角函数递推关系减少计算量 #

原函数中每次循环都要计算cos(angle*(i+1))和sin(angle*(i+1))，三角函数是相对耗时的操作。可以利用三角恒等式递推计算后续的三角函数值，避免重复调用cos/sin：

• cos((n+1)θ) = cos(nθ)cosθ - sin(nθ)sinθ
• sin((n+1)θ) = sin(nθ)cosθ + cos(nθ)sinθ

优化后的Numba函数示例：

import numba
import numpy as np

@numba.jit(nopython=True, fastmath=True)
def fourier_recursive(ai, bi, angle, nb_harmonique):
    serie = 0.0
    # 初始值：i=0对应(i+1)=1的三角函数值
    cos_prev = np.cos(angle)
    sin_prev = np.sin(angle)
    serie += ai[0] * cos_prev + bi[0] * sin_prev
    
    # 预存当前angle的cos和sin值，避免重复计算
    cos_theta = np.cos(angle)
    sin_theta = np.sin(angle)
    
    for i in range(1, nb_harmonique):
        # 递推计算第i+1次的三角函数值
        cos_curr = cos_prev * cos_theta - sin_prev * sin_theta
        sin_curr = sin_prev * cos_theta + cos_prev * sin_theta
        serie += ai[i] * cos_curr + bi[i] * sin_curr
        # 更新前一次的结果供下一轮使用
        cos_prev, sin_prev = cos_curr, sin_curr
    return serie

这个优化能将三角函数的调用次数从2*nb_harmonique次减少到2次，大幅降低单次调用的耗时。

2. 批量处理多个角度值 #

如果百万次调用是针对不同的angle，可以将所有angle整理成数组，结合批量并行计算，充分利用CPU的SIMD指令和多核资源：

@numba.jit(nopython=True, fastmath=True, parallel=True)
def fourier_batch(ai, bi, angles, nb_harmonique):
    results = np.empty_like(angles)
    # 预计算所有angle的cos和sin值
    cos_thetas = np.cos(angles)
    sin_thetas = np.sin(angles)
    
    # 并行处理每个angle
    for idx in numba.prange(len(angles)):
        serie = 0.0
        c_theta = cos_thetas[idx]
        s_theta = sin_thetas[idx]
        
        cos_prev = c_theta
        sin_prev = s_theta
        serie += ai[0] * cos_prev + bi[0] * sin_prev
        
        for i in range(1, nb_harmonique):
            cos_curr = cos_prev * c_theta - sin_prev * s_theta
            sin_curr = sin_prev * c_theta + cos_prev * s_theta
            serie += ai[i] * cos_curr + bi[i] * sin_curr
            cos_prev, sin_prev = cos_curr, sin_curr
        
        results[idx] = serie
    return results

使用时直接传入包含百万个角度的数组，比循环调用单次函数效率高得多，同时结合了递推优化的优势。

3. 降低数据精度（如果允许） #

如果业务场景对精度要求不高，可以将float64类型改为float32。CPU对32位浮点数的计算吞吐量远高于64位，能显著提升计算速度：

# 转换数据类型
ai_float32 = ai.astype(np.float32)
bi_float32 = bi.astype(np.float32)
angles_float32 = angles.astype(np.float32)

@numba.jit(nopython=True, fastmath=True)
def fourier_recursive_float32(ai, bi, angle, nb_harmonique):
    serie = 0.0
    cos_prev = np.cos(angle)
    sin_prev = np.sin(angle)
    serie += ai[0] * cos_prev + bi[0] * sin_prev
    
    cos_theta = np.cos(angle)
    sin_theta = np.sin(angle)
    
    for i in range(1, nb_harmonique):
        cos_curr = cos_prev * cos_theta - sin_prev * sin_theta
        sin_curr = sin_prev * cos_theta + cos_prev * sin_theta
        serie += ai[i] * cos_curr + bi[i] * sin_curr
        cos_prev, sin_prev = cos_curr, sin_curr
    return serie

注意：转换前需确认精度损失在可接受范围内。

4. 缓存重复角度的计算结果 #

如果百万次调用中存在重复的angle值，可以用字典缓存已计算的结果，避免重复计算：

from functools import lru_cache

# 注意：ai和bi如果是numpy数组，需要转成tuple才能被lru_cache缓存
def fourier_cached(ai, bi, angle, nb_harmonique):
    ai_tuple = tuple(ai)
    bi_tuple = tuple(bi)
    return _fourier_cached(ai_tuple, bi_tuple, angle, nb_harmonique)

@numba.jit(nopython=True, fastmath=True)
def _fourier_cached(ai_tuple, bi_tuple, angle, nb_harmonique):
    # 这里用递推实现的逻辑，和上面的fourier_recursive一致
    serie = 0.0
    cos_prev = np.cos(angle)
    sin_prev = np.sin(angle)
    serie += ai_tuple[0] * cos_prev + bi_tuple[0] * sin_prev
    
    cos_theta = np.cos(angle)
    sin_theta = np.sin(angle)
    
    for i in range(1, nb_harmonique):
        cos_curr = cos_prev * cos_theta - sin_prev * sin_theta
        sin_curr = sin_prev * cos_theta + cos_prev * sin_theta
        serie += ai_tuple[i] * cos_curr + bi_tuple[i] * sin_curr
        cos_prev, sin_prev = cos_curr, sin_curr
    return serie

这个方案仅适用于存在大量重复angle的场景，否则缓存的收益有限。

内容的提问来源于stack exchange，提问作者Clément