Introduces renormalized Hermite Polynomials

.github/CHANGELOG.md

@@ -2,6 +2,8 @@
 
 ### New features
 
+* Introduces the renormalized hermite polynomials. These new polynomials improve the speed and accuracy of `thewalrus.quantum.state_vector` and `thewalrus.quantum.density_matrix` and also `hafnian_batched` and `hermite_multimensional` when called with the optional argument `renorm=True`. [#108](https://github.com/XanaduAI/thewalrus/pull/108)
+
 ### Improvements
 
 * Updates the reference that should be used when citing The Walrus. [#102](https://github.com/XanaduAI/thewalrus/pull/102)

include/hermite_multidimensional.hpp

@@ -48,106 +48,108 @@ ullint vec2index(std::vector<int> &pos, int resolution) {
 
 }
 
+namespace libwalrus {
+
 /**
- * Returns the indices of the tensor corresponding to a given element
+ * Returns photon number statistics of a Gaussian state for a given covariance matrix `mat`.
+ * as described in *Multidimensional Hermite polynomials and photon distribution for polymode mixed light*
+ * [arxiv:9308033](https://arxiv.org/abs/hep-th/9308033).
+ *
+ * This implementation is based on the MATLAB code available at
+ * https://github.com/clementsw/gaussian-optics
  *
- * @param val
- * @param base
- * @param n
+ * @param mat a flattened vector of size \f$2n^2\f$, representing a
+ *       \f$2n\times 2n\f$ row-ordered symmetric matrix.
+ * @param d a flattened vector of size \f$2n\f$, representing the first order moments.
+ * @param resolution highest number of photons to be resolved.
  *
- * @return tensor index
  */
-std::vector<int> find_rep(int val, int base, int n) {
-    std::vector<int> x(n, 0);
-    int local_val = val;
+template <typename T>
+inline std::vector<T> hermite_multidimensional_cpp(std::vector<T> &R_mat, std::vector<T> &y_mat, int &resolution) {
+    int dim = std::sqrt(static_cast<double>(R_mat.size()));
 
-    x[0] = 1;
+    namespace eg = Eigen;
 
-    for (int i = 1; i < n; i++)
-        x[i] = x[i-1]*base;
+    eg::Matrix<T, eg::Dynamic, eg::Dynamic> R = eg::Map<eg::Matrix<T, eg::Dynamic, eg::Dynamic>, eg::Unaligned>(R_mat.data(), dim, dim);
+    eg::Matrix<T, eg::Dynamic, eg::Dynamic> y = eg::Map<eg::Matrix<T, eg::Dynamic, eg::Dynamic>, eg::Unaligned>(y_mat.data(), dim, dim);
 
-    std::vector<int> digits(n, 0);
+    ullint Hdim = pow(resolution, dim);
+    std::vector<T> H(Hdim, 0);
 
-    for (int i = 0; i < n; i++) {
-        digits[i] = local_val/x[n-i-1];
-        local_val = local_val - digits[i] * x[n-i-1];
-    }
 
-    return digits;
-}
+    H[0] = 1;
 
+    std::vector<int> nextPos(dim, 1);
+    std::vector<int> jumpFrom(dim, 1);
+    std::vector<int> ek(dim, 0);
+    std::vector<double> factors(resolution+1, 0);
+    int jump = 0;
 
-/**
- * Returns the sqrt of the factorial of an integer.
- *
- * @param nn input integer
- *
- * @return Square root of the factorial of \f$n\f$.
- */
-long double sqrtfactorial(int nn)
-{
-    long double n = static_cast<long double>(nn);
-
-    if(n > 1)
-        return std::sqrt(n) * sqrtfactorial(n - 1);
-    else
-        return 1;
-}
+    for (ullint jj = 0; jj < Hdim-1; jj++) {
 
+        if (jump > 0) {
+            jumpFrom[jump] += 1;
+            jump = 0;
+        }
 
-/**
- * Renormalizes an unnormalized photon number statistics of a Gaussian state.
- * Based on the MATLAB code available at: https://github.com/clementsw/gaussian-optics
- *
- * @param tn unnormalized flattened vector of size \f$res**nmodes$ representing unnormalized photon number statistics
- *       \f$2n\times 2n\f$ row-ordered symmetric matrix.
- * @param nmodes number of modes
- * @param res highest number of photons to be resolved.
- *
- * @return Renormalized photon number statistics
- */
-template <typename T>
-inline std::vector<T> renormalization(std::vector<T> tn, int nmodes, int res) {
-    std::vector<long double> invsqfacts(res, 0);
-    std::vector<int> digits(nmodes, 0);
 
-    ullint Hdim = pow(res, nmodes);
+        for (int ii = 0; ii < dim; ii++) {
+            std::vector<int> forwardStep(dim, 0);
+            forwardStep[ii] = 1;
 
-    for (int i = 0; i < res; i++)
-        invsqfacts[i] = sqrtfactorial(i);
+            if ( forwardStep[ii] + nextPos[ii] > resolution) {
+                nextPos[ii] = 1;
+                jumpFrom[ii] = 1;
+                jump = ii+1;
+            }
+            else {
+                jumpFrom[ii] = nextPos[ii];
+                nextPos[ii] = nextPos[ii] + 1;
+                break;
+            }
+        }
 
-    for (ullint i = 0; i < Hdim; i++) {
-        digits = find_rep(i, res, nmodes);
-        long double pref = 1;
-        for (int j = 0; j < nmodes; j++)
-            pref *= 1.0L/invsqfacts[digits[j]];
-        tn[i] = tn[i]*static_cast<double>(pref);
-    }
+        for (int ii = 0; ii < dim; ii++)
+            ek[ii] = nextPos[ii] - jumpFrom[ii];
+
+        int k = 0;
+        for(; k < static_cast<int>(ek.size()); k++) {
+            if(ek[k]) break;
+        }
+
+        ullint nextCoordinate = vec2index(nextPos, resolution);
+        ullint fromCoordinate = vec2index(jumpFrom, resolution);
+
+
+        H[nextCoordinate] = H[nextCoordinate] + y(k, 0);
+        H[nextCoordinate] = H[nextCoordinate] * H[fromCoordinate];
+
+        std::vector<int> tmpjump(dim, 0);
+
+        for (int ii = 0; ii < dim; ii++) {
+            if (jumpFrom[ii] > 1) {
+                std::vector<int> prevJump(dim, 0);
+                prevJump[ii] = 1;
+                std::transform(jumpFrom.begin(), jumpFrom.end(), prevJump.begin(), tmpjump.begin(), std::minus<int>());
+                ullint prevCoordinate = vec2index(tmpjump, resolution);
+                H[nextCoordinate] = H[nextCoordinate] - (static_cast<T>(jumpFrom[ii]-1))*static_cast<T>(R(k,ii))*H[prevCoordinate];
+
+            }
+        }
 
-    return tn;
+    }
+    return H;
 
 }
 
 
 
-namespace libwalrus {
 
-/**
- * Returns photon number statistics of a Gaussian state for a given covariance matrix `mat`.
- * as described in *Multidimensional Hermite polynomials and photon distribution for polymode mixed light*
- * [arxiv:9308033](https://arxiv.org/abs/hep-th/9308033).
- *
- * This implementation is based on the MATLAB code available at
- * https://github.com/clementsw/gaussian-optics
- *
- * @param mat a flattened vector of size \f$2n^2\f$, representing an
- *       \f$2n\times 2n\f$ row-ordered symmetric matrix.
- * @param d a flattened vector of size \f$2n\f$, representing the first order moments.
- * @param resolution highest number of photons to be resolved.
- *
- */
+
+
+
 template <typename T>
-inline std::vector<T> hermite_multidimensional_cpp(std::vector<T> &R_mat, std::vector<T> &y_mat, int &resolution, int &renorm) {
+inline std::vector<T> renorm_hermite_multidimensional_cpp(std::vector<T> &R_mat, std::vector<T> &y_mat, int &resolution) {
     int dim = std::sqrt(static_cast<double>(R_mat.size()));
 
     namespace eg = Eigen;
@@ -157,11 +159,12 @@ inline std::vector<T> hermite_multidimensional_cpp(std::vector<T> &R_mat, std::v
 
     ullint Hdim = pow(resolution, dim);
     std::vector<T> H(Hdim, 0);
-    std::vector<double> ren_factor(Hdim, 0);
 
     H[0] = 1;
-    ren_factor[0] = 1;
-
+    std::vector<double> intsqrt(resolution+1, 0);
+    for (int ii = 0; ii<=resolution; ii++) {
+        intsqrt[ii] = std::sqrt((static_cast<double>(ii)));
+    }
     std::vector<int> nextPos(dim, 1);
     std::vector<int> jumpFrom(dim, 1);
     std::vector<int> ek(dim, 0);
@@ -203,8 +206,7 @@ inline std::vector<T> hermite_multidimensional_cpp(std::vector<T> &R_mat, std::v
         ullint nextCoordinate = vec2index(nextPos, resolution);
         ullint fromCoordinate = vec2index(jumpFrom, resolution);
 
-
-	H[nextCoordinate] = H[nextCoordinate] + y(k, 0);
+        H[nextCoordinate] = H[nextCoordinate] + y(k, 0)/(static_cast<double>(intsqrt[nextPos[k]-1]));
         H[nextCoordinate] = H[nextCoordinate] * H[fromCoordinate];
 
         std::vector<int> tmpjump(dim, 0);
@@ -215,17 +217,12 @@ inline std::vector<T> hermite_multidimensional_cpp(std::vector<T> &R_mat, std::v
                 prevJump[ii] = 1;
                 std::transform(jumpFrom.begin(), jumpFrom.end(), prevJump.begin(), tmpjump.begin(), std::minus<int>());
                 ullint prevCoordinate = vec2index(tmpjump, resolution);
-                H[nextCoordinate] = H[nextCoordinate] - (static_cast<T>(jumpFrom[ii]-1))*static_cast<T>(R(k,ii))*H[prevCoordinate];
+                H[nextCoordinate] = H[nextCoordinate] - (intsqrt[jumpFrom[ii]-1]/intsqrt[nextPos[k]-1])*static_cast<T>(R(k,ii))*H[prevCoordinate];
 
             }
         }
 
     }
-
-    if (renorm) {
-        H = renormalization(H, dim, resolution);
-    }
-
     return H;
 
 }

tests/libwalrus_unittest.cpp

@@ -1272,7 +1272,6 @@ TEST(BatchHafnian, Clements) {
                                     2.91297712e-01, -6.63359067e-01, -1.47024726e-02,  2.58202829e+00,
                                     9.62872951e-01, -2.53217806e+00,  3.10967275e+00, -1.42559575e-14};
     int res = 4;
-    int renorm = 0;
     std::vector<std::complex<double>> d4c{std::complex<double>(0.0, 0.0), std::complex<double>(0.0, 0.0), std::complex<double>(0.0, 0.0), std::complex<double>(0.0, 0.0)};
     for(int i=0; i<4; i++) {
         for(int j=0; j<4; j++) {
@@ -1283,7 +1282,7 @@ TEST(BatchHafnian, Clements) {
     // Hermite polynomials, which means that we have to mat4 * d4 as the
     // second argument, this is precisely what is done in the previous two loops
 
-    out = libwalrus::hermite_multidimensional_cpp(mat4, d4c, res, renorm);
+    out = libwalrus::hermite_multidimensional_cpp(mat4, d4c, res);
 
     for (int i = 0; i < 256; i++) {
         EXPECT_NEAR(expected_re[i], std::real(out[i]), tol2);
@@ -1305,12 +1304,11 @@ TEST(BatchHafnian, UnitRenormalization) {
 
     std::vector<std::complex<double>> out(res*res, 0.0);
 
-    int renorm = 1;
 
     for (int i = 0; i < res; i++)
         expected_re[i*res+i] = pow(0.5, static_cast<double>(i)/2.0);
 
-    out = libwalrus::hermite_multidimensional_cpp(B, d, res, renorm);
+    out = libwalrus::renorm_hermite_multidimensional_cpp(B, d, res);
 
     for (int i = 0; i < res*res; i++) {
         EXPECT_NEAR(expected_re[i], std::real(out[i]), tol2);

thewalrus/_hermite_multidimensional.py

@@ -18,6 +18,10 @@
 
 from .libwalrus import hermite_multidimensional as hm
 from .libwalrus import hermite_multidimensional_real as hmr
+
+from .libwalrus import renorm_hermite_multidimensional as rhm
+from .libwalrus import renorm_hermite_multidimensional_real as rhmr
+
 from ._hafnian import input_validation
 
 # pylint: disable=too-many-arguments
@@ -77,9 +81,15 @@ def hermite_multidimensional(R, cutoff, y=None, renorm=False, make_tensor=True,
     yt = np.real_if_close(y)
 
     if Rt.dtype == np.float and yt.dtype == np.float:
-        values = np.array(hmr(Rt, yt, cutoff, renorm=renorm))
+        if renorm:
+            values = np.array(rhmr(Rt, yt, cutoff))
+        else:
+            values = np.array(hmr(Rt, yt, cutoff))
     else:
-        values = np.array(hm(R, y, cutoff, renorm=renorm))
+        if renorm:
+            values = np.array(rhm(np.complex128(R), np.complex128(y), cutoff))
+        else:
+            values = np.array(hm(np.complex128(R), np.complex128(y), cutoff))
 
     if make_tensor:
         shape = cutoff * np.ones([n], dtype=int)
@@ -131,6 +141,4 @@ def hafnian_batched(A, cutoff, mu=None, tol=1e-12, renorm=False, make_tensor=Tru
     return hermite_multidimensional(
         -A, cutoff, y=mu, renorm=renorm, make_tensor=make_tensor, modified=True
     )
-
-
 # Note the minus signs in the arguments. Those are intentional and are due to the fact that Dodonov et al. in PRA 50, 813 (1994) use (p,q) ordering instead of (q,p) ordering