# Unification in SymPy Enabling logical programming in Python

Unification is a way to ask questions by matching expressions against patterns. It is a powerful form of pattern matching found in logical programming languages like Prolog, Maude, and Datalog. It is the computational backbone behind the logical programming paradigm and is now a part of SymPy (in a pull request).

Consider the following example. Imagine that you want to find the name of the MatrixSymbol within the Transpose in the following expression (i.e. we’re looking for the string `'X'`

)

```
>>> X = MatrixSymbol('X', 3, 3)
>>> Y = MatrixSymbol('Y', 3, 3)
>>> expr = Transpose(X) + Y
```

Traditionally we could solve this toy problem with a simple function

```
def name_of_symbol_in_transpose_in_add(matadd):
for arg in matadd.args:
if isinstance(arg, Transpose) and isinstance(arg.arg, MatrixSymbol):
return arg.arg.name
```

We solve this task with unification by setting up a pattern and then unifying that pattern against a target expression

```
>>> A = MatrixSymbol('name', n, m)
>>> B = MatrixSymbol('B', m, n)
# Look for an expression tree like A.T + B
# Treat the leaves 'name', n, m, B as Wilds
>>> pattern = Transpose(A) + B
>>> wilds = 'name', n, m, B
>>> unify(pattern, expr, wilds=wilds).next()
{'name': 'X', m: 3, n: 3, B: Y}
```

We get back a matching for each of the wildcards (name, n, m, B) and see that `'name'`

was matched to the string `'X'`

. Is this better or worse than the straight Python solution? Given the relative number of users between Python and Prolog it’s a safe bet that the style of Python programs have some significant advantages over the logical programming paradigm. Why would we program in this strange way?

Unification allows a clean separation between *what we’re looking for* and *how we find it*. In the Python solution the mathematical definition of what we want is spread among a few lines and is buried inside of control flow.

```
for arg in matadd.args:
if isinstance(arg, Transpose) and isinstance(arg.arg, MatrixSymbol):
return arg.arg.name
```

In the unification solution the lines

```
pattern = Transpose(A) + B
wilds = 'name', n, m, B
```

expresse exactly *what* we’re looking for and gives no information on *how* it should be found. The how is wrapped up in the call to `unify`

`unify(pattern, expr, wilds=wilds).next()`

This separation of the *what* and *how* is what excites me about declarative programming. I think that this separation is useful when mathematical and algorithmic programmers need to work together to solve a large problem. This is often the case in scientific computing. Mathematical programmers think about *what* should be done while algorithmic programmers think about *how* it can be efficiently computed. Declarative techniques like unification enables these two groups to work independently.

## Multiple Matches

Lets see how unify works on a slightly more interesting expression

```
>>> expr = Transpose(X) + Transpose(Y)
>>> unify(pattern, expr)
<generator object unify at 0x548cb90>
```

In this situation because both matrices `X`

and `Y`

are inside transposes our pattern to match “the name of a symbol in a transpose” could equally well return the strings `'X'`

or `'Y'`

. The unification algorithm will give us both of these options

```
>>> for match in unify(pattern, expr):
... print match
{'name': 'Y', m: 3, n: 3, B: 'X'}
{'name': 'X', m: 3, n: 3, B: 'Y'}
```

Because expr is commutative we can match `{A: Transpose(X), B: Transpose(Y)}`

or `{A: Transpose(Y), B: Transpose(X)}`

with equal validity. Instead of choosing one `unify`

, returns an iterable of all possible matches.

## Combinatorial Blowup

In how many ways can we match the following pattern

```
w + x + y + z
```

to the following expression?

```
a + b + c + d + e + f
```

This is a variant on the standard “N balls in K bins” problem often given in a discrete math course. The answer is “quite a few.” How can we avoid this combinatorial blowup?

`unify`

produces matches lazily. It returns a Python generator which yields results only as you ask for them. You can ask for just one match (a common case) very quickly.

The bigger answer is that if you aren’t satisfied with this and want a better/stronger/faster way to find your desired match you could always *rewrite unify*. The `unify`

function is all about the *how* and is disconnected from the *what*. Algorithmic programmers can tweak unify without disrupting the mathematical code.

## Rewrites

Unification is commonly used in term rewriting systems. Here is an example

```
>>> sincos_to_one = rewriterule(sin(x)**2 + cos(x)**2, 1, wilds=[x])
>>> sincos_to_one(sin(a+b)**2 + cos(a+b)**2).next()
1
```

We were able to turn a mathematical identity `sin(x)**2 + cos(x)**2 => 1`

into a function very simply using unification. However unification only does exact pattern matching so we can only find the `sin(x)**2 + cos(x)**2`

pattern if that pattern is at the top node in the tree. As a result we’re not able to apply this simplification within a larger expression tree

```
>>> list(sincos_to_one(2 + c**(sin(a+b)**2 + cos(a+b)**2))) # no matches
[]
```

I will leave the solution of this problem to a future post. Instead, I want to describe why I’m working on all of this.

## Matrix Computations

My last post was about translating Matrix Expressions into high-performance Fortran Code. I ended this post with the following problem:

*So how can we transform a matrix expression like*

`(alpha*A*B).I * x`

…

*Into a graph of BLAS calls like one of the following?*

```
DGEMM(alpha, A, B, 0, B) -> DTRSV(alpha*A*B, x)
DTRMM(alpha, A, B) -> DTRSV(alpha*A*B, x)
```

This problem can be partially solved by unification and rewrite rules. Each `BLAS`

operation is described by a class

```
class MM(BLAS):
""" Matrix Multiply """
_inputs = (alpha, A, B, beta, C)
_outputs = (alpha*A*B + beta*C,)
```

The `_outputs`

and `_inputs`

fields mathematically define when `MM`

is appropriate. This is all we need to make a transformation

```
source = MM._outputs[0]
target = MM(*MM._inputs)
wilds = MM._inputs
rewriterule(source, target, wilds)
```

Unification allows us to describe `BLAS`

mathematically without thinking about
how each individual operation will be detected in an expression. The control
flow and the math are completely separated allowing us to think hard about each
problem in isolation.

## References

I learned a great deal from the following sources

- Artificial Intelligence: A Modern Approach by Stuart Russel and Peter Norvig (Particularly section 9.2 in the second edition)
- StackOverflow - Algorithms for Unification of list-based trees
- StackOverflow - Partition N items into K bins in Python lazily (Special thanks to Chris Smith who provided the best answer)
- Logic Programming
- Term Rewriting
- My favorite Prolog tutorial
- SymPy E-mail thread on this topic
- Pull Request

blog comments powered by Disqus