arrayfun with multiple inputs and a different size output

I'm trying to set arrayfun to compute Monte-Carlo simulations on a GPU. But first, I'd like to parameterize on a CPU (I don't have GPUs yet).

The thing is, arrayfun takes as inputs arrays of same size, and then returns a scalar as said here (second setup): https://fr.mathworks.com/help/matlab/ref/arrayfun.html

First issue: myfunc is not returning a scalar but a big matrix for each Monte-Carlo simulation. How must the setup be done ?

Second issue: My inputs are Px1 vectors, which are reused M times. I also need to pass an input matrix A (also reused for each simulation), which clearly won't be the same size as other inputs Px1. I'm not sure what to do. Any help would be very appreciated.

FYI: myfunc is returning a N x P matrix

The aim is to have M simulations, therefore I'd like arrayfun to run M times to have as a final result a N x P x M matrix. Final goal is to prepare the code for GPU enhancement.

Answer 1

[Guillaume]

Note: I'm not familiar enough with GPU computing to know if there are some restrictions that apply.

First issue: myfunc is not returning a scalar but a big matrix for each Monte-Carlo simulation. How must the setup be done

Simply tell arrayfun to put the non-scalar output in a cell array using 'UniformOutput', false:

arrayfun(@somefunc, in1, in2, in3, 'UniformOutput', false)

Second issue: I also need to pass an input matrix A

Create an aonymous function which passes A to the real function. This technique is called argument binding:

arrayfun(@(in1, in2, in3) somefunc(A, in1, in2, in3), in1, in2, in3, 'UniformOutput', false)

The anonymous function @(in1, in2, in3) somefunc(A, in1, in2, in3) binds A to the somefunc call. A must exists before the anonymous function is created and once the function is created changing the value of A (or clearing it) will not affect the A that is bound.

edit: accidently deleted my answer. Restored now

Answer 2

[Adam]

c = arrayfun( @(x,y) rand( x, y ), 1:4, 3:6, 'UniformOutput', false )

works fine to produce outputs that are matrices. If your matrix A does not change then you can just pass this as a bound argument to your function e.g.

A = rand(2,2);
c = arrayfun( @(x,y) rand( x, y ) + A(randi(4)), 1:4, 3:6, 'UniformOutput', false )

Apologies for the ridiculously contrived and meaningless example. I was just looking for the quickest way I could find to show the syntax!

If your A changes with each of the other inputs you would probably have to put all those matrices, A, into a cell array and use cellfun with 3 inputs instead.

Answer 3

[cedric W]

When using both your answers, I'm not sure this is doing what I want.

First, without A to start easy. Say I'm starting from:

c = arrayfun( @(x,y) rand( x, y ), 1:4, 3:6, 'UniformOutput', false )

So in that example my dimension 'P' above is equal to 4 (1:4 = 4 rows)

What I would like to have is a function, which gives me M times a result of size N (just a number passed to arrayfun as an input) x 4 (the dimension 4 is equivalent to the size of the other Px1 input vectors). So this is a 3D matrix/cell array. I can't see the 'M' component in your solutions.

If you are familiar with monte carlo, M is the number of simulations, N is the number of time steps, P is the number of correlated assets

My inputs are constant, they're not dependant on the simulation number. It's in myfunc that there're random variables which in the end output different N x P matrices

[Adam]

You can pass any constants you want into your function, as shown in both answers with A. It is only the variable arguments, those in the brackets at the front of the anonymous function handle that have size constraints on them as they will be filled by each pair of elements of your two input arrays, hence why they must obviously be the same size.

You can pass anything else you want (e.g. M) to the set of arguments in the brackets after the function name as bound arguments.

[cedric W]

Ok I understand a little more. Except that I have maybe 20 constants precalculated, this will be a long call.

One last point, is that I don't understand how to call the function M times ? Because I don't have any dependency on M in myfunc

[Adam]

You will likely (well, definitely by the sound of it) need to create your function as a standalone function, rather than an inline anonymous function. Then your function can contain for loops based on the input M if you wish.

With regard to your constants, I would put them in a struct if you have that many to pass in. An array would be neater in some senses if they are all scalars, but for referring to them afterwards in code that wouldn't be a friendly solution whereas with a struct you can pass one object but each constant can retain its meaningful name as a field of the struct.

Guillaume's example shows the usage of some function defined elsewhere. somefunc can be whatever you want it to be.

Though I would qualify all this advise, like Guillaume did, by saying that this works for standard CPU usage. I am also unfamiliar with what limitations there may be with GPU usage.

[cedric W]

The point is that arrayfun will then run on a GPU, so doing a for loop is not efficient. On a GPU it would be automatically dispatched, I guess.

myfunc so far is a standalone as you said. What if I pass a fake variable of size Mx1 to bluff the computer so that it thinks it has to run M times ?

[Guillaume]

arrayfun loops over the elements of the input array(s). Nothing more. Passing a fake M variable won't help. All the inputs to arrayfun must be the same size. Possibly, you could repmat all your input arrays by M but in my opinion you're trying to couple two indenpendant concepts in one arrayfun which is not a good idea. They should really stay separate.

Of course, you can always wrap your arrayfun into another arrayfun giving it a Mx1 input that is unused:

arrayfun(@(~) arrayfun(@(in1, in2, in3) somefunc(A, in1, in2, in3), in1, in2, in3, 'UniformOutput', false), 1:M, 'UniformOutput', false)

As said, I'm not familiar enough with GPU computing, don't have the toolbox, but I don't believe the above would result in significantly different code than using for loops which would be less contrived.

[Guillaume]

Note: I'm not familiar enough with GPU computing to know if there are some restrictions that apply.

First issue: myfunc is not returning a scalar but a big matrix for each Monte-Carlo simulation. How must the setup be done

Simply tell arrayfun to put the non-scalar output in a cell array using 'UniformOutput', false:

arrayfun(@somefunc, in1, in2, in3, 'UniformOutput', false)

Second issue: I also need to pass an input matrix A

Create an aonymous function which passes A to the real function. This technique is called argument binding:

arrayfun(@(in1, in2, in3) somefunc(A, in1, in2, in3), in1, in2, in3, 'UniformOutput', false)

The anonymous function @(in1, in2, in3) somefunc(A, in1, in2, in3) binds A to the somefunc call. A must exists before the anonymous function is created and once the function is created changing the value of A (or clearing it) will not affect the A that is bound.

[Guillaume]

As for the many constants another option to the struct or vector of constants is to wrap the function in a class (a function object):

classdef myfunc
   properties
      const1;
      const2;
      const3;
   end
   methods
      function this = myfunc(c1, c2, c3)  %constructor
         this.const1 = c1;
         this.const2 = c2;
         this.const3 = c3;
      end
      function out = run(this, A, in1, in2, in3)
          out = A * this.const1 * in1 + this.const2 * in2 + this.const3 * in3;
      end
   end

Used:

func = myfunc(constant1, constant2, constant3);  %set the constants of the function
arrayfun(@(A, in1, in2, in3) func.run(A, in1, in2, in3), in1, in2, in3, 'UniformOutput', false);
func.const3 = newvalue;  %change value of constant3
arrayfun(@(A, in1, in2, in3) func.run(A, in1, in2, in3), in1, in2, in3, 'UniformOutput', false);  %run with new value (others unchanged)

[cedric W]

Great, thanks a lot to both of you. I'll do something like that then.