## Simpler, more efficient, functional linear maps

A previous post described a data type of functional linear maps. As Andy Gill pointed out, we had a heck of a time trying to get good performance. This note describes a new representation that is *very simple* and much more efficient. It’s terribly obvious in retrospect but took me a good while to stumble onto.

The Haskell module described here is part of the vector-space library (version 0.5 or later) and requires ghc version 6.10 or better (for associated types).

**Edits**:

- 2008-11-09: Changed remarks about versions. The vector-space version 0.5 depends on ghc 6.10.
- 2008-10-21: Fixed the vector-space library link in the teaser.

### Linear maps

Semantically, a *linear map* is a function `f ∷ a → b`

such that, for all scalar values `s`

and "vectors" `u, v ∷ a`

, the following properties hold:

*f* (*s* ⋅ *u*) = *s* ⋅ *f* *u*

*f* (*u* + *v*) = *f* *u* + *f* *v*

By repeated application of these properties,

*f* (*s*_{1} ⋅ *u*_{1} + ⋯ + *s*_{n} ⋅ *u*_{n}) = *s*_{1} ⋅ *f* *u*_{1} + ⋯ + *s*_{n} ⋅ *f* *u*_{n}

Taking the *u _{i}* as basis vectors, this form implies that a linear function is determined by its behavior on any basis of its domain type.

Therefore, a linear function can be represented simply as a function from a basis, using the representation described in *Vector space bases via type families*.

`type u :-* v = Basis u → v`

The semantic function converts from `(u :-* v)`

to `(u → v)`

. It decomposes a source vector into its coordinates, applies the basis function to basis representations, and linearly combines the results.

`lapply ∷ ( VectorSpace u, VectorSpace v`

, Scalar u ~ Scalar v, HasBasis u ) ⇒

(u :-* v) → (u → v)

lapply lm = λ u → sumV [s *^ lm b | (b,s) ← decompose u]

or

`lapply lm = linearCombo ∘ fmap (first lm) ∘ decompose`

The reverse function is easier. Convert a function `f`

, presumed linear, to a linear map representation:

`linear ∷ (VectorSpace u, VectorSpace v, HasBasis u) ⇒`

(u → v) → (u :-* v)

It suffices to apply `f`

to basis values:

`linear f = f ∘ basisValue`

### Memoization

The idea of the linear map representation is to reconstruct an entire (linear) function out of just a few samples. In other words, we can make a very small sampling of function's domain, and re-use those values in order to compute the function's value at *all* domain values. As implemented above, however, this trick makes function application more expensive, not less. If `lm = linear f`

, then each use of `lapply lm`

can apply `f`

to the value of every basis element, and then linearly combine results.

A simple trick fixes this efficiency problem: *memoize* the linear map. We could do the memoization privately, e.g.,

`linear f = memo (f ∘ basisValue)`

If `lm = linear f`

, then no matter how many times `lapply lm`

is applied, the function `f`

can only get applied as many times as the dimension of the domain of `f`

.

However, there are several other ways to make linear maps, and it would be easy to forget to memoize each combining form. So, instead of the function representation above, I ensure that the function be memoized by representing it as a memo trie.

`type u :-* v = Basis u ↛ v`

The conversion functions `linear`

and `lapply`

need just a little tweaking. Split `memo`

into its definition `untrie ∘ trie`

, and then move the second phase (`untrie`

) into `lapply`

. We'll also have to add `HasTrie`

constraints:

`linear ∷ ( VectorSpace u, VectorSpace v`

, HasBasis u, HasTrie (Basis u) ) ⇒

(u → v) → (u :-* v)

linear f = trie (f ∘ basisValue)

lapply ∷ ( VectorSpace u, VectorSpace v, Scalar u ~ Scalar v

, HasBasis u, HasTrie (Basis u) ) ⇒

(u :-* v) → (u → v)

lapply lm = linearCombo ∘ fmap (first (untrie lm)) ∘ decompose

Now we can build up linear maps conveniently and efficiently by using the operations on memo tries shown in *Composing memo tries*. For instance, suppose that `h`

is a linear function of two arguments (linear in *both*, not it *each*) and `m`

and `n`

are two linear maps. Then `liftA2 h m n`

is the linear function that applies `h`

to the results of `m`

and `n`

.

`lapply (liftA2 h m n) a = h (lapply m a) (lapply n a)`

Exploiting the applicative functor instance for functions, we get another formulation:

`lapply (liftA2 h m n) = liftA2 h (lapply m) (lapply n)`

In other words, the meaning of a `liftA2`

is the `liftA2`

of the meanings, as discussed in *Simplifying semantics with type class morphisms*.

## Conal Elliott » Blog Archive » Comparing formulations of higher-dimensional, higher-order derivatives:

[...] implementation described in Higher-dimensional, higher-order derivatives, functionally and in Simpler, more efficient, functional linear maps. For the case of Rn, my formulation uses a trie with entries for n basis elements, while [...]

24 January 2009, 6:43 pm