High-Speed Resource Allocation Algorithm Using a Coherent Ising Machine for NOMA Systems

The high complexity of the RA problem in NOMA systems can be reduced using heuristic methods by dividing the problem into sub-problems [12, 9, 10, 11, 5]. In [5], the authors propose low complexity and efficient method for RA in a downlink NOMA system. In this study, user fairness, data rate, and energy efficiency are optimized. The RA problem is divided into the sub-problems of power allocation and channel allocation. The power allocation solution is derived from a theoretical analysis. The channel allocation is carried out based on the derived power allocation solution. In [9], the RA in full duplex-NOMA (FD-NOMA) systems is studied. In this study, the RA problem is successfully solved with low computational complexity by dividing the channel and power problems based on the block coordinate descent method. In [10], an algorithm for downlink NOMA systems that analytically examines the optimal power allocation (OPA) and ideal pairing search algorithm (IPSA) in terms of fairness was proposed. Simulation results indicate that the computational complexity is significantly reduced in comparison to the ES. In [11], NOMA-based device-to-device (D2D) communications are considered. It is shown that the communication performance can be improved to a significant extent using the NOMA technique. The RA solution is obtained with low computational complexity by applying an optimal power allocation after a channel allocation algorithm that assigns sub-carriers to the D2D groups. In [12], the RA in a NOMA system is optimized using the SA, which is a meta-heuristic algorithm. In this study, a joint optimization of the channel allocation, sub-channel power allocation, and inter-user power allocation is considered. Simulation results indicate that a sufficiently reliable solution in terms of system throughput can be obtained with low computational complexity in comparison to an ES.

As discussed above, we can see that the heuristic RA approaches can obtain a solution with low computational complexity. However, there is a fundamental trade-off between the performance of the solution and the computational complexity, and there have been no previous studies that can quickly optimize on the order of milliseconds.

With the development of ML applications in wireless communications, ML solutions for RA optimization problems in NOMA systems have been proposed [4, 13, 14]. In [4], deep reinforcement learning (DRL) based RA scheme is proposed for NOMA systems. Here, the framework of the DRL method is designed based on an attention-based neural network (ANN). Specifically, the joint optimization problem of channel and power allocation is formulated, and the channel allocation is learned according to the policy learned from the ANN. The simulation results indicate that the proposed method can achieve the same performance for six users with lower computational complexity in comparison to the ES. In [13], a DQN is adopted to achieve a joint RA for a large number of users, where reinforcement learning (RL) frames are used for the joint optimization of clustering, power allocation, and beamforming. Specifically, the clustering is adjusted using a DQN, and a backpropagation neural network (BPNN) is designed to learn the power coefficients for each cluster. Simulation results indicate that the performance of the proposed method is close to that of the ES. In [14], a method based on DRL is proposed to assign users to time-varying channels. Recently, DRL has been expected to be applied as a solution for such real-time dynamic decision-making problems in practical environments where the conventional RA algorithms are difficult to be applied. A complexity analysis and simulation results show that the proposed approach provides a superior performance with lower computational complexity in comparison to conventional optimization algorithms.

Therefore, ML-based RA methods can obtain near-optimal solutions with lower computational complexity than the ES. However, to obtain near-optimal solutions, a proper training time is required.

Optimization problems such as channel and power allocation are important not only in NOMA systems but also in a variety of other wireless communication systems. Optimization problems in wireless communication systems have recently become increasingly complex, and the demand for efficient solutions to meet the requirements of next-generation wireless communications has increased [28]. Hardware-based approaches, such as quantum computers and Ising machines, have been expected to be potential methods for solving optimization problems at high speed [15, 16, 17, 18]. In previous studies, several RA optimization problems in wireless communications have been studied using the CIM [22, 23, 24, 25]. In [23], the CIM is applied to RA problems in wireless communication systems such as IEEE802.11 wireless LAN systems and D2D communications. In addition, it has been shown that the CIM can approximately solve these optimization problems within the order of milliseconds. In [24], the CIM is applied to the optimization of scheduling problems for distributed antenna systems. In [25], the CIM is applied to the optimization problem of channel assignment in dense wireless LAN (WLAN) systems. In this study, the maximization of the system throughput in high-density WLAN systems using the CIM is considered. The performance of the proposed CIM-based channel assignment approach in comparison to other optimization approaches such as the SA and greedy algorithms is evaluated. Simulation results indicate that the CIM can obtain the optimal solution well with constant computational time, whereas other optimization algorithms increase the computational time and decrease the accuracy of the obtained solution as the number of users in the system increases.

In these studies, the CIM is assumed to work either as a base station (BS) or as a controller in the centralized systems. The ability of the CIM to operate at room temperature [29] and its simple components make it suitable for integration into many wireless transmission systems. Hence, the CIM is significantly more suitable for large-scale wireless networks if the size of the problem is up to the number of spins that can be achieved with the CIM.

In summary, the CIM may obtain optimal solutions to combinatorial optimization problems without the need for complex computer calculations. The principle of searching for optimal solutions using the CIM is based on the spontaneous energy reduction in physical phenomena. Existing research has demonstrated the higher speed of a CIM-based optimization. On the contrary, DRL or other ML algorithms require training samples to obtain the optimal solution, and trial and error to increase performance. In a real communication environment, such as a vehicle network, the environment changes frequently, and real-time optimization may not be possible for the ML algorithms because retraining is needed when the environment changes. In [4, 14], the complexity due to training in the ML algorithm was studied. In contrast to the complexity of training presented in these papers, the proposed CIM-based RA method can always solve the optimization problem at constant high speed, i.e., millisecond order [26]. Therefore, real-time RA can be achieved using our proposed method, which may improve the spectrum efficiency under dynamic changing environments.

The CIM is an Ising-type machine that can artificially reproduce an Ising model, a physical model of magnetic spins. The Ising model consists of two states of magnetic spins: upward and downward. Spins are mutually coupled and subject to interactions from other spins and effects from external magnetic fields. We consider Ising spins with an $N\times M$ two-dimensional structure. Here, the energy function of the Ising model, i.e., the Ising Hamiltonian is expressed as follows:

where $\sigma_{ij}\in\{-1,+1\}$ is the spin direction of the $(i,j)$th spin, $J_{ijkl}$ is the strength of the interaction between the $(i,j)$th and $(k,l)$th spins, and $\lambda_{ij}$ is the strength of the external magnetic field on the $(i,j)$th spin. Here, the CIM is the machine that can obtain the ground state of the Ising Hamiltonian at a high speed. In other words, by setting $J$ and $\lambda$ to correspond to the optimal solution of the optimization problem to the minimum of Eq. (15), we can approximately obtain a fast optimal solution using the CIM.

Initially, the CIM was a laser network consisting of one master laser and multiple slave lasers[27]. Currently, a measurement-feedback type CIM has been proposed to easily increase the number of couplings between spins to solve larger scale problems[16, 17]. A schematic diagram of the measurement feedback-type CIM is shown in Fig. 2. As shown in Fig. 2, parameters of the Ising model $J_{ijkl}$ and the $\lambda_{ij}$, parameters of the Ising model are pre-configured in a field programmable gate array (FPGA) module. The FPGA module calculates the feedback value using the pre-configured parameters and feeds it back to the original pulse. Therefore, interactions between spins can be programmatically reproduced, which makes it possible to fully connect spins. In the measurement-feedback CIM, the phase of the degenerate optical parametric oscillator (DOPO) pulse circulating on an optical fiber represents Ising spins $\sigma_{ij}$. In [19], it is shown that 100,000 DOPO pulses can be realized over 5 km of polarization-maintaining optical fiber and reach the reference score of the MAX-CUT problem for 100,000 nodes within a time of 593 µs.

The CIM is an Ising machine that can rapidly obtain the ground state of the Ising Hamiltonian by getting the combinations of the Ising spin directions. The set of the Ising spin directions that minimize the Ising Hamiltonian corresponds to the solution of the combinational optimization problem. The operation process of the CIM can be summarized as follows. As shown in Fig. 2, the initial pulses are first generated by a phase sensitive amplifier (PSA) on the optical fiber. The initial pulses are repeatedly amplified by PSA with a pump pulse each cycle and become a DOPO pulse. The phase of the DOPO pulse circulating on the fiber during oscillation is used to reproduce the Ising spin on the optical fiber in the CIM. Here, gradually increasing the pump pulse amplitude until the threshold is exceeded; the DOPO pulse will oscillate and take on a $\pi$ or 0 phase representing the orientation of the Ising spin $\sigma_{ij}$, where $\sigma_{ij}\in\{-1,+1\}$. 0 phase represents the up spin, i.e., $\sigma_{ij}=1$. Meanwhile, the $\pi$ phase represents the down spin, i.e., $\sigma_{ij}=-1$ While the DOPO pulse is being amplified, a part of the DOPO pulse on the fiber is separated by a coupler, and the in-phase amplitude is measured through a balanced homodyne detection (BHD). The measured in-phase amplitude $\sigma_{ij}$ is input into the FPGA module. A feedback signal $-\sum{J_{ijkl}\sigma_{ij}}+\lambda_{ij}$ is calculated using the mutual coupling between Ising spins $J_{ijkl}$ and external magnetic field $\lambda_{ij}$. Finally, the feedback signal is injected into the original DOPO pulse as a feedback pulse through the intensity modulator and push-pull modulator. By repeating the operation process of the CIM in this manner, we can eventually obtain a combination of $\sigma_{ij}$ that minimizes Eq. (15), which corresponds to the optimal solution to the optimization problem.

From the above, to obtain the optimal solution of optimization problem using the CIM, we need to transform the objective function into the form of an Ising Hamiltonian, and then derive the corresponding $J_{ijkl}$ and $\lambda_{ij}$ of the objective function. Then, by pre-configuring the derived these parameters into the FPGA module of the CIM, we can search for the ground state of the Ising Hamiltonian, which corresponds to the optimal solution to the optimization problem. In the next section, we introduce the problem transformation and the derivation of the corresponding $J_{ijkl}$ and $\lambda_{ij}$ for our formulated objective function in Eq. (11).

Let us consider a neural network with a two-dimensional structure of $N\times M$; the mutually connected neural network has the following energy structure:

where $x_{ij}\in\{0,1\}$ is the state of the $(i,j)$th neuron, $w_{ijkl}$ is the coupling weight between the $(i,j)$th and $(k,l)$th neurons, and $\theta_{ij}$ is the firing threshold of the $(i,j)$th neuron. From the similarity of the energy structure of the Ising Hamiltonian and that of the mutually connected neural network, i.e., Eqs. (15) and (16), the parameters of the CIM, $J_{ijkl}$ and $\lambda_{ij}$, can be derived using the parameters (i.e., $w_{ijkl}$ and $\theta_{ij}$) of the mutually coupled neural network. Here, the state of a neuron in a network is defined by the following equation:

where $u$ is the Heaviside step function, namely, $u[y]=0$ for $y\leq 0$, and $u[y]=1$ for $y>0$. The energy function of the mutually connected neural network, i.e., Eq. (16) always decreases and converges to a locally minimum value by updating the neurons iteratively according to Eq. (17) if the following conditions are satisfied. 1). The self-connections of all neurons are zero, i.e., $w_{ijij}=0$. 2). The mutual connections are symmetric, i.e., $w_{ijkl}=w_{klij}$.

To solve the optimization problem using a mutually connected neural network, we must formulate the optimization problem as an objective function using a neuron $x_{ij}$. Then, the formulated objective function is compared with Eq. (16) to derive $w_{ijkl}$ and $\theta_{ij}$. By iteratively updating the neurons according to Eq. (17) using the derived parameters $w_{ijkl}$ and $\theta_{ij}$, Eq. (16) converges to a local minimum value, and the states of the neurons at that time correspond to the approximate optimal solution to the optimization problem. In the next subsection, we define neurons $x_{ij}$ and derive $w_{ijkl}$ and $\theta_{ij}$ for optimizing the RA in a NOMA system.

To optimize the RA problem in a NOMA system using a mutually connected neural network, we define the binary variables in Eq. (12) as neurons. Then, our formulated object function Eq. (11) can be transformed into the energy function of the mutually connected neural network, which is expressed as follows:

where $\delta_{ij}$ is Kronecker’s delta, namely, $\delta_{ij}=1$ when $i=j$, and $\delta_{ij}=0$, otherwise. In addition, the constraints (8.c) and (8.d) can be formulated as constraint terms using neurons, which can be expressed as follows:

Note that the above constraints term can be fully satisfied even if the number of users in the cell $N_{u}$ is less than $2N_{c}$, since we define number of the dummy users $N_{d}$ to always satisfy $N_{u}+N_{d}=2N_{c}$. Using $E_{1}$, $E_{2}$, and $E_{3}$, the energy function of the RA optimization problem in a NOMA system can be obtained as follows:

where $\epsilon$, $\zeta$, and $\eta$ are parameters that adjust the scaling of each term, which are used to adjust the objective and constraint terms, causing the energy structure to more likely reach the optimal solution. By comparing Eq. (21) with Eq. (16), the coupling weight $w_{ijkl}$ and firing threshold $\theta_{ij}$ for solving the RA optimization problem in a NOMA system using mutually connected neural networks can be obtained as follows:

Using $w_{ijkl}$ and $\theta_{ij}$ obtained in Eqs. (22) and (23), the neurons were updated according to Eq. (17). Thereby, the combination of neurons minimizes the energy function in Eq. (16), i.e., the optimal solution to the NOMA RA problem can be obtained. However, owing to the monotonically decreasing nature of the energy function, this neural network is likely to fall into a local minimum of the network state. Hence, the optimal solutions may not be achieved using the mutually connected neural network. Further, it has been shown that the Ising Hamiltonian near the optimal solution can be reached with high probability by amplifying the DOPO pulse to near the oscillation threshold in CIM [31]. Therefore, in the next subsection, we will use the derived parameters $w_{ijkl}$ and $\theta_{ij}$ to derive $J_{ijkl}$ and $\lambda_{ij}$ which are the parameters used to solve the optimization problem in the CIM.

As previously noted, the Ising Hamiltonian of the CIM and the energy function of the mutually connected neural network have similar energy structures. As a difference, the neuron $x_{ij}\in\{0,1\}$ in the mutually connected neural network has states of 0 and 1, whereas the Ising-spin $\sigma_{ij}\in\{-1,+1\}$ in the Ising model has states of -1 and +1. To solve the RA problem using the CIM, we redefine the output of the neurons in a mutually connected neural network as -1 and +1. Using the neuron $\tilde{x}_{ij}\in\{-1,+1\}$ with a redefined output, the update equation for the neuron is expressed as follows:

where $r[y]=-1$ for $y\leq 0$, and $r[y]=+1$ for $y>0$, and $\tilde{w}_{ijkl}$ and $\tilde{\theta}_{ij}$ are the coupling weight and firing threshold of the neural network with output $\{-1,+1\}$, respectively. We then can obtain the following equation by transforming the internal state of Eq. (17) using $\tilde{x}_{ij}=2x_{ij}-1$:

By comparing Eqs. (16) and (25), we can obtain the following equations for $\tilde{w}_{ijkl}$ and $\tilde{\theta}_{ij}$, respectively.

By setting these derived parameters in the FPGA module, a fast RA optimization of the NOMA system using the CIM can be performed.

Then, we present the details of the optimization procedure of RA using the CIM in NOMA systems. In this study, the BS performs channel allocation using the CIM for the NOMA system to achieve the optimal total reachable data rate. A schematic of the optimization procedure of the NOMA system is shown in Fig. 3. For simplicity, we present a NOMA system with 4 users and 2 channels as an example. As shown in Fig. 3, each user transmits the data to the BS. Then, the channel assignment is determined by the CIM equipped in the BS. There are DOPO pulses on the optical fiber of the CIM, each of which behaves as an Ising spin. Because this problem involves 4 users and 2 channels, the RA problem can be solved using the CIM implemented with $2\times 4=8$ DOPO pulses on the fiber. When the CIM is operated, the states of the in-phase components shown on the left side of Fig. 3 are obtained. Here, the cycle number represents the number of times the DOPO pulse travels around the optical fiber of the CIM. As the cycle count increases, the in-phase component of the amplitude of the DOPO pulses increases and eventually splits into positive or negative values, i.e., 0 phase or $\pi$ phase. In the CIM, the spin orientation is determined by the in-phase components of the amplitude. Here, when user $i$ uses channel $j$, the spin is +1, i.e., an upward spin. Thus, the obtained spin combination corresponds to the channel assignment. Finally, based on this channel assignment, the BS sends a signal to the users after assigning the optimal power using Eq. (13). The pseudo-code for the algorithm is presented as Algorithm 1.

In this subsection, we discuss the time and computational complexity of the proposed method. In [26], it is demonstrated that the CIM takes a constant amount of time to obtain the ground state as long as the fiber length is constant. That is, the time complexity of the CIM is $\mathcal{O}(1)$ experimentally. The reason is that the computation time of one cycle for the CIM is related to the length of a fiber. For the CIM with 2000 spins and 1km fiber, 1000 cycles running of the CIM can be achieved within 5 ms while the spins can converge to the ground state within 1000 cycles by amplifying pump pulse amplitude appropriately. In our paper, the problems can be solved using the CIM with 2000 spins. Hence, the computation time is within 5 ms.

Then, we discuss the computation complexity of the main computation part of the CIM, i.e., FPGA. Ref. [16] shows that the memory resources required for matrix computation and power consumption in FPGAs are scaled by $\mathcal{O}(N^{2})$. A typical CPU runs 10-100Gflop per second, so its power efficiency is often below 1 Gop/J. Similarly, GPUs also suffer from high power consumption. Hence, CPUs and GPUs are difficult to apply to applications that require low power consumption [39]. In addition, although parallel processing is possible for modern CPUs with multi-core, unnecessary parts also need to be parallelized since the entire core is used. In the case of FPGA, only necessary processing needs to be programmed. Hence, parallel processing can be performed without wasting computation resources, increasing the computation efficiency. Moreover, owing to the merit of the parallel processing of the FPGA, the FPGA uses less circuit area when performing the same processing compared to the CPU/GPU and can perform processing without resource waste. It is natural with high power efficiency. In summary, our proposed method may be superior to the resource allocation method computed using the traditional computer in computation complexity.

In the simulation, a circular cell of a NOMA system with a radius of 500 m is considered. A BS is placed at the center of the cell. Users are randomly placed in the cell. Table II shows the detailed parameter setting used in our simulation.

We use the following simulation model of the CIM to evaluate our proposed method:

where $c_{ij}$ and $s_{ij}$ are the in-phase and quadrature-phase amplitude of the $(i,j)$th optical pulses, respectively, and $p$ is the pump pulse, which is used to amplify $c_{ij}$. By running the simulation with a gradual increase in $p$, the in-phase and quadrature components of the DOPO pulse can be obtained. As described in Section \@slowromancapiv@, the measurement feedback of the CIM reproduces the Ising spins using the amplitudes of DOPO pulses cycling over long fibers. In other words, the simulation model can be used to reproduce the behavior of the CIM and obtain the combination of spins that minimizes the Hamiltonian. Ising spin $\sigma_{ij}$ is implemented by $c_{ij}$: If $c_{ij}<0$, then $\sigma_{ij}=-1$, whereas if $c_{ij}>0$, then $\sigma_{ij}=+1$.

First, we evaluated the convergence of the proposed method in the NOMA system with 12 users and 6 subchannels, i.e., $N_{u}=12$ and $N_{c}=6$. Figure 4 shows the time evolution of the in-phase components of the DOPO pulses for each cycle and the corresponding energy of the Ising Hamiltonian at that time. Here, for the optimization time, we used the time required for 1000 cycles on an actual CIM machine as shown in [16].

The subfigures on the left side, in the center, on the right side of Fig. 4 are the in-phase amplitude of the spins, the corresponding Ising Hamiltonian, and the total data rate of the three RA methods in the NOMA system, respectively. From Fig. 4, we can see that when the in-phase amplitude of the spins converges, the minimum Ising Hamiltonian reaches, while the total data rate of the CIM-based RA achieves the optimal value. This means that in the process of the CIM operations, the spins spontaneously selected a combination of spins that exhibited the ground state of the Ising Hamiltonian, which corresponds to the solutions to the formulated optimization problem. Hence, the convergence of our proposed method has been verified.

In addition, it is shown in [19] that for certain MAX-CUT problems, the real CIM machines can provide a stable state approximately 1000 times faster than SA using state-of-the-art CPUs. Thus, our proposed method has a significant advantage in convergence speed. Moreover, in [40], it has been demonstrated that the phase of the DOPO pulse during oscillation converges to the ground state within 5 ms by amplifying the DOPO pulse with an appropriate pump pulse amplitude. The setting of the pump pulse amplitude in real CIM is out of the scope of this manuscript, while considering the space of the manuscript, we omit the detailed analytical proof process, which can be found in [19], [32], [40]. In this paper, we set the pump pulse amplitude, i.e., $p$ in equations (28) and (29), to an appropriate value by parameter search during the simulations. Thus, the in-phase amplitude of the DOPO pulse can converge within 5ms, as shown in Fig. 4.

Then, we evaluate the computation time of the SA, the ES, the DQN-based, and the CIM-based RA method under the varying number of users. Here, the computation time of the SA and the ES is evaluated using an Intel Xeon Gold 5222 CPU @ 3.80 GHz and C language, and the DQN-based RA is evaluated using the same CPU and Python with PyTorch, respectively. The number of users in the NOMA system $N_{u}$ is set as 6, 8, 10, 12, 14, 16, and 18 in the performance evaluation. The numbers of SA iterations under corresponding different $N_{u}$ were set as 10, 50, 100, 500, 1000, 5000, and 10000 times, respectively. In addition, the number of training episodes for the DQN-based RA under different $N_{u}$ is set as 100, 250, 400, 550, 700, 850, and 1000 times, respectively. Here, each episode consists of 120 steps. The BS (Base Station) collects information and trains the neural network for each step. The information consists of the CSI information, the channel assignment information, and the achieved data rate of the previous step of each user. The reason for the settings of the SA and DQN-based RA is that the larger the scale of the problem is, the larger number of iterations and training episodes required for SA and DQN-based RA methods. Here, owing to the enormous amount of computation time involved for the ES with the growth of the number of users, the method was evaluated for up to 14 users. Fig. 9 shows the evaluation results. The computation time of the real CIM in Fig. 9 is the time cost for RA using the CIM in the real world. In [16], it has been shown that the Ising problem can be solved using the CIM in 5 ms regardless of the problem size when the length of the fiber is 1-km-long. The computation time of DQN (train) and DQN (test) represents that for training the neural network and testing (making decisions) using the trained neural network, respectively. As shown in Fig. 9, the CIM-based and DQN-based RA for testing methods demonstrate almost constant computation time even as the number of users in the system increases. In contrast, other algorithms require a much longer time as the problem size increases. However, the DQN-based RA method requires training to obtain a good solution when the communication environment changes. The training time increases with the number of users in the NOMA system. Thus, as described above, the CIM is superior in speed compared to other methods.

Finally, we evaluate the performance of the proposed method when more realistic user activation is considered. In the simulation, the user activity follows the setup for the urban case in 3GPP TR 36.885 [34, 35]. In this scenario, the initial locations of the users are randomly generated according to the spatial Poisson process. Each user drives on the road at a constant speed. The speed of each user is generated randomly between 36km/h to 54km/h. The number of lanes is 2 in each direction and 4 in total in each street. The road grid size by the distance between intersections is 433m*250m. The simulation area size is 1299m*750m. When the user moves to the intersection, the direction is changed with a certain probability that follows the uniform distribution. In the performance evaluation, the number of users is set to 10, 12, 14, and 16. The time interval between each decision making is set to 20ms [14]. In this scenario, the wireless environment changes faster than the resource allocation using the ES, SA, and DQN (train) methods due to their long computation time, as shown in Fig. 9. These methods may no longer apply to mobile NOMA systems. Hence, we only evaluate the performance of the CIM-based and the DQN (test)-based RA methods under these settings. The evaluation results in terms of the ratio of the total data rate and the optimal solution are shown in Fig. 10. From Fig. 10, we can see that the CIM-based RA can achieve a higher data rate than the DQN (test)-based RA regardless of the number of users in the NOMA system. The reason can be summarized as follows. For the DQN (test)-based RA method, the time consumption for resource allocation is longer than the time interval between the decision-making. However, the user location and channel information are changed after the time interval. Hence, the information used by the DQN (test)-based RA method is the previous information that is different from the current information when allocating resources. On the other hand, the time consumption for the CIM-based RA method is shorter than the time interval between the decision-making. The CIM-based RA method can allocate resources using the current information. Hence, the optimal solution of the RA can almost be achieved by the CIM-based RA method. In summary, the proposed CIM-based RA method is applicable to the NOMA system with mobile users.

In this subsection, we discuss multiplexing more users (more than 2) in the same channel using the proposed CIM-based RA method for the NOMA system. This paper discusses NOMA systems where up to 2 users are multiplexed into the same channel. That is: in Eq. (11.d), we set a constraint term to limit the maximum number of users that are multiplexing in the same channel to up to 2. We believe that fast channel allocation in NOMA systems with more multiplexed users can be achieved by extending the constraint term in Eq. (11.d) as an arbitrary number of users $n$. Specifically, given the constraint, the energy function in Eq. (20) can be calculated. Finally, by recomputing the Ising model parameters, $J_{ijkl}$ and $\lambda_{ij}$, RA in the NOMA system with $n$ users multiplexed can be realized. Note that the number of dummy users $N_{d}$ to must be set as $N_{u}+N_{d}=nN_{c}$ to satisfy the above constraint. In addition, the number of channel assignment variables $x_{ij}$, $x_{kj}$ in Eq. (11) needs to be increased to correspond to the number of users multiplexed into the same channel.

In this subsection, we discuss the scalability of the proposed CIM-based RA methods in the practical settings with 2 or more users in the NOMA systems. As previously described, the CIM is the machine that can obtain the ground state of the Ising Hamiltonian at high speed by getting the optimal combination of the spin directions. The ground state of a quantum-mechanical system is its stationary state of lowest energy, which corresponds to the optimal solution of combinatorial optimization problems. In this paper, the combinatorial optimization problem is the RA problem in the NOMA system. The RA problem in the NOMA system can be solved by transforming the formulated RA combinatorial optimization problem into the form of the Ising Hamiltonian. The spin direction of the spin in the Ising Hamiltonian corresponds to the channel assignment variable of the RA problem in the NOMA system, i.e., whether the channel is allocated to the user or not. Each channel-user pair corresponds to one spin. The value of each spin direction is +1 or -1. If the spin direction is +1, the channel is allocated to the corresponding user. Otherwise, the channel is not allocated to the corresponding user. For example, channel $i$-user $j$ pair corresponding to the spin direction $\sigma_{ij}$. If $\sigma_{ij}$, the channel $i$ is allocated to user $j$. Otherwise, the channel $i$ is not allocated to user $j$. Hence, to solve the RA optimization problem in the NOMA system, $N_{u}*N_{c}$ spins are necessary, where $N_{u}$ and $N_{c}$ are the numbers of the users and channels, respectively. The related work [19] shows that the CIM can realize problems up to 100000 spins. This means that up to 446 users can be realized using the present CIM in our scenario, where up to 2 users are allocated to one channel simultaneously. In general, however, clustering with an arbitrary number of $n$ users of multiplexing is more attractive. Here, to achieve multiplexing at an arbitrary user $n$ using CIM, we need to change the definition of the assignment variable as described in the last subsection. Specifically, the number of indexes should be increased by considering which user $i$ is multiplexed with which user on channel $j$ in Eq. (12). Therefore, as n increases, the number of spins required to solve the problem also increases. For example, when $n=3$, the required number of spins is $N_{u}^{3}/3$, and RA for up to 67 users can be realized using the present CIM with 100000 spins. With the development of the CIM, much larger-scale optimization problems are expected to be solved by increasing the number of spins in the future.

In this study, the RA in the NOMA system is under the assumption of equal bandwidth allocation. According to [36], the performance of NOMA systems can be improved by performing unequal bandwidth and allocating higher power and wider bandwidth to users with stronger channel conditions. Therefore, we will consider unequal bandwidth allocation in our future work. We believe that our proposed CIM-based RA method can be easily extended to the NOMA system with unequal bandwidth allocation by assigning more than one resource block to one use, specifically, by considering extending the constraint in Eq. (11.c) to an arbitrary number n. Using these constraints, the energy function Eq. (19) can be computed. Finally, by recomputing $J_{ijkl}$ and $\lambda_{ij}$, more than one bandwidth can be allocated to one user using the CIM. Here, by adjusting the value of n, the upper limit of the bandwidth allocated to users with stronger channel conditions can be set. Thus, better system performance may be achieved with unequal bandwidth allocation compared to that with equal bandwidth allocation.

In this study, the channel allocation is optimized by the CIM, and then the power is assigned based on the optimizing methods in [5] using conventional computers. However, in large-scale NOMA systems, the time consumption for allocating power using conventional computers may reduce the speed of the optimization, which may bring a bottleneck to CIM-based optimization. Therefore, channel and power optimization will be jointly considered using the CIM in our future work.

In this study, we assumed that both of the BS and users have a single antenna, and the signal to each user is multiplexed in one channel at the same time. To improve the performance of the system, we will consider diversity schemes such as receiver diversity and antenna diversity in our future work. Specifically, the receiver diversity can improve the communication performance, such as a bit error rate, by combining the same data in the different channels to achieve a diversity gain [37]. Antenna diversity can be used to improve the outage performance of downlink NOMA systems combined with signal user multiple-input multiple-output [38].

Although quantum technologies as emerging technologies are expected to be deployed in 6G networks around 2030, it is still unclear when a practical quantum computer will be available. Hence, we will work with the communications and quantum computing industries in the future to try to achieve a practical implementation of the CIM in the next generation of wireless communication.