Joint EigenValue Decomposition for Quantum Information Theory and Processing

Abstract

The interest in quantum information processing has given rise to the development of programming languages and tools that facilitate the design and simulation of quantum circuits. However, since the quantum theory is fundamentally based on linear algebra, these high-level languages partially hide the underlying structure of quantum systems. We show that in certain cases of practical interest, keeping a handle on the matrix representation of the quantum systems is a fruitful approach because it allows the use of powerful tools of linear algebra to better understand their behavior and to better implement simulation programs. We especially focus on the Joint EigenValue Decomposition (JEVD). After giving a theoretical description of this method, which aims at finding a common basis of eigenvectors of a set of matrices, we show how it can easily be implemented on a Matrix-oriented programming language, such as Matlab (or, equivalently, Octave). Then, through two examples taken from the quantum information domain (quantum search based on a quantum walk and quantum coding), we show that JEVD is a powerful tool both for elaborating new theoretical developments and for simulation.

1. Introduction

The field of quantum information is experiencing a resurgence of interest due to the recent implementation of secure transmission systems [1] based on the teleportation of quantum states in metropolitan networks and in the context of satellite transmissions, further underscored by the development of quantum computers. A new path for intercontinental quantum communication opened up in 2017 when a source onboard a Chinese satellite made it possible to distribute entangled photons between two ground stations, separated by more than 1000 km [2, 3]. Experiments using optical fibers [4] and terrestrial free-space channels [5] have also proved that the use of quantum entanglement can be achieved over large distances.

Quantum programming languages, such as Q# [6] have been developed to facilitate the design and simulation of quantum circuits. The underlying quantum theory is quite complex and often counter-intuitive due to the fact that it relies on linear algebra and tensor products—for instance, the state of a set of three independent qubits (quantum bits) is not described by a 3-dimensional vector, as would be the case for classical bits, but by a $2^3$-dimensional vector which lives in a Hilbert space constructed by tensor products of lower-dimensional spaces. Therefore, these programming languages are helpful for people who do not need to bother with the underlying theory.

However, since the quantum theory is fundamentally based on linear algebra, there are cases of practical interest for the researcher in which keeping a handle on the matrix representation of the quantum systems is a fruitful approach because it allows the use of powerful tools of linear algebra to better understand their behavior and to better implement simulation programs.

In this chapter, our objective is to illustrate how the concept of Joint EigenValue Decomposition (JEVD) can provide interesting results in the domain of quantum information. The chapter is organized as follows. In Section 2, we give some mathematical background and in Section 3, we provide basic elements to understand quantum information. Then, in Section 4, we show an example of the application of JEVD to quantum coding, more precisely we propose an algorithm, based on JEVD, to identify a quantum encoder matrix from a collection of given Pauli errors. Finally, in Section 5, we show that JEVD is a powerful tool for the analysis of a quantum walk search. More precisely, we prove that, while the quantum walk operates in a huge state space, there exists a small subspace that captures all the essential elements of the quantum walk, and this subspace can be determined thanks to JEVD.

2. Mathematical background

2.1 Matrices and notations

We note $U^T$ the transpose of a matrix $U$ and $U^{\ast }$ the transpose conjugate of $U$.

$H$ is the normalized Hadamard $2\times 2$ matrix and $H_N$ the $N\times N$ Hadamard matrix obtained by the Kronecker product (defined in the next subsection):

$$H=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}1& 1\\ 1& -1\end{array}\right)\text{and}{H}_N=H^{\otimes n}\left(N=2^n\right)$$E1

$I_N$ is the $N\times N$ identity matrix (which will sometimes be noted $I$ when its dimension is implicit).

In the domain of quantum information processing, we mainly have unitary matrices. A square matrix $U$ is unitary [7] if $U^{\ast }U={\mathit{UU}}^{\ast }=I$. The columns of a unitary matrix are orthonormal and its eigenvalues are of norm 1. If the unitary matrix is real, its eigenvalues come by conjugate pairs.

We call “shuffle matrix” the permutation matrix $P_{a,b}$ which represents the permutation obtained when one writes elements row by row in an $a\times b$ matrix and reads them column by column. For instance, set $a=2$ and $b=3$. If one writes the elements $\mathrm{1,2,3,4,5,6}$ row by row in a $2\times 3$ matrix and reads them column by column, the order becomes $\mathrm{1,4,2,5,3,6}$. Then the shuffle matrix is the permutation matrix such that $\left(142536\right)=\left(\begin{array}{cccccc}1& 2& 3& 4& 5& 6\end{array}\right)P_{2,3}$. The inverse of $P_{a,b}$ is $P_{b,a}={\left(P_{a,b}\right)}^T$.

$G_n$ is the $n\times n$ Grover diffusion matrix defined by [8]:

$$G_n=-I_n+2{\theta }_n{\theta }_n^T$$E2

where ${\theta }_n$ the $\left(n\times 1\right)$ vector is defined by ${\theta }_n={\left[\begin{array}{cccc}1& 1& \cdots & 1\end{array}\right]}^T/\sqrt{n}$. It is easy to see that $G_n{\theta }_n={\theta }_n$. Therefore, ${\theta }_n$ is an eigenvector of $G_n$ with eigenvalue $+1$. We can also see that for any vector $v$ orthogonal to ${\theta }_n$ we have $G_{n}v=-v$. It follows that $G_n$ has two eigenvalues, $-1$ and $+1$, and the dimensions of the associated eigenspaces are $n-1$ and $1$.

2.2 Kronecker product

The Kronecker product, denoted by $\otimes$, is a bilinear operation on two matrices. If $A$ is a $k\times l$ matrix and $B$ is a $m\times n$ matrix, then the Kronecker product is the $\mathit{km}\times \mathit{ln}$ block matrix $C$ below:

$$C=A\otimes B=\left(\begin{array}{ccc}{a}_{11}B& \cdots & a_{1l}B\\ ⋮& \ddots & ⋮\\ a_{k1}B& \cdots & a_{\mathit{kl}}B\end{array}\right)$$E3

Assuming the sizes are such that one can form the matrix products $\mathit{AC}$ and $\mathit{BD}$, an interesting property, known as the mixed-product property, is:

$$\left(A\otimes B\right)\left(C\otimes D\right)=\left(\mathit{AC}\right)\otimes \left(\mathit{BD}\right)$$E4

The Kronecker product is associative, but not commutative. However, there exist permutation matrices (the shuffle matrices defined in the previous subsection) such that, if $A$ is an $a\times a$ square matrix and $B$ a $b\times b$ square matrix, then [9]:

$$\left(A\otimes B\right)P_{a,b}=\left(B\otimes A\right)P_{b,a}$$E5

2.3 Singular value decomposition, image, and kernel

The Singular Value Decomposition (SVD) of an $m\times n$ matrix $A$ is [7]:

$$A={\mathit{USV}}^{\ast }$$E6

where $U$ and $V$ are unitary matrices, and $S$ is diagonal. The diagonal of $S$ contains the Singular Values, which are real nonnegative numbers, ranked by decreasing order. The sizes of the matrices are $U\left(m\times m\right)$, $S\left(m\times n\right)$ and $V\left(n\times n\right)$. The SVD is a very useful linear algebra tool because it reveals a great deal about the structure of a matrix.

The image and the kernel of $A$ are defined by:

$$\mathit{Im}\left(A\right)=\left\{y\in {\mathbb{C}}^m:y=\mathit{Ax}\text{for\hspace{0.17em}some}x\in {\mathbb{C}}^n\right\}$$E7

$$ker\left(A\right)=\left\{x\in {\mathbb{C}}^n:\mathit{Ax}=0\right\}$$E8

When used in an algorithm, the notation $\mathit{\text{null}}$ will also be used for a procedure that computes a matrix whose columns are an orthonormal basis of the kernel of $A$.

The complement of a subspace $\mathcal{A}$ within a vector space $\mathcal{H}$ is defined by:

$${\left(\mathcal{A}\right)}^c=\left\{y\in \mathcal{H}:x^{\ast }y=0\text{for}\mathrm{all}x\in \mathcal{A}\right\}$$E9

In an algorithm, if the columns of $A$ are an orthonormal basis of $\mathcal{A}$ then the columns of $B=\mathit{\text{null}}\left(A^{\ast }\right)$ provide an orthonormal basis of ${\left(\mathcal{A}\right)}^c$.

The rank of $A$ is its number of nonzero singular values. When programmed on a computer determination of the rank must take into account finite precision arithmetic, which means that “zero” is replaced by “extremely small” (less than a given tolerance value). Let us note $r=\mathit{rank}\left(A\right)$. We have

$$dim\left(\mathit{Im}\left(A\right)\right)=r$$E10

$$dim\left(ker\left(A\right)\right)=n-r$$E11

An orthonormal basis of $ker\left(A\right)$ is obtained by taking the last $n-r$ columns of the matrix $V$.

2.4 Joint eigenspaces and joint eigenvalue decomposition (JEVD)

The eigenvalue decomposition of a unitary matrix $A$ is:

$$A={\mathit{VDV}}^{\ast }$$E12

where $D$ is a diagonal matrix, the diagonal of which contains the eigenvalues, and $V$ is a unitary matrix whose columns are the eigenvectors.

Let us note $E_{\lambda }^A$ the eigenspace of an operator $A$ associated with an eigenvalue $\lambda .$ The joint eigenspace $E_{\lambda ,\mu }^{A,B}$ is:

$$E_{\lambda ,\mu }^{A,B}=E_{\lambda }^A\cap E_{\mu }^B$$E13

A property of great interest in quantum information processing is that within $E_{\lambda ,\mu }^{A,B}$ (and even within any union of joint eigenspaces) the operators $A$ and $B$ commute.

Determination of the joint eigenspace on a computer may be determined through the complement, because:

$$E_{\lambda }^A\cap E_{\mu }^B={\left({\left(E_{\lambda }^A\right)}^c\cup {\left(E_{\mu }^B\right)}^c\right)}^c$$E14

Using Matrix-oriented programming languages, such as Matlab or Octave, this requires only a few lines. Let us note $A_{\lambda }$ and $B_{\mu }$ matrices whose columns are orthonormal bases of $E_{\lambda }^A$ and $E_{\mu }^B$ and $\left[.\right]$ the horizontal concatenation of matrices. The following computation procedure provides a matrix $C$ whose columns are an orthonormal basis of $E_{\lambda ,\mu }^{A,B}$:

$$C=\mathit{\text{null}}\left(\left[\begin{array}{cc}\mathit{\text{null}}\left(A_{\lambda }\right)& \mathit{\text{null}}\left(B_{\mu }\right)\end{array}\right]\right)$$E15

However, it is not efficient in terms of complexity and in the next sections we will propose faster computational procedures, adapted to each context.

A lower bound on the dimension of a joint eigenspace can be obtained as follows. Let us note $n$ the dimension of the full space. We have, obviously:

$$dim{E}_{\lambda }^A\cup E_{\mu }^B\le n$$E16

and we know that:

$$dim{E}_{\lambda }^A\cup E_{\mu }^B=dim{E}_{\lambda }^A+dim{E}_{\mu }^B-dim{E}_{\lambda }^A\cap E_{\mu }^B$$E17

Combining both equations, we obtain:

$$dim{E}_{\lambda ,\mu }^{A,B}\ge dim{E}_{\lambda }^A+dim{E}_{\mu }^B-n$$E18

3. Quantum information principles

A quantum system is described by a state vector $\left|\psi \right.⟩\in {\mathbb{C}}^N$, where $N$ is the dimension of the system. Since in the quantum formalism states $\left|\psi \right.⟩$ and $\gamma \left|\psi \right.⟩$ are equivalent, for any nonzero complex number $\gamma$, the state is usually represented by a normed vector and the global phase is considered irrelevant.

As long as it remains isolated, the evolution of a quantum system is driven by the Schrödinger equation. The latter is a first-order differential equation operating on the quantum state. Its integration shows that the quantum states at times $t_1$ and $t_2$ are linked by a unitary matrix $U$ such that $\left|{\psi }_2\right.⟩=U\left|{\psi }_1\right.⟩$. The norm is preserved because $U$ is unitary.

The second kind of evolution, called “measurement,” may occur if the system interacts with its environment. A measurement consists of the projection of the state onto a subspace of ${\mathbb{C}}^N$. When the measurement is controlled, it consists in defining a priori a decomposition of the state space into a direct sum of orthogonal subspaces $\underset{i}{\oplus }{\mathcal{H}}_i$. The measurement randomly selects one subspace. The result of the measurement is an identifier of the selected subspace (for instance, its index $i$). After measurement, the state is projected onto ${\mathcal{H}}_i$. If $P_i$ is the projection matrix onto ${\mathcal{H}}_i$, then the state becomes $P_i\left|\psi \right.⟩$ (which is then renormalized because the projection does not preserve the norm). The probability of ${\mathcal{H}}_i$ being selected is the square norm of $P_i\left|\psi \right.⟩$.

It is worth noting that a measurement may destroy a part of quantum information (because usually, a projection is an irreversible process), while the unitary evolution is reversible, and as such, preserves quantum information. Consequently, measurements must be used with extreme caution—how to design the system and the measurement device to measure only what is strictly required and not more is one of the difficult problems encountered in quantum information processing.

Quantum systems of special interest for quantum information processing are qubits (quantum bits) and qubit registers. A qubit belongs to a 2D quantum system with state a normed vector of ${\mathbb{C}}^2$. To highlight links with classical digital computation, it is convenient to note $\left|0\right.⟩$ and $\left|1\right.⟩$ for the orthonormal basis of ${\mathbb{C}}^2$. Physically any 2D quantum system can carry a quantum bit. For instance, the spin of an electron is a 2D quantum system, and the spins up and down can be associated with the basic states $\left|0\right.⟩$ and $\left|1\right.⟩$. A general qubit has an expression:

$$\left|\psi \right.⟩={\alpha }_0\left|0\right.⟩+{\alpha }_1\left|1\right.⟩$$E19

where ${\alpha }_0$ and ${\alpha }_1$ are complex numbers subject to ${\left|{\alpha }_0\right|}^2+{\left|{\alpha }_1\right|}^2=1$.

A qubit register is a $2^n$-D quantum system which, for convenience, is usually referred to as a standard orthonormal basis noted $\left\{\left|\mathrm{0\dots 00}\right.⟩,\left|\mathrm{0\dots 01}\right.⟩,\left|\mathrm{0\dots 10}\right.⟩,\dots ,\left|\mathrm{1\dots 11}\right.⟩\right\}$ and then, by analogy with classical digital processing, $n$ is the number of qubits. For instance, for $n=2$ the basis is $\left\{\left|00\right.⟩,\left|01\right.⟩,\left|10\right.⟩,\left|11\right.⟩\right\}$, where $\left|\mathit{ab}\right.⟩=\left|a\right.⟩\otimes \left|b\right.⟩$, and the quantum state of the register is:

$$\left|\psi \right.⟩=\sum _{\left(a,b\right)\in {\left\{0,1\right\}}^2}{\gamma }_{\mathit{ab}}\left|\mathit{ab}\right.⟩$$E20

Note that, contrary to classical digital registers, the qubits are usually not separable, hence the register must be considered as a whole. We say that the qubits are entangled. However in the special case where the coefficients ${\gamma }_{\mathit{ab}}$ can be decomposed in the form ${\gamma }_{\mathit{ab}}={\alpha }_a{\beta }_b$ the state can be written as a tensor product of the states of two qubits, which can be considered separately. Then, we have:

$$\left|\psi \right.⟩=\left({\alpha }_0\left|0\right.⟩+{\alpha }_1\left|1\right.⟩\right)\otimes \left({\beta }_0\left|0\right.⟩+{\beta }_1\left|1\right.⟩\right)$$E21

4. Application of JEVD to quantum coding

4.1 Principle of quantum coding

The objective of quantum coding is to protect quantum information [10]. In the classical domain, the information can be protected using redundancy—for instance, if we want to transmit bit 0 on a noisy communication channel, we can instead transmit 000 (and, similarly, transmit 111 instead of 1). On the receiver side, if one error has occurred on the channel, for instance, if the second bit is false, we receive 010 instead of 000, from which we can still guess that the most probable hypothesis is that the transmitted word was 000. Of course, if there were two errors the transmitted word could have been 111, but it is assumed that the probability of error is low, hence two errors are less likely than one error. More elaborated channel codes have been proposed, but fundamentally they are all based on the idea of adding redundancy and assuming that the probability of channel error is low.

In the quantum domain, it is impossible to use redundancy because it is impossible to copy a quantum state (this is due to the “no-cloning theorem” [11]). However, we can use entanglement to produce the quantum equivalent of classical redundancy. The principle of quantum coding is shown in Figure 1. Assume we want to protect the quantum state $\left|\psi \right.⟩$ of a $k$-qubit register. We add $r$ ancillary qubits initialized to $\left|0\right.⟩$ to form an $n$-qubit register ($n=k+r$). The encoder is represented by a unitary $2^n\times 2^n$ matrix $U$. Then, errors may occur on the encoded state: they are represented by a unitary matrix $E$. The decoder is represented by another unitary matrix $U^{\ast }$ which is the transpose conjugate of the encoding matrix. Finally, we measure the last $r$ qubits of the decoded state, and, depending on the result of the measurement, we apply the appropriate restoration matrix $U_c$ (which is a unitary matrix of size $2^k\times 2^k$) to the $k$-qubit register composed of the first $k$ qubits of the decoded state.

Figure 1.

Principle of quantum coding.

As an illustration, let us consider $n=2$, $k=1$ and the very simple quantum encoder shown in Figure 2. It is a basic quantum circuit known as the CNOT quantum gate, and it is represented by the unitary matrix below:

Figure 2.

CNOT quantum gate.

$$U=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 1\\ 0& 0& 1& 0\end{array}\right)$$E22

A quantum error on a qubit is described by a $2\times 2$ unitary matrix. It is convenient to decompose the error as a linear sum of the identity and the Pauli matrices below [12]:

$$\begin{array}{ccc}Z=\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right)& X=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)& Y=\left(\begin{array}{cc}0& i\\ -i& 0\end{array}\right)\end{array}$$E23

Let us consider that an error may appear on the first encoded qubit and that this error, if present, is represented by the unitary Pauli matrix $X$. Then, the error matrix which acts on the encoded state is:

$$E=X\otimes I=\left(\begin{array}{cccc}0& 1& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& 1\\ 0& 0& 1& 0\end{array}\right)$$E24

It is easy to check that:

$$F=U^{\ast }\mathit{EU}=\left(\begin{array}{cccc}0& 0& 0& 1\\ 0& 0& 1& 0\\ 0& 1& 0& 0\\ 1& 0& 0& 0\end{array}\right)=X\otimes X$$E25

The state at the input of the encoder is $\left[{\alpha }_0{\alpha }_1\right]{}^{\mathrm{T}}\otimes \left[10\right]{}^{\mathrm{T}}=\left[{\alpha }_{0}0{\alpha }_{1}0\right]{}^{\mathrm{T}}$. The state at the output of the decoder is, therefore, $\left[0{\alpha }_{1}0{\alpha }_0\right]{}^{\mathrm{T}}$_.

Measuring the second qubit on the output of the decoder consists in decomposing the state space into a direct sum ${\mathcal{H}}_0\oplus {\mathcal{H}}_1$ of two subspaces spanned by $\left\{\left|00\right.⟩,\left|10\right.⟩\right\}$ and $\left\{\left|01\right.⟩,\left|11\right.⟩\right\}$. The result of the measurement will be either 0 or 1 (index of the selected subspace), and by analogy with classical decoding, this result will be called the “syndrome.” The projections on these subspaces are $\left[00\right]{}^{\mathrm{T}}$ and $\left[{\alpha }_1{\alpha }_0\right]{}^{\mathrm{T}}$. Then the probability to obtain syndrome $1$ is 1.

The measurement then projects the state onto ${\mathcal{H}}_1$. Note that in this particular case, the information is preserved by the projection. Then, applying the operator $U_c=X$ to the projected state restores the initial state.

Similarly, if there is no error, we can see that $F$ is the identity matrix, then the projections on the subspaces are $\left[{\alpha }_0{\alpha }_1\right]{}^{\mathrm{T}}$ and $\left[00\right]{}^{\mathrm{T}}$. In that case, the syndrome is 0 and the state is projected onto ${\mathcal{H}}_0$. Correction is done by applying the operator $I$ to the projected state, which is equivalent to doing nothing.

The very simple code used above, as an illustration, cannot correct more complex errors (for instance, an error $Z$ on the first qubit). However, there exist efficient quantum codes, such as the Steane code [13], and the Shor code [14]. A remarkable result of quantum coding theory is that a linear combination of correctable errors is correctable [15].

Figure 3 shows the Steane Encoder, which is a $\left(n=7,k=1,t=1\right)$ quantum encoder. This means that it encodes $k=1$ qubits on $n=7$ qubits and it is able to correct any error occurring on $t=1$ encoded qubits. It is built with Hadamard (Eq. (1)) and CNOT (Eq. (22)) quantum gates. From this circuit description, it is possible to obtain the coding matrix $U$.

Figure 3.

Steane encoder.

4.2 Determination of encoder matrix using JEVD

The problem we address can be stated as follows (see Figure 1 for the notations)—given a list of $n$ independent Pauli errors $E_i$ with corresponding diagonal outer errors $F_i$, determine the unitary operator $U$ (quantum encoder) such that:

$$U^{\ast }{E}_{i}U=F_i\forall i\in \left\{1,\dots ,n\right\}$$E26

This equation shows that the columns of $U$ are the eigenvectors of $E_i$. Specification of the code by a small set of Pauli errors is very convenient and the interest of automatic determination of matrix $U$ is to allow further simulations of the behavior of the quantum code in various configurations.

To illustrate and validate the approach that will be developed below, let us consider the collection of $n=7$ Pauli errors shown in Table 1. Here, to be able to check the results, this collection has been chosen to correspond to the Steane encoder (Figure 3), while in a standard application of the method, it would be given a priori. The interest is that here we can compute the encoder matrix from the circuit and this will allow us to check that our method produces the correct encoder matrix.

$E_i$	$F_i$
$Z\otimes Z\otimes Z\otimes Z\otimes Z\otimes Z\otimes Z$	$Z\otimes I\otimes I\otimes I\otimes I\otimes I\otimes I$
$X\otimes X\otimes I\otimes I\otimes I\otimes X\otimes X$	$I\otimes Z\otimes I\otimes I\otimes I\otimes I\otimes I$
$X\otimes I\otimes X\otimes I\otimes X\otimes I\otimes X$	$I\otimes I\otimes Z\otimes I\otimes I\otimes I\otimes I$
$X\otimes I\otimes I\otimes X\otimes X\otimes X\otimes I$	$I\otimes I\otimes I\otimes Z\otimes I\otimes I\otimes I$
$Z\otimes Z\otimes I\otimes I\otimes I\otimes Z\otimes Z$	$I\otimes I\otimes I\otimes I\otimes Z\otimes I\otimes I$
$Z\otimes I\otimes Z\otimes I\otimes Z\otimes I\otimes Z$	$I\otimes I\otimes I\otimes I\otimes I\otimes Z\otimes I$
$Z\otimes I\otimes I\otimes Z\otimes Z\otimes Z\otimes I$	$I\otimes I\otimes I\otimes I\otimes I\otimes I\otimes Z$

Table 1.

Collection of Pauli errors.

We use $n$ independent equations in which each $F_i$ is a tensor product of $I$ and $Z$ only (including only one $Z$). Therefore, matrices $F_i$ are diagonal, and their diagonal elements are $+1$ and $-1$ in equal numbers.

Figure 4 shows the diagonals of the matrices $F_i$ (each row corresponding to one diagonal). Values $-1$ and $+1$ are represented, respectively, by black and white dots.

Figure 4.

Diagonals of matrices $F_i$.

Since matrix $U$ does not depend on $i$ in Eq. (26) its columns are joint eigenvectors of the $E_i$. For instance, in the example above, the ${20}^{\mathit{th}}$ column of $U$ is a joint eigenvector of $E_1$, $E_2$, …, $E_7$ associated to eigenvalues $+1$,$+1$,$-1$,$+1$,$+1$,$-1$,$-1$ (see Figure 4). In the general case, the set of $n$ eigenvalues corresponding to the column $c$ of $U$ is easily obtained by taking the binary representation of $c-1$ with the mapping $0\to +1$ and $1\to -1$.

Now, let us consider the determination of column $c$ of $U$. We know that it is a vector spanning a joint eigenspace of the $E_i$ corresponding to a given set of eigenvalues $\left\{{\lambda }_i,i=1,\dots ,n\right\}$. For each $E_i$ let $A_i$ denote the $2^n\times 2^{n-1}$ matrix whose orthonormal columns span the eigenspace associated to ${\lambda }_i$ and $B_i$ the $2^n\times 2^{n-1}$ matrix whose orthonormal columns span the kernel of $A_i$ (which corresponds to the eigenspace associated to $-{\lambda }_i$).

Let ${\mathcal{Y}}_k$ denote the joint eigenspace corresponding to eigenvalues $\left\{{\lambda }_j,j=1,\dots ,k\right\}$ with $k\in \left\{1,\dots ,n\right\}$. We propose Algorithm 1 to efficiently compute the column of $U$. It computes a series of matrices $Y_k$ whose columns are an orthonormal basis of ${\mathcal{Y}}_k$. Obviously, the searched column of $U$ is $Y_n$. For the moment, let us consider that $K\left(c\right)$=1 (the optimal value will be discussed later).

$$\begin{array}{l}\mathbf{\text{if}}K\left(c\right)=1\mathbf{\text{then}}\\ \mid Y_1=A_1\\ \mathbf{end}\\ \mathbf{\text{for}}k=K\left(c\right)+1\text{to}n\mathbf{do}\\ \left|\begin{array}{l}{C}_k=B_k^{\ast }{Y}_{k-1}\\ Z_k=\mathit{\text{null}}\left(C_k\right)\\ Y_k=Y_{k-1}{Z}_k\end{array}\right.\end{array}$$

Algorithm 1: Algorithm for determination of a joint eigenspace.

The sizes of the matrices are decreasing with $k$:

$C_k$: $2^{n-1}\times 2^{n-k+1}$ $Z_k$: $2^{n-k+1}\times 2^{n-k}{Y}_k$: $2^n\times 2^{n-k}$.

The intuitive ideas under the algorithm are the following:

• $C_k=B_k^{\ast }{Y}_{k-1}$: The orthonormal basis of ${\mathcal{Y}}_{k-1}$ is projected on the kernel of $A_k$. The components of the projected vectors are expressed in the orthonormal basis $B_k$ of that kernel. Consequently, $\mathit{Im}\left(C_k\right)$ is the projection of ${\mathcal{Y}}_{k-1}$ on the kernel of $A_k$, expressed in that kernel.

• A matrix $Z_k$ whose columns are an orthonormal basis of the complement of this projection is determined.

• Finally, the components of the basis vectors are restored to the original space by $Y_k=Y_{k-1}{Z}_k$

Let us prove that the matrices $Y_k$ have orthonormal columns. This is obviously the case for $k=1$. Then, by recursion, we have:

$$Y_k^{\ast }{Y}_k=Z_k^{\ast }{Y}_{k-1}^{\ast }{Y}_{k-1}{Z}_k=I$$E27

Now, let us prove, by recurrence, that $\mathit{Im}\left(Y_k\right)={\mathcal{Y}}_k$.

Obviously, this is the case for $k=1$. Assume this is the case for $k-1$. We have:

$$\mathit{Im}\left(Y_k\right)\subset \mathit{Im}\left(Y_{k-1}\right)={\mathcal{Y}}_{k-1}$$E28

We have also:

$$B_k^{\ast }{Y}_k=B_k^{\ast }\left(Y_{k-1}{Z}_k\right)=\left(B_k^{\ast }{Y}_{k-1}\right)Z_k=C_k{Z}_k=0$$E29

Then

$$\mathit{Im}\left(Y_k\right)\subset ker\left(B_k^{\ast }\right)=\mathit{Im}\left(A_k\right)$$E30

From $\mathit{Im}\left(Y_k\right)\subset {\mathcal{Y}}_{k-1}$ and $\mathit{Im}\left(Y_k\right)\subset \mathit{Im}\left(A_k\right)$ we deduce $\mathit{Im}\left(Y_k\right)\subset {\mathcal{Y}}_k$.

Conversely, assume that a vector $x$ belongs to ${\mathcal{Y}}_k$. Because ${\mathcal{Y}}_k\subset {\mathcal{Y}}_{k-1}$ there exists a vector $b$ such that $x=Y_{k-1}b$ and because $x\in \mathit{Im}\left(A_k\right)$ we have also $B_k^{\ast }x=0$

Then

$$B_k^{\ast }{Y}_{k-1}b=0\Rightarrow C_{k}b=0\Rightarrow \exists a:b=Z_{k}a$$

Therefore $x=Y_{k-1}b=Y_{k-1}{Z}_{k}a=Y_{k}a\Rightarrow x\in \mathit{Im}\left(Y_k\right)\Rightarrow {\mathcal{Y}}_k\subset \mathit{Im}\left(Y_k\right)$.

After execution of the algorithm to determine each column of $U$, there remains an indetermination because the joint eigenvectors (i.e., the columns of $U$) are determined up to a phase factor. This has no consequence on the performance of the quantum code. However, if we want to fix this residual indetermination, we proposed a fast and simple procedure in ref. [16]. The procedure requires an additional set of $n$ Pauli errors in which each additional $F_i$ is a tensor product of $I$ and $X$ only. As an example, for the Steane code, we use Table 2.

$E_i$	$F_i$
$X\otimes X\otimes X\otimes X\otimes X\otimes X\otimes X$	$X\otimes I\otimes I\otimes I\otimes I\otimes I\otimes I$
$Z\otimes I\otimes Z\otimes Z\otimes Z\otimes Z\otimes Z$	$I\otimes X\otimes I\otimes I\otimes I\otimes I\otimes I$
$Z\otimes Z\otimes I\otimes Z\otimes Z\otimes Z\otimes Z$	$I\otimes I\otimes X\otimes I\otimes I\otimes I\otimes I$
$Z\otimes Z\otimes Z\otimes I\otimes Z\otimes Z\otimes Z$	$I\otimes I\otimes I\otimes X\otimes I\otimes I\otimes I$
$X\otimes I\otimes I\otimes I\otimes X\otimes I\otimes I$	$I\otimes I\otimes I\otimes I\otimes X\otimes I\otimes I$
$X\otimes I\otimes I\otimes I\otimes I\otimes X\otimes I$	$I\otimes I\otimes I\otimes I\otimes I\otimes X\otimes I$
$X\otimes I\otimes I\otimes I\otimes I\otimes I\otimes X$	$I\otimes I\otimes I\otimes I\otimes I\otimes I\otimes X$

Table 2.

Additional Collection of Pauli errors.

After these remaining differences have been removed, we obtain an estimated matrix $U$ that is equal to the true matrix, up to a global phase (Figure 5). However, this remaining indetermination does not matter because, as said before, the global phase has no significance in quantum physics. Here we have chosen the global phase so that the encoder matrix is real.

Figure 5.

Estimated Matrix $U$ for the Steane encoder.

Figure 5 shows the matrix computed by our method. We have checked that it is equal to the matrix directly computed from the circuit description.

The programmer may speed up the computation by taking into account the fact that when computing columns $c$ of $U$, some matrices $Y_k$ have already been computed for other columns and can be reused. For instance, in Figure 4, we see that the joint eigenvalues corresponding to columns 19 and 20 are the same, except the last one. Then, when computing column 20, we can set $K\left(20\right)=n$ in Algorithm 1 instead of the default value $K\left(20\right)=1$, because the $Y_{n-1}$ for column 20 is the same as for column 19. More generally, Algorithm 2 written in pseudo-Octave code computes the optimal values of the $K\left(c\right)$.

$$\begin{array}{l}K=\left[\begin{array}{cc}1& 1\end{array}\right]\\ \mathbf{\text{for}}k=2\mathit{\text{to}}n\mathbf{do}\\ |K=\text{reshape}\left(\left[K;k\ast \text{ones}\left(1,2{}^{k-1}\right)\right],\left[12^k\right]\right)\\ \mathbf{end}\end{array}$$

Algorithm 2: Algorithm for computation of the optimal values $K\left(c\right)$.

For instance, for $n=3$ the algorithm produces $K=\left[\begin{array}{cccccccc}1& 3& 2& 3& 1& 3& 2& 3\end{array}\right]$.

5. Application of JEVD to quantum walk search

5.1 Principle of quantum walk search

Let us consider a particle that can move on a graph. In the classical world, at the time $t$ this particle is localized at a node of the graph. It can then randomly choose one of the edges linked to this node to reach one of the adjacent nodes at a time $t+1$. The repeated iteration of this process is the concept of classical random walk.

A quantum walk [17] relies also on a graph, but contrary to the classical walk, here the particle may be located at many nodes at the same time and can choose many edges simultaneously. At the time $t$, the state of the particle is then described by a state vector $\left|{\psi }_t\right.⟩$ and the evolution between times $t$ and $t+1$ is given by a unitary matrix $U=\mathit{SC}$ such that $\left|{\psi }_{t+1}\right.⟩=U\left|{\psi }_t\right.⟩$. The unitary matrices $C$ and $S$ represent, respectively, the choice of the edges and the movement to the adjacent nodes.

In the following, we will consider graphs associated with hypercubes [18]. We will note $n$ the dimension of the hypercube and $N=2^n$ the number of nodes. Figure 6 shows the graph corresponding to a hypercube of dimension $n=3$. It is convenient to label the nodes by binary words. In quantum language, these binary words $\kappa$ are also used to label the basis vectors of the so-called position space ${\mathcal{H}}^{\mathcal{S}}$.

Figure 6.

Hypercube for $n=3$.

The quantum state lives in a Hilbert space built by the tensor product of the position space ${\mathcal{H}}^{\mathcal{S}}$ (corresponding to the nodes) and the coin space ${\mathcal{H}}^{\mathcal{C}}$ (corresponding to the possible movements along the edges) $\mathcal{H}={\mathcal{H}}^{\mathcal{S}}\otimes {\mathcal{H}}^{\mathcal{C}}$. The dimensions of these state spaces are $N_e=\mathit{nN}$, $N$ and $n$.

It is usual to define $C$ as [19]:

$$C=I_N\otimes G_n$$E31

where $G_n$ is the $n\times n$ Grover diffusion matrix defined in Section 2. Matrix $C$ obtained for $n=3$ is shown on Figure 7.

Figure 7.

Matrices $C$ (left) and $O$ (right) for $n=3$.

The structure of $S$ is more complex. It is convenient to first define it in ${\mathcal{H}}^{\mathcal{C}}\otimes {\mathcal{H}}^{\mathcal{S}}$ and then to transpose it to $\mathcal{H}$ using the shuffle matrix $P=P_{n,N}$ (defined in subsection 2.2). Then:

$$S=P\widehat{S}{P}^T$$E32

where

$$\widehat{S}=\mathit{diag}\left({\widehat{S}}_1,\dots ,{\widehat{S}}_n\right)\mathit{\text{and}}{\widehat{S}}_d=I^{\otimes \left(n-d\right)}\otimes X\otimes I^{\otimes \left(d-1\right)}$$E33

The last equation just means that, because a movement along direction $d$ corresponds to an inversion of the $d^{\mathit{th}}$ bit in $\kappa$, the shift operator permutes the values associated to nodes that are adjacent along that dimension.

A quantum walk search can be described by repeated application of a unitary evolution operator $Q$, which can be written:

$$Q=\mathit{UO}$$E34

Here $O$ is the oracle, which aims at marking the solutions. An example of oracle structure is shown in Figure 7. It is a block-diagonal matrix, whose blocks are $-G_n$ when they correspond to a solution and $I_n$ otherwise. Denote $M$ the number of solutions and assume that $M\ll N$ (otherwise the quantum walk search would serve no purpose because the probability of rapidly finding a solution with a classical search would be high). In the example shown in the figure, there are $M=2$ solutions (located at positions 1 and 4 on Figure 6).

Let $t$ denote the number of iterations until a measurement is performed. Starting from an initial state $\left|{\psi }_0\right.⟩$, repeated iterations lead to the state $\left|{\psi }_t\right.⟩=Q^t\left|{\psi }_0\right.⟩$ which is then measured. The theory of quantum walk search [19] shows that the probability of success (that is the probability of obtaining a solution by measurement) oscillates as a function of $t$. This means that theoretical tools which help to understand and simulate quantum walk search lead to the development of methods to determine the optimal time of measurement.

In the sequel, we will show that JEVD is a fruitful tool in this context. Indeed, set $E$ to be the union of the joint eigenspaces of $U$ and $O$, and $\overline{E}$ its complement. Inside $E$_, the operators commute. So, if we note with index $E$ the restrictions of the operators to $E$, we have:

$$Q_E^2=\left(U_E{O}_E\right)\left(U_E{O}_E\right)=U_E{O}_E^2{U}_E=U_E^2$$E35

Then, inside $E$, there is no significant difference between the effective quantum walk $Q$ and the uniform quantum walk $U$, because, after each pair of successive iterations, the evolution is identical. Since the uniform quantum walk has no reason to converge to a solution, we deduce that the interesting part of the process lives in the complement of $E$, that is in $\overline{E}$.

After establishing results about the dimensions of the eigenspaces of $U$ and $O$, we will show that the concept of joint eigenspaces allows us to establish an upper bound on the dimension of the complement, with the remarkable result that this dimension grows only linearly with $n$. Then, we propose an algorithm for efficient computation of the joint eigenspaces and, finally, use it to check our theoretical upper bound.

5.2 Eigenspaces of $U$ and $O$

Set

$$F=H_N\otimes I_n$$E36

Then matrix $F$ diagonalizes $S$:

$$\mathit{FSF}=\left(H_N\otimes I_n\right)P_{N,n}\widehat{S}{P}_{n,N}\left(H_N\otimes I_n\right)$$E37

$$=P_{n,N}\left(I_n\otimes H_N\right)\widehat{S}\left(I_n\otimes H_N\right)P_{N,n}$$E38

$$=P^T\mathit{diag}\left(\dots ,H_N{\widehat{S}}_d{H}_N,\dots \right)P$$E39

The latter term is diagonal because the mixed product property, $H^2=I$ and $\mathit{HXH}=Z$, shows that:

$$H_N{\widehat{S}}_d{H}_N=I^{\otimes \left(n-d\right)}\otimes Z\otimes I^{\otimes \left(d-1\right)}$$E40

Once more, using the mixed product property, we can also prove that $F$ keeps $C$ unchanged, that is:

$$\mathit{FCF}=C$$E41

The diagonal of $\mathit{FSF}$ is the concatenation of the binary representation of the numbers 0 to $N-1$ with the mapping $0\to \left(+1\right)$ and $1\to \left(-1\right)$. That is:

$$\mathit{FSF}=\mathit{diag}\left(S_0,\dots ,S_{\kappa },\dots ,S_{N-1}\right)$$E42

Note that the diagonal of $S_{\kappa }$ contains $k$ times $-1$ and $n-k$ times $+1$ (where $k$ is the Hamming weight of $\kappa$).

Then, because $F^2=I$, $\mathit{FUF}$ is a block diagonal matrix:

$$\mathit{FUF}=\left(\mathit{FSF}\right)\left(\mathit{FCF}\right)$$E43

$$=\mathit{diag}\left(\dots ,S_{\kappa },\dots \right)C$$E44

Block $\kappa$ is then

$$U_{\kappa }=S_{\kappa }{G}_n$$E45

We have:

$$dim{E}_{-}^{U_{\kappa }}\ge dim{E}_{+,-}^{S_{\kappa },G_n}$$E46

$$\ge dim{E}_{+}^{S_{\kappa }}+dim{E}_{-}^{G_n}-n$$E47

$$\ge \left(n-k\right)+\left(n-1\right)-n$$E48

$$\ge n-k-1$$E49

and

$$dim{E}_{+}^{U_{\kappa }}\ge dim{E}_{-,-}^{S_{\kappa },G_n}$$E50

$$\ge dim{E}_{-}^{S_{\kappa }}+dim{E}_{-}^{G_n}-n$$E51

$$\ge k+\left(n-1\right)-n$$E52

$$\ge k-1$$E53

Then, there is only room left for at most 2 eigenvalues, specifically, at most a pair of conjugate ones.

Assume that this pair of eigenvalues exists. Since the diagonal of $G_n$ contains $-1+\frac{2}{n}$, the trace of $U_{\kappa }$ is:

$$\mathit{\text{trace}}\left(U_{\kappa }\right)=\left(\left(n-k\right)\left(+1\right)+k\left(-1\right)\right)\left(-1+\frac{2}{n}\right)$$E54

$$=-n+2k+2\left(1-2\frac{k}{n}\right)$$E55

The sum of the eigenvalues is equal to the trace and we already have eigenvalue $-1$ with multiplicity $n-k-1$ and eigenvalue $+1$ with multiplicity $k-1$. The sum of these $n-2$ eigenvalues is $-n+2k$. Then the sum of the two missing eigenvalues must be $2\left(1-2\frac{k}{n}\right)$. Let us denote them by ${\lambda }_k$ and ${\lambda }_k^{\ast }$. We must have $\mathit{Re}\left({\lambda }_k\right)=1-2\frac{k}{n}$. Then, since $\left|{\lambda }_k\right|=1$ we have

$${\lambda }_k=\left(1-2\frac{k}{n}\right)+i\frac{2}{n}\sqrt{k\left(n-k\right)}$$E56

Considering the eigenvalues of $-G_n$ and $I_n$ it is trivial to show that the dimensions of the eigenspaces of the oracle are:

$$dim{E}_{-}^O=M\text{and}dim{E}_{+}^O=N_e-M$$E57

5.3 Upper bound on the dimension of the complement

The eigenvalues of $U$ belong to $\left\{-1,+1,{\lambda }_k,{\lambda }_k^{\ast }\right\}$ where $k\in \left[1,n-1\right]$. Then, there are $2+2\left(n-1\right)=2n$ eigenspaces of $U$.

For $j\in \left[1,2n\right]$ let ${\alpha }_j$ be the dimensions of these eigenspaces and ${\beta }_j$ the dimensions of their intersections with $E_{+}^O$. An eigenvector of $U$ is in an intersection if and only if it is orthogonal to $E_{-}^O$. Then, because the dimension of $E_{-}^O$ is $M$, we have ${\beta }_j\ge {\alpha }_j-M$. Consequently

$$\sum _{j=1}^{2n}{\beta }_j\ge \sum _{j=1}^{2n}{\alpha }_j-2\mathit{nM}$$E58

Obviously, we have ${\sum }_{j=1}^{2n}{\alpha }_j=N_e$, so that

$$\sum _{j=1}^{2n}{\beta }_j\ge N_e-2\mathit{nM}$$E59

It follows that the dimension of the complement has an upper bound:

$$dim{E}_c\le 2\mathit{nM}$$E60

This is a remarkable result—despite the fact that the dimension of the Hilbert space grows exponentially ($N_e=n{2}^n$), the dimension of the complement grows only linearly with $n$.

5.4 Fast computation of the joint eigenspaces

5.4.1 Introduction

To check our theoretical upper bound, we propose an efficient algorithm for fast computation of the joint eigenspaces.

We have to compute orthonormal bases of joint eigenspaces of $U$ and $O$. The dimension of $E_{-}^O$ is small, hence, it makes sense to define it by an orthonormal basis generating the eigenspace. However, the dimension of $E_{+}^O$ is large (greater than $N_e/2$). Hence, it is computationally more efficient to define it by an orthonormal basis of its complement (which is $E_{-}^O$). Indeed $dim{E}_{-}^O\ll dim{E}_{+}^O$. We then have to design an algorithm adapted to each case.

5.4.2 Intersection of two eigenspaces defined by orthonormal bases

Let us consider a matrix $A$ whose columns are an orthonormal basis of an eigenspace of $U$, and a matrix $B$ whose columns are an orthonormal basis of $E_{-}^O$. Set $p$ and $q$ to be the number of columns of these matrices (their number of rows being $N_e$). We want to compute an $N_e\times r$ matrix $J$ whose columns are an orthonormal basis of the joint eigenspace (whose dimension we have set to be $r$). We propose the algorithm below, which is a straightforward adaptation of Theorem 1 in ref. [20].

First, we compute the $p\times q$ matrix $C$ below:

$$C=A^{\ast }B$$E61

Then, we compute the SVD of $C$:

$$C=U_c{S}_c{V}_c^{\ast }$$E62

Denote by $s_k$ the singular values (the diagonal elements of $S_c$) and determine $r$ such that $s_k\ge 1-\epsilon$ for $k=1,\dots ,r$, and $s_k<1-\epsilon$ for $k>r$. Here $\epsilon \ll 1$ is a very small positive value introduced to take into account the presence of small errors due to computer finite precision arithmetic. Finally:

$$J={\mathit{AU}}_c\left(:,1:r\right)$$E63

Or, equivalently, $J={\mathit{BV}}_c\left(:,1:r\right)$.

5.4.3 Intersection of two eigenspaces, one of them being defined by an orthonormal basis of its complement

Let us consider a matrix $A$ whose columns are an orthonormal basis of an eigenspace of $U$, and a matrix $B$ whose columns are an orthonormal basis of the complement of $E_{+}^O$ (that is $E_{-}^O$). First, we compute the $p\times q$ matrix $C$ (Eq. (61)). Then, we compute the $p\times r$ matrix ($r\le p$) $Z$ below:

$$Z=\mathit{\text{null}}\left(C^{\ast }\right)$$E64

and we obtain an $N_e\times r$ matrix $J$ whose columns are an orthonormal basis of $E_{\lambda ,+}^{U,O}$ from:

$$J=\mathit{AZ}$$E65

The justification of the algorithm is as follows. The $q$ columns of $C$ are a basis of the projection of $\mathit{Im}\left(B\right)$ into $\mathit{Im}\left(A\right)$, the components being expressed in the basis of $\mathit{Im}\left(A\right)$ The complement of $\mathit{Im}\left(C\right)$ in $\mathit{Im}\left(A\right)$ is the desired intersection (expressed in $\mathit{Im}\left(A\right)$). The columns of $Z$ are an orthonormal basis of this intersection. Finally, Eq. (65) restores the components in the original space.

5.5 Simulation results

Consider a hypercube of dimension $n=7$ with $M=3$ solutions located at nodes $\mathrm{2,8,9}$. The dimension of the state space is then $N_e=n{2}^n=896$. From the discussion above, we know that the dimension of the complement is upper bounded by $2\mathit{nM}=42$.

The algorithm gives us the dimensions of the joint eigenspaces of $U$ and $O$ (Table 3). The sum of the dimensions of the joint eigenspaces is then ${\sum }_{j=1}^{2n}{\beta }_j=858$, from which we obtain the dimension of the complement:

${\lambda }^O$	${\lambda }^U$	$dim{E}_{{\lambda }^O,{\lambda }^U}^{O,U}$
$-1$	any ${\lambda }_k$ or ${\lambda }_k^{\ast }$	0
$+1$	${\lambda }_0=+1$	321
$+1$	${\lambda }_1$ or ${\lambda }_1^{\ast }$	4
$+1$	${\lambda }_2$ or ${\lambda }_2^{\ast }$	18
$+1$	${\lambda }_3$ or ${\lambda }_3^{\ast }$	32
$+1$	${\lambda }_4$ or ${\lambda }_4^{\ast }$	32
$+1$	${\lambda }_5$ or ${\lambda }_5^{\ast }$	18
$+1$	${\lambda }_6$ or ${\lambda }_6^{\ast }$	4
$+1$	${\lambda }_7=-1$	321

Table 3.

Joint eigenspaces of $O$ and $U$ for $n=7$ and $M=3$ solutions located at nodes $\mathrm{2,8,9}$.

$$dim{E}_c=N_e-\sum _{j=1}^{2n}{\beta }_j=38$$E66

We can see that, as expected, this dimension ($dim{E}_c=38$) is much smaller than the dimension of the original state space ($N_e=896$). We can also check that it is less than the theoretical upper bound ($2\mathit{nM}=42$), as expected.

6. Conclusions

The recent growth of research on quantum communications and quantum information processing opens new challenges. In this chapter, we have shown that matrix theory concepts, such as JEVD, are powerful tools to propose new theoretical results as well as efficient simulation algorithms.

In the domain of quantum coding, we have shown how to determine the encoding matrix of a quantum code from a collection of Pauli errors. On a more speculative note to be part of future work concerning interception of quantum channels, it might also be useful to identify the quantum coder used by a noncooperative transmitter.

In the domain of quantum walk search, thanks to JEVD we have proved that there exists a small subspace of the whole Hilbert space which captures the essence of the search process, and we have given an algorithm that allows us to check this result by simulation.

Acknowledgments

The authors thank the IBNM (Institut Brestois du Numérique et des Mathématiques), CyberIoT Chair of Excellence, for its support.

Abbreviations

[qubit]	Quantum bit
[JEVD]	Joint EigenValue Decomposition
[SVD]	Singular Value Decomposition