Fast manipulation of multidimensional arrays in Matlab
Kevin P. Murphy
[email protected]
11 September 2002
1
Introduction
Probabilistic inference in graphical models with discrete random variables requires performing various oper
ations on multidimensional arrays (discrete potentials). This is true whether we use an exact algorithm like
junction tree [CDLS99, HD96] or an approximate algorithm like loopy belief propagation [AM00, KFL01].
These operations consist of elementwise multiplication/division of two arrays of potentially different sizes,
and summation (marginalization) over a subset of the dimensions. This report discusses efficient ways to
implement these operations in Matlab, with an emphasis on the implementation used in the Bayes Net
Toolbox (BNT).
2
Running example
We will introduce an example to illustrate the problems we want to solve.
Consider two arrays (tables),
Tbig and Tsmall, where Tbig represents a function over the variables
X
1
, X
2
, X
3
, X
4
and Tsmall represents
a function over
X
1
, X
3
. We will say Tbig’s domain is (1
,
2
,
3
,
4), Tsmall’s domain is (1
,
3), and the difference
in their domains is (2
,
4).
Let the size of variable
X
i
(i.e., the number of possible values it can have) be
denoted by
S
i
. Here is a straighforward implementation of elementwise multiplication:
for i1=1:S1
for i2=1:S2
for i3=1:S3
for i4=1:S4
Tbig[i1,i2,i3,i4] = Tbig[i1,i2,i3,i4] * Tsmall[i1,i3];
end
end
end
end
Similarly, here is a straightforward implementation of marginalizing the big array onto the small domain:
for i1=1:S1
for i3=1:S3
sum = 0;
for i2=1:S2
for i4=1:S4
sum = sum + Tbig[i1,i2,i3,i4];
end
end
Tsmall[i1,i3] = sum;
end
end
1
This preview has intentionally blurred sections. Sign up to view the full version.
View Full Document
Of course, these are not general solutions, because we have hardcoded the fact that Tbig is 4dimensional,
Tsmall is 2dimensional, and that we are marginalizing over dimensions 2 and 4. The general solution requires
mapping vectors of indices (or subscripts) to 1 dimensional indices (offsets into a 1D array), and vice versa.
We discuss these auxiliary functions next, before discussing a variety of different solutions based on these
functions.
3
Auxiliary functions
3.1
Converting from multidimensional indices to 1D indices
If a
d
dimensional array is stored in memory such that the leftmost indices toggle fastest (as in Matlab and
Fortran — C follows the opposite convention, toggling the rightmost indices fastest), then we can compute
the 1D index from a vector of indices, subs, as follows:
ndx = 1 + (i11) + (i21)*S1 + (i31)*S1*S2 + ... + (id1)*S1*S2*...*S(d1)
where the vector of indices is
subs
= (
i
1
, . . . , i
d
), and the size of the dimensions are
sz
= (
S
1
, . . . , S
d
).
(We only need to add and subtract 1 if we use 1based indexing, as in Matlab; in C, we would omit this
step.)
This can be written more concisely as
ndx = 1 + sum((subs1) .* [1 cumprod(sz(1:end1))])
If we pass the subscripts separately, we can use the builtin Matlab function
ndx = sub2ind(sz, i1, i2, ...)
BNT offers a similar function, except the subscripts are passed as a vector:
ndx = subv2ind(sz, subs)
The BNT function is vectorized, i.e., if each row of subs contains a set of subscripts. then ndx will be a
vector of 1D indices. This can be computed by simple matrix multiplication. For example, if
This is the end of the preview.
Sign up
to
access the rest of the document.
 Winter '08
 STAFF
 matlab, Multiplication, Indices, tsmall, Tbig

Click to edit the document details