基于单遍算法的MATLAB偏度计算代码加速及C++扩展咨询

问题描述 #

我正在用单遍算法计算时间序列的偏度，MATLAB实现代码如下：

n = 10000000;
timeseries = randi(1000, 1, n);

tic
sk = skewness_onepass(timeseries);
toc

function skewness = skewness_onepass(x)
    
    N = length(x);
    M1 = 0;
    M2 = 0;
    M3 = 0;

    % 单遍算法基于高阶统计量递推逻辑
    for n = 1: N
        delta = x(n) - M1;
        M1 = M1 + delta / n;
        M3 = M3 + delta * delta * delta * (n - 1) * (n - 2) / (n * n) - 3 * delta * M2 / n;
        M2 = M2 + delta * delta * (n - 1) / n;
    end
    
    std_dev = sqrt(M2 / N);

    % 计算偏度
    skewness = (M3 / N) / (std_dev * std_dev * std_dev);
    
end

我的需求是尽可能加速这段代码，应用对时间要求极高。同时我不熟悉C++，如果用C++扩展MATLAB可行的话，需要具体实现代码和执行步骤，且希望这套方法能复用在其他类似统计函数上。

一、MATLAB原生优化（无需C++） #

先试试不用C++的优化，足够应对大部分场景：

• 减少循环内重复计算：把重复的乘法、除法提前计算，减少运算量，修改后的核心循环如下：

  for n = 1: N
      delta = x(n) - M1;
      delta_sq = delta * delta;
      delta_cu = delta_sq * delta;
      inv_n = 1 / n;
      inv_n_sq = inv_n * inv_n;
      
      M1 = M1 + delta * inv_n;
      M3 = M3 + delta_cu * (n-1)*(n-2)*inv_n_sq - 3*delta*M2*inv_n;
      M2 = M2 + delta_sq*(n-1)*inv_n;
  end
  

  这样避免了重复计算delta*delta、1/n等，减少了除法次数（除法比乘法慢）。
• 对比MATLAB内置函数：MATLAB自带的skewness函数是高度优化的底层实现，速度远快于手写循环，直接调用skewness(timeseries)就能获得最优性能，除非你必须用自定义的单遍逻辑。
• 启用JIT编译：确保MATLAB的JIT编译器开启（默认开启），把函数保存为独立的.m文件，不要写在脚本里，JIT会自动优化循环。

二、C++扩展实现（MEX文件） #

如果原生优化仍达不到要求，用MEX文件把核心逻辑改成C++，速度能提升数倍：

1. 编写C++代码 #

创建名为skewness_onepass.cpp的文件，内容如下：

#include "mex.h"
#include <cmath>

// MATLAB MEX接口函数
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    // 检查输入输出合法性
    if (nrhs != 1) {
        mexErrMsgIdAndTxt("skewness:badInput", "必须输入一个时间序列数组");
    }
    if (nlhs > 1) {
        mexErrMsgIdAndTxt("skewness:badOutput", "只能输出一个偏度值");
    }
    
    // 获取输入数组的指针和长度（注意MATLAB是1索引，C++是0索引）
    double *x = mxGetPr(prhs[0]);
    mwSize N = mxGetNumberOfElements(prhs[0]);
    
    // 初始化一阶、二阶、三阶矩
    double M1 = 0.0;
    double M2 = 0.0;
    double M3 = 0.0;
    
    // 单遍循环计算矩
    for (mwSize n = 1; n <= N; ++n) {
        double delta = x[n-1] - M1;
        double inv_n = 1.0 / n;
        double inv_n_sq = inv_n * inv_n;
        double delta_sq = delta * delta;
        double delta_cu = delta_sq * delta;
        
        M1 += delta * inv_n;
        M3 += delta_cu * (n-1)*(n-2)*inv_n_sq - 3*delta*M2*inv_n;
        M2 += delta_sq * (n-1)*inv_n;
    }
    
    // 计算标准差和偏度
    double std_dev = sqrt(M2 / N);
    double skewness_val = (M3 / N) / (std_dev * std_dev * std_dev);
    
    // 设置输出结果
    plhs[0] = mxCreateDoubleScalar(skewness_val);
}

2. 编译MEX文件 #

在MATLAB命令行中切换到代码所在目录，执行：

mex skewness_onepass.cpp

• Windows系统：需要提前安装Visual Studio（勾选C++开发组件）；
• Linux/macOS：需要安装GCC或Clang编译器，MATLAB会自动检测。

3. 调用编译后的函数 #

和调用普通MATLAB函数完全一致：

n = 10000000;
timeseries = randi(1000, 1, n);

tic
sk = skewness_onepass(timeseries);
toc

三、复用性说明 #

这套C++ MEX模板可以直接复用在其他单遍统计量计算（比如峰度、高阶矩）：

1. 保留mexFunction的框架，只修改统计量的初始化（比如计算峰度时增加M4 = 0.0）；
2. 替换循环内的矩更新公式（对应不同统计量的递推式）；
3. 最后修改输出的统计量计算逻辑即可。

这样核心的接口代码无需改动，只需要聚焦统计计算的核心逻辑。

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