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