CODEWITHSUNDEEP
arraysmatlabloopsconditional-statementsvectorization
green diod
green diod

Reputation: 1499

Vectorizing a function involving a while loop or if-clause in a loop (Matlab)

Let's say I have a function that can compute one output from one input, e.g.

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

I'd like to be able to apply this function to a vector input say sqrt_newton(2:9) like with built-in functions. What is the best way to achieve this with a condition at the beginning of the loop or some if-clause inside? I'd like to avoid writing an extra function as a wrapper just to loop through the input vector, if possible at all.

My current cumbersome solution

What I do up until now is:

It seems that the numel function alleviates the need of all this heavy lifting but extra comments would be most welcome.

Upvotes: 3

Views: 463

Answers (3)

Ian Hincks
Ian Hincks

Reputation: 4118

There's always arrayfun. You can keep the code you have, putting it into an inner function.

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

Edit: The above solution has the advantage of being trivial to implement after you have it working with a 1x1 input, but the loops in the other answers are way faster for large inputs. For example, on my computer, the code

tic; sqrt_newton(rand(500)); toc

runs in ~1.24 seconds with my code, 0.06 seconds with @Ramashalanka's code, and 0.28 seconds with @GuntherStruyf's code.

Upvotes: 2

Ramashalanka
Ramashalanka

Reputation: 8854

Well, you could vectorize it as follows (with 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

Then you get:

>> 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

as expected. However, I wouldn't recommend it, since you're doing unnecessary operations on elements that have already converged. I'd just use a for loop over x at the start of the function:

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

Set the size yall at the start to avoid it increasing in size throughout the loop.

Upvotes: 2

Gunther Struyf
Gunther Struyf

Reputation: 11168

I think in general you'll have to use a loop, because of the unknown character of the operation of the function. In case it's a linear operation, vectorization is possible.

For your example I'd use the following:

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

I use numel, instead of size, so it can handle any array I throw at it

Upvotes: 2

Related Questions