Repository: PennyLaneAI/pennylane
Pull request #1331 (2021-05-18 07:51:59)

[unitaryhack] Change the order of the covariance matrix and the vector of means internally in `default.gaussian` for `gaussian_state`
**Context:**

Change the order of the covariance matrix and the vector of means internally in default.gaussian for gaussian_state

**Description of the Change:**

Changed the order of mu and cov in gaussian_state and adapted the documentation

**Benefits:**

Improved internal code quality for `default.gaussian`.

**Possible Drawbacks:**
N/A

**Related GitHub Issues:**
Closes PennyLaneAI#1257


l45k: Hey @josh146 and @antalszava , 
I am not completyl shure if changeing the order of mu and cov is enough. I searched for other parts of the code that need to be changed, however I found none. So please guide me into the right direction if I missed something.
antalszava: Hi @l45k, what we would like is to always have the covariance matrix first and the vector of means second. The signature of the `gaussian_state` function should be good as is, what we would like to change is the return statement: from `return mu, cov` we'd like to go to `return cov, mu`.

This applies to the entire `pennylane/devices/default_gaussian.py` file and related tests: we would always like the covariance matrix to comes first and the vector of means to come second. It will likely also include subtleties when each function is being used because the code for `default.gaussian` assumes at the moment that the vector of means is first and the covariance matrix is second.

Examples where a change is required:
* `fock_prob` function
* `displacement` function (the way of how a `state` is represented)
* `vacuum_state` function and other state preparation functions (the way of how a `state` is represented)
* `photon_number` function,
* etc.

Let us know if you have further questions. :slightly_smiling_face: 
MaldoAlberto: For this problem my team and I saw that the situation of bad value conversion is found in this repository https://github.com/PennyLaneAI/pennylane-sf 
or is it necessarily in the Pennylane repository?
antalszava: Hi @MaldoAlberto, could you point me to which part of PennyLane-SF you've found?

The `default.gaussian` device is independent of PennyLane-SF. and we would like to switch to the convention that Straweberry Fields and PennyLane-SF use (cov mat first, vector of means second) within `default.gaussian`.

We did have previously an issue that was assumed to be coming from PennyLane-SF, but turned out that it was due to the documentation in this repo.
MaldoAlberto: @antalszava 

from the example code the problem is this command qml.GaussianState(mean_vector, cov_mat, wires=[0, 1]) the values that are passed are the r, V following this in the file simulator.py in the command line 111 you have

op = self._operation_map[operation](*sf_par)

that sf_par passes the r,V commands to the strawberryfields file ops.py, but the arguments must be V and r, then to add the condition to modify any method like GaussianState must be considered in the lines of code 98-104, in the case of else if we add a reverse it will work GaussianState with r and V respectively :

        # convert PennyLane parameter conventions to
        # Strawberry Fields conventions
        if operation == "CatState":
            sf_par = (par[0] * np.exp(par[1] * 1j), par[2])
        else:
            sf_pair = par
            sf_pair.reverse()

But if there are more methods, other conditions for the _operation_map method should be considered.

I hope my explanation is correct.

antalszava: The `qml.GaussianState(mean_vector, cov_mat, wires=[0, 1])` is incorrect, this is what the issue https://github.com/PennyLaneAI/pennylane-sf/issues/67 boiled down to. It should be `qml.GaussianState(cov_mat, mean_vector, wires=[0, 1])`. This is also what the current documentation for `GaussianState` reflects:
[`GaussianState(V, r, wires)`](https://pennylane.readthedocs.io/en/latest/code/api/pennylane.GaussianState.html) where `V` is the covariance matrix and `r` is the vector of means. This problem was solved with #1255. The original confusion came from the fact, that the documentation was incorrect.

Currently, the #1257 issue focuses only on the `default.gaussian` device. This device is not related to PennyLane-SF, but it's stand-alone and can execute operations correctly, including `qml.GaussianState(cov_mat, mean_vector, wires=[0, 1])`.

Internally, however, the `default.gaussian` device still uses the old order of the covariance matrix and vector of means. We would like this to be adjusted, so that everywhere, even in the internal code of `default.gaussian` we always have the cov mat coming before the vector of means for example in function signatures, when passing arguments, when returning these quantities, etc.

Hope this helps, let us know if you have further questions! :)
l45k: @antalszava  thank you for clearifying! I will change it accordingly. 
l45k: Just to make sure. The intended bahvior is to have the cov_matrix first and the mean_vector second. Do we expect the same for the return values? For example the `gaussian_state`

```
def gaussian_state(cov, mu, hbar=2.0):
    r"""Returns a Gaussian state.

    This is simply a bare wrapper function,
    since the means vector and covariance matrix
    can be passed via the parameters unchanged.

    Note that both the means vector and covariance
    matrix should be in :math:`(\x_1,\dots, \x_N, \p_1, \dots, \p_N)`
    ordering.

    Args:
        cov (array): covariance matrix. Must be dimension :math:`2N\times 2N`,
            where N is the number of modes
        mu (array): vector means. Must be length-:math:`2N`,
            where N is the number of modes
        hbar (float): (default 2) the value of :math:`\hbar` in the commutation
            relation :math:`[\x,\p]=i\hbar`

    Returns:
        tuple: the mean and the covariance matrix of the Gaussian state
    """
    ...
    return mu, cov
```

returns the mean first and covariance second.  From my understanding this should be changed as well?

Thanks in advance!
antalszava: @l45k yep, that's right! The reason behind this change would be to eliminate any confusion that might arise from the difference of the ordering in the code, so indeed, the return statement would also have to be changed.
l45k: I commited the changes. However, when I run the tests I get a lot of Failed test. I didn't had the time to check why this happens, but I gues this might be due to some other code that needs to be adjustments to run with the changed parameters.
mariaschuld: Hey @l45k, I only had a superficial look at the tests, but the values for the gradients are far off the expected values (and it is not just an ordering issue), so there still seems to be some fundamental problem with the code changes! Let us know when you have questions, or if you solved the failing tests and are ready for a re-review. 
crazy4pi314: Hey @l45k  and the PennyLane team! I wanted to let folks know I am extending the unitaryHACK deadline till the end of this week for PRs that were already in progress to give some time to iterate with maintainers. Let me know if you have questions, or if I can be of help in getting this issue closed!
l45k: Hey @mariaschuld and @antalszava, I finally manage to find the error. I missed to change the order in one case of the poly_quad_expectations. This is now corrected and all test finish successfully. So I would like to give it anothter try.
antalszava: Hi @l45k, that's great to hear! :slightly_smiling_face: Could you resolve the conflict in the `.github/CHANGELOG.md` file please?

It should involve the following locally:
1. Checking out the master branch: `git checkout master`.
2. Pulling the latest changes from there: `git pull`
3. Checking out this feature branch: `git checkout issue_1257`
4. Merging the master branch into the feature branch: `git merge master`
5. Merging should indicate that there are conflicts in the `.github/CHANGELOG.md` file. Those can be resolved by opening the file and comparing the difference between the `<<<<<<< issue_1257` and `>>>>>>> master` markers.
l45k: @antalszava sorry, I used the wrong origin when I tried to merge the newest changes. I hope it works now. Please let me know if there is anything else I can do.
antalszava: Hi @l45k, it seems that multiple checks are failing due to parts of the codebase being removed (e.g., Projector operation).
l45k: I wasn't able to get the source right. Therfore, I reseted the branch to the current HEAD of the master branch. I don't know why this closed the PR. I implemented the changes again and the test should now run without errors. Can you reopen the PR or do I have to open a new one?
codecov[bot]: # [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1331?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
> Merging [#1331](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1331?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (7fa1fe3) into [master](https://codecov.io/gh/PennyLaneAI/pennylane/commit/8e69187d7c90a43d7dbf261dc9249096d4667cee?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (8e69187) will **decrease** coverage by `9.35%`.
> The diff coverage is `97.50%`.

[![Impacted file tree graph](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1331/graphs/tree.svg?width=650&height=150&src=pr&token=9S3YzlUgGk&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1331?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)

```diff
@@            Coverage Diff             @@
##           master    #1331      +/-   ##
==========================================
- Coverage   98.21%   88.86%   -9.36%     
==========================================
  Files         160      155       -5     
  Lines       11888    11360     -528     
==========================================
- Hits        11676    10095    -1581     
- Misses        212     1265    +1053     
```


| [Impacted Files](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1331?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Coverage Δ | |
|---|---|---|
| [pennylane/devices/default\_gaussian.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1331/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2RldmljZXMvZGVmYXVsdF9nYXVzc2lhbi5weQ==) | `99.64% <97.50%> (+0.35%)` | :arrow_up: |
| [pennylane/beta/\_\_init\_\_.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1331/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2JldGEvX19pbml0X18ucHk=) | `0.00% <0.00%> (-100.00%)` | :arrow_down: |
| [pennylane/devices/tf\_ops.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1331/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2RldmljZXMvdGZfb3BzLnB5) | `0.00% <0.00%> (-100.00%)` | :arrow_down: |
| [pennylane/devices/jax\_ops.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1331/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2RldmljZXMvamF4X29wcy5weQ==) | `0.00% <0.00%> (-100.00%)` | :arrow_down: |
| [pennylane/interfaces/torch.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1331/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2ludGVyZmFjZXMvdG9yY2gucHk=) | `0.00% <0.00%> (-100.00%)` | :arrow_down: |
| [pennylane/beta/devices/\_\_init\_\_.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1331/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2JldGEvZGV2aWNlcy9fX2luaXRfXy5weQ==) | `0.00% <0.00%> (-100.00%)` | :arrow_down: |
| [pennylane/beta/devices/numpy\_ops.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1331/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2JldGEvZGV2aWNlcy9udW1weV9vcHMucHk=) | `0.00% <0.00%> (-100.00%)` | :arrow_down: |
| [pennylane/interfaces/tf.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1331/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2ludGVyZmFjZXMvdGYucHk=) | `0.00% <0.00%> (-97.44%)` | :arrow_down: |
| [pennylane/beta/devices/default\_tensor.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1331/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2JldGEvZGV2aWNlcy9kZWZhdWx0X3RlbnNvci5weQ==) | `0.00% <0.00%> (-95.23%)` | :arrow_down: |
| [pennylane/devices/default\_qubit\_jax.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1331/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2RldmljZXMvZGVmYXVsdF9xdWJpdF9qYXgucHk=) | `0.00% <0.00%> (-94.37%)` | :arrow_down: |
| ... and [44 more](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1331/diff?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1331?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1331?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Last update [8e69187...7fa1fe3](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1331?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).

antalszava: Hi @l45k, it looks the formatting check failed, as some sources were not formatted. Could you format the `tests/devices/test_default_gaussian.py` file?

Also, there seems to be a conflict in the changelog. This can be resolved locally for example by merging `master` into this feature branch and then carefully modifying the changelog file.
antalszava: Hi @l45k, it looks the formatting is still off. Note that `black-21.6b0` should be used and the command would be `black -l 100 tests/devices/test_default_gaussian.py`.
co9olguy: Hi @l45k, just wanted to check the status of this contribution. Is it ready for the team to give another round of review? Have you taken into considerations the comment from @antalszava?
l45k: Hey @antalszava, sorry to bother you again. I think I will need some guidance for this one. I only changed the order of the arguments and nothing within the code. So this error should have been there already before I touched the code. Do you have a suggestion how to solve it?
l45k: Hey @antalszava I implemented the required changes. I hope it is now all right. 

Patch:
diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md
index c3b3a0696..021870159 100644
--- a/.github/CHANGELOG.md
+++ b/.github/CHANGELOG.md
@@ -141,6 +141,9 @@
 
 <h3>Improvements</h3>
 
+* Change the order of the covariance matrix and the vector of means internally
+  in `default.gaussian`. [(#1331)](https://github.com/PennyLaneAI/pennylane/pull/1331)
+
 * Added the `id` attribute to templates, which was missing from 
   PR [(#1377)](https://github.com/PennyLaneAI/pennylane/pull/1377).
   [(#1438)](https://github.com/PennyLaneAI/pennylane/pull/1438)
@@ -183,7 +186,7 @@
 
 This release contains contributions from (in alphabetical order):
 
-Olivia Di Matteo, Josh Izaac, Romain Moyard, Ashish Panigrahi, Maria Schuld.
+Olivia Di Matteo, Josh Izaac, Leonhard Kunczik, Romain Moyard, Ashish Panigrahi, Maria Schuld.
 
 
 # Release 0.16.0 (current release)
diff --git a/pennylane/devices/default_gaussian.py b/pennylane/devices/default_gaussian.py
index 542e34342..ea1f4eb5e 100644
--- a/pennylane/devices/default_gaussian.py
+++ b/pennylane/devices/default_gaussian.py
@@ -79,7 +79,7 @@ def partitions(s, include_singles=True):
                 yield ((item_partition),) + p
 
 
-def fock_prob(mu, cov, event, hbar=2.0):
+def fock_prob(cov, mu, event, hbar=2.0):
     r"""Returns the probability of detection of a particular PNR detection event.
 
     For more details, see:
@@ -93,8 +93,8 @@ def fock_prob(mu, cov, event, hbar=2.0):
       <https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.119.170501>`_
 
     Args:
-        mu (array): length-:math:`2N` means vector
         cov (array): :math:`2N\times 2N` covariance matrix
+        mu (array): length-:math:`2N` means vector
         event (array): length-:math:`N` array of non-negative integers representing the
             PNR detection event of the multi-mode system.
         hbar (float): (default 2) the value of :math:`\hbar` in the commutation
@@ -180,17 +180,17 @@ def displacement(state, wire, alpha, hbar=2):
     """Displacement in the phase space.
 
     Args:
-        state (tuple): contains means vector and covariance matrix
+        state (tuple): contains covariance matrix and means vector
         wire (int): wire that the displacement acts on
         alpha (float): complex displacement
 
     Returns:
-        tuple: contains the vector of means and covariance matrix
+        tuple: contains the covariance matrix and the vector of means
     """
-    mu = state[0]
+    mu = state[1]
     mu[wire] += alpha.real * math.sqrt(2 * hbar)
     mu[wire + len(mu) // 2] += alpha.imag * math.sqrt(2 * hbar)
-    return mu, state[1]
+    return state[0], mu
 
 
 def squeezing(r, phi):
@@ -358,7 +358,7 @@ def vacuum_state(wires, hbar=2.0):
     """
     means = np.zeros((2 * wires))
     cov = np.identity(2 * wires) * hbar / 2
-    state = [means, cov]
+    state = [cov, means]
     return state
 
 
@@ -376,7 +376,7 @@ def coherent_state(a, phi=0, hbar=2.0):
     alpha = a * cmath.exp(1j * phi)
     means = np.array([alpha.real, alpha.imag]) * math.sqrt(2 * hbar)
     cov = np.identity(2) * hbar / 2
-    state = [means, cov]
+    state = [cov, means]
     return state
 
 
@@ -393,7 +393,7 @@ def squeezed_state(r, phi, hbar=2.0):
         array: the squeezed state
     """
     means = np.zeros((2))
-    state = [means, squeezed_cov(r, phi, hbar)]
+    state = [squeezed_cov(r, phi, hbar), means]
     return state
 
 
@@ -413,7 +413,7 @@ def displaced_squeezed_state(a, phi_a, r, phi_r, hbar=2.0):
     """
     alpha = a * cmath.exp(1j * phi_a)
     means = np.array([alpha.real, alpha.imag]) * math.sqrt(2 * hbar)
-    state = [means, squeezed_cov(r, phi_r, hbar)]
+    state = [squeezed_cov(r, phi_r, hbar), means]
     return state
 
 
@@ -429,7 +429,7 @@ def thermal_state(nbar, hbar=2.0):
         array: the thermal state
     """
     means = np.zeros([2])
-    state = [means, (2 * nbar + 1) * np.identity(2) * hbar / 2]
+    state = [(2 * nbar + 1) * np.identity(2) * hbar / 2, means]
     return state
 
 
@@ -458,25 +458,25 @@ def gaussian_state(cov, mu, hbar=2.0):
     # pylint: disable=unused-argument
 
     # Note: the internal order of mu and cov is different to the one used in Strawberry Fields
-    return mu, cov
+    return cov, mu
 
 
-def set_state(state, wire, mu, cov):
+def set_state(state, wire, cov, mu):
     r"""Inserts a single mode Gaussian into the
     state representation of the complete system.
 
     Args:
-        state (tuple): contains means vector
-            and covariance matrix of existing state
+        state (tuple): contains covariance matrix
+            and means vector of existing state
         wire (Wires): wire corresponding to the new Gaussian state
-        mu (array): vector of means to insert
         cov (array): covariance matrix to insert
+        mu (array): vector of means to insert
 
     Returns:
         tuple: contains the vector of means and covariance matrix.
     """
-    mu0 = state[0]
-    cov0 = state[1]
+    cov0 = state[0]
+    mu0 = state[1]
     N = len(mu0) // 2
 
     # insert the new state into the means vector
@@ -488,7 +488,7 @@ def set_state(state, wire, mu, cov):
     cols = ind.reshape(1, -1)
     cov0[rows, cols] = cov
 
-    return mu0, cov0
+    return cov0, mu0
 
 
 # ========================================================
@@ -496,12 +496,12 @@ def set_state(state, wire, mu, cov):
 # ========================================================
 
 
-def photon_number(mu, cov, params, hbar=2.0):
+def photon_number(cov, mu, params, hbar=2.0):
     r"""Calculates the mean photon number for a given one-mode state.
 
     Args:
-        mu (array): length-2 vector of means
         cov (array): :math:`2\times 2` covariance matrix
+        mu (array): length-2 vector of means
         params (None): no parameters are used for this expectation value
         hbar (float): (default 2) the value of :math:`\hbar` in the commutation
             relation :math:`[\x,\p]=i\hbar`
@@ -528,7 +528,7 @@ def homodyne(phi=None):
     """
     if phi is not None:
 
-        def _homodyne(mu, cov, params, hbar=2.0):
+        def _homodyne(cov, mu, params, hbar=2.0):
             """Arbitrary angle homodyne expectation."""
             # pylint: disable=unused-argument
             rot = rotation(phi)
@@ -538,7 +538,7 @@ def homodyne(phi=None):
 
         return _homodyne
 
-    def _homodyne(mu, cov, params, hbar=2.0):
+    def _homodyne(cov, mu, params, hbar=2.0):
         """Arbitrary angle homodyne expectation."""
         # pylint: disable=unused-argument
         rot = rotation(params[0])
@@ -549,13 +549,13 @@ def homodyne(phi=None):
     return _homodyne
 
 
-def poly_quad_expectations(mu, cov, wires, device_wires, params, hbar=2.0):
+def poly_quad_expectations(cov, mu, wires, device_wires, params, hbar=2.0):
     r"""Calculates the expectation and variance for an arbitrary
     polynomial of quadrature operators.
 
     Args:
-        mu (array): vector of means
         cov (array): covariance matrix
+        mu (array): vector of means
         wires (Wires): wires to calculate the expectation for
         device_wires (Wires): corresponding wires on the device
         params (array): a :math:`(2N+1)\times (2N+1)` array containing the linear
@@ -601,12 +601,12 @@ def poly_quad_expectations(mu, cov, wires, device_wires, params, hbar=2.0):
     return ex, var
 
 
-def fock_expectation(mu, cov, params, hbar=2.0):
+def fock_expectation(cov, mu, params, hbar=2.0):
     r"""Calculates the expectation and variance of a Fock state probability.
 
     Args:
-        mu (array): length-:math:`2N` vector of means
         cov (array): :math:`2N\times 2N` covariance matrix
+        mu (array): length-:math:`2N` vector of means
         params (Sequence[int]): the Fock state to return the expectation value for
         hbar (float): (default 2) the value of :math:`\hbar` in the commutation
             relation :math:`[\x,\p]=i\hbar`
@@ -615,7 +615,7 @@ def fock_expectation(mu, cov, params, hbar=2.0):
         tuple: the Fock state expectation and variance
     """
     # pylint: disable=unused-argument
-    ex = fock_prob(mu, cov, params[0], hbar=hbar)
+    ex = fock_prob(cov, mu, params[0], hbar=hbar)
 
     # var[|n><n|] = E[|n><n|^2] -  E[|n><n|]^2 = E[|n><n|] -  E[|n><n|]^2
     var = ex - ex ** 2
@@ -720,7 +720,7 @@ class DefaultGaussian(Device):
         if operation == "GaussianState":
             if len(device_wires) != self.num_wires:
                 raise ValueError(
-                    "GaussianState means vector or covariance matrix is "
+                    "GaussianState covariance matrix or means vector is "
                     "the incorrect size for the number of subsystems."
                 )
             self._state = self._operation_map[operation](*par, hbar=self.hbar)
@@ -728,9 +728,9 @@ class DefaultGaussian(Device):
 
         if "State" in operation:
             # set the new device state
-            mu, cov = self._operation_map[operation](*par, hbar=self.hbar)
+            cov, mu = self._operation_map[operation](*par, hbar=self.hbar)
             # state preparations only act on at most 1 subsystem
-            self._state = set_state(self._state, device_wires[:1], mu, cov)
+            self._state = set_state(self._state, device_wires[:1], cov, mu)
             return  # we are done here
 
         # get the symplectic matrix
@@ -740,11 +740,11 @@ class DefaultGaussian(Device):
         S = self.expand(S, device_wires)
 
         # apply symplectic matrix to the means vector
-        means = S @ self._state[0]
+        means = S @ self._state[1]
         # apply symplectic matrix to the covariance matrix
-        cov = S @ self._state[1] @ S.T
+        cov = S @ self._state[0] @ S.T
 
-        self._state = [means, cov]
+        self._state = [cov, means]
 
     def expand(self, S, wires):
         r"""Expands a Symplectic matrix S to act on the entire subsystem.
@@ -779,13 +779,13 @@ class DefaultGaussian(Device):
     def expval(self, observable, wires, par):
 
         if observable == "PolyXP":
-            mu, cov = self._state
+            cov, mu = self._state
             ev, var = self._observable_map[observable](
-                mu, cov, wires, self.wires, par, hbar=self.hbar
+                cov, mu, wires, self.wires, par, hbar=self.hbar
             )
         else:
-            mu, cov = self.reduced_state(wires)
-            ev, var = self._observable_map[observable](mu, cov, par, hbar=self.hbar)
+            cov, mu = self.reduced_state(wires)
+            ev, var = self._observable_map[observable](cov, mu, par, hbar=self.hbar)
 
         if self.shots is not None:
             # estimate the ev
@@ -797,14 +797,14 @@ class DefaultGaussian(Device):
 
     def var(self, observable, wires, par):
 
-        mu, cov = self.reduced_state(wires)
+        cov, mu = self.reduced_state(wires)
 
         if observable == "PolyXP":
             _, var = self._observable_map[observable](
-                mu, cov, wires, self.wires, par, hbar=self.hbar
+                cov, mu, wires, self.wires, par, hbar=self.hbar
             )
         else:
-            _, var = self._observable_map[observable](mu, cov, par, hbar=self.hbar)
+            _, var = self._observable_map[observable](cov, mu, par, hbar=self.hbar)
         return var
 
     def sample(self, observable, wires, par):
@@ -839,7 +839,7 @@ class DefaultGaussian(Device):
                 "default.gaussian does not support sampling {}".format(observable)
             )
 
-        mu, cov = self.reduced_state(wires)
+        cov, mu = self.reduced_state(wires)
         rot = rotation(phi)
 
         muphi = rot.T @ mu
@@ -855,14 +855,14 @@ class DefaultGaussian(Device):
         self._state = vacuum_state(self.num_wires, self.hbar)
 
     def reduced_state(self, wires):
-        r"""Returns the vector of means and the covariance matrix of the specified wires.
+        r"""Returns the covariance matrix and the vector of means of the specified wires.
 
         Args:
             wires (Wires): requested wires
 
         Returns:
-            tuple (means, cov): means is an array containing the vector of means,
-            and cov is a square array containing the covariance matrix
+            tuple (cov, means): cov is a square array containing the covariance matrix,
+            and means is an array containing the vector of means
         """
         if len(wires) == self.num_wires:
             # reduced state is full state
@@ -876,7 +876,7 @@ class DefaultGaussian(Device):
         rows = ind.reshape(-1, 1)
         cols = ind.reshape(1, -1)
 
-        return self._state[0][ind], self._state[1][rows, cols]
+        return self._state[0][rows, cols], self._state[1][ind]
 
     @property
     def operations(self):


Test patch:
diff --git a/tests/devices/test_default_gaussian.py b/tests/devices/test_default_gaussian.py
index 7e61b8ad7..696f121ac 100644
--- a/tests/devices/test_default_gaussian.py
+++ b/tests/devices/test_default_gaussian.py
@@ -16,15 +16,14 @@ Unit tests for the :mod:`pennylane.plugin.DefaultGaussian` device.
 """
 # pylint: disable=protected-access,cell-var-from-loop,no-self-use
 
-import pytest
-from scipy.special import factorial as fac
-from scipy.linalg import block_diag
 import numpy as np
 import numpy.random
+import pytest
+from scipy.linalg import block_diag
+from scipy.special import factorial as fac
 
 import pennylane as qml
 from pennylane import DeviceError
-from pennylane.wires import Wires
 from pennylane.devices.default_gaussian import (
     fock_prob,
     rotation,
@@ -41,7 +40,7 @@ from pennylane.devices.default_gaussian import (
     thermal_state,
     DefaultGaussian,
 )
-
+from pennylane.wires import Wires
 
 U = np.array(
     [
@@ -50,7 +49,6 @@ U = np.array(
     ]
 )
 
-
 U2 = np.array(
     [
         [
@@ -80,10 +78,8 @@ U2 = np.array(
     ]
 )
 
-
 H = np.array([[1.02789352, 1.61296440 - 0.3498192j], [1.61296440 + 0.3498192j, 1.23920938 + 0j]])
 
-
 hbar = 2
 
 
@@ -171,7 +167,7 @@ class TestAuxillaryFunctions:
         probs = [0.430461524043, 0.163699407559, 0.0582788388927, 0.00167706931355]
 
         for idx, e in enumerate(events):
-            res = fock_prob(mu, cov, e, hbar=hbar)
+            res = fock_prob(cov, mu, e, hbar=hbar)
             assert res == pytest.approx(probs[idx], abs=tol)
 
 
@@ -308,14 +304,14 @@ class TestStates:
     def test_vacuum_state(self, tol):
         """Test the vacuum state is correct."""
         wires = 3
-        means, cov = vacuum_state(wires, hbar=hbar)
+        cov, means = vacuum_state(wires, hbar=hbar)
         assert means == pytest.approx(np.zeros([2 * wires]), abs=tol)
         assert cov == pytest.approx(np.identity(2 * wires) * hbar / 2, abs=tol)
 
     def test_coherent_state(self, tol):
         """Test the coherent state is correct."""
         a = 0.432 - 0.123j
-        means, cov = coherent_state(a, hbar=hbar)
+        cov, means = coherent_state(a, hbar=hbar)
         assert means == pytest.approx(np.array([a.real, a.imag]) * np.sqrt(2 * hbar), abs=tol)
         assert cov == pytest.approx(np.identity(2) * hbar / 2, abs=tol)
 
@@ -323,7 +319,7 @@ class TestStates:
         """Test the squeezed state is correct."""
         r = 0.432
         phi = 0.123
-        means, cov = squeezed_state(r, phi, hbar=hbar)
+        cov, means = squeezed_state(r, phi, hbar=hbar)
 
         # test vector of means is zero
         assert means == pytest.approx(np.zeros([2]), abs=tol)
@@ -340,7 +336,7 @@ class TestStates:
         phi_a = np.angle(alpha)
         r = 0.432
         phi_r = 0.123
-        means, cov = displaced_squeezed_state(a, phi_a, r, phi_r, hbar=hbar)
+        cov, means = displaced_squeezed_state(a, phi_a, r, phi_r, hbar=hbar)
 
         # test vector of means is correct
         assert means == pytest.approx(
@@ -355,7 +351,7 @@ class TestStates:
     def thermal_state(self, tol):
         """Test the thermal state is correct."""
         nbar = 0.5342
-        means, cov = thermal_state(nbar, hbar=hbar)
+        cov, means = thermal_state(nbar, hbar=hbar)
         assert means == pytest.approx(np.zeros([2]), abs=tol)
         assert np.all((cov.diag * 2 / hbar - 1) / 2 == nbar)
 
@@ -415,12 +411,12 @@ class TestDefaultGaussianDevice:
                     mu = np.array([0.432, 0.123, 0.342, 0.123])
                     p = [cov, mu]
                     w = list(range(2))
-                    expected_out = [mu, cov]
+                    expected_out = [cov, mu]
                 elif gate_name == "Interferometer":
                     w = list(range(2))
                     p = [U]
                     S = fn(*p)
-                    expected_out = S @ gaussian_dev._state[0], S @ gaussian_dev._state[1] @ S.T
+                    expected_out = S @ gaussian_dev._state[0] @ S.T, S @ gaussian_dev._state[1]
             else:
                 # the parameter is a float
                 p = [0.432423, -0.12312, 0.324, 0.751][: op.num_params]
@@ -428,19 +424,19 @@ class TestDefaultGaussianDevice:
                 if gate_name == "Displacement":
                     alpha = p[0] * np.exp(1j * p[1])
                     state = gaussian_dev._state
-                    mu = state[0].copy()
+                    mu = state[1].copy()
                     mu[w[0]] += alpha.real * np.sqrt(2 * hbar)
                     mu[w[0] + 2] += alpha.imag * np.sqrt(2 * hbar)
-                    expected_out = mu, state[1]
+                    expected_out = state[0], mu
                 elif "State" in gate_name:
-                    mu, cov = fn(*p, hbar=hbar)
+                    cov, mu = fn(*p, hbar=hbar)
                     expected_out = gaussian_dev._state
-                    expected_out[0][[w[0], w[0] + 2]] = mu
+                    expected_out[1][[w[0], w[0] + 2]] = mu
 
                     ind = np.concatenate([np.array([w[0]]), np.array([w[0]]) + 2])
                     rows = ind.reshape(-1, 1)
                     cols = ind.reshape(1, -1)
-                    expected_out[1][rows, cols] = cov
+                    expected_out[0][rows, cols] = cov
                 else:
                     # if the default.gaussian is an operation accepting parameters,
                     # initialise it using the parameters generated above.
@@ -451,7 +447,7 @@ class TestDefaultGaussianDevice:
                         # reorder from symmetric ordering to xp-ordering
                         S = block_diag(S, np.identity(2))[:, [0, 2, 1, 3]][[0, 2, 1, 3]]
 
-                    expected_out = S @ gaussian_dev._state[0], S @ gaussian_dev._state[1] @ S.T
+                    expected_out = S @ gaussian_dev._state[0] @ S.T, S @ gaussian_dev._state[1]
 
             gaussian_dev.apply(gate_name, wires=Wires(w), par=p)
 
@@ -791,7 +787,7 @@ class TestDefaultGaussianIntegration:
                 return (alpha + a).real * np.sqrt(2 * hbar)
 
             if "State" in g:
-                mu, _ = qop(*x, hbar=hbar)
+                _, mu = qop(*x, hbar=hbar)
                 return mu[0]
 
             S = qop(*x)