Pseudoinverses.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pairs of functions $f$ and $g$ from finite sets to themselves having the property that $fgf=f$ and $gfg=g$ occur in chip firing (Todo ref), rotor routing (Todo ref) and Bulgarian solitaire.  Counting the number of ordered pairs $(f,g)$ that are pseudoinverses of each other is straightforward.\n",
    "\n",
    "Let us assume wlog that our underlying set is $[n]=\\{1, 2, \\ldots, n\\}$.\n",
    "\n",
    "First let us count the number of pairs that satisfy the first condition $fgf=f$.  The condition only constrains $g$ on the image of $f$, so if $y=f(x)$ is in the image of $f$, then $f(g(y))=y$ and so $g(y)\\in f^{-1}(y)$ for all $y$ in the image of $f$.  In particular, $g$ is injective on the image of $f$.  Let $k$ be the cardinality of the image of $f$.  There are $n!/(n-k)!$ ways of specifying $g$ on the image of $f$, as we need to specify $k$ distinct ordered values out of $[n]$.  This in turn specifies $f$ on the $k$ values of $f$, as $f(g(y)) = y$.  There are $k$ choices for each of remaining $n-k$ values of $f$, and $n$ choices of the remaining values of $g$.\n",
    "\n",
    "Combining all of these, along with the fact that there are $\\binom{n}{k}$ ways of choosing the image of $f$ gives us\n",
    "$$ \\sum_{k=0}^n \\binom{n}{k} \\frac{n!}{(n-k)!} (nk)^{n-k} $$ \n",
    "\n",
    "pairs $(f,g)$ such that $fgf=f$. Including the zero index has no effect except in the $n=0$ case where the trivial function has an image of cardinality 0.\n",
    "\n",
    "This gives us [A239768](https://oeis.org/A239768)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def count_one_sided_inverse(n):\n",
    "    return sum(binomial(n,k) * factorial(n) * (n*k)^(n-k) /factorial(n-k) for k in srange(n+1))\n",
    "print \"A239768 \", [count_one_sided_inverse(n) for n in srange(10)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Counting the pairs of functions that satisfy both constraints takes only a slight modification of the above analysis.  For values $y$ in the image of $f$ we have that $f(g(y))=y$ so $g(f(g(y)))=g(y)$ so we only need to look at values outside the image of $f$.  Here if $g(x)=z$ we need $g(f(z))=z$.  This happens whenever $z = g(y)$ for some $y$ in the image of $f$.  As $g$ is injective on the image of $f$, there are $k$ choices for each value of $g(z)$ for $z$  not in the image of $f$.  So the total count becomes\n",
    "$$ \\sum_{k=0}^n \\binom{n}{k} \\frac{n!}{(n-k)!} k^{2(n-k)}. $$\n",
    "\n",
    "This is [A277337](https://oeis.org/A277337)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def count_pseudoinverse(n):\n",
    "    return sum(binomial(n,k) * factorial(n) * k^(2*(n-k)) /factorial(n-k) for k in srange(n+1))\n",
    "print \"A239768 \", [count_pseudoinverse(n) for n in srange(10)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On a finite set a pair of pseudoinverse functions $(f,g)$ has the property that $(fg)^2=fg$ and $(gf)^2 = gf$.  Is the reverse true?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import itertools as it\n",
    "def count_coincidences(n, test):\n",
    "    count = 0\n",
    "    for f in it.product(range(n), repeat=n):\n",
    "        for g in it.product(range(n), repeat=n):\n",
    "            good = test(n,f,g)\n",
    "            if good:\n",
    "                count += 1\n",
    "    return count\n",
    "\n",
    "def list_coincidences(n, test):\n",
    "    l = []\n",
    "    for f in it.product(range(n), repeat=n):\n",
    "        for g in it.product(range(n), repeat=n):\n",
    "            good = test(n,f,g)\n",
    "            if good:\n",
    "                l.append(([x+1 for x in f], [x+1 for x in g]))\n",
    "    return l\n",
    "\n",
    "\n",
    "def check_idempotent(n, f, g):\n",
    "    for i in range(n):\n",
    "        if f[g[i]] != f[g[f[g[i]]]] or g[f[i]] != g[f[g[f[i]]]]:\n",
    "            return False\n",
    "    return True\n",
    "\n",
    "print [count_coincidences(n, check_idempotent) for n in range(1,5)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nope.\n",
    "\n",
    "So let us first look at the pairs $(f,g)$ where $fg$ is idempotent.  Let there be $k$ points $y$ where $f(g(y))=y$ then $f$ maps every point in the image of $g$ to one of these points, and $f$ is arbitrary on the points not in the image of $g$.\n",
    "\n",
    "We can divide things up\n",
    "\n",
    "- $j$ points $z$ not fixed by $fg$ where $g(z)=g(y)$ with $y$ fixed by $fg$.\n",
    "- $l$ points in the image of $g$ that are not in the image of $g(Y)$.\n",
    "\n",
    "So $l \\le j \\le n-k$.  There are $n!/(n-k)!$ ways of defining $g$ on $Y$.  There are $l!\\left\\{{j}\\atop {l}\\right\\}$ ways of assigning $g$ on the $j$ points and $k^{n-k-j}$ ways of assigning the other points. Finally there are $k^l$ ways of assigning $f$ on the $l$ points, and $n^{n-k-l}$ ways of assigning $f$ on the remainder.  So, combining everything we get\n",
    "$$\\sum_{k=1}^n \\frac{n!}{(n-k)!}\\binom{n}{k}\\sum_{j=0}^{n-k}\\binom{n-k}{j}\\sum_{l=0}^j \\binom{n-k}{l}\n",
    "l!\\left\\{{j} \\atop {l}\\right\\} k^{n-k-j}k^l n^{n-k-l}$$\n",
    "\n",
    "expanding and canceling\n",
    "\n",
    "$$\\sum_{k=1}^n \\frac{(n!)^2 {(nk)}^{n-k}}{k!}\\sum_{j=0}^{n-k}\\frac{1}{j!(n-k-j)!k^j}\\sum_{l=0}^j\\frac{k^l}{(n-k-l)!n^l}\\left\\{{j}\\atop {l}\\right\\}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def compute_idempotents(n):\n",
    "    return sum(factorial(n)^2 * (k*n)^(n-k) /factorial(k) * \n",
    "               sum((1/(factorial(j)*factorial(n-k-j)*k^j)) * \n",
    "                   sum(k^l/(factorial(n-k-l) * n^l) *\n",
    "                                    stirling_number2(j,l) \n",
    "                       for l in srange(j+1)) \n",
    "                   for j in srange(n-k+1)) \n",
    "               for k in srange(1,n+1))\n",
    "def one_sided_test(n,f,g):\n",
    "    for i in srange(n):\n",
    "        if f[g[i]] != f[g[f[g[i]]]]:\n",
    "            return False\n",
    "    return True\n",
    "\n",
    "#takes a long time for n=6 try something smaller if you are in a hurry\n",
    "print [count_coincidences(n, one_sided_test) for n in srange(1,5)]\n",
    "print [compute_idempotents(n) for n in srange(1,20)]\n",
    "print list_coincidences(2, one_sided_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It works.\n",
    "\n",
    "Now let us try to count with both constraints. \n",
    "\n",
    "The points $z$ where $g(f(z))=z$ are just the image under $g$ of the points $y$ where $f(g(y))=y$.  The only modification that we need to make are to restrict the values of $f$ on the set $n-k-l$ to $n-j$ distinct values that are either fixed by $gf$ or are mapped to one of those points by $g$ giving us\n",
    "\n",
    "$$\\sum_{k=1}^n \\frac{(n!)^2 k^{n-k}}{k!}\\sum_{j=0}^{n-k}\\frac{{(n-j)}^{n-k}}{j!(n-k-j)!k^j}\\sum_{l=0}^j\\frac{k^l}{(n-k-l)!{(n-j)}^l}\\left\\{{j}\\atop {l}\\right\\} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def compute_double_idempotents(n):\n",
    "    return sum(factorial(n)^2 * k^(n-k) / factorial(k) * \n",
    "               sum((n-j)^(n-k) / (factorial(j) * factorial(n-k-j) * k^j) * \n",
    "                   sum(stirling_number2(j,l) *  k^l / (factorial(n-k-l) * (n-j)^l ) \n",
    "                    for l in srange(j+1)) \n",
    "                for j in srange(n-k+1)) \n",
    "            for k in srange(1,n+1))\n",
    "\n",
    "\n",
    "\n",
    "print [count_coincidences(n, check_idempotent) for n in srange(1,5)]\n",
    "print [compute_double_idempotents(n) for n in srange(1,10)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This works as well. There should be some way of simplifying the idempotent counting functions.\n",
    "\n",
    "TODO:\n",
    "\n",
    "+ compute asymptotic approximations for the counts.\n",
    "+ See if there is some way of fitting this into the Moore-Penrose framework.\n",
    "+ Look at further restrictions -- e.g. $f$ and $g$ have the same limit cycles.\n",
    "+ Look at when f and g are conjugate - $g = a^{-1}fa$ for some permutation or involution $a$.\n",
    "+ Look at functions f such that any g so fg and gf is idempotent is also a pseudoinverse.  This is true for permutations "
   ]
  }
 ],
 "metadata": {
  "git": {
   "suppress_outputs": true
  },
  "kernelspec": {
   "display_name": "SageMath 7.3",
   "language": "",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}