Packed Arrays

An efficient internal storage format in Mathematica 4

Rob Knapp
Numerical Algorithms Development
Wolfram Research, Inc.
rknapp@wolfram.com

While packed arrays may sound quite mysterious and deep, the idea is really quite simple. When it is clearly appropriate to represent a list of a single type machine numbers (integer, real, or complex) internally as an array, Mathematica does so. The basic technology is actually there in Mathematica 3.0, and is visible in the List operations that were added to the compiler. Basically what has been done for Mathematica 4 is to extend this representation to kernel functions where it is appropriate. The internal representation is essentially invisible to you as a user (you have to use functions in the Developer` context to even see it) so you do not have to worry about it. However, understanding the packed array representation can help you write faster, more memory-efficient code.

Representations of data in Mathematica

If you enter a list in Mathematica,

list  =  {1, 2, 3, 4}

its internal representation may be different from the OutputForm. The basic common representation for lists (and all expressions) is given by the FullForm.

FullForm[list]

Each integer you see in the list is itself an expression. Internally in the Mathematica C code, the data stored carries information about the type of expression (in this case, list, or number which is a machine integer).

If you enter the data a different way,

packed  =  Range[4]

it is internally stored in a different way. The basic print forms are all the same. For all intents and purposes it is the same.

packed  ===  list

You can only tell that it is stored in a different way by using a Developer` context function:

Developer`PackedArrayQ  /@  {list,  packed}

Advantages of packed arrays

There are two advantages of packed arrays, memory and speed.

The memory is the most obvious of these.

ByteCount[a = N[Range[1000]]]

If you change the first element to be a different type, the result cannot be stored as a packed array.

ByteCount[b = a; b[[1]] = 1; 
      b]

\(N[%\/%%]\)

The speed advantage typically comes in common operations.

Timing[Do[a + a, {1000}]]

Timing[Do[b + b, {1000}]]

Understanding why the timing difference occurs is the key to knowing when packed arrays will help your code run faster.

In the first operation, a is stored as a packed array, so when Plus sees it, it knows all of the elements are machine real numbers, so that internally, a function which adds arrays of real numbers is used. Except for overflow and underflow checking, which Mathematica does, this is done as fast as you could get by writing a loop to do the addition in a C-program for two arrays of machine numbers.

In the second operation, b is stored as List[…], so what happens first is Plus is threaded over the lists, producing (roughly) the following.

{b[[1]] + b[[1]], 
      b[[2]] + b[[2]], ....}

Plus processes 1000 expressions. For each pair of numbers, Plus uses the information attached to those numbers (i.e., real, integer, machine, bignum, and so forth) to determine what function to use to do the actual adding. This all takes time and accounts for nearly all of the difference here.

There are different ways that packed arrays can speed things up, but this exemplifies the type which typically gives the greatest increase in speed.

When do you get packed arrays?

When trying to write efficient programs, this is a very important question. There are three important situations that commonly produce packed arrays.

1. Listable and structural operations with packed array input (where the 
output can be a packed array.)

Listable operations mean those with the attribute Listable.

Select[Names["System`*"],  
      MemberQ[Attributes[#],  Listable] &]

The majority of which will produce packed arrays given packed array arguments.

Structural operations mean commands like Part, Join, Transpose, and so forth, which act on the structure of an expression.

2. Functions that use an internal array format for machine numbers.

Fourier is a good example of this:

Developer`PackedArrayQ[Fourier[{1. , 2. , 3. , 4. }]]

Whether the input is a packed array or not, if Fourier does its computation with machine numbers, it internally uses arrays already, so it is natural to give the output in terms of a packed array.

3. Functions which can quickly determine from simple input that the output 
could be a packed array.

You have already used one example in this class, Range. In Range[4], it is easy to see before generating any numbers that the result can be represented as a packed array.

For something like the following, it is a little more difficult to determine.

Timing[
      Developer`PackedArrayQ[
        data  =  Table[Random[], {10000}, {2}]]]

Commands like this are handled by auto-compilation. However, since trying the compiler does take some time, there is a (user modifiable) threshold length that must be exceeded before compilation is attempted.

Can packed arrays ever be slower?

Yes. Here is an example.

f[{a_, b_}]  :=  a^2  -  b^2

Timing[Developer`PackedArrayQ[Map[f, data]]]

Timing[udata  =  Developer`FromPackedArray[data];]

Timing[Map[f, udata];]

However, if you inline the function, the auto-compilation can work.

Timing[
      Developer`PackedArrayQ[
        Map[((#[[1]])^2  +  (#[[2]])^2) &,  
          data]]]

It is not always the case that things are so easily patched up. Some commands are incompatible with packed arrays. However, before complaining too loudly about these cases, please take into account the whole time going into your computation. In the above example, you were not taking into account the time it took to generate the random list. Without auto-compilation and packed arrays, it takes quite a bit longer to compute.

Developer`SetSystemOptions["TableCompileLength" -> 
          Infinity];
Timing[
      Developer`PackedArrayQ[
        data  =  Table[Random[], {10000}, {2}]]]

How do I write code that takes advantage of packed arrays?

The most important thing to do is to make sure that you do operations on lists of elements all at once or in groups whenever possible.

The good news is that this programming principle has not really changed from Mathematica 1.0. Efficient code in Mathematica has always used functional and list based constructs.

The one other thing to keep in mind is that type needs to be consistent for there to even be a packed array representation.

Some aids are provided that can help you in putting together your code.

Developer`PackedArrayForm[Range[1000]]

Developer`SetSystemOptions["UnpackMessage" -> 
          True];

Timing[Apply[Plus,  N[Range[1000]]]]

Timing[Tr[N[Range[1000]]]]

Sometimes, you know your code will run more efficiently if you pack the input beforehand.

SeedRandom[1234];
data  =  Table[Random[Complex], {240}];

one[a_]  :=  
        Max[Abs[1.   -  ( Sin[a]^2  -  Cos[a]^2)]];

Timing[one[data]]

One[a_]  :=  
      Module[{a1  =  Developer`ToPackedArray[a]}, one[a1]]

Timing[One[data]]

Probably, the best thing to do is show you an example.

Example: LU Decomposition

Here is the code for actually doing an LU decomposition (factorization) from the package LinearAlgebra`GaussianElimination. Since this functionality exists in the kernel, the real purpose of this package is pedagogical, but the original author seems to have gone too far with the intent to stick to standard C-looking code so that people familiar only with procedural code can understand it. The result is a prime example of how you should not write code in Mathematica. Let’s see why and how you can improve it.

lufactor[aa_]:=Module[{a = aa, pivot, ii, iip, i, ip, j, k, mpiv, m, 
            n = Length[aa], tmp},
pivot = Table[i, {i, n}];
For[ii = 1,ii ≤ n,ii++,  
  
mpiv = Abs[a[[pivot[[ii]], ii]]];
k = ii;
For[i = ii + 1,i ≤ n,i++,tmp = 
                          Abs[a[[pivot[[i]], 
                                ii]]];

                      If[tmp > mpiv, mpiv = tmp;  k = i]];
tmp = pivot[[ii]];
pivot[[ii]] = iip = 
                      pivot[[k]];
pivot[[k]] = tmp;
mpiv = a[[iip, ii]];
  
For[i = ii + 1,i ≤ n,i++,
ip = pivot[[i]];
m = a[[ip, ii]] = 
                            a[[ip, ii]]/mpiv;
For[j = ii + 1,j ≤ n,j++,a[[ip, j]] -= m*a[[iip, j]]]];];
LU[a, pivot]]

Here are two test examples so you can compare timings of different versions.

SeedRandom[1234];
machine  =  
        Table[Random[], {80}, {80}];
big  =  Table[Random[Real, {0, 1}, 20], {40}, {40}];

Timing[luf  =  lufactor[machine];]

Timing[LUDecomposition[machine];]

Timing[lufb  =  lufactor[big];]

Timing[LUDecomposition[big];]

Ouch! You are a factor of 300 slower in the machine case and a factor of 5 slower in the bignum example. Before screaming an exasperated “Mathematica is just too slow,” lets try to improve the code.

The first thing noticed is that this code is using For to handle the loops.  In Mathematica, For is slower than the alternative iterator constructs. Try to avoid using For in innermost loops unless it would make for more complicated code. In this example simply changing from For to Do makes a substantial difference.

lufactor1[aa_]:=Module[{a = aa, pivot, ii, iip, i, ip, j, k, mpiv, m, 
            n = Length[aa], tmp},
pivot = Table[i, {i, n}];
Do[  
  
mpiv = Abs[a[[pivot[[ii]], ii]]];
k = ii;
Do[tmp = 
                          Abs[a[[pivot[[i]], 
                                ii]]];

                      If[tmp > mpiv, mpiv = tmp;  k = i], {i, ii + 1, n}];
tmp = pivot[[ii]];
pivot[[ii]] = iip = pivot[[k]];
pivot[[k]] = tmp;
mpiv = 
                    a[[iip, ii]];
  
Do[ip = pivot[[i]];
m = a[[ip, ii]] = 
                            a[[ip, ii]]/mpiv;
Do[a[[ip, j]] -= m*a[[iip, j]],{j, ii + 1, n}],
{i, ii + 1, n}],
{ii, n}];
LU[a, pivot]]

Timing[lufactor1[machine]  ===  luf]

Timing[lufactor1[big]  ===  lufb]

The next change will help make the code more readable. At this point, it won’t seem to make much of a speed difference, but it does help. The idea is to replace the lines of code finding the correct pivot to use with a function.

SetAttributes[FindPivot,  HoldFirst];

FindPivot[pivots_,  a_,  
          c_]  :=  
Module[{k = 
              First[MaxPosition[
                    Abs[a[[Take[pivots, {c, -1}], c]]]]] + c - 1}, 
          pivots[[{c, k}]] = pivots[[{k, c}]]; 
          pivots[[c]]];

This function finds which of the possible pivots (in column c, rows c and above) is largest in absolute value, and swaps with the one in row c.

This, in turn, uses the simple function

MaxPosition[x_]  :=  First[Position[x,  Max[x]]];

and puts it into the main function.

lufactor2[aa_]:=Module[{a = aa, pivot, ii, iip, i, ip, j, k, mpiv, m, 
            n = Length[aa], tmp},
pivot = Table[i, {i, n}];
Do[  
iip  =  FindPivot[pivot,  a,  ii];
mpiv = a[[iip, ii]];
  
Do[ip = pivot[[i]];
m = a[[ip, ii]] = 
                              a[[ip, ii]]/mpiv;
Do[
                            a[[ip, j]] -= m*a[[iip, j]], {j, ii + 1, 
                              n}],
{i, ii + 1, n}],
{ii, n}];
LU[a, pivot]]

Timing[lufactor2[machine]  ===  luf]

Timing[lufactor2[big]  ===  lufb]

You are still going in the right direction, but there is no big savings. The place to look for big savings is the innermost loop. The next change replaces the innermost loop by giving Mathematica commands that operate on all the elements of a list:

lufactor3[aa_]:=Module[{a = aa, pivot, ii, iip, i, ip, j, mpiv, m, 
            n = Length[aa]},
pivot = Table[i, {i, n}];
Do[  
iip  =  FindPivot[pivot,  a, ii];
mpiv = a[[iip, ii]];
  
rc  =  
                    Table[j, {j, ii + 1, n}];
Do[ip = pivot[[i]];
m = -(a[[ip, ii]] /= mpiv);
a[[ip, rc]] += m*a[[iip, rc]],
{i, ii + 1, n}],
{ii, n}];
LU[a, pivot]]

Timing[lufactor3[machine]  ===  luf]

luf3  =  lufactor3[machine];
{Max[Abs[First[luf3]  -  First[luf]]], 
      luf3[[2]]  ===  luf[[2]]}

Timing[lufactor3[big]  ===  lufb]

So replacing the innermost loop was spectacularly successful. There is still a set of nested loops that can be worked on though.

lufactor4[aa_]:=Module[{a = aa, pivot, ii, iip, i, ip, j, k, mpiv, m, 
            n = Length[aa], row},
pivot = Table[i, {i, n}];
Do[  
iip  =  FindPivot[pivot,  a,  ii];
mpiv = a[[iip, ii]];
  
rc  =  
                    Table[j, {j, ii + 1, n}];
prc  =  pivot[[rc]];
m  =  -(a[[prc, ii]]  /=  
                            mpiv);
row  =  a[[iip, rc]];
a[[prc, rc]]  +=  
                      Map[(row #) &,  m],
{ii, n}];
LU[a, pivot]]

Timing[lufactor4[machine]  ===  luf3]

Timing[lufactor4[big]  ===  lufb]