MATLAB中如何向量化自定义函数脚本？含实现代码求助

首先，你的核心问题是原函数不支持向量/矩阵输入，导致bsxfun调用时出现维度不匹配。我们从底层到顶层逐步改造，让函数能处理批量输入，同时摆脱嵌套循环：

步骤1：先让Integrand完全支持向量化 #

原Integrand是标量级运算，直接修改为支持任意维度的输入（只要输入维度一致），利用MATLAB原生的向量化算术运算：

function ret = Integrand(r0, sigma, lambda, t, tau, x, y)
    % 预计算重复出现的项，减少冗余计算
    lambda_t = lambda .* t;
    lambda_tau = lambda .* tau;
    lambda_t_tau = lambda .* (t - tau);
    exp_lambda_t = exp(-lambda_t);
    exp_lambda_tau = exp(-lambda_tau);
    exp_2lambda_t_tau = exp(-2 * lambda_t_tau);
    
    sinh_term = sinh(lambda .* (tau - t));
    cosh_term = cosh(lambda .* (tau - t));
    
    % 拆分分子部分
    part1 = (lambda .* r0 .* exp_lambda_t) ./ 2;
    part2 = ((x - r0 .* exp_lambda_t) ./ 2) .* lambda .* cosh_term ./ sinh_term;
    part3 = ((y - r0 .* exp_lambda_tau) ./ 2) .* lambda ./ sinh_term;
    numerator = part1 + part2 - part3;
    
    % 拆分高斯分布部分
    sqrt_term = sqrt(pi .* sigma.^2 ./ lambda .* (1 - exp_2lambda_t_tau));
    exponent_arg = ((x - r0 .* exp_lambda_t - exp(-lambda_t_tau) .* (y - r0 .* exp_lambda_tau)).^2) ...
        ./ (sigma.^2 ./ lambda .* (1 - exp_2lambda_t_tau));
    gauss_part = exp(-exponent_arg) ./ sqrt_term;
    
    ret = numerator .* gauss_part;
end

修改后，Integrand可以直接处理向量/矩阵输入，无需依赖arrayfun。

步骤2：替换m循环的权重计算（向量化条件判断） #

原函数中m循环是逐个计算权重，我们把m做成向量，用逻辑索引替代逐次if判断，完全去掉这个循环：

% 替换原m循环的代码块：
k_current = k;
m_vec = 1:(k_current-1);

% 预计算所有条件变量
R1 = mod(k_current, 2);
Q1 = floor(k_current / 2);
R2 = mod(m_vec, 2);
Q2 = floor(m_vec / 2);

% 初始化权重数组并通过逻辑索引赋值
weight = zeros(1, length(m_vec));
weight(R1==0 & R2==1 & Q2<=Q1-1) = 4/3;
weight(R1==0 & R2==0 & Q2<=Q1-2) = 2/3;
weight(R1==1 & R2==1 & Q2<=Q1-2) = 4/3;
weight(R1==1 & R2==0 & Q2<=Q1-3 & Q1>=3) = 2/3;
weight(R1==1 & R2==0 & m_vec==2*(Q1-1)) = 17/24;
weight(R1==1 & m_vec==2*Q1-1 & Q1>=1) = 9/8;
weight(R1==1 & m_vec==2*Q1 & Q1>=1) = 9/8;

这种向量操作的效率远高于逐元素循环。

步骤3：改造k循环为多维度矩阵操作 #

原k循环是逐列更新g，我们把g改成3维数组（对应r0元素、t0元素、k迭代次数），利用MATLAB的广播机制处理不同维度的输入：

function density = HittingDensityOUlevel0(beta, r0, sigma, lambda, t0, i, p, MinOrMax)
    % 统一输入维度，方便广播
    r0 = reshape(r0, [], 1);
    t0 = reshape(t0, 1, []);
    
    % 初始化3维g数组：行=r0数量，列=t0数量，层=k迭代次数
    g = zeros(numel(r0), numel(t0), i);
    
    % 处理k=1的初始情况
    g(:,:,1) = -MinOrMax .* 2 .* Integrand(r0, sigma, lambda, t0+p, t0, beta, r0);
    
    % 处理k>=2的迭代
    for k = 2:i
        m_vec = 1:(k-1);
        R1 = mod(k, 2);
        Q1 = floor(k / 2);
        R2 = mod(m_vec, 2);
        Q2 = floor(m_vec / 2);
        
        % 生成权重数组
        weight = zeros(1, k-1);
        weight(R1==0 & R2==1 & Q2<=Q1-1) = 4/3;
        weight(R1==0 & R2==0 & Q2<=Q1-2) = 2/3;
        weight(R1==1 & R2==1 & Q2<=Q1-2) = 4/3;
        weight(R1==1 & R2==0 & Q2<=Q1-3 & Q1>=3) = 2/3;
        weight(R1==1 & R2==0 & m_vec==2*(Q1-1)) = 17/24;
        weight(R1==1 & m_vec==2*Q1-1 & Q1>=1) = 9/8;
        weight(R1==1 & m_vec==2*Q1 & Q1>=1) = 9/8;
        
        % 生成TimeTicker，支持批量t0
        TimeTicker = t0 + p * (1:(k-1));
        
        % 批量计算Integrand值
        integrand_vals = Integrand(r0, sigma, lambda, t0 + k*p, TimeTicker, beta, beta);
        
        % 计算Sum，利用3维数组求和
        Sum = 2 * p * sum(weight .* g(:,:,1:(k-1)) .* integrand_vals, 3);
        
        % 更新当前k层的g值
        g(:,:,k) = -MinOrMax .* 2 .* Integrand(r0, sigma, lambda, t0 + k*p, t0, beta, r0) + MinOrMax .* Sum;
    end
    
    % 取最后一层作为结果，返回r0与t0所有组合的密度值矩阵
    density = g(:,:,end);
end

步骤4：调用改造后的函数（无需bsxfun） #

现在函数已经支持批量r0和t0输入，直接传入向量即可，MATLAB会自动广播处理所有组合：

% 替换为你的实际参数
beta = ...;
sigma = ...;
lambda = ...;
i = ...;
p = ...;
MinOrMax = ...;

t0 = [1,2,3];
r0 = [4,5];

% 调用函数，返回2x3矩阵：density(j,k)对应r0(j)与t0(k)的组合结果
density_matrix = HittingDensityOUlevel0(beta, r0, sigma, lambda, t0, i, p, MinOrMax);

额外优化建议 #

• 如果你的MATLAB版本低于R2016b（不支持隐式广播），可以用bsxfun辅助处理维度扩展，确保输入维度匹配。
• 若i（迭代次数）很大，可考虑将k循环进一步向量化，但通常i是数值近似的阶数，不会特别大，当前循环的开销可接受。
• 预计算所有重复的指数、双曲函数项，能进一步提升计算效率。

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