Finding the Closest Unitary for a Given Matrix

clostest_unitary.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:f39c9d11e4dbdd82092c81ff101e701c93459b99cda4a8f4117429a03c6bf852"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Finding the Closest Unitary for a Given Matrix"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Author: Michael Goerz <goerz@physik.uni-kassel.de>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import scipy.linalg"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$\\newcommand{\\trace}{\\operatorname{tr}}$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For a given matrix $\\hat{A}$, we want to find the closest unitary matrix\n",
      "$\\hat{U}$, in the sense that the [operator norm][WP:OpNorm] (aka\n",
      "[2-norm][WP:MatrixNorm]) of their difference should be minimal.\n",
      "\n",
      "$$\\hat{U} = \\arg \\min_{\\tilde{U}} \\Vert \\hat{A}-\\tilde{U} \\Vert$$\n",
      "\n",
      "[WP:OpNorm]:     http://en.wikipedia.org/wiki/Operator_norm\n",
      "[WP:MatrixNorm]: http://en.wikipedia.org/wiki/Matrix_norm#Induced_norm"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For example, let's consider the following matrix:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "A = np.matrix([\n",
      "[-1.75900766E-02-1.15354406E-01j, 6.10816904E-03+9.49971160E-03j, 1.79090787E-02+1.33311069E-02j,-1.82163102E-03-8.77682357E-04j],\n",
      "[-9.77987203E-03+5.01950535E-03j, 8.74085180E-04+3.25580543E-04j,-6.74874670E-03-5.82800747E-03j,-1.95106265E-03+9.84138284E-04j],\n",
      "[-6.11175534E-03+2.26761191E-02j,-5.04355339E-03-2.57661178E-02j, 2.15674643E-01+8.36337993E-01j, 1.76098908E-02+1.74391779E-02j],\n",
      "[ 1.51473418E-03+1.07178057E-03j, 6.40793740E-04-1.94176372E-03j,-1.28408900E-02+2.66263921E-02j, 4.84726807E-02-3.84341687E-03j]\n",
      "])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In cases where $\\hat{A}$ of dimension $n$ is obtained from projecting down a\n",
      "unitary matrix $\\hat{U}_m$ of a larger dimension ($\\hat{A} = \\hat{P}\n",
      "\\hat{U}_m \\hat{P}$, where $\\hat{P}$ is the projector from dimension\n",
      "$m$ to dimension $n$) , one possible way to quantify the \"distance from unitarity\" is\n",
      "$$d_u(\\hat{A}) = 1 - \\frac{1}{n} \\trace[\\hat{A}^\\dagger \\hat{A}].$$\n",
      "\n",
      "This situation is common in quantum information, where $\\hat{A}$ is the\n",
      "projection of the unitary evolution of a large Hilbert space of dimension $m$\n",
      "into a small logical subspace of dimension $n$. The quantity $d_u(\\hat{A})$ is\n",
      "then simply the population lost from logical subspace.\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def delta_uni(A):\n",
      "    \"\"\" Assuming that A is a matrix obtained from projecting a unitary matrix\n",
      "        to a smaller subspace, return the loss of population of subspace, as a\n",
      "        distance measure of A from unitarity.\n",
      "        \n",
      "        Result is in [0,1], with 0 if A is unitary.\n",
      "    \"\"\"\n",
      "    return 1.0 - (A.H * A).trace()[0,0].real / A.shape[0]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For the given matrix, we lose about 80% of the population in the logical subspace:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "delta_uni(A)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "0.80861752310942581"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "A $\\hat{U}$ that minimizes $\\Vert \\hat{A} - \\hat{U}\\Vert$ can \n",
      "be calculated via a [singular value decomposition (SVD)][WP:SVD] [Reich]:\n",
      "$$\\hat{A} = \\hat{V} \\hat{\\Sigma}  \\hat{W}^\\dagger,$$\n",
      "$$ \\hat{U} = \\hat{V} \\hat{W}^\\dagger.$$\n",
      "\n",
      "[WP:SVD]: http://en.wikipedia.org/wiki/Singular_value_decomposition"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def closest_unitary(A):\n",
      "    \"\"\" Calculate the unitary matrix U that is closest with respect to the\n",
      "        operator norm distance to the general matrix A.\n",
      "\n",
      "        Return U as a numpy matrix.\n",
      "    \"\"\"\n",
      "    V, __, Wh = scipy.linalg.svd(A)\n",
      "    U = np.matrix(V.dot(Wh))\n",
      "    return U\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "U = closest_unitary(A)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The SVD also allows to calculate the distance that $\\hat{A}$ has from $\\hat{U}$.\n",
      "$$ d(\\hat{A}, \\hat{U}) = \\max_i \\vert \\sigma_i - 1\\vert, $$\n",
      "where $\\sigma_i$ are the diagonal entries of $\\hat{\\Sigma}$ from the SVD.  This\n",
      "is a more general measure of \"distance from unitarity\" than `delta_uni`"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def deltaU(A):\n",
      "    \"\"\" Calculate the operator norm distance \\Vert\\hat{A} - \\hat{U}\\Vert between\n",
      "        an arbitrary matrix $\\hat{A}$ and the closest unitary matrix $\\hat{U}$\n",
      "    \"\"\"\n",
      "    __, S, __ = scipy.linalg.svd(A)\n",
      "    d = 0.0\n",
      "    for s in S:\n",
      "        if abs(s - 1.0) > d:\n",
      "            d = abs(s-1.0)\n",
      "    return d"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For matrices obtained from projecting down from a larger Hilbert space, the maximum distance is 1. For general matrices, it can be larger."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "deltaU(A) # should be close to 1, we are *very* non-unitary"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "0.99902145939850873"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "deltaU(U) # should be zero, within machine precision"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 9,
       "text": [
        "6.6613381477509392e-16"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can double check this with the implemenation of the 2-norm in SciPy:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "scipy.linalg.norm(A-U, ord=2) # should be equal to deltaU(A)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 10,
       "text": [
        "0.99902145939850939"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Reference"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "[Reich] D.M.Reich. \"Characterisation and Identification of Unitary Dynamics Maps in Terms of Their Action on Density Matrices\" (unpublished)"
     ]
    }
   ],
   "metadata": {}
  }
 ]
}