Estimating Entanglement Entropy via Variational Quantum Circuits with
  Classical Neural Networks
Entropy plays a crucial role in both physics and information science,
encompassing classical and quantum domains. In this work, we present the
Quantum Neural Entropy Estimator (QNEE), a novel approach that combines
classical neural network (NN) with variational quantum circuits to estimate the
von Neumann and Renyi entropies of a quantum state. QNEE provides accurate
estimates of entropy while also yielding the eigenvalues and eigenstates of the
input density matrix. Leveraging the capabilities of classical NN, QNEE can
classify different phases of quantum systems that accompany the changes of
entanglement entropy. Our numerical simulation demonstrates the effectiveness
of QNEE by applying it to the 1D XXZ Heisenberg model. In particular, QNEE
exhibits high sensitivity in estimating entanglement entropy near the phase
transition point. We expect that QNEE will serve as a valuable tool for quantum
entropy estimation and phase classification.
Sangyun Lee, Hyukjoon Kwon, Jae Sung Lee
# Estimating Entanglement Entropy via Variational Quantum Circuits with Classical Neural Networks

###### Abstract

Entropy plays a crucial role in both physics and information science, encompassing classical and quantum domains. In this work, we present the Quantum Neural Entropy Estimator (QNEE), a novel approach that combines classical neural network (NN) with variational quantum circuits to estimate the von Neumann and Renyi entropies of a quantum state. QNEE provides accurate estimates of entropy while also yielding the eigenvalues and eigenstates of the input density matrix. Leveraging the capabilities of classical NN, QNEE can classify different phases of quantum systems that accompany the changes of entanglement entropy. Our numerical simulation demonstrates the effectiveness of QNEE by applying it to the 1D XXZ Heisenberg model. In particular, QNEE exhibits high sensitivity in estimating entanglement entropy near the phase transition point. We expect that QNEE will serve as a valuable tool for quantum entropy estimation and phase classification.

## I Introduction

Entropy, originally introduced in classical thermodynamics and statistical physics, has evolved into a fundamental concept used to quantify information across various scientific disciplines, ranging from computer science to classical and quantum information theory. In particular, the Shannon entropy establishes a profound connection between thermodynamics and information theory, which is further extended to quantum systems through the introduction of the von Neumann entropy. This extension has paved the way for the exploration of intricate relationships between information and physical reality, as demonstrated by the discovery of information fluctuation theorems [1; 2] and Landauer's principle [3; 4; 5; 6] in both classical and quantum domains. Moreover, in the realm of quantum information theory, entropy has been employed to unveil non-classical characteristics inherent to quantum systems. Specifically, the entanglement entropy captures non-classical correlations between distinct parties of a quantum system and plays a pivotal role in distinguishing different phases of quantum materials [7; 8; 9; 10].

In spite of its recognized importance as a valuable theoretical tool, the practical determination of entropy values in physical systems presents substantial challenges. The difficulty arises primarily from the requirement of accurately evaluating the probability distribution, which becomes impractical in high-dimensional systems [11; 12]. In the realm of quantum systems, estimating entropy becomes even more formidable due to the dependence on measurement choices. Specifically, the von Neumann entropy of a quantum state is determined as the minimum entropy among all potential measurement bases, which is achieved by taking the state's eigenbasis. This additional layer of complexity introduces difficulties in obtaining the eigendecomposition of a quantum density matrix, with computational costs exponentially increasing with the size of the quantum system. While several indirect methods have been proposed [13; 14; 15; 16] and experimentally demonstrated [17; 18; 19; 20] to estimate or bound entanglement of quantum states, limitations in accuracy persist. Recently, variational quantum algorithms (VQA) [21], such as the variational quantum state diagonalization (VQSD) [22] and the variational quantum state eigensolver (VQSE) [23] have emerged as potential approaches for estimating quantum state eigenvalues. However, it still remains uncertain whether these methods can efficiently estimate quantum entropy.

The challenging nature associated with quantum entropy estimation prompts the pursuit of a novel and efficient algorithm that can be effectively integrated with a variational quantum circuit, leading to improved performance. To this end, one potential candidate is classical neural networks (NNs), which have demonstrated exceptional capabilities in optimizing objective functions. Notably, this machine-learning technique has recently found application in estimating entropy production for classical thermodynamic systems, named as neural estimator of entropy production (NEEP) algorithm [24], by optimizing the objective function defined by the Donsker-Varadhan inequality [25]. This utilization of NNs in classical thermodynamics provides inspiration for their potential application in the context of quantum entropy estimation. By leveraging the strengths of classical NNs in optimizing von Neumann entropy, it would be possible to enhance the accuracy and efficiency of quantum entropy estimation, paving the way for further advancements in the field.
In this work, we present a quantum-classical hybrid algorithm, Quantum Neural Entropy Estimator (QNEE), for estimating the von Neumann and Renyi entropies of a quantum state. Our method capitalizes on the efficient entropy estimation capabilities of NNs from the outcomes of parametrized quantum circuits while utilizing variational parameters to align the quantum circuit in the correct measurement basis. We highlight that our methodology can be readily extended for any quantum circuit ansatz as the classical NN is constructed independent of the variational quantum circuit's structure.

Classical NN using the cost function based on the Donsker-Varadhan inequality has some useful properties, such as robustness against small-size data and rare events, which have already been verified in classical stochastic processes [12; 24]. Another advantage of this type of cost function is that it can provide a reliable upper bound of the entropy even if the measurement basis is not fully aligned correctly. This opens a possible application of the proposed approach as a phase classifier of quantum systems even without the accurate estimation of entanglement value.

To showcase the effectiveness of our proposed protocol, we apply it to estimate the ground state entanglement entropy of the Heisenberg XXZ model and to distinguish its phase relying on the external magnetic field. Our numerical results demonstrate that QNEE exhibits good performance in estimating entanglement entropy as well as the subsystem's eigenvalues. We also provide a detailed comparison between QNEE and the existing scheme of VQSE.

This paper is organized as follows. In Sec. II, we construct QNEE by combining a classical NN with a variational quantum circuit. In Sec. III, we apply QNEE to estimate the ground state entanglement entropy of the 1D XXZ Heisenberg model and its performance as a phase classifier. We conclude the paper in Sec. IV.

**Concurrent work:** In the course of finalizing this manuscript, we became aware of related research that has emerged in arXiv [26; 27]. It is essential to highlight that our work draws inspiration from the authors' earlier contributions [12; 24; 28] as its quantum extension. Therefore, our cost function for NN is distinct from that employed in the concurrent works. Furthermore, we must emphasize that the primary focus of our work lies in the application of phase classification to a physical Hamiltonian. In contrast, the concurrent works under examination delve into the investigation of random quantum states as illustrative examples. These distinctions in scope and approach further underscore the unique contributions and significance of our study.

## II Quantum entropy estimator

### Cost function of QNEE

In this section, we establish QNEE's cost function to estimate the quantum entropy of a given density matrix \(\hat{\rho}\), describing the state of a system. The cost function is based on the following Gibbs variational principle:

\[S(\hat{\rho})\equiv-\mathrm{tr}\{\hat{\rho}\ln\hat{\rho}\}\leq- \langle\hat{O}\rangle_{\hat{\rho}}+\ln\mathrm{tr}\{e^{\hat{O}}\} \tag{1}\]

for any Hermitian operator \(\hat{O}\), whose expectation value is defined as \(\langle\hat{O}\rangle_{\hat{\rho}}=\mathrm{tr}\{\hat{O}\hat{\rho}\}\). We present a simple derivation of Eq. (1) in Appendix A. Equation (1) was utilized as a quantum entropy estimator in Ref. [27]. However, as explained in Sec. II.3 and Appendix A, a cost function, including the logarithmic function, can be biased and malfunction for large data if additional manipulation is not employed. To avoid this complication, we further linearize Eq. (1) using the inequality \(\ln\mathrm{tr}\{e^{\hat{O}}\}\leq\mathrm{tr}\{e^{\hat{O}}\}-1\) as

\[S(\hat{\rho})\leq-\langle\hat{O}\rangle_{\hat{\rho}}+\mathrm{tr} \{e^{\hat{O}}\}-1. \tag{2}\]

Equation (2) is usually a more appropriate cost function for the stochastic (mini-batch) gradient method [29], which is typically used for large-system NN optimization [30]. The both bounds (1) and (2) imply that we can indirectly estimate the von Neumann entropy using the operator \(\hat{O}\)'s expectation value and its eigenvalues. However, the operator saturating the both bounds in Eqs. (1) and (2) is given by \(\hat{O}^{*}=\ln\hat{\rho}\). In other words, the exact estimation using these bounds requires the information of \(\hat{\rho}\). Nevertheless, an optimization method enables us to estimate the von Neumann entropy without accessing full information of the density matrix \(\hat{\rho}\). For classical stochastic processes, various optimization methods have been employed to estimate the entropy production without detailed information [12; 31; 32; 33; 34; 35; 36; 37; 38; 39; 40; 41; 42; 43; 44; 45; 46; 47; 48; 49; 50; 51; 52; 53; 54; 55; 56; 57; 58; 59; 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; 85; 86; 87; 88; 89; 90; 82; 84; 86; 87; 88; 89; 91; 92; 93; 94; 95; 96; 97; 98; 99; 100; 101; 102; 103; 104; 105; 106; 107; 108; 109; 110; 111; 112; 113; 114; 115; 116; 117; 118; 119; 120; 121; 122; 123; 124; 125; 126; 127; 128; 129; 130; 131; 132; 133; 134; 135; 136; 137; 138; 139; 140; 141; 142; 143; 144; 145; 146; 147; 148; 149; 150; 151; 152; 153; 154; 155; 156; 157; 158; 159; 160; 161; 162; 163; 164; 165; 166; 167; 168; 169; 170; 171; 172; 173; 174; 175; 176; 1777; 178; 179; 180; 181; 182; 183; 184; 185; 186; 187; 188; 189; 190; 188; 189; 189; 191; 200; 201; 202; 203; 204; 205; 206; 207; 208; 209; 210; 211; 213; 214; 215; 216; 217; 218; 219; 22; 231; 232; 233; 234; 235; 236; 237; 238; 239; 240; 251; 252; 253; 254; 255; 256; 257; 258; 26; 279; 28; 281; 282; 283; 284; 285; 286; 287; 288; 289; 290; 291; 288; 289; 292; 293; 294; 295; 296; 297; 300; 301; 302; 303; 304; 305; 306; 307; 308; 309; 310; 311; 329; 333; 341; 342; 358; 359; 360; 371; 372; 373; 381; 382; 383; 398; 399; 40; 41; 42; 43; 44; 45; 46; 47; 48; 49; 50; 51; 52; 53; 54; 55; 56; 57; 58; 59; 60; 61; 62; 63; 64; 65; 66; 67; 68; 68; 69; 70; 72; 74; 75; 76; 777; 78; 79; 80; 81; 83; 84; 85; 86; 87; 88; 89; 92; 88; 89; 93; 941; 88; 87; 88; 89; 95; 96; 97; 98; 99; 99; 101; 11; 12; 13; 14; 15; 16; 17; 18; 199; 192; 193; 194; 195; 196; 197; 198; 199; 202; 203; 204; 205; 206; 207; 208; 209; 211; 224; 215; 216; 217; 218; 219; 232; 241; 219; 25; 26; 27; 28; 299; 301; 312; 32; 333; 342; 35; 361; 373; 38; 39; 40; 410; 411; 42; 43; 45; 46; 47; 48; 49; 51; 50; 53; 54; 56; 57; 58; 59; 61; 70; 71; 72; 73; 74; 75; 76; 78; 89; 931; 80; 81; 82; 83; 84; 85; 86; 87; 88; 89; 94; 95; 96; 97; 101; 11; 12; 13; 14; 15; 16; 17; 18; 198; 199; 203; 21; 222; 23; 24; 25; 26; 27; 28; 29; 31; 32; 334; 35; 37; 39; 50; 51; 538; 59; 60; 61; 62; 63; 640; 65; 67; 68; 69; 70; 73; 71; 72; 74; 76; 78; 79; 81; 80; 82; 84; 85; 86; 87; 88; 89; 99; 90; 911; 92; 933; 94; 95; 96; 97; 103; 98; 104; 105; 106; 107; 108; 107; 108; 119; 111; 139; 120; 133; 134; 14; 15; 16; 17; 18; 199; 193; 18; 194; 195; 196; 197; 198; 204; 21; 22; 24; 26; 27; 28; 299; 303; 31; 32; 34; 36; 37; 38; 39; 40; 41; 43; 44; 45; 46; 47; 48; 49; 51; 52; 53; 54; 56;
where \(\hat{\rho}_{\hat{V}}=\hat{V}\hat{\rho}\hat{V}^{\dagger}\). The optimization over \(\hat{O}\) then becomes finding the optimal diagonal matrix \(\hat{D}\) and the unitary operator \(\hat{V}\). By taking \(\hat{D}(\Theta_{N})=\sum_{i=1}^{2^{n}}h(s_{i};\Theta_{N})|s_{i}\rangle\langle s_{ i}|\) in the computational basis for \(n\) qubits \(\{|s_{i}\rangle\}=\{|00\cdots 0\rangle,\cdots,|11\cdots 1\rangle\}\), we encode \(h(s_{i};\Theta_{N})\) into the output function of a classical NN with a weight set \(\Theta_{N}\). We also encode the unitary operator \(\hat{V}(\Theta_{Q})\) into a variational quantum circuit parameterized by a vector \(\Theta_{Q}\). The QNEE's cost function \(C(\Theta_{N},\Theta_{Q})\) is then defined as

\[S(\hat{\rho})\leq-\sum_{i=1}^{2^{n}}h(s_{i};\Theta_{N})P_{\hat{V}}(s_{i}; \Theta_{Q})+\sum_{i=1}^{2^{n}}e^{h(s_{i};\Theta_{N})}-1\equiv C(\Theta_{N}, \Theta_{Q}), \tag{4}\]

where \(P_{\hat{V}}(s_{i};\Theta_{Q})=\langle s_{i}|\hat{V}(\Theta_{Q})\hat{\rho}\hat {V}^{\dagger}(\Theta_{Q})|s_{i}\rangle\) is the probability of the quantum circuit parameterized by \(\Theta_{Q}\) to output the string \(s_{i}\). Here we define \(\hat{V}_{\rm D}\) as the unitary matrix rendering \(\hat{\rho}\) diagonalized as \(\rho_{\hat{V}_{\rm D}}=\hat{V}_{\rm D}\hat{\rho}\hat{V}_{\rm D}^{\dagger}= \sum_{i}\lambda_{i}|s_{i}\rangle\langle s_{i}|\), where \(\lambda_{i}=\langle s_{i}|\hat{\rho}_{\hat{V}_{\rm D}}|s_{i}\rangle\) is the probability for measuring \(s_{i}\) from \(\hat{\rho}_{\hat{V}_{\rm D}}\). When \(\hat{V}(\Theta_{Q})=\hat{V}_{\rm D}\) and \(h(s_{i};\Theta_{N})=\ln\lambda_{i}\), the inequality in Eq. (4) is saturated, in which case \(P_{\hat{V}}(s_{i};\Theta_{Q})=\lambda_{i}\).

Equation (4) is reminiscent of the cost function of NEEP: a neural estimator for entropy production [24]. Both the cost functions of QNEE and NEEP can be derived from the \(f\)-divergence [33; 34]. The unsupervised learning with a variational cost function has been established as an effective method to estimate the divergence of data rather than other methods, especially when the number of samples is inadequate or the data include rare events [12; 24]. Because the entropy can be considered as the divergence that one distribution is given as a maximally mixed state, we expect that QNEE would be effective when the given density matrix has a few dominant eigenvalues while the remainings are nearly zero, or the number of given shots is insufficient, which is the usual situation for many-qubit systems. The von Neumann entropy is small in general under these conditions.

The \(f\)-divergence yields a diverse set of information measures and their variational functional forms. Consequently, this enables the formulation of cost functions for estimating various physical quantities. One prominent example of such cost functions is the Renyi entropy, defined as

\[S_{\alpha}(\hat{\rho})=\frac{1}{1-\alpha}\ln\mathrm{tr}\{\hat{\rho}^{\alpha}\} \tag{5}\]

for \(\alpha>0\) and \(\alpha\neq 1\). The corresponding cost function of Eq. (5) is

\[\frac{e^{(1-\alpha)S_{\alpha}(\hat{\rho})}-1}{\alpha(1-\alpha)}\leq \sum_{i=1}^{2^{n}}P_{\hat{V}}(s_{i};\Theta_{Q})\frac{e^{(\alpha-1 )h(s_{i};\Theta_{N})}-1}{(1-\alpha)}+\sum_{i=1}^{2^{n}}\frac{e^{\alpha h(s_{i} ;\Theta_{N})}-1}{\alpha}\equiv C_{\alpha}(\Theta_{N},\Theta_{Q}). \tag{6}\]

The detailed derivation of Eq. (6) from the \(f\)-divergence is presented in Appendix B. It is also straightforward to show that Eq. (6) is saturated under the same conditions as Eq. (4), that is, \(h(s_{i};\Theta_{N})\rightarrow\ln\lambda_{i}\) and \(P_{\hat{V}}(s_{i};\Theta_{Q})\rightarrow\lambda_{i}\). With the optimal parameters saturating the equality, the cost function can be easily converted into the Renyi entropy. Note that the cost function for Renyi entropy used in Ref. [27] consists of the two logarithmic terms, which has the same complication as the Gibbs variational principle (1), not straightforwardly applicable to the stochastic gradient method to handle large data without additional data manipulation. These two logarithmic terms of the cost function in Ref. [27] cannot be simply linearized as done in Eq. (2) either since it is not trivial to determine which one is larger between the original one and the linearized one. Therefore, it is not easy to find the variational form appropriate for the stochastic gradient method from the cost function in Ref. [27]. In contrast, our cost function has no such problem, thus, it would be more appropriate for manipulating larger quantum system via NN.

Equation (6) is reduced to the bound for the von Neumann entropy (4) in the limit of \(\alpha\to 1\). Moreover, Eq. (6), originally derived for estimating the Renyi entropy, can also be applied to estimate the von Neumman entropy even for general \(\alpha(\neq 1)\) cases. This is due to the fact that the NN output with optimal parameters \(\Theta_{N}^{*}\) for the Renyi entropy is the same as that of the von Neumann entropy (4), i.e., \(h(s_{i};\Theta_{N}^{*})=\ln\lambda_{i}\) for any \(\alpha\). This freedom for choosing \(\alpha\) arbitrarily opens another possibility to render the NN performance enhanced. A similar technique was used to estimate stochastic entropy production more accurately by tuning \(\alpha\)[33].

### Overview of the QNEE algorithm

Here we present an outline of QNEE's estimation procedure. Figure 1 shows the schematic diagram of the comprehensive QNEE algorithm. As depicted in Fig. 1(a), we prepare an \(n\)-qubit quantum system described by the density matrix \(\hat{\rho}\), for which the von Neumann entropy needs to be estimated. Additionally, we randomly select the initial values of \(\Theta_{N}\) and \(\Theta_{Q}\), which serve as the parameters for the NN and the quantum circuit, respectively. This setup
allows us to commence the repetitive procedure for updating \(\Theta_{N}\) and \(\Theta_{Q}\) until the cost function in Eq. (4) is optimized. The one unit, as illustrated in Figs. 1(b) and (c), of this repetitive procedure can be divided into two processes in terms of the updating parameters: (i) the updating process of \(\Theta_{N}\) via the optimization of the cost function using the NN and (ii) the updating process of \(\Theta_{Q}\) obtained through the calculation of the gradient of the cost function.

The brief overview of this one-unit process is as follows. First, training and test data sets, consisting of \(N_{s}\) strings each, are produced from the quantum circuit with a given \(\Theta_{Q}\). With the training set, we perform unsupervised learning to update \(\Theta_{N}\) using the classical NN, along with minimizing the cost function. In the course of the \(\Theta_{N}\) updating process, we evaluate the cost function by inputting the test data set into the NN several times and finding the lowest one among them. This lowest value of the cost function denoted by \(C^{\rm NN}\) and the corresponding parameters \(\Theta_{N}\) are recorded. Details on this NN optimization process and the structure of the NN are presented in Sec. II.3. For updating \(\Theta_{Q}\), we need to evaluate the gradient of the cost function with respect to \(\Theta_{Q}\), i.e., \(\nabla_{\Theta_{Q}}C^{\rm NN}\). To do this, it is necessary to obtain \(C^{\rm NN}\) for small perturbation \(\delta\) from \(\Theta_{Q}\) in various directions through the NN optimization method mentioned above. Employing this gradient descent method, \(\Theta_{Q}\) is updated slightly toward the direction that minimizes \(C^{\rm NN}\). Details on the structure of the quantum circuit and the \(\Theta_{Q}\) updating process are explained in Sec. II.4.

Upon conducting an adequate number of repetitions of this unit process to update \(\Theta_{N}\) and \(\Theta_{Q}\), we anticipate the optimization of parameters and cost functions. Utilizing the obtained optimal parameters \(\Theta_{N}^{*}\) and \(\Theta_{Q}^{*}\), QNEE produces three essential outputs: the von Neumann entropy (\(S(\hat{\rho})\)), the eigenvalues of the quantum state (\(\lambda_{i}\)), and their corresponding eigenvectors (\(|\lambda_{i}\rangle\)) as follows:

\[S(\hat{\rho}) =C(\Theta_{N}^{*},\Theta_{Q}^{*}), \tag{7}\] \[\lambda_{i} =\exp{[h(s_{i};\Theta_{Q}^{*})]},\] (8) \[|\lambda_{i}\rangle =\hat{V}^{\dagger}(\Theta_{Q}^{*})|s_{i}\rangle=\hat{V}_{D}^{ \dagger}|s_{i}\rangle. \tag{9}\]

We note that the same process can be employed to estimate the Renyi entropy with the cost function in Eq. (6).

Figure 1: Architecture of QNEE. (a) QNEE requires three inputs: density matrix \(\hat{\rho}\), initial parameter \(\Theta_{Q}\) for a quantum circuit, and \(\Theta_{N}\) for a NN. (b) From the given \(\Theta_{Q}\) and density matrix \(\hat{\rho}\), two string data sets are generated: One for training and the other one for testing. (c) We update \(\Theta_{N}\) by minimizing the cost function evaluated with the train data, \(\min_{\Theta_{N}}C^{\rm train}\). During the learning process, the cost function evaluated with the test data \(C^{\rm test}\) is recorded. Then, the NN outputs the lowest value of \(C^{\rm test}\) from the recorded values, denoted as \(C^{\rm NN}\). The quantum circuit’s parameter \(\Theta_{Q}\) is updated by minimizing the NN’s output \(\min_{\Theta_{Q}}C^{\rm NN}\). On successful training, QNEE finds the optimal parameters \(\Theta_{N}^{*}\) and \(\Theta_{Q}^{*}\), and QNEE generates three outputs. (d) The cost function itself becomes von Neumann entropy with optimal parameters. (e) The trained NN can generate the eigenvalues of \(\hat{\rho}\). (f) With a conjugated quantum circuit \(\hat{V}^{\dagger}(\Theta_{Q}^{*})\), eigenvectors can be generated.
### Structure of the classical NN and optimization of \(\Theta_{N}\)

In this section, we elaborate the structure and training process of the classical NN. We design the NN as a combination of an embedding layer and fully connected three hidden layers parameterized by \(\Theta_{N}\). The input string is converted into a vector via the embedding layer and is subsequently processed through the hidden layers with 256 units each along with ReLU activation functions [35]. Therefore, the NN can transform the string data generated via the quantum circuit \(\hat{V}(\Theta_{Q})\) into the scalar value \(h(s;\Theta_{N})\). To optimize \(\Theta_{N}\) of the classical NN, we first generate training and test sets consisting of \(N_{s}\) strings each from the quantum circuit parameterized by \(\Theta_{Q}\). Using the training set, we evaluate the cost function in Eq. (4) as

\[C^{\text{train}}(\Theta_{N},\Theta_{Q})=-\frac{1}{N_{s}}\sum_{s\in\text{train set}}h(s;\Theta_{N})+\sum_{i=1}^{2^{n}}e^{h(s_{i};\Theta_{N})}-1. \tag{10}\]

Over a predetermined number of iterations \(N_{\text{iter}}\), \(\Theta_{N}\) is updated using ADAM optimizer [36] for minimizing the cost function.

During these training iterations, we evaluate the cost function several times using the test set with the updated \(\Theta_{N}\) as

\[C^{\text{test}}(\Theta_{N},\Theta_{Q})=-\frac{1}{N_{s}}\sum_{s\in\text{test set}}h(s;\Theta_{N})+\sum_{i=1}^{2^{n}}e^{h(s_{i};\Theta_{N})}-1. \tag{11}\]

Among these several evaluated \(C^{\text{test}}\), the lowest one is recorded as \(C^{\text{NN}}\) and the corresponding \(\Theta_{N}\) is selected for the NN parameters at a given \(\Theta_{Q}\). Note that the learning curves shown in Figs. 4(d-f) are the plot of the recorded \(C^{\text{NN}}\).

When the NN is ideally trained, the optimal parameter for the NN is given as \(\tilde{\Theta}_{N}(\Theta_{Q})=arg\,\min_{\Theta_{N}}C(\Theta_{N},\Theta_{Q})\). Then, the output of the NN and the cost function are written as

\[h(s_{i};\tilde{\Theta}_{N})= \ln P_{\hat{V}}(s_{i};\Theta_{Q}), \tag{12}\] \[C(\tilde{\Theta}_{N},\Theta_{Q})= -\sum_{i=1}^{2^{n}}P_{\hat{V}}(s_{i};\Theta_{Q})\ln\{P_{\hat{V}}( s_{i};\Theta_{Q})\}, \tag{13}\]

respectively. Equation (13) is the Shannon entropy \(S[P_{\hat{V}}(s_{i};\Theta_{Q})]\) in terms of the diagonal elements of \(\hat{\rho}_{\hat{V}}\). Therefore, we can write down the following variational principle with the ideal-string distribution \(P_{\hat{V}}(s_{i};\Theta_{Q})\) as

\[S[P_{\hat{V}}(s_{i};\Theta_{Q})]\leq-\sum_{i=1}^{2^{n}}h(s_{i};\Theta_{N})P_{ \hat{V}}(s_{i};\Theta_{Q})+\sum_{i=1}^{2^{n}}e^{h(s_{i};\Theta_{N})}-1. \tag{14}\]

This bound (14) can also be derived from the Donsker-Varadhan inequality \(D_{\text{KL}}[P\|Q]\geq\langle h(s)\rangle_{P}-\ln\langle e^{h(s)}\rangle_{Q}\) which is the variational form of the Kullback-Leibler divergence \(D_{\text{KL}}[P\|Q]=\sum_{i}P(s_{i})\ln\left[P(s_{i})/Q(s_{i})\right]\) for two distributions \(P(s_{i})\) and \(Q(s_{i})\)[25]. One can easily derive Eq. (14) from the Donsker-Varadahan inequality by choosing the probability \(Q\) as a uniform distribution and applying the inequality \(\ln x\leq x-1\) to the logarithm function. Note that Eq. (12) reveals the role of the last two terms in \(C^{\text{train}}\); these terms \(\sum_{i=1}^{2^{n}}e^{h(s_{i};\Theta_{N})}-1\) ensure that \(e^{h}\) represents a normalized probability.

Finally, it is worth noting that although Eq. (10) is currently formulated as a full-batch summation, it can be adapted to utilize mini-batch evaluation for calculating the cost function. This implies that the stochastic gradient descent method [29] can be directly applied to the QNEE's cost function without any modifications, presenting an advantage in using the linearization of the logarithmic function, \(\ln x\leq x-1\). The stochastic gradient descent method is efficient in handling large-scale data and facilitating optimizer escape from local minima. Given the exponential increase in data size with the number of qubits, the practical necessity of stochastic gradient descent for managing such vast amounts of data becomes evident. Moreover, we can also extend the application of stochastic gradient descent to the cost function for Renyi entropy estimation, as given in Eq. (6).

### Structure of the quantum circuit and optimization of \(\Theta_{Q}\)

With the cost function calculated by the NN, we update the quantum circuit. Here we explain our choice of ansatz and the optimization process of the quantum circuit. Optimizing \(\Theta_{Q}\) can be viewed as a VQA, the accuracy of which
heavily depends on the structure of the target quantum state and quantum circuit ansatz. In this work, we choose the layered hardware-efficient ansatz [37], which is widely adopted when intricate information of \(\hat{\rho}\) is inaccessible. The ansatz consists of two-qubit gates of one kind and parameterized one-qubit gates, as described in Fig. 2. \(\hat{V}(\Theta_{Q})\) is constructed from \(N_{l}\) layered two-qubit gates \(\bar{V}(\Theta_{Q}^{(i)})\) that are drawn as yellow paths. \(\bar{V}(\Theta_{Q}^{(i)})\) contains four \(\hat{R}_{Y}\) gates and one controlled-Z gate, which is the same choice as the quantum circuit in Ref. [23]. Each quantum circuit layer can be an identity gate with certain parameters to guarantee that increasing \(N_{l}\) improves the accuracy of QNEE [23].

In this study, a gradient descent method is chosen to optimize the cost function via the quantum circuit. To this end, it is necessary to estimate the cost function at differently perturbed \(\Theta_{Q}\)s for obtaining its gradient. For example, if the ansatz is parameterized by one single scalar, we estimate \(C^{\rm NN}(\Theta_{Q})\) and \(C^{\rm NN}(\Theta_{Q}+\delta)\) via the NN, and calculate the gradient \(\nabla_{\Theta_{Q}}C(\Theta_{Q})=[C^{\rm NN}(\Theta_{Q}+\delta)-C^{\rm NN}( \Theta_{Q})]/\delta\), where \(\delta\) is a small value. Then, the parameter is updated as \(\Theta_{Q}-\eta_{Q}\nabla_{\Theta_{Q}}C(\Theta_{Q})\) with a learning rate \(\eta_{Q}\). Therefore, two training processes of the NN are required to calculate the gradient for a scalar \(\Theta_{Q}\) case. When the dimension of \(\Theta_{Q}\) increases, the number of necessary directions for perturbation to obtain the gradient also increases. Using the calculated gradient and a learning rate that we choose, \(\Theta_{Q}\) is updated. After an adequate number of iterations of this \(\Theta_{Q}\) updating process, we take the lowest value of the cost function as the estimated von Neumann entropy among all \(C^{\rm NN}\) values evaluated during the iteration procedure.

When the quantum circuit is ideally trained, the optimal parameter is found as \(\Theta_{Q}^{*}=arg\min_{\Theta_{Q}}C(\bar{\Theta}_{N},\Theta_{Q})\) that diagonalizes the input density matrix with the eigenvalue \(\lambda_{i}=P_{\hat{V}_{D}}(s_{i};\Theta_{Q}^{*})\). Therefore, the cost function is saturated to the von Neumann entropy, and the output of the NN with \(\Theta_{N}^{*}=\arg\min_{\Theta_{N}}C(\Theta_{N},\Theta_{Q}^{*})\) is the logarithm of the eigenvalue as shown in Eqs. (7) and (8), respectively. Additionally, the eigenvectors of \(\hat{\rho}\) can be constructed by the conjugated quantum circuit \(\hat{V}^{\dagger}(\Theta_{Q}^{*})\) as Eq. (9).

Usually, many qubit states or highly entangled states require an ansatz with deep depth. However, the increased dimension of \(\Theta_{Q}\) requires more data to compute the gradient of the cost function. Moreover, a large number of layers can induce a flat landscape of a cost function, which leads to difficulty in training. This is named the Barren plateau [38]. Although various methods are suggested to mitigate the Barren Plateau, being free of this effect seems uneasy. Even with a shallow quantum circuit, QNEE provides useful applications, which will be discussed in Section III in more detail.

### Other methods for entropy estimation

Various VQAs are applicable to estimate quantum entropy [22; 23; 39]. These algorithms can be classified into two categories. One category is the VQAs based on purity minimization. VQSD is one of these algorithms. VQSD employs the SWAP (swapping two qubits) test to compute the purity \(\mathrm{tr}\{\hat{\rho}^{2}\}\). As the result of minimization, the input density matrix is diagonalized, from which the von Neumann entropy can be estimated. The SWAP test requires \(\hat{\rho}\otimes\hat{\rho}\) as its input, which means that the width of the quantum circuit must be twice the number of qubits of \(\hat{\rho}\). In addition, another quantum circuit is necessary to compute the cost function.

The other category exploits majorization [23]; we say that \(\vec{x}\) majorizes \(\vec{y}\), when \(\sum_{i}^{k}x_{i}^{\downarrow}\geq\sum_{i}^{k}y_{i}^{\downarrow}\) for arbitrary \(k\) and \(\sum_{i}^{m}x_{i}=\sum_{i}^{m}y_{i}\)[40]. Here \(m\) is the dimension of the vector, and \(\vec{x}^{\downarrow}\) is a decreasing ordered vector, i.e., \(x_{i}^{\downarrow}\geq x_{i+1}^{\downarrow}\). Because the eigenvalues of \(\hat{\rho}\) majorize the diagonal elements of \(\hat{\rho}\) and the inner product with an ordered energy-level vector is a Schur convex function, the expectation value of an artificially ordered Hamiltonian can be a cost function

Figure 2: Schematic diagrams of (a) layered hardware efficient ansatzs \(V(\Theta_{Q})\) and (b) its block gate \(\bar{V}\). \(V(\Theta_{Q})\) consists of \(N_{l}\) layers of \(\bar{V}\) block gates. In this research, we choose \(\bar{V}\) as the combination of four \(R_{Y}\) gates and one controlled-Z gate. For \(n\) qubits, the dimension of \(\Theta_{Q}\) or the total number of parameters is given as \(2nN_{l}\).
for diagonalizing \(\hat{\rho}\)'s subspace (refer [23] for detailed information). After optimization, the \(\ell\)-largest eigenvalues are estimated using the quantum circuit with the optimal parameter \(\Theta_{Q}^{*}\). From the string data, we can estimate the eigenvalue \(\lambda_{i}\sim N_{i}/N_{s}\) where \(N_{s}\) denotes the number of total shots, and \(N_{i}\) denotes the number of \(s_{i}\) data. One benefit of using the cost function is that the requirement for the width of the quantum circuit only needs to match the number of qubits of \(\hat{\rho}\), which is the same requirement as the QNEE. Note that the VQSE focuses on finding a few largest eigenvalues, which makes the VQSE to be well saturated by expanding the allowed solution space of \(V(\Theta_{Q})\). We note that QNEE's intermediate cost function (see Eq. (13)) also exploits majorization (see Appendix B). Due to the same structure of the quantum circuit ansatz \(V(\Theta_{Q})\) used in QNEE and VQSE, we can fairly compare QNEE with VQSE under the same conditions. We provide this comparison in Sec. III.

## III Application of QNEE to a physical model

### Applications of QNEE

QNEE is a combination of the classical NN and the quantum circuit to estimate the von Neumann entropy. Hence, QNEE shares some advantages of using classical NNs for entropy estimation. For instance, NN such as NEEP has been reported to show high precision in estimating large-valued \(f\)-divergences, which corresponds to low entropy (see Appendix B), with small-size data or data with rare events [24; 12]. In this sense, one can expect that QNEE provides high precision for the small-valued von Neumann entropy. While estimating the quantum entropy additionally requires the diagonalization of the density matrix, we note that this does not necessitate a deep depth of the circuit ansatz for low-entropy quantum states. It is also worth noting that QNEE provides the upper bound Eq. (4) of the von Neumann entropy even when \(\Theta_{Q}\) is not fully optimized throughout the training process. This guarantees that QNEE efficiently characterizes quantum states with low entropy.

Based on these properties originating from the classical NN, QNEE can be utilized as a phase classifier for quantum states even with a shallow-depth ansatz and a small number of iterations. Particularly, it has been actively studied in many-body physics that phase transitions in various many-body quantum systems accompany the drastic change of the entanglement entropy, given by the von Neumann entropy of the subsystem [7; 8; 9; 10; 41; 42; 43]. As an illustrative example, we test our protocol using the XXZ Heisenberg model in a 1D chain (XXZ chain). This model shows the Pokrovsky-Talapov transition that accompanies a sudden change of entanglement entropy [41]. The entanglement entropy abruptly decreases after passing the critical point.

In the following subsections, we present numerical results of classifying the phase of the XXZ chain using QNEE, and compare its performances with the existing methods of VQSE [23]. We employ PyTorch [44], Qiskit [45], and QuTiP [44] to implement QNEE.

### Classifying the phase of XXZ chain

In order to verify the efficiency of QNEE as a phase classifier, we consider the ground state of the XXZ chain [41; 46; 47] Hamiltonian,

\[H_{XXZ}=\sum_{l=0}^{L-1}(\hat{\sigma}_{l}^{x}\hat{\sigma}_{l+1}^{x}+\hat{ \sigma}_{l}^{y}\hat{\sigma}_{l+1}^{y}+\Delta\hat{\sigma}_{l}^{z}\hat{\sigma}_ {l+1}^{z}-\lambda\hat{\sigma}_{l}^{z}), \tag{15}\]

where \(\hat{\sigma}_{l}^{k}\) denotes \(l\)-th spin operator in \(k\in\{x,y,z\}\) direction and \(\Delta\), \(\lambda\) are real numbers. \(L\) is the total number of qubits in the chain. The last term \(\lambda\hat{\sigma}_{l}^{z}\) can be interpreted as the effects of external magnetic field or chemical potential [46]. Here, we impose the periodic boundary condition.

The ground state \(|\Psi_{G}\rangle\) of the XXZ chain shows phase transition at \(\lambda_{c}=2(1-\Delta)\), which is named Pokrovsky-Talapov transiton [41]. For \(0\leq\lambda<\lambda_{c}\), the system is in the critical or massless phase. In the critical regime, the entanglement entropy shows log scaling behavior \(S\propto\ln n+c\), where \(n\) is the number of qubits of the reduced matrix, and \(c\) is a constant [9]. In the case of the finite-size XXZ chain, the scaling relation is modified as \(S\propto\ln\{\frac{L}{n}\sin\frac{n}{L}\}+c^{\prime}\). For \(\lambda>\lambda_{c}\), the system is in the non-critical or ferromagnetic phase, where the entanglement entropy vanishes. This implies that a sudden decrease of the entanglement entropy can be observed as \(\lambda\) increases.

We estimate the entanglement entropy of the XXZ chain with 8 qubits by estimating the von Neumann entropy of the ground state's reduced density matrix \(\rho=\text{tr}_{B}\{|\Psi_{G}\rangle\langle\Psi_{G}|\}\), where \(\text{tr}_{B}\) denotes partial tracing out the sub system \(B\). We set \(\Delta=0.05\) and vary magnetic field parameter \(\lambda\) from 0 to 3. Near the critical point \(\lambda\sim 2\), the entanglement entropy undergoes the phase transition. The estimation results of QNEE for the three and four-qubit reduced density
matrices are presented in Fig. 3. We emulate the variational quantum circuit using a classical computer by collecting the outcome strings from the exact probability distribution of the quantum circuit's outcome. The number of layer \(N_{l}\) is taken as 8 for the three-qubit case and 10 for the four-qubit case. The learning rate of gradient-descent method for the quantum circuit is set as \(\eta_{Q}=0.01\) for both cases. For the NN training process, the learning rate and weight decay of ADAM optimizer are set as 0.00001 and 0.00005, respectively. We train the NN over \(N_{\text{iter}}=10^{4}\) iterations at an initial \(\Theta_{Q}\) for sufficient optimization. In the subsequent NN training processes, the NN is trained over \(N_{\text{iter}}=100\) iterations at each \(\Theta_{Q}\). The entanglement entropy is estimated from five different initial \(\Theta_{Q}\), and the minimum value among them is taken. In the insets of Fig. 3, we plot the mean and standard deviation of the five different trials. For each \(\Theta_{Q}\), 30,000 shots are used to find the optimal parameter \(\bar{\Theta}_{N}\).

The numerical results demonstrate that QNEE well estimates the entanglement entropy of the XXZ chain ground state for a wide range of the parameter regime. In particular, QNEE provides an accurate estimation near \(\lambda_{c}\) and the non-critical regime. This is expected as these are low-entropy states that a classical NN can estimate well. In particular, the standard deviation of QNEE is smaller than the size of markers when \(\lambda\geq\lambda_{c}\), which means that a single trial is enough to classify those phases with QNEE. We also discuss the learning rate of QNEE.

Another output of QNEE is the eigenvalues of the reduced density matrix, the so-called entanglement spectrum. In Fig. 4 (a - c), we plot the estimated eigenvalues for three different values of magnetic field \(\lambda=0.5,2.0,3.0\) for the three-qubit reduced density matrix using the same data as in Fig. 3. Even though QNEE mainly focuses on entropy estimation, its eigenvalue estimation is comparable to VQSE, which is designed for this particular task. A more detailed comparison between QNEE and VQSE will be discussed in the following subsection.

In Fig. 4 (d-f), we plot the learning curves of QNEE for the case of three qubits. The blue line represents the learning curve with the cost function calculated by the classical NN using a finite number of samples. The red curve represents the learning curve with the cost function obtained from the exact wave function of the ground state. As shown in Fig. 4, 100 \(\sim\) 200 iterations are enough to classify the phases. A similar number of iterations are also enough for the four-qubit case.

### Comparison between QNEE and VQSE

We compare the performance of QNEE with VQSE in more detail. We first clarify that the cost functions of QNEE and VQSE are designed for different purposes; the former is for entropy estimation, while the latter is for estimating \(\ell+1\) lowest eigenvalues (see Appendix C and Ref. [23] for more details). Therefore, it is expected that QNEE will provide better performance in the entropy estimation task, as observed in Fig. 3. On the other hand, the eigenvalue estimation by VQSE is, in general, more accurate than QNEE.

In our XXZ chain model, both QNEE and VQSE show good agreement in entropy/eigenvalue estimation in the non-critical regime (\(\lambda>\lambda_{c}\)), where the entropy is vanishing. Near the transition point (\(\lambda\sim\lambda_{c}\)), QNEE produces a

Figure 3: Plots of estimation results on the reduced density matrices of the XXZ chain. Total number of qubits is 8. The reduced-density matrix is prepared by partially tracing out (a) five qubits and (b) four qubits. The entanglement entropy is estimated with 5 different randomly initialized quantum circuits. ‘\(\times\)’ symbol represents the lowest value of cost functions among 5 trials. Insets show the mean and standard deviations of 5 trials. Orange-colored symbols represent the results of QNEE. Blue-colored symbols represent the results of VQSE. The exact value of entanglement entropy is denoted by a black line.
more reliable entropy estimation than VQSE, while the accuracies of the eigenvalue estimations are comparable. In the critical regime, (\(\lambda<\lambda_{c}\)) QNEE estimates the entropy more accurately than VQSE. In particular, for \(\lambda=0.5\), VQSE estimates the eigenvalues (see Fig. 4(a)) better than QNEE, but QNEE yields a higher accuracy of the estimated entropy than VQSE (see Fig. 3(a)). This could lead to an intriguing implication that, in certain cases, the von Neumann entropy can be accurately estimated using QNEE, even when the estimated eigenvalues do not exactly match the correct values. On the other hand, the main advantage of VQSE comes from its low computational cost. As VQSE does not involve a training process, it runs much faster than QNEE to estimate the quantum state's entropy and eigenvalues while yielding a comparable performance.

We also note that QNEE is less sensitive to the initial choices of \(\Theta_{Q}\) than VQSE. This can be captured by comparing the variance of estimated entropy by taking different initial points of \(\Theta_{Q}\) (see the insets of Fig. 3). The estimator's variance for QNEE ranges up to 0.09. In contrast, the variance of VQSE ranges up to 0.16, which is significantly higher than that of QNEE. From the simulation data, we further analyze the behavior of errors depending on the exact entropy value. In Fig. 5, we plot the absolute error of estimation results of QNEE \(|S(\hat{\rho})-C(\Theta_{N},\Theta_{Q})|\) as a function of entanglement entropy for 3- and 4-qubit reduced density matrices of the XXZ chain. QNEE results show a clear tendency where the error and its fluctuation decrease as the entanglement entropy decreases (see Fig. 5). While VQSE shares a similar tendency, this is not as clear as the QNEE case.

Figure 4: (a-c) Plots of estimated eigenvalues and (d-f) learning curves for the reduced density matrix with three qubits. These results are from three \(\lambda\): (a, d) \(\lambda=0.5\), (b, e) \(\lambda=2.0\), and (c, f) \(\lambda=3.0\). In the top plots, orange (blue) color denotes the estimation results of QNEE (VQSE). \(\times\) represents the exact (analytic) result. In the bottom plots, blue line represents cost function \(C^{\rm NN}\) trained with finite samples. The red line represents ideally trained \(C(\bar{\Theta}_{N},\Theta_{Q})\) without finite-sample effect. Here we utilize the same data that is used for Fig 3(a).

Figure 5: Scatter plot of the absolute errors of QNEE (orange symbol) and VQSE (blue symbol) as a function of entanglement entropy. For this plot, we take data from reduced density matrices of three and four qubits of the XXZ chain.
To summarize, QNEE performs better in entropy estimation, while VQSE is more specialized in finding the first few largest eigenvalues. We expect that such a comparison will also be valid for general cases with a larger number of qubits, where a considerable number of non-zero eigenvalues are involved in the quantum state's entropy. However, we must note that the observed tendency between the two methods is not definite and may depend on the quantum state of interest and the choice of variational parameters in quantum circuit ansatz and classical optimizers.

## IV Discussion

We have proposed a novel quantum entropy estimator, QNEE, which combines variational quantum circuits with classical NNs. Our method allows us to estimate the von Neumann entropy and the Renyi entropy from finite samples of circuit outcomes without accessing full information of the density matrix. By leveraging the Donsker-Varadhan inequality to optimize the cost function, classical NNs in the QNEE inherit valuable characteristics from recently developed classical entropy estimators, including NEEP [24]. We emphasize that classical NNs operate independently of the structure of variational quantum circuits, offering a universal application to any form of quantum circuit ansatz.

As a physical model, we have applied QNEE to estimate the entanglement entropy of the ground state. We have demonstrated that QNEE effectively classifies the massless and ferromagnetic phases of the 1D XXZ Heisenberg model. In particular, numerical evidence shows that QNEE can sensitively estimate the entanglement entropy near the critical point even before the estimated eigenvalues converge to the exact ones. We also highlight that such phase classification can still work even without accurately estimating entanglement entropy when one of the phases has low entanglement, as the estimator is designed to upper bound the exact entropy value.

On the other hand, there are some cautions when applying QNEE. Although the cost function provides an upper bound on the exact entropy value in an ideal scenario, the estimated value obtained from insufficient samples can be underestimated compared to the exact value. Such a case has already been reported when using classical stochastic entropy production estimators, including NEEP [24]. However, it has also been studied that the NN-based estimator still provides efficient estimation rather than other methods though [12]. Also, the performance of QNEE is sensitive to the learning rate of the gradient-descent method for the quantum circuit. A lower learning rate is recommended as the number of qubits increases or the number of shots decreases. This can be a potential issue when applying the QNEE to a high dimensional quantum state, as a high learning rate can induce fluctuations in the learning curve.

Our methodology establishes a pathway for investigating quantum state properties using fewer iterations of quantum circuits, which are considerably more resource-intensive than classical processing. We anticipate that QNEE will serve as a vital tool for addressing a broad spectrum of estimation and classification challenges within quantum systems.

###### Acknowledgements.

 The authors acknowledge the Korea Institute for Advanced Study for providing computing resources (KIAS Center for Advanced Computation Linux Cluster System). This research was supported by the KIAS Individual Grant Nos. PG081801 (S.L.), CG085301 (H.K.), and PG064901 (J.S.L.) at the Korea Institute for Advanced Study. S. L. thanks Dong Kyum Kim for helpful discussions about machine learning.

## Appendix A Derivation of the upper bound of the von Neumann entropy

Here, we derive the upper bound of the von Neumann entropy, i.e., Eq. (2). In terms of the eigenstate \(|\lambda_{m}\rangle\) of a density matrix \(\hat{\rho}\), that is \(\hat{\rho}=\sum_{m}\lambda_{m}|\lambda_{m}\rangle\langle\lambda_{m}|\), one can write an inequality for an observable \(\hat{O}\) as follows:

\[\ln\mathrm{tr}\{e^{\hat{O}}\}= \ln\sum_{m}\langle\lambda_{m}|e^{\hat{O}}|\lambda_{m}\rangle\] \[\geq \ln\sum_{m}e^{\langle\lambda_{m}|\hat{O}|\lambda_{m}\rangle}\] \[= \ln\sum_{m}\lambda_{m}e^{\langle\lambda_{m}|\hat{O}|\lambda_{m} \rangle-\ln\lambda_{m}}\] \[\geq \ln e^{\sum_{m}\lambda_{m}\left((\lambda_{m}|\hat{O}|\lambda_{m} \rangle-\ln\lambda_{m}\right)}\] \[= \langle\hat{O}\rangle_{\hat{\rho}}+S(\hat{\rho}). \tag{10}\]
Convexity of the exponential function is used for deriving the inequalities in the second and the fourth lines of Eq. (17). We note that the second-line inequality is saturated when the observable operator is diagonalized by the basis set \(\{|\lambda_{m}\rangle\}\). Equation (17) can be rearranged as

\[S(\hat{\rho})\leq-\langle\hat{O}\rangle_{\hat{\rho}}+\ln\mathrm{tr}\{e^{\hat{O }}\}, \tag{18}\]

which represents the Gibbs variational principle. In general, the Gibbs variational principle focuses on finding an optimal density matrix, enabling Eq. (18) to be saturated, which is given by the Gibbs state \(\hat{\rho}^{*}=e^{\hat{O}}/\mathrm{tr}\{e^{\hat{O}}\}\). In contrast, QNEE focuses on finding an optimal observable achieving the equality of Eq. (18), which is given by \(\hat{O}^{*}=\ln\hat{\rho}\).

For implementing the variational method to obtain the optimal entropy value via NN, one might be able to use the right-hand side (rhs) of Eq. (18) as a cost function. However, when employing the stochastic (mini-batch) gradient descent method for optimizing the parameters of the NN, the cost function, including a logarithmic function such as the rhs of Eq. (18) can malfunction and yield an incorrect result without additional data manipulation [30]. To avoid this complication, we linearize the logarithmic function by applying the additional inequality, \(\mathrm{tr}\{e^{\hat{O}}\}-1\geq\ln\mathrm{tr}\{e^{\hat{O}}\}\), which is also saturated when \(\hat{O}^{*}=\ln\hat{\rho}\). Application of this additional inequality to Eq. (18) leads to Eq. (2) as follows:

\[S(\hat{\rho})\leq-\langle\hat{O}\rangle_{\hat{\rho}}+\mathrm{tr}\{e^{\hat{O}} \}-1. \tag{19}\]

## Appendix B Derivation of the upper bound of the Renyi entropy

Here, we derive the upper bound of the Renyi entropy in Eq. (6). For a convex and twice-differentiable real-valued function \(f(x)\) with arbitrary real numbers \(p\) and \(q\) (\(\neq 0\)), the following inequality

\[-qf(p/q)\leq-pf^{\prime}(x)+q\left\{xf^{\prime}(x)-f(x)\right\} \tag{20}\]

holds [33]. We can demonstrate this by differentiating the rhs of Eq. (20) with respect to \(x\), which yields \(f^{\prime\prime}(x)(-p+qx)\). As \(f^{\prime\prime}(x)>0\) by convexity, the rhs has the unique minimum \(-qf(p/q)\) at \(x=p/q\). As Eq. (20) is valid for arbitrary \(p\) and \(q\), we can replace \(p\) and \(q\) with some distributions \(P(s_{i})\) and \(Q(s_{i})\), where \(\sum_{i}P(s_{i})=\sum_{i}Q(s_{i})=1\). By considering that \(x=x(s_{i};\theta)\) is a function of \(s_{i}\) and parameterized by \(\theta\), and summing over all \(i\), we can write Eq. (20) as the following variational form:

\[-D_{f}(P\|Q)\leq\sum_{i}\left[-P(s_{i})f^{\prime}\left(x(s_{i};\theta)\right) +Q(s_{i})\left\{x(s_{i};\theta)f^{\prime}(x(s_{i};\theta))-f(x(s_{i};\theta) )\right\}\right], \tag{21}\]

where \(D_{f}(P\|Q)\equiv\sum_{i}Q(s_{i})f\left(P(s_{i})/Q(s_{i})\right)\) is the \(f\)-divergence. To connect this inequality to the Renyi entropy, we choose the convex function \(f(x)\) as

\[f(x)=\frac{x^{\alpha}-\alpha x-(1-\alpha)}{\alpha(\alpha-1)} \tag{22}\]

for \(\alpha\neq 0,1\). Then, the \(f\)-divergence can be rewritten as

\[D_{f}(P\|Q)=\frac{1}{\alpha(\alpha-1)}\left[\sum_{i}\frac{P^{\alpha}(s_{i})} {Q^{\alpha-1}(s_{i})}-1\right]. \tag{23}\]

We can make a relation between the \(f\)-divergence and the Renyi-divergence \(D_{\alpha}(P\|Q)\) as

\[D_{\alpha}(P\|Q)\equiv \frac{1}{\alpha-1}\ln\left[\sum_{i}\frac{P^{\alpha}(s_{i})}{Q^{ \alpha-1}(s_{i})}\right]\] \[= \frac{1}{\alpha-1}\ln\left[\alpha(\alpha-1)D_{f}(P\|Q)+1\right]. \tag{24}\]

By choosing \(Q\) as the uniform distribution, i.e., \(Q(s_{i})=1/d\) for all \(i\) with \(d=\sum_{i}1\), the Renyi divergence is reduced to the Renyi entropy as

\[D_{\alpha}(P\|Q)=-H_{\alpha}(P)+\ln d, \tag{25}\]
where the Renyi entropy is defined as \(H_{\alpha}(P)\equiv\frac{1}{1-\alpha}\ln\sum_{i}P^{\alpha}(s_{i})\). Plugging Eqs. (3), (4), and (5) into Eq. (4), and substituting \(x(s_{i};\theta)\) with \(Ne^{h(s_{i};\theta)}\), Eq. (4) can be arranged as

\[\frac{e^{(1-\alpha)H_{\alpha}(P)}}{\alpha(1-\alpha)}\leq\sum_{i}P(s_{i})\frac{ e^{(\alpha-1)h(s_{i};\theta)}}{1-\alpha}+\sum_{i}\frac{e^{\alpha h(s_{i};\theta)}} {\alpha}. \tag{6}\]

The quantum generalization of the Renyi entropy for a given density matrix \(\hat{\rho}\) is defined as \(S_{\alpha}(\hat{\rho})\equiv\frac{1}{1-\alpha}\ln\mathrm{tr}\{\hat{\rho}^{ \alpha}\}=\frac{1}{1-\alpha}\ln\left(\sum_{i}\lambda_{i}^{\alpha}\right)\). By noting that the probability distribution \(P_{i}=\mathrm{tr}\{\hat{\rho}\Pi_{i}\}\) for a set of measurement operators \(\{\Pi_{i}\}\) is always majorized by the distribution of eigenvalues \(\lambda_{i}\) and that the Renyi entropy is a Shcur-concave function, we obtain

\[H_{\alpha}(P)\geq S_{\alpha}(\hat{\rho}).\]

In particular, the inequality is saturated when the measurement basis is the eigenbasis of the density matrix \(\hat{\rho}\). By taking \(P_{\hat{V}}(s_{i};\Theta_{Q})=\langle s_{i}|\hat{V}(\Theta_{Q})\hat{\rho}\hat {V}^{\dagger}(\Theta_{Q})|s_{i}\rangle=\mathrm{tr}\{\hat{\rho}_{\hat{V}}\Pi_{i}\}\) and subtracting \(1/[\alpha(1-\alpha)]\) from both sides of Eq. (6), we can reach the final result as

\[\frac{e^{(1-\alpha)S_{\alpha}(\hat{\rho})}-1}{\alpha(1-\alpha)}\leq\frac{e^{( 1-\alpha)H_{\alpha}(P_{\hat{V}}(s_{i};\Theta_{Q}))}-1}{\alpha(1-\alpha)}\leq \sum_{i}P_{\hat{V}}(s_{i};\Theta_{Q})\frac{e^{(\alpha-1)h(s_{i};\Theta_{N})}- 1}{(1-\alpha)}+\sum_{i}\frac{e^{\alpha h(s_{i};\Theta_{N})}-1}{\alpha}, \tag{7}\]

which is Eq. (6). Note that Eq. (7) is reduced to the inequality for the von Neumann entropy (4) in the \(\alpha\to 1\) limit.

## Appendix C details of VQSE

In this section, we illustrate the details of VQSE [23]. VQSE exploits the relationship between majorization and diagonalization. The cost function is given as the expectation value of an artificial Hamiltonian

\[\hat{H}\equiv (1-t)\hat{H}_{L}+t\hat{H}_{G}, \tag{8}\]

where \(\hat{H}_{L}\equiv\hat{I}-\frac{1}{2}\sum_{j=1}^{\ell}r_{j}\hat{\sigma}_{j}^{z}\) is local Hamiltonian that is designed to mitigate Barren plateau, and \(\hat{H}_{G}\equiv\hat{I}-\sum_{s_{j}\in S}q_{j}|s_{j}\rangle\langle s_{j}|\) is a global Hamiltonian. \(t\), \(r_{j}\), and \(q_{j}\) are real number parameters. \(S\) is the set of \(m\) most numerous strings from the quantum circuit, where \(m\) is a positive integer less than \(2^{\ell}\). At the start of the optimization process, \(t=0\) and \(t\) increases at every certain iteration up to \(t=1\). The set \(S\) is also updated with the same period with \(t\). \(\hat{I}\) is the identity operator, and \(\ell\) is the number of qubits. There exist many choices for \(q_{j}\) and \(r_{j}\). The choice used in reference [23] is \(r_{j}=r_{1}+(j-1)\delta\), and the eigenvalues of \(H_{G}\) are set to be the same value as \(\ell+1\) lowest eigenvalues of \(\hat{H}_{L}\). The other eigenvalues of \(\hat{H}_{G}\) are set to be 1. This choice of \(r_{j}\) ensures \(\ell+1\) lowest non-degenerate eigenvalues of \(H_{L}\). Therefore, \(\ell+1\) eigenvalues can be estimated. In our research, we set \(r_{1}=0.2\), \(\delta=0.01\), and \(\ell=n\). We set the learning rate for training VQSE to be 0.05 and update \(t\) and \(S\) at every 25 iterations.
