APL inspired arrays for Haskell

​

Orthotope

Lennart Augustsson

(Epic Games)

What’s in a name?

From Wikipedia:

​

In geometry, an orthotope^[2] (also called a hyperrectangle or a box) is the generalization of a rectangle to higher dimensions. It is formally defined as the Cartesian product of orthogonal intervals. 

Since all the good names like array and tensor were taken I picked orthotope.

APL history

• APL is an array language
• Designed in the 1960s by Ken Iverson
• Much ahead of its time; very functional
• Most known for being cryptic and using special symbols
• (~R∊R∘.×R)/R←1↓⍳R�(this computes prime numbers up to R)

Orthotope

• A Haskell array library
• Provides mostly structural (i.e., no computation) operations on arrays
• Many operations in constant time
• Here presented by examples
• And a thanks to Google for making it open source

Array notation

Array notation

Rank 0 arrays are allowed, often called scalars.

Array type

• The orthotope library contains variations of how strongly typed the functions are.
• The examples will use the version that has the shape as part of the type.
• This requires various Haskell extensions
• Some examples

Array [2,3] Int -- 2x3 array of Int

Array [6] Float -- 6 element array (i.e. vector) of Int

Array [4,4,4,4] Bool — 4x4x4x4 array of Bool

Some examples: fromList, shape, rank, size, iota

> a = fromList [1..6] :: Array [2,3] Int

> pp a -- pp pretty prints

1 2 3

4 5 6

> :t fromList

fromList :: Shape sh => [a] -> Array sh a

> a = fromList @[2,3] [1..6]

> pp a

1 2 3

4 5 6

> shape a

[2,3]

> rank a

2

> size a

6

> pp $ iota @10

0 1 2 3 4 5 6 7 8 9

Some examples: index, unScalar

> pp a

1 2 3

4 5 6

> pp $ index a 1

4 5 6

> pp $ a `index` 1 `index` 0

4

> :t index

index :: KnownNat s => Array (s : sh) a -> Int -> Array sh a

> print $ a `index` 1 `index` 0

fromList @[] [4]

> print $ unScalar $ a `index` 1 `index` 0

4

> :t unScalar

unScalar :: Array '[] a -> a

Some examples: transpose

> pp a

1 2 3

4 5 6

> pp $ transpose @[1,0] a

1 4

2 5

3 6

> transpose @[2,0] a

​

<interactive>:27:1: error:

• Could not deduce (Elem 2 '[]) arising from a use of ‘transpose’

from the context: (Num a, Enum a)

bound by the inferred type of

it :: (Num a, Enum a) => Array (Permute '[2, 0] '[2, 3]) a

at <interactive>:27:1-18

• In the expression: transpose @[2, 0] a

In an equation for ‘it’: it = transpose @[2, 0] a

Some examples: transpose

> b = fromList @[2,3,4] [1..24]

> pp b

1 2 3 4

5 6 7 8

9 10 11 12

​

13 14 15 16

17 18 19 20

21 22 23 24

​

> pp $ transpose @[2,1,0] b

1 13

5 17

9 21

​

2 14

6 18

10 22

​

3 15

7 19

11 23

​

4 16

8 20

12 24

​

​

Some examples: transpose

> :t transpose

transpose

:: (Permutation is, Shape sh, Shape is,

KnownNat (Rank sh), Rank is <= Rank sh) =>

Array sh a -> Array (Permute is sh) a

Some examples: reshape

> pp $ reshape @[3,2] a

1 2

3 4

5 6

> pp $ reshape @[3,3,3] iota

0 1 2

3 4 5

6 7 8

​

9 10 11

12 13 14

15 16 17

​

18 19 20

21 22 23

24 25 26

> :t reshape

reshape :: (Shape sh, Shape sh', Size sh ~ Size sh') =>

Array sh a -> Array sh' a

Some examples: broadcast

> b = reshape @[2,3] $ transpose @[1,0] a

> pp b

1 3 5

2 4 6

> :t a

Array [3,2] Integer

> pp $ (broadcast @[0,2] a :: Array [3,4,2] Integer)

1 2

1 2

1 2

1 2

3 4

3 4

3 4

3 4

5 6

5 6

5 6

5 6

Some examples: arithmetic

> pp $ a + a

2 4 6

8 10 12

> pp $ a + 1

2 3 4

5 6 7

> pp $ sqrt a

1.0 1.4142135623730951 1.7320508075688772

2.0 2.23606797749979 2.449489742783178

> c = fromList @'[2] [100,1000]

> pp c

100 1000

> pp $ a + broadcast @'[0] c

101 102 103

1004 1005 1006

Some examples: reduce

> pp a

1 2 3

4 5 6

> pp $ reduce (+) 0 a

21

> :t reduce

reduce :: Shape sh =>

(a -> a -> a) -> a -> Array sh a -> Array '[] a

> sumA = reduce (+) 0

> :t sumA

sumA :: (Shape sh, Num a) => Array sh a -> Array '[] a

> pp $ sumA a

21

Some examples: rerank

> pp a

1 2 3

4 5 6

> pp $ rerank @1 sumA a

6 15

> pp $ rerank @1 sumA $ transpose @[1,0] a

5 7 9

> pp b

1 2 3 4

5 6 7 8

9 10 11 12

​

13 14 15 16

17 18 19 20

21 22 23 24

> pp $ rerank @2 sumA b

10 26 42

58 74 90

Matrix multiply

• To multiply two matrices you dot product all rows in the first matrix with all columns in second matrix.

​

​

​

1 * 11 + 2 * 13 + 3 * 15 = 82

1 * 12 + 2 * 14 + 3 * 16 = 88

4 * 11 + 5 * 13 + 6 * 15 = 199

4 * 12 + 5 * 14 + 6 * 16 = 214

*

=

Some examples: matrix multiply

> pp a

1 2 3

4 5 6

> b = fromList @[3,2] [11..16]

> pp b

11 12

13 14

15 16

> pp $ transpose @[1,0] b

11 13 15

12 14 16

​

​

Some examples: matrix multiply

> as = broadcast @[0,2] a :: Array [2,2,3] Int

> pp as

1 2 3

1 2 3

​

4 5 6

4 5 6

> bs = broadcast @[1,2] (transpose @[1,0] b) :: Array [2,2,3] Int

> pp bs

11 13 15

12 14 16

​

11 13 15

12 14 16

> cs = as * bs

> pp cs

11 26 45

12 28 48

​

44 65 90

48 70 96

> ds = rerank @2 sumA cs

> pp ds

82 88

199 214

Some examples: matrix multiply

matMul :: (KnownNat n, KnownNat o, KnownNat m, Num a) =>

Array [m,n] a -> Array [n,o] a-> Array [m,o] a

matMul a b =

let as = broadcast @[0,2] a

bs = broadcast @[1,2] $ transpose @[1,0] b

in rerank @2 sumA (as * bs)

​

sumA :: (Shape sh, Num a) =>

Array sh a -> Array '[] a

sumA = reduce (+) 0

​

Some examples: 1-D convolution

• A 1-dimensional convolution takes a vector and a kernel and produces the dot product of the subvectors with the kernel.

​

Vector

​

Kernel

​

Subvectors

​

Mul by kern

​

Sum

​

Some examples: 1-D convolution

> pp a

1 2 3

4 5 6

> pp k

10 11 2

> pp (window @'[3] a)

1 2 3

2 3 4

3 4 5

4 5 6

> pp (broadcast @'[1] k :: Array [4,3] Int)

10 11 12

10 11 12

10 11 12

10 11 12

> pp (window @'[3] a * broadcast @'[1] k)

10 22 36

20 33 48

30 44 60

40 55 72

> pp $ rerank @1 sumA (window @'[3] a * broadcast @'[1] k)

68 101 134 167

​

Some examples: 1-D convolution

conv1D :: forall sk sa a .

(KnownNat sk, KnownNat sa,

sk <= sa, Num a) =>

Array '[sk] a -> Array '[sa] a -> Array '[sa+1-sk] a

conv1D k a =

rerank @1 asum (window @'[sk] a * broadcast @'[1] k)

​

​

​

Complexity

• Are all these operations expensive? No, thanks to a clever representation.
• Not my idea, has been in use since the 60s, “dope vectors”
• Allows any operation that is an affine transformation on the index space in constant time.
• Constant time operations:
• fromVector, stretch, constant, transpose, window, stride, slice, rev, index
• Linear or constant time operations:
• reshape, toVector, rerank
• Linear time operations:
• fromList, toList, map, zipWith, append, reduce

Variations, typing

• The library comes in three flavors with varying degree of strong typing:
• Shaped, the shape is part of the type, e.g., Array [3,2] Int
• Ranked, the rank is part of the type, e.g., Array 2 Int
• Dynamic, the shape is not part of the type, e.g., Array Int
• The less strongly typed versions do runtime type checks.
• The APIs are as similar as possible, e.g.,
• Shaped, fromList @[3,2] xs
• Dynamic, fromList [3,2] xs

Variations, storage

• To get efficient array operations it can be important to use unboxed numbers.
• The library comes in four variations for the underlying storage.
• Regular: boxed from Data.Vector
• Unboxed: unboxed from Data.Vector.Unboxed
• Storable: unboxed, pinned from Data.Vector.Storable
• Generic: any storage
• This means that there is a total of 3*4 = 12 variations of the library.
• This code duplication is really a flaw in the abstraction power of Haskell.

Implementation

• The implementation uses a vector for storing the elements and
• for each dimension, its size
• for each dimension, the step size in the vector
• an offset into the vector

Array

Sizes

Steps

Offset

Vector

d₁

d₂

d₃

.

.

.

s₁

s₂

s₃

.

.

.

o

• The position (in the vector) of an element with indices i₁, i₂, i₃, … is�o + i₁*s₁ + i₂*s₂ + i₃*s₃ + …

Implementation, example

• The same old example, the 2x3 matrix
• 1 2 3�4 5 6

Array

Sizes

Steps

Offset

Vector

2

3

3

1

0

Implementation, example

• The constant 3x4 matrix
• 7 7 7 7�7 7 7 7�7 7 7 7

Array

Sizes

Steps

Offset

Vector

3

4

0

0

0

Implementation, example, transpose

• The 2x3 matrix
• 1 2 3�4 5 6

Transpose

Sizes

Steps

Offset

Vector

3

2

1

3

1

• The transpose of it
• 1 4�2 5�3 6

Implementation, example, slice

• The 2x3 matrix
• 1 2 3�4 5 6

Slice

Sizes

Steps

Offset

Vector

2

2

3

1

1

• A 2x2 slice of it
• 2 3�5 6

Implementation, drawbacks

Main problem: loss of locality.

• After many of the operations "scramble" the step sizes and so elements that logically are adjacent are not physically adjacent.
• The main culprit is reshape

Complexity again

• All operations that can be expressed as an affine transformation on the indices can be done in constant time.
• stretch, constant, transpose, window, stride, slice, rev, index, (some) reshape, ...
• Affine transformations are those that can be expressed as
• A * I + B�where I is the vector of indices and A, B are matrices.
• Multiple affine transforms can be combined because matrix multiplication is associative.

Conclusion

• Orthotope is available at https://hackage.haskell.org/package/orthotope
• Matrix multiplication is available at https://hackage.haskell.org/package/orthotope-hmatrix
• Uses hmatrix, which in turn uses BLAS and LAPACK.
• Orthotope together with matrix multiplication can be used for reasonably efficient implementation of neural nets.