# Vectorization tricks for cell arrays in MATLAB

## Motivations

As is known, MATLAB is notoriously slow in executing `for`

loops. Looping through a large array is usually a nightmare, even more so if we add `if/then`

within, and/or for sparse matrices. Nevertheless, MATLAB is highly optimized in vectorized array and matrix operations using the LAPACK/BLAS backend, and as an interpreted/scripting language, MATLAB
operates in 4k by 4k matrices hundreds even thousands of times faster than direct implementations in compiled languages like C/C++ and Java.

Recently I am learning to code “virtual element method” for Long Chen’s $i$FEM. As it requires a polytopal data structure, the biggest difference with traditional finite element is:

The number of indices stored in each element is of variable lengths.

This poses some great challenges when we want to write a general purpose code to retrieve sparse matrices based on the topological structures.
In `numpy`

, the solution is usually adding paddings. Yet MATLAB has a more neat solution to this, which is the cell array:

```
a = {[2,3], [], [2, 3, 1]}
```

then accessing the element of `a`

can either be done through `{}`

or `()`

indexing, with the latter being the preferred way of manipulating it as `a(i)`

returns the $i$-th entry as a cell, while `a{i}`

will yield a matrix or vector whenever applicable.

### First example

Consider the following scripts to generate something like a variable length dice rolling simulation:

```
N = 100000;
result = cell(N,1);
nRoll = randi(10, [N,1]); % num of rolls each trial is between 1 to 10
for i = 1:N
result{i} = randi(6, [1, nRoll(i)]);
end
```

the first 5 elements of `result`

may look like

```
result
ans =
[6,6,5,5,1,4]
[2,1,3]
[1,4,4,6,3,4]
[4,6,1,1,1,6,6,2,6]
[5,6,2,4,4,4,6,4,4,6]
```

Now we want to retrieve the number of rolling in each sequence (i.e., recovering `nRoll`

from `result`

), the first way is through a `for`

loop:

```
avg1 = zeros(N,1);
for i = 1:N
avg1(i) = length(result{i});
end
```

On my old Retina Macbook Pro 2013, 10 runs yield an average of `0.03`

seconds to execute. However, if we apply a built-in function through `cellfun`

:

```
avg2 = cellfun('length', result);
```

the average time is only `0.004`

seconds for 10 runs. While the speed-up may not sound so significant, let us move to the next example.

### Element-wise matrix vector multiplication

Now suppose we have a bunch of matrices $A$ stored in a 3D array, we want to compute element by element matrix vector multiplication $b = A\,\mathbf{x}$. An example setup can be:

```
N = 100000;
n = 10;
A = rand(N,n,n);
x = rand(N,n);
```

We want to obtain another array `b`

such that `b(i,:) = A(i,:,:)*x(i,:)`

, of course the simplest is using a `for`

loop to execute this literally:

```
b = zeros(N,n);
for i = 1:N
b(i,:) = A(i,:,:)*x(i,:);
end
```

It takes about `1.5`

seconds. Now we vectorize this routine using `bsxfun`

:

```
b = sum(bsxfun(@times, A, x), 2);
```

This takes about `0.03`

seconds, which is about 50 times speed-up.

#### Remark on tensor prod

Note that this is NOT a tensor product between these two arrays, as

\([A_1, A_2] \otimes [\mathbf{x}_1, \mathbf{x}_2] = \begin{pmatrix} A_1 \mathbf{x}_1 & A_1 \mathbf{x}_2 \\ A_2 \mathbf{x}_1 & A_2 \mathbf{x}_2 \end{pmatrix}\) computes more off-diagonal entries than the element-wise matrix vector product. The method here can be viewed as a vectorization for retrieving the diagonal blocks of a general tensor product between 3D arrays.

## Use accumarray as np.append for cells

MATLAB has a built-in function dedicated to vectorization which is `accumarray`

which adds up `VALS`

in `SUBS`

by default in the expression of:

```
A = accumarray(subs,val,sz)
```

For example, if you have say

```
subs = [1 2 3 8 1 3 5 6 1 5 5]';
ix = accumarray(subs,1);
```

and `ix`

will be `[3 1 2 0 3 1 0 1]`

, by letting `VALS`

be 1, we are just recording the number of times each number appearing in the `subs`

array.

Now say we do not want to add up the times, instead we just want to append 1’s each time an index appears in the array, i.e., we want to get a cell array like the following:

```
[1,1,1]
1
[1,1]
[]
[1,1,1]
1
[]
1
```

`1`

is appended three times as `1`

appears three times in `subs`

, and empty cells represent that number has not appeared once in the `subs`

. There is a single line of code can do this:

```
ix = accumarray(subs,1,[],@(x) {x});
```

which takes advantage of the fourth variable in `accumarray`

can be a function handle.

## mat2cell is cell array’s reshape

`reshape`

function plays somewhat a central role in MATLAB vectorization for matrices and vectors. For example, to convert back and forth between `:`

(linear indexing) and the `[row,col]`

indices, sometimes there are too many pointers involved the MATLAB built-in `sub2ind`

does not meed our needs; or sometimes we would like to do some row-wise reshape instead of the default column-based linear indexing, `reshape`

are handy.

In the context of cell array, `reshape`

cannot be used. There is a function `num2cell`

but can only output uniform result (no variable length arrays within each cell). `mat2cell`

is here to rescue.

Suppose we have the following scenerio, we have a cell array recording how the connection of a node with neighbors in a graph:

```
C = {
[2, 3, 4, 6], ...
[1, 5, 4], ...
[2, 4], ...
[2], ...
[4, 6], ...
[1, 2, 3, 4, 5]
}
```

Say we would like to check which nodes are connected with node 2. If `C`

has consistent sizes, a simple `C == 2`

yields a logical matrix that gives you in this information; if an indicator is more preferred, `sum(C==2, 2)`

or `any(C==2,2)`

will do the trick. With more pointers involved, you may want to use `reshape(C(:) == 2, size(C))`

as well. However, none of these simple elegant code works for cell array.

Nevertheless, there does exist a simple elegant code to achieve something like `reshape(C(:) == 2, size(C))`

for cell array. First we can **flatten** the cell array by `:`

syntax:

```
C_flat = [C{:}];
```

or

```
C_flat = horzcat(C{:});
```

Then we obtain the logical array we want for this flattened array:

```
isConnected = (C_flat == 2);
```

Now we would like to reshape `isConnected`

to have the shape of `C`

, the first row has 4 indicators, the second row has 3, etc. We can do this by a trick used earlier:

```
sizeC = cellfun(@length, C);
isConnected = mat2cell(isConnected, sizeC);
```

Now `isConnected`

has the same shape as `C`

. To get an indicator array for each cell, `cellfun`

can be used again:

```
cellfun(@any, isConnected)
```

### A less-convoluted example

**(Updated July 27, 2021)**
In the comment, Jan Valdman asked how we can convert the following code to a one-liner: for the `C`

defined above

```
v=[(1:8)' (2:9)']
v =
1 2
2 3
3 4
4 5
5 6
6 7
7 8
8 9
```

Now say we want to vectorize the following routine:

```
vC={v(C{1},:), v(C{2},:), v(C{3},:),v(C{4},:), v(C{5},:), v(C{6},:)}
```

I think this is actually a more illustrative example than the one in my original post. So let us give a try. Using the `mat2cell`

trick above again, the full code should be:

```
sizeC = cellfun(@length, C); % get the length of each cell
C_flat = [C{:}]; % flatten the cell array
vC = mat2cell(v(C_flat,:), sizeC); % reshape the flattened array to the shape of C
```

### Final remark on cellfun

Notice that if we use some customized function handle on `cellfun`

, for example,

```
cellfun(@(x) yourOwnFunction(x), yourCellArray,'UniformOutput',false)
```

we will lose some of the charm in speedup as MATLAB reverts back to essentially executing `for`

loops.

## Linear indexing for cell array

Continuing from the previous section, we may encounter some scenarios such as accessing $t$-th cell’s $i$-th element in a cell array, with $t$ and $i$ both changing. In the context of matrix vectorization, using linear indexing with `sub2ind`

and `ind2sub`

are extremely useful, however, there is no such thing for cell array. We can construct a linear indexing for cell array again using some of the tricks above.

Say we have:

```
T = {
[1,10,11,2],
[2,11,12,7,3],
[11,28,29,21,18,12],
[10,27,28,11],
[27,34,35,28],
[28,35,36,37,29],
[35,40,41,36],
[36,41,42,37],
[41,45,46,42]
}
```

and we want to build an array `A`

such that `A(j) = T{t(j)}(i(j))`

, a `for`

loop can be used as follows:

```
t = [1 1 3 4 9];
i = [2 3 3 1 4];
for j = 1:length(t)
A(j) = T{t(j)}(i(j));
end
```

and

```
A
ans =
[10 11 29 10 42]
```

BTW here we assume `t`

is sorted in ascended order, and `i`

is well-posed in that it does not exceed the size of the corresponding cell.

The vectorization of the routine above can benefit from the linear indexing of the cell array `T`

: if we can convert `T`

to its linear indexed counterpart `T_flat`

, and retrieving the linear `ind`

for `(t,i)`

, then constructing `A`

is as simple as `A = T_linear(ind)`

.

The routine can be as follows: first we flatten `T`

as before:

```
T_flat = [T{:}]';
```

To construct the linear indexing we need the length array again and using `cumsum`

:

```
sZ = cellfun('length', T);
ix = [0; cumsum(sZ(1:end-1))]';
```

The array `ix`

has same dimension with `T`

and acts in place of the `mod`

operator used in `ind2sub`

:

```
ix
ans =
[0 4 9 15 19 23 28 32 36]
```

Now if we want to retrieve the say the 4th cell’s 2nd entry, `ix(4)+2`

would yield the correct linear index for `T_flat`

: now a single line of code yield the same result as before.

```
T_flat(ix(t)+i)'
ans =
[10 11 29 10 42]
```

## Applications

Some routines above are used in functions on retrieving topological relation and obtaining topological indicators related to geometries for edges in a polygonal mesh here. The journal article references is (Cao, 2021).

## Comments