Loren on the Art of MATLAB

Turn ideas into MATLAB

Use nested functions to memoize costly functions 27

Posted by Loren Shure,

Some languages provide ways to store results of costly calculations so if the calculation had been performed before, it need not be repeated when requested later. There is a way to accomplish this in MATLAB as well, using nested functions. Consider this function:
function f = memoize1(F)
%one-arg F, inputs testable with ==

x = [];
y = [];
f = @inner;
    function out = inner(in)
        ind = find(in == x);
        if isempty(ind)
            out = F(in);
            x(end+1) = in;
            y(end+1) = out;
            out = y(ind);
Let's see what happens when we use this function, admittedly with not such an expensive function to calculate. Also, you could substitute any other function handle here, including one referring to an anonymous function, provided it has one input and one output.
>> f = memoize1(@sin)
f =
>> x = f(3)
x =
>> y = f(7)
y =
>> z = f(3)
z =
>> f(.2)
ans =
>> f(pi/8)
ans =
Now let's look into what's happening with our function f using the functions function. This should be used primarily for debugging and gaining insight, not as a programmatic tool.

Get information about the function handle:

>> s = functions(f)
s =
function: 'memoize1/inner'
type: 'nested'
file: 'H:DocumentsLORENv7memoize1.m'
workspace: {[1x1 struct]}

Get info about the nested function including its workspace:

>> s.workspace{1}
ans =

f: @memoize1/inner
F: @sin
x: [3 7 0.2000 0.3927]
y: [0.1411 0.6570  0.1987 0.3827]
So, if we were to put a breakpoint at inner and issued the command f(3) from the prompt, we'd be able to follow the logic about indexing into the saved values rather than doing the calculation again. In addition, if we wanted to quit MATLAB and come back at a later time, we could save our function handle in a MAT-file and reload it later. We would not need to save anything in addition to the handle, since the associated workspace would be also be saved. memoize1 can clearly be improved by vectorizing it (probably using ismember). Here is the latest version of this variant of the code:

function f = memoize2(F)
%one-arg F, inputs testable with ==
%allow nonscalar input.

x = [];
y = [];
f = @inner;
    function out = inner(in)
        out = zeros(size(in));  % preallocate output
        [tf,loc] = ismember(in,x);  % find which in's already computed in x
        ft = ~tf;  % ones to be computed
        out(ft) = F(in(ft));  % get output values for ones not already in
        % place new values in storage
        x = [x in(ft(:).')];
        y = [y reshape(out(ft),1,[])];
        out(tf) = y(loc(tf));  % fill in the rest of the output values

Variants of this memoize code pattern might be useful for functions that have recurrence relationships, recursive functions, or refining underlying meshes. It's possible that you would want to have a buffer limit so the size of the stored data in the nested function didn't grow beyond some bound. But you get the idea.


Comments are closed.

27 CommentsOldest to Newest

Markus replied on : 3 of 27
Hi! I am a little bit confused why your function memoize1 does what it does. I was used that I need a "persistent" declaration when I want a function to keep variable values from one call to another. Is there a special behaviour when using function handles to nested functions? It does not seem to be consistent/straightforward to me at the moment. Regards Markus
Loren replied on : 4 of 27
With persistent, you only have one instance because all invocations of a function with a persistent variable share the value (assuming regular and not nested functions). With nested functions, each instance of the function handle carries along its own workspace. I recommend you read the user documentation carefully and examine the demos and examples to understand more (especially nesteddemo).
Rod replied on : 5 of 27
I have tried to use your memoize2 function and it doesn't work for me. The problem comes from the function ismember. Below are three easy tests which should give the same answer: >> ismember(10.*(0:0.1:1),10.*(0:0.1:2)) ans = 1 1 1 1 1 1 0 0 1 1 1 >> ismember(0:10,0:20) ans = 1 1 1 1 1 1 1 1 1 1 1 >> ismember(0:10,10.*(0:0.1:2)) ans = 1 1 1 0 1 1 0 0 1 1 1 Whereas the two sets of data are the same (or should be same) for the three examples I don't get the same answer! Any idea what's wrong here? Thanks
Markus replied on : 7 of 27
The winning contribution of the last Matlab programming contest is a bad example for using nested functions. The functions solveri and CalculateMoves are nested functions (this is better visible after correct identing). Both functions are quite large and use quite a number of variables. After starting to modify the code I quickly ran into the problem that I erroneously used variables twice and the code did not work any more. The use of nested functions requires a lot of discipline and due to that is very error-prone, especially when different people work on the same code.
Alex replied on : 8 of 27
Hi Loren, Calling nested functions which access cell arrays seems to take time proportional to the size of the array, while nested functions using numeric arrays do not. I asked about this in the newsgroup (thread "matlab and copy-on-write" -- I thought it was due to a memory copy, but that seems to be incorrect) but have yet to find an explanation or work-around. Any ideas?
Ged replied on : 9 of 27
Hi Loren, > With persistent, you only have one instance This is an obvious advantage for nested functions when multiple copies are useful, but how do the two options compare when both are feasible? In terms of speed, etc? Instead of memoising, I have lots of data which I run a computationally intensive function on, parts of which are invariant to some aspects of the particular data. I can use either persistent variables or nested functions to store partial results for the invariant parts, speeding up the function for further calls on remaining data. Alternatively, I could use other functions to compute the invariant parts and simply pass these in to a different function. I expected the first option (persistence) to be much quicker than the third (passing) since no memory needs to be copied back and forth, but I note on profiling that the persistent line itself seems to be quite a slow part of the code (when repeated over many data). I hadn't considered the second option (nesting) at all; I wonder if it might be quicker than persistence... If so, why? Persistent arrays seem like a concept that should map efficiently to the lower level implementation, while nested functions keeping track of their own workspaces sounds more complex. Many thanks, Ged
Ged replied on : 10 of 27
Hello again, Had a quick look at this myself... It seems that with largish arrays (e.g. 1000x1000), the options are much the same, but that with smaller arrays (e.g. 100x100) persistent is slower than nested, with the persistent declaration line itself being one of the slower lines in the profile. I can't think of any good reasons for this, but it looks to me like nested beats persistent in any case, as well as allowing you to have multiple instances. Hope that's of interest to someone (and that maybe Loren or someone else can suggest reasons why persistent isn't better) Ged
Loren replied on : 11 of 27
Ged- Performance aside (though it is a legitimate question), persistent doesn't easily allow multiple instantiations, for example, 2 gui screens, while nested functions allow this naturally. --Loren
David replied on : 12 of 27
Hi, I was just wondering: the function 'inner' retains the variable values between runs. When does memory allocated by it get freed?
Loren replied on : 13 of 27
David- inner itself doesn't retain the values. The function handle to the nested function does. All related get cleared up when the last reference to the nested function is cleared. However, local variables in inner get cleared right after each time it is run. Only x and y are retained with the nested function handle. --Loren
Yi Cao replied on : 14 of 27
Loren It is an interesting usage of nested function. Two questions: 1. After complete all calls and save f to a mat file, then clear whole workspace. Now, load the mat file and check with whos, the function handle, f only occupies 16 bytes. But check with s=functions(f), and s.workspace{1} all information is still there. I worder where the information is stored? 2. In the inner function, x and y are not pre-allocated. Will this affect the performance when the sizes of x and y becomes large? -Yi
Loren replied on : 15 of 27
Yi- The rest of the information is also stored in the mat-file. What you are seeing is the size of the function handle itself, not everything it encapsulates. As for pre-allocation, that would be nice, and it might help performance, especially if I chunked it up - but it can't solve everything as we don't know how big to make it. Whether or not that slows things down materially depends on the very picky details of how it gets called and how much and how often the array is grown. For the example, it was going to make things muddier and harder to understand, I think, to write the chunked memory version. --Loren
Ljubomir Josifovski replied on : 16 of 27
Does anyone know of a way to define a function that does not compute its arguments at the point of the function call, but rather "on-demand" so to say, as they are needed/referenced inside it? I am thinking about a typical "verbose printing" which is something like if (verbosityLevel(...)) fprintf(...). If this construct is moved away in its own function, at the point of call the (maybe costly) arguments to fprintf are computed, only to be ignored if false verbosityLevel(). The half-solution seems to be using the operators cross-circuiting like verbosityLevel(...) && fprintf(...). Thanks, Ljubomir
Loren replied on : 17 of 27
Ljubomir- Instead of sending in expressions, perhaps rearchitect your function so it takes in function handles and original data. Then you can decide whether to evaluate functions as needed. You could manage this even with anonymous functions perhaps, e.g., something like this input1 = @() a'*a; Here's an experiment I just ran.
a = rand(1000);
tic, b = a'*a; toc
tic, f = @() a'*a, toc
tic, c = f(); toc
Elapsed time is 0.000887 seconds.
Elapsed time is 0.081510 seconds.
>> tic, c = f(); toc
Elapsed time is 0.074612 seconds.
Setting up the function handle takes little time, but the evaluation is similar to the time for the expression itself. --Loren
Ljubomir Josifovski replied on : 18 of 27
Thanks Loren. Do you maybe know of a way of applying this conveniently to varargin-s? For example in
  pv = @(v, varargin) (v<=1) && fprintf(varargin{:});
  slow = @() fprintf('slow called ...\n');
  pv(1, '   when needed\n', slow());
  pv(2, '   even when not needed\n', slow());
would like to avoid the call to slow() in the course of the 2nd call to pv(), pv(2, ...). Thanks, Ljubomir
Loren replied on : 19 of 27
Ljubomir- I can't think of a good way around it right now without explicitly writing a function instead of using an anonymous function. Note that your second pv statement itself doesn't print, though slow executes. Do you have an example that is more numerically oriented? I bet the logic would be easier (for me) to trace then. You might ask on the newsgroup. Someone there may know how to do this off the top of their heads. --Loren
Loren replied on : 20 of 27
Ljubomir- slow() is getting called each timein yours because the arguments to pv get evaluated before making the call to pv. You need to pass in slow, and not slow(). And you need to rewrite pv to execute slow when required. You are always causing slow to be executed right now. --Loren
Loren replied on : 21 of 27
Ljubomir- I've got something to work, but I am sure you will find it too specific. Perhaps you can work with the idea, though.
pv = @(v, varargin) (v<=1) && fprintf(varargin{1},varargin{2}());
slow = @() fprintf('slow called ...\n');
pv(1, '   when needed\n', slow);
pv(2, '   even when not needed\n', slow);
>> pv(1, '   when needed\n', slow);
slow called ...
   when needed
>> pv(2, '   even when not needed\n', slow);
Ljubomir Josifovski replied on : 22 of 27
Many thanks for your help Loren. Yep, will try to extend your idea. Matlab would not have a general "reference of any expression" operator, right? By this I mean something akin to what a function handle is to a function, but aplicable to any expression (not just functions). Thanks again, Ljubomir
Ljubomir Josifovski replied on : 23 of 27
Just to report back: the (half-?)solution (for avoiding arguments evaluation when print is skipped) I settled for is to use logical ops short circuiting directly and turn the print calls into
   is_verbose_level(...) && fprintf(...);

Bit ugly but does it. Thanks, Ljubomir                
David Young replied on : 24 of 27
Although this blog entry has been around for a while, I only read it recently. I thought a version for functions with multiple arguments and multiple results might be useful, and I've put one in the file exchange as part of this contribution: https://www.mathworks.com/matlabcentral/fileexchange/33068-a-multidimensional-map-class . It follows your pattern quite closely, though the implementation ultimately uses containers.Map rather than == and find.
Rafael replied on : 26 of 27
Hi Loren. I have been converting a lot of our instrument control and data acquisition code to use nested functions and "memoize". I like it because it allows us to think of the code as if it was object-oriented, and has made it a lot more re-usable. But I have been worrying lately about memory usage when I have lots of these function handles associated with persistent workspaces. How can I make sure that those workspaces disappear from memory when we don't need them anymore? In your example, if I have f = memoize(@sin) and I use it for a while, can I clear the persistent workspace by just saying "clear f"? Of is there anything else I need to do? Thanks!
Loren Shure replied on : 27 of 27
Rafael- Clearing every reference to a particular handle should remove the memory associated with that instance. As long as there are no other copies of the reference that f is pointing to (i.e., you didn't say g = f...), then clear f should remove the nested workspace associated with f. In other words, as long as there are no outstanding references, the memory should be returned. --Loren