A method and accelerator for MCMC acceleration based on adaptive data subsampling

Abstract

The application discloses an MCMC acceleration method and accelerator based on adaptive data subsampling, and is applied to Bayesian modeling. i‑1 The acceleration method comprises the following steps: generating a candidate parameter sample, calculating the upper bound of the energy difference boundary between the candidate parameter sample θ' and the current parameter sample θ, generating the number B of small batch subsampling according to the upper bound of the energy difference boundary and the current hyperparameter χ, sampling according to the number B, and updating the small batch data set and the acceptance probability likelihood LH of θ' according to the acceptance probability of the data point, calculating the acceptance probability α by using the finally accumulated LH, accepting or rejecting θ' as the parameter sample of the current iteration of the model according to α, dynamically adjusting the hyperparameter χ every fixed iteration time, and repeating the above process to obtain the final parameter sample θNs of the Bayesian model. The application can realize the best balance between the sampling speed and the sample efficiency of the parameter in the model establishment process.

Description

TECHNICAL FIELD

[0001] The application belongs to the field of artificial intelligence and information processing, and particularly relates to an MCMC acceleration method and accelerator based on adaptive data subsampling. BACKGROUND

[0002] The standard Markov Chain Monte Carlo (MCMC) operation pseudocode is shown in Table 1. A starting point is first selected from the parameter space as the initial state. Then, in each loop iteration, a candidate state θ' is generated near the current state θ i-1 by using a Gaussian random walk (line 2). The acceptance probability a is calculated according to the Metropolis Hastings (MH) formula (line 3). The candidate state is accepted with probability a and used as the next state; otherwise, the current state is kept unchanged with probability 1-a (lines 4-8).

[0003]

[0004] The main workload of MCMC is concentrated in the calculation of the acceptance probability. For each proposed sample, the likelihood of the entire data set needs to be checked one by one when calculating the acceptance probability. Assuming that each data point is independent, the standard MH method requires the multiplication of the likelihoods p(x | θ) of all data points when calculating the posterior probability distribution n , that is:

[0005]

[0006] where p(θ) is the prior distribution, which represents the initial probability distribution given the parameter θ. p(x n | θ) refers to the likelihood of observing the data point x n given the parameter θ, that is, the probability of the occurrence of the data point x n . Therefore, the MCMC method requires more computing resources, especially when the MCMC algorithm encounters a large data set. Running MCMC on such a large data set is usually too computationally intensive to produce any practical use [1] . In recent years, a lot of effort has been devoted to accelerating the MCMC algorithm. These methods can be roughly divided into two categories: divide-and-conquer methods and data subsampling-based algorithms.

[0007] The divide-and-conquer method divides the initial data set into batches and then runs MCMC in parallel on different data partitions, and then combines the results to obtain an approximation of the posterior. This method has two problems: first, the MCMC method does not return the posterior as a mathematical function, so an additional step is needed to combine the batch posteriors according to some strategy[2] Another problem is that it is not possible to efficiently combine non-Gaussian batch posteriors. When the batch posteriors are not Gaussian, we need to consider how to effectively combine these non-Gaussian batch posteriors. [3] These issues undoubtedly significantly increase the algorithm's complexity and computational cost. Furthermore, processing and combining non-Gaussian batch posteriors can introduce approximation errors, especially when using methods such as Gaussian approximation. This can lead to inaccurate estimates of the true posterior distribution, affecting subsequent inference results. Consequently, the divide-and-conquer approach struggles to accurately perform Bayesian inference.

[0008] Data subsampling methods are popular for their efficiency and aim to reduce the number of data points that need to be processed in each iteration when calculating the Metropolis acceptance probability, thereby speeding up the calculation. This method does not use the entire dataset to calculate the posterior probability distribution, but instead uses subsamples or small batches of MH to calculate an approximation. [4] Some current data subsampling methods introduce approximate MH rules to control the bias by accepting or rejecting samples based on partial data. [5] Some effectively subsample the data based on the contribution of each likelihood term, which can be transformed into an unbiased estimator of the likelihood function through bias correction. [6] Others achieve the exact posterior distribution by replacing each product term in the likelihood probability with a simpler function. [7] .

[0009] Existing data subsampling methods either require a strict lower bound on the likelihood function or introduce high additional computational overhead. Moreover, existing data subsampling MCMC optimization methods do not consider designing specialized hardware devices to accelerate the MCMC algorithm, thus limiting their promotion in practical applications.

[0010] In summary, existing MCMC methods and acceleration methods based on data subsampling MCMC have the following main problems:

[0011] 1) MCMC methods are slow when applied to big data processing and cannot meet the requirements of real-time processing;

[0012] 2) Existing MCMC acceleration methods based on data subsampling do not support real-time adjustment of the sub-dataset size and cannot achieve the optimal balance between speed and sampling efficiency.

[0013] 3) Existing MCMC acceleration devices based on data subsampling need to read and write from off-chip memory when reading sub-data sets, which results in large communication overhead and does not maximize the use of on-chip storage resources.

[0014] 4) Existing MCMC acceleration devices based on data subsampling require a large amount of hardware resources due to the high algorithm complexity, which affects the hardware execution efficiency. Summary of the Invention

[0015] The present invention provides an MCMC acceleration method and accelerator based on adaptive data subsampling, which can achieve an optimal balance between parameter sampling speed and sample efficiency during the model establishment process.

[0016] In order to achieve the above technical objectives, the present invention adopts the following technical solutions:

[0017] An MCMC acceleration method based on adaptive data subsampling for Bayesian modeling, including:

[0018] Step 1: Given the initial parameter sample θ0 and initial hyperparameter χ0 of the Bayesian model;

[0019] Step 2: Generate candidate parameter samples θ' for the i-th iteration, and calculate the difference between the candidate parameter samples θ' and the parameter samples θ for the i-1-th iteration. i-1 The energy difference between them is bounded by the upper bound M(θ i-1 ,θ');

[0020] Step 3: According to the upper bound of energy difference M(θ i-1 ,θ') and the current hyperparameter χ, generate the number of mini-batch subsampling B;

[0021] Step 4: Sample the training data set according to the number B of mini-batch subsampling, and update the mini-batch data set according to the acceptance probability of each sampled data point, as well as the likelihood LH of the acceptance probability of the candidate parameter sample θ';

[0022] Step 5: Use the mini-batch data set to finally update the likelihood LH of the acceptance probability and calculate the acceptance probability α of the candidate parameter sample θ';

[0023] Step 6: If the acceptance probability α of the candidate parameter sample θ' is greater than the random number in the given range, then the candidate parameter sample θ' is accepted as the parameter sample θ of the Bayesian model for the i-th iteration. i , that is, θ i =θ'; otherwise, reject the candidate parameter sample θ' and replace the parameter sample θ i-1 As the parameter sample of the i-th iteration of the Bayesian model, θ i =θ i-1 ;

[0024] Step 7: For every given number of iterations, the hyperparameter X is dynamically adjusted based on the number of candidate parameter samples accepted within the interval.

[0025] Step 8: Let i = i + 1, return to step 2, and continue to generate the parameter sample θ for the i + 1th iteration. i+1 , until the maximum number of iterations is reached, and the final parameter sample θ of the Bayesian model is obtained Ns .

[0026] Furthermore, step 2 uses a Gaussian random walk with a step size of s to generate candidate parameter samples θ', which can be expressed as:

[0027] θ'~θ i-1 +Normal(0,s 2 I D )

[0028] Where, Normal(0,s 2 I D ) represents a normally distributed random variable with a mean vector of 0 and a covariance matrix of s 2 I D ; Among them, I D is a D×D identity matrix.

[0029] Furthermore, the upper bound of the energy difference M(θ i-1 ,θ') is calculated as: M(θ i-1 ,θ')=||θ'-θ||.

[0030] Furthermore, the number of mini-batch subsamplings B is sampled from a Poisson distribution and is expressed as:

[0031] B~Poisson(χC 2 M 2 (θ i-1 ,θ')+CM(θ i-1 ,θ'))

[0032] Where Poisson represents Poisson distribution, and C is the sum of the weights of the training data set.

[0033] Furthermore, step 4 specifically includes:

[0034] Step 4.1: Sample the training data set according to the weight of each training data, and the index n of the sampled data point b satisfy:

[0035] P(n b =n)=c n / C;

[0036] In the formula, P(n b =n) represents the probability that the nth training data is sampled, n is the index of each training data in the training data set, c n is the weight of the nth training data, and C is the sum of the weights of the training data set;

[0037] Step 4.2, calculate two intermediate variables φ1 and φ2:

[0038]

[0039] Where U(θ) is the component energy function of the Bayesian model, and Indicates index n b The parameters of the Bayesian model for the data points are θ' and θ i-1 The component energy, Indicates index n b The weight of the data point;

[0040] Step 4.3, calculate the acceptance probability of each data point:

[0041]

[0042] Where, Indicates index n b The acceptance probability of the sampled data points;

[0043] Step 4.4, determine the acceptance probability of each data point Is it greater than or equal to the random number u1, u1 ~ Uniform (0,1); if Then add the data point to the mini-batch dataset and update the likelihood LH of the acceptance probability of the candidate parameter sample θ';

[0044]

[0045] Furthermore, the acceptance flag accumulator is set, and its value is initialized to 0 in step 1; in the loop iteration process from step 2 to step 8:

[0046] If the candidate parameter sample θ' is accepted as the parameter sample θ of the Bayesian model in step 6 i , then the value of the acceptance flag accumulator is increased by 1;

[0047] Step 7 is specifically as follows: at each interval of a given number of iterations, determine whether the value of the acceptance flag accumulator is less than the ideal number of acceptances: if so, adjust the hyperparameter χ = (1 + δ)χ, and set the value of the acceptance flag accumulator to 0; otherwise, adjust the hyperparameter χ = (1 - δ)χ, and set the value of the acceptance flag accumulator to 0; where δ is the adjustment ratio.

[0048] A hardware accelerator for implementing any of the above-mentioned MCMC acceleration methods based on adaptive data subsampling on an FPGA, comprising: a sample proposal evaluator, a mini-batch sampler, a likelihood probability calculation module, a parameter sample generation module, a random number generator, a data register, and a hyperparameter adjustment module;

[0049] The sample proposal evaluator is used to generate a candidate parameter sample θ' for the i-th iteration, calculate the difference between the candidate parameter sample θ' and the parameter sample θ for the i-1-th iteration i-1 The energy difference between them is bounded by the upper bound M(θ i-1 ,θ'), and the calculated energy difference is bounded by the upper bound M(θ i-1 ,θ') is stored in the data register;

[0050] The sample proposal evaluator is further configured to: read the upper bound M(θ i-1 ,θ') and the current hyperparameter χ, and according to the read M(θ i-1 ,θ') and χ generate the number of mini-batch subsampling B;

[0051] The mini-batch sampler is used to: sample in the training data set according to the number B of mini-batch sub-sampling;

[0052] The likelihood probability calculation module is used to: update the small batch data set according to the acceptance probability of each sampled data point, and update the likelihood LH of the acceptance probability of the candidate parameter sample θ';

[0053] The parameter sample generation module is used to: calculate the acceptance probability α of the candidate parameter sample θ' using the likelihood LH of the final updated acceptance probability of the mini-batch data set; and determine whether to accept the candidate parameter sample θ' according to the acceptance probability: if the acceptance probability α of the candidate parameter sample θ' is greater than the random number u2 within a given range, then accept the candidate parameter sample θ' as the parameter sample θ of the i-th iteration of the Bayesian model i , that is, θ i =θ'; otherwise, reject the candidate parameter sample θ' and replace the parameter sample θ i-1 As the parameter sample of the i-th iteration of the Bayesian model, θ i =θ i-1 ;

[0054] The hyperparameter adjustment module is used to dynamically adjust the hyperparameter χ according to the number of candidate parameter samples accepted in each interval of a given number of iterations;

[0055] The random number generator is used to generate a random number u2 within a given range.

[0056] Further, the hardware accelerator further comprises an on-chip memory and an off-chip memory, the on-chip memory is used to save training data with higher weights and the weights thereof in the training data set, the off-chip memory is used to save training data with lower weights and the weights thereof in the training data set, and the process of allocating the memory to the two kinds of training data is completed during the preprocessing of the acceleration method.

[0057] Further, the hardware accelerator further comprises an on-chip memory and an off-chip memory, the on-chip memory is used to save training data with higher weights and the weights thereof in the training data set, the off-chip memory is used to save training data with lower weights and the weights thereof in the training data set, and the process of allocating the memory to the two kinds of training data is completed during the preprocessing of the acceleration method.

[0058] Further, the likelihood probability calculation module comprises two likelihood calculation sub-modules, which are respectively used for calculating the likelihood of the acceptance probability of the candidate parameter sample θ' and the parameter sample θ i-1 .

[0059] Compared with the prior art, the application has the beneficial technical effects that:

[0060] 1. The application uses an improved data subsampling method, improves the model parameter sampling speed of the MCMC method in processing a large data set, and can meet the requirements of real-time processing.

[0061] 2. By improving the existing data subsampling MCMC method, an adaptive data subsampling method is proposed, which can reduce the size of the sub-data set while ensuring the calculation accuracy of the acceptance rate of the candidate parameter sample in the model parameter generation process, so as to achieve the best balance between the sampling speed of the model parameter and the sample efficiency; further, the speed of model establishment can be improved, and the accuracy of the model established based on the generated parameter sample in performing respective tasks is ensured.

[0062] 3. The hardware accelerator of the adaptive data subsampling-based MCMC acceleration method is realized on the FPGA, and the overall performance problem of the memory-limited system is solved by designing an intelligent storage system.

[0063] 4. The hardware accelerator of the adaptive data subsampling-based MCMC acceleration method is realized on the FPGA, the data processing flow and resource allocation are optimized, the calculation efficiency is significantly improved, and the hardware resources are saved. BRIEF DESCRIPTION OF DRAWINGS

[0064] Figure 1 is the overall architecture of the hardware accelerator described in the embodiments of the application;

[0065] Figure 2 is the principle diagram of the memory subsystem and the small-batch sampler module in the embodiments of the application;

[0066] Figure 3 is the principle diagram of the likelihood probability calculation module in the embodiments of the application;

[0067] Figure 4 This is a risk comparison chart of the present invention and other methods in an example;

[0068] Figure 5 This is a comparison chart of the test accuracy of the present invention and other methods on examples.

[0069] Figure 6 This is a comparison chart of probability distribution of the present invention and other methods in examples. DETAILED DESCRIPTION

[0070] The following is a detailed description of an embodiment of the present invention. This embodiment is based on the technical solution of the present invention, provides a detailed implementation method and a specific operation process, and further explains the technical solution of the present invention.

[0071] The embodiment of the present invention proposes an MCMC acceleration method and hardware accelerator based on adaptive data subsampling, which is applied to Bayesian modeling. The main idea is to improve the scalability and efficiency of the MH algorithm on large-scale datasets by using dynamically adjustable minibatch data while maintaining the accuracy of the MH algorithm. The hardware implementation is based on full consideration of the architecture of modern hardware. The following problems can be solved: 1) Traditional MCMC methods have high time and computational costs when processing large datasets; 2) Although other data subsampling MCMC methods can improve efficiency by approximating the acceptance rate, this approximation may introduce bias, resulting in inaccurate sampling results; 3) Under traditional design methods, the batch size is usually fixed, and they may not be able to fully utilize the local characteristics of the data to optimize the computational amount of each iteration; 4) Under traditional design methods, cross-memory access of large data leads to increased memory access latency, thereby limiting the overall performance of the system.

[0072] Example 1

[0073] This embodiment provides an MCMC acceleration method based on adaptive data subsampling, namely Self-ShrinkingMCMC, referred to as SS-MCMC acceleration method in this paper, which is applied to Bayesian modeling. and the Bayesian parameterized model θ, the goal is to calculate the posterior distribution With the idea of ​​factorization, the formula can be expressed as:

[0074]

[0075] Referring to the SS-MCMC pseudocode shown in Table 2, the MCMC acceleration method based on adaptive data subsampling in this embodiment includes the following steps:

[0076] Step 1: Given the initial parameter sample θ0 and initial hyperparameter χ0 of the Bayesian model.

[0077] In the specific implementation process, the input includes not only the initial parameter sample θ0 and the initial hyperparameter χ0, but also the number of parameter samples Ns, the amount of training data N, the local limit parameter c n , C, M, ideal acceptance number target accept, adjustment interval adjust interval and adjustment ratio δ.

[0078] In addition, this embodiment sets an acceptance flag accumulator current accept, whose value is initialized to 0 in step 1, that is, current accept=0.

[0079] According to the above input data, the parameter sample θ corresponding to the i-th iteration is generated by repeating steps 2 to 8 in a loop, with each loop being one iteration. i .

[0080] Step 2: Generate candidate parameter samples θ' for the i-th iteration, and calculate the difference between the candidate parameter samples θ' and the parameter samples θ for the i-1-th iteration. i-1 The energy difference between them is bounded by the upper bound M(θ i-1 ,θ').

[0081] In this embodiment, a Gaussian random walk with a step size of s is used to generate candidate parameter samples θ', which can be expressed as:

[0082] θ'~θ i-1 +Normal(0,s 2 I D )

[0083] Where, Normal(0,s 2 I D ): This is a multivariate normal distribution, which has a mean vector of 0 (i.e., the mean of each dimension is 0) and a covariance matrix of s 2 I D , here I D is a D×D identity matrix, s 2 is a scalar that represents the variance in each dimension. Therefore, s 2 I D It is a diagonal matrix, and the elements on the diagonal are all s 2 , indicating that the variance in each dimension is s 2 , and the covariance between different dimensions is 0, that is, the dimensions are independent. θ'~θ i-1 +Normal(0,s 2 I D) This expression means that the parameter θ' is obtained by i-1 Add a mean of 0 and a variance of s to the 2 The update is done by independent normal random variables.

[0084] The calculation formula of the upper bound of the energy difference in this embodiment is: M(θ i-1 ,θ')=||θ'-θ||.

[0085] Step 3: According to the upper bound of energy difference M(θ i-1 ,θ') and the current hyperparameter χ, generate the number B of mini-batch subsampling.

[0086] In this embodiment, the number B of small batch sub-sampling is sampled from the Poisson distribution, which is expressed as:

[0087] B~Poisson(χC 2 M 2 (θ i-1 ,θ')+CM(θ i-1 ,θ'))

[0088] Where Poisson represents Poisson distribution, and C is the sum of the weights of the training data set.

[0089] Step 4: Sample the training dataset according to the number B of mini-batch subsampling, and update the mini-batch dataset based on the acceptance probability of each sampled data point, as well as the likelihood LH of the acceptance probability of the candidate parameter sample θ'. Specifically, it includes:

[0090] Step 4.1: Sample the training data set according to the weight of each training data, and the index n of the sampled data point b satisfy:

[0091] P(n b =n)=c n / C;

[0092] In the formula, P(n b =n) represents the probability that the nth training data is sampled, n is the index of each training data in the training data set, c n is the weight of the nth training data, and C is the sum of the weights of the training data set;

[0093] Step 4.2, calculate two intermediate variables φ1 and φ2:

[0094]

[0095] Where U(θ) is the component energy function of the Bayesian model, and Indicates index n bThe parameters of the Bayesian model for the data points are θ' and θ i-1 The component energy, Indicates index n b The weight of the data point;

[0096] Step 4.3, calculate the acceptance probability of each data point:

[0097]

[0098] Where, Indicates index n b The acceptance probability of the sampled data points;

[0099] Step 4.4, determine the acceptance probability of each data point Is it greater than or equal to the random number u1, u1 ~ Uniform (0,1); if Then add the data point to the mini-batch dataset and update the likelihood LH of the acceptance probability of the candidate parameter sample θ';

[0100]

[0101] The likelihood LH of the acceptance probability is set to its initial value of 0 before sampling in step 4.

[0102] Step 5: Use the mini-batch data set to finally update the likelihood LH of the acceptance probability and calculate the acceptance probability α = exp(LH) of the candidate parameter sample θ'.

[0103] Step 6: If the acceptance probability α of the candidate parameter sample θ' is greater than the random number u2 within the given range, then the candidate parameter sample θ' is accepted as the parameter sample θ of the Bayesian model for the i-th iteration. i , that is, θ i =θ', and add 1 to the value of the acceptance flag accumulator current accept; otherwise, reject the candidate parameter sample θ' and replace the parameter sample θ i-1 As the parameter sample of the i-th iteration of the Bayesian model, θ i =θ i-1 .

[0104] The random number u2 in this embodiment is given in the range of (0, 1), that is, u2~Uniform(0, 1) is generated by the random number generator.

[0105] In step 7, every interval with a given number of iterations, the hyperparameter χ is dynamically adjusted according to the number of candidate parameter samples accepted in the interval.

[0106] In this embodiment, according to the set adjustment interval, if the current number of iterations is an integer multiple of the adjustment interval, it is determined whether the value of the current acceptance flag accumulator is less than the ideal number of acceptances:

[0107] If so, adjust the hyperparameter χ = (1 + δ)χ and set the value of the acceptance flag accumulator to 0;

[0108] Otherwise, adjust the hyperparameter χ = (1-δ)χ ​​and set the value of the acceptance flag accumulator to 0.

[0109] Step 8: Let i = i + 1, return to step 2, and continue to generate the parameter sample θ for the i + 1th iteration. i+1 , until the maximum number of iterations is reached.

[0110] Step 9: Output the parameter sample θ of the Bayesian model i , i=1,…,Ns.

[0111] In this embodiment, is the component energy function, n represents the index of the data point in the training dataset D, and θ represents the parameters of the Bayesian model. This function measures the "energy" or "cost" of the observed data point under the given parameters θ, which can be used to measure the fit between the parameters and the data. U n The smaller the value of (θ), the closer the parameter θ is to the observed data x. n In the MCMC acceleration method, the present invention focuses on calculating the acceptance rate α used to decide whether to update the parameters. At this time, the acceptance rate α can be calculated using the following formula:

[0112]

[0113] In order to allow the use of small batches of data in each iteration to approximate the posterior distribution of the entire dataset without significantly changing the target distribution, the present invention limits the energy difference between data points. Specifically, this embodiment uses positive constants c1,…,c N and C, which satisfies and symmetric functions For any θ,θ′∈Θ, the energy difference is bounded to be |U n (θ)-U n (θ')|≤c n ·M(θ,θ'). By this restriction, the energy function difference ΔU n Can be locally c n and M(θ,θ'), the algorithm can only focus on a small batch of data points, thus greatly reducing the amount of computation.

[0114] In this embodiment, it is assumed that each energy function U n It's Ln -Lipschitz continuous, set c n =L n And M(θ,θ')=‖θ'-θ||. Among them, L n Is a positive constant representing the energy function U n The Lipschitz constant, which limits the maximum rate of change of the function value, and ||θ'-θ|| represents the L2 norm between θ' and θ. The batch size B is sampled from the Poisson distribution (line 5 of the pseudo code shown in Table 2). Then, according to c n The data points are sampled with replacement (line 8 of the pseudo code shown in Table 2), c n can be considered as the weight of the data point, c n The larger the value of is, the more likely it is to be included in the mini-batch of size B, and each sampled point is likely to be deleted from the mini-batch (line 12 of the pseudo code shown in Table 2).

[0115] When calculating the probability of whether the selected data point needs to be eliminated (line 10 of the pseudocode shown in Table 2), the probability in line 9 of the pseudocode shown in Table 2 is updated through a single calculation, so that the SS-MCMC method of the present invention can automatically update the likelihood of the MH acceptance probability after the data point is selected, thereby greatly reducing the complexity of the calculation. However, compared with the standard MCMC method that only calculates the probability of the proposal and keeps the previous proposal unchanged, due to the change of the subset, the SS-MCMC method of the present invention needs to recalculate the probability of the proposal and the current parameters in each iteration. Although this introduces additional calculations, the use of subsets reduces the overall computational cost.

[0116] The hyperparameter χ can be used to balance the algorithm's batch size and convergence rate. Larger values ​​of χ may increase the batch size, improving the convergence rate, but may also result in higher computational cost. Smaller values ​​of χ may reduce the computational effort per iteration, reducing the convergence rate. One of the key innovations of the SS-MCMC acceleration method of this invention is the batch size auto-shrinking mechanism, as shown in lines 20 to 24 of the pseudocode shown in Table 2. This mechanism dynamically adjusts the required subset size by adjusting the hyperparameter based on the sample acceptance rate, ensuring an optimal balance between sample efficiency and computational cost.

[0117]

[0118]

[0119] Example 2

[0120] This embodiment provides a hardware accelerator for implementing the MCMC acceleration method based on adaptive data subsampling described in Example 1 on an FPGA. Figure 1As shown, it includes: sample proposal evaluator, mini-batch sampler, likelihood probability calculation module, parameter sample generation module, random number generator, data register, hyperparameter adjustment module, acceptance flag accumulator, etc.

[0121] The sample proposal evaluator is used to generate a candidate parameter sample θ' for the i-th iteration, calculate the difference between the candidate parameter sample θ' and the parameter sample θ for the i-1-th iteration i-1 The energy difference between them is bounded by the upper bound M(θ i-1 ,θ'), and the calculated energy difference is bounded by the upper bound M(θ i-1 ,θ') is stored in the data register for subsequent operations.

[0122] The sample proposal evaluator is further configured to: read the upper bound M(θ i-1 ,θ') and the current hyperparameter χ, and according to the read M(θ i-1 ,θ') and χ generate the number B of mini-batch subsampling.

[0123] The mini-batch sampler samples the training dataset according to the number B of mini-batch sub-sampling.

[0124] In big data scenarios, observed data sets are usually initially stored in off-chip DDR memory and then processed by the MCMC accelerator. In a more preferred embodiment of the present invention, in order to reduce the communication overhead caused by cross-memory access, a data prefetching technique that can make optimal use of on-chip memory is introduced. As shown in Table 2, the SS-MCMC pseudo code is used according to the corresponding c n To optimize this process, this embodiment stores data points with high inclusion probabilities (data with higher weights have higher inclusion probabilities) in the faster on-chip BRAM, while data points with lower probabilities are stored in the slower off-chip memory. This allocation is determined during preprocessing before the SS-MCMC accelerator starts the sample proposal evaluator, effectively meeting the needs of real-time data transmission. The customized memory subsystem and mini-batch sampler module are as follows: Figure 2 shown.

[0125] During the preprocessing process, the weights of each point in the data set are calculated, and the data points with high weights are stored in the on-chip BRAM, and their corresponding weights are also stored in the BRAM of the adjacent memory address. Specifically, the address of the data point with the highest weight in the BRAM is first determined, and used as the head address of the on-chip BRAM to quickly access key data. At the same time, a starting address is set for the off-chip memory as its head address to facilitate data reading. In addition, the address of the last data point in the BRAM, that is, the tail address of the BRAM, needs to be determined to manage the entire storage area. In order to improve efficiency, an address mapping table can also be maintained to map data points to their actual locations in the BRAM or off-chip memory, and this mapping table can be stored in the cache to reduce search time. Through these carefully designed strategies, the acceleration method of the present invention can significantly improve the data processing speed and meet the needs of real-time data processing.

[0126] The likelihood probability calculation module is used to update the small batch data set according to the acceptance probability of each sampled data point, and update the likelihood LH of the acceptance probability of the candidate parameter sample θ'.

[0127] The likelihood probability calculation of the embodiment of the present invention has the following architecture: Figure 3 This module uses a pipelined design to efficiently calculate the likelihood of each data point. The likelihood calculation module completes the data removal probability and the likelihood calculation in line 9 of pseudocode 2. If the data is removed, the likelihood is filled with 0. If the data is selected, its likelihood is calculated and output.

[0128] The algorithm's i-1th iteration generates parameter samples θ i-1 In the process of , the training data contained in the small batch data set used is not guaranteed to be the same as the training data in the small batch data set used in the i-th iteration (the training data used in each iteration is randomly selected according to the data weight and is random). Therefore, in the i-th iteration, the parameter sample θ must be recalculated. i-1 The likelihood of the training data used in the i-th iteration. Given that this stage represents the main computational load of the SS-MCMC acceleration method of the present invention, the embodiment of the present invention uses parallel computing to process multiple likelihood calculation submodules simultaneously. The likelihood probability calculation of this embodiment has two likelihood calculation submodules, namely Figure 3 Likelihood module 1 and likelihood module 2 in the , respectively calculate the candidate parameter sample θ 'and the current parameter θ in each iteration i-1 The likelihood is based on the mini-batch used in iteration i. This is required by SS-MCMC.

[0129] The parameter sample generation module is configured to: update the likelihood LH of the obtained acceptance probability by using the mini-batch dataset, calculate the acceptance probability a of the candidate parameter sample θ'; and determine whether to accept the candidate parameter sample θ' according to the acceptance probability: if the acceptance probability a of the candidate parameter sample θ' is greater than a random number u2 in a given range, accept the candidate parameter sample θ' as the parameter sample θ of the i th iteration of the Bayesian model, that is, θ = θ'; otherwise, reject the candidate parameter sample θ', and set the parameter sample θ as the parameter sample θ of the i th iteration of the Bayesian model, that is, θ = θ i , that is, θ i = θ'. i-1 i i-1 .

[0130] The acceptance flag accumulator is configured to count the number of times of accepting the proposed parameter sample within a given number of iteration times per interval.

[0131] The hyperparameter adjustment module is configured to: dynamically adjust the hyperparameter χ according to the number of times of accepting the candidate parameter sample within the interval, every given number of iteration times.

[0132] The random number generator is configured to: generate a random number u2 in a given range.

[0133] The hardware accelerator is implemented on the FPGA to realize the MCMC acceleration method based on adaptive data subsampling in Embodiment 1, and the specific implementation steps are as follows:

[0134] S1: In the preprocessing stage, the weights of the observable data stored in the off-chip DDR are calculated one by one, and the data and the corresponding weights are respectively stored in the on-chip BRAM and the off-chip DDR according to the high and low of the weights, and an address mapping table is sequentially formed and stored in the cache, waiting for subsequent operations.

[0135] S2: The sample proposal evaluator reads the value of the current sample, generates a candidate sample on this basis by random walk, and calculates the energy difference between the two and stores it in the parameter register and generates a mini-batch size B, waiting for the next operation.

[0136] ​​S3: After S1 is executed, the small batch sampler starts sampling, and does not need to rely on the execution result of S2. When the small batch sampler is initially run, the weight of the first data is read from the head address of the BRAM to the accumulation counter, and the calculated accumulated weight is stored in the adjacent position of the data; then the weight of the next data is input to the accumulation counter, and the result is also stored in the adjacent position of the data; until the accumulation count reaches the last data. After the small batch sampler is initialized, the random number generator generates a random number u between 0 and 1 multiplied by the local limit parameter C, and then sequentially compares it with the accumulated weight of the data until the value of the accumulated weight is not less than the product value, and finally outputs the data and weight value corresponding to the accumulated weight to the FIFO in the chip. The process is continuously executed to ensure that the data flow meets the design requirements.

[0137] S4: After S2 is executed, the likelihood evaluation module reads the parameter register and the data of the FIFO, completes a likelihood probability calculation and inputs it to the accumulator.

[0138] S5: When the number of times of S4 execution exceeds the size of B, the flag of S4 waiting is triggered, and the value of the likelihood probability accumulator is read out and compared with the random number to determine whether to accept the candidate sample. After the determination is completed, the flag of S2 execution is triggered.

[0139] Example verification:

[0140] The performance of the SS-MCMC accelerator of the present application is evaluated using logistic regression on the MNIST dataset of handwritten digits. And compared with several different MCMC methods (other methods are also deployed on FPGA devices). The remaining several methods are:

[0141] MH: represents the traditional Metropolis-Hastings MCMC method.

[0142] flyMCMC: a commonly used subsampling MCMC method used for comparison, proposed in the paper "Firefly Monte Carlo: Exact MCMC with Subsets of Data".

[0143] TunaMH: the latest subsampling MCMC method, proposed in the paper "Asymptotically Optimal Exact Minibatch Metropolis-Hastings".

[0144] First, the standard MH algorithm is run to determine the true prediction mean, and then the risk is calculated as the mean square error (MSE) of the estimated values generated by the four evaluation algorithms. Figure 4The results show the risk of estimating the first dimension of the generated parameter samples on the MNIST dataset. The results show that the proposed SS-MCMC achieves the fastest risk reduction, reducing the risk faster within the same execution time. This rapid convergence is mainly attributed to the computational acceleration brought by its dynamic batch adjustment and optimized hardware design, which has a significant advantage over MH and flyMCMC methods in reducing risk. In addition, due to its adaptive hyperparameter adjustment, SS-MCMC outperforms TunaMH in reducing risk within the same time frame, which can effectively optimize the batch size.

[0145] The sample parameters generated by the four different methods are used as parameters in the logistic regression model, and reasoning is performed on the MNIST test dataset to compare the classification accuracy of different methods to evaluate the test accuracy. Figure 5 The test accuracy of our method is compared to other methods on examples. The line widths in the figures represent the standard deviation observed across these runs (each algorithm was run for 100 iterations). The results confirm that our SS-MCMC method consistently outperforms both MH and flyMCMC in terms of test accuracy. While our SS-MCMC and TunaMH achieve similarly high accuracy, our SS-MCMC achieves this level much faster, highlighting its efficiency in producing accurate results.

[0146] Figure 6 The performance of the proposed method, the conventional MH method, and two other different methods in completing the classification task of handwritten digits 7 and 9 using logistic regression on the MNIST dataset is given. The parameter samples generated by each method are multidimensional, and their performance is usually evaluated by observing the probability distribution of these algorithms in a certain dimension. Specifically, Figure 6 Focusing on the first-dimensional probability distributions of the parameter samples generated by these methods, experimental results show that the mean of the proposed SS-MCMC acceleration method is closest to that of the conventional MH method, experimentally demonstrating the unbiased nature of the proposed SS-MCMC acceleration method. This finding further validates the effectiveness and accuracy of the SS-MCMC acceleration method.

[0147] The above embodiments are preferred embodiments of the present application. Ordinary technicians in this field can also make various changes or improvements on this basis. Without departing from the overall concept of the present application, these changes or improvements should fall within the scope of protection required by the present application.

[0148] References:

[0149] [1]Korattikara A,Chen Y,Welling M.Austerity in MCMC land:Cutting theMetropolis-Hastings budget[C] / / International conference on machinelearning.PMLR,2014:181-189.

[0150] [2]Huang Z,Gelman A.Sampling for Bayesian computation with largedatasets[J].Available at SSRN 1010107,2005.

[0151] [3]Neiswanger W,Wang C,Xing E.Asymptotically exact,embarrassinglyparallel MCMC[J].arXivpreprint arXiv:1311.4780,2013.

[0152] [4]Bardenet R,Doucet A,Holmes C.On Markov chain Monte Carlo methodsfor tall data[J].Journal ofMachine Learning Research,2017,18(47).

[0153] [5]Rémi Bardenet,Arnaud Doucet,and Chris Holmes.Towards scaling upMarkov chain Monte Carlo:an adaptive subsampling approach.In InternationalConference on Machine Learning,2014.

[0154] [6]Quiroz M,Kohn R,Villani M,et al.Speeding up MCMC by efficient datasubsampling[J].Journal ofthe American Statistical Association,2018.

[0155] [7]Maclaurin D,Adams R P.Firefly Monte Carlo:Exact MCMC with subsetsofdata[J].arXiv preprint arXiv:1403.5693,2014.

Claims

1. A method for accelerating MCMC based on adaptive data subsampling on FPGA, characterized in that: Applications to Bayesian modeling, including: Step 1: Given the initial parameter sample of the Bayesian model and initial hyperparameters ; Step 2: Generate Candidate parameter samples for iterations , calculate the candidate parameter sample With the Parameter samples for iterations The energy difference between ; Step 3: According to the upper bound of energy difference and the current hyperparameters , the number of mini-batch subsamplings generated ; Step 4: Subsample by the number of mini-batches Sampling in the training dataset and updating the mini-batch dataset and candidate parameter samples based on the acceptance probability of each sampled data point The likelihood of acceptance probability ; Step 5: Finally update the likelihood of the obtained acceptance probability using the mini-batch dataset , calculate the candidate parameter sample The acceptance probability ; Step 6: If the candidate parameter sample The acceptance probability If the random number is greater than the given range, the candidate parameter sample is accepted As a Bayesian model Parameter samples for iterations ,Right now ; Otherwise reject the candidate parameter sample , the parameter sample As a Bayesian model The parameter sample of the iteration, that is, ; Step 7: For every given number of iterations, the hyperparameters are dynamically adjusted based on the number of candidate parameter samples received within the interval. ; Step 8, let , return to step 2 and continue to generate Parameter samples for iterations , until the maximum number of iterations is reached, and the final parameter sample of the Bayesian model is obtained ; S1: In the preprocessing stage, the weights of observable data stored in the off-chip DDR are calculated one by one, and the data and the corresponding weights are stored in the on-chip BRAM and the off-chip DDR respectively according to the weight. An address mapping table is sequentially formed and stored in the cache, waiting for subsequent operations; S2: The sample proposal evaluator reads the value of the current sample, and based on this, randomly walks to generate candidate samples, calculates the energy difference between the two, stores it in the parameter register, and generates a minibatch size B, waiting for the next operation; where minibatch stands for small batch; S3: After S1 is executed, the small batch sampler starts sampling without relying on the execution result of S2; When the mini-batch sampler runs for the first time, it reads the weight from the BRAM header address into the cumulative counter and stores the calculated cumulative weight in the adjacent position of the data; then it takes the weight of the next data and inputs it into the cumulative counter, and also stores the result in the adjacent position of the data; until the cumulative count reaches the last data; After the mini-batch sampler is initialized, the random number generator generates a random number u between 0 and 1, multiplies it by the local limit parameter C, and then compares it with the accumulated weight of the data one by one until the value of the accumulated weight is not less than the product value. Finally, the data and weight value corresponding to the accumulated weight are output to the FIFO on the chip. S4: After S2 is executed, the likelihood evaluation module reads the data from the parameter register and FIFO, completes a likelihood probability calculation and inputs it into the accumulator; S5: When the number of times S4 is executed exceeds the size of B, the S4 waiting flag is triggered, and the value of the likelihood probability accumulator is read out and compared with the random number to determine whether to accept the candidate sample. After the judgment is completed, the S2 execution flag is triggered.

2. The MCMC acceleration method based on adaptive data subsampling implemented on FPGA according to claim 1, characterized in that: Step 2 uses a step size of Generate candidate parameter samples using Gaussian random walk , expressed as: ； Where, represents a normally distributed random variable with a mean vector of 0 and a covariance matrix of ;in, is a The identity matrix of .

3. The MCMC acceleration method based on adaptive data subsampling implemented on FPGA according to claim 1, characterized in that: Energy difference upper bound The calculation formula is: .

4. The MCMC acceleration method based on adaptive data subsampling implemented on FPGA according to claim 1, characterized in that: The number of mini-batch subsamples B is sampled from a Poisson distribution and is expressed as: ； Where, represents the Poisson distribution, is the sum of the weights of the training dataset.

5. The MCMC acceleration method based on adaptive data subsampling implemented on FPGA according to claim 1, characterized in that: Step 4 specifically includes: Step 4.1: Sample the training data set according to the weight of each training data point, and the index of the sampled data point satisfy: ； Where, Indicates the The probability that a training data set is sampled is the index of each training data in the training dataset, For the The weight of the training data, is the sum of the weights of the training data set; Step 4.2, calculate two intermediate variables 、 : ； ； Where, is the component energy function of the Bayesian model, and Indicates that the index is The parameters of the Bayesian model for the data points are and The component energy, Indicates that the index is The weight of the data point; Step 4.3, calculate the acceptance probability of each data point: ； Where, Indicates that the index is The acceptance probability of the sampled data points; Step 4.4, determine the acceptance probability of each data point Is it greater than or equal to the random number , ;like , then add the data point to the mini-batch dataset and update the candidate parameter sample The likelihood of acceptance probability ; 。 6. The MCMC acceleration method based on adaptive data subsampling implemented on FPGA according to claim 1, characterized in that: Set the acceptance flag accumulator, whose value is initialized to 0 in step 1; in the loop iteration process from step 2 to step 8: If the candidate parameter sample is accepted in step 6 As a Bayesian model Parameter samples for iterations , then the value of the acceptance flag accumulator is increased by 1; Step 7 is as follows: At each given number of iterations, determine whether the value of the acceptance flag accumulator is less than the ideal acceptance number: If so, adjust the hyperparameters , and set the value of the acceptance flag accumulator to 0; otherwise, adjust the hyperparameters , and set the value of the acceptance flag accumulator to 0; To adjust the proportions.

7. A hardware accelerator for implementing the MCMC acceleration method based on adaptive data subsampling according to any one of claims 1 to 6 on an FPGA, characterized in that: include: Sample proposal evaluator, mini-batch sampler, likelihood probability calculation module, parameter sample generation module, random number generator, data register, hyperparameter adjustment module; The sample proposal evaluator is used to generate Candidate parameter samples for iterations , calculate the candidate parameter sample With the Parameter samples for iterations The energy difference between , and the calculated energy difference is bounded by Store to data register; The sample proposal evaluator is further configured to: read an upper bound of the energy difference from a data register and the current hyperparameters , and according to the read and The number of subsamples to generate a mini-batch ; The mini-batch sampler is used to: Sampling in the training dataset; The likelihood probability calculation module is used to update the small batch data set according to the acceptance probability of each sampled data point and update the candidate parameter samples The likelihood of acceptance probability ; The parameter sample generation module is used to finally update the likelihood of the acceptance probability using the mini-batch data set , calculate the candidate parameter sample The acceptance probability ; and determine whether to accept the candidate parameter sample based on the acceptance probability :If the candidate parameter sample The acceptance probability A random number greater than the given range , then accept the candidate parameter sample As a Bayesian model Parameter samples for iterations ,Right now ; Otherwise reject the candidate parameter sample , the parameter sample As a Bayesian model The parameter sample of the iteration, that is, ; The hyperparameter adjustment module is used to dynamically adjust the hyperparameters according to the number of candidate parameter samples received in each interval of a given number of iterations. ; The random number generator is used to generate random numbers within a given range. .

8. The hardware accelerator according to claim 7, characterized in that: It also includes on-chip memory and off-chip memory, the on-chip memory is used to store training data with higher weights and their weights in the training data set, and the off-chip memory is used to store training data with lower weights and their weights in the training data set, and the process of allocating memory to the two types of training data is completed during the preprocessing of the acceleration method.

9. The hardware accelerator according to claim 7, characterized in that: The invention also includes an address mapping table stored in the cache, for mapping each training data to their actual location in the on-chip memory or the off-chip memory.

10. The hardware accelerator according to claim 7, wherein: The likelihood probability calculation module includes two likelihood calculation submodules, which are used to calculate candidate parameter samples. and parameter samples Likelihood of acceptance probability.

Images