## 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 (s1 ⋅ u1 + ⋯ + sn ⋅ un) = s1 ⋅ f u1 + ⋯ + sn ⋅ f un

Taking the ui 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.

### One Comment

1. #### 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 […]