Introduction

Quantum memories are a fundamental of any global-scale quantum Internet1,2,3,4,5,6. However, while quantum repeaters can be realized without the necessity of quantum memories1,3, these units, in fact, are required for guaranteeing an optimal performance in any high-performance quantum networking scenario3,4,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32. Therefore, the utilization of quantum memories still represents a fundamental problem in the quantum Internet33,34,35,36,37,38,39,40,41,42, since the near-term quantum devices (such as quantum repeaters5,6,8,32,43,44,45,46,47) and gate-model quantum computers48,49,50,51,52,53,54,55,56,57,58,59 have to store the quantum states in their local quantum memories43,44,45,46,47,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84. The main problem here is the efficient readout of the stored quantum systems and the low retrieval efficiency of these systems from the quantum registers of the quantum memory. Currently, no general solution to this problem is available, since the quantum register evolves the stored quantum systems via an unknown operation, and the input quantum system is also unknown, in a general scenario4,5,7,8,9,11,12. The optimization of the readout procedure is therefore a hard and complex problem. Several physical implementations have been developed in the last few years85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105. However, these experimental realizations have several drawbacks, in general because the output signal-to-noise ratio (SNR) values are still not satisfactory for the construction of a powerful, global-scale quantum communication network. As another important application field in quantum communication, the methods of quantum secure direct communication106,107,108,109 also require quantum memory.

Here, we define a novel quantum memory called high-retrieval-efficiency (HRE) quantum memory for near-term quantum devices. An HRE quantum memory unit integrates local unitary operations on its hardware level for the optimization of the readout procedure. An HRE quantum memory unit utilizes the advanced techniques of quantum machine learning to achieve a significant improvement in the retrieval efficiency110,111,112. We define the integrated unitary operations of an HRE quantum memory, prove the learning procedure, and evaluate the achievable output SNR values. The local unitaries of an HRE quantum memory achieve the optimization of the readout procedure in an unsupervised manner without the use of any labeled data or any training sequences. The readout procedure of an HRE quantum memory is realized in a completely blind manner. It requires no information about the input quantum system or about the quantum operation of the quantum register. (It is motivated by the fact that this information is not accessible in any practical setting).

The proposed model assumes that the main challenge is the recovery the stored quantum systems from the quantum register of the quantum memory unit, such that both the input quantum system and the transformation of the quantum memory are unknown. The optimization problem of the readout process also integrates the efficiency of the write-in procedure. In the proposed model, the noise and uncertainty added by the write-in procedure are included in the unknown transformation of the QR quantum register of the quantum memory that results in a σQR mixed quantum system in QR.

The novel contributions of our manuscript are as follows:

  1. 1.

    We define a novel quantum memory called high-retrieval-efficiency (HRE) quantum memory.

  2. 2.

    An HRE quantum memory unit integrates local unitary operations on its hardware level for the optimization of the readout procedure and utilizes the advanced techniques of quantum machine learning.

  3. 3.

    We define the integrated unitary operations of an HRE quantum memory, prove the learning procedure, and evaluate the achievable output signal-to-noise ratio values. We prove that local unitaries of an HRE quantum memory achieve the optimization of the readout procedure in an unsupervised manner without the use of any labeled data or training sequences.

  4. 4.

    We evaluate the retrieval efficiency of an HRE quantum memory and the output SNR.

  5. 5.

    The proposed results are convenient for gate-model quantum computers and near-term quantum devices.

This paper is organized as follows. Section 2 defines the system model and the problem statement. Section 3 evaluates the integrated local unitary operations of an HRE quantum memory. Section 4 proposes the retrieval efficiency in terms of the achievable output SNR values. Finally, Section 5 concludes the results. Supplemental material is included in the Appendix.

System Model and Problem Statement

System model

Let \({\rho }_{in}\) be an unknown input quantum system formulated by n unknown density matrices,

$${\rho }_{in}=\mathop{\sum }\limits_{i=1}^{n}\,{\lambda }_{i}^{(in)}|{\psi }_{i}\rangle \langle {\psi }_{i}|,$$
(1)

where \({\lambda }_{i}^{(in)}\ge 0\), and \({\sum }_{i=1}^{n}\,{\lambda }_{i}^{(in)}=1\).

The input system is received and stored in the QR quantum register of the HRE quantum memory unit. The quantum systems are d-dimensional systems (\(d=2\) for a qubit system). For simplicity, we focus on \(d=2\) dimensional quantum systems throughout the derivations.

The UQR unknown evolution operator of the QR quantum register defines a mixed state σQR as

$$\begin{array}{rcl}{\sigma }_{QR} & = & {U}_{QR}{\rho }_{in}{U}_{QR}^{\dagger }\\ & = & \mathop{\sum }\limits_{i=1}^{n}\,{\lambda }_{i}|{\varphi }_{i}\rangle \langle {\varphi }_{i}|,\end{array}$$
(2)

where \({\lambda }_{i}\ge 0\), \({\sum }_{i=1}^{n}\,{\lambda }_{i}=1\).

Let us allow to rewrite (2) for a particular time t, \(t=1,\ldots ,T\), where T is a total evolution time, via a mixed system \({\sigma }_{QR}^{(t)}\), as

$$\begin{array}{rcl}{\sigma }_{QR}^{(t)} & = & {U}_{QG}^{(t)}{\rho }_{in}{({U}_{QG}^{(t)})}^{\dagger }\\ & = & \mathop{\sum }\limits_{i=1}^{n}\,{\lambda }_{i}^{(t)}|{\varphi }_{i}^{(t)}\rangle \langle {\varphi }_{i}^{(t)}|\\ & = & \mathop{\sum }\limits_{i=1}^{n}\,(\sqrt{{\lambda }_{i}^{(t)}}|{\varphi }_{i}^{(t)}\rangle )(\sqrt{{\lambda }_{i}^{(t)}}\langle {\varphi }_{i}^{(t)}|)\\ & = & \mathop{\sum }\limits_{i=1}^{n}\,{X}_{i}^{(t)}{({X}_{i}^{(t)})}^{\dagger }\\ & = & {X}^{(t)}{({X}^{(t)})}^{\dagger },\end{array}$$
(3)

where \({U}_{QR}^{(t)}\) is an unknown evolution matrix of the QR quantum register at a given t, with a dimension

$${\rm{\dim }}({U}_{QR}^{(t)})={d}^{n}\times {d}^{n},$$
(4)

with \(0\le {\lambda }_{i}^{(t)}\le 1\), \({\sum }_{i}\,{\lambda }_{i}^{(t)}=1\), while \({X}_{i}^{(t)}\in {\mathbb{C}}\) is an unknown complex quantity, defined as

$${X}_{i}^{(t)}=\sqrt{{\lambda }_{i}^{(t)}}|{\varphi }_{i}^{(t)}\rangle $$
(5)

and

$${X}^{(t)}=\mathop{\sum }\limits_{i=1}^{n}\,{X}_{i}^{(t)}.$$
(6)

Then, let us rewrite \({\sigma }_{QR}^{(t)}\) from (3) as

$${\sigma }_{QR}^{(t)}={\rho }_{in}+{\zeta }_{QR}^{(t)},$$
(7)

where \({\rho }_{in}\) is as in (1), and \({\zeta }_{QR}^{(t)}\) is an unknown residual density matrix at a given t.

Therefore, (7) can be expressed as a sum of M source quantum systems,

$${\sigma }_{QR}^{(t)}=\mathop{\sum }\limits_{m=1}^{M}\,{\rho }_{m},$$
(8)

where \({\rho }_{m}\) is the m-th source quantum system and \(m=1,\ldots ,M\), where

$$M=2,$$
(9)

in our setting, since

$${\rho }_{1}={\rho }_{in}$$
(10)

and

$${\rho }_{2}={\zeta }_{QR}^{(t)}.$$
(11)

In terms of the M subsystems, (3) can be rewritten as

$$\begin{array}{rcl}{\sigma }_{QR}^{(t)} & = & \mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{i=1}^{n}\,{\lambda }_{i}^{(m,t)}|{\varphi }_{i}^{(m,t)}\rangle \langle {\varphi }_{i}^{(m,t)}|\\ & = & \mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{i=1}^{n}\,\sqrt{{\lambda }_{i}^{(m,t)}}|{\varphi }_{i}^{(m,t)}\rangle \sqrt{{\lambda }_{i}^{(m,t)}}\langle {\varphi }_{i}^{(m,t)}|\\ & = & \mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{i=1}^{n}\,{X}_{i}^{(m,t)}{({X}_{i}^{(m,t)})}^{\dagger }\\ & = & \mathop{\sum }\limits_{m=1}^{M}\,{X}^{(m,t)}{({X}^{(m,t)})}^{\dagger },\end{array}$$
(12)

where \({X}_{i}^{(m,t)}\) is a complex quantity associated with an m-th source system,

$${X}_{i}^{(m,t)}=\sqrt{{\lambda }_{i}^{(m,t)}}|{\varphi }_{i}^{(m,t)}\rangle ,$$
(13)

with \(0\le {\lambda }_{i}^{(m,t)}\le 1\), \({\sum }_{m}\,{\sum }_{i}\,{\lambda }_{i}^{(m,t)}=1\), and

$${X}^{(m,t)}=\mathop{\sum }\limits_{i=1}^{n}\,{X}_{i}^{(m,t)}.$$
(14)

The aim is to find the VQG inverse matrix of the unknown evolution matrix UQR in (2), as

$${V}_{QG}={U}_{QG}^{-1},$$
(15)

that yields the separated readout quantum system of the HRE quantum memory unit for \(t=1,\ldots ,T\), such that for a given t,

$${\sigma }_{out}^{(t)}={V}_{QG}^{(t)}{\sigma }_{QR}^{(t)}{({V}_{QG}^{(t)})}^{\dagger },$$
(16)

where

$${V}_{QG}^{(t)}={({U}_{QG}^{(t)})}^{-1}.$$
(17)

For a total evolution time T, the target σout density matrix is yielded at the output of the HRE quantum memory unit, as

$${\sigma }_{out}\approx \mathop{\sum }\limits_{i=1}^{n}\,{\lambda }_{i}^{(in)}|{\psi }_{i}\rangle \langle {\psi }_{i}|$$
(18)

with a sufficiently high SNR value,

$${\rm{SNR}}({\sigma }_{out})\ge x,$$
(19)

where x is an SNR value that depends on the actual physical layer attributes of the experimental implementation.

The problem is therefore that both the input quantum system (1) and the transformation matrix UQR in (2) of the quantum register are unknown. As we prove, by integrating local unitaries to the HRE quantum memory unit, the unknown evolution matrix of the quantum register can be inverted, which allows us to retrieve the quantum systems of the quantum register. The retrieval efficiency will be also defined in a rigorous manner.

Problem statement

The problem statement is as follows.

Let M be the number of source systems in the QR quantum register such that the sum of the M source systems identifies the mixed state of the quantum register. Let m be the index of the source system, \(m=1,\ldots ,M\), such that \(m=1\) identifies the unknown input quantum system stored in the quantum register (target source system), while \(m=2,\ldots ,M\) are some unknown residual quantum systems. The input quantum system, the residual systems, and the transformation operation of the quantum register are unknown. The aim is then to define local unitary operations to be integrated on the HRE quantum memory unit for an HRE readout procedure in an unsupervised manner with unlabeled data.

The problems to be solved are summarized in Problems 1–4.

Problem 1.

Find an unsupervised quantum machine learning method, UML, for the factorization of the unknown mixed quantum system of the quantum register via a blind separation of the unlabeled quantum register. Decompose the unknown mixed system state into a basis unitary and a residual quantum system.

Problem 2.

Define a unitary operation for partitioning the bases with respect to the source systems of the quantum register.

Problem 3.

Define a unitary operation for the recovery of the target source system.

Problem 4.

Evaluate the retrieval efficiency of the HRE quantum memory in terms of the achievable SNR.

The resolutions of the problems are proposed in Theorems 1–4.

The schematic model of an HRE quantum memory unit is depicted in Fig. 1.

Figure 1
figure 1

The schematic model of a high-retrieval-efficiency (HRE) quantum memory unit. The HRE quantum memory unit contains a QR quantum register and integrated local unitary operations. The \(n\) input quantum systems, \({\rho }_{1}\ldots {\rho }_{n}\), are received and stored in the quantum register. The state of the QR quantum register defines a mixed state, \({\sigma }_{QR}={\sum }_{i}\,{\lambda }_{i}{\rho }_{i}\), where \({\sum }_{i}\,{\lambda }_{i}=1\). The stored density matrices of the \(QR\) quantum register are first transformed by a \({U}_{ML}\), a quantum machine learning unitary (depicted by the orange-shaded box) that implements an unsupervised learning for a blind separation of the unlabeled input, and decomposable as \({U}_{ML}={U}_{F}{U}_{CQT}{U}_{P}{U}_{CQT}^{\dagger }\), where \({U}_{F}\) is a factorization unitary, \({U}_{CQT}\) is the quantum constant \(Q\) transform with a windowing function \({f}_{W}\) for the localization of the wave functions of the quantum register, \({U}_{P}\) is a basis partitioning unitary, while \({U}_{CQT}^{\dagger }\) is the inverse of \({U}_{CQT}\). The result of \({U}_{ML}\) is processed further by the \({\tilde{U}}_{{\rm{DSTFT}}}^{\dagger }\) unitary (depicted by the green-shaded box) that realizes the inverse quantum discrete short-time Fourier transform (DSTFT) operation (depicted by the yellow-shaded box), and by the \({U}_{DFT}\) (quantum discrete Fourier transform) unitary to yield the desired output \({\rho ^{\prime} }_{1}\ldots {\rho ^{\prime} }_{n}\).

The procedures realized by the integrated unitary operations of the HRE quantum memory are depicted in Fig. 2.

Figure 2
figure 2

Detailed procedures of an HRE quantum memory. The unknown input quantum system is stored in the \(QR\) quantum register that realizes an unknown transformation. The density matrix of the quantum register is the sum of \(M=2\) source systems, where source system \(m=1\) identifies the valuable unknown input quantum system stored in the quantum register, while \(m=2\) identifies an unknown undesired residual quantum system. The \({U}_{F}\) unitary evaluates \(K\) bases for the source system and defines a \(W\) auxiliary quantum system. The \({U}_{CQT}\) unitary is a preliminary operation for the partitioning of the \(K\) bases onto \(M\) clusters via unitary \({U}_{P}\). The \({U}_{P}\) unitary regroups the bases with respect to the \(M=2\) source systems. The results are then processed by the \({\tilde{U}}_{{\rm{DSTFT}}}^{\dagger }\) and \({U}_{DFT}\) unitaries to extract the source system \(m=1\) on the output of the memory unit.

Experimental implementation

An experimental implementation of an HRE quantum memory in a near-term quantum device52 can integrate standard photonics devices, optical cavities and other fundamental physical devices. The quantum operations can be realized via the framework of gate-model quantum computations of near-term quantum devices52,53,54,55,56, such as superconducting units53. The application of a HRE quantum memory in a quantum Internet setting1,2,4,5,6 can be implemented via noisy quantum links between the quantum repeaters8,32,43,44,45,46,47 (e.g., optical fibers7,62,113, wireless quantum channels27,28, free-space optical channels114) and fundamental quantum transmission protocols24,115,116,117.

Integrated Local Unitaries

This section defines the local unitary operations integrated on an HRE quantum memory unit.

Quantum machine learning unitary

The UML quantum machine learning unitary implements an unsupervised learning for a blind separation of the unlabeled quantum register. The UML unitary is defined as

$${U}_{ML}={U}_{F}{U}_{CQT}{U}_{P}{U}_{CQT}^{\dagger },$$
(20)

where UF is a factorization unitary, UCQT is the quantum constant Q transform, UP is a partitioning unitary, while \({U}_{CQT}^{\dagger }\) is the inverse of \({U}_{CQT}\).

Factorization unitary

Theorem 1.

(Factorization of the unknown mixed quantum system of the quantum register). The UF unitary factorizes the unknown \({\sigma }_{QR}\) mixed quantum system of the QR quantum register into a unitary \({u}_{mk}={e}^{-i{H}_{mk}\tau /\hslash }\), with a Hamiltonian \({H}_{mk}\) and application time \(\tau \), and into a system \({w}_{kt}\), where \(t=1,\ldots ,T\), \(m=1,\ldots ,M\), and \(k=1,\ldots ,K\), and where T is the evolution time, M is the number of source systems of \({\sigma }_{QR}\), and K is the number of bases.

Proof. The aim of the UF factorization unitary is to factorize the mixed quantum register (2) into a basis matrix UB and a quantum system \({\overrightarrow{\rho }}_{W}\), as

$$\begin{array}{rcl}{U}_{F}{\overrightarrow{\sigma }}_{QR}{U}_{F}^{\dagger } & = & {U}_{F}({U}_{QR}{\rho }_{in}{U}_{QR}^{\dagger }){U}_{F}^{\dagger }\\ & = & {U}_{B}{\overrightarrow{\rho }}_{W}{U}_{B}^{\dagger },\end{array}$$
(21)

where UB is a complex basis matrix, defined as

$${U}_{B}=\{{u}_{mk}\}\in {{\mathbb{C}}}^{M\times K},$$
(22)

and \({\overrightarrow{\rho }}_{W}\in {{\mathbb{C}}}^{K\times T}\) is a complex matrix, defined as

$${\overrightarrow{\rho }}_{W}={\{{\rho }_{W}^{(t)}\}}_{t=1}^{T},$$
(23)

where

$$\begin{array}{rcl}{\rho }_{W}^{(t)} & = & \mathop{\sum }\limits_{k=1}^{K}\,{v}_{k}^{(t)}|{\phi }_{k}\rangle \langle {\phi }_{k}|\\ & = & \mathop{\sum }\limits_{k=1}^{K}\,\sqrt{{v}_{k}^{(t)}}|{\phi }_{k}\rangle \sqrt{{v}_{k}^{(t)}}\langle {\phi }_{k}|\\ & = & \mathop{\sum }\limits_{k=1}^{K}\,{W}_{k}^{(t)}{({W}_{k}^{(t)})}^{\dagger },\end{array}$$
(24)

where \(0\le {v}_{k}^{(t)}\le 1\), and \({\sum }_{k=1}^{K}\,{v}_{k}^{(t)}=1\), while \(K\) is the total number of bases of \({U}_{B}\), while \({W}_{k}^{(t)}\in {\mathbb{C}}\) is a complex quantity, as

$${W}_{k}^{(t)}=\sqrt{{v}_{k}^{(t)}}|{\phi }_{k}\rangle .$$
(25)

The first part of the problem is therefore to find (22), where \({u}_{mk}\) is a unitary that sets a computational basis for \({W}_{k}^{(t)}\) in (25), defined as

$${u}_{mk}={e}^{-i{H}_{mk}\tau /\hslash },$$
(26)

where Hmk is a Hamiltonian, as

$${H}_{mk}={G}_{mk}|{k}_{m}\rangle \langle {k}_{m}|,$$
(27)

where Gmk is the eigenvalue of basis \(|{k}_{m}\rangle \), \({H}_{mk}|{k}_{m}\rangle ={G}_{mk}|{k}_{m}\rangle \), while \(\tau \) is the application time of umk.

The second part of the problem is to determine W, as

$$W=\{{W}_{k}^{(t)}={w}_{kt}\}\in {{\mathbb{C}}}^{K\times T},$$
(28)

where \({W}_{k}^{(t)}={w}_{kt}\) is a system state, that formulates \({\tilde{X}}^{(m,t)}\) as

$$\begin{array}{rcl}{\tilde{X}}^{(m,t)} & = & {[{U}_{B}W]}_{mt}\\ & = & \mathop{\sum }\limits_{k=1}^{K}\,{u}_{mk}{w}_{kt},\end{array}$$
(29)

where \({\tilde{X}}^{(m,t)}\) is an approximation of \({X}^{(m,t)}\),

$${\tilde{X}}^{(m,t)}\approx {X}^{(m,t)},$$
(30)

where \({X}^{(m,t)}\) is defined in (14).

As follows, for the total evolution time T, \(\overrightarrow{X}\in {{\mathbb{C}}}^{M\times T}\) can be defined as

$$\overrightarrow{X}={\{{X}^{(1,t)},\ldots ,{X}^{(M,t)}\}}_{t=1}^{T},$$
(31)

and the challenge is to evaluate (31) as a decomposition

$$\begin{array}{rcl}\tilde{X} & = & {U}_{B}W\\ & = & {e}^{-i{H}_{\Sigma }\tau /\hslash }W\\ & = & \mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{t=1}^{T}\,\mathop{\sum }\limits_{k=1}^{K}\,{u}_{mk}{W}_{k}^{(t)}\\ & = & \mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{t=1}^{T}\,\mathop{\sum }\limits_{k=1}^{K}\,{e}^{-i{H}_{mk}\tau /\hslash }\sqrt{{v}_{k}^{(t)}}|{\phi }_{k}^{(m,t)}\rangle .\end{array}$$
(32)

Thus, by applying of the umk unitaries for the total evolution time T, \(\tilde{X}\in {{\mathbb{C}}}^{M\times T}\) is as

$$\begin{array}{rcl}\tilde{X} & = & {U}_{B}W\\ & = & \mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{k=1}^{K}\,{\ell }_{m,k}(\tau )|{k}_{m}\rangle ,\\ & = & \alpha (\mathop{\sum }\limits_{{k}_{1}=1}^{{K}_{1}}\,|{k}_{1}\rangle +\ldots +\mathop{\sum }\limits_{{k}_{M}=1}^{{K}_{M}}\,|{k}_{M}\rangle ),\end{array}$$
(33)

where Km is the number of bases associated with the m-th source system,

$$\mathop{\sum }\limits_{m=1}^{K}\,{K}_{m}=K,$$
(34)

and \(0\le {|{\ell }_{m,k}(\tau )|}^{2}\le 1\), \({\sum }_{m=1}^{M}\,{\sum }_{k=1}^{K}\,{|{\ell }_{m,k}(\tau )|}^{2}=1\).

In our setting \(M=2\), and our aim is to get the system state \(m=1\) on the output of the HRE quantum memory, thus a \(|{\Phi }^{\ast }\rangle \) target output system state is defined as

$$|{\Phi }^{\ast }\rangle =\frac{1}{\sqrt{{K}_{1}}}\,\mathop{\sum }\limits_{{k}_{1}=1}^{{K}_{1}}\,|{k}_{1}\rangle ,$$
(35)

where K1 is the number of bases for source system \(m=1\), \({k}_{1}=1,\ldots ,{K}_{1}\).

Let rewrite the system state \(\tilde{X}\) (32) as

$$\tilde{X}={\{{\tilde{X}}^{(1,t)},\ldots ,{\tilde{X}}^{(M,t)}\}}_{t=1}^{T},$$
(36)

and let

$${X}^{(t)}=\mathop{\sum }\limits_{m=1}^{M}\,{X}^{(m,t)},$$
(37)

and

$${\tilde{X}}^{(t)}=\mathop{\sum }\limits_{m=1}^{M}\,{\tilde{X}}^{(m,t)}.$$
(38)

Then, let \({\rho }_{\overrightarrow{X}}\) be a density matrix associated with \(\overrightarrow{X}\), defined as

$${\rho }_{\overrightarrow{X}}=\mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{t=1}^{T}\,{\overrightarrow{X}}^{(m,t)}{({\overrightarrow{X}}^{(m,t)})}^{\dagger }$$
(39)

and let

$${\rho }_{\tilde{X}}=\mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{t=1}^{T}\,{\tilde{X}}^{(m,t)}{({\tilde{X}}^{(m,t)})}^{\dagger }$$
(40)

be the density matrix associated with (36).

The aim of the estimation is to minimize the \(D(\cdot \parallel \cdot )\) quantum relative entropy function taken between \({\rho }_{\overrightarrow{X}}\) and \({\rho }_{\tilde{X}}\), thus an \(f({U}_{F})\) objective function for \({U}_{F}\) is defined via (37) and (38) as

$$\begin{array}{rcl}f({U}_{F}) & = & \mathop{{\rm{\min }}}\limits_{\tilde{X}}\,D({\rho }_{\overrightarrow{X}}\parallel {\rho }_{\tilde{X}})\\ & = & \mathop{{\rm{\min }}}\limits_{\tilde{X}}\,{\rm{Tr}}({\rho }_{\overrightarrow{X}}\,\log \,({\rho }_{\overrightarrow{X}}))-{\rm{Tr}}({\rho }_{\overrightarrow{X}}\,\log \,({\rho }_{\tilde{X}})).\end{array}$$
(41)

To achieve the objective function \(f({U}_{F})\) in (41), a factorization method is defined for \({U}_{F}\) that is based on the fundamentals of Bayesian nonnegative matrix factorization118,119,120,121,122,123,124,125,126,127 (Footnote: The \({U}_{F}\) factorization unitary applied on the mixed state of the quantum register is analogous to a Poisson-Exponential Bayesian nonnegative matrix factorization118,119,120,121 process). The method adopts the Poisson distribution as \( {\mathcal L} (\,\cdot \,)\) likelihood function and the exponential distribution for the control parameters118,119,120,121 \({\alpha }_{mk}\) and \({\beta }_{kt}\) defined for the controlling of \({u}_{mk}\) and \({w}_{kt}\).

Let \({u}_{mk}\) and \({w}_{kt}\) from (29) be defined via the control parameters \({\alpha }_{mk}\) and \({\beta }_{kt}\) as exponential distributions

$${u}_{mk}\simeq {\alpha }_{mk}{e}^{-{\alpha }_{mk}{u}_{mk}},$$
(42)

with mean \({\alpha }_{mk}^{-1}\), and

$${w}_{kt}\simeq {\beta }_{kt}{e}^{-{\beta }_{kt}{w}_{kt}},$$
(43)

with mean \({\beta }_{kt}^{-1}\).

Using (41), (42) and (43), a \( {\mathcal L} (\,\cdot \,)\) log likelihood function

$$- {\mathcal L} (x,y|z)=-\,\log \,{\rm{\Pr }}(x,y|z)$$
(44)

can be defined as

$$\begin{array}{l}- {\mathcal L} ({U}_{B},W|\overrightarrow{X})\\ \begin{array}{rcl} & = & D({\rho }_{\overrightarrow{X}}\parallel {\rho }_{\tilde{X}})+\mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{k=1}^{K}\,{\alpha }_{mk}{u}_{mk}{({\alpha }_{mk}{u}_{mk})}^{\dagger }\\ & & +\,\mathop{\sum }\limits_{k=1}^{K}\,\mathop{\sum }\limits_{t=1}^{T}\,{\beta }_{kt}{w}_{kt}{({\beta }_{kt}{w}_{kt})}^{\dagger },\end{array}\end{array}$$
(45)

thus the objective function \(f({U}_{F})\) can be rewritten via as (45)

$$f({U}_{F})=\mathop{{\rm{\min }}}\limits_{\tilde{X}}\,(- {\mathcal L} ({U}_{B},W|\overrightarrow{X}))\mathrm{}.$$
(46)

The problem is therefore can be reduced to determine the model parameters

$$\zeta =\{{U}_{B},W\}$$
(47)

that are treated as latent variables for the estimation of the control parameters118,119,120,121,125,126,127

$${\tau }_{mk}^{(t)}=\{{\alpha }_{mk},{\beta }_{kt}\}.$$
(48)

A maximum likelihood estimation \(\tilde{\zeta }\) of (47) is as

$$\tilde{\zeta }={\rm{\arg }}\,\mathop{{\rm{\max }}}\limits_{\zeta }\, {\cal{D}} (\overrightarrow{X}|\zeta ),$$
(49)

where \( {\cal{D}} (\,\cdot \,)\) is some distribution, that identifies an incomplete estimation problem.

The estimation of (47) can also be yielded from a maximization of a marginal likelihood function \( {\mathcal L} (\overrightarrow{X}|\zeta )\) as

$$ {\mathcal L} (\overrightarrow{X}|\zeta )=\int \int \,\sum _{\overrightarrow{\kappa }}\, {\cal{D}} (\overrightarrow{X}|\overrightarrow{\kappa })\, {\cal{D}} (\overrightarrow{\kappa }|{U}_{B},W)\, {\cal{D}} ({U}_{B},W|\zeta )\,d{U}_{B}dW,$$
(50)

where \(\overrightarrow{\kappa }\) is a complex matrix, \(\overrightarrow{\kappa }\in {{\mathbb{C}}}^{M\times T}\),

$$\overrightarrow{\kappa }={\{{\kappa }^{(1,t)},\ldots ,{\kappa }^{(M,t)}\}}_{t=1}^{T},$$
(51)

where

$${\kappa }^{(m,t)}={({\kappa }_{k=1}^{(m,t)},\ldots ,{\kappa }_{k=K}^{(m,t)})}^{T},$$
(52)

with

$${\kappa }_{k}^{(m,t)}={\kappa }_{mkt}$$
(53)

where

$${\kappa }_{mkt}={u}_{mk}{w}_{kt}.$$
(54)

The quantity in (54) can be estimated via (42) and (43) as

$${\kappa }_{mkt}\approx {\alpha }_{mk}{e}^{-{\alpha }_{mk}{u}_{mk}}{\beta }_{kt}{e}^{-{\beta }_{kt}{w}_{kt}}.$$
(55)

Using (54), \({\tilde{X}}^{(m,t)}\) in (29) can be rewritten as

$${\tilde{X}}^{(m,t)}=\mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{k=1}^{K}\,{\kappa }_{mkt}.$$
(56)

However, since the exact solution does not exists118,119,120,121, since it would require the factorization of \( {\cal{D}} (\overrightarrow{\kappa },{U}_{B},W|\overrightarrow{X},\zeta )\), such that \(\zeta ,{U}_{B},W\) are unknown.

This problem can be solved by a variational Bayesian inference procedure118,119,120,121,125,126,127, via the maximization of the lower bound of a likelihood function \({ {\mathcal L} }_{{ {\cal{D}} }_{v}}\)

$$\begin{array}{rcl}{ {\mathcal L} }_{{ {\cal{D}} }_{v}} & = & \iint \,\sum _{\overrightarrow{\kappa }}\,{ {\cal{D}} }_{v}(\overrightarrow{\kappa },{U}_{B},W)\,\log \,\tfrac{ {\cal{D}} (\overrightarrow{X},\overrightarrow{\kappa },{U}_{B},W|\zeta )}{{ {\cal{D}} }_{v}(\overrightarrow{\kappa },{U}_{B},W)}d{U}_{B}dW\\ & = & {\mathbb{E}}(\log \, {\cal{D}} (\overrightarrow{X},\overrightarrow{\kappa },{U}_{B},W|\zeta ))+H({ {\cal{D}} }_{v}(\overrightarrow{\kappa },{U}_{B},W)),\end{array}$$
(57)

where \({ {\cal{D}} }_{v}\) is a variational distribution, while \(H({ {\cal{D}} }_{v}(\overrightarrow{\kappa },{U}_{B},W))\) is the entropy of variational distribution \({ {\cal{D}} }_{v}(\overrightarrow{\kappa },{U}_{B},W)\),

$$H({ {\cal{D}} }_{v}(\overrightarrow{\kappa },{U}_{B},W))=\mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{t=1}^{T}\,H({\kappa }^{(m,t)})+\mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{k=1}^{K}\,H({u}_{mk})+\mathop{\sum }\limits_{k=1}^{K}\,\mathop{\sum }\limits_{t=1}^{T}\,H({w}_{kt}),$$
(58)

and where \({ {\cal{D}} }_{v}(\overrightarrow{\kappa },{U}_{B},W)\) is a joint variational distribution, as

$$\begin{array}{rcl}{ {\cal{D}} }_{v}(\overrightarrow{\kappa },{U}_{B},W) & = & { {\cal{D}} }_{v}(\overrightarrow{\kappa }){ {\cal{D}} }_{v}({U}_{B}){ {\cal{D}} }_{v}(W)\\ & = & \prod _{m}\,\prod _{t}\,\prod _{k}\,{ {\cal{D}} }_{v}({\kappa }_{mkt}){ {\cal{D}} }_{v}({u}_{mk}){ {\cal{D}} }_{v}({w}_{kt}),\end{array}$$
(59)

from which distribution \( {\cal{D}} (\overrightarrow{\kappa },{U}_{B},W|\overrightarrow{X},\zeta )\) can be approximated as118,119,120,121

$$ {\cal{D}} (\overrightarrow{\kappa },{U}_{B},W|\overrightarrow{X},\zeta )\approx \prod _{m}\,\prod _{t}\,\prod _{k}\,{ {\cal{D}} }_{v}({\kappa }_{mkt}){ {\cal{D}} }_{v}({u}_{mk}){ {\cal{D}} }_{v}({w}_{kt}).$$
(60)

The function \({ {\mathcal L} }_{{ {\cal{D}} }_{v}}\) in (57) is related to (50) as

$$ {\mathcal L} (\overrightarrow{X}|\zeta )\ge { {\mathcal L} }_{{ {\cal{D}} }_{v}}.$$
(61)

The result in (59) therefore also determines the number \(K\) of bases selected for the factorization unitary \({U}_{F}\). The \({ {\cal{D}} }_{v}\) variational distributions \({ {\cal{D}} }_{v}({\kappa }_{mkt})\), \({ {\cal{D}} }_{v}({u}_{k})\) and \({ {\cal{D}} }_{v}({w}_{kt})\) are determined for the unitary UF as follows.

Let \({ {\cal{D}} }_{v}(\Phi )\) refer to the variational distribution of a given \(\Phi \),

$$\Phi \in \{\overrightarrow{\kappa },{U}_{B},W\}.$$
(62)

Since only the joint (posterior) distribution \( {\cal{D}} (\overrightarrow{X},\overrightarrow{\kappa },{U}_{B},W|\zeta )\) is obtainable, the variational distributions have to be evaluated as

$${{\mathbb{E}}}_{{ {\cal{D}} }_{v}(i\ne \Phi )}(\log \, {\cal{D}} (\overrightarrow{X},\overrightarrow{\kappa },{U}_{B},W|\zeta ))=\,\log \,{ {\cal{D}} }_{v}(\Phi ),$$
(63)

where \({{\mathbb{E}}}_{{ {\cal{D}} }_{v}(i\ne \Phi )}(\,\cdot \,)\) is the expectation function of the \({ {\cal{D}} }_{v}(i)\) variational distribution of i, such that \(i\ne \Phi \), where \(\Phi \) is as in (62), with

$${{\mathbb{E}}}_{a}(f(a)+g(a))={{\mathbb{E}}}_{a}(f(a))+{{\mathbb{E}}}_{a}(g(a)),$$
(64)

for some functions \(f(a)\) and \(g(a)\), and

$${{\mathbb{E}}}_{a}(bf(a))=b{{\mathbb{E}}}_{a}(f(a))$$
(65)

for some constant b, (note: for simplicity, we use \({\mathbb{E}}(\,\cdot \,)\) for the expectation function), while

$$\begin{array}{l}\log \, {\cal{D}} (\overrightarrow{X},\overrightarrow{\kappa },{U}_{B},W|\zeta )\\ \begin{array}{rcl} & = & \mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{t=1}^{T}\,\log \,{f}_{\delta }({X}^{(m,t)}-\mathop{\sum }\limits_{k=1}^{K}\,{\kappa }_{mkt})+\mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{k=1}^{K}\,\mathop{\sum }\limits_{t=1}^{T}\,({\kappa }_{mkt}\,\log ({u}_{mk}{w}_{kt})\\ & & -\,{u}_{mk}{w}_{kt}-\,\log \,{f}_{\Gamma }({\kappa }_{mkt}+1))+\mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{k=1}^{K}\,(\log \,{\alpha }_{mk}-{\alpha }_{mk}{u}_{mk})\\ & & +\,\mathop{\sum }\limits_{k=1}^{K}\,\mathop{\sum }\limits_{t=1}^{T}\,(\log \,{\beta }_{kt}-{\beta }_{kt}{w}_{kt}),\end{array}\end{array}$$
(66)

where \({f}_{\delta }(\,\cdot \,)\) is the Dirac delta function, while \({f}_{\Gamma }(\,\cdot \,)\) is the Gamma function,

$${f}_{\Gamma }(x)={\int }_{0}^{\infty }\,{t}^{x-1}{e}^{-t}dt.$$
(67)

By utilizing a variational Poisson–Exponential Bayesian learning118,119,120,121, these variational distributions can be evaluated as follows.

The \({ {\cal{D}} }_{v}({\kappa }_{mkt})\) variational distribution is as

$${ {\cal{D}} }_{v}({\kappa }_{mkt})= {\mathcal M} ({\kappa }_{mkt}|{\eta }_{mkt})$$
(68)

where \( {\mathcal M} \) is a multinomial distribution, while \({\eta }_{mkt}\) is a multinomial parameter

$${\eta }_{mkt}=\frac{{e}^{{\mathbb{E}}(\log {u}_{mk})+{\mathbb{E}}(\log {w}_{kt})}}{{\sum }_{j}\,{e}^{{\mathbb{E}}(\log {u}_{mj})+{\mathbb{E}}(\log {w}_{jt})}},$$
(69)

while the \({ {\cal{D}} }_{v}({\kappa }^{(m,t)})\) variational distribution is as

$$\begin{array}{l} {\mathcal M} ({\kappa }^{(m,t)}|{X}^{(m,t)},{\eta }_{k}^{(m,t)})\\ \,=\,{f}_{\delta }({X}^{(m,t)}-\mathop{\sum }\limits_{k=1}^{K}{\kappa }_{mkt}){X}^{(m,t)}!\prod _{k}\,\tfrac{{({\eta }_{mkt})}^{{\kappa }_{mkt}}}{{\kappa }_{mkt}!},\end{array}$$
(70)

where \({\eta }_{k}^{(m,t)}\) is a multinomial parameter vector

$${\eta }_{k}^{(m,t)}={({\eta }_{k=1}^{(m,t)},\ldots ,{\eta }_{k=K}^{(m,t)})}^{T},$$
(71)

such that

$$\mathop{\sum }\limits_{k=1}^{K}\,{\eta }_{k}^{(m,t)}=1.$$
(72)

The \({ {\cal{D}} }_{v}({u}_{mk})\) variational distribution is as

$$\begin{array}{l}{ {\cal{D}} }_{v}({u}_{mk})\\ \,=\,{e}^{(\mathop{\sum }\limits_{t=1}^{T}{\mathbb{E}}({\kappa }_{mkt})\log {u}_{mk}-(\mathop{\sum }\limits_{t=1}^{T}{\mathbb{E}}({w}_{kt})+{\alpha }_{mk}){u}_{mk})}\\ \,=\,{\mathscr G}({u}_{mk}|{\tilde{\alpha }}_{mk}(A),{\tilde{\alpha }}_{mk}(B)),\end{array}$$
(73)

where \({\mathscr G}(\,\cdot \,)\) is a Gamma distribution,

$${\mathscr G}(x;a,b)={e}^{(a-1)\log x-\frac{x}{b}-\log {f}_{\Gamma }(a)-a\log b},$$
(74)

where a is a shape parameter, while b is a scale parameter, \({f}_{\Gamma }(\,\cdot \,)\) is the Gamma function (67). The entropy of (74) is as

$$H({\mathscr G}(x;a,b))=-\,(a-1){\partial }_{{{\mathscr G}}_{\log }}(a)+\,\log \,b+a+\,\log \,{f}_{\Gamma }(a),$$
(75)

where \({\partial }_{{{\mathscr G}}_{\log }}(\,\cdot \,)\) is the derivative of the log gamma function (digamma function),

$${\partial }_{{{\mathscr G}}_{\log }}(x)=\frac{d\,\log \,{f}_{\Gamma }(x)}{dx},$$
(76)

while \({\mathbb{E}}({\kappa }_{mkt})\) is evaluated as

$${\mathbb{E}}({\kappa }_{mkt})={X}^{(m,t)}{\eta }_{mkt},$$
(77)

while \({\tilde{\alpha }}_{mk}(A)\) and \({\tilde{\alpha }}_{mk}(B)\) are control parameters for \({U}_{B}\), defined as

$${\tilde{\alpha }}_{mk}(A)=1+\mathop{\sum }\limits_{t=1}^{T}\,{\mathbb{E}}({\kappa }_{mkt}),$$
(78)

while \({\tilde{\alpha }}_{mk}(B)\) is defined as

$${\tilde{\alpha }}_{mk}(B)=\frac{1}{{\sum }_{t=1}^{T}\,{\mathbb{E}}({w}_{kt})+{\alpha }_{mk}}.$$
(79)

The \({ {\cal{D}} }_{v}({w}_{kt})\) variational distribution is as

$$\begin{array}{l}{ {\cal{D}} }_{v}({w}_{kt})\\ \,=\,{e}^{(\mathop{\sum }\limits_{m=1}^{M}{\mathbb{E}}({\kappa }_{mkt})\log {w}_{kn}-(\mathop{\sum }\limits_{m=1}^{M}{\mathbb{E}}({u}_{mk})+{\beta }_{kt}){w}_{kt})}\\ \,=\,{\mathscr G}({w}_{kt}|{\tilde{\beta }}_{kt}(A),{\tilde{\beta }}_{kt}(B)),\end{array}$$
(80)

where \({\tilde{\beta }}_{kt}(A)\) and \({\tilde{\beta }}_{kt}(B)\) are control parameters for W, defined as

$${\tilde{\beta }}_{kt}(A)=1+\mathop{\sum }\limits_{m=1}^{M}\,{\mathbb{E}}({\kappa }_{mkt}),$$
(81)

and

$${\tilde{\beta }}_{kt}(B)=\frac{1}{{\sum }_{m=1}^{M}\,{\mathbb{E}}({u}_{mk})+{\beta }_{kt}}.$$
(82)

Given the variational parameters \({\tilde{\alpha }}_{mk}(A)\), \({\tilde{\alpha }}_{mk}(B)\), \({\tilde{\beta }}_{kt}(A)\) and \({\tilde{\beta }}_{kt}(B)\) in (78), (79), (81) and (82), the estimates of \({U}_{B}\) and \(W\) are realized by the determination of the Gamma means \({\mathbb{E}}({u}_{mk})\) and \({\mathbb{E}}({w}_{kt})\)118,119,120,121. It can be verified that the mean \({\mathbb{E}}({w}_{kt})\) in (73), (79) and (80) can be evaluated via (81) and (82) as a mean of a Gamma distribution

$${\mathbb{E}}({w}_{kt})={\tilde{\beta }}_{kt}(A){\tilde{\beta }}_{kt}(B),$$
(83)

while \({\mathbb{E}}(\log \,{w}_{kt})\) is as

$${\mathbb{E}}(\log \,{w}_{kt})={\partial }_{{{\mathscr G}}_{\log }}({\tilde{\beta }}_{kt}(A))+\,\log \,{\tilde{\beta }}_{kt}(B),$$
(84)

where \({\partial }_{{{\mathscr G}}_{\log }}(\,\cdot \,)\) digamma function (76).

The mean \({\mathbb{E}}({u}_{mk})\) in (80) and (82) can be evaluated via (78) and (79), as a mean of a Gamma distribution

$${\mathbb{E}}({u}_{mk})={\tilde{\alpha }}_{mk}(A){\tilde{\alpha }}_{mk}(B),$$
(85)

and \({\mathbb{E}}(\log \,{u}_{mk})\) is yielded as

$${\mathbb{E}}(\log \,{u}_{mk})={\partial }_{{{\mathscr G}}_{\log }}({\tilde{\alpha }}_{mk}(A))+\,\log \,{\tilde{\alpha }}_{mk}(B).$$
(86)

As the \({ {\cal{D}} }_{v}({\kappa }_{mkt})\), \({ {\cal{D}} }_{v}({u}_{mk})\) and \({ {\cal{D}} }_{v}({w}_{kt})\) variational distributions are determined via (68), (73) and (80) the evaluation of (59) is straightforward.

Using the defined terms, the term \({\mathbb{E}}(\log \, {\cal{D}} (\overrightarrow{X},\overrightarrow{\kappa },{U}_{B},W|\zeta ))\) from (57) can be evaluated as

$$\begin{array}{l}{\mathbb{E}}(\log \, {\cal{D}} (\overrightarrow{X},\overrightarrow{\kappa },{U}_{B},W|\zeta ))\\ \begin{array}{rcl} & = & \mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{t=1}^{T}\,{\mathbb{E}}(\log \,{f}_{\delta }({X}^{(m,t)}-\mathop{\sum }\limits_{k=1}^{K}\,{\kappa }_{mkt}))\\ & & +\,\mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{k=1}^{K}\,{\mathbb{E}}(\log \,{u}_{mk})\,\mathop{\sum }\limits_{t=1}^{T}\,{\mathbb{E}}({\kappa }_{mkt})\\ & & +\,\mathop{\sum }\limits_{k=1}^{K}\,\mathop{\sum }\limits_{t=1}^{T}\,{\mathbb{E}}(\log \,{w}_{kt})\,\mathop{\sum }\limits_{m=1}^{M}\,{\mathbb{E}}({\kappa }_{mkt})-\mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{k=1}^{K}\,\mathop{\sum }\limits_{t=1}^{T}\,{\mathbb{E}}({u}_{mk})\,{\mathbb{E}}({w}_{kt})\\ & & -\,\mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{k=1}^{K}\,\mathop{\sum }\limits_{t=1}^{T}\,{\mathbb{E}}(\log \,{f}_{\Gamma }({\kappa }_{mkt}+1))\\ & & +\,\mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{k=1}^{K}\,(\log \,{\alpha }_{mk}-{\alpha }_{mk}{\mathbb{E}}({u}_{mk}))\\ & & +\,\mathop{\sum }\limits_{k=1}^{K}\,\mathop{\sum }\limits_{t=1}^{T}\,(\log \,{\beta }_{kt}-{\beta }_{kt}{\mathbb{E}}({w}_{kt})),\end{array}\end{array}$$
(87)

while the \(H({ {\cal{D}} }_{v}(\overrightarrow{\kappa },{U}_{B},W))\) entropy of the variational distribution from (58) can be evaluated as

$$\begin{array}{l}H({ {\cal{D}} }_{v}(\overrightarrow{\kappa },{U}_{B},W))\\ \begin{array}{rcl} & = & \mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{t=1}^{T}\,(\,-\,\log \,{f}_{\Gamma }({X}^{(m,t)}+1)-\mathop{\sum }\limits_{k=1}^{K}\,{\mathbb{E}}({\kappa }_{mkt})\,\log \,{\eta }_{mkt})\\ & & +\,\mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{k=1}^{K}\,\mathop{\sum }\limits_{t=1}^{T}\,{\mathbb{E}}(\log \,{f}_{\Gamma }({\kappa }_{mkt}+1))\\ & & -\,\mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{t=1}^{T}\,{\mathbb{E}}(\log \,{f}_{\delta }({X}^{(m,t)}-\mathop{\sum }\limits_{k=1}^{K}\,{\kappa }_{mkt}))\\ & & +\,\mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{k=1}^{K}\,(-({\tilde{\alpha }}_{mk}(A)-1){\partial }_{{{\mathscr G}}_{\log }}({\tilde{\alpha }}_{mk}(A))+\,\log ({\tilde{\alpha }}_{mk}(B))\\ & & +\,{\tilde{\alpha }}_{mk}(A)+\,\log \,{f}_{\Gamma }({\tilde{\alpha }}_{mk}(A)))\\ & & +\,\mathop{\sum }\limits_{k=1}^{K}\,\mathop{\sum }\limits_{t=1}^{T}\,(-({\tilde{\beta }}_{kt}(A)-1){\partial }_{{{\mathscr G}}_{\log }}({\tilde{\beta }}_{k}(A))+\,\log ({\tilde{\beta }}_{k}(B))\\ & & +\,{\tilde{\beta }}_{k}(A)+\,\log \,{f}_{\Gamma }({\tilde{\beta }}_{k}(A))).\end{array}\end{array}$$
(88)

Thus, from (87) and (88), the lower bound \({ {\mathcal L} }_{{ {\cal{D}} }_{v}}\)in (57) is as

$$\begin{array}{lcl}{ {\mathcal L} }_{{{\mathscr{D}}}_{v}} & = & -\,\mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{t=1}^{T}\,\mathop{\sum }\limits_{k=1}^{K}\,{\mathbb{E}}({u}_{mk})\,{\mathbb{E}}({w}_{kt})\\ & & +\,\mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{t=1}^{T}\,(\,-\,\log \,{f}_{\Gamma }({X}^{(m,t)}+1)-\mathop{\sum }\limits_{k=1}^{K}\,{\mathbb{E}}({\kappa }_{mkt})\,\log \,{\eta }_{mkt})\\ & & +\,\mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{k=1}^{K}\,{\mathbb{E}}(\log \,{u}_{mk})\,\mathop{\sum }\limits_{t=1}^{T}\,{\mathbb{E}}({\kappa }_{mkt})\\ & & +\,\mathop{\sum }\limits_{k=1}^{K}\,\mathop{\sum }\limits_{t=1}^{T}\,{\mathbb{E}}(\log \,{w}_{kt})\,\mathop{\sum }\limits_{m=1}^{M}\,{\kappa }_{mkt}\\ & & +\,\mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{k=1}^{K}\,(\log \,{\alpha }_{mk}-{\alpha }_{mk}{\mathbb{E}}({u}_{mk}))\\ & & +\,\mathop{\sum }\limits_{k=1}^{K}\,\mathop{\sum }\limits_{t=1}^{T}\,(\log \,{\beta }_{kt}-{\beta }_{kt}{\mathbb{E}}({w}_{kt}))\\ & & +\,\mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{k=1}^{K}\,(-({\tilde{\alpha }}_{mk}(A)-1){\partial }_{{{\mathscr G}}_{\log }}({\tilde{\alpha }}_{mk}(A))+\,\log \,{\tilde{\alpha }}_{mk}(B)\\ & & +\,{\tilde{\alpha }}_{mk}(A)+\,\log \,{f}_{\Gamma }({\tilde{\alpha }}_{mk}(A)))\\ & & +\,\mathop{\sum }\limits_{k=1}^{K}\,\mathop{\sum }\limits_{t=1}^{T}\,(-({\tilde{\beta }}_{kt}(A)-1){\partial }_{{{\mathscr G}}_{\log }}({\tilde{\beta }}_{kt}(A))+\,\log \,{\tilde{\beta }}_{kt}(B)\\ & & +\,{\tilde{\beta }}_{kt}(A)+\,\log \,{f}_{\Gamma }({\tilde{\beta }}_{kt}(A))).\end{array}$$
(89)

The next problem is the \({\tilde{\tau }}_{k}^{(t)}\) estimation of the control parameters \({\alpha }_{mk},{\beta }_{kt}\) in (48) as

$${\tilde{\tau }}_{mk}^{(t)}=\{{E}_{mk},{F}_{kt}\},$$
(90)

such that \({E}_{mk}\) is a basis estimation

$${E}_{mk}\approx {\alpha }_{mk}$$
(91)

and \({F}_{kt}\) is a system estimation

$${F}_{kt}\approx {\beta }_{kt},$$
(92)

such that the variational lower bound \({ {\mathcal L} }_{{ {\cal{D}} }_{v}}\) in (89) is maximized118,119,120,121. It is achieved for the unitary \({U}_{F}\) as follows. The maximization problem can be formalized via the \(\partial ({ {\mathcal L} }_{{ {\cal{D}} }_{v}})\) derivative of \({ {\mathcal L} }_{{ {\cal{D}} }_{v}}\)

$$\frac{\partial ({ {\mathcal L} }_{{ {\cal{D}} }_{v}})}{\partial {\alpha }_{mk}}=\frac{1}{{\alpha }_{mk}}-{\mathbb{E}}({u}_{mk})+\frac{\partial (\log \,({\tilde{\alpha }}_{mk}(B)))}{\partial {\alpha }_{mk}}=0,$$
(93)

and

$$\frac{\partial ({ {\mathcal L} }_{{ {\cal{D}} }_{v}})}{\partial {\beta }_{kt}}=\frac{1}{{\beta }_{kt}}-{\mathbb{E}}({w}_{kt})+\frac{\partial (\log ({\tilde{\beta }}_{kt}(B)))}{\partial {\beta }_{kt}}=0,$$
(94)

which is solvable via118,120

$${({\alpha }_{mk})}^{2}+\mathop{\sum }\limits_{t=1}^{T}\,{\mathbb{E}}({w}_{kt}){\alpha }_{mk}-\frac{{\sum }_{t=1}^{T}\,{\mathbb{E}}({w}_{kt})}{{\mathbb{E}}({u}_{mk})}=0,$$
(95)

and

$${({\beta }_{kt})}^{2}+\mathop{\sum }\limits_{m=1}^{M}\,{\mathbb{E}}({u}_{mk}){\beta }_{kt}-\frac{{\sum }_{m=1}^{M}\,{\mathbb{E}}({u}_{mk})}{{\mathbb{E}}({w}_{kt})}=0.$$
(96)

After some calculations, \({E}_{mk}\) and \({F}_{kt}\) from (90) are as

$${E}_{mk}=\frac{1}{2}(\,-\,\mathop{\sum }\limits_{t=1}^{T}\,{\mathbb{E}}({w}_{kt})+{({(\mathop{\sum }\limits_{t=1}^{T}{\mathbb{E}}({w}_{kt}))}^{2}+4\frac{{\sum }_{t=1}^{T}{\mathbb{E}}({w}_{kt})}{{\mathbb{E}}({u}_{mk})})}^{\frac{1}{2}}),$$
(97)

and

$${F}_{kt}=\frac{1}{2}(\,-\,\mathop{\sum }\limits_{m=1}^{M}\,{\mathbb{E}}({u}_{mk})+{({(\mathop{\sum }\limits_{m=1}^{M}{\mathbb{E}}({u}_{mk}))}^{2}+4\frac{{\sum }_{m=1}^{M}{\mathbb{E}}({u}_{mk})}{{\mathbb{E}}({w}_{kt})})}^{\frac{1}{2}}),$$
(98)

respectively.

From (97) and (98), the \({\tilde{\tau }}_{mk}^{(t)}\) estimation in (90) is therefore straightforwardly yielded. Therefore, using the parameters \({\tilde{\alpha }}_{mk}(B),{\tilde{\alpha }}_{mk}(A),{\tilde{\beta }}_{kt}(A),{\tilde{\beta }}_{kt}(B)\) and \({\eta }_{mkt}\), the optimal variational distributions \({ {\cal{D}} }_{v}({\kappa }_{mkt})\), \({ {\cal{D}} }_{v}({u}_{mk})\) and \({ {\cal{D}} }_{v}({w}_{kt})\) can be substituted to estimate \({\tilde{\tau }}_{mk}^{(t)}\).

Using (97) and (98), the estimation of terms \({u}_{k}\) (42), \({w}_{kt}\) (43) and \({\kappa }_{kt}\) (55) are yielded as

$${\tilde{u}}_{mk}={E}_{mk}{e}^{-{E}_{mk}{\tilde{u}}_{mk}},$$
(99)
$${\tilde{w}}_{kt}={F}_{kt}{e}^{-{F}_{kt}{\tilde{w}}_{kt}},$$
(100)

and

$${\tilde{\kappa }}_{mkt}={E}_{mk}{e}^{-{E}_{mk}{\tilde{u}}_{mk}}{F}_{kt}{e}^{-{F}_{kt}{\tilde{w}}_{kt}}.$$
(101)

The evaluation of (97) and (98) therefore is yielded in an iterative manner through the \({\tilde{\alpha }}_{mk}(B)\), \({\tilde{\alpha }}_{mk}(A)\), \({\tilde{\beta }}_{kt}(A)\), \({\tilde{\beta }}_{kt}(B)\) and \({\eta }_{mkt}\), and the K* optimal number of bases, K, is determined with respect to (89) such that

$${K}^{\ast }={\rm{\arg }}\,\mathop{{\rm{\max }}}\limits_{K}\,{ {\mathcal L} }_{{ {\cal{D}} }_{v}}(K),$$
(102)

where \({ {\mathcal L} }_{{ {\cal{D}} }_{v}}(K)\) refers to \({ {\mathcal L} }_{{ {\cal{D}} }_{v}}\) from (89) at a particular base number \(K\).

The proof is concluded here. ■

The schematic representation of unitary \({U}_{F}\) is depicted in Fig. 3.

Figure 3
figure 3

Representation of the \({U}_{F}\) unitary over a total evolution time \(t\), with \(K\) factored bases and \(M\) source systems (\(M=2\) in our setting). The factorization is represented by the solid-line arrows. At a given \(t\), \(t=1,\ldots ,T\), the input system of \({U}_{F}\) subject of factorization is \({X}^{(m,t)}=\sqrt{{\lambda }_{i}^{(m,t)}}|{\varphi }_{i}^{(m,t)}\rangle \), \(m=1,\ldots ,M\). Term \({\kappa }_{mkt}\) is expressed as \({\kappa }_{mkt}={u}_{mk}{w}_{kt}\), where \({u}_{mk}={e}^{-i{H}_{mk}\tau /\hslash }\) is a unitary, \({u}_{mk}\in {\mathbb{C}}\), \(k=1,\ldots ,K\), which sets a computational basis for \({w}_{kt}\), \({w}_{kt}={W}_{k}^{(t)}=\sqrt{{v}_{k}^{(t)}}|{\phi }_{k}\rangle \). The basis matrix is \({U}_{B}=\{{u}_{mk}\}\in {{\mathbb{C}}}^{M\times K}\) with \(K\) bases, \({H}_{mk}={G}_{mk}|{k}_{m}\rangle \langle {k}_{m}|\) is a Hamiltonian, and \(W=\{{W}_{k}^{(t)}={w}_{kt}\}\in {{\mathbb{C}}}^{K\times T}\), \({w}_{kt}\in {\mathbb{C}}\). The factorization decomposes \({X}^{(m,t)}\) into \({X}^{(m,t)}={[{U}_{B}\overrightarrow{W}]}_{mt}\), and for the total evolution \(\overrightarrow{X}={U}_{B}\overrightarrow{W}\), where \(\overrightarrow{X}={\{{X}^{(1,t)},\ldots ,{X}^{(M,t)}\}}_{t=1}^{T}\), while \({\kappa }_{kt}\) is as \({\kappa }_{mkt}={u}_{mk}{w}_{kt}\). Terms \({\alpha }_{mk}\) and \({\beta }_{kt}\) are control parameters for \({u}_{mk}\) and \({w}_{kt}\) (controlling is depicted by the dashed-line arrows) to evaluate the parameters as \({u}_{mk}\simeq {\alpha }_{mk}{e}^{-{\alpha }_{mk}{u}_{mk}}\) and \({w}_{kt}\simeq {\beta }_{kt}{e}^{-{\beta }_{kt}{w}_{kt}}\), estimated by \({E}_{mk}\) and \({F}_{kt}\) as \({\tilde{u}}_{mk}={E}_{mk}{e}^{-{E}_{mk}{\tilde{u}}_{mk}}\) and \({\tilde{w}}_{kt}={F}_{kt}{e}^{-{F}_{kt}{\tilde{w}}_{kt}}\).

Quantum constant Q transform

As the \(\{{\tilde{u}}_{mk}\}\) basis estimations (99) are determined via \(\{{E}_{mk}\}\) (97), the next problem is the partitioning of the \(K\) bases with respect to \(M\), see (8). To achieve the partitioning, first the bases of \({U}_{B}\) are transformed by the \({U}_{CQT}\) is the quantum constant \(Q\) transform128. The \({U}_{CQT}\) operation is similar to the discrete QFT (quantum Fourier transform) transform117, and defined in the following manner.

The \({U}_{CQT}\) transform is defined as

$${U}_{CQT}(|k\rangle ,m)=\frac{1}{\sqrt{N}}\,\mathop{\sum }\limits_{j=0}^{N-1}\,{f}_{W}(j-m){e}^{2\pi ijQ/m}|j\rangle =|{\varphi }_{k}\rangle ,$$
(103)

where \(|k\rangle \) is a quantum state of the computational basis \(B\), and in the current setting

$$N=K,$$
(104)

and

$$|k\rangle ={E}_{mk},$$
(105)

thus \(B\) is as

$$B:\{|0\rangle ,\ldots ,|K-1\rangle \},$$
(106)

while \(h\) is selected such that

$$0\le (j-h)\le N-1=K-1$$
(107)

holds, and \(Q\) is defined via the following relation

$$\frac{2\pi k}{K}=\frac{2\pi Q}{h},$$
(108)

from which \(Q\) is yielded at a given \(h\), \(k\) and \(K\), as

$$Q=\frac{hk}{K},$$
(109)

while \({f}_{W}(\,\cdot \,)\) is a windowing function129 that localizes the wavefunctions of the quantum register, defined via parameter \(h\) as

$${f}_{W}(j-h)=\frac{1}{2}(1-\,\cos \,(\frac{2\pi (h-m)}{K-1})).$$
(110)

(Footnote: The function in (110) is the so-called Hanning window129).

The \(|{\varphi }_{k}\rangle \) output states of \({U}_{CQT}\) therefore identify a set \({\mathscr S}_{\varphi }\) of states, as

$${\mathscr S}_{\varphi }:\{|{\varphi }_{k}\rangle :k=0,\ldots ,K-1\}$$
(111)

that formulates an orthonormal basis.

The \({U}_{CQT}^{\dagger }\) inverse of \({U}_{CQT}\) will be processed as the \({U}_{P}\) partitioning is completed, with the same \({f}_{W}(\,\cdot \,)\) windowing function, defined as

$${U}_{CQT}^{\dagger }(|k\rangle ,h)=\frac{1}{\sqrt{K}}\,\mathop{\sum }\limits_{j=0}^{K-1}\,{f}_{W}(j-h){e}^{-2\pi ijQ/h}|j\rangle .$$
(112)

Applying (103) on the \(K\) estimated bases \(\{{E}_{mk}\}\) yields the \({C}_{B}\) transformed bases, as

$$\begin{array}{rcl}{C}_{B} & = & {U}_{CQT}({U}_{B})\\ & = & \{{C}_{mk}\}\in {{\mathbb{C}}}^{M\times K},\end{array}$$
(113)

where \({C}_{mk}\) is as,

$${C}_{mk}={U}_{CQT}({E}_{mk}).$$
(114)

After the application of (113), the resulting system is therefore as

$${C}_{B}W=({U}_{CQT}{U}_{B})W,$$
(115)

where \({C}_{B}W\in {{\mathbb{C}}}^{M\times T}\).

Basis partitioning unitary

Theorem 2.

(Partitioning the bases of source systems). The \(Q\) transformed bases can be partitioned to \(M\) partitions via the \({U}_{P}\) partitioning unitary operation.

Proof. As the \({U}_{CQT}\) transforms of the \(\{{E}_{mk}\}\) basis estimations (99) are determined via \({C}_{B}\) (113), the \(Q\) transformed bases are partitioned to \(M\) partitions via the \({U}_{P}\) unitary operation, as follows.

Let the system state from (115) be denoted by

$$S={C}_{B}W$$
(116)

and let \(\tilde{S}\) be the estimation of \(S\)130, defined as

$$\tilde{S}=\langle \langle {\mathcal R} {\mathcal E} \rangle {\mathcal H} \rangle ,$$
(117)

where

$${\mathscr T}\in \{ {\mathcal R} , {\mathcal E} , {\mathcal H} \}$$
(118)

is a tensor (multidimensional array)131,132 with dimension \({\rm{\dim }}({\mathscr T})\), and size

$$s({\mathscr T})=\mathop{\prod }\limits_{i=1}^{{\rm{\dim }}({\mathscr T})}\,|{d}_{i}({\mathscr T})|,$$
(119)

where \(|{d}_{i}({\mathscr T})|\) is the size of the i-th dimension \({d}_{i}({\mathscr T})\).

Let

$$ {\mathcal R} ={\mathscr A}\circ {\mathcal B} $$
(120)

be a translation tensor of size

$$\begin{array}{rcl}s( {\mathcal R} ) & = & \mathop{\prod }\limits_{i=1}^{{\rm{\dim }}( {\mathcal R} )}\,|{d}_{i}( {\mathcal R} )|\\ & = & \mathop{\prod }\limits_{i=1}^{{\rm{\dim }}({\mathscr A})}\,|{d}_{i}({\mathscr A})|\times \mathop{\prod }\limits_{i=1}^{{\rm{\dim }}( {\mathcal B} )}\,|{d}_{i}( {\mathcal B} )|,\end{array}$$
(121)

with

$${\rm{\dim }}( {\mathcal R} )=3,$$
(122)

as

$$|{d}_{1}( {\mathcal R} )|=M,$$
(123)
$$|{d}_{2}( {\mathcal R} )|=1,$$
(124)

and

$$|{d}_{3}( {\mathcal R} )|=M$$
(125)

and let

$$ {\mathcal E} ={\mathscr A}\circ {\mathscr C}$$
(126)

be a tensor of size

$$\begin{array}{rcl}s( {\mathcal E} ) & = & \mathop{\prod }\limits_{i=1}^{{\rm{\dim }}( {\mathcal E} )}\,|{d}_{i}( {\mathcal E} )|\\ & = & \mathop{\prod }\limits_{i=1}^{{\rm{\dim }}({\mathscr A})}\,|{d}_{i}({\mathscr A})|\times \mathop{\prod }\limits_{i=1}^{{\rm{\dim }}({\mathscr C})}\,|{d}_{i}({\mathscr C})|,\end{array}$$
(127)

with

$${\rm{\dim }}( {\mathcal E} )=2$$
(128)

as

$$|{d}_{1}( {\mathcal E} )|=M,$$
(129)
$$|{d}_{2}( {\mathcal E} )|=K,$$
(130)

and with

$${\rm{\dim }}( {\mathcal H} )=3,$$
(131)

as

$$|{d}_{1}( {\mathcal R} )|=1,$$
(132)
$$|{d}_{2}( {\mathcal R} )|=K$$
(133)

and

$$|{d}_{3}( {\mathcal R} )|=T,$$
(134)

thus

$${\rm{\dim }}({\mathscr A})=M$$
(135)

and

$${\rm{\dim }}( {\mathcal B} )=M$$
(136)

while

$${\rm{\dim }}({\mathscr C})=K.$$
(137)

The term \(\langle {\mathcal R} {\mathcal E} \rangle \) is evaluated as

$$\begin{array}{l}{\langle {\mathcal R} {\mathcal E} \rangle }_{\{1:{\rm{\dim }}({\mathscr A}),1:{\rm{\dim }}({\mathscr A})\}}({j}_{1},\ldots ,{j}_{{\rm{\dim }}( {\mathcal B} )},{k}_{1},\ldots ,{k}_{{\rm{\dim }}({\mathscr C})})\\ \begin{array}{rcl} & = & \mathop{\sum }\limits_{{i}_{1}=1}^{{d}_{1}({\mathscr A})}\cdots \mathop{\sum }\limits_{{i}_{{\rm{\dim }}({\mathscr A})}=1}^{{d}_{{\rm{\dim }}({\mathscr A})}({\mathscr A})}\, {\mathcal R} ({i}_{1},\ldots ,{i}_{{\rm{\dim }}({\mathscr A})},{j}_{1},\ldots ,{j}_{{\rm{\dim }}( {\mathcal B} )})\\ & & \times \, {\mathcal E} ({i}_{1},\ldots ,{i}_{{\rm{\dim }}({\mathscr A})},{k}_{1},\ldots ,{k}_{{\rm{\dim }}({\mathscr C})}),\end{array}\end{array}$$
(138)

where \( {\mathcal R} (i,j)\) is the indexing for the elements of the tensor.

Let \( {\mathcal E} (\forall m,k)\) refer to the j-th column of \( {\mathcal E} \), and let \( {\mathcal H} (1,k,\forall t)\) refer to the j-th lateral slice of \( {\mathcal H} \). Then, let be a \({U}_{P}\) unitary operation that achieves the decomposition of (117) with respect to a given \(k\), \(k=1,\ldots ,K\), as

$${[S]}_{k}=\langle \langle {\mathcal R} {\mathcal E} (\forall m,k)\rangle {\mathcal H} (1,k,\forall t)\rangle $$
(139)

with a particular cost function \(f({U}_{P})\) of the \({U}_{P}\) unitary defined via the quantum relative entropy function, as

$$\begin{array}{rcl}f({U}_{P}) & = & \mathop{{\rm{\min }}}\limits_{\tilde{S}}\,D({\rho }_{S}\parallel \tilde{S}{\tilde{S}}^{\dagger })\\ & = & \mathop{{\rm{\min }}}\limits_{\tilde{S}}\,{\rm{Tr}}({\rho }_{S}\,\log \,({\rho }_{S}))-{\rm{Tr}}({\rho }_{S}\,\log \,(\tilde{S}{\tilde{S}}^{\dagger })),\end{array}$$
(140)

where \({\rho }_{S}\) is the density matrix associated with \(S\) is as in (116),

$${\rho }_{S}={U}_{P}{U}_{CQT}{U}_{F}(\mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{t=1}^{T}\,{\overrightarrow{X}}^{(m,t)}{({\overrightarrow{X}}^{(m,t)})}^{\dagger }),$$
(141)

while \(\tilde{S}\) is given in (117).

Using (139), the \(Q\)-transformed bases are partitioned into \(M\) classes, the partition Ω outputted by \({U}_{P}\) is evaluated as

$$\Omega ={\rm{\arg }}\,\mathop{{\rm{\max }}}\limits_{k}\,({[{\rm{Q}}]}_{k}),$$
(142)

where Q is a \(1\times K\) size matrix, such that

$${[{\rm{Q}}]}_{k}=\mathop{\sum }\limits_{m=1}^{M}\,\langle \langle {\mathcal R} (\forall ,1,\forall ) {\mathcal E} (\forall ,k)\rangle {\mathcal H} (1,k,\forall )\rangle .$$
(143)

Since \(M=2\) in our setting, the partition (142) can be rewritten as

$$\Omega ={\Omega }_{Q}^{(1)}+{\Omega }_{Q}^{(2)},$$
(144)

where \({\Omega }_{Q}^{(m)}\) identifies a cluster of \({K}_{m}\) Q-transformed bases for m-th system state,

$${\Omega }_{Q}^{(m)}={\{{\Omega }_{Q}^{(m,{k}_{m})}\}}_{{k}_{m}=1}^{{K}_{m}},$$
(145)

of

$$|{\Omega }_{Q}^{(m)}|={K}_{m}$$
(146)

bases formulated via the base estimations (99) for the m-th system state in (8), such that

$$\mathop{\sum }\limits_{m=1}^{M}\,{K}_{m}=K.$$
(147)

Since the partitioning is made over the Q transformed bases, the output of \({U}_{P}\) is then transformed by the \({U}_{CQT}^{\dagger }\) inverse transformation (112). ■

Inverse quantum constant Q transform

Applying the \({U}_{CQT}^{\dagger }\) inverse transformation (112) on the partitions (143) of the Q transformed bases yields the decomposition of the bases of \({U}_{B}\) onto \(M\) classes, as

$${U}_{CQT}^{\dagger }(\Omega )=\theta =\mathop{\sum }\limits_{m=1}^{M}\,{\gamma }^{(m)},$$
(148)

and since \(M=2\)

$$\theta ={\gamma }^{(1)}+{\gamma }^{(2)},$$
(149)

where \({\gamma }^{(m)}\) identifies a cluster of \({K}_{m}\) bases for m-th system state.

Therefore, the resulting system state is as

$$\begin{array}{l}{U}_{CQT}^{\dagger }({U}_{P}({C}_{B}W))\\ \begin{array}{rcl} & = & {U}_{CQT}^{\dagger }({U}_{P}{U}_{CQT}{U}_{B})W\\ & = & \chi W.\end{array}\end{array}$$
(150)

The next problem is therefore the evaluation of the estimations of the \(M=2\) source systems \({\rho }_{in}\) and \({\zeta }_{QR}^{(t)}\), as given in (7) from \(\chi W\). Using the system state (150), the system separation is produced by the \({U}_{{\rm{DSTFT}}}^{\dagger }\) unitary that realizes the inverse quantum DSTFT (discrete short-time Fourier transform)129.

Inverse quantum DSTFT and quantum DFT

The result of unitary \({U}_{ML}\) is evaluated further by the \({U}_{{\rm{DSTFT}}}^{\dagger }\) unitary.

Theorem 3.

(Target source system recovery). Source system \(m=1\) can be extracted by the \({U}_{{\rm{DSTFT}}}^{\dagger }\) and \({U}_{DFT}\) discrete quantum Fourier transform on the output of an HRE quantum memory.

Proof. The \({U}_{{\rm{DSTFT}}}^{\dagger }\) inverse quantum DSTFT transformation applied to a state \(|k\rangle \) of the computational basis

$$B:\{|0\rangle ,\ldots ,|K-1\rangle \},$$
(151)

is defined as

$${U}_{{\rm{DSTFT}}}^{\dagger }(|k\rangle ,h)=\frac{1}{\sqrt{K}}\,\mathop{\sum }\limits_{j=0}^{K-1}\,{f}_{W}(j-h){e}^{-2\pi ijk/K}|j\rangle =|{\psi }_{k}\rangle ,$$
(152)

where \(h\) is selected such that

$$0\le (j-h)\le K-1$$
(153)

holds, set

$${\mathscr S}_{\psi }:\{|{\psi }_{k}\rangle :k=0,\ldots ,K-1\}$$
(154)

formulates an new orthonormal basis, while \({f}_{W}(\,\cdot \,)\) is a windowing function129.

Using system state \(\chi W\) in (150), let \({\gamma }^{(m,k)}\) be a k-th basis of cluster \({\gamma }^{(m)}\), and let \({(\chi W)}^{(m,t)}\) be defined as

$${(\chi W)}^{(m,t)}={[\chi W]}_{mk}=\mathop{\sum }\limits_{k=1}^{K}\,{\gamma }^{(m,k)}{W}_{k}^{(m,t)}$$
(155)

and let system \(|\chi W\rangle \) identify (33) as

$$|\chi W\rangle =\alpha \,\mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{{k}_{m}=1}^{{K}_{m}}\,|{k}_{m}\rangle ,$$
(156)

where \(|{k}_{m}\rangle \) is the eigenvector of the Hamiltonian of \({\gamma }^{(m,{k}_{m})}\), \({K}_{m}\) is the cardinality of cluster \({\gamma }^{(m)}\), while \({\sum }_{m\mathrm{=1}}^{M}{\sum }_{{k}_{m}\mathrm{=1}}^{{K}_{m}}\alpha \mathrm{=1}\).

Since the \(|{k}_{1}\rangle \) values are some parameters of \({U}_{ML}\), we can redefine (156) as

$$|\chi W\rangle =\alpha \,\mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{{k}_{m}=1}^{{K}_{m}}\,|{k}_{1}+{x}_{m,{k}_{m}}\rangle ,$$
(157)

where

$${x}_{m,{k}_{m}}=\{\begin{array}{ll}0, & {\rm{if}}\,m=1\\ !0, & {\rm{otherwise}}\end{array},$$
(158)

and

$$\alpha =\frac{1}{\sqrt{K}}.$$
(159)

In our setting, using \({k}_{m=1}\) as input parameter available from the \({U}_{ML}\) block, we redefine the formula of (152) via a unitary \({\tilde{U}}_{{\rm{DSTFT}}}^{\dagger }\), as

$${\tilde{U}}_{{\rm{DSTFT}}}^{\dagger }(|{k}_{m}\rangle ,h)=\frac{1}{\sqrt{K}}\,\mathop{\sum }\limits_{j=0}^{K-1}\,{f}_{W}(j-h){e}^{-2\pi ij{k}_{1}/K}|j\rangle =|{\psi }_{{k}_{m}}\rangle ,$$
(160)

where we set \({f}_{W}(j-h)\) to unity,

$${f}_{W}(j-h)=1.$$
(161)

Thus, applying (160) on (157) yields

$$\begin{array}{l}{\tilde{U}}_{{\rm{DSTFT}}}^{\dagger }(\alpha \,\mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{{k}_{m}=1}^{{K}_{m}}\,|{k}_{1}+{x}_{m,{k}_{m}}\rangle )\\ \begin{array}{rcl} & = & \tfrac{1}{\sqrt{K}}\,\mathop{\sum }\limits_{j=0}^{K-1}\,((\alpha \,\mathop{\sum }\limits_{m=0}^{M-1}\,\mathop{\sum }\limits_{{k}_{m}=0}^{{K}_{m}-1}\,{e}^{-2\pi ij{k}_{1}/K}){e}^{-2\pi ij{x}_{m,{k}_{m}}/K})|j\rangle \\ & = & \tfrac{1}{\sqrt{K}}(\mathop{\sum }\limits_{j=0}^{K-1}\,{e}^{-2\pi ij{x}_{m,{k}_{m}}/K})(\alpha \,\mathop{\sum }\limits_{m=0}^{M-1}\,\mathop{\sum }\limits_{{k}_{m}=0}^{{K}_{m}-1}\,{e}^{-2\pi ij{k}_{1}/K})|j\rangle ,\end{array}\end{array}$$
(162)

where

$$j=\frac{K}{{k}_{m}},$$
(163)

and \({\sum }_{j=0}^{K-1}\,{e}^{-2\pi ij{x}_{m,{k}_{m}}/K}=1\), thus (162) can be rewritten as

$$\begin{array}{l}{\tilde{U}}_{{\rm{DSTFT}}}^{\dagger }(\alpha \mathop{\sum }\limits_{m=1}^{M}\mathop{\sum }\limits_{{k}_{m}=1}^{{K}_{m}}|{k}_{1}+{x}_{m,{k}_{m}}\rangle )\\ \,=\,\tfrac{1}{\sqrt{K}}(\alpha \,\mathop{\sum }\limits_{m=0}^{M-1}\,\mathop{\sum }\limits_{{k}_{m}=0}^{{K}_{m}-1}\,{e}^{-2\pi i(\tfrac{K}{{k}_{m}}){k}_{1}/K})|\tfrac{K}{{k}_{m}}\rangle .\end{array}$$
(164)

As follows, if

$$j=\frac{K}{{k}_{1}},$$
(165)

then, the resulting \({\rm{\Pr }}(j)\) probability is

$$\begin{array}{rcl}{\rm{\Pr }}(j) & = & \tfrac{1}{K}{|\alpha \mathop{\sum }\limits_{k=0}^{{K}_{1}-1}{e}^{-2\pi ij{k}_{1}/K}|}^{2}\\ & = & \tfrac{1}{K}{|\alpha \mathop{\sum }\limits_{k=0}^{{K}_{1}-1}{e}^{-2\pi i\tfrac{K}{{k}_{1}}{k}_{1}/K}|}^{2}\\ & = & \tfrac{1}{K}{|\alpha |}^{2}{K}_{1}^{2}\\ & = & \tfrac{1}{{K}^{2}}{K}_{1}^{2},\end{array}$$
(166)

while for the remaining j-s, the probabilities are vanished out, thus

$${\rm{\Pr }}(j)=\mathrm{0,}$$
(167)

if

$$j\ne \frac{K}{{k}_{1}}.$$
(168)

Therefore, applying the \({U}_{DFT}\) discrete quantum Fourier transform on the resulting system state (164), defined in our setting as

$${U}_{DFT}(|k\rangle )=\frac{1}{\sqrt{{K}_{1}}}\,\mathop{\sum }\limits_{j=0}^{{K}_{1}-1}\,{e}^{2\pi ijk/{K}_{1}}|j\rangle ,$$
(169)

yields the source system \(m=1\) in terms of the K1 bases, as

$$\begin{array}{rcl}{U}_{DFT}{\tilde{U}}_{{\rm{DSTFT}}}^{\dagger }(\alpha \,\mathop{\sum }\limits_{m=1}^{M}\,\mathop{\sum }\limits_{{k}_{m}=1}^{{K}_{m}}\,|{k}_{1}+{x}_{m,{k}_{m}}\rangle ) & = & \tfrac{1}{\sqrt{{K}_{1}}}\,\mathop{\sum }\limits_{{k}_{m}=1}^{{K}_{1}}\,|{k}_{1}\rangle \\ & = & |{\Phi }^{\ast }\rangle ,\end{array}$$
(170)

that identifies the target system from (35).

The proof is concluded here. ■

The state of the \(QR\) quantum register after the \({\tilde{U}}_{{\rm{CQT}}}^{\dagger }\) operation and after the \({\tilde{U}}_{{\rm{DSTFT}}}^{\dagger }\) operation is depicted in Fig. 4.

Figure 4
figure 4

(a) The state of the \(QR\) quantum register after the \({\tilde{U}}_{{\rm{CQT}}}^{\dagger }\) operation. The quantum register contains \(K={\sum }_{m}\,{K}_{m}\) states, \(|{k}_{1}+{x}_{m,{k}_{m}}\rangle \), each with probability \({|\alpha |}^{2}=1/K\), with a unit distance between the states (depicted by the red dots). (b) The state of the \(QR\) quantum register after the \({\tilde{U}}_{{\rm{DSTFT}}}^{\dagger }\) operation. The quantum register contains \({K}_{1}\) quantum states, \(|\frac{K}{{k}_{1}}\rangle \), \({k}_{1}=0,\ldots ,{K}_{1}-1\), each with probability \(\frac{1}{K}{|\alpha |}^{2}{K}_{1}^{2}=\frac{1}{{K}^{2}}{K}_{1}^{2}\), with a distance \(\frac{K}{{K}_{1}^{2}-{K}_{1}}\) between the states (depicted by the red dots; the vanished-out states of the quantum register are depicted by the black dots).

Retrieval Efficiency

This section evaluates the retrieval efficiency of an HRE quantum memory in terms of the achievable output SNR values.

Theorem 4.

(Retrieval efficiency of an HRE quantum memory). The SNR of the output quantum system of an HRE quantum memory is evolvable from the difference of the wave function energy ratios taken between the input system, the quantum register system, and the output quantum system.

Proof. Let \(|{\psi }_{in}\rangle \) be an arbitrary quantum system fed into the input of an HRE quantum memory unit,

$$|{\psi }_{in}\rangle =\sum _{i}\,{a}_{i}|i\rangle ,$$
(171)

and let \(|\varphi \rangle \) be the state outputted from the QR quantum register,

$$|\varphi \rangle ={U}_{QR}|{\psi }_{in}\rangle ,$$
(172)

where UQG is an unknown transformation.

Let \(|{\Phi }^{\ast }\rangle \) be the output system of as given in (170), that can be rewritten as

$$|{\Phi }^{\ast }\rangle =U|\varphi \rangle =U({U}_{QR}|{\psi }_{in}\rangle ),$$
(173)

where U is the operator of the integrated unitary operations of the HRE quantum memory, defined as

$$U={U}_{ML}{\tilde{U}}_{{\rm{DSTFT}}}^{\dagger }{U}_{DFT}={U}_{F}{U}_{CQT}{U}_{P}{U}_{CQT}^{\dagger }{\tilde{U}}_{{\rm{DSTFT}}}^{\dagger }{U}_{DFT}.$$
(174)

Then, let \({{\mathcal O}}_{V}\) be a verification oracle that computes the energy E of a wavefunction \(|\psi \rangle ={\sum }_{i}\,{c}_{i}|{\phi }_{i}\rangle \)133 as

$$E(\psi )=\frac{\int \,\langle \psi |\hat{H}|\psi \rangle }{\int \,\langle \psi |\psi \rangle }=\frac{{\sum }_{ij}\,{c}_{i}^{\ast }{c}_{j}\int \,\langle {\phi }_{i}|\hat{H}|{\phi }_{j}\rangle }{{\sum }_{ij}\,{c}_{i}^{\ast }{c}_{j}\,\int \,\langle {\phi }_{i}|{\phi }_{j}\rangle },$$
(175)

where \(\hat{H}\) is a Hamiltonian.

Then, let evaluate the corresponding energies of wavefunctions \(|{\psi }_{in}\rangle \), \(|\varphi \rangle \) and \(|{\Phi }^{\ast }\rangle \) via \({{\mathcal O}}_{V}\), as

$$S=E({\psi }_{in}),$$
(176)
$$X=E(\varphi ),$$
(177)

and

$$T=E({\Phi }^{\ast }).$$
(178)

Then, let Δ be the difference of the ratios of wavefunction energies, defined as

$$\Delta =R(S,T)-R(S,X)$$
(179)

where

$$R(S,T)=\frac{S}{T},$$
(180)

and

$$R(S,X)=\frac{S}{X}.$$
(181)

From the quantities of (176)–(178), let \({\rm{SNR}}(|{\Phi }^{\ast }\rangle )\) be the SNR of the output system \(|{\Phi }^{\ast }\rangle \), defined as

$$\begin{array}{rcl}{\rm{SNR}}(|{\Phi }^{\ast }\rangle ) & = & 10\,{\log }_{10}R(S,T)\\ & = & {\log }_{10}\Delta +\tfrac{1}{10}\,{\rm{SNR}}(|X\rangle ),\end{array}$$
(182)

where

$${\rm{SNR}}(|X\rangle )=10\,{\log }_{10}\,R(S,X),$$
(183)

while Δ is as given in (179).

Therefore, the SNR of the output system can be evolved from the difference of the ratios of the wavefunction energies as

$$\begin{array}{rcl}{\rm{SNR}}(|{\Phi }^{\ast }\rangle ) & = & 10{\log }_{10}R(S,T)\\ & = & 10({\log }_{10}\Delta +{\log }_{10}R(S,X))\\ & = & 10({\log }_{10}(R(S,T)-R(S,X))+{\log }_{10}R(S,X))\\ & = & 10({\log }_{10}\tfrac{R(S,T)}{R(S,X)}+{\log }_{10}R(S,X)).\end{array}$$
(184)

It also can be verified that Δ from (179) can be rewritten as

$$\Delta ={10}^{{\Delta }_{{\rm{SNR}}}/10},$$
(185)

where ΔSNR is an SNR difference, defined as

$${\Delta }_{{\rm{SNR}}}={\rm{SNR}}(|{\Phi }^{\ast }\rangle )-{\rm{SNR}}(|X\rangle ).$$
(186)

The high SNR values are reachable at moderate values of wavefunction energy ratio differences (179), therefore a high retrieval efficiency (high SNR values) can be produced by the local unitaries of the memory unit (see also Fig. 5).

Figure 5
figure 5

Verification of the retrieval efficiency of an HRE quantum memory unit via an \({{\mathcal O}}_{V}\) verification oracle. In the verification procedure, an unknown quantum system \(|\psi \rangle \) is stored in the \(QR\) quantum register that is evolved by an unknown operation \({U}_{QR}\) of the \(QR\) quantum register. The output of \(QR\) is an unknown quantum system \(|\varphi \rangle \) that is processed further by the \(U\) integrated unitary operations of the HRE quantum memory. The output system of the HRE quantum memory is \(|{\Phi }^{\ast }\rangle \) (170). The \({{\mathcal O}}_{V}\) oracle evaluates the SNR of the readout quantum system \(|{\Phi }^{\ast }\rangle \).

The proof is concluded here. ■

The verification of the retrieval efficiency of the output of an HRE quantum memory unit is depicted in Fig. 5.

Figure 6
figure 6

The output SNR values, \({\rm{SNR}}(|{\Phi }^{\ast }\rangle )=10\,{\log }_{10}\,R(S,T)\), of an HRE quantum memory in the function of \(\Delta =R(S,T)-R(S,X)\), where \(R(S,T)=\frac{S}{T}\), \(R(S,X)=\frac{S}{X}\), \(S=E({\psi }_{in})\), \(X=E(\varphi )\), and \(T=E({\Phi }^{\ast })\).

The output SNR values in the function of the Δ wave function energy ratio difference are depicted in Fig. 6.

Conclusions

Quantum memories are a cornerstone of the construction of quantum computers and a high-performance global-scale quantum Internet. Here, we defined the HRE quantum memory for near-term quantum devices. We defined the unitary operations of an HRE quantum memory and proved the learning procedure. We showed that the local unitaries of an HRE quantum memory integrates a group of quantum machine learning operations for the evaluation of the unknown quantum system, and a group of unitaries for the target system recovery. We determined the achievable output SNR values. The HRE quantum memory is a particularly convenient unit for gate-model quantum computers and the quantum Internet.

Ethics statement

This work did not involve any active collection of human data.