Quantum Algorithm for DOA Estimation in Hybrid Massive MIMO

Throughout, we obey the following conventions. $\mathop{I}\nolimits_{N}$ refers to an $N\times N$ identity matrix. For an arbitrary matrix $A\in C^{M\times N}$, let $A=U\Lambda V^{H}$ be the singular value decomposition (SVD) of $A$. Here, $U=\left[{\mathop{u}\nolimits_{1},\mathop{u}\nolimits_{2}\ldots\mathop{u}\nolimits_{M}}\right]$ and $V=\left[{\mathop{v}\nolimits_{1},\mathop{v}\nolimits_{2}\ldots\mathop{v}\nolimits_{N}}\right]$ are unitary matrices with the left and right singular vectors of $X$ as columns, respectively. $\Lambda$ is a diagonal matrix with a singular value ${\mathop{\sigma}\nolimits_{i}}$ as a diagonal element. The economy-sized SVD can be $A=\sum\limits_{i=1}^{r}{\mathop{\sigma}\nolimits_{i}\mathop{u}\nolimits_{i}\mathop{v}\nolimits_{i}^{H}}$, where $r$ is the rank of $A$. The form of eigenvalue decomposition of $AA^{H}$ can be expressed as $AA^{H}{\rm{=}}\sum\limits_{i=1}^{r}{\mathop{\sigma}\nolimits_{i}^{2}\mathop{u}\nolimits_{i}\mathop{u}\nolimits_{i}^{H}}$, where ${\mathop{\sigma}\nolimits_{i}^{2}}$ is the eigenvalue of $AA^{H}$. $\left\|\cdot\right\|_{F}$ refers to Frobenius norm, and $\otimes$ denotes Kronecker product. The vectorization of a matrix $A$ can be denoted as ${\rm{vec}}\left(A\right)=\left(x_{11},x_{21}\ldots x_{M1},x_{12},x_{22}\ldots x_{M2},\ldots,x_{MN}\right)^{T}$ and ${\rm{vec}}\left(A\right)=\sum\limits_{i=1}^{r}{\mathop{\sigma}\nolimits_{i}\mathop{u}\nolimits_{i}\otimes\mathop{v}\nolimits_{i}}$.

Without loss of generality, we assume that there exists a set of predetermined DOA angles $\left\{\theta^{\left(i\right)}\right\}_{i=1}^{Q},{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}\theta^{\left(1\right)}<\theta^{\left(2\right)}<\ldots<\theta^{\left(Q\right)}$. The beam direction is switched by the analog beamformer to the predetermined DOA angle in turn. Then, the combination of the received signals in the receiver can be represented by $c_{q}\left(t\right)=a^{H}\left({\theta^{\left(q\right)}}\right)y\left(t\right)$, which is shown in Fig. 1.

Thus, a sample of the signal combination can be given by

$$c_{q}\left[n\right]=c_{q}\left({nT_{s}}\right)=a^{H}\left({\theta^{q}}\right)y\left[n\right]$$

(1)

where ${T_{s}}$ denotes the sample period and $y\left[n\right]=y\left({n\mathop{T}\nolimits_{s}}\right)$. The steering vector $a\left({\theta^{\left(q\right)}}\right)$ is defined as follows

$$a\left({{\theta^{\left(q\right)}}}\right)={\left[{1,{e^{-j\frac{{2\pi d}}{\lambda}\sin\left({{\theta^{\left(q\right)}}}\right)}},\ldots,{e^{-j\frac{{2\pi d\left({M-1}\right)}}{\lambda}\sin\left({{\theta^{\left(q\right)}}}\right)}}}\right]^{T}}$$

(2)

where the distance $d$ between consecutive elements is equal to half of the signal wavelength $\lambda$. Let $P_{q}$ be the average power of $c\left[n\right]$, then we can obtain

$$\begin{array}[]{l}P_{q}=\frac{1}{N}\sum_{n=0}^{N-1}{{\left|{c_{q}\left[n\right]}\right|}^{2}}\\ =a^{H}\left({\theta^{\left(q\right)}}\right)\frac{1}{N}\sum_{n=0}^{N-1}{y\left[n\right]y^{H}\left[n\right]}a\left({\theta^{\left(q\right)}}\right)\end{array}$$

(3)

When the number of samples is large enough, the sample average in Eq. (3) can be replaced by the statistical average, and we can obtain

$$\mathop{P}\nolimits_{q}=\mathop{a}\nolimits^{H}\left({\mathop{\theta}\nolimits^{\left(q\right)}}\right)Ra\left({\mathop{\theta}\nolimits^{\left(q\right)}}\right)$$

(4)

Using the vectorization definition to Eq. (4),

$$\begin{array}[]{l}vec\left\{{\mathop{a}\nolimits^{H}\left({\mathop{\theta}\nolimits^{\left(q\right)}}\right)Ra\left({\mathop{\theta}\nolimits^{\left(q\right)}}\right)}\right\}\\ =\left[{\mathop{a}\nolimits^{T}\left({\mathop{\theta}\nolimits^{\left(q\right)}}\right)\otimes\mathop{a}\nolimits^{H}\left({\mathop{\theta}\nolimits^{\left(q\right)}}\right)}\right]^{T}vec\left(R\right)\end{array}$$

(5)

Let $a_{q}=a\left({\theta^{\left(q\right)}}\right)\otimes a^{*}\left({\theta^{\left(q\right)}}\right)$ and $r=vec\left(R\right)$, then Eq. (4) can be rewritten as

$$\mathop{a}\nolimits_{q}^{T}r=\mathop{P}\nolimits_{q}$$

(6)

Given the group of predetermined DOA angles, we can obtain

$$Ar=P$$

(7)

where $A={\left({a_{1},a_{2},\ldots a_{Q}}\right)}^{T}$ is a $Q\times M^{2}$ matrix and $P={\left({P_{1},P_{2},\ldots P_{Q}}\right)}^{T}$ is a $Q\times 1$ vector. To avoid the ill-conditioned result, diagonal loading is adopted, and then we can obtain the vector form of the spatial covariance matrix as follows

$$\hat{r}={\left({A^{H}A+\sigma^{2}I_{M^{2}}}\right)}^{-1}A^{H}P$$

(8)

where $\hat{R}=unvec\left({\hat{r}}\right)$ is the desired spatial covariance matrix.

We assume that there are $L$ narrowband source signals incident upon an array of $M$ antenna elements. These $M$ elements are linearly spaced with equal distance between consecutive elements shown in Fig. 2.

Without loss of generality, we assume that the number of antenna elements is greater than the number of signals, i.e., $\left({M\gg L}\right)$. The distance between consecutive elements is equal to half of the signal wavelength $\lambda$; namely, $d_{1}=d_{2}=d_{3}=\cdots=d_{M}=\frac{\lambda}{2}$.

Let $y\left(k\right)={\left[{y_{1}\left(k\right),y_{2}\left(k\right),\ldots y_{M}\left(k\right)}\right]}^{T}$ be the $k$th observation data from $L$ source signals impinging on an array of $M$ elements; that is,

$$y\left(k\right)=\sum\limits_{i=0}^{L-1}{\mathop{s}\nolimits_{i}\left(k\right)a\left({\mathop{\theta}\nolimits_{i}}\right)}+n\left(k\right)=As\left(k\right)+n\left(k\right)$$

(9)

where $A=\left[{a\left({{\theta_{0}}}\right),a\left({{\theta_{1}}}\right),\ldots a\left({{\theta_{L-1}}}\right)}\right]$ is an $M\times L$ array matrix, signal vector and noise vector are independent. $s\left(k\right)=\mathop{\left[{\mathop{s}\nolimits_{0}\left(k\right),\mathop{s}\nolimits_{1}\left(k\right),\ldots,\mathop{s}\nolimits_{L-1}\left(k\right)}\right]}\nolimits^{T}$ is an incident signal vector with zero mean value, $n\left(k\right)$ is the Gaussian noise vector with zero mean value and $\sigma^{2}I_{N}$ covariance matrix. Then, the $M\times M$ spatial covariance matrix of the received signal vector in Eq. (9) can be obtained as

$$R=E\left\{{y\left(k\right)y^{H}\left(k\right)}\right\}$$

(10)

As the spatial covariance matrix is Hermitian, the eigen-decomposition form can be represented as

$$R{\rm{=}}{U_{s}}{\Lambda_{s}}{\rm{}}U_{s}^{H}+{\rm{}}{U_{n}}{\Lambda_{n}}{\rm{}}U_{n}^{H}$$

(11)

where ${U_{s}}=\left[u_{1},u_{2},\dots,u_{L}\right]$ is an $M\times L$ matrix representing signal subspace, ${U_{n}}=\left[u_{L+1},u_{L+2},\dots,u_{M}\right]$ is a group of the orthogonal base vectors of noise subspace.

The object of the direction angle estimation is

$$\begin{array}[]{c}{\theta_{MUSIC}}=\mathop{\arg\min}\limits_{\theta}{\rm{}}{a^{H}}\left(\theta\right){\rm{}}{U_{n}}{\rm{}}U_{n}^{H}a\left(\theta\right)\end{array}$$

(12)

Analogously, direction angle estimation can also be represented in terms of its reciprocal to obtain peaks; that is,

$$\begin{array}[]{c}\mathop{P}\nolimits_{MUSIC}=\frac{1}{{\mathop{a}\nolimits^{H}\left(\theta\right)\mathop{U}\nolimits_{n}\mathop{U}\nolimits_{n}^{H}a\left(\theta\right)}}\end{array}$$

(13)

In conventional MUSIC algorithm, the distribution of received signal vector ${y\left(k\right)}$ is unknown. Therefore, the spatial covariance matrix in Eq. (10) can be estimated by the sample average, that is $R\approx\frac{1}{N}\sum_{k=1}^{N}{y\left(k\right)y^{H}\left(k\right)}$ where $N$ indicates the number of samples.

In this subroutine, to obtain the vector form of the spatial covariance matrix in Eq. (8), we assume that $P$ is stored in quantum random access memory Kerenidis and Prakash (2017) with a suitable data structure of binary trees. Then, the quantum state $\left|P\right\rangle$ can be prepared as

$$U_{P}:\left|0\right\rangle\to\frac{1}{{{\left\|P\right\|}_{2}}}\sum\nolimits_{i=0}^{Q-1}{P_{i}}\left|i\right\rangle$$

(14)

Next, we design an efficient quantum circuit to prepare any row of $A$. As every row of $A$ can be presented as Kronecker product of the steering vector and its conjugate $a_{q}=a\left({\theta^{\left(q\right)}}\right)\otimes a^{*}\left({\theta^{\left(q\right)}}\right),q=1,2,\dots,Q$. Therefore, we first design a quantum circuit to prepare $\left|a\left(\theta^{q}\right)\right\rangle$, the detailed circuit can be shown in Fig. 3.

Here, we shown an example of preparing a $4\times 1$ steering vector to demonstrate the quantum circuit implementation.

$\displaystyle\begin{array}[]{l}\frac{1}{2}\left|{\frac{{2\pi d\sin\theta}}{\lambda}}\right\rangle\left({\left|0\right\rangle+\left|1\right\rangle}\right)\left({\left|0\right\rangle+\left|1\right\rangle}\right)\left|0\right\rangle=\frac{1}{2}\left|{\frac{{2\pi d\sin\theta}}{\lambda}}\right\rangle\left({\left|0\right\rangle\left|0\right\rangle+\left|0\right\rangle\left|1\right\rangle+\left|1\right\rangle\left|0\right\rangle+\left|1\right\rangle\left|1\right\rangle}\right)\left|0\right\rangle\\ =\frac{1}{2}\left|{\frac{{2\pi d\sin\theta}}{\lambda}}\right\rangle\left({\left|0\right\rangle\left|0\right\rangle\left|0\right\rangle+\left|0\right\rangle\left|1\right\rangle\mathop{e}\nolimits^{-j\frac{{2\pi d\sin\theta}}{\lambda}}\left|0\right\rangle+\left|1\right\rangle\left|0\right\rangle\mathop{e}\nolimits^{-j\frac{{2\pi 2d\sin\theta}}{\lambda}}\left|0\right\rangle+\left|1\right\rangle\left|1\right\rangle\mathop{e}\nolimits^{-j\frac{{2\pi 3d\sin\theta}}{\lambda}}\left|0\right\rangle}\right)\\ =\frac{1}{2}\left|{\frac{{2\pi d\sin\theta}}{\lambda}}\right\rangle\left({\left|0\right\rangle\left|0\right\rangle+\mathop{e}\nolimits^{-j\frac{{2\pi d\sin\theta}}{\lambda}}\left|0\right\rangle\left|1\right\rangle+\mathop{e}\nolimits^{-j\frac{{2\pi 2d\sin\theta}}{\lambda}}\left|1\right\rangle\left|0\right\rangle+\mathop{e}\nolimits^{-j\frac{{2\pi 3d\sin\theta}}{\lambda}}\left|1\right\rangle\left|1\right\rangle}\right)\left|0\right\rangle\\ =\left|{\frac{{2\pi d\sin\theta}}{\lambda}}\right\rangle\left|{a\left(\theta\right)}\right\rangle\left|0\right\rangle\end{array}$

(19)

Thus, any row of $A$ can be prepared in following quantum circuit Fig. 4, and then we define two unitary mappings $\mathop{U}\nolimits_{\rm M}$ and $\mathop{U}\nolimits_{\rm N}$ as follows

$$\begin{array}[]{l}U_{\rm M}:\left|i\right\rangle\left|0\right\rangle\to\frac{1}{{{\left\|{A_{i}}\right\|}_{2}}}\sum\nolimits_{j=0}^{M^{2}-1}{A_{ij}\left|i\right\rangle\left|j\right\rangle}=\sum\nolimits_{j=0}^{M^{2}-1}{\frac{A_{ij}}{M}\left|i\right\rangle\left|j\right\rangle}\end{array}$$

(20)

$$\begin{array}[]{l}\mathop{U}\nolimits_{\rm N}:\left|0\right\rangle\left|j\right\rangle\to\frac{1}{{\mathop{\left\|A\right\|}\nolimits_{F}}}\sum\nolimits_{i=0}^{Q-1}{\mathop{\left\|{\mathop{A}\nolimits_{i}}\right\|}\nolimits_{2}\left|i\right\rangle\left|j\right\rangle}\\ =\frac{1}{{M\sqrt{Q}}}\sum\nolimits_{i=0}^{Q-1}{M\left|i\right\rangle\left|j\right\rangle}\end{array}$$

(21)

where ${A_{i}}$ is the $i$th row of $A$. Obviously, $U_{\rm N}$ can be implemented with the complexity ${\rm O}\left(\log Q\right)$ , and based on circuit Fig. 3 and 4, the implementation of $U_{\rm M}$ can be performed in the complexity ${\rm O}\left(\mathrm{poly}\log MQ\right)$. Moreover, we can obtain the following normalization form of $A$,

$$\left\langle i\right|\left\langle 0\right|\mathop{U}\nolimits_{\rm M}^{H}\mathop{U}\nolimits_{\rm N}\left|0\right\rangle\left|j\right\rangle=\frac{{\mathop{A}\nolimits_{ij}}}{{\mathop{\left\|A\right\|}\nolimits_{F}}}=\frac{{\mathop{A}\nolimits_{ij}}}{{M\sqrt{Q}}}$$

(22)

Then, the quantum state form of ${\hat{r}}$ can be represented

$\displaystyle\left|{\hat{r}}\right\rangle\propto{\left({A^{H}A+\sigma^{2}I_{M^{2}}}\right)}^{-1}A^{H}\left|P\right\rangle=\sum\nolimits_{i=1}^{\gamma}{\frac{{{\alpha_{i}\sigma}_{i}}}{{\sigma_{i}^{2}+\sigma^{2}}}}\left|{v_{i}}\right\rangle$

(23)

where $\alpha_{i}=\left\langle{u_{i}}\right|\left.P\right\rangle$ and $\gamma$ is the rank of $A$, and $\left|{v_{i}}\right\rangle$, $\left|{u_{i}}\right\rangle$ and $\sigma_{i}$ are the right, left singular vectors and corresponding singular value, respectively. By considering $U_{\rm M}$ and $U_{\rm N}$, we prove the following singular vector transition theorem,

See Appendix A for details of the singular vector transition theorem. Therefore, based on the definition of $U_{\rm M}$, $U_{\rm N}$ and the above theorem, the detailed quantum algorithm is as follows:
(1) Prepare the state as follows

$$\left|{\hat{\varphi}}\right\rangle{\rm{=}}\left|0\right\rangle\left|P\right\rangle{\rm{=}}\sum\nolimits_{i=1}^{Q}{\mathop{\alpha}\nolimits_{i}\left|0\right\rangle\left|{\mathop{u}\nolimits_{i}}\right\rangle}$$

(25)

(2) Apply the quantum singular value estimation (QSVE) technique Kerenidis and Prakash (2017), we have the state

$$\left|{\hat{\varphi}}\right\rangle{\rm{=}}\sum\nolimits_{i=1}^{Q}{\mathop{\alpha}\nolimits_{i}\left|{\mathop{\sigma}\nolimits_{i}}\right\rangle\left|{\mathop{u}\nolimits_{i}}\right\rangle}$$

(26)

(3) Add a register with ${\left|0\right\rangle}$ and perform the controlled rotation operation, controlled by the register storing the singular value, we have the following state

$$\begin{array}[]{l}\left|{\hat{\varphi}}\right\rangle=\\ \sum\nolimits_{i=1}^{Q}{\mathop{\alpha}\nolimits_{i}\left({\frac{{C\mathop{\sigma}\nolimits_{i}}}{{\mathop{\sigma}\nolimits_{i}^{2}+\mathop{\sigma}\nolimits^{2}}}\left|0\right\rangle+\sqrt{1-\mathop{\left({\frac{{C\mathop{\sigma}\nolimits_{i}}}{{\mathop{\sigma}\nolimits_{i}^{2}+\mathop{\sigma}\nolimits^{2}}}}\right)}\nolimits^{2}}\left|1\right\rangle}\right)\left|{\mathop{\sigma}\nolimits_{i}}\right\rangle\left|{\mathop{u}\nolimits_{i}}\right\rangle}\end{array}$$

(27)

where $C$ is a constant, that is $C={\raise 3.01385pt\hbox{$1$}\!\mathord{\left/{\vphantom{1{\mathop{\max}\limits_{i}\frac{{\mathop{\sigma}\nolimits_{i}}}{{\mathop{\sigma}\nolimits_{i}^{2}+\mathop{\sigma}\nolimits^{2}}}}}}\right.\kern-1.2pt}\!\lower 3.01385pt\hbox{${\mathop{\max}\limits_{i}\frac{{\mathop{\sigma}\nolimits_{i}}}{{\mathop{\sigma}\nolimits_{i}^{2}+\mathop{\sigma}\nolimits^{2}}}}$}}$.
(4) Uncomputing the second register and measure the first register $\left|0\right\rangle$ with the probability $P\left({\left|0\right\rangle}\right)$, we have the state

$$\left|{\hat{\varphi}}\right\rangle{\rm{=}}\frac{1}{{\sqrt{\sum\nolimits_{i=1}^{\gamma}{\mathop{\left|{\mathop{\alpha}\nolimits_{i}}\right|}\nolimits^{2}\mathop{\left({\frac{{\mathop{\sigma}\nolimits_{i}}}{{\mathop{\sigma}\nolimits_{i}^{2}+\mathop{\sigma}\nolimits^{2}}}}\right)}\nolimits^{2}}}}}\sum\nolimits_{i=1}^{\gamma}{\frac{{\mathop{\alpha}\nolimits_{i}\mathop{\sigma}\nolimits_{i}}}{{\mathop{\sigma}\nolimits_{i}^{2}+\mathop{\sigma}\nolimits^{2}}}}\left|{\mathop{u}\nolimits_{i}}\right\rangle$$

(28)

where $P\left({\left|0\right\rangle}\right)=\sum\nolimits_{i=1}^{\gamma}{\mathop{\left|{\mathop{\alpha}\nolimits_{i}}\right|}\nolimits^{2}\mathop{\left({\frac{{\mathop{C\sigma}\nolimits_{i}}}{{\mathop{\sigma}\nolimits_{i}^{2}+\mathop{\sigma}\nolimits^{2}}}}\right)}\nolimits^{2}}$.
(5) Perform the unitary transformation ${\rm T}:\left|{\mathop{u}\nolimits_{i}}\right\rangle\to\left|{\mathop{v}\nolimits_{i}}\right\rangle$, we can obtain

$$\left|{\hat{r}}\right\rangle=\left|{\hat{\varphi}}\right\rangle{\rm{=}}\frac{1}{{\sqrt{\sum\nolimits_{i=1}^{\gamma}{\mathop{\left|{\mathop{\alpha}\nolimits_{i}}\right|}\nolimits^{2}\mathop{\left({\frac{{\mathop{\sigma}\nolimits_{i}}}{{\mathop{\sigma}\nolimits_{i}^{2}+\mathop{\sigma}\nolimits^{2}}}}\right)}\nolimits^{2}}}}}\sum\nolimits_{i=1}^{\gamma}{\frac{{\mathop{\alpha}\nolimits_{i}\mathop{\sigma}\nolimits_{i}}}{{\mathop{\sigma}\nolimits_{i}^{2}+\mathop{\sigma}\nolimits^{2}}}}\left|{\mathop{v}\nolimits_{i}}\right\rangle$$

(29)

where more details of ${\rm T}$ is deferred to Appendix A. Thus, we can obtain the quantum state $\left|{\hat{\varphi}}\right\rangle$, which is proportional to the vector form of the spatial covariance matrix. Subsequently, for the goal of eigen-decomposition in MUSIC, we transform $\left|{\hat{\varphi}}\right\rangle$ into a density matrix based on Schmidt decomposition theorem,

By Schmidt decomposition theorem, $\left|{\hat{\varphi}}\right\rangle$ can be rewritten as follows

$$\begin{array}[]{c}\left|{\hat{\varphi}}\right\rangle{\rm{=}}\frac{1}{{\sqrt{\sum\nolimits_{i=1}^{\gamma}{{{\left|{{\alpha_{i}}}\right|}^{2}}{{\left({\frac{{{\sigma_{i}}}}{{\sigma_{i}^{2}+{\sigma^{2}}}}}\right)}^{2}}}}}}\sum\nolimits_{i=1}^{\gamma}{\frac{{{\alpha_{i}}{\sigma_{i}}}}{{\sigma_{i}^{2}+{\sigma^{2}}}}}\left|{{\rm{}}{v_{i}}}\right\rangle\\ {\rm{=}}\sum\nolimits_{j=1}^{\gamma}{\sigma_{j}^{\hat{r}}}{{\left|{{\rm{}}u_{j}^{\hat{r}}}\right\rangle}_{A}}{{\left|{{\rm{}}v_{j}^{\hat{r}}}\right\rangle}_{B}}\end{array}$$

(31)

where $\left|u_{j}^{{\hat{r}}}\right\rangle$ and $\left|v_{j}^{\hat{r}}\right\rangle$ are left and right singular vectors of $\left|{\hat{r}}\right\rangle$, $\sigma_{j}^{\hat{r}}$ is corresponding singular value, respectively. As the spatial covariance matrix is positive definite, its singular vectors are the same as eigenvectors, namely, $\left|u_{j}^{\hat{r}}\right\rangle_{A}=\left|v_{j}^{\hat{r}}\right\rangle_{B}$. Therefore, we take the partial trace on the $B$ register of $\left|{\hat{\varphi}}\right\rangle\left\langle{\hat{\varphi}}\right|$ and the density matrix representation $\mathop{\rho}\nolimits_{\hat{R}}$ can be denoted as

$$\rho_{\hat{R}}={{\rm tr}}_{2}\left({\left|{\hat{\varphi}}\right\rangle\left\langle{\hat{\varphi}}\right|}\right)=\sum\nolimits_{i=1}^{\gamma}{{\left({\sigma_{j}^{\hat{r}}}\right)}^{2}}{\left|{u_{j}^{\hat{r}}}\right\rangle}_{A}\left\langle{u_{j}^{\hat{r}}}\right|$$

(32)

where $\rho_{\hat{R}}$ have the same eigenvectors as the spatial covariance matrix.

Here, we propose the variational quantum eigensolver for the density matrix, which variationally learns the largest nozero eigenvalues of the density matrix as well as a gate sequence $V\left(\theta\right)$ that prepares the corresponding eigenvectors. Before elaborately introuducing our variational quantum eigensolver, we first briefly review following von Neumann theorem,

Based on von Neumann theorem and Subspace-search variational quantum eigensolver Nakanishi et al. (2019), our variational quantum eigensolver can be depicted in Algorithm 1

In step 1, given a series of parameter vectors $\overrightarrow{\theta}=\left(\theta_{1},\theta_{2},\dots\theta_{m}\right)$, the quantum circuit $V$ is defined as

$$V\left(\overrightarrow{\theta}\right)=V_{m}\left(\theta_{m}\right)V_{m-1}\left(\theta_{m-1}\right)\dots V_{1}\left(\theta_{1}\right)$$

(34)

where the ansatz circuit can be shown in Fig. 5. Note that the number of parameters is logarithmically proportional to the dimension of the density matrix. These parameterized quantum circuits have shown significant potential power in quantum neural network and quantum circuit Born machines Marcello et al. (2019); Killoran et al. (2019); Perdomo et al. (2019).

In step 3, $C\left(\overrightarrow{\theta}\right)$ can be denoted as

$$C\left(\overrightarrow{\theta}\right)=\sum\nolimits_{i=1}^{L}\frac{{q_{i}}}{\left({\textstyle\sum_{i=1}^{L}}q_{i}\right)}\left\langle{\varphi_{i}}\right|V^{H}\left(\overrightarrow{\theta}\right)\rho_{\hat{R}}V\left(\overrightarrow{\theta}\right)\left|{\varphi_{i}}\right\rangle$$

(35)

Further, we can transform $C\left(\overrightarrow{\theta}\right)$ into a trace form as follows

$$C\left(\overrightarrow{\theta}\right)=\mathrm{tr}\left(V^{H}\left(\overrightarrow{\theta}\right)\rho_{\hat{R}}V\left(\overrightarrow{\theta}\right)\rho_{f}\right)$$

(36)

where $\rho_{f}={\textstyle\sum_{i=1}^{L}}\frac{q_{i}}{{\textstyle\sum_{i=1}^{L}}q_{i}}\left|\varphi_{i}\right\rangle\left\langle\varphi_{i}\right|$. Unlike existing variational quantum eigensolvers, the density matrix $\rho_{\hat{R}}$ can not be represented as a linear combination of unitary matrices. This is because the density matrix is the output of the quantum algorithm in the previous quantum subroutine, we do not have any prior information on $\rho_{\hat{R}}$ without quantum tomography. However, we transform the expection form in the objective function into the trace of two density matrices. Thus, in our quantum algorithm, a linear combination of unitary matrices is not required and the objective function $C\left(\overrightarrow{\theta}\right)$ can effectively be implemented in Destructive Swap Test, which is shown in Fig. 6. Moreover, the detailed quantum implementation scheme is depicted in Fig. 7, and for the implementation, we consider the following a $4\times 4$ density matrix (using two qubit).

Example: An random density matrix $\rho_{in}$ with four fixed eigenvalues $\lambda_{1}=0.4,\lambda_{2}=0.3,\lambda_{3}=0.2,\lambda_{4}=0.1$, a group of positive real weights are $q_{1}=4,q_{2}=3,q_{3}=2,q_{4}=1$ and $\left|\varphi_{1}\right\rangle=\left|1\right\rangle,\left|\varphi_{2}\right\rangle=\left|2\right\rangle,\left|\varphi_{3}\right\rangle=\left|3\right\rangle,\left|\varphi_{4}\right\rangle=\left|4\right\rangle$. The experiment’s results of the VQDME implementation are shown in Fig. 8. As shown in the above figures, we can obtain four eigenvalues and eigenvectors within 200 quantum-classical hybrid iterations, which are very close to the exact ones.

To implement the direction angles estimation, we can rewrite the object (12) as

$$\begin{array}[]{l}\mathop{\theta}\nolimits_{MUSIC}=\mathop{\min}\limits_{\theta}\mathop{a}\nolimits^{H}\left(\theta\right)\mathop{U}\nolimits_{n}\mathop{U}\nolimits_{n}^{H}a\left(\theta\right)=\\ \mathop{\min}\limits_{\theta}\mathop{a}\nolimits^{H}\left(\theta\right)\left({I-\mathop{U}\nolimits_{s}\mathop{U}\nolimits_{s}^{H}}\right)a\left(\theta\right)=\\ \mathop{\min}\limits_{\theta}\left({M-\mathop{a}\nolimits^{H}\left(\theta\right)\mathop{U}\nolimits_{s}\mathop{U}\nolimits_{s}^{H}a\left(\theta\right)}\right)=\\ \mathop{\max}\limits_{\theta}\mathop{a}\nolimits^{H}\left(\theta\right)\mathop{U}\nolimits_{s}\mathop{U}\nolimits_{s}^{H}a\left(\theta\right)=\mathop{\max}\limits_{\theta}{\textstyle\sum_{i=1}^{L}}\left(\alpha_{i}^{\theta}\right)^{2}\end{array}$$

(37)

where $\alpha_{i}^{\theta}=a^{H}\left(\theta\right)u_{i},i=1,2,\dots,L$ are the projections on the signal subspace. Although we have implemented the optimal quantum network $V\left(\overrightarrow{\theta}^{\ast}\right)$, it is hard to construct the controlled $V\left(\overrightarrow{\theta}^{\ast}\right)$. Namely, we can not prepare a coherent state of $\left\{\left|u_{i}\right\rangle\right\}_{i=1}^{L}$ and the projections on the signal subspace can also be calculated by the quantum parallelism. Fortunately, the projections on entire space can be obtained by the quantum computing, and a quantum labeling algorithm is required to extract the projections on the signal subspace from ones on entire space. Here, we design a quantum labeling operation based on the set $\left\{\left|\varphi_{i}\right\rangle\right\}_{i=1}^{L}$ and the detailed implementation is as follows:
(1) Prepare the coherent quantum state $\left|{\mathop{\phi}\nolimits_{S}}\right\rangle$ for the searching space vector set $\left|{a\left({\mathop{\theta}\nolimits_{n}}\right)}\right\rangle$, $n=1,2,\ldots,K$, as

$$\begin{array}[]{l}\left|{\mathop{\phi}\nolimits_{S}}\right\rangle=\frac{1}{{\sqrt{K}}}\sum\nolimits_{n=1}^{K}{\left|n\right\rangle\left|{a\left({\mathop{\theta}\nolimits_{n}}\right)}\right\rangle}=\\ \frac{1}{{\sqrt{K}}}\sum\nolimits_{n=1}^{K}{\left|n\right\rangle\sum\nolimits_{i=1}^{M}{\mathop{\alpha}\nolimits_{i}^{n}\left|{\mathop{u}\nolimits_{i}^{\hat{r}}}\right\rangle}}\end{array}$$

(38)

where ${\left|{\mathop{u}\nolimits_{i}^{\hat{r}}}\right\rangle}$ is the eigenvector of $\mathop{\rho}\nolimits_{\hat{R}}$ and $\mathop{\alpha}\nolimits_{i}^{n}=\left\langle{\mathop{u}\nolimits_{i}^{\hat{r}}}\right.\left|{a\left({\mathop{\theta}\nolimits_{n}}\right)}\right\rangle,i=1,2,\dots,M$ is the projections on the entire space of searching space vector.

(2) Perform the inverse optimal quantum network $\mathop{V}\nolimits^{H}\left({\mathop{\theta}\nolimits^{*}}\right)$ on the second register, we have the following state

$$\left|{\mathop{\phi}\nolimits_{S}}\right\rangle=\frac{1}{{\sqrt{K}}}\sum\nolimits_{n=1}^{K}{\left|n\right\rangle\sum\nolimits_{i=1}^{M}{\mathop{\alpha}\nolimits_{i}^{n}\left|{\mathop{\varphi}\nolimits_{i}}\right\rangle}}$$

(39)

(3) Append a single qubit label register with the initial $\left|0\right\rangle$. As we select some certain orthogonal basis vectors as ${\left|{\mathop{\varphi}\nolimits_{i}}\right\rangle}$, without loss of generality let $\left\{\left|{\mathop{\varphi}\nolimits_{i}}\right\rangle\right\}_{i=1}^{L}$ be the computational basis set $\left\{\left|i\right\rangle\right\}_{i=1}^{L}$, then we can construct the labeling mapping as

$${\textstyle\sum_{k=1}^{L}}\left|k\right\rangle\left\langle k\right|\otimes X+\left(I_{M}-{\textstyle\sum_{k=1}^{L}}\left|k\right\rangle\left\langle k\right|\right)\otimes I_{2}$$

(40)

and the operation on $\left\{\left|{\mathop{\varphi}\nolimits_{i}}\right\rangle\right\}_{i=1}^{L}$ can be denoted as

$$CU:\left|{\mathop{\varphi}\nolimits_{i}}\right\rangle\mathop{\left|0\right\rangle}\nolimits^{{\mathop{\rm label}\nolimits}}\to\left|{\mathop{\varphi}\nolimits_{i}}\right\rangle\mathop{\left|1\right\rangle}\nolimits^{{\mathop{\rm label}\nolimits}},{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}i=1,2,\ldots,L$$

(41)

$$CU:\left|{\mathop{\varphi}\nolimits_{i}}\right\rangle\mathop{\left|0\right\rangle}\nolimits^{{\mathop{\rm label}\nolimits}}\to\left|{\mathop{\varphi}\nolimits_{i}}\right\rangle\mathop{\left|0\right\rangle}\nolimits^{{\mathop{\rm label}\nolimits}},{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}i=L+1,\ldots,M$$

(42)

Thus, we have the state as follows

$$\begin{array}[]{l}CU:\left|{\mathop{\phi}\nolimits_{S}}\right\rangle\to\\ \sum\nolimits_{n=1}^{K}{\frac{{\left|n\right\rangle}}{{\sqrt{K}}}\left({\sum\nolimits_{i=1}^{L}{\mathop{\alpha}\nolimits_{i}^{n}\left|{\mathop{\varphi}\nolimits_{i}}\right\rangle\mathop{\left|1\right\rangle}\nolimits^{{\mathop{\rm label}\nolimits}}+\sum\nolimits_{i=L+1}^{M}{\mathop{\alpha}\nolimits_{i}^{n}\left|{\mathop{\varphi}\nolimits_{i}}\right\rangle\mathop{\left|0\right\rangle}\nolimits^{{\mathop{\rm label}\nolimits}}}}}\right)}\end{array}$$

(43)

Therefore, the projections of $\left|{a\left({\mathop{\theta}\nolimits_{n}}\right)}\right\rangle$ on the signal subspace labeled ${\mathop{\left|1\right\rangle}\nolimits^{{\mathop{\rm label}\nolimits}}}$ are separated from projections on the noise subspace with label ${\mathop{\left|0\right\rangle}\nolimits^{{\mathop{\rm label}\nolimits}}}$.
(4) Measure the label register in ${\mathop{\left|1\right\rangle}\nolimits^{{\mathop{\rm label}\nolimits}}}$ with the probability $\mathop{P}\nolimits_{S}\left({\left|1\right\rangle}\right)$, then we have the state

$$\left|{\mathop{\phi}\nolimits_{S}}\right\rangle=\frac{1}{{\sqrt{\sum\nolimits_{n=1}^{K}{\sum\nolimits_{i=1}^{L}{\mathop{\left({\mathop{\alpha}\nolimits_{i}^{n}}\right)}\nolimits^{2}}}}}}\sum\nolimits_{n=1}^{K}{\left|n\right\rangle\sum\nolimits_{i=1}^{L}{\mathop{\alpha}\nolimits_{i}^{n}\left|{\mathop{\varphi}\nolimits_{i}}\right\rangle}}$$

(44)

(5) Perform a small amount of sampling on the first register of $\left|{\mathop{\phi}\nolimits_{S}}\right\rangle$ with the probability $P\left({\left|n\right\rangle}\right)=\frac{{\sum\nolimits_{i=1}^{L}{\mathop{\left({\mathop{\alpha}\nolimits_{i}^{n}}\right)}\nolimits^{2}}}}{{\sqrt{\sum\nolimits_{n=1}^{K}{\sum\nolimits_{i=1}^{L}{\mathop{\left({\mathop{\alpha}\nolimits_{i}^{n}}\right)}\nolimits^{2}}}}}}$, and satisfy

$$\frac{{\sum\nolimits_{i=1}^{L}{\mathop{\left({\mathop{\alpha}\nolimits_{i}^{n}}\right)}\nolimits^{2}}}}{{\sqrt{\sum\nolimits_{n=1}^{K}{\sum\nolimits_{i=1}^{L}{\mathop{\left({\mathop{\alpha}\nolimits_{i}^{n}}\right)}\nolimits^{2}}}}}}\propto\mathop{a}\nolimits^{H}\left(\theta\right)\mathop{U}\nolimits_{s}\mathop{U}\nolimits_{s}^{H}a\left(\theta\right)$$

(45)

As sample can be more easily selected with higher probability, the most frequently selected samples can represent the direction estimation values that satisfy the object (33).

Time complexity for $\left|{\hat{\varphi}}\right\rangle$: The state $\left|P\right\rangle$ can be prepared with the complexity ${\rm O}\left({{\mathop{\rm poly}\nolimits}\left({\log Q}\right)}\right)$, $\mathop{U}\nolimits_{\rm N}$ and $\mathop{U}\nolimits_{\rm M}$ can be implemented in the runtime ${\rm O}\left({\log Q}\right)$ and ${\rm O}\left({{\mathop{\rm poly}\nolimits}\left({\log QM}\right)}\right)$, respectively. Let the error in SVE be $\mathop{\varepsilon}\nolimits_{p}\mathop{\left\|A\right\|}\nolimits_{F}$, where $\mathop{\varepsilon}\nolimits_{p}$ is the error in the phase estimation. We define the function $h\left(\sigma\right)$ as

$$h\left(\lambda\right)=\frac{\lambda}{{\mathop{\lambda}\nolimits^{2}+\mathop{\sigma}\nolimits^{2}}},{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}\lambda\in\left[{\frac{1}{{\mathop{\kappa}\nolimits_{A}}},1}\right]$$

(46)

where ${\mathop{\kappa}\nolimits_{A}}$ is the condition number of $A$. Then, we take the derivative of $h\left(\sigma\right)$ as

$$h^{\prime}\left(\lambda\right)=\left\{{\begin{array}[]{*{20}{c}}{\frac{{\mathop{\sigma}\nolimits^{2}-\mathop{\lambda}\nolimits^{2}}}{{\left({\mathop{\lambda}\nolimits^{2}+\mathop{\sigma}\nolimits^{2}}\right)}}<0,\lambda>\sigma{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}}\\ {\frac{{\mathop{\sigma}\nolimits^{2}-\mathop{\lambda}\nolimits^{2}}}{{\left({\mathop{\lambda}\nolimits^{2}+\mathop{\sigma}\nolimits^{2}}\right)}}\geq 0,0\leq\lambda\leq\sigma{\kern 1.0pt}{\kern 1.0pt}}\end{array}}\right.$$

(47)

As the parameter ${\mathop{\sigma}\nolimits^{2}}$ is required to be a number closed to 0 in the reconstruction of the spatial covariance matrix, we can deduce $\sigma\leq\frac{1}{{\mathop{\kappa}\nolimits_{A}}}$. Therefore, $h\left(\sigma\right)$ is demonstrated to be monotonic decreasing with $\left[{\frac{1}{{\mathop{\kappa}\nolimits_{A}}},1}\right]$. Thus, we can estimate

$$\frac{1}{{1+\mathop{\sigma}\nolimits^{2}}}\leq\sqrt{\sum\nolimits_{i=1}^{Q}{\mathop{\left|{\mathop{\alpha}\nolimits_{i}}\right|}\nolimits^{2}\mathop{\left({\frac{{\mathop{\sigma}\nolimits_{i}}}{{\mathop{\sigma}\nolimits_{i}^{2}+\mathop{\sigma}\nolimits^{2}}}}\right)}\nolimits^{2}}}<\mathop{\kappa}\nolimits_{A}$$

(48)

, $C=\frac{1}{{\mathop{\kappa}\nolimits_{A}}}+\mathop{\kappa}\nolimits_{A}\mathop{\sigma}\nolimits^{2}$ and

$\displaystyle P\left({\left|0\right\rangle}\right)=\sum\nolimits_{i=1}^{Q}{\mathop{\left|{\mathop{\alpha}\nolimits_{i}}\right|}\nolimits^{2}\mathop{\left({\frac{{\mathop{C\sigma}\nolimits_{i}}}{{\mathop{\sigma}\nolimits_{i}^{2}+\mathop{\sigma}\nolimits^{2}}}}\right)}\nolimits^{2}}\geq\mathop{C}\nolimits^{2}\mathop{\left({\frac{1}{{1+\mathop{\sigma}\nolimits^{2}}}}\right)}\nolimits^{2}=\mathop{\left({\frac{1}{{\mathop{\kappa}\nolimits_{A}}}+\mathop{\kappa}\nolimits_{A}\mathop{\sigma}\nolimits^{2}}\right)}\nolimits^{2}\mathop{\left({\frac{1}{{1+\mathop{\sigma}\nolimits^{2}}}}\right)}\nolimits^{2}>\frac{1}{4}\mathop{\left({\frac{1}{{\mathop{\kappa}\nolimits_{A}}}+\mathop{\kappa}\nolimits_{A}\mathop{\sigma}\nolimits^{2}}\right)}\nolimits^{2}$

(49)

Moreover, we can estimate $\left|{h\left(\lambda\right)-h\left({\tilde{\lambda}}\right)}\right|$ as

$$\begin{array}[]{l}\left|{h\left(\lambda\right)-h\left({\tilde{\lambda}}\right)}\right|\approx\frac{{\left|{h\left(\lambda\right)-h\left({\tilde{\lambda}}\right)}\right|}}{{\left|{\lambda-\tilde{\lambda}}\right|}}\left|{\lambda-\tilde{\lambda}}\right|\approx\left|{h^{\prime}\left(\lambda\right)}\right|\left|{\lambda-\tilde{\lambda}}\right|<\\ \frac{{\mathop{\kappa}\nolimits_{A}\mathop{\varepsilon}\nolimits_{p}\mathop{\left\|A\right\|}\nolimits_{F}}}{{1+\mathop{\mathop{\sigma}\nolimits^{2}\kappa}\nolimits_{A}^{2}}}\end{array}$$

(50)

We define the ideal state $\left|{\mathop{\varphi}\nolimits^{*}}\right\rangle$, then we have

$$\begin{array}[]{l}\mathop{\left\|{\left|{\hat{\varphi}}\right\rangle-\left|{\mathop{\varphi}\nolimits^{*}}\right\rangle}\right\|}\nolimits^{2}\leq\mathop{\left\|{\sum\nolimits_{i=1}^{Q}{\frac{{\mathop{\sigma}\nolimits_{i}}}{{\mathop{\sigma}\nolimits_{i}^{2}+\mathop{\sigma}\nolimits^{2}}}\left|{\mathop{v}\nolimits_{i}}\right\rangle-\sum\nolimits_{i=1}^{Q}{\frac{{\mathop{\tilde{\sigma}}\nolimits_{i}}}{{\mathop{\tilde{\sigma}}\nolimits_{i}^{2}+\mathop{\sigma}\nolimits^{2}}}\left|{\mathop{v}\nolimits_{i}}\right\rangle}}}\right\|}\nolimits^{2}\\ \leq\mathop{\left\|{\sum\nolimits_{i=1}^{Q}{\left({\frac{{\mathop{\sigma}\nolimits_{i}}}{{\mathop{\sigma}\nolimits_{i}^{2}+\mathop{\sigma}\nolimits^{2}}}-\frac{{\mathop{\tilde{\sigma}}\nolimits_{i}}}{{\mathop{\tilde{\sigma}}\nolimits_{i}^{2}+\mathop{\sigma}\nolimits^{2}}}}\right)\left|{\mathop{v}\nolimits_{i}}\right\rangle}}\right\|}\nolimits^{2}\leq\mathop{\left({\frac{{\mathop{\kappa}\nolimits_{A}\mathop{\varepsilon}\nolimits_{p}\mathop{\left\|A\right\|}\nolimits_{F}}}{{1+\mathop{\mathop{\sigma}\nolimits^{2}\kappa}\nolimits_{A}^{2}}}}\right)}\nolimits^{2}\end{array}$$

(51)

Let the final error be $\varepsilon$, that is $\left\|{\left|{\hat{\varphi}}\right\rangle-\left|{\mathop{\varphi}\nolimits^{*}}\right\rangle}\right\|<\varepsilon$, we have

$$\mathop{\varepsilon}\nolimits_{p}=\frac{{\varepsilon\left({1+\mathop{\mathop{\sigma}\nolimits^{2}\kappa}\nolimits_{A}^{2}}\right)}}{{\mathop{\mathop{\kappa}\nolimits_{A}\left\|A\right\|}\nolimits_{F}}}$$

(52)

Considering the complexity of the singular vector transition, applying the amplitude amplification technique Brassard et al. (2002), the complexity of implementing $\left|{\hat{\varphi}}\right\rangle$ is estimated as

$$\begin{array}[]{l}{\rm O}\left({\frac{{\mathop{\mathop{\kappa}\nolimits_{A}^{2}\left\|A\right\|}\nolimits_{F}}}{{\varepsilon\left({1+\mathop{\mathop{\sigma}\nolimits^{2}\kappa}\nolimits_{A}^{2}}\right)}}{\mathop{\rm poly}\nolimits}\left({\log MQ}\right)}\right)=\\ {\rm O}\left({\frac{{\mathop{\mathop{\kappa}\nolimits_{A}^{2}\left\|A\right\|}\nolimits_{F}{\mathop{\rm poly}\nolimits}\left({\log MQ}\right)}}{\varepsilon}}\right)\end{array}$$

(53)

Time complexity for VQDME: Let the error of VQDME algorithm be $\varepsilon$. Based on Wang et al. (2020a), we can approximately estimate the complexity of VQDME as

$${\rm O}\left({\frac{{\mathop{T}\nolimits_{VQDME}\mathrm{poly}\left(\log M\right)}}{{\mathop{\varepsilon}\nolimits^{2}}}}\right)$$

(54)

where ${\mathop{T}\nolimits_{VQDME}}$ is the iteration number of the optimization. Therefore, we can obtain the signal subspace with ${\rm O}\left({\frac{{\mathop{T}\nolimits_{VQDME}\mathrm{poly}\left(\log M\right)}}{{\mathop{\varepsilon}\nolimits^{2}}}}\right)$ times quantum-classical hybrid operations.

Time complexity for quantum labeling: First, the searching space can be prepared with the complexity ${\rm O}\left({{\rm{poly}}\log KM}\right)$, and then we can estimate

$$\mathop{P}\nolimits_{S}\left({\left|1\right\rangle}\right)=\frac{1}{K}\sum\nolimits_{n=1}^{K}{\sum\nolimits_{i=1}^{L}{\mathop{\left|{\mathop{\alpha}\nolimits_{i}^{n}}\right|}\nolimits^{2}}}=\Theta\left(1\right)$$

(55)

Moreover, a small amount of sampling are required for the direction estimation satisfying the object (33).

Therefore, Combining the above three quantum subroutines, our MUSIC-based DOA estimation algorithm can be implemented with the complexity

$${\rm O}\left({\frac{{\mathop{\mathop{\kappa}\nolimits_{A}^{2}T_{VQDME}\left\|A\right\|}\nolimits_{F}{\rm{poly}}\log MKQ}}{{\mathop{\varepsilon}\nolimits^{3}}}}\right){\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}$$

(56)

Without loss of generality, the amplitude of the element in $A$ is less than or equal to $1$ and $Q\gg M$, thus $\left\|A\right\|_{F}$ is approximate ${\rm O}\left(Q\right)$.