Arrays Everywhere

Rob Knapp, Wolfram Research

Abstract:  Working with arrays is fundamental to efficient use of the Wolfram Language
and appear in almost any computation.  The basic representation for an array
is a nested list, but there are other representations such as PackedArray,
SparseArray, and structured arrays that are more effective for some applications.
This presentation will show how to use the various representations of arrays
combined with compilation to make your computations more efficient.

Arrays

... come in all shapes and sizes

In this presentation, I will focus on a fairly specific definition of array:

An regular arrangement of data that has a depth d  with dimensions of length d and such that each element of the data can be accessed via part indexes of length d

This includes vectors, matrices, and more.

In the Wolfram language, for an array 𝒜 we have

ArrayDepth[𝒜] gives d

Dimensions[𝒜] gives the dimensions

Part[𝒜,(Subscript[p, 1]), … ,(Subscript[p, d])] gives an element when 1<=(Subscript[p, k])<=(Subscript[dims, k]), for 1<=k<=d

Testing Arrays

Two expressions:

expr1 = {{1,2,3},{4,5}};
expr2 = {{1,2,3},{4,5,6}};

Test if they are arrays

ArrayQ[expr1]

ArrayQ[expr2]

Dimensions[expr1]

Dimensions[expr2]

each element is accessible via part with 2 indexes

expr2[[2,1]]

You can test the elements at the same time

ArrayQ[expr2, _, IntegerQ]

You can also test for the depth

ArrayQ[expr2, 2, IntegerQ]

Note VectorQ and MatrixQ are just ArrayQ with specific depths

MatrixQ[expr2]

The elements by default must not themselves be arrays, but we can change this

ArrayQ[expr1, 1, True&]

ArrayQ[expr1, 1,MatchQ[#,{_Integer?Positive..}]&]

Constructing Arrays

Here is the most obvious way:

Array[f,{2,3}]

A convenient way to generate one with symbolic elements

syms = Array[Indexed[a,{##}]&,{2,3,4}]

A constant array is an extension of a scalar ...

ConstantArray[0.,{3,1}]

... or arrays

ConstantArray[syms,{2,2}]

Dimensions[%]

Arrays of random values are convienient for testing operations

rarr = RandomInteger[{-9,9},{3,2,1}]

... Constructing Arrays

A very common constructor is Table

Table[Indexed[a,{i,j}],{i,2},{j,3}]

This is the same as

Table[Table[Indexed[a,{i,j}],{j,3}],{i,2}]

chebNodes[n_] := Table[Cos[π (i-1/2)/n],{i,n,0,-1}]

n5 = chebNodes[5]

You can check with Table again

Table[FullSimplify[ChebyshevT[5,x]],{x,n5}]

Suppose you wanted to make a grid of these points;
this is easily done with an outer product.

out = Outer[List, n5,n5];
ListPlot[out]

Dimensions[out]

Arrays Lurking in the Shadows

Arrays really are (almost) everywhere in the Wolfram Language, even if you don’t always see them.

image=;

ArrayQ[ImageData[image]]

Dimensions[ImageData[image]]

Spikey hardly looks like an array

spikey=PolyhedronData["RhombicHexecontahedron"]

ArrayQ[coordinates=spikey[[1,1]]]

Max[N[coordinates]]

cells = spikey[[1,2]];
ArrayQ[cells[[1]]]

Head[cells]

mr = TriangulateMesh[BoundaryDiscretizeGraphics[spikey]];

coordinates=MeshCoordinates[mr];
cells = MeshCells[mr,3];

You can operate on the array data (within limits!) to effect a transformation

morphed = MeshRegion[Map[#*Norm[#-3]&,coordinates],cells]

Another example:

ifun = NDSolveValue[{x'[t]==x[t],x[0]==1},x,{t,0,1}]

Map[ArrayQ, ifun]

InputForm[ifun]

Representations of Arrays

The default fundamental representation of arrays in the Wolfram Language is a List of Lists ...

FullForm[{{1,2,3},{4,5,6}}]

FullForm[{{{1,2},{3,4},{5,6}},{{7,8},{9,10},{11,12}}}]

There are other implicit and explicit representations:

PackedArray (implicit)

SparseArray

Structured representations

QuantityArray

Structured Matrices (e.g. PermutationMatrix, BlockDiagonalMatrix)

NumericArray

PackedArray representations

This is an implicit “under the hood” way of representing arrays that have consitently typed machine sized numbers that works automatically

pa = Range[47]

Developer`PackedArrayQ[pa]

One big advantage of PackedArrays is the memory required

ByteCount[pa]

ua = Table[i,{i,47}];
ByteCount[ua]

Because the representation is implicit, even though the storage is different the objects are the same

SameQ[pa, ua]

You may be asking why Table is not producing a PackedArray.
Here is a cyptic answer

ListPlot[Table[{n,ByteCount[Table[i,{i,n}]]},{n,500}]]

"TableCompileLength" /.("CompileOptions" /.SystemOptions["CompileOptions"])

PackedArray conversions

You can convert something with consistently typed elements to a PackedArray

Developer`PackedArrayQ[Developer`ToPackedArray[ua]]

The conversion is not done when the types are not consistent

mixed = ReplacePart[ua,1.,1];
Developer`PackedArrayQ[Developer`ToPackedArray[mixed]]

You can specify coercion to machine reals

Developer`PackedArrayQ[rmixed = Developer`ToPackedArray[mixed,Real]]

Conversion from a PackedArray always works:

Developer`PackedArrayQ[Developer`FromPackedArray[rmixed]]

PackedArray operations

When inputs to an operation are PackedArrays and the return type is consistent, it is often the case that the operation will return a PackedArray.  When this happens, often the operation can be done quite a bit faster too

pa = Range[1663];
ua = Developer`FromPackedArray[pa];

RepeatedTiming[Developer`PackedArrayQ[Outer[List,pa,pa]]]

RepeatedTiming[Developer`PackedArrayQ[Outer[List,ua,ua]]]

There are messages, off by default, that you can turn off to see if an operation has to unpack to work

On["Packing"]

pa = Range[10];
ua = Developer`FromPackedArray[pa];

Map[f,pa]

These messages can be useful in figuring out what operations may not be playing nicely with PackedArrays, but will be irritating as I go further, so back off they go:

Off["Packing"]

SparseArray representations

Here is a sparse array

n = 47;
sparse = SparseArray[{{i_,i_}:>i,{n,1}->1,{1,n}->-1},{n,n}]

It can be converted to the normal representation with Normal

dense = Normal[sparse]

Both represent the same data and the basic structural operations reflect this:

Map[Function[{op}, SameQ[op[sparse], op[dense]]],{ArrayDepth, Dimensions, Part[#,7,7]&}]

The objects are not the same since they are different ways of representing the data

SameQ[dense, sparse]

The sparse representation saves memory:

{ByteCount[sparse], ByteCount[dense]}

Sparse arrays are useful for many tasks.
One particular case is numerical linear algebra

b = RandomReal[1,n];
sparse = N[sparse];
dense = N[dense];

{ts,xs} = RepeatedTiming[LinearSolve[sparse, b]];
{td,xd} = RepeatedTiming[LinearSolve[dense, b]];
{ts, td,td/ts,Max[Abs[xs-xd]]}

With bigger size, the timing difference gets larger

n = 1663;
dense = DiagonalMatrix[Range[n]];
dense[[n,1]]=1;
dense[[1,n]]=-1;
dense = N[dense];
sparse = SparseArray[dense];
b = RandomReal[1,n];

{ts,xs} = RepeatedTiming[LinearSolve[sparse, b]];
{td,xd} = RepeatedTiming[LinearSolve[dense, b]];
{ts, td,td/ts,Max[Abs[xs-xd]]}

(Sparse) Representation contract

Even though the sparse representation is different from the list of lists, it is treated as structurally the same and generally you should expect the following identity to hold for most operations op:

Normal[op[sparse]] == op[Normal[sparse]]

When there is not a sparse algorithm to do the operation, then the sparse representation will be converted with normal and the dense algorithm will be used.

Examples:

Normal[LinearSolve[sparse,b]] == LinearSolve[Normal[sparse],b]

Normal[Total[sparse]]==Total[Normal[sparse]]

Structured representations

There are additional representations based on particular structures.
I will give two examples of these, QuantityArray and PermutationMatrix.

These generally satisfy the same properties that I mentioned for the sparse representation:

Structural operations like ArrayQ, ArrayDepth, Dimensions give the same results

Generally, Normal[op[representation]] == op[Normal[representation]]

When there is a special algorithm for an operation, it is used, otherwise the representation is converted to another representation for which there is an algorithm.

QuantityArray

Construct a 2-column matrix that has meters and seconds for the units of the columns:

RepeatedTiming[qa = QuantityArray[RandomReal[1,{100,2}],{"Meters", "Seconds"}]]

Here is the normal representation

Short[nqa=Normal[qa]]

Huge memory savings

{ByteCount[qa],ByteCount[nqa]}

When an special algorithm is defined, there is often a huge difference in execution times

{tqa,sqa}=RepeatedTiming[Total[qa]];
{tn,sn} = RepeatedTiming[Total[nqa]];
Column@{tqa,tn,tn/tqa,sqa, sn,Normal[sqa]==sn}

PermutationMatrix

This is a simple example of a set of structured matrices that are currently under development.

Create a random permutation:

perm = RandomSample[1;;1663,1663];

Construct a permutation matrix from this.

pm = PermutationMatrix[perm]

m = Normal[pm];
{ByteCount[pm], ByteCount[m]}

Compare the timings of matrix products to permute rows and columns of a random matrix

a = RandomInteger[{-9,9},{1663,1663}];
{tp,dp} = RepeatedTiming[pm.a.Transpose[pm]];
{tm,dn} = RepeatedTiming[m.a.Transpose[m]];
{tp,tm,tm/tp,Normal[dp]===dn}

Note that in this case the algorithm is quite simple

RepeatedTiming[partp=a[[perm,perm]];]

partp===dn

NumericArray

NumericArray is an array that has elements that are machine sized numbers that can be different from the defaults that the Wolfram language uses.  These are useful for situations where memory savings can be very important or in interfaces to code that supports single precision floating point operations (e.g. GPUs, machine learning)

Because conversion to normal arrays loses this type information, they are not considered as direct representations.

Basic structural operations and Normal work this this as for the representations, but many other operations do not.  The best way to do operations on these is using the compiler with the appropriate types.

na = NumericArray[RandomInteger[9,1234],"Integer8"]

Developer`PackedArrayQ[pa = Normal[na]]

Since packed arrays use 64 bit integers, there is a significant memory savings

{ByteCount[na], ByteCount[pa]}

Map[#^2&,na]

Array Arithmetic (and elementary functions)

Basic arithmetic in the Wolfram language is done elementwise.

ax = Array[Indexed[x,{##}]&,{3,4}];
ay =Array[Indexed[y,{##}]&,{3,4}];
v3 = Array[Indexed[v,{##}]&,{3}];
v4 = Array[Indexed[v,{##}]&,{4}];

MatrixForm[ax+ay]

When the top dimensions are consistent, then arrays of different depths can be combined

MatrixForm[ax+v3]

This does not work

ax+v4

If you want to add the vector to each of the columns, then you can use Threaded:

MatrixForm[ax+Threaded[v4]]

Listable

The elementwise nature of arithmetic works through a very general mechanism triggered by the Listable attribute.

SetAttributes[f,Listable];
MatrixForm[f[ax,ay]]

MatrixForm[f[ax,Threaded[v4]]]

There a a lot of system commands that are

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

Listable speed

It is almost always the case that using the elementwise arithmetic directly is faster (and usually more convenient) than iterating over the elements.

I’ll use the Chebyshev nodes from before as a simple example
(I’ve changed it here to use machine numbers)

chebNodes[n_] := Table[Cos[π (i-.5)/n],{i,N[n],0.,-1.}]

chebNodesListable[n_] := Cos[π (Range[n,0,-1]-.5)/n];

TableForm[Table[
n = 10^k;
{tt,nt} = RepeatedTiming[chebNodes[n]];
{tl,nl} = RepeatedTiming[chebNodesListable[n]];
{tt,tl,tt/tl,Max[Abs[nt-nl]]},
{k,5}],TableHeadings->{10^Range[5],{"chebNodes", "chebNodesListable", "ratio", "difference"}}]

Here is another example -- shift coordinates in 2D so that their mean is zero.

coordinates = RandomReal[1,{47,2}];
mean = Mean[coordinates]

{tl,cl} = RepeatedTiming[coordinates-Threaded[mean]];
{tm,cm} = RepeatedTiming[Map[(#-mean)&, coordinates]];
{tl,tm,tm/tl,Mean[cl],Max[Abs[cl-cm]]}

Listable operations with array representations

The listable operations work with the array representations I have shown too, converting to other representations if needed.

n = 7;
s = SparseArray[{{i_,i_}:>i,{n,1}->1,{1,n}->1},{n,n}];
p = PermutationMatrix[RotateRight[Range[n]]];

MatrixForm[s + p]

When both operands are sparse, the result is sparse

s+SparseArray[p]

What if your operation is not listable

If your function does not consist of combining listable operations, often compilation can make things fast

SetAttributes[tent1, Listable];
tent1[x_] := Which[
x<=-1,0.,
x<=0,1+x,
x<=1,1-x,
True,0.];
Plot[tent1[x],{x,-2,2}]

This is slow

data = RandomReal[1,10^6];
RepeatedTiming[tent1[data];]

tentDeclaration = FunctionDeclaration[tentFun, Typed[ForAllType[{"a"},{Element["a","FloatingPointReals"]},{"a"}->"a"]]@
Function[{x},
Module[{zero = Typed[0,TypeOf[x]], one = Typed[1,TypeOf[x]]},
Which[
x<=-one,zero,
x<=zero,one+x,
x<=one,one-x,
True,zero]
]
]
];

tentR64 = FunctionCompile[tentDeclaration, Function[{Typed[x,"Real64"]},tentFun[x]]]

SetAttributes[tent2, Listable];
With[{cf = tentR64}, tent2[x_] := cf[x]]

This is a little bit faster

RepeatedTiming[tent2[data];]

Now compile to work on a PackedArray

tentP64 = FunctionCompile[tentDeclaration, Function[{Typed[x,"PackedArray"["Real64",1]]},Table[tentFun[xi],{xi,x}]]]

This is a lot faster

RepeatedTiming[tentP64[data];]

The advantage of writing the code polymorphically is that you can make it work for NumericArray too

tentN32 = FunctionCompile[tentDeclaration, Function[{Typed[x,"NumericArray"["Real32",1]]},Module[{res = x},
Do[res[[i]]=tentFun[x[[i]]],{i,Length[x]}];
res
]]]

Slightly faster because 32 bit floats are faster than double (64 bit) floats

nadata = NumericArray[data, "Real32"];
RepeatedTiming[tentN32[nadata];]

You do have to compile a function for each type (a lot of the speed comes from not having to change for various types).  There are also ways to set up DownValues that do some type dispatching to make it more convenient.
Here is a (not totally optimal) definition that handles machine real arrays of any depth.

With[{sf =  tentR64,af = tentP64}, 
tent3[x_?NumericQ] := sf[x];
tent3[x_?ArrayQ]:= Module[{d = ArrayDepth[x],dims,y=x},
If[d > 1,
dims=Dimensions[x];
y = Flatten[x]
];
y = af[y];
If[d > 1,
ArrayReshape[y,dims],
y
]
]
];

tent3[RandomReal[1,{2,3,4}]]

This is a bit slower -- it is worth compiling separately for commonly used depths

data3 = ArrayReshape[data,{100, 100, 100}];
RepeatedTiming[tent3[data3];]

Array Valued Variables

The default assumption for variables in the Wolfram language is that they are scalar complex valued.
Array valued variables are a powerful concept in modelling, making descriptions both more efficient and more readable.

NDSolve and many optimization functions can currently handle array valued variables
and we will be expanding the capabilities in the future.

A nice example is finding the distance between two convex regions.

ℛ1 = ConvexHullMesh[RandomReal[1,{10,3}]];
ℛ2 = Ball[{-1,-1,-1}];

{minval, minxy} = NMinimize[Norm[x-y],{x∈ℛ1 ,y∈ℛ2}]

Show[ℛ1,Graphics3D[{ℛ2,Red, Thick, Line[{x,y}] /. minxy}]]

Another example is minimizing the norm of a symmetric matrix given some fixed elements and trace

NMinimize[{Norm[x],{(IndexedDefault[x, 4, 1])==1,(IndexedDefault[x, 1, 4])==1, Tr[x]==7}},Element[x,Matrices[4,Reals, Symmetric]]]

⚠:  Be very careful when combining array valued variables with array constants

NMinimize[{Norm[x-Reverse[IdentityMatrix[4]]],{(IndexedDefault[x, 4, 1])==1,(IndexedDefault[x, 1, 4])==1, Tr[x]==7}},Element[x,Matrices[4,Reals, Symmetric]]]

The problem is the following

x - Reverse[IdentityMatrix[4]]

When x is a matrix, the array depth of this is 4, and we wanted a matrix norm

You can use a constant parameter

NMinimize[{Norm[x-rid],{(IndexedDefault[x, 4, 1])==1,(IndexedDefault[x, 1, 4])==1, Tr[x]==7,rid ==Reverse[IdentityMatrix[4]]}},Element[x,Matrices[4,Reals, Symmetric]]]

You can use Inactive

NMinimize[Norm[Inactive[Plus][x,-rid]],{(IndexedDefault[x, 4, 1])==1,(IndexedDefault[x, 1, 4])==1, Tr[x]==7,rid ==Reverse[IdentityMatrix[4]]},Element[x,Matrices[4,Reals, Symmetric]]]

or, someday in the future, you might be able to do

x= SymbolicMatrix["x",4,Reals,Symmetric];
NMinimize[Norm[x-Reverse[IdentityMatrix[4]]],{(IndexedDefault[x, 4, 1])==1,(IndexedDefault[x, 1, 4])==1, Tr[x]==7},x]

When x is set to (effectively declared as) a symbolic matrix, the listable arithmetic does not 
evaluate fully as it does when it considers x to be a scalar.

Summary

Arrays are just about everywhere

Arrays can be represented in different ways

List of lists

PackedArray (implicit)

SparseArray

Structured representations

NumericArray

Using appropriate representations for all or part of computations can give you some enormous speedups.

(Listable) Arithmetic on arrays is elementwise and very efficient

Use array operations when you can

If your code works with machine numbers of consistent type check if it is running with PackedArrays

The compiler can help speed up code that is not already listable

Array valued variables are expressive and efficient