Merge pull request #27 from XanaduAI/new_docs_hermite

New docs hermite

docs/algorithms.rst

@@ -27,6 +27,11 @@ which states that a bipartite graph with two parts having :math:`n/2` elements c
 It should be noted that improving over Ryser's algorithm is a well-known open problem: e.g. Knuth :cite:`TAOCP` asks for an arithmetic circuit for the permanent with less than :math:`2^n` operations. Also note that since the exact calculation of the permanent of :math:`(0,1)` matrices is in the \#P complete class :cite:`valiant1979complexity` the above identity shows that deciding if the hafnian of a complex matrix is larger than a given value is also in the \#P complete class.
 
 
+.. tip::
+
+   Permanents can be calculated directly using Ryser's algorithm via :func:`hafnian.perm`.
+
+
 Finally, note that the approximate methods developed for counting perfect matchings are aimed at (weighted-)graphs with real or positive entries  :cite:`2016arXiv160107518B,10.1007/3-540-36494-3_38`. Of particular note is the approximate algorithm introduced by Barvinok for matrices with non-negative entries :cite:`barvinok1999polynomial` further analyzed in Ref. :cite:`rudelson2016hafnians`.
 
 
@@ -140,11 +145,34 @@ Compare this with Björklund's algorithm, which requires :math:`O\left((A n)^3 \
 
 .. tip::
 
-   *Implemented as* :func:`hafnian.hafnian_repeated`. The vector :math:`\mathbf{m}` is passed in the variable ``rpt``. The loop hafnian calculation can be done by passing the variable ``mu`` with the values of the vector :math:`mathbf{u}` and the option ``loops=True``.
+   *Implemented as* :func:`hafnian.hafnian_repeated`. The vector :math:`\mathbf{m}` is passed in the variable ``rpt``. The loop hafnian calculation can be done by passing the variable ``mu`` with the values of the vector :math:`\mathbf{u}` and the option ``loops=True``.
+
+Batched algorithm
+-----------------
+Using the Hermite polynomials, and their connection to the matrix elements of Gaussian states and hafnians discussed in the next section, one can calculate the hafnians of all reductions of a matrix :math:`\mathbf{A}` up to a given cutoff. The reduction of matrix :math:`\mathbf{B}` is precisely the matrix :math:`\mathbf{B}_{\mathbf{m}}` obtained by repeating (or removing) the :math:`i^{\text{th}}` row and column :math:`m_i` times. Thus given a cutoff :math:`m_{\max}`, one can use the batched algorithm to calculate
+
+.. math::
+	\text{haf}\left( \mathbf{B}_\mathbf{m} \right)
+
+for all :math:`0\leq m_i < m_{\max}`, thus this function returns a tensor with :math:`n^{m_{\max}}` components.
+
+One can also use this function to calculate the same loop hafnians that Kan's algorithm returns
+
+.. math::
+	\text{lhaf}\left( \text{vid}(\mathbf{B}_\mathbf{m},\mathbf{u}_\mathbf{m}) \right)
+
+if provided also with a vector :math:`\mathbf{u}`. Note that this parameter is optional.
+
+Internally, these hafnians are calculated by using the recursion relation of the multidimensional Hermite polynomials discussed in the next section.
+
+.. tip::
+
+   *Implemented as* :func:`hafnian.hafnian_batched`. The loop hafnian calculation can be done by passing the variable ``mu`` with the values of the vector :math:`\mathbf{u}`.
+
 
 
 Approximate algorithm
--------------------------
+---------------------
 In 1999 Barvinok :cite:`barvinok1999polynomial` provided a surprisingly simple algorithm to approximate the hafnian of a symmetric matrix with positive entries. Let the matrix have entries :math:`A_{i,j}` and define the antisymmetric stochastic matrix with entries that distribute according to :math:`W_{i,j} = -W_{i,j}   \sim \mathcal{N}(0,A_{i,j})`, where :math:`\mathcal{N}(\mu,\sigma^2)` is the normal distribution with mean :math:`\mu` and variance :math:`\sigma^2`. The following now holds:
 
 .. math::

docs/gbs.rst

@@ -53,18 +53,79 @@ where :math:`{Z} = \left( \begin{smallmatrix} \mathbb{I} & 0\\ 0& -\mathbb{I} \e
 where :math:`\Gamma` is hermitian and positive definite, while :math:`C` is symmetric.
 
 
-An important property of Gaussian states is that reduced (or marginal) states of a global Gaussian state are also Gaussian. This implies that the reduced covariance matrix of a subsystem of a Gaussian state together with a reduced vector of means fully characterize a reduced Gaussian state. The reduced covariance matrix for modes :math:`S = i_1,\ldots,i_n` is obtained from the covariance matrix of the global state :math:`sigma` by keeping the columns and rows  :math:`i_1,\ldots,i_n` and :math:`i_1+\ell,\ldots,i_n+\ell` of the original covariance matrix :math:`\sigma`. Similarly one obtains the vector of means by keeping only entries :math:`i_1,\ldots,i_n` and :math:`i_1+\ell,\ldots,i_n+\ell` of the original vector of means :math:`\vec \beta`. Using the :ref:`notation <notation>` previously introduced, one can succinctly write the covariance matrix of modes :math:`S=i_1,\ldots,i_m` as :math:`\sigma_{(S)}`, and similarly the vector of means as :math:`\vec{\beta}_{(S)}`.
+
+Gaussian states in the quadrature basis
+***************************************
+
+Historically, Gaussian states are parametrized not in terms of the covariance matrix :math:`\sigma` of the complex amplitudes :math:`\alpha_j`, but rather in terms of its quadrature components, the canonical positions :math:`q_j` and canonical momenta :math:`p_j`,
+
+.. math::
+	\alpha_j = \frac{1}{\sqrt{2 \hbar}} \left( q_j+ i p_j \right),
+
+where :math:`\hbar` is a positive constant. There are at least three conventions for the value of this constant, :math:`\hbar \in \{1/2,1,2 \}`.
+
+It is convenient to write the following vector :math:`\mathbf{r} = (q_0,\ldots,q_{\ell-1},p_0,\ldots,p_{\ell-1})` whose mean is related to the vector of means :math:`\vec \beta` as
+
+.. math::
+	\vec \beta &= \frac{1}{\sqrt{\hbar}}\mathbf{R} \mathbf{r}, \\
+	\mathbf{\bar{r}} &= \sqrt{\hbar} \mathbf{R}^\dagger \vec \beta, \\
+	\mathbf{R} &= \frac{1}{\sqrt{2}}\begin{bmatrix}
+		\mathbb{I} & i \mathbb{I}  \\
+		\mathbb{I} & -i \mathbb{I}
+		\end{bmatrix}.
+
+Similarly the complex normal covariance matrix :math:`\sigma` of the variables :math:`(\alpha_0,\ldots,\alpha_{\ell-1})` is related to the normal covariance matrix :math:`\mathbf{V}` of the variables :math:`\mathbf{r} = (q_0,\ldots,q_{\ell-1},p_0,\ldots,p_{\ell-1})` as
+
+.. math::
+	\sigma &= \frac{1}{\hbar} \ \mathbf{R} \mathbf{V} \mathbf{R}^\dagger \\
+	\mathbf{V} &= {\hbar} \ \mathbf{R}^\dagger \sigma \mathbf{R}.
+
+.. tip::
+
+   * To inter convert between the complex covariance matrix :math:`\sigma` and the quadrature covariance matrix :math:`\mathbf{V}` use the functions :func:`hafnian.quantum.Qmat` and :func:`hafnian.quantum.Covmat`
+
+
+An important property of Gaussian states is that reduced (or marginal) states of a global Gaussian state are also Gaussian. This implies that the reduced covariance matrix of a subsystem of a Gaussian state together with a reduced vector of means fully characterize a reduced Gaussian state. The reduced covariance matrix for modes :math:`S = i_1,\ldots,i_n` is obtained from the covariance matrix of the global state :math:`sigma` or :math:`\mathbf{V}` by keeping the columns and rows  :math:`i_1,\ldots,i_n` and :math:`i_1+\ell,\ldots,i_n+\ell` of the original covariance matrix :math:`\sigma`. Similarly one obtains the vector of means by keeping only entries :math:`i_1,\ldots,i_n` and :math:`i_1+\ell,\ldots,i_n+\ell` of the original vector of means :math:`\vec \beta` or :math:`\mathbf{\bar{r}}`. Using the :ref:`notation <notation>` previously introduced, one can succinctly write the covariance matrix of modes :math:`S=i_1,\ldots,i_m` as :math:`\sigma_{(S)}` or :math:`\mathbf{V}_{(S)}` , and similarly the vector of means as :math:`\vec{\beta}_{(S)}` or :math:`\mathbf{\bar{r}}_{(S)}`.
 
 .. tip::
 
    * To obtain the reduced covariance matrix and vector of means for a certain subset of the modes use :func:`hafnian.quantum.reduced_gaussian`.
 
 
+Note that for :math:`\mathbf{V}` to be a valid **quantum** covariance matrix it needs to be symmetric and satisfy the uncertainty relation
+
+.. math::
+	V + i \frac{\hbar}{2} \Omega \geq 0.
+
+
+.. tip::
+
+   * To verify if a given quadrature covariance matrix is a valid quantum covariance matrix use the function :func:`hafnian.quantum.is_valid_cov`
+
+A Gaussian state is pure :math:`\varrho = \ket{\psi} \bra{\psi}` if and only if :math:`\text{det}(\mathbf{V}/\tfrac{\hbar}{2}) = 1`.
+
+.. tip::
+
+   * To verify if a given quadrature covariance matrix is a valid quantum covariance matrix and corresponds to a pure state use the function :func:`hafnian.quantum.is_pure_cov`
+
+Finally, there is a special subset of Gaussian states called **classical** whose covariance matrix satisfies
+
+.. math::
+	\mathbf{V} \geq \tfrac{\hbar}{2}\mathbb{I}.
+
+This terminology is explained in the next section when sampling is discussed.
+
+.. tip::
+
+   * To verify if a given quadrature covariance matrix is a valid quantum covariance matrix and corresponds to a classical state use the function :func:`hafnian.quantum.is_classical_cov`
+
 
 
 Gaussian states in the Fock basis
 *********************************
 In this section we use a generalization :cite:`quesada2019franck,quesada2019simulating` of the results of Hamilton et al. :cite:`hamilton2017gaussian` by providing an explicit expression for Fock basis matrix elements :math:`\langle \mathbf{m} | \rho | \mathbf{n} \rangle`, :math:`\mathbf{n} = (n_0,\ldots, n_{\ell-1}), \mathbf{m} = (m_0,\ldots, m_{\ell-1})`, of an :math:`\ell`-mode Gaussian state :math:`\rho` with covariance matrix :math:`\mathbf{\sigma}` and displacement vector :math:`\vec \beta`.
+Note that these matrix elements can also be calculated using multidimensional Hermite polynomials as shown by Dodonov et al. :cite:`dodonov1994multidimensional`. Depending on how many of these elements are required one can prefer to calculate loop hafnians or multidimensional Hermite polynomials. In particular if one only needs a few matrix elements it is more advantageous to use the formulas derived below. On the other hand if one requires **all** the matrix elements up to a certain Fock occupation cutoff it is more efficient to use the methods of Dodonov et al., which are also implemented in this library.
+
 
 We first define the following useful quantities:
 
@@ -91,6 +152,14 @@ Note that one can also obtain the probability of detecting a certain photon numb
 .. math:: p(\mathbf{n}|\varrho) = \langle \mathbf{n} | \varrho | \mathbf{n} \rangle.
 
 
+.. tip::
+
+   * To obtain the matrix element of gaussian state with quadrature covariance matrix :math:`\mathbf{V}` and vector of means :math:`\mathbf{r}` use the function :func:`hafnian.quantum.density_matrix_element`.
+
+.. tip::
+
+   * To obtain the Fock space density matrix of gaussian state with quadrature covariance matrix :math:`\mathbf{V}` and vector of means :math:`\mathbf{r}` use the function :func:`hafnian.quantum.density_matrix`.
+
 
 In the case where the Gaussian state :math:`\varrho = |\psi \rangle \langle \psi|` is pure then the matrix element
 
@@ -100,6 +169,14 @@ factorizes into a product of two amplitudes. In Ref. :cite:`quesada2019franck` i
 
 
 
+.. tip::
+
+   * To obtain the overlap of a *pure* gaussian state with quadrature covariance matrix :math:`\mathbf{V}` and vector of means :math:`\mathbf{r}` and a given Fock state :math:`\langle \mathbf{n}|` use the function :func:`hafnian.quantum.pure_state_amplitude`.
+
+.. tip::
+
+   * To obtain the Fock space state vector (ket) of a pure gaussian state with quadrature covariance matrix :math:`\mathbf{V}` and vector of means :math:`\mathbf{r}` use the function :func:`hafnian.quantum.state_vector`.
+
 
 
 
@@ -131,62 +208,5 @@ The torontonian can be thought of as a generating function for hafnians (cf. the
 
 
 
-Gaussian states in the quadrature basis
-***************************************
-
-Historically, Gaussian states are parametrized not in terms of the covariance matrix :math:`\sigma` of the complex amplitudes :math:`\alpha_j` but rather in terms of its quadrature components, the canonical positions :math:`q_j` and canonical momenta :math:`p_j` as follows
-
-.. math::
-	\alpha_j = \frac{1}{\sqrt{2 \hbar}} \left( q_j+ i p_j \right),
-
-where :math:`\hbar` is a positive constant. There are at least three conventions for the value of this constant, :math:`\hbar \in \{1/2,1,2 \}`.
-
-It is convenient to write the following vector :math:`\mathbf{r} = (q_0,\ldots,q_{\ell-1},p_0,\ldots,p_{\ell-1})` whose mean is related to the vector of means :math:`\vec \beta` as
-
-.. math::
-	\vec \beta &= \frac{1}{\sqrt{\hbar}}\mathbf{R} \mathbf{r}, \\
-	\mathbf{r} &= \sqrt{\hbar} \mathbf{R}^\dagger \vec \beta, \\
-	\mathbf{R} &= \frac{1}{\sqrt{2}}\begin{bmatrix}
-		\mathbb{I} & i \mathbb{I}  \\
-		\mathbb{I} & -i \mathbb{I}
-		\end{bmatrix}.
-
-Similarly the complex normal covariance matrix :math:`\sigma` of the variables :math:`(\alpha_0,\ldots,\alpha_{\ell-1})` is related to the normal covariance matrix :math:`\mathbf{V}` of the variables :math:`\mathbf{r} = (q_0,\ldots,q_{\ell-1},p_0,\ldots,p_{\ell-1})` as
-
-.. math::
-	\sigma &= \frac{1}{\hbar} \ \mathbf{R} \mathbf{V} \mathbf{R}^\dagger \\
-	\mathbf{V} &= {\hbar} \ \mathbf{R}^\dagger \sigma \mathbf{R}.
-
-Finally, note that for :math:`\mathbf{V}` to be a valid **quantum** covariance matrix it needs to be symmetric and satisfy the uncertainty relation
-
-.. math::
-	V + i \frac{\hbar}{2} \Omega \geq 0.
-
-.. tip::
-
-   * To interconvert between the complex covariance matrix :math:`\sigma` and the quadrature covariance matrix :math:`\mathbf{V}` use the functions :func:`hafnian.quantum.Qmat` and :func:`hafnian.quantum.Covmat`
-
-.. tip::
-
-   * To verify if a given quadrature covariance matrix is a valid quantum covariance matrix use the function :func:`hafnian.quantum.is_valid_cov`
-
-.. tip::
-
-   * To verify if a given quadrature covariance matrix is a valid quantum covariance matrix and corresponds to a pure state use the function :func:`hafnian.quantum.is_pure_cov`
-
-.. tip::
-
-   * To obtain the matrix element of gaussian state with quadrature covariance matrix :math:`\mathbf{V}` and vector of means :math:`\mathbf{r}` use the function :func:`hafnian.quantum.density_matrix_element`.
-
-.. tip::
-
-   * To obtain the Fock space density matrix of gaussian state with quadrature covariance matrix :math:`\mathbf{V}` and vector of means :math:`\mathbf{r}` use the function :func:`hafnian.quantum.density_matrix`.
 
-.. tip::
-
-   * To obtain the overlap of a *pure* gaussian state with quadrature covariance matrix :math:`\mathbf{V}` and vector of means :math:`\mathbf{r}` and a given Fock state :math:`\langle \mathbf{n}|` use the function :func:`hafnian.quantum.pure_state_amplitude`.
-
-.. tip::
-
-   * To obtain the Fock space state vector (ket) of a pure gaussian state with quadrature covariance matrix :math:`\mathbf{V}` and vector of means :math:`\mathbf{r}` use the function :func:`hafnian.quantum.state_vector`.
 

docs/gbs_sampling.rst

@@ -57,7 +57,9 @@ where :math:`p(N_{k-1}=n_{k-1},\ldots,N_0=n_0)` has already been calculated from
 
 .. tip::
 
-   * To generate samples from a gaussian state speficied by a quadrature covariance matrix use :func:`hafnian.samples.generate_hafnian_sample`.
+   * To generate samples from a gaussian state specified by a quadrature covariance matrix use :func:`hafnian.samples.generate_hafnian_sample`.
+
+      Note that the above algorithm can also be generalized to states with finite means for which one only needs to provide the mean with the optional argument ``mean``.
 
 
 Threshold detection samples
@@ -66,4 +68,15 @@ Note the arguments presented in the previous section can also be generalized to
 
 .. tip::
 
-   * To generate samples from a gaussian state speficied by a quadrature covariance matrix use :func:`hafnian.samples.generate_torontonian_sample`.
+   * To generate threshold samples from a gaussian state specified by a quadrature covariance matrix use :func:`hafnian.samples.generate_torontonian_sample`.
+
+
+Sampling of classical states
+****************************
+
+In the previous section it was mentioned that states whose covariance matrix satisfies :math:`\mathbf{V} \geq \mathbb{I}` are termed classical. These designation is due to the fact that for these states it is possible to obtain a polynomial (cubic) time algorithm to generate photon number or threshold samples :cite:`rahimi2015can`.
+
+.. tip::
+
+   * To generate photon number or threshold samples from a classical gaussian state specified by a quadrature covariance matrix use :func:`hafnian.samples.hafnian_sample_classical_state` or :func:`torontonian_sample_classical_state`.
+

docs/hermite.rst

@@ -0,0 +1,60 @@
+.. role:: raw-latex(raw)
+   :format: latex
+
+.. role:: html(raw)
+   :format: html
+.. _hermite:
+
+
+Multidimensional Hermite polynomials
+====================================
+.. sectionauthor:: Nicolás Quesada <nicolas@xanadu.ai>
+
+In this section we study the multidimensional Hermite polynomials. These polynomials where originally introduced by C. Hermite in 1865. Mizrahi :cite:`mizrahi1975generalized` provides an exhaustive reference on these subject. We however, follow the succinct treatment of Kok and Braunstein :cite:`kok2001multi`.
+
+In the next section, where we discuss Gaussian states we will explain how these polynomials relate to hafnians and loop hafnians. For the moment just let us introduce them and study their formal properties.
+
+Generating function definition
+******************************
+Given two complex vectors :math:`\alpha,\beta \in \mathbb{C}^\ell` and a symmetric matrix :math:`\mathbf{B} = \mathbf{B}^T \in \mathbb{C}^{\ell \times \ell}`,
+
+.. math::
+    G_B(\alpha,\beta) = \exp\left( \alpha \mathbf{B} \beta^T - \tfrac{1}{2}\beta \mathbf{B} \beta^T\right) = \sum_{\mathbf{m} \geq \mathbf{0}} \prod_{i=1}^{\ell} \frac{\beta_i^{n_i}}{n_i} H_{\mathbf{m}}^{(\mathbf{B})}(\alpha),
+
+where the notation :math:`\mathbf{m} \geq 0` is used to indicate that the sum goes over all vectors  in :math:` \mathbb{N}^{\ell}_0` (the set of vectors of nonnegative integers of size :math:`\ell`). This generating function provides an implicit definition of the multidimensional Hermite polynomials.
+It is also straightforward to verify that :math:`H_{\mathbf{0}}^{(\mathbf{B})}(\alpha) = 1`.
+
+In the one dimensional case, :math:`\ell=1`, one can compare the generating function above with the ones for the "probabilists' Hermite polynomials" :math:`He_n(x)` and "physicists' Hermite polynomials" :math:`H_n(x)` to find
+
+.. math::
+    He_n(x) = H_{n}^{([1])}(x), \\
+    H_n(x) = H_{n}^{([2])}(x).
+
+.. tip::
+   The multidimensional Hermite polynomials are *Implemented as* :func:`hafnian.hermite_multidimensional`.
+
+
+Recursion relation
+******************
+Based on the generating function introduced in the previous section one can derive the following recursion relation
+
+.. math::
+    H_{\mathbf{m}+\mathbf{e}_i}^{(\mathbf{B})}(\alpha) - \sum_{j=1}^\ell B_{i,j} \alpha_j H_{\mathbf{m}}^{(\mathbf{B})}(\alpha) + \sum_{j=1}^\ell B_{i,j} m_j H_{\mathbf{m}-\mathbf{e}_j}^{(\mathbf{B})}(\alpha) = 0,
+
+
+where :math:`\mathbf{e}_j` is a vector with zeros in all its entries except in the :math:`i^{\text{th}}` element.
+
+
+
+
+From this recursion relation, or by Taylor expanding the generating function, one easily finds
+
+.. math::
+    H_{\mathbf{e}_i}^{(\mathbf{B})}(\alpha) = \sum_{j=1}^\ell B_{ij} \alpha_j.
+
+
+Using this recursion relation one can calculate all the multidimensional Hermite polynomials up to a given cutoff.
+
+
+The connection between the multidimensional Hermite polynomials and **pure** Gaussian states was reported by Wolf :cite:`wolf1974canonical`, and later by Kramer, Moshinsky and Seligman :cite:`kramer1975group`. This same connection was also pointed out by Doktorov, Malkin and Man'ko in the context of vibrational modes of molecules :cite:`doktorov1977dynamical`.
+Furthermore, this connection was later generalized to **mixed** Gaussian states by Dodonov, Man'ko and Man'ko :cite:`dodonov1994multidimensional`.