Symbolic Reasoning about Quantum Circuits in Coq

A Hilbert space $\mathcal{H}$ is a complete vector space equipped with an inner product

$$\langle\cdot|\cdot\rangle:\mathcal{H}\times\mathcal{H}\rightarrow\mathbb{C}$$

such that

1. 1.

   $\langle\psi|\psi\rangle\geq 0$ for any $|\psi\rangle\in\mathcal{H}$, with equality if and only if $|\psi\rangle=0$;

2. 2.

   $\langle\phi|\psi\rangle=\langle\psi|\phi\rangle^{\ast}$;

3. 3.

   $\langle\phi|\sum_{i}c_{i}|\psi_{i}\rangle=\sum_{i}c_{i}\langle\phi|\psi_{i}\rangle$,

where $\mathbb{C}$ is the set of complex numbers, and for each $c\in\mathbb{C}$, $c^{\ast}$ stands for the complex conjugate of $c$. A vector $|\psi\rangle\in\mathcal{H}$ is normalised if its length $\sqrt{\langle\psi|\psi\rangle}$ is equal to $1$. Two vectors $|\psi\rangle$ and $|\phi\rangle$ are orthogonal if $\langle\psi|\phi\rangle=0$. An orthonormal basis of a Hilbert space $\mathcal{H}$ is a basis $\{|i\rangle\}$ where each $|i\rangle$ is normalised and any pair of them is orthogonal.

Let $\mathcal{L(H)}$ be a set of linear operators on $\mathcal{H}$. For any $A\in\mathcal{L(H)}$, $A$ is Hermitian if $A^{\dagger}=A$ where $A^{\dagger}$ is the adjoint operator of $A$ such that $\langle\psi|A^{\dagger}|\phi\rangle=\langle\phi|A|\psi\rangle^{*}$ for any $|\psi\rangle,|\phi\rangle\in\mathcal{H}$. The fundamental spectral theorem states that the set of all normalised eigenvectors of a Hermitian operator in $\mathcal{L(H)}$ constitutes an orthonormal basis for $\mathcal{H}$. That is, there exists a so-called spectral decomposition for each Hermitian $A$ such that

$$A~{}=~{}\sum_{i}\lambda_{i}|i\rangle\langle i|,$$

where the set $\{|i\rangle\}$ constitutes an orthonormal basis of $\mathcal{H}$, $\{\lambda_{i}\}$ denotes the set of eigenvalues of $A$, and $|i\rangle\langle i|$ is the projector to the corresponding eigenspace of $\lambda_{i}$. A linear operator $A\in\mathcal{L(H)}$ is unitary if $A^{\dagger}A=AA^{\dagger}=I_{\mathcal{H}}$ where $I_{\mathcal{H}}$ is the identity operator on $\mathcal{H}$. The trace of $A$ is defined as ${\rm tr}(A)=\sum_{i}\langle i|A|i\rangle$ for some given orthonormal basis $\{|i\rangle\}$ of $\mathcal{H}$. It is worth noting that the trace function is actually independent of the orthonormal basis selected. It is also easy to check that the trace function is linear and ${\rm tr}(AB)={\rm tr}(BA)$ for any $A,B\in\mathcal{L(H)}$.

Let $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$ be two Hilbert spaces. Their tensor product $\mathcal{H}_{1}\otimes\mathcal{H}_{2}$ is defined as a vector space consisting of linear combinations of the vectors $|\psi_{1}\psi_{2}\rangle=|\psi_{1}\rangle|\psi_{2}\rangle=|\psi_{1}\rangle\otimes|\psi_{2}\rangle$ with $|\psi_{1}\rangle\in\mathcal{H}_{1}$ and $|\psi_{2}\rangle\in\mathcal{H}_{2}$. Here the tensor product of two vectors is defined by a new vector such that

$$\left(\sum_{i}\lambda_{i}|\psi_{i}\rangle\right)\otimes\left(\sum_{j}\mu_{j}|\phi_{j}\rangle\right)=\sum_{i,j}\lambda_{i}\mu_{j}|\psi_{i}\rangle\otimes|\phi_{j}\rangle.$$

Then $\mathcal{H}_{1}\otimes\mathcal{H}_{2}$ is also a Hilbert space where the inner product is defined in the following way: for any $|\psi_{1}\rangle,|\phi_{1}\rangle\in\mathcal{H}_{1}$ and $|\psi_{2}\rangle,|\phi_{2}\rangle\in\mathcal{H}_{2}$,

$$\langle\psi_{1}\otimes\psi_{2}|\phi_{1}\otimes\phi_{2}\rangle=\langle\psi_{1}|\phi_{1}\rangle_{\mathcal{H}_{1}}\langle\psi_{2}|\phi_{2}\rangle_{\mathcal{H}_{2}},$$

where $\langle\cdot|\cdot\rangle_{\mathcal{H}_{i}}$ is the inner product of $\mathcal{H}_{i}$. For any $A_{1}\in\mathcal{L}(\mathcal{H}_{1})$ and $A_{2}\in\mathcal{L}(\mathcal{H}_{2})$, $A_{1}\otimes A_{2}$ is defined as a linear operator in $\mathcal{L}(\mathcal{H}_{1}\otimes\mathcal{H}_{2})$ such that for each $|\psi_{1}\rangle\in\mathcal{H}_{1}$ and $|\psi_{2}\rangle\in\mathcal{H}_{2}$,

$$(A_{1}\otimes A_{2})|\psi_{1}\psi_{2}\rangle=A_{1}|\psi_{1}\rangle\otimes A_{2}|\psi_{2}\rangle.$$

According to von Neumann’s formalism of quantum mechanics vN55 , an isolated physical system is associated with a Hilbert space which is called the state space of the system. A pure state of a quantum system is a normalised vector in its state space, and a mixed state is represented by a density operator on the state space. Here a density operator $\rho$ on Hilbert space $\mathcal{H}$ is a positive linear operator such that ${\rm tr}(\rho)=1$. Another equivalent representation of a density operator is an ensemble NC11 of pure states. In particular, given an ensemble $\{(p_{i},|\psi_{i}\rangle)\}$ where $p_{i}\geq 0$, $\sum_{i}p_{i}=1$, and $|\psi_{i}\rangle$ are pure states, then $\rho=\sum_{i}p_{i}|\psi_{i}\rangle\langle\psi_{i}|$ is a density operator. Conversely, each density operator can be generated by an ensemble of pure states in this way. Finally, a pure state can be regarded as a special mixed state.

The state space of a composite system (for example, a quantum system consisting of many qubits) is the tensor product of the state spaces of its components. Note that in general, the state of a composite system cannot be decomposed into a tensor product of the reduced states on its component systems. A well-known example is the 2-qubit state

$$|\Psi\rangle=\frac{1}{\sqrt{2}}(|00\rangle+|11\rangle).$$

This kind of state is called an entangled state. Entanglement is an important feature of quantum mechanics which has no counterpart in the classical world, and is the key to many quantum information processing tasks.

Let $|\psi\rangle$ be a state vector and $\theta$ a real number. In quantum mechanics, the state $e^{i\theta}|\psi\rangle$ is considered to be equal to $|\psi\rangle$, up to the global phase factor $e^{i\theta}$. The reason is that from an observational point of view, global phases are irrelevant to the observed properties of the physical system under consideration and can thus be ignored as far as quantum states are concerned NC11 .

The evolution of a closed quantum system is described by a unitary operator on its state space. If the states of the system at times $t_{1}$ and $t_{2}$ are $|\psi_{1}\rangle$ and $|\psi_{2}\rangle$, respectively, then $|\psi_{2}\rangle=U|\psi_{1}\rangle$ for some unitary operator $U$ which depends only on $t_{1}$ and $t_{2}$. A convenient way to understand unitary operators is in terms of their matrix representations. In fact, the unitary operator and matrix viewpoints turn out to be completely equivalent. An $m$ by $n$ complex unitary matrix $U$ with entries $U_{ij}$ can be considered as a unitary operator sending vectors in the vector space $\mathbb{C}^{n}$ to the vector space $\mathbb{C}^{m}$, under matrix multiplication of the matrix $U$ by a vector in $\mathbb{C}^{n}$.

We often denote a single qubit as a vector $|\psi\rangle=\alpha|0\rangle+\beta|1\rangle$ parameterized by two complex numbers satisfying the condition $|\alpha|^{2}+|\beta|^{2}=1$. A unitary operator for a qubit is then described by a $2\times 2$ unitary matrix. Quantum circuits are a popular model for quantum computation, where quantum gates usually stand for basic unitary operators whose mathematical meanings are given by appropriate unitary matrices. Some commonly used quantum gates to appear in the current work include the $1$-qubit Hadamard gate $H$, the Pauli gates $I_{2},X,Y,Z$, the controlled-NOT gate $CX$ performed on two qubits, and the $3$-qubit Toffoli gate. Their matrix representations are given below:

$$I_{2}=\left(\begin{array}[]{cc}1&0\\ 0&1\\ \end{array}\right),\qquad X=\left(\begin{array}[]{cc}0&1\\ 1&0\\ \end{array}\right),\qquad Y=\left(\begin{array}[]{cc}0&-i\\ i&0\\ \end{array}\right),\qquad Z=\left(\begin{array}[]{cc}1&0\\ 0&-1\\ \end{array}\right),$$

$$H=\frac{1}{\sqrt{2}}\left(\begin{array}[]{cc}1&1\\ 1&-1\\ \end{array}\right),\qquad CX=\left(\begin{array}[]{cccc}1&0&0&0\\ 0&1&0&0\\ 0&0&0&1\\ 0&0&1&0\end{array}\right),\qquad TOF=\left(\begin{array}[]{cccccccc}1&0&0&0&0&0&0&0\\ 0&1&0&0&0&0&0&0\\ 0&0&1&0&0&0&0&0\\ 0&0&0&1&0&0&0&0\\ 0&0&0&0&1&0&0&0\\ 0&0&0&0&0&1&0&0\\ 0&0&0&0&0&0&0&1\\ 0&0&0&0&0&0&1&0\end{array}\right).$$

For example, in Figure 1 we display a circuit that can generate the $3$-qubit Greenberger–Horne–Zeilinger state (GHZ state) GHZ , which is $\frac{|000\rangle+|111\rangle}{\sqrt{2}}$. In the circuit, a Hadamard gate is applied on the first qubit, then two controlled-NOT gates are used, with the first qubit controlling the second, which in turn controls the third.

A quantum measurement is described by a collection $\{M_{m}\}$ of measurement operators, where the indices $m$ refer to the measurement outcomes. It is required that the measurement operators satisfy the completeness equation $\sum_{m}M_{m}^{{\dagger}}M_{m}=I_{\mathcal{H}}$. If the state of a quantum system is $|\psi\rangle$ immediately before the measurement, then the probability that result $m$ occurs is given by

$$p(m)=\langle\psi|M_{m}^{\dagger}M_{m}|\psi\rangle,$$

and the state of the system after the measurement is $\frac{M_{m}|\psi\rangle}{\sqrt{p(m)}}.$ If the states of the system at times $t_{1}$ and $t_{2}$ are mixed, say $\rho_{1}$ and $\rho_{2}$, respectively, then $\rho_{2}=U\rho_{1}U^{{\dagger}}$ after the unitary operation $U$ is applied on the system. For the same measurement $\{M_{m}\}$ as above, if the system is in the mixed state $\rho$, then the probability that measurement result $m$ occurs is given by

$$p(m)={\rm tr}(M_{m}^{{\dagger}}M_{m}\rho),$$

and the state of the post-measurement system is $\frac{M_{m}\rho M_{m}^{{\dagger}}}{p(m)}$ provided that $p(m)>0$.

Our symbolic reasoning is based on terms constructed from scalars and basic vectors using some constructors:

• •

  Scalars are complex numbers. We write $\mathbb{C}$ for the set of complex numbers. Our formal treatment of complex numbers is based on the definitions and lemmas in the library of Paykin et al. PRZ17 for complex numbers, which, in turn, is adapted from the Coquelicot library BLM .

• •

  Basic vectors are the base states of a qubit, i.e., $|0\rangle$ and $|1\rangle$ in the Dirac notation. Mathematically, $|0\rangle$ stands for the vector $[1\ 0]^{T}$ and $|1\rangle$ for $[0\ 1]^{T}$.

• •

  Constructors include the scalar product $\cdot$, matrix product $\times$, matrix addition $+$, tensor product $\otimes$, and the conjugate transpose $A^{\dagger}$ of a matrix $A$.

In Dirac notations, $\langle 0|$ represents the dual of $|0\rangle$, i.e. $|0\rangle^{\dagger}$; similarly for $\langle 1|$. The term $\langle j|\times|k\rangle$ is abbreviated into $\langle j|k\rangle$, for any $j,k\in\{0,1\}$. This notation introduces an intuitive explanation of quantum operation. For example, the effect of the $X$ operator is to map $|0\rangle$ into $|1\rangle$ and $|1\rangle$ into $|0\rangle$. Thus we define $X$ in Coq as $|0\rangle\langle 1|+|1\rangle\langle 0|$, instead of $\begin{bmatrix}0\,1\\ 1\,0\end{bmatrix}$. Then it is obvious that $X|0\rangle=|1\rangle\langle 0|0\rangle+|0\rangle\langle 1|0\rangle=|1\rangle$ and similarly for $X|1\rangle$.

Some commonly used vectors and gates can be derived from the basic terms. For example, we define the vectors $|+\rangle$ and $|-\rangle$ as follows.

$$|+\rangle~{}=~{}\frac{1}{\sqrt{2}}\cdot|0\rangle+\frac{1}{\sqrt{2}}\cdot|1\rangle\qquad\qquad|-\rangle~{}=~{}\frac{1}{\sqrt{2}}\cdot|0\rangle+(-\frac{1}{\sqrt{2}})\cdot|1\rangle$$

We also define four simple matrices $B_{0},...,B_{3}$.

$$B_{0}~{}=~{}|0\rangle\times\langle 0|\quad\quad B_{1}~{}=~{}|0\rangle\times\langle 1|\quad\quad B_{2}~{}=~{}|1\rangle\times\langle 0|\quad\quad B_{3}~{}=~{}|1\rangle\times\langle 1|$$

(1)

The Hadamard matrix $H$ is a combination of the four matrices above.

$$H~{}=~{}\frac{1}{\sqrt{2}}\cdot B_{0}+\frac{1}{\sqrt{2}}\cdot B_{1}+\frac{1}{\sqrt{2}}\cdot B_{2}+(-\frac{1}{\sqrt{2}})\cdot B_{3}$$

Similarly, the Pauli-X gate and the controlled-NOT gate $CX$ are given as follows.

$$X~{}=~{}B_{1}+B_{2}\qquad\qquad CX~{}=~{}B_{0}\otimes I_{2}+B_{3}\otimes X$$

(2)

Notice that tensor products of matrices from $\{B_{j}:j=0,1,2,3\}$ constitute an orthonormal basis for the space of all matrices of dimension $2^{n}$, for all $n\geq 1$. Thus any matrix that appeared in the computation of a quantum circuit can be represented as a term.

Suppose the state of a quantum system is represented by a vector. The central idea of our symbolic reasoning is to employ the laws in Table 1 to rewrite terms, trying to put together the basic vectors and simplify them using the laws in L1. Technically, we design a series of strategies for that purpose.

Let $T_{1}=T_{2}$ be a law that relates two terms $T_{1}$ and $T_{2}$. We say that the law is sound if both terms can be reduced to exactly the same matrix, after their metavariables (e.g., $\ket{0}$) are instantiated by concrete matrices (e.g., $[1\ 0]^{T}$). We define a strategy called orthogonal$\_$reduce to verify that the laws in L1 are sound. In this case, we explicitly represent $|0\rangle$ and $|1\rangle$ as matrices. Both of them are $2\times 1$ matrices, so we can use matrix multiplication to prove that the corresponding elements in the matrices on both sides of the equation are the same. For example, the law $\langle 0|0\rangle=1$ is actually shown via $[1\ 0]\begin{bmatrix}1\\ 0\end{bmatrix}=1$. We add the laws in L1 to the set named S_db. The laws in L2-16 are also proved through explicit matrix representation.

We have given a formal soundness proof in Coq for each of the laws in Table 1. We collect the soundness properties in a library that contains many useful properties about matrices.

The four matrices $B_{0},...,B_{3}$ defined in (1) are the basic building blocks to represent $2\times 2$ matrices, and are intensively used in our symbolic reasoning. We design a strategy called base$\_$reduce to prove some equations about them acting on the base states $|0\rangle$ and $|1\rangle$. For example, let us consider the equation $B_{0}\times|0\rangle=|0\rangle$. We first represent $B_{0}$ by $|0\rangle\times\langle 0|$, then use the associativity of matrix multiplication to form the subterm $\langle 0|\times|0\rangle$. Now we can use the proved laws in L1 to rewrite $\langle 0|\times|0\rangle$ into $1$. The last step is to deal with scalar multiplications. We add these equations to the set named B_db.

The Pauli and Hadamard gates are probably the most widely used single-qubit gates. We introduce a strategy called gate$\_$reduce to prove some equations about the matrices $I_{2},X,Y,Z,H$ acting on base states. For example, consider the equation $X\times|0\rangle=|1\rangle$. We first expand $X$ into $B_{1}+B_{2}$. In order to prove the equation $(B_{1}+B_{2})\times|0\rangle=|1\rangle$, we use the distributivity of matrix multiplication over addition to rewrite the left hand side of the equation into the sum of $B_{1}\times|0\rangle$ and $B_{2}\times|0\rangle$. Then we employ the proved laws in B_db to rewrite them. Eventually, we deal with scalar multiplications and cancel zero matrices. We add these equations to the set named G_db.

Furthermore, we add into S_db and G_db some commonly used equations about the matrices $B_{0},...,B_{3}$ and $I_{2},X,Y,Z,H$ acting on states $|+\rangle$ and $|-\rangle$. For example, we have $H|+\rangle=|0\rangle$ and $H|-\rangle=|1\rangle$, etc.

We propose a strategy called operate$\_$reduce that puts together all the results above to reason about circuits applied to input states represented in the vector form. We will have a close look at the strategy by using an example. Let us revisit the 3-qubit GHZ state. The state can be generated by applying the circuit in Figure 1 to the initial state $|0\rangle\otimes|0\rangle\otimes|0\rangle$. We would like to verify that the output state is indeed what we expect by establishing the following equation.

$$\begin{array}[]{rl}&(I_{2}\otimes CX)\times(CX\otimes I_{2})\times(H\otimes I_{2}\otimes I_{2})\times(|0\rangle\otimes|0\rangle\otimes|0\rangle)\\ =&\frac{1}{\sqrt{2}}\cdot(|0\rangle\otimes|0\rangle\otimes|0\rangle)+\frac{1}{\sqrt{2}}\cdot(|1\rangle\otimes|1\rangle\otimes|1\rangle)\end{array}$$

(3)

We first use the right associativity of tensor product and matrix multiplication to change the expression on the left hand side in (3) into

$$\begin{array}[]{rl}&(I_{2}\otimes CX)\times((CX\otimes I_{2})\times((H\otimes(I_{2}\otimes I_{2}))\times(|0\rangle\otimes(|0\rangle\otimes|0\rangle)))).\end{array}$$

(4)

Then we calculate the inner layer matrix multiplications in sequence. We first calculate $(H\otimes(I_{2}\otimes I_{2}))\times(|0\rangle\otimes(|0\rangle\otimes|0\rangle))$. We exploit the law in L13 to change a matrix product of tensored terms into a tensor product of matrix multiplications. Then we employ the laws established in G_db and B_db to rewrite terms. This transformation is as follows.

$$\begin{array}[]{rl}&(H\otimes(I_{2}\otimes I_{2}))\times(|0\rangle\otimes(|0\rangle\otimes|0\rangle))\\ =&(H\times|0\rangle)\otimes((I_{2}\times|0\rangle)\otimes(I_{2}\times|0\rangle))\\ =&|+\rangle\otimes(|0\rangle\otimes|0\rangle)\\ \end{array}$$

Correspondingly, the expression in (4) turns into

$$\begin{array}[]{rl}&(I_{2}\otimes CX)\times((CX\otimes I_{2})\times(|+\rangle\otimes(|0\rangle\otimes|0\rangle))).\end{array}$$

The inner layer matrix multiplication to be calculated next is the following expression: $(CX\otimes I_{2})\times(|+\rangle\otimes(|0\rangle\otimes|0\rangle))$. Different from the calculation of the previous layer, we have to expand the multiple-qubit quantum gates first. Here is the CX gate. After this step, we use the distributivity of tensor product and matrix product over addition to rewriting it. We also exploit the law in L13 as well as the equations in G_db and B_db. So the inference goes as follows.

$$\begin{array}[]{rl}&(CX\otimes I_{2})\times(|+\rangle\otimes(|0\rangle\otimes|0\rangle))\\ =&((B_{0}\otimes I_{2}+B_{3}\otimes X)\otimes I_{2})\times(|+\rangle\otimes(|0\rangle\otimes|0\rangle))\\ =&((B_{0}\otimes I_{2}\otimes I_{2}+B_{3}\otimes X\otimes I_{2})\times(|+\rangle\otimes(|0\rangle\otimes|0\rangle))\\ =&((B_{0}\otimes(I_{2}\otimes I_{2}))\times(|+\rangle\otimes(|0\rangle\otimes|0\rangle))+((B_{3}\otimes(X\otimes I_{2})\times(|+\rangle\otimes(|0\rangle\otimes|0\rangle))\\ =&((B_{0}\times|+\rangle)\otimes(I_{2}\times|0\rangle)\otimes(I_{2}\times|0\rangle))+((B_{3}\times|+\rangle)\otimes(X\times|0\rangle)\otimes(I_{2}\times|0\rangle))\\ =&\frac{1}{\sqrt{2}}\cdot(|0\rangle\otimes(|0\rangle\otimes|0\rangle))+\frac{1}{\sqrt{2}}\cdot(|1\rangle\otimes(|1\rangle\otimes|0\rangle))\end{array}$$

Other laws such as the associativity and distributivity of scalar multiplication may also be used in the inference, though they are not demonstrated in this example.

We continue in this way until no more matrix multiplication is possible. It is worth noting that if there is a zero matrix in the summands of a certain layer of calculation, it can be eliminated immediately, without being carried into the next layer of calculation. Finally, we simplify the expressions about complex numbers.

The above steps appear a bit complex, but they can be fully automated in Coq, fortunately. The script for implementing the strategy operate$\_$reduce is as follows.

Ltac inner_reduce :=
  unfold_operator;
  kron_plus_distr;
  isolate_scale;
  assoc_right;
  try repeat rewrite Mmult_plus_distr_l;
  try repeat rewrite Mmult_plus_distr_r;
  repeat rewrite <- Mscale_kron_dist_r;
  repeat mult_kron;
  repeat rewrite Mscale_mult_dist_r;
  repeat (autorewrite with G_db;
  repeat cancel_0;
  repeat rewrite Mscale_kron_dist_r);
  repeat rewrite <- Mmult_plus_distr_l.

Ltac operate_reduce :=
  assoc_right;
  repeat inner_reduce;
  reduce_scale;
  unified_base.

In summary, using the strategy operate$\_$reduce, we can formally prove the equation in (3) automatically.

As we can see, our general framework is to formalize circuits symbolically as terms, and simplify terms involving matrices by (semi-)automated term rewriting instead of actually computing the matrices. The strategies in Sections 3.2-3.5 are designed with the common goal: to reduce matrix multiplications in the form of $\braket{i}{j}$ into scalars, and simplify the matrix representation by absorbing ones and eliminating zeros. This symbolic approach of reasoning about circuits turns out to be effective; see Section 6.7 for more detailed discussion.

Although it is very intuitive to represent pure quantum states by vectors, there is an inconvenience. In quantum mechanics, the global phase of a qubit is often ignored. For example, we would not distinguish $|\psi\rangle$ from $e^{i\theta}|\psi\rangle$ for any $\theta$. However, when written in vector form, $|\psi\rangle$ and $e^{i\theta}|\psi\rangle$ may be different because of the coefficient $e^{i\theta}$ present in the latter but not in the former. Therefore, we use the symbol $\approx$ to denote such an equivalence, i.e. $e^{i\theta}|\psi\rangle\approx|\psi\rangle$. As a matter of fact, we can be more general and define an observational equivalence for matrices, as given below.

Definition obs_equiv {m n : nat} (A B : Matrix m n) : Prop :=
   exists c : C, Cmod c = R1 /$\backslash$  c .* A = B.

Infix "≈" := obs_equiv.

In the above definition, the condition Cmod c = R1 means that the norm of the complex number c is 1 and c .* A = B says that the matrix A is equal to B after a scalar product with the coefficient c. See Section 6.1 for more concrete examples that use the relation $\approx$.

Note that if quantum states are represented by density matrices, we have

$$(e^{i\theta}|\psi\rangle)(e^{i\theta}|\psi\rangle)^{\dagger}~{}=~{}(e^{i\theta}|\psi\rangle)(e^{-i\theta}\langle\psi|)~{}=~{}|\psi\rangle\langle\psi|.$$

Therefore, the discrepancy entailed by the global phase disappears: the two vectors $e^{i\theta}|\psi\rangle$ and $|\psi\rangle$ correspond to the same density matrix $|\psi\rangle\langle\psi|$. Representing states by density matrices can thus omit unnecessary details and in some cases simplify our reasoning. This is a small but useful trick in formal verification of quantum circuits, which does not seem to have been exploited in the literature. In Sections 6.2 and 6.3, we give two examples where the input and output quantum states of the circuits are given in terms of density matrices.

If the state of a quantum system is represented by a density matrix, the reasoning strategies discussed above can still be used. For instance, suppose a system is in the initial state given by density matrix $\rho$. After the execution of a quantum circuit implementing some unitary transformation $U$, the system changes into the new state $\rho^{\prime}=U\rho U^{\dagger}$. Let $\rho=\sum_{j}\lambda_{j}|j\rangle\langle j|$ be its spectral decomposition, where $\lambda_{j}$ are eigenvalues of $\rho$ and the vectors $|j\rangle$ the corresponding eigenvectors. It follows that

$$\rho^{\prime}~{}=~{}U(\sum_{j}\lambda_{j}|j\rangle\langle j|)U^{\dagger}~{}=~{}\sum_{j}\lambda_{j}U|j\rangle(U|j\rangle)^{\dagger}\ .$$

(5)

Therefore, we can first simplify $U|j\rangle$ into a vector, take its dual and then obtain $\rho^{\prime}$ easily. Our approach to symbolic reasoning also applies in this setting.

We define two functions density and super in advance. The former converts states in the vector form into corresponding states in the density matrix form. The latter formalizes the transformation process between states in the density matrix form.

    
Definition density {n} ($\psi$ : Matrix n 1) : Matrix n n := $\psi$ × $\psi$†.
Definition super {m n} (M : Matrix m n) : Matrix n n -> Matrix m m :=
   fun $\rho$ => M × $\rho$ × M†.
        

We introduce the simplification strategy called super$\_$reduce for states in the density matrix form.

Ltac super_reduce:=
  unfold super,density;
(* Expand super and density *)
  match goal with
(* Match the pattern of target
   with U × $\psi$ × $\psi$† × U† *)
   ||-context [ (?A × ?B) × ?A† ] =>
     match B with
      | (?C × ?C†) =>
         transitivity ((A × C) × (C† × A†)
(* Cast uniform types *)
     end
  end;
  [repeat rewrite <- Mmult_assoc; reflexivity|..];
  rewrite <- Mmult_adjoint;
(* Extract adjoint *)
  let Hs := fresh "Hs" in
    match goal with
    ||-context [ (?A × ?B) × ?C† ) = ?D × ?D†]=>
      match C with
      | ?A × ?B=> assert (A × B = D) as Hs
      end
    end;
(* Use operate_reduce to prove vector states
   and rewrite it in density matrix form *)
  [try reflexivity; try operate_reduce |
   repeat rewrite Hs; reflexivity].

In the above strategy, we first expand the density and super functions in the target. Next, we match the pattern of the target to see whether it is in the form $U\times|\psi\rangle\times\langle\psi|\times U^{\dagger}$ (the middle of the equation in (5)) and cast uniform types, for the reasons to be discussed in Section 5. Then we exploit the law in L14 to extract adjoint of multiplication terms so the target becomes $U\times|\psi\rangle\times(U\times|\psi\rangle)^{\dagger}$ as in the right hand side of the equation in (5). Finally, we use the strategy operate$\_$reduce to conduct the proof for states in vector form and rewrite it back in density matrix form. Note that we omit dimensions for presentation purpose, in practice we need to specify these implicit arguments.

As mentioned in Section 2.2, the mixed state $\rho=\sum_{i}p_{i}|\psi_{i}\rangle\langle\psi_{i}|$ is an ensemble $\{(p_{i},|\psi_{i}\rangle)\}$ of pure states $|\psi_{i}\rangle$, where $p_{i}$ is the probability of the pure state $|\psi_{i}\rangle$ with $\sum_{i}p_{i}=1$. In Coq, we use a list of pairs of real numbers and density operators to define mixed states.

    
Definition Pure n := (R * (Matrix n n)).
Definition Mix n := (list (Pure n)).

Recall that the application of unitary operator $U$ on the mixed state $\rho$ is in the following form:

$$\rho=\sum_{i}p_{i}|\psi_{i}\rangle\langle\psi_{i}|~{}\xrightarrow{U}~{}U\rho U^{\dagger}=\sum_{i}p_{j}U|\psi_{i}\rangle\langle\psi_{i}|U^{\dagger}.$$

(6)

We formalize (6) in Coq as follows.

Fixpoint UnitMix {n} (A : Matrix n n) (m : Mix n): Mix n :=
  match m with
  | [] => []
  | a :: b => (match a with
              |(x , y) => (x , super A y)
              end) :: (UnitMix A b)
  end.
    

After a measurement, results can be expressed with mixed states, as introduced in Section 2.2. For simplicity, we use projection measurements, so for each measurement operator $M_{m}$ we have $M_{m}^{{\dagger}}M_{m}=M_{m}$. We first define the measurement operators of any dimension as follows, where the two parameters $n$ and $k$ stand for the space dimension and the position of an active qubit.

    
Definition Mea0 (n k : N) := (I (2^k) ⊗ ∣0⟩⟨0∣ ⊗ I (2^(n-k))).
Definition Mea1 (n k : N) := (I (2^k) ⊗ ∣1⟩⟨1∣ ⊗ I (2^(n-k))).
Definition Mea (n k : N) := Mea0 n k .+ Mea1 n k.
    

Then we formalize measurements on the mixed state in Coq as follows. Note that a pure state $\rho$ can be regarded as a special mixed state $[(1,\rho)]$.

    
Fixpoint MeaMix {n} (m k : N) (l : Mix n) : Mix n :=
  match l with
  | [] => []
  | a :: b => match a with
              | (x , y) =>
  [((x * (trace((Mea0 m k) × y))), /(trace ((Mea0 m k)× y)).* super(Mea0 m k) y);
  ((x * (trace((Mea1 m k) × y))), /(trace ((Mea1 m k)× y)).* super(Mea1 m k) y)]
              end ++ (MeaMix m k b)
  end.
  

A natural way of interpreting a quantum circuit without measurements is to consider each quantum gate as a unitary matrix and the whole circuit as a composition of matrices that eventually reduces to a single matrix. From this viewpoint, two circuits are equivalent if they denote the same unitary matrix, that is, matrix equivalence $=$ suffices to stand for circuit equivalence.

Directly showing that two matrices are equivalent requires to inspect their elements and compare them component-wisely. Instead, we can take a functional view of matrix equivalence. Let $A,B$ be two $2^{m}\times 2^{m}$ matrices, then $A=B$ if and only if $A|v\rangle=B|v\rangle$ for any vector $|v\rangle\in\mathcal{H}$, where $\mathcal{H}$ is the space of all $m$-qubit states.

Lemma MatrixEquiv_spec: forall {n} (A B: Matrix n n),
   A = B <-> (forall v: Vector n, A × v = B × v).

At first sight, it appears difficult to verify whether $A|v\rangle=B|v\rangle$ for all vectors $|v\rangle$, since there are infinitely many vectors in the state space $\mathcal{H}$. However, both matrices $A$ and $B$ represent linear operators, which means that it suffices to consider the vectors in an orthonormal basis of $\mathcal{H}$, where there are only $2^{m}$ vectors.

In Figure 2 we list some laws that are often useful in simplifying circuits before showing that they are equivalent. Let us verify the validity of the laws. Take the first one as an example. Its validity is stated in Lemma unit$\_$X. In order to prove that lemma, we apply MatrixEquiv$\_$spec and reduce it to Lemma unit$\_$X’, which can be easily proved by the strategy operate$\_$reduce.

Lemma unit_X : X × X = I_2.

Lemma unit_X’ : forall v : Vector 2,  X × X × v = I_2 × v. 

In the right column of Figure 2, the subscripts of $X,Y,Z$ and $H$ indicate on which qubits the quantum gates are applied. For example, $X_{2}$ means that the Pauli-X gate is applied on the second qubit. Thus, the operation $Y_{1}\times X_{2}$ actually stands for $(Y\otimes I_{2})\times(I_{2}\otimes X)$.

In Figure 3, we display some equivalent circuits. In diagram (a), on the right of = is a schematic specification of swapping two qubits, which is implemented by the circuit on the left. Mathematically, the equality is described by Lemma Eq15 below.

    
Definition SWAP :=  B0 ⊗ B0 .+ B1 ⊗ B2 .+ B2 ⊗ B1 .+ B3 ⊗ B3.
Definition XC :=  X ⊗ B3 .+ I_2 ⊗ B0.

Lemma Eq15 : SWAP =  CX × XC × CX.

In diagram (b), there is a controlled operation performed on the second qubit, conditioned on the first qubit being set to zero. It is equivalent to a $CX$ gate enclosed by two Pauli-X gates on the first qubit. The equality is specified by Lemma Eq16 below.

 
    
Definition not_CX := B0 ⊗ X .+ B3 ⊗ I_2.

Lemma Eq16 : not_CX = (X ⊗ I_2) × CX × (X ⊗ I_2).
    

In diagram (c), the controlled phase shift gate on the left is equivalent to a circuit for two qubits on the right. Lemma Eq17 gives a description of this equality.

 
    
Definition CE (u: R) := B0 ⊗ I_2 .+ B3 ⊗ (Cexp u .* B0 .+ Cexp u .* B3).

Lemma Eq17 : CE u = (B0 .+ Cexp u .* B3) ⊗ I_2.

In diagram (d), a controlled gate with two targets is equivalent to the concatenation of two $CX$ gates. This is formalized by Lemma Eq18 below.

 
    
Definition CXX := B0 ⊗ I_2 ⊗ I_2 .+ B3 ⊗ X ⊗ X.
Definition CIX := B0 ⊗ I_2 ⊗ I_2 .+ B3 ⊗ I_2 ⊗ X.

Lemma Eq18 : CXX = CIX × (CX ⊗ I_2).

The previous four lemmas can all be proved by using the strategy operate$\_$reduce in conjunction with MatrixEquiv$\_$spec.

In Section 3 we have formalized the preparation of the 3-qubit GHZ state (cf. Figure 1). Now let us have a look at the Bell states. Depending on the input states, the circuit in Figure 4 gives four possible output states. The correctness of the circuit is validated by the four lemmas below, where the states are given in terms of density matrices and the circuit is described by a super-operator. It is easy to prove them by using our strategy super$\_$reduce.

Definition bl00 := /√2 .* (∣0,0⟩) .+ /√2 .* (∣1,1⟩).
Definition bl01 := /√2 .* (∣0,1⟩) .+ /√2 .* (∣1,0⟩).
Definition bl10 := /√2 .* (∣0,0⟩) .+ (-/√2) .* (∣1,1⟩).
Definition bl11 := /√2 .* (∣0,1⟩) .+ (-/√2) .* (∣1,0⟩).

Lemma pb00 : super (CX × (H ⊗ I_2)) (density ∣0,0⟩) = density bl00.
Lemma pb01 : super (CX × (H ⊗ I_2)) (density ∣0,1⟩) = density bl01.
Lemma pb10 : super (CX × (H ⊗ I_2)) (density ∣1,0⟩) = density bl10.
Lemma pb11 : super (CX × (H ⊗ I_2)) (density ∣1,1⟩) = density bl11.

An alternative way of interpreting a quantum circuit without measurements is to consider it as an operator that changes input quantum states to output states and abstracts away unobservable details. Therefore, we reuse the notion of matrix equivalence $\approx$ introduced in Section 3 and call it observational equivalence for quantum circuits.

The rationale of using $\approx$ is that from an observational point of view, global phases are irrelevant to the observed properties of the physical system under consideration and can be ignored as far as quantum states are concerned.

The following lemma provides a functional view of observational equivalence. It is a counterpart of MatrixEquiv$\_$spec given in Section 4.1. Let $A,B$ be two operators, we have $A\approx B$ if and only if $A|\psi\rangle\approx B|\psi\rangle$ for any state $|\psi\rangle$.

Lemma ObsEquiv_operator: forall {n} (A B: Matrix n n),
   A ≈ B <-> (forall $\psi$: Matrix n 1, A × $\psi$ ≈  B × $\psi$).

Furthermore, two states are equal modulo a global phase, i.e. $|\psi\rangle\approx|\phi\rangle$ if and only if their density matrices are exactly the same, i.e. $|\psi\rangle\langle\psi|=|\phi\rangle\langle\phi|$. We formally prove this property as it motivated us to introduce density matrices to represent quantum states in Section 3.6.

Lemma ObsEquiv_state: forall {n} ($\psi$ $\phi$: Matrix n 1),
   $\psi$ ≈ $\phi$ <-> $\psi$ × ($\psi$†) = $\phi$ × ($\phi$†) .

Although both matrix equivalence $=$ and observational equivalence $\approx$ can be used for relating circuits, the former is strictly finer than the latter in the sense that $A=B$ implies $A\approx B$ but not the other way around. Therefore, in the rest of the paper, we relate circuits by $=$ whenever possible, as they are also related by $\approx$.

Moreover, it is not difficult to see that $=$ is a congruence relation. For example, if $A,B$ are two quantum gates and $A=B$, then we can add a control qubit to form controlled-$A$ and controlled-$B$ gates, which are still identified by $=$. However, the relation $\approx$ does not satisfy such kind of congruence property. To see this, note that $I\approx-I$. But the controlled-$I$ gate is quite different from controlled-($-I$): the latter transforms $\frac{1}{\sqrt{2}}(|00\rangle+|11\rangle)$ into $\frac{1}{\sqrt{2}}(|00\rangle-|11\rangle)$ whereas the former keeps it unchanged. In Section 6.1 we will see a concrete example of using $\approx$, where quantum states are identified by purposefully ignoring their global phases.