CODEWITHSUNDEEP
matlabvectorization
romanbird
romanbird

Reputation: 178

Matlab function with two vector arguments

I have a function defined as follows:

@(A,N0,Tb,lambda,p) p.*erfc(1.0./sqrt(N0).*sqrt(Tb).*(A-lambda)).*(1.0./2.0)-erfc(1.0./sqrt(N0).*sqrt(Tb).*(A+lambda)).*(p-1.0).*(1.0./2.0)

I have N0 and lambda as vectors. If I pass N0 as a vector and lambda as a scalar, I get a vector output.

For each lambda, I want to evaluate the function with the full N0 vector. (i.e., I want to return an array of vectors instead of just one vector). Is there a way to do this without looping?

Upvotes: 1

Views: 185

Answers (3)

knedlsepp
knedlsepp

Reputation: 6084

I will assume your data is given by:

  • A, Tb, p scalar values
  • N0, lambda column-vectors.

Version 1 - (fast and specific to the problem)

As the two vectors you are trying to combine are combined via multiplication, you can make use of linear algebra. If you multiply a column vector by a row vector, you will get the matrix of pairwise products:

>> (1:4)'*(10:12)
ans =
    10    11    12
    20    22    24
    30    33    36
    40    44    48

So your code would result in:

f = @(A,N0,Tb,lambda,p) p*erfc(sqrt(Tb./N0)*(A-lambda).')*0.5 ...
                   +(1-p)*erfc(sqrt(Tb./N0)*(A+lambda).')*0.5

Version 2 - (more general and reusing code)

If your formula were to combine the values by something else than multiplication, you couldn't go with linear algebra. For this there is the marvelous function bsxfun, which combines the column and row vectors in a similar way, but you can specify the function to use instead of multiplication. You can use it here in the following way:

Let's call your initial function f, you could then define a new function g, that handles lambda and N0 pairwise with this:

g = @(A,N0,Tb,lambda,p) bsxfun(@(N0,lambda) f(A,N0,Tb,lambda,p), N0, lambda')

Upvotes: 2

Divakar
Divakar

Reputation: 221514

Some discussion and code

Proposed here is another bsxfun based approach tailor-made for the specific function handle used in the problem -

%// This could could be reused in the function definition, 
%// so that pre-calulcating it makes sense for performance
f1vals = 1.0./sqrt(N0).*sqrt(Tb);

%// Finally get the output using f1vals and other inputs
out = p.*erfc(bsxfun(@times,f1vals,A - lambda')).*0.5 - ...
                     erfc(bsxfun(@times,f1vals,A + lambda')).*(p-1.0).*0.5;

Please note that this approach also assumes N0 and lambda as column-vectors and others as scalars.

Benchmarking

The proposed solution code was run against @knedlsepp's Version -1 solution code and @rayryeng's solution code for various datasizes. The benchmarking code is listed here on ideone.com. The runtimes obtained on my system were -

1) N0 and lambda of lengths 500:

--------- With Proposed solution
Elapsed time is 0.068731 seconds.
--------- With Knedlsepp Version - 1 solution
Elapsed time is 0.082629 seconds.
--------- With Rayryeng solution
Elapsed time is 0.132094 seconds.

2) N0 and lambda of lengths 1000:

--------- With Proposed solution
Elapsed time is 0.220694 seconds.
--------- With Knedlsepp Version - 1 solution
Elapsed time is 0.244323 seconds.
--------- With Rayryeng solution
Elapsed time is 0.334007 seconds.

3) N0 and lambda of lengths 3500:

--------- With Proposed solution
Elapsed time is 2.571570 seconds.
--------- With Knedlsepp Version - 1 solution
Elapsed time is 3.009692 seconds.
--------- With Rayryeng solution
Elapsed time is 3.747122 seconds.

Upvotes: 1

rayryeng
rayryeng

Reputation: 104464

Yes you can. You'll need to do a bit of pre-processing before you use your function handle. I'm going to assume that N0 and lambda are both column vectors where N0 is M elements long and lambda is N elements long. If you want to create a 2D matrix such that each row is the result of applying a single value of lambda while N0 remains constant, what you'll need to do is make both N0 and lambda a matrix such that they are both N x M long. N0 will be transposed such that it will be a row vector, and stacked row-wise N times to create this matrix. lambda will also be N x M, where the lambda vectors gets stacked column-wise M times to create this matrix.

You would then submit these matrices as input into your function handle. The result would be a N x M matrix, where each column is the result of applying one lambda to the vector N0. As proof of this concept, consider the first part of your function handle:

p.*erfc(1.0./sqrt(N0).*sqrt(Tb) ...

The result of this step will have M identical columns produced from the erfc function. We then multiply this with (A - lambda). lambda will be shaped in such a way where each row reflects a single lambda to be applied to every single value of N0. The same can be said about the other side of the equation with (A + lambda). As such, try doing this:

M = numel(N0);
N = numel(lambda);
N0_matrix = repmat(N0.', N, 1);
lambda_matrix = repmat(lambda, 1, M);
f = @(A,N0,Tb,lambda,p) p.*erfc(1.0./sqrt(N0).*sqrt(Tb).*(A-lambda)).*(1.0./2.0)-erfc(1.0./sqrt(N0).*sqrt(Tb).*(A+lambda)).*(p-1.0).*(1.0./2.0);
out = f(A, N0_matrix, Tb, lambda_matrix, p);

out will be that N x M matrix we talked about.

Upvotes: 2

Related Questions