Efficient tree implementation in MATLAB

Question

Tree class in MATLAB

I am implementing a tree data structure in MATLAB. Adding new child nodes to the tree, assigning and updating data values related to the nodes are typical operations that I expect to execute. Each node has the same type of data associated with it. Removing nodes is not necessary for me. So far, I've decided on a class implementation inheriting from the handle class to be able to pass references to nodes around to functions that will modify the tree.

Edit: December 2nd

First of all, thanks for all the suggestions in the comments and answers so far. They have already helped me to improve my tree class.

Someone suggested trying digraph introduced in R2015b. I have yet to explore this, but seeing as it does not work as a reference parameter similarly to a class inheriting from handle, I am a bit sceptical how it will work in my application. It is also at this point not yet clear to me how easy it will be to work with it using custom data for nodes and edges.

Edit: (Dec 3rd) Further information on the main application: MCTS

Initially, I assumed the details of the main application would only be of marginal interest, but since reading the comments and the answer by @FirefoxMetzger, I realise that it has important implications.

I am implementing a type of Monte Carlo tree search algorithm. A search tree is explored and expanded in an iterative manner. Wikipedia offers a nice graphical overview of the process:

In my application I perform a large number of search iterations. On every search iteration, I traverse the current tree starting from the root until a leaf node, then expand the tree by adding new nodes, and repeat. As the method is based on random sampling, at the start of each iteration I do not know which leaf node I will finish at on each iteration. Instead, this is determined jointly by the data of nodes currently in the tree, and the outcomes of random samples. Whatever nodes I visit during a single iteration have their data updated.

Example: I am at node n which has a few children. I need to access data in each of the children and draw a random sample that determines which of the children I move to next in the search. This is repeated until a leaf node is reached. Practically I am doing this by calling a search function on the root that will decide which child to expand next, call search on that node recursively, and so on, finally returning a value once a leaf node is reached. This value is used while returning from the recursive functions to update the data of the nodes visited during the search iteration.

The tree may be quite unbalanced such that some branches are very long chains of nodes, while others terminate quickly after the root level and are not expanded further.

Current implementation

Below is an example of my current implementation, with example of a few of the member functions for adding nodes, querying the depth or number of nodes in the tree, and so on.

classdef stree < handle
    %   A class for a tree object that acts like a reference
    %   parameter.
    %   The tree can be traversed in both directions by using the parent
    %   and children information.
    %   New nodes can be added to the tree. The object will automatically
    %   keep track of the number of nodes in the tree and increment the
    %   storage space as necessary.

    properties (SetAccess = private)
        % Hold the data at each node
        Node = { [] };
        % Index of the parent node. The root of the tree as a parent index
        % equal to 0.
        Parent = 0;
        num_nodes = 0;
        size_increment = 1;
        maxSize = 1;
    end

    methods
        function [obj, root_ID] = stree(data, init_siz)
            % New object with only root content, with specified initial
            % size
            obj.Node = repmat({ data },init_siz,1);
            obj.Parent = zeros(init_siz,1);
            root_ID = 1;
            obj.num_nodes = 1;
            obj.size_increment = init_siz;
            obj.maxSize = numel(obj.Parent);
        end

        function ID = addnode(obj, parent, data)
            % Add child node to specified parent
            if obj.num_nodes < obj.maxSize
                % still have room for data
                idx = obj.num_nodes + 1;
                obj.Node{idx} = data;
                obj.Parent(idx) = parent;
                obj.num_nodes = idx;
            else
                % all preallocated elements are in use, reserve more memory
                obj.Node = [
                    obj.Node
                    repmat({data},obj.size_increment,1)
                    ];

                obj.Parent = [
                    obj.Parent
                    parent
                    zeros(obj.size_increment-1,1)];
                obj.num_nodes = obj.num_nodes + 1;

                obj.maxSize = numel(obj.Parent);

            end
            ID = obj.num_nodes;
        end

        function content = get(obj, ID)
            %% GET  Return the contents of the given node IDs.
            content = [obj.Node{ID}];
        end

        function obj = set(obj, ID, content)
            %% SET  Set the content of given node ID and return the modifed tree.
            obj.Node{ID} = content;
        end

        function IDs = getchildren(obj, ID)
            % GETCHILDREN  Return the list of ID of the children of the given node ID.
            % The list is returned as a line vector.
            IDs = find( obj.Parent(1:obj.num_nodes) == ID );
            IDs = IDs';
        end
        function n = nnodes(obj)
            % NNODES  Return the number of nodes in the tree.
            % Equal to root + those whose parent is not root.
            n = 1 + sum(obj.Parent(1:obj.num_nodes) ~= 0);
            assert( obj.num_nodes == n);
        end

        function flag = isleaf(obj, ID)
            % ISLEAF  Return true if given ID matches a leaf node.
            % A leaf node is a node that has no children.
            flag = ~any( obj.Parent(1:obj.num_nodes) == ID );
        end

        function depth = depth(obj,ID)
            % DEPTH return depth of tree under ID. If ID is not given, use
            % root.
            if nargin == 1
                ID = 0;
            end
            if obj.isleaf(ID)
                depth = 0;
            else
                children = obj.getchildren(ID);
                NC = numel(children);
                d = 0; % Depth from here on out
                for k = 1:NC
                    d = max(d, obj.depth(children(k)));
                end
                depth = 1 + d;
            end
        end
    end
end

However, performance at times is slow, with operations on the tree taking up most of my computation time. What specific ways would there be to make the implementation more efficient? It would even be possible to change the implementation to something else than the handle inheritance type if there are performance gains.

Profiling results with current implementation

As adding new nodes to the tree is the most typical operation (along with updating the data of a node), I did some profiling on that. I ran the profiler on the following benchmarking code with Nd=6, Ns=10.

function T = benchmark(Nd, Ns)
% Tree benchmark. Nd: tree depth, Ns: number of nodes per layer
% Initialize tree
T = stree(rand, 10000);
add_layers(1, Nd);
    function add_layers(node_id, num_layers)
        if num_layers == 0
            return;
        end
        child_id = zeros(Ns,1);
        for s = 1:Ns
            % add child to current node
            child_id(s) = T.addnode(node_id, rand);

            % recursively increase depth under child_id(s)
            add_layers(child_id(s), num_layers-1);
        end
    end
end

Results from the profiler:

R2015b performance

---

It has been discovered that R2015b improves the performance of MATLAB's OOP features. I redid the above benchmark and indeed observed an improvement in performance:

So this is already good news, although further improvements are of course accepted ;)

Reserving memory differently

It was also suggested in the comments to use

obj.Node = [obj.Node; data; cell(obj.size_increment - 1,1)];

to reserve more memory rather than the current approach with repmat. This improved performance slightly. I should note that my benchmark code is for dummy data, and since the actual data is more complicated this is likely to help. Thanks! Profiler results below:

Questions on even further increasing performance

1. Perhaps there is an alternative way to maintain memory for the tree that is more efficient? Sadly, I typically don't know ahead of time how many nodes there will be in the tree.
2. Adding new nodes and modifying the data of existing nodes are the most typical operations I do on the tree. As of now, they actually take up most of the processing time of my main application. Any improvements on these functions would be most welcome.

Just as a final note, I would ideally like to keep the implementation as pure MATLAB. However, options such as MEX or using some of the integrated Java functionalities may be acceptable.

Answer 1

TL:DR You deep copy the entire data stored on each insertation, initialize the parent and Node cell bigger then what you expect to need.

Your data does have a tree structure, however you are not utilising this in your implementation. Instead the implemented code is a computational hungry version of a look up table (actually 2 tables), that stores the data and the relational data for the tree.

The reasons I am saying this are the following:

• To insert you call stree.addnote(parent, data), which will store all data in the tree object stree's fields Node = {} and Parent = []
• you seem to know prior which element in your tree you want to access as the search code is not given (if you use stree.getchild(ID) for it I have some bad news)
• once you processed a node you trace it back using find() which is a list search

By no means does that mean the implementation is clumsy for the data, it may even be the best, depending on what you are doing. However it does explain your memory allocation issues and gives hints on how to resolve them.

---

Keep Data as lookup table

One of the ways to store data is to keep the underlying look up table. I would only do this, if you know the ID of the first element you want to modify without searching for it. This case allows you to make your structure more efficient in two steps.

First initialise your arrays bigger then what you expect you need to store the data. If the look up table's capacity is exceeded, a new one is initialized, which is X fields larger, and a deep-copy of the old data is made. If you need to expand capcity once or twice (during all insertations) it might not be an issue, but in your case a deep copy is made for ever insertation!

Second I would change the internal structure and merge the two tables Node and Parent. The reason for this is that back-propagation in your code takes O(depth_from_root * n), where n is the number of nodes in your table. This is because find() will iterate over the entire table for each parent.

Instead you can implement something similar to

table = cell(n,1) % n bigger then expected value
end_pointer = 1 % simple pointer to the first free value

function insert(data,parent_ID)
    if end_pointer < numel(table)
        content.data = data;
        content.parent = parent_ID;
        table{end_pointer} = content;
        end_pointer = end_pointer + 1;
    else
        % need more space, make sure its enough this time
        table = [table cell(end_pointer,1)];
        insert(data,parent_ID);
    end
end

function content = get_value(ID)
    content = table(ID);
end

This instantly gives you access to the parent's ID without the need to find() it first, saving n iterations each step, so afford becomes O(depth). If you do not know your initial node, then you have to find() that one, which costs O(n).

Note that this structure has no need for is_leaf(), depth(), nnodes() or get_children(). If you still need those I need more insight in what you want to do with your data, as this highly influences a proper structure.

---

Tree Structure

This structure makes sense, if you never know the first node's ID and thus allways have to search for it.

The benefit is that the search for an arbitrary note works with O(depth), so searching is O(depth) instead of O(n) and back propagating is O(depth^2) instead of O(depth + n). Note that depth can be anything from log(n) for a perfectly balanced tree, that may be possible depending on your data, to n for the degenerated tree, that just is a linked list.

However to suggest something proper I would require more insight, as every tree structure kind of has its own nich. From what I can see so far, I'd suggest an unbalanced tree, that is 'sorted' by the simple order given by a nodes wanted parent. This may be further optimized depending on

• is it possible to define a total order on your data
• how do you treat double values (same data appearing twice)
• what scale is your data (thousands, millions, ...)
• is a lookup / search allways paired with back propagation
• how long are the chains of 'parent-child' on your data (or how balanced and deep will the tree be using this simple order)
• is there allways just one parent or is the same element inserted twice with different parents

I'll happily provide example code for above tree, just leave me a comment.

EDIT: In your case an unbalanced tree (that is construted paralell to doing the MCTS) seems to be the best option. The code below assumes that the data is split in state and score and further that a state is unique. If it isn't this will still work, however there is a possible optimisation to increase MCTS preformance.

classdef node < handle
    % A node for a tree in a MCTS
    properties
        state = {}; %some state of the search space that identifies the node
        score = 0;
        childs = cell(50,1);
        num_childs = 0;
    end
    methods
        function obj = node(state)
            % for a new node simulate a score using MC
            obj.score = simulate_from(state); % TODO implement simulation state -> finish
            obj.state = state;
        end
        function value = update(obj)
            % update the this node using MC recursively
            if obj.num_childs == numel(obj.childs)
                % there are to many childs, we have to expand the table
                obj.childs = [obj.childs cell(obj.num_childs,1)];
            end
            if obj.do_exploration() || obj.num_childs == 0
                % explore a potential state
                state_to_explore = obj.explore();

                %check if state has already been visited
                terminate = false;
                idx = 1;
                while idx <= obj.num_childs && ~terminate
                    if obj.childs{idx}.state_equals(state_to_explore)
                        terminate = true;
                    end
                    idx = idx + 1;
                end

                %preform the according action based on search
                if idx > obj.num_childs
                    % state has never been visited
                    % this action terminates the update recursion 
                    % and creates a new leaf
                    obj.num_childs = obj.num_childs + 1;
                    obj.childs{obj.num_childs} = node(state_to_explore);
                    value = obj.childs{obj.num_childs}.calculate_value();
                    obj.update_score(value);
                else
                    % state has been visited at least once
                    value = obj.childs{idx}.update();
                    obj.update_score(value);
                end
            else
                % exploit what we know already
                best_idx = 1;
                for idx = 1:obj.num_childs
                    if obj.childs{idx}.score > obj.childs{best_idx}.score
                        best_idx = idx;
                    end
                end
                value = obj.childs{best_idx}.update();
                obj.update_score(value);
            end
            value = obj.calculate_value();
        end
        function state = explore(obj)
            %select a next state to explore, that may or may not be visited
            %TODO
        end
        function bool = do_exploration(obj)
            % decide if this node should be explored or exploited
            %TODO
        end
        function bool = state_equals(obj, test_state)
            % returns true if the nodes state is equal to test_state
            %TODO
        end
        function update_score(obj, value)
            % updates the score based on some value
            %TODO
        end
        function calculate_value(obj)
            % returns the value of this node to update previous nodes
            %TODO
        end
    end
end

A few comments on the code:

• depending on the setup the obj.calculate_value() might not be needed. E.g. if it is some value that can be computed by evaluating the child's scores alone
• if a state can have multiple parents it makes sense to reuse the note object and cover it in the structure
• as each node knows all its children a subtree can be easily generated using node as root node
• searching the tree (without any update) is a simple recursive greedy search
• depending on the branching factor of your search, it might be worth to visit each possible child once (upon initialization of the node) and later do randsample(obj.childs,1) for exploration as this avoids copying / reallocating of the child array
• the parent property is encoded as the tree is updated recursively, passing value to the parent upon finishing the update for a node
• The only time I reallocate memory is when a single node has more then 50 childs any I only do reallocation for that individual node

This should run a lot faster, as it just worries about whatever part of the tree is chosen and does not touch anything else.

Answer 2

I know that this might sound stupid... but how about keeping the number of free nodes instead of total number of nodes? This would require comparison against a constant (which is zero), which is single property access.

One other voodoo improvement would be moving .maxSize near .num_nodes, and placing both those before the .Node cell. Like this their position in memory won't change relative to the beginning of the object because of the growth of .Node property (the voodoo here being me guessing the internal implementation of objects in MATLAB).

Later Edit When I profiled with the .Node moved at the end of the property list, the bulk of the execution time was consumed by extending the .Node property, as expected (5.45 seconds, compared to 1.25 seconds for the comparison you mentioned).

Answer 3

You can try to allocate a number of elements that is proportional to the number of elements you have actually filled: this is the standard implementation for std::vector in c++

obj.Node = [obj.Node; data; cell(q * obj.num_nodes,1)];

I don't remember exactly but in MSCC q is 1 while it is .75 for GCC.

---

This is a solution using Java. I don't like it very much, but it does its job. I implemented the example you extracted from wikipedia.

import javax.swing.tree.DefaultMutableTreeNode

% Let's create our example tree
top = DefaultMutableTreeNode([11,21])
n1 = DefaultMutableTreeNode([7,10])
top.add(n1)
n2 = DefaultMutableTreeNode([2,4])
n1.add(n2)
n2 = DefaultMutableTreeNode([5,6])
n1.add(n2)
n3 = DefaultMutableTreeNode([2,3])
n2.add(n3)
n3 = DefaultMutableTreeNode([3,3])
n2.add(n3)
n1 = DefaultMutableTreeNode([4,8])
top.add(n1)
n2 = DefaultMutableTreeNode([1,2])
n1.add(n2)
n2 = DefaultMutableTreeNode([2,3])
n1.add(n2)
n2 = DefaultMutableTreeNode([2,3])
n1.add(n2)
n1 = DefaultMutableTreeNode([0,3])
top.add(n1)

% Element to look for, your implementation will be recursive
searching = [0 1 1];
idx = 1;
node(idx) = top;
for item = searching,
    % Java transposes the matrices, remember to transpose back when you are reading
    node(idx).getUserObject()'
    node(idx+1) = node(idx).getChildAt(item);
    idx = idx + 1;
end
node(idx).getUserObject()'

% We made a new test...
newdata = [0, 1]
newnode = DefaultMutableTreeNode(newdata)
% ...so we expand our tree at the last node we searched
node(idx).add(newnode)

% The change has to be propagated (this is where your recursion returns)
for it=length(node):-1:1,
    itnode=node(it);
    val = itnode.getUserObject()'
    newitemdata = val + newdata
    itnode.setUserObject(newitemdata)
end

% Let's see if the new values are correct
searching = [0 1 1 0];
idx = 1;
node(idx) = top;
for item = searching,
    node(idx).getUserObject()'
    node(idx+1) = node(idx).getChildAt(item);
    idx = idx + 1;
end
node(idx).getUserObject()'