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 [...]

Leave a comment