3

我被分配了一个任务，我应该在其中编写一个算法，该算法通过重心公式执行多项式插值。公式指出：

p(x) = (SIGMA_(j=0 到 n) w(j)*f(j)/(x - x(j)))/(SIGMA_(j=0 到 n) w(j)/(x - x(j)))

我编写了一个运行良好的算法，并且得到了我想要的多项式输出。然而，这需要使用一些相当长的循环，并且对于大的网格数，将不得不进行许多讨厌的循环操作。因此，如果有人对我如何改进这一点有任何提示，我将不胜感激，这样我就可以避免所有这些循环。

在算法中，x代表f我们应该插值的给定点。 w代表重心权重，在运行算法之前已经计算过。并且grid是应该进行插值的 linspace：

function p = barycentric_formula(x,f,w,grid)

%Assert x-vectors and f-vectors have same length.
if length(x) ~= length(f)
    sprintf('Not equal amounts of x- and y-values. Function is terminated.')
    return;
end

n = length(x);
m = length(grid);
p = zeros(1,m);

% Loops for finding polynomial values at grid points.  All values are
% calculated by the barycentric formula.
for i = 1:m
    var = 0;
    sum1 = 0;
    sum2 = 0;
    for j = 1:n
        if grid(i) == x(j)
            p(i) = f(j);
            var = 1;
        else
            sum1 = sum1 + (w(j)*f(j))/(grid(i) - x(j));
            sum2 = sum2 + (w(j)/(grid(i) - x(j)));
        end
    end
    if var == 0
        p(i) = sum1/sum2;
    end    
end
4

1 回答 1

4

这是matlab“矢量化”的经典案例。我会说 - 只需删除循环。几乎就是这么简单。首先，看一下这段代码：

function p = bf2(x, f, w, grid)

m = length(grid);
p = zeros(1,m);

for i = 1:m
    var = grid(i)==x;
    if any(var)
        p(i) = f(var);
    else
        sum1 = sum((w.*f)./(grid(i) - x));
        sum2 = sum(w./(grid(i) - x));
        p(i) = sum1/sum2;
    end
end
end

我已经删除了内部循环j。实际上，我在这里所做的只是删除(j)索引并将算术运算符从/to./和 from *to 更改.*- 相同，但前面有一个点表示该操作是逐个元素执行的。与普通的矩阵运算符相比，这称为数组运算符。另请注意，处理网格点所在的特殊情况x与您在原始实现中的情况非常相似，只使用了var这样的向量x(var)==grid(i)。

现在，您还可以删除最外层的循环。这有点棘手，在 MATLAB 中有两种主要方法可以做到这一点。我会用更简单的方法来做，它可能效率较低，但阅读起来更清晰 - 使用repmat：

function p = bf3(x, f, w, grid)

% Find grid points that coincide with x.
% The below compares all grid values with all x values
% and returns a matrix of 0/1. 1 is in the (row,col)
% for which grid(row)==x(col)

var  = bsxfun(@eq, grid', x);

% find the logical indexes of those x entries
varx = sum(var, 1)~=0;

% and of those grid entries
varp = sum(var, 2)~=0;

% Outer-most loop removal - use repmat to
% replicate the vectors into matrices.
% Thus, instead of having a loop over j
% you have matrices of values that would be
% referenced in the loop

ww = repmat(w, numel(grid), 1);
ff = repmat(f, numel(grid), 1);
xx = repmat(x, numel(grid), 1);
gg = repmat(grid', 1, numel(x));

% perform the calculations element-wise on the matrices
sum1 = sum((ww.*ff)./(gg - xx),2);
sum2 = sum(ww./(gg - xx),2);
p    = sum1./sum2;

% fix the case where grid==x and return
p(varp) = f(varx);

end

完全矢量化的版本可以使用bsxfun而不是repmat. 这可能会更快一些，因为矩阵不是明确形成的。但是，对于小型系统，速度差异可能不会很大。

此外，第一个带有一个循环的解决方案在性能方面也不算太差。我建议你测试一下，看看哪个更好。也许完全矢量化不值得？第一个代码看起来更具可读性..

于 2012-10-24T18:37:50.980 回答