向量化涉及循环中的while循环或if子句的函数（Matlab）答案

【问题标题】：Vectorizing a function involving a while loop or if-clause in a loop (Matlab)向量化涉及循环中的while循环或if子句的函数（Matlab）
【发布时间】：2012-03-22 11:29:12
【问题描述】：

假设我有一个可以从一个输入计算一个输出的函数，例如

function y = sqrt_newton(x)
    y = x ./ 2;
    yo = y;
    y = 0.5.*(y + x ./ y);
    while abs(y - yo) > eps * abs(y)
        yo = y;
        y = 0.5.*(y + x ./ y);
    end
end

我希望能够将此函数应用于矢量输入，例如 sqrt_newton(2:9)，就像使用内置函数一样。在循环开始时使用条件或内部某些 if 子句来实现此目的的最佳方法是什么？如果可能的话，我想避免编写一个额外的函数作为包装器来循环输入向量。

我目前的繁琐解决方案

到目前为止我所做的是：

我必须首先将输入扩展为相同的大小（使用金融工具箱中的 finargsz，但如果您知道另一个执行相同功能的核心函数，那就太好了）
使用size记录形状
deal 输入
遍历所有输入元素
reshape输出

似乎numel 函数减轻了所有这些繁重工作的需要，但额外的 cmets 将是最受欢迎的。

【问题讨论】：

标签： arrays matlab loops conditional-statements vectorization

【解决方案1】：

总是有arrayfun。您可以保留已有的代码，将其放入内部函数中。

function y = sqrt_newton(z)

    y = arrayfun(@inner, z);

    function y= inner(x)
        y = x ./ 2;
        yo = y;
        y = 0.5.*(y + x ./ y);
        while abs(y - yo) > eps * abs(y)
            yo = y;
            y = 0.5.*(y + x ./ y);
        end
    end
end

编辑：上述解决方案的优点是在使用 1x1 输入后实现起来很简单，但其他答案中的循环对于大型输入来说要快得多。比如在我的电脑上，代码

tic; sqrt_newton(rand(500)); toc

使用我的代码在~1.24 seconds 运行，0.06 seconds 使用@Ramashalanka 的代码运行，0.28 seconds 使用@GuntherStruyf 的代码运行。

【讨论】：

嗯，我喜欢您的解决方案，因为它相对于原始函数的开销最小。太糟糕了，它是如此缓慢。

【解决方案2】：

我认为通常你必须使用循环，因为函数操作的未知字符。如果是线性运算，则可以进行矢量化。

对于您的示例，我将使用以下内容：

function y = sqrt_newton(x)
    y = x ./ 2;
    yo = y;
    y = 0.5.*(y + x ./ y);
    for i=1:numel(x)
        while abs(y(i) - yo(i)) > eps * abs(y(i))
            yo(i) = y(i);
            y(i) = 0.5*(y(i) + x(i) / y(i));
        end
    end
end

我使用 numel 而不是 size，所以它可以处理我扔给它的任何数组

【讨论】：

事实上，我使用了类似的解决方案，但我期待其他方式可以巧妙地做到这一点。请参阅我在问题下的额外评论。

【解决方案3】：

好吧，您可以将其矢量化如下（使用any）：

function y = sqrt_newton(x)
    y = x / 2;
    yo = y;
    y = 0.5*(y + x ./ y);
    while any(abs(y - yo) > eps * abs(y))
        yo = y;
        y = 0.5*(y + x ./ y);
    end
end

然后你得到：

>> sqrt_newton(2:9)
ans =
    1.4142    1.7321    2.0000    2.2361    2.4495    2.6458    2.8284    3.0000

>> ans.^2-(2:9)
ans =
   1.0e-14 *
   -0.0444   -0.0444         0    0.0888   -0.0888    0.0888   -0.1776         0

正如预期的那样。但是，我不推荐它，因为您正在对已经收敛的元素进行不必要的操作。我只是在函数的开头使用for 循环x：

function yall = sqrt_newton(xall)
yall = zeros(size(xall));
for xn=1:numel(xall)
    x = xall(xn);
    y = x / 2;
    yo = y;
    y = 0.5*(y + x ./ y);
    while abs(y - yo) > eps * abs(y)
        yo = y;
        y = 0.5*(y + x ./ y);
    end
    yall(xn)=y;
end
end

在开始时设置大小yall，以避免它在整个循环中增大。

【讨论】：