axiom-developer
[Top][All Lists]

## [Axiom-developer] 20070927.01.tpd.patch

 From: daly Subject: [Axiom-developer] 20070927.01.tpd.patch Date: Thu, 27 Sep 2007 22:16:22 -0500

This patch adds Martin Rubey's implementation of Gunter Rote's Pfaffian
algorithm as well as some documentation and regression test cases.

========================================================================
diff --git a/changelog b/changelog
index d085c48..c254ce6 100644
--- a/changelog
+++ b/changelog
@@ -1,3 +1,5 @@
+20070927 tpd src/input/Makefile pfaffian regression test added
+20070927 mxr src/input/pfaffian.input regression test added
20070920 tpd src/input/Makefile add bug101.input regression test
20070920 tpd src/input/bug101.input test laplace(log(z),z,w) bug 101
20070920 wxh src/algebra/laplace.spad fix laplace(log(z),z,w) bug 101
diff --git a/src/input/Makefile.pamphlet b/src/input/Makefile.pamphlet
index b2d5fe9..b9e5374 100644
--- a/src/input/Makefile.pamphlet
+++ b/src/input/Makefile.pamphlet
@@ -338,9 +338,9 @@ REGRES= algaggr.regress algbrbf.regress  algfacob.regress
alist.regress  \
oct.regress       ode.regress      odpol.regress    op1.regress \
opalg.regress     operator.regress op.regress       ovar.regress \
-    patch51.regress \
+    patch51.regress   page.regress \
patmatch.regress  pat.regress      perman.regress   perm.regress \
-    pfr1.regress      pfr.regress      page.regress     pmint.regress \
+    pfaffian.regress  pfr1.regress     pfr.regress      pmint.regress \
poly1.regress     polycoer.regress poly.regress     psgenfcn.regress \
quat1.regress     quat.regress     r20abugs.regress r20bugs.regress \
@@ -592,8 +592,9 @@ FILES= ${OUT}/algaggr.input${OUT}/algbrbf.input
${OUT}/algfacob.input \${OUT}/pascal1.input \
${OUT}/pascal.input \${OUT}/patch51.input \
-       ${OUT}/patmatch.input${OUT}/perman.input \
-       ${OUT}/perm.input${OUT}/pfr.input        ${OUT}/pfr1.input \ +${OUT}/patmatch.input ${OUT}/perman.input \ +${OUT}/perm.input     ${OUT}/pfaffian.input \ +${OUT}/pfr.input      ${OUT}/pfr1.input \${OUT}/pinch.input    ${OUT}/plotfile.input${OUT}/pollevel.input \
${OUT}/pmint.input${OUT}/polycoer.input \
${OUT}/poly1.input${OUT}/psgenfcn.input   \
@@ -869,7 +870,8 @@ DOCFILES= \
${DOC}/pat.input.dvi${DOC}/patch51.input.dvi   \
${DOC}/patmatch.input.dvi \${DOC}/pdecomp0.as.dvi       ${DOC}/perman.input.dvi \ -${DOC}/perm.input.dvi        ${DOC}/pfr1.input.dvi \ +${DOC}/perm.input.dvi        ${DOC}/pfaffian.input.dvi \ +${DOC}/pfr1.input.dvi       \
${DOC}/pfr.input.dvi${DOC}/pinch.input.dvi      \
${DOC}/plotfile.input.dvi${DOC}/plotlist.input.dvi   \
${DOC}/pollevel.input.dvi${DOC}/poly1.input.dvi      \
diff --git a/src/input/pfaffian.input.pamphlet
b/src/input/pfaffian.input.pamphlet
new file mode 100755
--- /dev/null
+++ b/src/input/pfaffian.input.pamphlet
@@ -0,0 +1,483 @@
+\documentclass{article}
+\usepackage{amsfonts}
+\usepackage{axiom}
+\begin{document}
+\title{\$SPAD/src/input pfaffian} +\author{Timothy Daly, Gunter Rote, Martin Rubey} +\maketitle +\begin{abstract} +In mathematics, the determinant of a skew-symmetric matrix can always +be written as the square of a polynomial in the matrix entries. This +polynomial is called the Pfaffian of the matrix. The Pfaffian is +nonvanishing only for$2n \times 2n$skew-symmetric matrices, in which +case it is a polynomial of degree$n$. +\end{abstract} +\eject +\tableofcontents +\eject +\section{Examples} +$$+{\rm Pfaffian}\left[ +\begin{array}{cc} +0 & a\\ +-a & 0 +\end{array} +\right] = a +$$ + +$$+{\rm Pfaffian}\left[ +\begin{array}{cccc} +0 & a & b & c\\ +-a & 0 & d & e\\ +-b & -d & 0 & f\\ +-c & -e & -f & 0 +\end{array} +\right] = af -be + dc +$$ + +$$+{\rm Pfaffian}\left[ +\begin{array}{cccc} +\begin{array}{cc} +0 & \lambda_1\\ +-\lambda_1 & 0 +\end{array} & 0 & \cdots & 0\\ +0 & +\begin{array}{cc} +0 & \lambda_2\\ +-\lambda_2 & 0 +\end{array} & & 0\\ +\vdots & & \ddots & \vdots\\ +0 & 0 & \cdots & +\begin{array}{cc} +0 & \lambda_n\\ +-\lambda_n & 0 +\end{array} +\end{array} +\right] = \lambda_1\lambda_2\cdots\lambda_n +$$ +\section{Formal definition} + +Let$A = \{a_{i,j}\}$be a$2n \times 2n$skew-symmetric matrix. +The Pfaffian of A is defined by the equation + +$$Pf(A) = \frac{1}{s^n n!}\sum_{\sigma\in{}S_{2n}}sgn(\sigma) +\prod_{i=1}^n{a_{\sigma(2i-1),\sigma(2i)}}$$ + +where$S_{2n}$is the symmetric group and$sgn(\sigma)$+is the signature of$\sigma$. + +One can make use of the skew-symmetry of$A$to avoid summing over all +possible permutations. Let$\Pi$be the set of all partitions of +$\{1, 2, \ldots, 2n\}$into pairs without regard to order. +There are$(2n - 1)!!$such partitions. An element +$\alpha \in \Pi$, can be written as + +$$\alpha=\{(i_1,j_1),(i_2,j_2),\cdots,(i_n,j_n)\}$$ + +with$i_k < j_k$and$i_1 < i_2 < \cdots < i_n$. Let + +$$+\pi= +\left[ +\begin{array}{cccccc} +1 & 2 & 3 & 4 & \cdots & 2n \\ +i_1 & j_1 & i_2 & j_2 & \cdots & j_{n} +\end{array} +\right] +$$ + +be a corresponding permutation. This depends only on the partition$\alpha$+and not on the particular choice of$\Pi$. Given a partition$\alpha$as above +define + +$$A_\alpha =sgn(\pi)a_{i_1,j_1}a_{i_2,j_2}\cdots a_{i_n,j_n}$$ + +The Pfaffian of$A$is then given by + +$$Pf(A)=\sum_{\alpha\in\Pi} A_\alpha$$ + +The Pfaffian of a$n\times n$skew-symmetric matrix for n odd is +defined to be zero. + +\subsection{Alternative definition} + +One can associate to any skew-symmetric$2n \times 2n$+matrix$A=\{a_{ij}\}$a bivector + +$$\omega=\sum_{i<j} a_{ij}~e^i\wedge e^j$$ + +where$\{e^1, e^2, \ldots,  e^{2n}\}$is the standard basis of +$\mathbb{R}^{2n}$. The Pfaffian is then defined by the equation + +$$\frac{1}{n!}\omega^n = \mbox{Pf}(A)~e^1\wedge e^2\wedge\cdots\wedge e^{2n}$$ + +here$\omega^n$denotes the wedge product of$n$copies of +$\omega^n$with itself. + +\subsection{Derivation from Determinant} + +The Pfaffian can be derived from the determinant for a skew-symmetric +matrix$A$as follows. Using Laplace's formula we can write the +determinant as + +$$\det(A)=(-1)^{p+1}a_{p1}A_{p1} + +(-1)^{p+2}a_{p2}A_{p2}+\cdots+(-1)^{n+p}a_{pn}A_{pn}$$ + +where$A_{pi}$is the$p,i$minor matrix of the matrix$A$. We further +use Laplace's formula to note that + +$$\det(A[A_{ij}]) = \vert A \vert^n$$ + +since this determinant is that of an$n \times n$matrix whose only +non-zero elements are the diagonals (each with value det(A)) and +$[A_{ij}]$is a matrix whose$i,j$th component is the corresponding$i,j$+minor matrix. In this way, following a proof by Parameswaran, we can +write the determinant as, + +$$\det(A)\equiv\Delta_n= +\left| +\begin{array}{cccc} +a_{11}&a_{12}&\cdots&a_{1n}\\ +a_{21}&a_{22}&\cdots&a_{2n}\\ +\vdots&&&\vdots\\ +a_{n1}&a_{n2}&\cdots&a_{nn} +\end{array} +\right| +$$ + +The minor of +$$+\left| +\begin{array}{cc} +a_{11}&a_{12}\\ +a_{21}&a_{22} +\end{array} +\right| +$$ +would be$\Delta_{n - 2}$. With this notation + +$$\left| +\begin{array}{cccc} +1&0&\cdots&0\\ +0&1&\cdots&0\\ +a_{31}&a_{32}&\cdots&a_{3n}\\ +\cdots&\cdots&\cdots&\cdots\\ +a_{n1}&a_{n2}&\cdots&a_{nn} +\end{array} +\right| +\times +\left| +\begin{array}{cccc} +A_{11}&A_{21}&\cdots&A_{n1}\\ +A_{12}&A_{22}&\cdots&A_{n2}\\ +A_{13}&A_{23}&\cdots&A_{n3}\\ +\cdots&\cdots&\cdots&\cdots\\ +A_{1n}&A_{2n}&\cdots&A_{nn} +\end{array} +\right| = +\left| +\begin{array}{cc} +A_{11}&A_{21}\\ +A_{12}&a_{22} +\end{array} +\right| +\Delta_n^{n-2} +$$ + +Thus + +$$\Delta_{n-2}\Delta_n^{n-1}= +\left| +\begin{array}{cc} +A_{11}&A_{21}\\ +A_{12}&a_{22} +\end{array} +\right| +\Delta_n^{n-2} +$$ + +Of course, it was only arbitrarily that we chose to remove the first +two rows, and more generically we can write + +$$A_{rr}A_{ss} - A_{rs}A_{sr} = \Delta_{rs,rs}\Delta_n$$ + +where$\Delta_{rs,rs}$is the determinant of the original matrix with +the rows$r$and$s$, as well as the columns$r$and$s$removed. The +equation above simplifies in the skew-symmetric case to + +$$A_{rs} = \sqrt{\Delta_{rs,rs}\Delta_n}$$ + +We now plug this back into the original formula for the determinant, + +$$\Delta_n= +(-1)^{p+1}a_{p1}\sqrt{\Delta_{p1,p1}\Delta_n} + +(-1)^{p+2}a_{p2}\sqrt{\Delta_{p2,p2}\Delta_n} + +\cdots + +(-1)^{n+p}a_{pn}\sqrt{\Delta_{pn,pn}\Delta_n} +$$ + +or with slight manipulation, + +$$\sqrt{\Delta_n}=(-1)^{p+1} +\left( \ a_{p1}\sqrt{\Delta_{p1,p1}} - +a_{p2}\sqrt{\Delta_{p2,p2}} + +\cdots + +(-1)^{n-1}a_{pn}\sqrt{\Delta_{pn,pn}} \ +\right) +$$ + +The determinant is thus the square of the right hand side, and so we +identify the right hand side as the Pfaffian. + +\section{Identities} + +For a$2n \times 2n$skew-symmetric matrix$A$and an arbitrary +$2n \times 2n$matrix B, +\begin{itemize} +\item$Pf(A)^2 = \det(A)$+\item$Pf(BAB^T) = \det(B)Pf(A)$+\item$Pf(\lambda{}A) = \lambda^nPf(A)$+\item$Pf(A^T) = ( - 1)^nPf(A)$+\item For a block-diagonal matrix +$$A_1\oplus A_2= +\left[ +\begin{array}{cc} +A_1 & 0 \\ +0 & A_2 +\end{array} +\right] +$$ +$$Pf(A1\oplus A2) = Pf(A_11)Pf(A_2)$$ +\item For an arbitrary$n \times n$matrix$M$: +$$\mbox{Pf} +\left[ +\begin{array}{cc} +0 & M \\ +-M^T & 0 +\end{array} +\right] = +(-1)^{n(n-1)/2}\det M +$$ +\end{itemize} +\section{Applications} + +The Pfaffian is an invariant polynomial of a skew-symmetric matrix +(note that it is not invariant under a general change of basis but +rather under a proper orthogonal transformation). As such, it is +important in the theory of characteristic classes. In particular, it +can be used to define the Euler class of a Riemannian manifold which +is used in the generalized Gauss-Bonnet theorem. + +The number of perfect matchings in a planar graph turns out to be the +absolute value of a Pfaffian, hence is polynomial time +computable. This is surprising given that for a general graph, the +problem is very difficult (so called \#P-complete). This result is used +to calculate the partition function of Ising models of spin glasses in +physics, respectively of Markov random fields in machine learning +(Globerson and Jaakkola, 2007), where the underlying graph is +planar. Recently it is also used to derive efficient algorithms for +some otherwise seemingly intractable problems, including the efficient +simulation of certain types of restricted quantum computation. + +The calculation of the number of possible ways to tile a standard +chessboard or 8-by-8 checkerboard with 32 dominoes is a simple example +of a problem which may be solved through the use of the Pfaffian +technique. There are 12,988,816 possible ways to tile a chessboard in +this manner. Specifically, 12988816 is the number of possible ways to +cover an 8-by-8 square with 32 1-by-2 rectangles. 12988816 is a square +number:$12988816 = 3604^2$). Note that 12988816 can be written in the +form:$2\times 1802^2 + 2\times 1802^2$, +where all the numbers have a digital root of 2. + +More generally, the number of ways to cover a$2n \times 2n$square with +$2n^2\$ dominoes (as calculated independently by Temperley and M.E. Fisher and
+Kasteleyn in 1961) is given by
+
+$$+\prod_{j=1}^N +\prod_{k=1}^N +\left ( 4\cos^2 \frac{\pi j}{2n + 1} + +4\cos^2 \frac{\pi k}{2n + 1} \right ) +$$
+
+This technique can be applied in many mathematics-related subjects,
+for example, in the classical, 2-dimensional computation of the
+dimer-dimer correlator function in quantum mechanics.
+
+\section{History}
+
+The term Pfaffian was introduced by Arthur Cayley, who used the term
+in 1852: The permutants of this class (from their connection with the
+researches of Pfaff on differential equations) I shall term
+Pfaffians.'' The term honors German mathematician Johann Friedrich
+Pfaff.
+
+\section{Axiom code}
+I have hacked together an algorithm to compute a Pfaffian, using an algorithm
+of Gunter Rote.  Currently it's only an .input script, but if it's useful for
+somebody else than myself, we could make it a little more professional.
+
+Martin
+
+<<*>>=
+)spool pfaffian.output
+)set message test on
+)set message auto off
+)clear all
+
+--S 1 of 12
+B0 n == matrix [[(if i=j+1 and odd? j then -1 else _
+                   if i=j-1 and odd? i then 1 else 0) _
+                     for j in 1..n] for i in 1..n]
+--R
+--R                                                                   Type:
Void
+--E 1
+
+--S 2 of 12
+PfChar(lambda, A) ==
+    n := nrows A
+    (n = 2) => lambda^2 + A.(1,2)
+    M := subMatrix(A, 3, n, 3, n)
+    r := subMatrix(A, 1, 1, 3, n)
+    s := subMatrix(A, 3, n, 2, 2)
+
+    p := PfChar(lambda, M)
+    d := degree(p, lambda)
+
+    B := B0(n-2)
+    C := r*B
+    g := [(C*s).(1,1), A.(1,2), 1]
+    if d >= 4 then
+        B := M*B
+        for i in 4..d by 2 repeat
+            C := C*B
+            g := cons((C*s).(1,1), g)
+    g := reverse! g
+
+    res := 0
+    for i in 0..d by 2 for j in 2..d+2 repeat
+        c := coefficient(p, lambda, i)
+        for e in first(g, j) for k in 2..-d by -2 repeat
+            res := res +  c * e * lambda^(k+i)
+
+    res
+--R
+--R                                                                   Type:
Void
+--E 2
+
+--S 3 of 12
+pfaffian A == eval(PfChar(l, A), l=0)
+--R
+--R                                                                   Type:
Void
+--E 3
+
+--S 4 of 12
+m:Matrix(Integer):=[[0,15],[-15,0]]
+--R
+--R        + 0    15+
+--R   (4)  |        |
+--R        +- 15  0 +
+--R                                                         Type: Matrix
Integer
+--E 4
+
+--S 5 of 12
+pfaffian m
+--R
+--R   (5)  15
+--R                                                     Type: Polynomial
Integer
+--E 5
+
+--S 6 of 12
+m1:Matrix(Polynomial(Integer)):=[[0,a,b,c],[-a,0,d,e],[-b,-d,0,f],[-c,-e,-f,0]]
+--R
+--R
+--R        + 0    a    b   c+
+--R        |                |
+--R        |- a   0    d   e|
+--R   (6)  |                |
+--R        |- b  - d   0   f|
+--R        |                |
+--R        +- c  - e  - f  0+
+--R                                              Type: Matrix Polynomial
Integer
+--E 6
+
+--S 7 of 12
+pfaffian m1
+--R
+--R   Compiling function B0 with type PositiveInteger -> Matrix Integer
+--R
+--R   (7)  a f - b e + c d
+--R                                                     Type: Polynomial
Integer
+--E 7
+
+--S 8 of 12
+(a,b,c,d,e,f):=(3,5,7,11,13,17)
+--R
+--R   (8)  17
+--R                                                        Type:
PositiveInteger
+--E 8
+
+--S 9 of 12
+m1
+--R
+--R
+--R        + 0    a    b   c+
+--R        |                |
+--R        |- a   0    d   e|
+--R   (9)  |                |
+--R        |- b  - d   0   f|
+--R        |                |
+--R        +- c  - e  - f  0+
+--R                                              Type: Matrix Polynomial
Integer
+--E 9
+
+--S 10 of 12
+a*f-b*e+d*c
+--R
+--R
+--R   (10)  63
+--R                                                        Type:
PositiveInteger
+--E 10
+
+--S 11 of 12
+n:=pfaffian m1
+--R
+--R   (11)  a f - b e + c d
+--R                                                     Type: Polynomial
Integer
+--E 9
+
+--S 12 of 12
+eval(n,['a,'b,'c,'d,'e,'f]::List(Symbol),[a,b,c,d,e,f])
+--R
+--R   (12)  63
+--R                                                     Type: Polynomial
Integer
+--E 12
+)spool
+)lisp (bye)
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} {\bf http://en.wikipedia.org/wiki/Pfaffian''}
+\bibitem{2} The statistics of dimers on a lattice, Part I'', Physica,
+vol.27, 1961, pp.1209-25, P.W. Kasteleyn.
+\bibitem{3} Propp, James (2004),
+Lambda-determinants and domino-tilings'', arXiv:math.CO/0406301.
+\bibitem{4} Globerson, Amir and Tommi Jaakkola (2007),
+Approximate inference using planar graph decomposition'',
+Advances in Neural Information Processing Systems 19, MIT Press.
+\bibitem{5} The Games and Puzzles Journal'',
+No.14, 1996, pp.204-5, Robin J. Chapman, University of Exeter
+\bibitem{6} Domino Tilings and Products of Fibonacci and Pell numbers'',
+Journal of Integer Sequences, Vol.5, 2002, Article 02.1.2,
+James A. Sellers, The Pennsylvania State University
+\bibitem{7} The Penguin Dictionary of Curious and Interesting Numbers'',
+revised ed., 1997, ISBN 0-14-026149-4, David Wells, p.182.
+\bibitem{8} A Treatise on the Theory of Determinants'',
+1882, Macmillan and Co., Thomas Muir, Online
+\bibitem{9} Skew-Symmetric Determinants'',
+The American Mathematical Monthly, vol. 61, no.2., 1954, p.116,
+S. Parameswaran Online-Subscription
+\end{thebibliography}
+\end{document}