## Notes on Finally Tagless

For reading group this week we read finally tagless partially evaluated http://okmij.org/ftp/tagless-final/JFP.pdf. It took me a couple minutes to regain my footing on my preferred explanation of what tagless is. This is why we write blog posts, to remember that kind of thing.

One thing that is very confusing about finally tagless as presented is that people tend to be talking about dsls with binding forms, like embedded lambda calculi, or tensor summation and things. This is complicated and I think to some degree orthogonal to the core of the the idea. Instead I’ll use Bool as my running example, which is so simple that it perhaps obscures the idea in the opposite direction.

When you define a data type, you define constructors. Constructors are just functions. This is more readily apparent using GADT syntax.

What makes constructor feel like more than just ordinary functions is that you can pattern match out of them too. Applying constructors and pattern matching out of them is a completely lossless process. The two processes are dual in some sense. In some sense, it seems like programming is one big shuffling game. In some sense. In some sense. In some sense.

In some sense. iN SoME SeNSe

Anyway, pattern matching is it’s own thing that doesn’t feel like other piece of the language. But pattern matching can be captured as a first class object with the notion of an eliminator / recursor function. If you think about it, what pattern matching is is a thing that takes that data type and then gives you the stuff inside the data type. So pattern matching is the same as function that takes in a functions that tell me what to do with that stuff for each case.

The bohm-berarducci encoding of data types makes the pattern matching function the data type itself.

data BoolI where
TrueI  :: BoolI
FalseI :: BoolI

type BoolC = forall s. s -> s -> s
truec :: BoolC
truec = \x y -> x
falsec :: BoolC
falsec = \x y -> y

to :: BoolI -> BoolC
to TrueI = truec
to FalseI = falsec

from :: BoolC -> BoolI
from f = f TrueI FalseI


In the final encoding of the datatype, we replace the data keyword with the class keyword. We can witness the isomorphism with an instance for BoolI and an intepretation function from BoolI to BoolF

class BoolF a where
truef :: a
falsef :: a

instance BoolF BoolI where
truef = TrueI
falsef = FalseI

interpf :: BoolF a => BoolI -> a
interpf TrueI = truef
interpf FalseI = falsef

However, there are some very nice features of this encoding. Typeclass resolution happens completely at compile time. This means that you can write something once and have it run many ways, with no runtime overhead. This is useful for dsls, especially ones you intend to just immediately interpret out of.

Once way you can have something run many ways is by having a syntax tree for the thing you want to do. Then you can write different intepreters. But then you have the runtime cost of interpretation.

interpi :: BoolI -> Int
interpi TrueI = 40
interpi FalseF = 27

interps :: BoolI -> String
interps TrueI = "hi"
interps FalseF = "fred"

instance BoolF Int where
truef = 40
falsef = 27

instance BoolF String where
truef = "hi"
falsef = "fred" 

A second feature is the openness of typeclasses compared to data types. Suppose you wanted to add another field to BoolI. Now you need to correct all your functions. However, you can make the new field a new typeclass and all your old functions still work. You can require the power you need.

A third thing is that finally tagless does get you some of the type restriction available with GADTs in a language without them. GADTs are IN SOME SENSE just constructors without the most general inferred type. But they also let you recover the type information you hid away upon pattern matching.

We can see the correspondence in a different way. A typeclass constraint corresponds to the implicit supplying of a dictionary with fields correspond to the typeclass.

s -> s -> s ~ (s,s) -> s ~ {truef :: s, falsef :: s} -> s ~ BoolF s => s

What is finally tagless not so good for? Brains mostly. It is quite a mind bending style to use. If you want to do deep pattern matching in some simplifier, it is possible, yet rather difficult to achieve. I’ve seen this done in some Oleg papers somewhere (on SQL query optimization I think?)

Here’s another example on list

class ListF f where
cons :: a -> f a -> f a
nil :: f a

instance ListF [] where
cons = (:)
nil = []

interpl :: LiftF f => [a] -> f a
interpl (x : xs) = cons x (interpl xs)
interpl [] = nil

type ListC a = forall s. (a -> s -> s) -> s -> s

Going the other direction from finally tagless is interesting as well. If you take a typeclass and replace the keyword class with data, you get something like the “free” version of that class. Two cases in mind are that of the free monoid and free monad. The usual presentation of these looks different though. That is because they are canonized. These data types need to be thought of as “modulo” the laws of the typeclass, which probably shows up in a custom Eq instance. I’m a little hazy on exactly how to explain the Pure constructors, but you do need them I think.

data FreeMonoid a where
Mappend :: FreeMonoid a -> FreeMonoid a -> FreeMonoid a
Mempty :: FreeMonoid a
Pure :: a -> Freemonoid a
Bind :: FreeMond f a -> (a -> FreeMonad f b) -> FreeMonad f b
Return :: a -> FreeMonad f a
Pure' :: f a -> FreeMonad f a


http://okmij.org/ftp/tagless-final/JFP.pdf – tagless final paper. Also some very interesting things related to partial evaluation

https://oleg.fi/gists/posts/2019-06-26-linear-church-encodings.html – interesting explanation of bohm-berarducci

http://okmij.org/ftp/tagless-final/course/lecture.pdf – oleg’s course

## Solving the Laplace Equations with Linear Relations

The Laplace equation is ubiquitous in physics and engineering.

$latex \nabla^2 \phi = 0 = \partial_x^2 \phi + \partial_y^2 \phi = 0 It and slight variants of it describes electrostatics, magnetostatics, steady state heat flow, elastic flex, pressure, velocity potentials. There are a couple reasons for that. • It results from minimizing the squared gradient of a field $|\nabla \phi |^2$ which can make sense from an energy minimization perspective. • Similarly it results from the combination of a flow conservation law and a linear constitutive relation connecting flow and field (such as Ohm’s law, Fick’s law, or Hooke’s law). • It also gets used even if not particularly appropriate because we know how to mathematically deal with it, for example in image processing. There are a couple of questions we may want to ask about a Laplace equation system • Given the field on the boundary, determine the field in the interior (Dirchlet problem) • Given the normal derivative of the field on the boundary determine the field in the interior (Neumann problem) • Given sources in the interior and 0 boundary condition, determine the field. The Laplace equation is called the Poisson equation when you allow a source term on the right hand side. $\nabla^2 \phi = \rho$. • Given the field at the boundary, determine the derivative at the boundary. Dirichlet-to-Neumann map or Poincare-Steklov operator. Given the Dirichlet to Neumann map, you do not have to consider the interior of a region to use it. The Dirichlet to Neumann map is sort of the same thing as an effective resistance or scattering matrix. It gives you a black box representation of a region based solely on the variables at its boundary. This linear relation algebra is useful for many things that I’d have considered a use case for the Schur complement. The Schur complement arises when you do Gaussian elimination on a blocked matrix. It is good that this pattern has a name, because once you know about it, you’ll see it in many places. Domain decomposition, marginalized gaussian distributions, low-rank update, Scattering matrices. By composing the linear relations corresponding to the Dirchlet-Neumann relations of regions, we can build the Dirichlet-Neumann relations of larger regions. To make this more concrete, let us take the example of electrical circuits like before. A grid of resistors is a finite difference approximation to the continuous problem $-\nabla \phi = E$ Electric field is gradient of potential $E = \rho j$ continuum ohm’s law $\nabla\cdot j = 0$ current conservation In this post, I mentioned how you can make reasonable 2 dimensional diagrams out of a monoidal category, by sort of arbitrarily flipping one wire up and one wire down as in the diagram below. This defines a horizontal and vertical composition which have to do the required book-keeping (associations) to keep an arrow in canonical form. I had considered this for the management method of weights in neural networks, but it is way more natural as the actual geometrical layout of a finite difference grid of a laplace equation. So we can reuse our categorical circuit combinators to build a finite difference Laplace equation. This can be implemented in Haskell doing the following. Neato.  type HLinRel2D u d l r = HLinRel (Either u l) (Either d r) {- A stencil of 2d resistors for tiling u / \ / | l -/\/\/----/\/\/\- r | / \ / | d -} stencil :: HLinRel2D IV IV IV IV stencil = (hpar r10 r10) <<< short <<< (hpar r10 r10) where r10 = resistor 10 horicomp :: forall w w' w'' w''' a b c. (BEnum w, BEnum w', BEnum a, BEnum w''', BEnum b, BEnum w'', BEnum c ) => HLinRel2D w' w'' b c -> HLinRel2D w w''' a b -> HLinRel2D (Either w' w) (Either w'' w''') a c horicomp f g = hcompose f' g' where f' :: HLinRel (Either (Either w' w''') b) (Either (Either w'' w''') c) f' = (first hswap) <<< hassoc' <<< (hpar hid f) <<< hassoc <<< (first hswap) g' :: HLinRel (Either (Either w' w) a) (Either (Either w' w''') b) g' = hassoc' <<< (hpar hid g) <<< hassoc rotate :: (BEnum w, BEnum w', BEnum a, BEnum b) => HLinRel2D w w' a b -> HLinRel2D a b w w' rotate f = hswap <<< f <<< hswap vertcomp :: (BEnum w, BEnum w', BEnum a, BEnum d, BEnum b, BEnum w'', BEnum c ) => HLinRel2D w' w'' c d -> HLinRel2D w w' a b -> HLinRel2D w w'' (Either c a) (Either d b) vertcomp f g = rotate (horicomp (rotate f) (rotate g) ) view raw laplace.hs hosted with ❤ by GitHub ### Bits and Bobbles • Not the tightest post, but I needed to get it out there. I have a short attention span. • Homology and defining simplices as categories. One way of describing Homology is about linear operators that are analogues of finite difference operators (or better yet, discrete differential geometry operators / exterior derivatives). To some degree, it is analyzing the required boundary conditions to fully define differential equations on weirdo topological surfaces, which correspond to geometrical loops/holes. You can figure this out by looking at subspaces and quotients of the difference operators. Here we have here a very category theory way of looking at partial differential equations. How does it all connect? • Continuous circuit models – https://en.wikipedia.org/wiki/Distributed-element_model Telegrapher’s equation is classic example. • Cody mentioned that I could actually build circuits and measure categorical identities in a sense. That’s kind of cool. Or I could draw conductive ink on carbon paper and actually make my string diagrams into circuits? That is also brain tickling • Network circuits • I really want to get coefficients that aren’t just doubles. allowing rational functions of a frequency $\omega$ would allow analysis of capacitor/inductor circuits, but also tight binding model systems for fun things like topological insulators and the Haldane model http://www.philipzucker.com/topologically-non-trivial-circuit-making-haldane-model-gyrator/ . I may need to leave Haskell. I’m not seeing quite the functionality I need. Use Sympy? https://arxiv.org/abs/1605.02532 HLinear. Flint bindings for haskell? Looks unmaintained. Could also use a grobner basis package as dynamite for a mouse. • This is relevant for the boundary element method. Some really cool other stuff relevant here. http://people.maths.ox.ac.uk/martinsson/2014_CBMS/ Blah Blah blah: The subtlest aspect of differential equations is that of boundary conditions. It is more correct usually to consider the system of the interior differential equation and the boundary conditions as equally essential parts of the statement of the problem you are considering. ## Categorical LQR Control with Linear Relations Last time, I described and built an interesting expansion of describing systems via linear functions, that of using linear relations. We remove the requirement that mappings be functional (exactly one output vector for any given input vector), which extends our descriptive capability. This is useful for describing “open” pieces of a linear system like electric circuits. This blog post is about another application: linear dynamical systems and control. Linear relations are an excellent tool for reachability/controllability. In a control system, it is important to know where it is possible to get the system to. With linear dynamics $x_{t+1} = Ax_{t} + B u_{t}$, an easy controllability question is “can my system reach into some subspace”. This isn’t quite the most natural question physically, but it is natural mathematically. My first impulse would be to ask “given the state starts in this little region over here, can it get to this little region over here”, but the subspace question is a bit easier to answer. It’s a little funky but it isn’t useless. It can detect if the control parameter only touches 1 of 2 independent dynamical systems for example. We can write the equations of motion as a linear relation and iterate composition on them. See Erbele for more. There is a set of different things (and kind of more interesting things) that are under the purview of linear control theory, LQR Controllers and Kalman filters. The LQR controller and Kalman filter are roughly (exactly?) the same thing mathematically. At an abstract mathematical level, they both rely on the fact that optimization of a quadratic objective $x^T Q x + r^T x + c$ with a linear constraints $Ax=b$ is solvable in closed form via linear algebra. The cost of the LQR controller could be the squared distance to some goal position for example. When you optimize a function, you set the derivative to 0. This is a linear equation for quadratic objectives. It is a useful framework because it has such powerful computational teeth. The Kalman filter does a similar thing, except for the problem of state estimation. There are measurements linear related to the state that you want to match with the prior knowledge of the linear dynamics of the system and expected errors of measurement and environment perturbations to the dynamics. There are a couple different related species of these filters. We could consider discrete time or continuous time. We can also consider infinite horizon or finite horizon. I feel that discrete time and finite horizon are the most simple and fundamental so that is what we’ll stick with. The infinite horizon and continuous time are limiting cases. We also can consider the dynamics to be fixed for all time steps, or varying with time. Varying with time is useful for the approximation of nonlinear systems, where there are different good nonlinear approximation (taylor series of the dynamics) depending on the current state. There are a couple ways to solve a constrained quadratic optimization problem. In some sense, the most conceptually straightforward is just to solve for a span of the basis (the “V-Rep” of my previous post) and plug that in to the quadratic objective and optimize over your now unconstrained problem. However, it is commonplace and elegant in many respects to use the Lagrange multiplier method. You introduce a new parameter for every linear equality and add a term to the objective $\lambda^T (A x - b)$. Hypothetically this term is 0 for the constrained solution so we haven’t changed the objective in a sense. It’s a funny slight of hand. If you optimize over x and $\lambda$, the equality constraint comes out as your gradient conditions on $\lambda$. Via this method, we convert a linearly constrained quadratic optimization problem into an unconstrained quadratic optimization problem with more variables. Despite feeling very abstract, the value these Lagrangian variables take on has an interpretation. They are the “cost” or “price” of the equality they correspond to. If you moved the constraint by an amount $Ax - b + \epsilon$, you would change the optimal cost by an amount $\lambda \epsilon$. (Pretty sure I have that right. Units check out. See Boyd https://www.youtube.com/watch?v=FJVmflArCXc) The Lagrange multipliers enforcing the linear dynamics are called the co-state variables. They represent the “price” that the dynamics cost, and are derivatives of the optimal value function V(s) (The best value that can be achieved from state s) that may be familiar to you from dynamic programming or reinforcement learning. See references below for more. Let’s get down to some brass tacks. I’ve suppressed some terms for simplicity. You can also add offsets to the dynamics and costs. A quadratic cost function with lagrange multipliers. Q is a cost matrix associated with the x state variables, R is a cost matrix for the u control variables. $C = \lambda_0 (x_0 - \tilde{x}_0) + \sum_t x_t^T Q x_t + u_t^T R u_t + \lambda_{t+1}^T (x_{t+1} - A x_t - B u_t )$ Equations of motion results from optimizing with respect to $\lambda$ by design. $\nabla_{\lambda_{t+1}} C = x_{t+1} - A x_t - B u_t = 0$. The initial conditions are enforced by the zeroth multiplier. $\nabla_{\lambda_0} C = x_i - x_{0} = 0$ Differentiation with respect to the state x has the interpretation of backwards iteration of value derivative, somewhat related to what one finds in the Bellman equation. $\nabla_{x_t} C = Q x_{t} + A^T \lambda_{t+1} - \lambda_{t} = 0 \Longrightarrow \lambda_{t} = A^T \lambda_{t+1} + Q x_{t}$ The final condition on value derivative is the one that makes it the most clear that the Lagrange multiplier has the interpretation of the derivative of the value function, as it sets it to that. $\nabla_{x_N} C = Q x_N - \lambda_{N} = 0$ Finally, differentiation with respect to the control picks the best action given knowledge of the value function at that time step. $\nabla_{u_t} C = R u_{t} - B^T \lambda_t = 0$ Ok. Let’s code it up using linear relations. Each of these conditions is a separate little conceptual box. We can build the optimal control step relation by combining these updates together The following is a bit sketchy. I am confident it makes sense, but I’m not confident that I have all of the signs and other details correct. I also still need to add the linear terms to the objective and offset terms to the dynamics. These details are all kind of actually important. Still, I think it’s nice to set down in preliminary blog form.  -- state of an oscillator data SHOState = X | P deriving (Show, Enum, Bounded, Eq, Ord) data Control = F deriving (Show, Enum, Bounded, Eq, Ord) -- Costate newtype wrapper newtype Co a = Co a deriving (Show, Enum, Bounded, Eq, Ord) type M = Matrix Double dynamics :: forall x u. (BEnum x, BEnum u) => Matrix Double -> Matrix Double -> HLinRel (Either x u) x dynamics a b = HLinRel (a ||| b ||| -i ) (vzero cx) where cx = card @x cu = card @u i = ident cx initial_cond :: forall x. BEnum x => Vector Double-> HLinRel Void x initial_cond x0 = HLinRel i x0 where cx = card @x i = ident cx valueUpdate :: forall x l. (BEnum x, BEnum l) => M -> M -> HLinRel (Either x l) l valueUpdate a q = HLinRel ((tr a) ||| q ||| i) (vzero cl) where cl = card @l i = ident cl optimal_u :: forall u l. (BEnum u, BEnum l) => M -> M -> HLinRel u l optimal_u r b = HLinRel (r ||| tr b) (vzero cu) where cu = card @u step :: forall x u l. (BEnum x, BEnum u, BEnum l) => M -> M -> M -> M -> HLinRel (Either x l) (Either x l) step a b r q = f5 <<< f4 <<< hassoc' <<< f3 <<< f2 <<< hassoc <<< f1 where f1 :: HLinRel (Either x l) (Either (Either x x) l) f1 = first hdup f2 :: HLinRel (Either x (Either x l)) (Either x l) f2 = second (valueUpdate a q) f3 :: HLinRel (Either x l) (Either x (Either l l)) f3 = second hdup f4 :: HLinRel (Either (Either x l) l) (Either (Either x u) l) f4 = first (second (hconverse (optimal_u r b))) f5 :: HLinRel (Either (Either x u) l) (Either x l) f5 = first (dynamics a b) -- iterate (hcompose (step a b r q)) :: [HLinRel (Either x l) (Either x l)] view raw optimal_relation.hs hosted with ❤ by GitHub ### Bits and Bobbles • The code in context https://github.com/philzook58/ConvexCat/blob/master/src/LinRel.hs • Some of the juicier stuff is nonlinear control. Gotta walk before we can run though. I have some suspicions that a streaming library may be useful here, or a lens-related approach. Also ADMM. • Reminiscent of Keldysh contour. Probably a meaningless coincidence. • In some sense H-Rep is the analog of (a -> b -> Bool) and V-rep [(a,b)] • Note that the equations of motion are relational rather than functional for a control systems. The control parameters $u$ describe undetermined variables under our control. • Loopback (trace) the value derivative for infinite horizon control. • Can solve the Laplace equation with a Poincare Steklov operator. That should be my next post. • There is something a touch unsatisfying even with the reduced goals of this post. Why did I have to write down the quadratic optimization problem and then hand translate it to linear relations? It’s kind of unacceptable. The quadratic optimization problem is really the most natural statement, but I haven’t figured out to make it compositional. The last chapter of Rockafellar? ### References https://arxiv.org/abs/1405.6881 Baez and Erbele – Categories in Control http://www.scholarpedia.org/article/Optimal_control http://www.argmin.net/2016/05/18/mates-of-costate/ – interesting analogy to backpropagation. Lens connection? https://courses.cs.washington.edu/courses/cse528/09sp/optimality_chapter.pdf ## Linear Relation Algebra of Circuits with HMatrix Oooh this is a fun one. I’ve talked before about relation algebra and I think it is pretty neat. http://www.philipzucker.com/a-short-skinny-on-relations-towards-the-algebra-of-programming/. In that blog post, I used finite relations. In principle, they are simple to work with. We can perform relation algebra operations like composition, meet, and join by brute force enumeration. Unfortunately, brute force may not always be an option. First off, the finite relations grow so enormous as to be make this infeasible. Secondly, it is not insane to talk about relations or regions with an infinite number of elements, such as some continuous blob in 2D space. In that case, we can’t even in principle enumerate all the points in the region. What are we to do? We need to develop some kind of finite parametrization of regions to manipulate. This parametrization basically can’t possibly be complete in some sense, and we may choose more or less powerful systems of description for computational reasons. In this post, we are going to be talking about linear or affine subspaces of a continuous space. These subspaces are hyperplanes. Linear subspaces have to go through the origin, while affine spaces can have an offset from the origin. In the previous post, I mentioned that the finite relations formed a lattice, with operations meet and join. These operations were the same as set intersection and union so the introduction of the extra terminology meet and join felt a bit unwarranted. Now the meet and join aren’t union and intersection anymore. We have chosen to not have the capability to represent the union of two vectors, instead we can only represent the smallest subspace that contains them both, which is the union closed under vector addition. For example, the join of a line and point will be the plane that goes through both. Linear/Affine stuff is great because it is so computational. Most questions you cant to ask are answerable by readily available numerical linear algebra packages. In this case, we’ll use the Haskell package HMatrix, which is something like a numpy/scipy equivalent for Haskell. We’re going to use type-level indices to denote the sizes and partitioning of these spaces so we’ll need some helper functions. In case I miss any extensions, make typos, etc, you can find a complete compiling version here https://github.com/philzook58/ConvexCat/blob/master/src/LinRel.hs type BEnum a = (Enum a, Bounded a) -- cardinality. size was already taken by HMatrix :( card :: forall a. (BEnum a) => Int card = (fromEnum (maxBound @a)) - (fromEnum (minBound @a)) + 1 In analogy with sets of tuples for defining finite relations, we partition the components of the linear spaces to be “input” and “output” indices/variables $\begin{bmatrix} x_1 & x_2 & x_3 & ... & y_1 & y_2 & y_3 & ... \end{bmatrix}$. This partition is somewhat arbitrary and easily moved around, but the weakening of strict notions of input and output as compared to functions is the source of the greater descriptive power of relations. Relations are extensions of functions, so linear relations are an extension of linear maps. A linear map has the form $y = Ax$. A linear relation has the form $Ax + By = 0$. An affine map has the form $y = Ax + b$ and an affine relation has the form $Ax + By = b$. There are at least two useful concrete representation for subspaces. 1. We can write a matrix $A$ and vector $b$ down that corresponds to affine constraints. $Ax = b$. The subspace described is the nullspace of $A$ plus a solution of the equation. The rows of A are orthogonal to the space. 2. We can hold onto generators of subspace. $x = A' l+b$ where l parametrizes the subspace. In other words, the subspace is generated by / is the span of the columns of $A'$. It is the range of $A'$. We’ll call these two representations the H-Rep and V-Rep, borrowing terminology from similar representations in polytopes (describing a polytope by the inequalities that define it’s faces or as the convex combination of it’s vertices). https://inf.ethz.ch/personal/fukudak/lect/pclect/notes2015/PolyComp2015.pdf These two representations are dual in many respects. -- HLinRel holds A x = b constraint data HLinRel a b = HLinRel (Matrix Double) (Vector Double) deriving Show -- x = A l + b. Generator constraint. data VLinRel a b = VLinRel (Matrix Double) (Vector Double) deriving Show It is useful to have both reps and interconversion routines, because different operations are easy in the two representations. Any operations defined on one can be defined on the other by sandwiching between these conversion functions. Hence, we basically only need to define operations for one of the reps (if we don’t care too much about efficiency loss which, fair warning, is out the window for today). The bulk of computation will actually be performed by these interconversion routines. The HMatrix function nullspace performs an SVD under the hood and gathers up the space with 0 singular values. -- if A x = b then x is in the nullspace + a vector b' solves the equation h2v :: HLinRel a b -> VLinRel a b h2v (HLinRel a b) = VLinRel a' b' where b' = a <\> b -- least squares solution a' = nullspace a -- if x = A l + b, then A' . x = A' A l + A' b = A' b because A' A = 0 v2h :: VLinRel a b -> HLinRel a b v2h (VLinRel a' b') = HLinRel a b where b = a #> b' -- matrix multiply a = tr$ nullspace (tr a') -- orthogonal space to range of a.
-- tr is transpose and not trace? A little bit odd, HMatrix.

These linear relations form a category. I’m not using the Category typeclass because I need BEnum constraints hanging around. The identity relations is $x = y$ aka $Ix - Iy = 0$.

hid :: forall a. BEnum a => HLinRel a a
hid =  HLinRel (i ||| (- i)) (vzero s) where
s = card @a
i = ident s

Composing relations is done by combining the constraints of the two relations and then projecting out the interior variables. Taking the conjunction of constraints is easiest in the H-Rep, where we just need to vertically stack the individual constraints. Projection easily done in the V-rep, where you just need to drop the appropriate section of the generator vectors. So we implement this operation by flipping between the two.

hcompose :: forall a b c. (BEnum a, BEnum b, BEnum c) => HLinRel b c -> HLinRel a b -> HLinRel a c
hcompose (HLinRel m b) (HLinRel m' b') = let a'' = fromBlocks [[       ma',           mb' ,    0       ],
[         0 ,    mb,        mc          ]] in
let b'' = vjoin [b', b] in
let (VLinRel q p) = h2v (HLinRel a'' b'') in -- kind of a misuse
let q' = (takeRows ca q)  -- drop rows belonging to @b
===
(dropRows (ca + cb) q) in
let [x,y,z] =  takesV [ca,cb,cc] p in
let p'=  vjoin [x,z] in -- rebuild without rows for @b
v2h (VLinRel q' p') -- reconstruct HLinRel
where
ca = card @a
cb = card @b
cc = card @c
sb = size b -- number of constraints in first relation
sb' = size b' -- number of constraints in second relation
ma' = takeColumns ca m'
mb' = dropColumns ca m'
mb = takeColumns cb m
mc = dropColumns cb m

(<<<) :: forall a b c. (BEnum a, BEnum b, BEnum c) => HLinRel b c -> HLinRel a b -> HLinRel a c
(<<<) = hcompose

We can implement the general cadre of relation operators, meet, join, converse. I feel the converse is the most relational thing of all. It makes inverting a function nearly a no-op.

hjoin :: HLinRel a b -> HLinRel a b -> HLinRel a b
hjoin v w = v2h $vjoin' (h2v v) (h2v w) -- hmatrix took vjoin from me :( -- joining means combining generators and adding a new generator -- Closed under affine combination l * x1 + (1 - l) * x2 vjoin' :: VLinRel a b -> VLinRel a b -> VLinRel a b vjoin' (VLinRel a b) (VLinRel a' b') = VLinRel (a ||| a' ||| (asColumn (b - b'))) b -- no constraints, everything -- trivially true htop :: forall a b. (BEnum a, BEnum b) => HLinRel a b htop = HLinRel (vzero (1,ca + cb)) (konst 0 1) where ca = card @a cb = card @b -- hbottom? hconverse :: forall a b. (BEnum a, BEnum b) => HLinRel a b -> HLinRel b a hconverse (HLinRel a b) = HLinRel ( (dropColumns ca a) ||| (takeColumns ca a)) b where ca = card @a cb = card @b  Relational inclusion is the question of subspace inclusion. It is fairly easy to check if a VRep is in an HRep (just see plug the generators into the constraints and see if they obey them) and by using the conversion functions we can define it for arbitrary combos of H and V.  -- forall l. A' ( A l + b) == b' -- is this numerically ok? I'm open to suggestions. vhsub :: VLinRel a b -> HLinRel a b -> Bool vhsub (VLinRel a b) (HLinRel a' b') = (naa' <= 1e-10 * (norm_2 a') * (norm_2 a) ) && ((norm_2 ((a' #> b) - b')) <= 1e-10 * (norm_2 b') ) where naa' = norm_2 (a' <> a) hsub :: HLinRel a b -> HLinRel a b -> Bool hsub h1 h2 = vhsub (h2v h1) h2 heq :: HLinRel a b -> HLinRel a b -> Bool heq a b = (hsub a b) && (hsub b a) instance Ord (HLinRel a b) where (<=) = hsub (>=) = flip hsub instance Eq (HLinRel a b) where (==) = heq It is useful the use the direct sum of the spaces as a monoidal product. hpar :: HLinRel a b -> HLinRel c d -> HLinRel (Either a c) (Either b d) hpar (HLinRel mab v) (HLinRel mcd v') = HLinRel (fromBlocks [ [mab, 0], [0 , mcd]]) (vjoin [v, v']) where hleft :: forall a b. (BEnum a, BEnum b) => HLinRel a (Either a b) hleft = HLinRel ( i ||| (- i) ||| (konst 0 (ca,cb))) (konst 0 ca) where ca = card @a cb = card @b i = ident ca hright :: forall a b. (BEnum a, BEnum b) => HLinRel b (Either a b) hright = HLinRel ( i ||| (konst 0 (cb,ca)) ||| (- i) ) (konst 0 cb) where ca = card @a cb = card @b i = ident cb htrans :: HLinRel a (Either b c) -> HLinRel (Either a b) c htrans (HLinRel m v) = HLinRel m v hswap :: forall a b. (BEnum a, BEnum b) => HLinRel (Either a b) (Either b a) hswap = HLinRel (fromBlocks [[ia ,0,0 ,-ia], [0, ib,-ib,0]]) (konst 0 (ca + cb)) where ca = card @a cb = card @b ia = ident ca ib = ident cb hsum :: forall a. BEnum a => HLinRel (Either a a) a hsum = HLinRel ( i ||| i ||| - i ) (konst 0 ca) where ca = card @a i= ident ca hdup :: forall a. BEnum a => HLinRel a (Either a a) hdup = HLinRel (fromBlocks [[i, -i,0 ], [i, 0, -i]]) (konst 0 (ca + ca)) where ca = card @a i= ident ca hdump :: HLinRel a Void hdump = HLinRel 0 0 hlabsorb ::forall a. BEnum a => HLinRel (Either Void a) a hlabsorb = HLinRel m v where (HLinRel m v) = hid @a  A side note: Void causes some consternation. Void is the type with no elements and is the index type of a 0 dimensional space. It is the unit object of the monoidal product. Unfortunately by an accident of the standard Haskell definitions, actual Void is not a BEnum. So, I did a disgusting hack. Let us not discuss it more. ### Circuits Baez and Fong have an interesting paper where they describe building circuits using a categorical graphical calculus. We have the pieces to go about something similar. What we have here is a precise way in which circuit diagrams can be though of as string diagrams in a monoidal category of linear relations. An idealized wire has two quantities associated with it, the current flowing through it and the voltage it is at. -- a 2d space at every wire or current and voltage. data IV = I | V deriving (Show, Enum, Bounded, Eq, Ord) When we connect wires, the currents must be conserved and the voltages must be equal. hid and hcompose from above still achieve that. Composing two independent circuits in parallel is achieve by hpar. We will want some basic tinker toys to work with. A resistor in series has the same current at both ends and a voltage drop proportional to the current resistor :: Double -> HLinRel IV IV resistor r = HLinRel ( (2><4) [ 1,0,-1, 0, r, 1, 0, -1]) (konst 0 2)  Composing two resistors in parallel adds the resistance. (resistor r1) <<< (resistor r2) == resistor (r1 + r2)) A bridging resistor allows current to flow between the two branches bridge :: Double -> HLinRel (Either IV IV) (Either IV IV) bridge r = HLinRel ( (4><8) [ 1,0, 1, 0, -1, 0, -1, 0, -- current conservation 0, 1, 0, 0, 0, -1 , 0, 0, --voltage maintained left 0, 0, 0, 1, 0, 0, 0, -1, -- voltage maintained right r, 1, 0,-1, -r, 0, 0, 0 ]) (konst 0 4)  Composing two bridge circuits is putting the bridge resistors in parallel. The conductance $G=\frac{1}{R}$ of resistors in parallel adds. hcompose (bridge r1) (bridge r2) == bridge 1 / (1/r1 + 1/r2). An open circuit allows no current to flow and ends a wire. open ~ resistor infinity open :: HLinRel IV Void open = HLinRel (fromList [[1,0]]) (konst 0 1) At branching points, the voltage is maintained, but the current splits. cmerge :: HLinRel (Either IV IV) IV cmerge = HLinRel (fromList [[1, 0, 1, 0, -1, 0], [0,1,0,0,0 ,-1 ], [0,0,0,1, 0, -1]]) (konst 0 3) This cmerge combinator could also be built using a short == bridge 0 , composing a branch with open, and then absorbing the Void away. We can bend wires up or down by using a composition of cmerge and open. cap :: HLinRel (Either IV IV) Void cap = hcompose open cmerge cup :: HLinRel Void (Either IV IV) cup = hconverse cap ground :: HLinRel IV Void ground = HLinRel ( (1><2) [ 0 , 1 ]) (vzero 1)  Voltage and current sources enforce current and voltage to be certain values vsource :: Double -> HLinRel IV IV vsource v = HLinRel ( (2><4) [ 1,0,-1, 0, 0, 1, 0, -1]) (fromList [0,v]) isource :: Double -> HLinRel IV IV isource i = HLinRel (fromList [ [1,0, -1, 0], -- current conservation [1, 0, 0, 0]]) (fromList [0,i])  Measurements of circuits proceed by probes. type VProbe = () vprobe :: HLinRel IV VProbe vprobe = HLinRel ( (2><3) [1,0,0, 0,1,-1]) (konst 0 2)  Inductors and capacitors could be included easily, but would require the entries of the HMatrix values to be polynomials in the frequency $\omega$, which it does not support (but it could!). We'll leave those off for another day. We actually can determine that the rules suggested above are being followed by computation. r20 :: HLinRel IV IV r20 = resistor 20 main :: IO () main = do print (r20 == (hid <<< r20)) print (r20 == r20 <<< hid) print (r20 == (hmeet r20 r20)) print$ resistor 50 == r20 <<< (resistor 30)
print $(bridge 10) <<< (bridge 10) == (bridge 5) print$ v2h (h2v r20) == r20
print $r20 <= htop print$ hconverse (hconverse r20) == r20
print $(open <<< r20) == open ### Bits and Bobbles • Homogenous systems are usually a bit more elegant to deal with, although a bit more unfamiliar and abstract. • Could make a pandas like interface for linear relations that uses numpy/scipy.sparse for the computation. All the swapping and associating is kind of fun to design, not so much to use. Labelled n-way relations are nice for users. • Implicit/Lazy evaluation. We should let the good solvers do the work when possible. We implemented our operations eagerly. We don't have to. By allowing hidden variables inside our relations, we can avoid the expensive linear operations until it is useful to actually compute on them. • Relational division = quotient spaces? • DSL. One of the beauties of the pointfree/categorical approach is that you avoid the need for binding forms. This makes for a very easily manipulated DSL. The transformations feel like those of ordinary algebra and you don't have to worry about the subtleties of index renaming or substitution under binders. • Sparse is probably really good. We have lots of identity matrices and simple rearrangements. It is very wasteful to use dense operations on these. • Schur complement https://en.wikipedia.org/wiki/Schur_complement are the name in the game for projecting out pieces of linear problems. We have some overlap. • Linear relations -> Polyhedral relations -> Convex Relations. Linear is super computable, polyhedral can blow up. Rearrange a DSL to abuse Linear programming as much as possible for queries. • Network circuits. There is an interesting subclass of circuits that is designed to be pretty composable. https://en.wikipedia.org/wiki/Two-port_network Two port networks are a very useful subclass of electrical circuits. They model transmission lines fairly well, and easily composable for filter construction. It is standard to describe these networks by giving a linear function between two variables and the other two variables. Depending on your choice of which variables depend on which, these are called the z-parameters, y-parameters, h-parameters, scattering parameters, abcd parameters. There are tables of formula for converting from one form to the others. The different parameters hold different use cases for composition and combining in parallel or series. From the perspective of linear relations this all seems rather silly. The necessity for so many descriptions and the confusing relationship between them comes from the unnecessary and overly rigid requirement of have a linear function-like relationship rather than just a general relation, which depending of the circuit may not even be available (there are degenerate configurations where two of the variables do not imply the values of the other two). A function relationship is always a lie (although a sometimes useful one), as there is always back-reaction of new connections. -- voltage divider divider :: Double -> Double -> HLinRel (Either IV IV) (Either IV IV) divider r1 r2 = hcompose (bridge r2) (hpar (resistor r1) hid)  The relation model also makes clearer how to build lumped models out of continuous ones. https://en.wikipedia.org/wiki/Lumped-element_model null • Because the type indices have no connection to the actual data types (they are phantom) it is a wise idea to use smart constructors that check that the sizes of the matrices makes sense.  -- smart constructors hLinRel :: forall a b. (BEnum a, BEnum b) => Matrix Double -> Vector Double -> Maybe (HLinRel a b) hLinRel m v | cols m == (ca + cb) && (size v == rows m) = Just (HLinRel m v) | otherwise = Nothing where ca = card @a cb = card @b  • Nonlinear circuits. Grobner Bases and polynomial relations? • Quadratic optimization under linear constraints. Can't get it to come out right yet. Clutch for Kalman filters. Nice for many formulations like least power, least action, minimum energy principles. Edit: I did more in this direction here http://www.philipzucker.com/categorical-lqr-control-with-linear-relations/ • Quadratic Operators -> Convex operators. See last chapter of Rockafellar. • Duality of controllers and filters. It is well known (I think) that for ever controller algorithm there is a filter algorithm that is basically the same thing. • LQR - Kalman • Viterbi filter - Value function table • particle filter - Monte Carlo control • Extended Kalman - iLQR-ish? Use local approximation of dynamics • unscented kalman - ? #### Refs ## Neural Networks with Weighty Lenses (DiOptics?) I wrote a while back how you can make a pretty nice DSL for reverse mode differentiation based on the same type as Lens. I’d heard some interesting rumblings on the internet around these ideas and so was revisiting them. type Lens s t a b = s -> (a, b -> t) type AD x dx y dy = x -> (y, dy -> dx) Composition is defined identically for reverse mode just as it is for lens. After chewing on it a while, I realized this really isn’t that exotic. How it works is that you store the reverse mode computation graph, and all necessary saved data from the forward pass in the closure of the (dy -> dx). I also have a suspicion that if you defunctionalized this construction, you’d get the Wengert tape formulation of reverse mode ad. Second, Lens is just a nice structure for bidirectional computation, with one forward pass and one backward pass which may or may not be getting/setting. There are other examples for using it like this. It is also pretty similar to the standard “dual number” form type FAD x dx y dy = (x,dx)->(y,dy) for forward mode AD. We can bring the two closer by a CPS/Yoneda transformation and then some rearrangement.  x -> (y, dy -> dx) ==> x -> (y, forall s. (dx -> s) -> (dy -> s)) ==> forall s. (x, dx -> s) -> (y, dx -> s)  and meet it in the middle with (x,dx) -> (y,dy) ==> forall s. (x, s -> dx) -> (y, s -> dy) I ended the previous post somewhat unsatisfied by how ungainly writing that neural network example was, and I called for Conal Elliot’s compiling to categories plugin as a possible solution. The trouble is piping the weights all over the place. This piping is very frustrating in point-free form, especially when you know it’d be so trivial pointful. While the inputs and outputs of layers of the network compose nicely (you no longer need to know about the internal computations), the weights do not. As we get more and more layers, we get more and more weights. The weights are in some sense not as compositional as the inputs and outputs of the layers, or compose in a different way that you need to maintain access to. I thought of a very slight conceptual twist that may help. The idea is we keep the weights out to the side in their own little type parameter slots. Then we define composition such that it composes input/outputs while tupling the weights. Basically we throw the repetitive complexity appearing in piping the weights around into the definition of composition itself. These operations are easily seen as 2 dimensional diagrams. Here’s the core reverse lens ad combinators import Control.Arrow ((***)) type Lens'' a b = a -> (b, b -> a) comp :: (b -> (c, (c -> b))) -> (a -> (b, (b -> a))) -> (a -> (c, (c -> a))) comp f g x = let (b, dg) = g x in let (c, df) = f b in (c, dg . df) id' :: Lens'' a a id' x = (x, id) relu' :: (Ord a, Num a) => Lens'' a a relu' = \x -> (frelu x, brelu x) where frelu x | x > 0 = x | otherwise = 0 brelu x dy | x > 0 = dy | otherwise = 0 add' :: Num a => Lens'' (a,a) a add' = \(x,y) -> (x + y, \ds -> (ds, ds)) dup' :: Num a => Lens'' a (a,a) dup' = \x -> ((x,x), \(dx,dy) -> dx + dy) sub' :: Num a => Lens'' (a,a) a sub' = \(x,y) -> (x - y, \ds -> (ds, -ds)) mul' :: Num a => Lens'' (a,a) a mul' = \(x,y) -> (x * y, \dz -> (dz * y, x * dz)) recip' :: Fractional a => Lens'' a a recip' = \x-> (recip x, \ds -> - ds / (x * x)) div' :: Fractional a => Lens'' (a,a) a div' = (\(x,y) -> (x / y, \d -> (d/y,-x*d/(y * y)))) sin' :: Floating a => Lens'' a a sin' = \x -> (sin x, \dx -> dx * (cos x)) cos' :: Floating a => Lens'' a a cos' = \x -> (cos x, \dx -> -dx * (sin x)) pow' :: Num a => Integer -> Lens'' a a pow' n = \x -> (x ^ n, \dx -> (fromInteger n) * dx * x ^ (n-1)) --cmul :: Num a => a -> Lens' a a --cmul c = lens (* c) (\x -> \dx -> c * dx) exp' :: Floating a => Lens'' a a exp' = \x -> let ex = exp x in (ex, \dx -> dx * ex) fst' :: Num b => Lens'' (a,b) a fst' = (\(a,b) -> (a, \ds -> (ds, 0))) snd' :: Num a => Lens'' (a,b) b snd' = (\(a,b) -> (b, \ds -> (0, ds))) -- some monoidal combinators swap' :: Lens'' (a,b) (b,a) swap' = (\(a,b) -> ((b,a), \(db,da) -> (da, db))) assoc' :: Lens'' ((a,b),c) (a,(b,c)) assoc' = \((a,b),c) -> ((a,(b,c)), \(da,(db,dc)) -> ((da,db),dc)) assoc'' :: Lens'' (a,(b,c)) ((a,b),c) assoc'' = \(a,(b,c)) -> (((a,b),c), \((da,db),dc)-> (da,(db,dc))) par' :: Lens'' a b -> Lens'' c d -> Lens'' (a,c) (b,d) par' l1 l2 = l3 where l3 (a,c) = let (b , j1) = l1 a in let (d, j2) = l2 c in ((b,d) , j1 *** j2) first' :: Lens'' a b -> Lens'' (a, c) (b, c) first' l = par' l id' second' :: Lens'' a b -> Lens'' (c, a) (c, b) second' l = par' id' l labsorb :: Lens'' ((),a) a labsorb (_,a) = (a, \a' -> ((),a')) labsorb' :: Lens'' a ((),a) labsorb' a = (((),a), \(_,a') -> a') rabsorb :: Lens'' (a,()) a rabsorb = comp labsorb swap' And here are the two dimensional combinators. I tried to write them point-free in terms of the combinators above to demonstrate that there is no monkey business going on. We type WAD' w w' a b = Lens'' (w,a) (w',b) type WAD'' w a b = WAD' w () a b -- terminate the weights for a closed network {- For any monoidal category we can construct this composition? -} -- horizontal composition hcompose :: forall w w' w'' w''' a b c. WAD' w' w'' b c -> WAD' w w''' a b -> WAD' (w',w) (w'',w''') a c hcompose f g = comp f' g' where f' :: Lens'' ((w',r),b) ((w'',r),c) f' = (first' swap') comp assoc'' comp (par' id' f) comp assoc' comp (first' swap') g' :: Lens'' ((r,w),a) ((r,w'''),b) g' = assoc'' comp (par' id' g) comp assoc' rotate :: WAD' w w' a b -> WAD' a b w w' rotate f = swap' comp f comp swap' -- vertical composition of weights vcompose :: WAD' w' w'' c d -> WAD' w w' a b -> WAD' w w'' (c, a) (d, b) vcompose f g = rotate (hcompose (rotate f) (rotate g) ) -- a double par. diagpar :: forall w w' a b w'' w''' c d. WAD' w w' a b -> WAD' w'' w''' c d -> WAD' (w,w'') (w',w''') (a, c) (b, d) diagpar f g = t' comp (par' f g) comp t where t :: Lens'' ((w,w''),(a,c)) ((w,a), (w'',c)) -- yikes. just rearrangements. t = assoc'' comp (second' ((second' swap') comp assoc' comp swap')) comp assoc' t' :: Lens'' ((w',b), (w''',d)) ((w',w'''),(b,d)) -- the tranpose of t t' = assoc'' comp (second' ( swap' comp assoc'' comp (second' swap'))) comp assoc' id''' :: WAD' () () a a id''' = id' -- rotate:: WAD' w a a w -- rotate = swap' liftIO :: Lens'' a b -> WAD' w w a b liftIO = second' liftW :: Lens'' w w' -> WAD' w w' a a liftW = first' wassoc' = liftW assoc' wassoc'' = liftW assoc'' labsorb'' :: WAD' ((),w) w a a labsorb'' = first' labsorb labsorb''' :: WAD' w ((),w) a a labsorb''' = first' labsorb' wswap' :: WAD' (w,w') (w',w) a a wswap' = first' swap' -- and so on we can lift all combinators I wonder if this is actually nice? I asked around and it seems like this idea may be what davidad is talking about when he refers to dioptics http://events.cs.bham.ac.uk/syco/strings3-syco5/slides/dalrymple.pdf Perhaps this will initiate a convo. Edit: He confirms that what I’m doing appears to be a dioptic. Also he gave a better link http://events.cs.bham.ac.uk/syco/strings3-syco5/papers/dalrymple.pdf He is up to some interesting diagrams ### Bits and Bobbles • Does this actually work or help make things any better? • Recurrent neural nets flip my intended role of weights and inputs. • Do conv-nets naturally require higher dimensional diagrams? • This weighty style seems like a good fit for my gauss seidel and iterative LQR solvers. A big problem I hit there was getting all the information to the outside, which is a similar issue to getting the weights around in a neural net. ## Functors, Vectors, and Quantum Circuits Vectors are dang useful things, and any energy you put into them seems to pay off massive dividends. Vectors and Linear Algebra are useful for: • 2D, 3D, 4D geometry stuff. Computer graphics, physics etc. • Least Squares Fitting • Solving discretized PDEs • Quantum Mechanics • Analysis of Linear Dynamical Systems • Probabilistic Transition Matrices There are certain analogies between Haskell Functors and Vectors that corresponds to a style of computational vector mathematics that I think is pretty cool and don’t see talked about much. Due to the expressivity of its type system, Haskell has a first class notion of container that many other languages don’t. In particular, I’m referring to the fact that Haskell has higher kinded types * -> * (types parametrized on other types) that you can refer to directly without filling them first. Examples in the standard library include Maybe, [], Identity, Const b, and Either b. Much more vector-y feeling examples can be found in Kmett’s linear package V0, V1, V2, V3, V4. For example, the 4 dimensional vector type V4 data V4 a = V4 a a a a This really isn’t such a strange, esoteric thing as it may appear. You wouldn’t blink an eye at the type struct V4 { double x, y, z, w; } in some other language. What makes Haskell special is how compositional and generic it is. We can build thousand element structs with ease via composition. What we have here is an alternative to the paradigm of computational vectors ~ arrays. Instead we have computational vectors ~ structs. In principle, I see no reason why this couldn’t be as fast as arrays, although with current compiler expectations it probably isn’t. Monoidal categories are a mathematical structure that models this analogy well. It has been designed by mathematicians for aesthetic elegance, and it seems plausible that following its example leads us to interesting, useful, and pleasant vector combinators. And I personally find something that tickles me about category theory. So to get started, let’s talk a bit about functors. ## The Algebra of Functors Functors in Haskell are a typeclass for containers. They allow you to map functions over all the items in the container. They are related to the categorical notion of functor, which is a mapping between categories. type Container = * -> * -- Note: This is actually a kind signature. -- Kinds and types are the same thing in Haskell. You can lift the product and sum of types to the product and sum of Functors which you may find in Data.Functor.Product and Data.Functor.Sum. This is analogous to the lifting of ordinary addition and multiplication to the addition and multiplication of polynomials, which are kind of like numbers with a “hole”. newtype Product f g a = Pair (f a , g a) newtype Sum f g a = Sum (Either (f a) (g a)) Functors also compose. A container of containers of a is still a container of a. We can form composite containers by using the Compose newtype wrapper. newtype Compose f g a = Compose (f (g a)) When you use this Compose newtype, instead of having to address the individual elements by using fmap twice, a single application of fmap will teleport you through both layers of the container. Product, Sum, and Compose are all binary operator on functors. The type constructor has kind {- Enter into ghci -} :kind Compose ==> (* -> *) -> (* -> *) -> (* -> *) :kind Product ==> (* -> *) -> (* -> *) -> (* -> *) :kind Sum ==> (* -> *) -> (* -> *) -> (* -> *) Some important other functors from the algebra of types perspective are Const Void a, Const () a, and Identity a. These are identity elements for Sum, Product, and Compose respectively. You can define mappings between containers that don’t depend on the specifics of their contents. These mappings can only rearrange, copy and forget items of their contained type. This can be enforced at the type level by the polymorphic type signature type f ~> g = forall a. f a -> g a These mappings correspond in categorical terminology to natural transformations between the functors f and g. There is a category where objects are Functors and morphisms are natural transformations. Sum, Product, and Compose all obeys the laws necessary to be a monoidal product on this category. How the lifting of functions works for Compose is kind of neat. mon_prod :: Functor f' => (f ~> f') -> (g ~> g') -> (Compose f g ~> Compose f' g') mon_prod ntf ntg (Compose fg) = Compose (fmap ntg (ntf fg)) -- or equvalently Compose (ntf (fmap ntg fg)) with a (Functor f) typeclass requirement. Because the natural transformations require polymorphic types, when you apply ntf to fg the polymorphic variable a in the type of ntf restricts to a ~ g a'. Product and Sum have a straight forward notion of commutativity ( (a,b) is isomorphic to (b,a)) . Compose is more subtle. sequenceA from the Traversable typeclass can swap the ordering of composition. sequenceA . sequenceA may or may not be the identity depending on the functors in question, so it has some flavor of a braiding operation. This is an interesting post on that topic https://parametricity.com/posts/2015-07-18-braids.html Combinators of these sorts are used arise in at least the following contexts • Data types a la carte – A systematic way of building extensible data types • GHC Generics – A system for building generic functions that operate on data types that can be described with sums, products, recursion, and holes. • In and around the Lens ecosystem Also see the interesting post by Russell O’Connor and functor oriented programming http://r6.ca/blog/20171010T001746Z.html. I think the above is part of that to which he is referring. ### Vector Spaces as Shape Vector spaces are made of two parts, the shape (dimension) of the vector space and the scalar. Just as a type of kind * -> * can be thought of as a container modulo it’s held type, it can also be a vector modulo its held scalar type. The higher kinded type for vector gives an explicit slot to place the scalar type. type Vect = * -> * The standard Haskell typeclass hierarchy gives you some of the natural operations on vectors if you so choose to abuse it in that way. • Functor ~> Scalar Multiplication: smul s = fmap (* s) • Applicative ~> Vector Addition: vadd x y = (+) <$> x <*> y
• Traversable ~> Tranposition. sequenceA has the type of transposition and works correctly for the linear style containers like V4.

The linear library does use Functor for scalar multiplication, but defines a special typeclass for addition, Additive. I think this is largely for the purposes for bringing Map like vectors into the fold, but I’m not sure.

Once we’ve got the basics down of addition and scalar multiplication, the next thing I want is operations for combining vector spaces. Two important ones are the Kronecker product and direct sum. In terms of indices, the Kronecker product is a space that is indexed by the cartesian product (,) of its input space indices and the direct sum is a space indexed by the Either of its input space indices. Both are very useful constructs. I use the Kronecker product all the time when I want to work on 2D or 3D grids for example. If you’ll excuse my python, here is a toy 2-D finite difference Laplace equation example. We can lift the 1D second derivative matrix $K = \partial_x^2$ using the kronecker product $K2 = K \otimes I + I \otimes K$. The direct sum is useful as a notion of stacking matrices.

import numpy as np
N = 10 # We're making a 10x10 grid
row = np.zeros(N)
row[0] = -2
row[1] = 1
K = np.toeplitz(row,row) #toeplitz makes constant diagonal matrices
I = np.eye(N) #identity matrix
K2 = np.kron(K,I) + np.kron(I,K)

The following is perhaps the most important point of the entire post.

type Kron = Compose
type DSum = Product

Compose of vector functors gives the Kronecker product, and Product gives the direct sum (this can be confusing but its right. Remember, the sum in direct sum refers to the indices).

We can form the Kronecker product of vectors given a Functor constraint.

kron :: (Num a, Functor f, Functor g) => f a -> g a -> Kron f g a
liftI2 f (Compose x) (Compose y) = Compose $liftI2 (liftI2 f) x y ## Monoidal Categories The above analogy can be put into mathematical terms by noting that both vectors and functor are monoidal categories. I talked a quite a bit about monoidal categories in a previous post http://www.philipzucker.com/a-touch-of-topological-computation-3-categorical-interlude/ . Categories are the combo of a collection of objects and arrows between the objects. The arrows can compose as long as the head of one is on the same object as the tail of the other. On every object, there is always an identity arrow, which when composed will do nothing. We need a little extra spice to turn categories into monoidal categories. One way of thinking about it is that monoidal categories have ordinary category composition and some kind of horizontal composition, putting things side to side. Ordinary composition is often doing something kind of sequentially, applying a sequence of functions, or a sequence of matrices. The horizontal composition is often something parallel feeling, somehow applying the two arrows separately to separate pieces of the system. ### Why are they called Monoidal? There is funny game category people play where they want to lift ideas from other fields and replace the bits and pieces in such a way that the entire thing is defined in terms of categorical terminology. This is one such example. A monoid is a binary operations that is associative and has an identity. Sometimes people are more familiar with the concept of a group. If not, ignore the next sentence. Monoids are like groups without requiring an inverse. Numbers are seperately monoids under both addition, multiplication and minimization (and more), all of which are associative operations with identity (0, 1, and infinity respectively). Exponentiation is a binary operation that is not a monoid, as it isn’t associative. The Monoid typeclass in Haskell demonstrates this http://hackage.haskell.org/package/base-4.12.0.0/docs/Data-Monoid.html class Semigroup a => Monoid a where mempty :: a mappend :: a -> a -> a A common example of a monoid is list, where mempty is the empty list and mappend appends the lists. There are different set-like intuitions for categories. One is that the objects in the category are big opaque sets. This is the case for Hask, Rel and Vect. A different intuitiion is that the category itself is like a set, and the objects are the elements of that set. There just so happens to be some extra structure knocking around in there: the morphisms. This is the often more the feel for the examples of preorders or graphs. The word “monoidal” means that they we a binary operation on the objects. But in the category theory aesthetic, you also need that binary operation to “play nice” with the morphisms that are hanging around too. Functors are the first thing that has something like this. It has other properties that come along for the ride. A Functor is a map that takes objects to objects and arrows to arrows in a nice way. A binary functor takes two objects to and object, and two arrows to one arrow in a way that plays nice (commutes) with arrow composition. ### String diagrams String diagrams are a graphical notation for monoidal categories. Agin I discussed this more here. Morphisms are denoted by boxes. Regular composition is shown by plugging arrows together vertically. Monoidal product is denoted by putting the arrows side to side. When I was even trying to describe what a monoidal category was, I was already using language evocative of string diagrams. You can see string diagrams in the documentation for the Arrow library. Many diagrams that people use in various fields can be formalized as the string diagrams for some monoidal category. This is big chunk of Applied Category Theory. This is the connection to quantum circuits, which are after all a graphical notation for very Kroneckery linear operations. type Qubit = V2 type C = Complex Double assoc :: Functor f => (Kron (Kron f g) h) ~> (Kron f (Kron g h)) assoc = Compose . (fmap Compose) . getCompose . getCompose assoc' :: Functor f => (Kron f (Kron g h)) ~> (Kron (Kron f g) h) assoc' (Compose x) = Compose$ Compose $(fmap getCompose x) kron'' :: (Additive f, Additive g, Additive k, Additive f', Additive g') => LinOp1 f f' a -> LinOp1 g g' a -> Kron (Kron f g) k a -> Kron (Kron f' g') k a kron'' lf lg fgk = let v = (assoc fgk) in assoc' (Compose$ fmap lg $getCompose (lf v)) sigx' :: LinOp1 Qubit Qubit C sigx' (Compose (V2 up down)) = Compose$ V2 down up

sigz' :: LinOp1 Qubit Qubit C
sigz' (Compose (V2 up down)) = Compose  $V2 up ((-1) *^ down) sigy' :: LinOp1 Qubit Qubit C sigy' (Compose (V2 up down)) = Compose$ V2 ((-i) *^ down) (i *^ up) where i = 0 :+ 1

swap' :: (Traversable f, Applicative g) => LinOp1 (Kron f g) (Kron g f) a
swap' (Compose (Compose fgk)) = Compose $Compose$ sequenceA fgk

cnot :: LinOp1 (Kron Qubit Qubit) (Kron Qubit Qubit) a
cnot (Compose (Compose (V2 (V2 up1 down1) v))) = Compose $Compose$ V2 (V2 down1 up1) v

phase :: Double -> LinOp1 Qubit Qubit C
phase phi (Compose (V2 up down)) = Compose $V2 up ((cis phi) *^ down) lefting :: (Additive f, Additive k, Additive g) => LinOp1 f g a -> LinOp1 (Kron f k) (Kron g k) a lefting l = kron'' l id -- Qubit.assoc' . l . Qubit.assoc righting :: (Additive k, Additive f, Additive g) => LinOp1 f g a -> LinOp1 (Kron k f) (Kron k g) a righting l = kron'' id l -- (Compose (Compose fkk)) = Compose$ Compose $(fmap (getCompose . l . Compose) fkk) example :: LinOp1 (Kron Qubit Qubit) (Kron Qubit Qubit) C example = (lefting sigx') . (lefting sigy') . (righting sigz') . swap'  There is an annoying amount of stupid repetitive book keeping with the associative structure of Kron. This can largely be avoided hopefully with coerce, but I’m not sure. I was having trouble with roles when doing it generically. ### Bit and Bobbles • Woof. This post was more draining to write than I expected. I think there is still a lot left to say. Sorry about the editing everyone! Bits and pieces of this post are scattered in this repo • How would you go about this in other languages? C, Rust, OCaml, C++, Agda • The discussion of Vect = * -> * is useful for discussion of 2-Vect, coming up next. What if we make vectors of Vect? Wacky shit. • Metrics and Duals vectors. type Dual f a = f a -> a. type Dual1 f a = forall k. Additive k => Kron f k a -> k a • Adjunction diagrams have cups and caps. Since we have been using representable functors, they actually have a right adjunction that is tupling with the vector space index type. This gives us something that almost feels like a metric but a weirdly constrained metric. • LinOp1 form is yoneda? CPS? Universally quantified k is evocative of forall c. (a -> c) -> (b -> c) ### References ## Appendix ### Representable/Naperian Functors Containers that are basically big product types are also known as representable, Naperian, or logarithmic. Representable places emphasis on the isomorphism between such a container type and the type (->) i which by the algebra of types correspond is isomorphic to $a^i$ (i copies of a). They are called Naperian/Logarithmic because there is a relationship similar to exponentiation between the index type a and the container type f. If you take the Product f g, this container is indexed by (a + b) = Either a b. Compose f g is indexed by the product (a,b). (f r) ~ r^a The arrow type is written as an exponential b^a because if you have finite enumerable types a and b, that is the number of possible tabulations available for f. The Sum of two representable functors is no longer representable. Regular logarithms of sums Log(f + g) do not have good identities associated with them. See Gibbons article. There is a good argument to be made that representable functors are a good match for vectors/well typed tensor programming. But note that there is a reasonable interpretation for container types with sum types in them. These can be thought of as subspaces, different bases, or as choices of sparsity patterns. When you define addition, you’ll need to say how these subspaces reconcile with each other. — two bases at 45 degrees to each other. data V2_45 a = XY a a | XY' a a data Maybe a = Just a | Notinhg -- a 1-d vectro space with a special marker for the zero vector. data Maybe2 a = Just2 a a | Nothing2 -- a 2d vector space with special zero marker ### Monoidal Products on Hask Hask is a name for the category that has objects as Haskell types and morphisms as Haskell functions. Note that it’s a curious mixing of type/value layers of Haskell. The objects are types whereas the function morphisms are Haskell values. Composition is given by (.) and the identity morphisms are given by id. For Haskell, you can compose functions, but you can also smash functions together side by side. These combinators are held in Control.Arrow. You can smash together types with tuple (,) or with Either. Both of these are binary operators on types. The corresponding mapping on morphisms are given by (***) :: a b c -> a b' c' -> a (b, b') (c, c') (+++) :: a b c -> a b' c' -> a (Either b b') (Either c c')  These are binary operators on morphisms that play nice with the composition structure of Haskell. ### Monoidal Combinators of Functors A monoidal category also has unit objects. This is given by the Identity functor rightUnitor :: Functor f => Compose f Identity a -> f a rightUnitor (Compose f) = fmap runIdentity f rightUnitor' :: f ~> Compose f Identity rightUnitor' = Compose . fmap Identity leftUnitor' :: f ~> Compose Identity f leftUnitor' = Compose . Identity leftUnitor :: Identity *** f ~> f leftUnitor = runIdentity . getCompose There is also a sense of associativity. It is just newtype rearrangement, so it can also be achieved with a coerce (although not polymorphically?). assoc :: Functor f => ((f *** g) *** h) ~> (f *** (g *** h)) assoc = Compose . (fmap Compose) . getCompose . getCompose assoc' :: Functor f => (f *** (g *** h)) ~> ((f *** g) *** h) assoc' (Compose x) = Compose$ Compose $(fmap getCompose x) Similarly, we can define a monoidal category structure using Product or Sum instead of Compose. These are all actually just newtype rearrangement, so they should all just be instances of coerce, but I couldn’t get the roles to go through generically? ## Concolic Weakest Precondition is Kind of Like a Lens That’s a mouthful. Lens are described as functional getters and setters. The simple lens type is type Lens a b = a -> (b, b -> a) . The setter is a->b and the getter is a -> b -> a This type does not constrain lenses to obey the usual laws of getters and setters. So we can use/abuse lens structures for nontrivial computations that have forward and backwards passes that share information. Jules Hedges is particular seems to be a proponent for this idea. I’ve described before how to encode reverse mode automatic differentiation in this style. I have suspicions that you can make iterative LQR and guass-seidel iteration have this flavor too, but I’m not super sure. My attempts ended somewhat unsatisfactorily a whiles back but I think it’s not hopeless. The trouble was that you usually want the whole vector back, not just its ends. I’ve got another example in imperative program analysis that kind of makes sense and might be useful though. Toy repo here: https://github.com/philzook58/wp-lens In program analysis it sometimes helps to run a program both concretely and symbolically. Concolic = CONCrete / symbOLIC. Symbolic stuff can slowly find hard things and concrete execution just sprays super fast and can find the dumb things really quick. We can use a lens structure to organize a DSL for describing a simple imperative language The forward pass is for the concrete execution. The backward pass is for transforming the post condition to a pre condition in a weakest precondition analysis. Weakest precondition semantics is a way of specifying what is occurring in an imperative language. It tells how each statement transforms post conditions (predicates about the state after the execution) into pre conditions (predicates about before the execution). The concrete execution helps unroll loops and avoid branching if-then-else behavior that would make the symbolic stuff harder to process. I’ve been flipping through Djikstra’s book on this. Interesting stuff, interesting man. I often think of a state machine as a function taking s -> s. However, this is kind of restrictive. It is possible to have heterogenous transformations s -> s’. Why not? I think I am often thinking about finite state machines, which we really don’t intend to have a changing state size. Perhaps we allocated new memory or something or brought something into or out of scope. We could model this by assuming the memory was always there, but it seems wasteful and perhaps confusing. We need to a priori know everything we will need, which seems like it might break compositionally. We could model our language making some data type like data Imp = Skip | Print String | Assign String Expr | Seq Imp Imp | ... and then build an interpreter interp :: Imp -> s -> s' But we can also cut out the middle man and directly define our language using combinators. type Stmt s s' = s ->s' To me this has some flavor of a finally tagless style. Likewise for expressions. Expressions evaluate to something in the context of the state (they can lookup variables), so let’s just use type Expr s a = s -> a And, confusingly (sorry), I think it makes sense to use Lens in their original getter/setter intent for variables. So Lens structure is playing double duty. type Var s a = Lens' s a With that said, here we go.  type Stmt s s' = s -> s' type Lens' a b = a -> (b, b -> a) set l s a = let (_, f) = l s in f a type Expr s a = s -> a type Var s a = Lens' s a skip :: Stmt s s skip = id sequence :: Stmt s s' -> Stmt s' s'' -> Stmt s s'' sequence = flip (.) assign :: Var s a -> Expr s a -> Stmt s s assign v e = \s -> set v s (e s) (===) :: Var s a -> Expr s a -> Stmt s s v === e = assign v e ite :: Expr s Bool -> Stmt s s' -> Stmt s s' -> Stmt s s' ite e stmt1 stmt2 = \s -> if (e s) then stmt1 s else stmt2 s while :: Expr s Bool -> Stmt s s -> Stmt s s while e stmt = \s -> if (e s) then ((while e stmt) (stmt s)) else s assert :: Expr s Bool -> Stmt s s assert e = \s -> if (e s) then s else undefined abort :: Stmt s s' abort = const undefined  Weakest precondition can be done similarly, instead we start from the end and work backwards Predicates are roughly sets. A simple type for sets is type Pred s = s -> Bool  Now, this doesn’t have much deductive power, but I think it demonstrates the principles simply. We could replace Pred with perhaps an SMT solver expression, or some data type for predicates, for which we’ll need to implement things like substitution. Let’s not today. A function a -> b  is equivalent to forall c. (b -> c) -> (a -> c) . This is some kind of CPS / Yoneda transformation thing. A state transformer s -> s' to predicate transformer (s' -> Bool) -> (s -> Bool) is somewhat evocative of that. I’m not being very precise here at all. Without further ado, here’s how I think a weakest precondition looks roughly.  type Lens' a b = a -> (b, b -> a) set l s a = let (_, f) = l s in f a type Expr s a = s -> a type Var s a = Lens' s a type Pred s = s -> Bool type Stmt s s' = Pred s' -> Pred s skip :: Stmt s s skip = \post -> let pre = post in pre -- if sequence :: Stmt s s' -> Stmt s' s'' -> Stmt s s'' sequence = (.) assign :: Var s a -> Expr s a -> Stmt s s assign v e = \post -> let pre s = post (set v s (e s)) in pre (===) :: Var s a -> Expr s a -> Stmt s s v === e = assign v e ite :: Expr s Bool -> Stmt s s' -> Stmt s s' -> Stmt s s' ite e stmt1 stmt2 = \post -> let pre s = if (e s) then (stmt1 post) s else (stmt2 post) s in pre abort :: Stmt s s' abort = \post -> const False assert :: Expr s Bool -> Stmt s s assert e = \post -> let pre s = (e s) && (post s) in pre {- -- tougher. Needs loop invariant while :: Expr s Bool -> Stmt s s -> Stmt s s while e stmt = \post -> let pre s = if (e s) then ((while e stmt) (stmt post)) s else in pre -}  Finally here is a combination of the two above that uses the branching structure of the concrete execution to aid construction of the precondition. Although I haven’t expanded it out, we are using the full s t a b parametrization of lens in the sense that states go forward and predicates come back.  type Lens' a b = a -> (b, b -> a) set l s a = let (_, f) = l s in f a type Expr s a = s -> a type Var s a = Lens' s a type Pred a = a -> Bool type Stmt s s' = s -> (s', Pred s' -> Pred s) -- eh. Screw the newtype skip :: Stmt s s skip = \x -> (x, id) sequence :: Stmt s s' -> Stmt s' s'' -> Stmt s s'' sequence f g = \s -> let (s', j) = f s in let (s'', j') = g s' in (s'', j . j') assign :: Var s a -> Expr s a -> Stmt s s assign v e = \s -> (set v s (e s), \p -> \s -> p (set v s (e s))) --if then else ite :: Expr s Bool -> Stmt s s' -> Stmt s s' -> Stmt s s' ite e stmt1 stmt2 = \s -> if (e s) then let (s', wp) = stmt1 s in (s', \post -> \s -> (e s) && (wp post s)) else let (s', wp) = stmt2 s in (s', \post -> \s -> (not (e s)) && (wp post s)) assert :: Pred s -> Stmt s s assert p = \s -> (s, \post -> let pre s = (post s) && (p s) in pre) while :: Expr s Bool -> Stmt s s -> Stmt s s while e stmt = \s -> if e s then let (s' , wp) = (while e stmt) s in (s', \post -> let pre s'' = (post s'') && (wp post s'') in pre) else (s, \p -> p) {- -- declare and forget can change the size and shape of the state space. -- These are heterogenous state commpands declare :: Iso (s,Int) s' -> Int -> Stmt s s' declare iso defalt = (\s -> to iso (s, defalt), \p -> \s -> p$ to iso (s, defalt))

forget :: Lens' s s' -> Stmt s s' -- forgets a chunk of state

declare_bracket :: Iso (s,Int) s' -> Int ->  Stmt s' s' -> Stmt s s
declare_bracket iso defalt stmt = (declare iso default) . stmt . (forget (_1 . iso))


Neat. Useful? Me dunno.

## Linear Algebra of Types

It gives my brain a pleasant thrum to learn new mathematics which mimics the algebra I learned in middle school. Basically this means that the new system has operations with properties that match those of regular numbers as much as possible. Two pretty important operations are addition and multiplication with the properties of distributivity and associativity. Roughly this corresponds to the mathematical notion of a Semiring.

Some examples of semirings include

• And-Or
• Min-plus
• Matrices.
• Types

I have written before about how types also form a semiring, using Either for plus and (,) for times. These constructions don’t obey distributivity or associativity “on the nose”, but instead are isomorphic to the rearranged type, which when you squint is pretty similar to equality.

Matrices are grids of numbers which multiply by “row times column”. You can form matrices out of other semirings besides just numbers. One somewhat trivial but interesting example is block matrices, where the elements of the matrix itself are also matrices. Another interesting example is that of relations, which can be thought of as matrices of boolean values. Matrix multiplication using the And-Or semiring on the elements corresponds to relational composition.

What if we put our type peanut butter in our matrix chocolate and consider matrices of types, using the Either(,) semiring?

The simplest implementation to show how this could go can be made using the naive list based implementation of vectors and matrices. We can directly lift this representation to the typelevel and the appropriate value-level functions to type families.

type a :*: b = (a,b)
type a :+: b = Either a b

type family Dot v v' where
Dot '[x] '[y] = x :*: y
Dot (x : xs) (y : ys) = (x :*: y) :+: (Dot xs ys)

type family MVMult m v where
MVMult '[r] v = '[Dot r v]
MVMult (r : rs) v = (Dot r v) : (MVMult rs v)

type family VMMult m v where
VMMult v ' = '[Dot v c]
VMMult v (c : cs) = (Dot v c) : (VMMult v cs)

type family MMMult' m m' where
MMMult' '[r] cs = '[VMMult r cs]
MMMult' (r : rs) cs = (VMMult r cs) : (MMMult' rs cs)

type family MMMult m m' where
MMMult m m' = MMMult' m (Transpose m')

type family Transpose m where
Transpose ((r1 : rs') : rs) = (r1 : (Heads rs)) : (Conss rs' (Transpose (Tails rs)))
Transpose '[] = '[]

-- some mapped helper functions
-- verrrrrry ugly. Eh. Get 'er dun
type family Tails v where
Tails ((v : vs) : xs) = vs : (Tails xs)
Tails '[] = '[]
type family Conss v vs where
Conss (x : xs) (y : ys) = (x : y) : (Conss xs ys)
Conss '[] '[] = '[]

type family Index v i where
Index (x : xs) 0 = x
Index (x : xs) n = Index xs (n-1)


This was just for demonstration purposes. It is not my favorite representation of vectors. You can lift a large fraction of possible ways to encode vector spaces at the value level up to the type level, such as the linear package, or using dual vectors type V2 a = a -> a -> a. Perhaps more on that another day.

### What is the point?

Ok. That’s kind of neat, but why do it? Well, one way to seek an answer to that question is to ask “what are matrices useful for anyway?”

One thing they can do is describe transition systems. You can write down a matrix whose entire $a_{ij}$ describes something about the transition from state $i$ to state $j$. For example the entry could be:

• The cost of getting from $i$ to $j$ (min-plus gives shortest path),
• The count of ways to get from $i$ to $j$ (combinatorics of paths)
• The connectivity of the system from $i$ to $j$ using boolean values and the and-or semiring
• The probability of transition from $i$ to $j$
• The quantum amplitude of going from $i$ to $j$ if we’re feeling saucy.

If we form a matrix describing a single time step, then multiplying this matrix by itself gives 2 time steps and so on.

Lifting this notion to types, we can build a type exactly representing all the possible paths from state $i$ to $j$.

Concretely, consider the following humorously bleak transition system: You are going between home and work. Every 1 hour period you can make a choice to do a home activity, commute, or work. There are different options of activities at each.

data Commute = Drive
data Home = Sleep | Eat
data Work = TPSReport | Bitch | Moan

This is described by the following transition diagram

The transitions are described by the following matrix.type:

type T = '[ '[Home    ,  Commute ],
'[Commute ,  Work    ]]


What is the data type that describe all possible 4-hour day? You’ll find the appropriate data types in the following matrix.

type FourHour = MMMult T (MMMult T (MMMult T T))

Now, time to come clean. I don’t think this is necessarily the best way to go about this problem. There are alternative ways of representing it.

Here are two data types that describe an indefinite numbers of transition steps.

data HomeChoice = StayHome Home HomeChoice | GoWork Commute WorkChoice
data WorkChoice = StayWork Work WorkChoice | GoHome Commute HomeChoice

Another style would hold the current state as a type parameter in the type using a GADT.

data Path state where
StayWork :: Work -> Path Work -> Path Work
CommuteHome :: Commute -> Path Home ->  Path Work
StayHome :: Home -> Path Home -> Path Home
CommuteWork :: Commute -> Path Work ->  Path Home

We could construct types that are to the above types as Vec n  is to [] by including an explicit step size parameter.

Still, food for thought.

### Further Thoughts

The reason i was even thinking about this is because we can lift the above construction to perform a linear algebra of vectors spaces. And I mean the spaces, not the vectors themselves. This is a confusing point.

Vector spaces have also have two natural operations on them that act like addition and multiplication, the direct sum and kronecker product. These operations do form a semiring, although again not on the nose.

This is connected to the above algebra of types picture by considering the index types of these vector spaces. The simplest way to denote this in Haskell is using the free vector space construction as shown in this Dan Piponi post. The Kronecker product makes tuples of the indices and the direct sum has an index that is the Either of the original index types.

type Vec b r = [(b, a)]
-- Example 2D vector space type
type V2D = Vec Bool Double

This is by far not the only way to go about it. We can also consider using the Compose-Product semiring on functors (Compose is Kron, Product is DSum) to get a more index-free kind of feel and work with dense vectors.

Going down this road (plus a couple layers of mathematical sophistication) leads to a set of concepts known as 2Vect. Dan Roberts and James Vicary produced a Mathematica package for 2Vect which is basically incomprehensible to me. It seems to me that typed functional programming is a more appropriate venue for something of this kind of pursuit, given how evocative/ well modeled by category theory it can be. These mathematical ideas are applicable to describing anyonic vector spaces. See my previous post below. It is not a coincidence that the Path data type above is so similar to FibTree data type. The root type variable takes the place of the work/home state, and the tuple structure take the place of a Vec-like size parameter n .

More to on this to come probably as I figure out how to explain it cleanly.

Edit: WordPress, your weird formatting is killing me.

Edit: Hoo Boy. This is why we write blog posts. Some relevant material was pointed out to me that I was not aware of. Thanks @DrEigenbastard.

http://blog.sigfpe.com/2010/08/constraining-types-with-regular.html

http://blog.sigfpe.com/2010/08/divided-differences-and-tomography-of.html

## Relational Algebra with Fancy Types

Last time, I tried to give a primer of relations and relational algebra using the Haskell type type Rel a b = [(a,b)]. In this post we’re going to look at these ideas from a slightly different angle. Instead of encoding relations using value level sets, we’ll encode relations in the type system. The Algebra of Programming Agda repo and the papers quoted therein are very relevant, so if you’re comfortable wading into those waters, give them a look. You can find my repo for fiddling here

At this point, depending on what you’ve seen before, you’re either thinking “Yeah, sure. That’s a thing.” or you’re thinking “How and why the hell would you do such a ridiculous thing.”

Most of this post will be about how, so let’s address why first:

1. Examining relations in this style illuminates some constructions that appear around the Haskell ecosystem, particularly some peculiar fellows in the profunctor package.
2. More syntactic approach to relations allows discussion of larger/infinite domains. The finite enumerations of the previous post is nice for simplicity, but it seems you can’t get far that way.
3. Mostly because we can – It’s a fun game. Maybe a useful one? TBD.

With that out of the way, let’s go on to how.

### Translating Functions to Relation GADTs

We will be using some Haskell extensions in this post, at the very least GADTs and DataKinds. For an introduction to GADTs and DataKinds, check out this blog post. DataKinds is an extension that reflects every data constructor of data types to a type constructor. Because there are values True and False there are corresponding types created'True and 'False. GADTs is an extension of the type definition mechanism of standard Haskell. They allow you to declare refined types for the constructors of your data and they infer those refined type when you pattern match out of the data as well, such that the whole process is kind of information preserving.

We will use the GADT extension to define relational datatypes with the kind

a -> b -> *
. That way it has a slot a for the “input” and b for the “output” of the relation. What will goes in these type slots will be DataKind lifted types like 'True, not ordinary Haskell types like Int. This is a divergence from from the uses of similar kinds you see in Category, Profunctor, or Arrow. We’re doing a more typelevel flavored thing than you’ll see in those libraries. What we’re doing is clearly a close brother of the singleton approach to dependently typed programming in Haskell.

Some examples are in order for what I mean. Here are two simple boolean functions, not and and defined in ordinary Haskell functions, and their equivalent GADT relation data type.

not True = False
not False = True

data Not a b where
NotTF :: Not 'True 'False
NotFT :: Not 'False 'True

and True True = True
and False _ = False
and _ False = False

data And a b where
AndTT :: And '( 'True, 'True) 'True
AndFU :: And '( 'False, a) 'False
AndUF :: And '( a, 'False) 'False

You can already start to see how mechanical the correspondence between the ordinary function definition and our new fancy relation type. A function is often defined via cases. Each case corresponds to a new constructor of the relation and each pattern that occurs in that case is the pattern that appears in the GADT. Multiple arguments to the relations are encoded by uncurrying everything by default.

Any function calls that occur on the right hand side of a function definition becomes fields in the constructor of our relation. This includes recursive calls and external function calls. Here are some examples with a Peano style natural number data type.

data Nat = S Nat | Z

plus Z x = x
plus (S x) y = S (plus x y)

data Plus a b where
PZ :: Plus '( 'Z, a) a
PS :: Plus '( a,b) c -> Plus '( 'S a, b) c


We can also define things that aren’t functions. Relations are a larger class of things than functions are, which is part of their utility. Here is a “less than equal” relation LTE.

data LTE a b where
LTERefl :: LTE n n
LTESucc :: LTE n m -> LTE n ('S m)

You can show that elements are in a particular relation by finding a value of that relational type. Is ([4,7], 11) in the relation Plus? Yes, and I can show it with with the value PS (PS (PS (PS PZ))) :: Plus (4,7) 11 . This is very much the Curry-Howard correspondence. The type R a b corresponds to the proposition/question is $(a,b) \in R$ .

### The Fun Stuff : Relational Combinators

While you need to build some primitive relations using new data type definitions, others can be built using relational combinators. If you avoid defining too many primitive relations like the above and build them out of combinators, you expose a rich high level manipulation algebra. Otherwise you are stuck in the pattern matching dreck. We are traveling down the same road we did in the previous post, so look there for less confusing explanations of the relational underpinnings of these constructions, or better yet some of the references below.

Higher order relational operators take in a type parameters of kind

a -> b -> *
and produce new types of a similar kind. The types appearing in these combinators is the AST of our relational algebra language.

The first two combinators of interest is the composition operator and the identity relation. An element $(a,c)$ is in $R \cdot Q$ if there exists a $b$ such that $(a,b) \in R$ and $(b,c) \in Q$. The fairly direct translation of this into a type is

{- rcompose :: Rel b c -> Rel a b -> Rel a c  -}

data RCompose k1 k2 a c where
RCompose :: k1 b c -> k2 a b -> RCompose k1 k2 a c

type k <<< k' = RCompose k k'
type k >>> k' = RCompose k' k

The type of the composition is the same as that of Profunctor composition found in the profunctors package.

type RCompose = Procompose

Alongside a composition operator, it is a knee jerk to look for an identity relation and we do have one

data Id a b where
Refl :: Id a a

-- monomorphic identity. Leave this out?
data IdBool a b where
ReflTrue :: IdBool 'True 'True
ReflFalse :: IdBool 'False 'False

This is also a familiar friend. The identity relation in this language is the Equality type.

-- identity function is the same as Equality
type Id a b = Id (a :~: b)

We can build an algebra for handling product and sum types by defining the appropriate relational combinators. These are very similar to the combinators in the Control.Arrow package.

-- Product types

data Fan k k' a b where
Fan :: k a b -> k' a c -> Fan k k' a '(b,c)

type k &&& k' = Fan k k'

data Fst a b where
Prj1 :: Fst '(a, b) a

data Snd a b where
Prj2 :: Snd '(a, b) b

-- Sum type

data Split k k' a b where
CaseLeft :: k a c -> Split k k' ('Left a) c
CaseRight :: k' b c -> Split k k' ('Right b) c

type k ||| k' = Split k k'

data Inj1 a b where
Inj1 :: Inj1 a ('Left a)
data Inj2 a b where
Inj2 :: Inj2 a ('Right a)

-- some derived combinators
type Par f g = Fan (f <<< Fst) (g <<< Snd)
type Dup  = Fan Id Id
type Swap = Fan Snd Fst


The converse of relations is very interesting operation and is the point where relations really differ from functions. Inverting a function is tough. Conversing a relation always works. This data type has no analog in profunctor to my knowledge and probably shouldn't.

data RConverse k a b where
RConverse :: k a b -> RConverse k b a
-- Shorter synonym
type RCon = RConverse

Relations do not have a notion of currying. The closest thing they have is

data Trans k a b where
Trans :: k '(a,b) c -> Trans k a '(b,c)

### Lattice Operators

For my purposes, lattices are descriptions of sets that trade away descriptive power for efficiency. So most operations you'd perform on sets have an analog in the lattice structure, but it isn't a perfect matching and you're forced into approximation. It is nice to have the way you perform these approximation be principled, so that you can know at the end of your analysis whether you've actually really shown anything or not about the actual sets in question.

The top relation holds all values. This is represented by making no conditions on the type parameters. They are completely phantom.

newtype Top a b = Top ()

Bottom is a relation with no inhabitants.

newtype Bottom a b = Bottom Void

The meet is basically the intersection of the relations, the join is basically the union.

newtype RMeet k k' a b = RMeet (k a b, k' a b)
type k /\ k' = RMeet k k'

newtype RJoin k k' a b = RJoin (Either (k a b) (k' a b))
type k \/ k' = RJoin k k' 

A Lattice has an order on it. This order is given by relational inclusion. This is the same as the :-> combinator can be found in the profunctors package.

type (:->) p q = forall a b. p a b -> q a b
type RSub p q = p :-> q

Relational equality can be written as back and forth inclusion, a natural isomorphism between the relations. There is also an interesting indirect form.

data REq k k' = REq {to' :: k :-> k', from' :: k' :-> k }

#### Relational Division

If we consider the equation (r <<< p) :-> q with p and q given, in what sense is there a solution for r? By analogy, this looks rather like r*p = q, so we're asking a kind of division question. Well, unfortunately, this equation may not necessarily have a solution (neither do linear algebraic equations for that matter), but we can ask for the best under approximation instead. This is the operation of relational division. It also appears in the profunctor package as the right Kan Extension. You'll also find the universal property of the right division under the name curryRan and uncurryRan in that module.

newtype Ran p q a b = Ran { runRan :: forall x. p x a -> q x b }
type RDiv = Ran

One formulation of Galois connections can be found in the adjunctions file. Galois Connections are very slick, but I'm running out of steam, so let's leave that one for another day.

### Properties and Proofs

We can prove many properties about these relational operations. Here a a random smattering that we showed using quickcheck last time.

prop_ridleft ::  (k <<< Id) :-> k
prop_ridleft (RCompose k IdRefl) = k

prop_ridright ::  (Id <<< k) :-> k
prop_ridright (RCompose IdRefl k) = k

prop_meet :: p /\ q :-> p
prop_meet (RMeet (p, q)) = p

prop_join :: p :-> p \/ q
prop_join p = RJoin (Left p)

meet_assoc :: RMeet k (RMeet k' k'') a b -> RMeet (RMeet k k') k'' a b
meet_assoc (RMeet (k, (RMeet (k',k'')))) = RMeet (RMeet (k,k'), k'')

prop_top :: k :-> Top
prop_top _ = top

prop_bottom :: Bottom :-> k
prop_bottom (Bottom x) = absurd x

bottom_compose :: REq (k <<< Bottom) Bottom
bottom_compose = REq (\(RCompose k (Bottom b)) -> absurd b) prop_bottom

data Iso a b = Iso {to :: a -> b, from :: b -> a}
type a <-> b = Iso a b

meet_universal :: (p ::-> RMeet k k') <-> (p ::-> k, p ::-> k')
meet_universal = Iso to from where
to (RSub f) = (RSub $\p -> case f p of RMeet (k,k') -> k , RSub$ \p -> case f p of RMeet (k,k') -> k')
from (RSub f,RSub g) = RSub $\p -> RMeet (f p, g p) prop_con :: RCon (RCon k) :-> k prop_con (RConverse (RConverse k)) = k ### Odds and Ends • Recursion Schemes - Recursion schemes are a methodology to talk about recursion in a point free style and where the rubber meets the road in the algebra of programming. Here is an excellent series of articles about them. Here is a sample of how I think they go: data MapMaybe k a b where MapJust :: k a b -> MapMaybe k ('Just a) ('Just b) MapNothing :: MapMaybe k 'Nothing 'Nothing data Cata map k a b where Cata :: k fa a -> map (Cata map k) x fa -> Cata map k ('Fix x)  • Higher Order Relations? • Examples of use. Check out the examples folder in the AoP Agda repo. These are probably translatable into Haskell. • Interfacing with Singletons. Singletonized functions are a specialized case or relations. Something like? •  newtype SFun a b = SFun (Sing a -> Sing b)  • A comment to help avoid confusion. What we've done here feels confusingly similar to profunctor, but it is in fact distinct I think. Profunctors are described as a categorical generalization of relations , but to be honest, I kind of don't get it. Despite many of our constructions appearing in the profunctor package, the profunctor typeclass itself appears to not play a role in our formulation. There just isn't a good way to dimap under our relations as written, unless you construct free profunctors. Converse at the least is a wrench in the works. • Star and graphs. Transition relations are a powerful methodology. A transition relation is in some respects the analog of a square matrix. We can iteratively compose it with itself. -- Check out "term rewriting and all that" -- This is also the reflection without remorse data type -- TSequence http://okmij.org/ftp/Haskell/zseq.pdf -- this is also a free instance of Category data Star k a b where Done :: Star k a a Roll :: k b c -> Star k a b -> Star k a c data KPlus k a b where PDone :: k a b -> KPlus k a b PRoll :: k b c -> KPlus k a b -> KPlus k a c type SymClos k a b = RJoin k (RCon k) a b type RefClos k a b = RJoin k Id a b {- n-fold composition -} -- similar to Fin. -- This is also the Vec n is to list and this is to reflection without remorse. Kind of interesting data NFold n k a b where One :: k a b -> NFold ('S n) k a b More :: k b c -> NFold n k a b -> NFold ('S n) k a b ### References ## Doing Basic Ass Shit in Haskell Haskell has lots of fancy weird corners, but you want to get rippin’ and runnin’ The Haskell phrase book is a new useful thingy. Nice and terse. https://typeclasses.com/phrasebook This one is also quite good https://lotz84.github.io/haskellbyexample/ I also like what FP complete is up to. Solid set of useful stuff, although a bit more emphasis towards their solutions than is common https://haskell.fpcomplete.com/learn I was fiddling with making some examples for my friends a while ago, but I think the above do a similar better job. https://github.com/philzook58/basic-ass-shit Highlights include: Makin a json request -# LANGUAGE OverloadedStrings, DeriveGeneric #-} module JsonRequest where import Data.Aeson import Network.Wreq import GHC.Generics import Control.Lens data ToDo = ToDo { userId :: Int, id :: Int, title :: String, completed :: Bool } deriving (Generic, Show) instance ToJSON ToDo instance FromJSON ToDo my_url = "https://jsonplaceholder.typicode.com/todos/1" main = do r <- get my_url print$ ((decode $r ^. responseBody) :: Maybe ToDo) -- ((decode$ r ^. responseBody) :: Maybe ToDo)

Showing a plot of a sine function

module Plot where

import Graphics.Rendering.Chart.Easy
import Graphics.Rendering.Chart.Backend.Cairo -- Chart-cairo
import Graphics.Image as I -- hip

filename = "example1_big.png"

main = do
toFile def filename $plot (line "a sine" [[ (x :: Double, sin x) | x <- [0, 0.1 .. 2 * pi]]]) plotimg <- readImageRGB VU filename -- yeah,I want the plot to pop up displayImage plotimg print "Press Enter to Quit" getLine Doing a least squares fit of some randomly created data module LeastSquares where import Numeric.LinearAlgebra n = 20 x = linspace n (-3,7::Double) y0 = 3 * x main = do noise <- randn 1 n let y = (flatten noise) + y0 let sampleMatrix = (asColumn x) ||| (konst 1 (n,1)) let sol = (sampleMatrix <\> y) print$ "Best fit is y = " ++ show (sol ! 0) ++ " * x + " ++ (show (sol ! 1))

I love Power Serious. https://www.cs.dartmouth.edu/~doug/powser.html Infinite power series using the power of laziness in something like 20 lines

http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.42.8903&rep=rep1&type=pdf Jerzy Karczmarczuk doing automatic differentiation in Haskell before it was cool. Check out Conal Elliott's stuff after.

Very simple symbolic differentiation example. When I saw this in SICP for the first time, I crapped my pants.

data Expr = X | Plus Expr Expr | Times Expr Expr | Const Double
deriv :: Expr -> Expr
deriv X = Const 1
deriv (Const _) = Const 0
deriv (Plus x y) = Plus (deriv x) (deriv y)
deriv (Times x y) = (Times (deriv x) y) Plus (Times x (deriv y))

https://www.cs.kent.ac.uk/people/staff/dat/miranda/whyfp90.pdf Why functional Programming Matters by John Hughes

https://www.cs.cmu.edu/~crary/819-f09/Backus78.pdf John Backus emphasizing escaping the imperative mindset in his 1978 Turing Award speech. A call to arms of functional programming

https://www.cs.tufts.edu/~nr/cs257/archive/richard-bird/sudoku.pdf Richard Bird defining sudoku solutions and then using equation reasoning to build a more efficient solver

#### Here's how I find useful Haskell packages:

I google. I go to hackage (if I'm in a subpage, click on "contents" in the upper right hand corner). Click on a category that seems reasonable (like "web" or something) and then sort by Downloads (DL). This at least tells me what is popular-ish. I look for tutorials if I can find them. Sometimes there is a very useful getting started snippet in the main subfile itself. Some packages are overwhelming, others aren't.

The Real World Haskell book is kind of intimidating although a lovely resource.

The wiki has a pretty rockin set of tutorials. Has some kind soul been improving it?

I forgot learn you a Haskell has a chapter on basic io

When you're ready to sit down with Haskell more, the best intro is currently the Haskell Book

You may also be interested in https://www.edx.org/course/introduction-functional-programming-delftx-fp101x-0 this MOOC

https://github.com/data61/fp-course or this Data61 course

Then there is a fun infinitude of things to learn after that.

______

More ideas for simple examples?

This post is intentionally terse.

IO is total infective poison.

standard output io

main = do
x <- getStrLn
putStrLn "Hello"
print [1,2,3]
print (Just 19022.32)
print x



mutation & loops. You probably don't want these. They are not idiomatic Haskell, and you may be losing out on some of the best lessons Haskell has to offer.

file IO

web requests

http://www.serpentine.com/wreq/tutorial.html

web serving - scotty

image processing

basic data structures

command line arguments

plotting

Parallelism and Concurrency