MCMC acceleration method and accelerator based on adaptive data subsampling

Abstract

The invention discloses an MCMC acceleration method and accelerator based on adaptive data sub-sampling, which are applied to Bayesian modeling, and the acceleration method comprises the following steps: generating a candidate parameter sample, and calculating the upper bound of the energy difference between the candidate parameter sample theta'and the current parameter sample theta i-1; generating the number B of small-batch quantum sampling according to the upper bound of the energy difference bound and the current hyper-parameter x; sampling according to the number B, and updating the acceptance probability likelihood LH of the small-batch data set and the theta'according to the acceptance probability of the data points; calculating an acceptance probability alpha by using the finally accumulated LH; accepting or rejecting theta'as a parameter sample of the current iteration of the model according to alpha; dynamically adjusting the hyper-parameter x every fixed iteration times; and repeating the process to obtain a final parameter sample theta Ns of the Bayesian model. According to the method, the optimal balance between the sampling speed and the sample efficiency of the parameters in the model building process can be realized.

Description

Technical Field

[0001] The present invention belongs to the field of artificial intelligence and information processing, and in particular relates to an MCMC acceleration method and an accelerator based on adaptive data subsampling. Background Art

[0002] The pseudo code of the standard Markov chain Monte Carlo (MCMC) operation is shown in Table 1. First, a starting point is selected from the parameter space as the initial state. Then, in each loop iteration, the first step is to use a Gaussian random walk (line 2) to find the current state θ. i-1 Generate a candidate state θ' nearby. The second step is to calculate the acceptance probability α according to the Metropolis Hastings (MH) formula (line 3). The third step is to accept the candidate state with probability α and use it as the next state; otherwise, keep the current state unchanged with probability 1-α (lines 4-8).

[0003]

[0004] The main workload of MCMC is concentrated in calculating the acceptance probability. For each proposed sample, the possibility 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 calculates the posterior probability distribution in each iteration. When the probability of all data points p(x n |θ) multiplied, that is:

[0005]

[0006] where p(θ) is the prior distribution, which represents the initial probability distribution under given parameters θ. n │θ) refers to the observation of data point x when the parameter θ is known n The likelihood of the data point x n Therefore, MCMC methods require 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 be of any practical use. [1] In recent years, a lot of effort has been invested in accelerating MCMC algorithms. These methods can be roughly divided into two categories: divide-and-conquer methods and algorithms based on data subsampling.

[0007] The divide-and-conquer approach divides the initial dataset into batches, runs MCMC in parallel on different data partitions, and then combines these results to obtain an approximation of the posterior. There are two problems with this approach: First, MCMC methods do not return the posterior as a mathematical function, so an additional step is required to combine the batches of 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 problems undoubtedly greatly increase the complexity and computational cost of the algorithm. At the same time, approximation errors may be introduced when processing and combining non-Gaussian batch posteriors, especially when using methods such as Gaussian approximation. This may lead to inaccurate estimates of the true posterior distribution, affecting subsequent inference results. As a result, the divide-and-conquer approach is difficult to accurately complete the work of Bayesian reasoning.

[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 an approximate MH rule to control the bias by accepting or rejecting samples based on part of the data. [5] Some effectively subsample the data based on the contribution of each likelihood term, and can transform it 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 effort. In addition, 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, the existing MCMC methods and acceleration methods based on data subsampling MCMC have the following main problems:

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

[0012] 2) The 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, so there is a large communication overhead and the on-chip storage resources are not utilized to the maximum extent.

[0014] 4) The 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 an accelerator based on adaptive data subsampling, which can achieve an optimal balance between the sampling speed and sample efficiency of parameters in the model building 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 sample θ' for the i-th iteration, and calculate 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 ,θ');

[0020] Step 3: According to the upper bound of energy difference M(θ i-1 ,θ') and the current hyperparameter χ, generate the number of mini-batch sub-sampling 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, and update 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 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, that is, θ i =θ i-1 ;

[0024] Step 7: For every given number of iterations, the hyperparameter X is dynamically adjusted according to the number of candidate parameter samples accepted in the interval;

[0025] Step 8, let i = i + 1, return to step 2, and continue to generate the parameter sample θ of 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 θ', expressed as:

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

[0028] In the formula, Normal(0,s 2 I D ) represents a normally distributed random variable, and the mean vector of the normal distribution is 0, and the covariance matrix is ​​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 sub-samplings B is sampled from a Poisson distribution, 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 point, 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 of the nth training data being 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 is n b The parameters of the Bayesian model for the data points are θ' and θ i-1 The component energy at Indicates index is n b The weight of the data point;

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

[0041]

[0042] In the formula, Indicates index is 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 data set and update the likelihood LH of the acceptance probability of the candidate parameter sample θ';

[0044]

[0045] Further, 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 i-th iteration 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 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 small 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 θ' of the i-th iteration, calculate the candidate parameter sample θ' and the parameter sample θ of 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 used to: read the energy difference upper bound M(θ i-1 ,θ') and the current hyperparameter χ, and according to the read M(θ i-1 ,θ') and χ generate the number of mini-batch sub-sampling 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 acceptance probability finally updated by the small 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, that is, θ i =θ i-1 ;

[0054] The hyperparameter adjustment module is used to dynamically adjust the hyperparameter x according to the number of times candidate parameter samples are 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] Furthermore, the hardware accelerator also includes an on-chip memory and an 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.

[0057] Furthermore, the method further includes an address mapping table stored in the cache, which is used to map the training data to their actual locations in the on-chip memory or the off-chip memory.

[0058] Furthermore, the likelihood probability calculation module includes two likelihood calculation submodules, which are used to calculate the candidate parameter sample θ' and the parameter sample θ i-1 The likelihood of acceptance probability.

[0059] Compared with the prior art, the beneficial technical effects of the present invention are:

[0060] 1. The present invention uses an improved data subsampling method to improve the model parameter sampling speed of the MCMC method when processing large data sets, which can meet the requirements of real-time processing;

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

[0062] 3. The present invention implements a hardware accelerator of the MCMC acceleration method based on adaptive data subsampling on FPGA, and solves the overall performance problem of the memory-limited system by designing an intelligent storage system;

[0063] 4. The present invention implements a hardware accelerator of the MCMC acceleration method based on adaptive data subsampling on FPGA, optimizes the data processing flow and resource allocation, achieves a significant improvement in computing efficiency, and saves hardware resources. BRIEF DESCRIPTION OF THE DRAWINGS

[0064] Figure 1 It is the overall architecture of the hardware accelerator described in the embodiment of the present invention;

[0065] Figure 2 is a schematic diagram of a memory subsystem and a small batch sampler module of an embodiment of the present invention;

[0066] Figure 3 is a schematic diagram of a likelihood probability calculation module according to an embodiment of the present invention;

[0067] Figure 4 It is a risk comparison chart of the present invention and other methods in examples;

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

[0069] Figure 6 It is a probability distribution comparison diagram of the present invention and other methods on 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, and provides a detailed implementation method and a specific operation process to further explain 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 data sets by using dynamically adjustable minibatch data while maintaining the accuracy of the MH algorithm, and implement the hardware on the basis of fully considering the architecture of modern hardware. The following problems can be solved: 1) The traditional MCMC method has high time cost and high computational cost when processing large data sets; 2) Although other data subsampling MCMC methods can improve efficiency by approximating the acceptance rate, this approximation may introduce bias and lead to inaccurate sampling results; 3) Under the traditional design method, the batch size is usually fixed, and they may not be able to fully utilize the local characteristics of the data to optimize the amount of calculation for each iteration; 4) Under the traditional design method, cross-memory access of big 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, this formula can be expressed as:

[0074]

[0075] Referring to the SS-MCMC pseudo code 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, and 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 accept 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, by repeating steps 2 to 8 in a loop, each loop is 1 iteration, and the parameter sample θ corresponding to the i-th iteration is generated. i .

[0080] Step 2: Generate candidate parameter sample θ' for the i-th iteration, and calculate 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 ,θ').

[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] In the formula, 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 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 adding the previous parameter θ i-1 Add a mean of 0 and a variance of s to the 2 The update is based on 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 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 in the training data set according to the number B of mini-batch sub-sampling, and update the mini-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 θ'. Specifically, it includes:

[0090] Step 4.1: Sample the training data set according to the weight of each training data point, 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 of the nth training data being 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 is n bThe parameters of the Bayesian model for the data points are θ' and θ i-1 The component energy at Indicates index is n b The weight of the data point;

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

[0097]

[0098] In the formula, Indicates index is 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 data set 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 have an 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 i-th iteration of the Bayesian model. i , that is, θ i = θ', and let the value of the acceptance flag accumulator current accept increase by 1; 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, that is, θ 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, for every 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 adjust 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 acceptance number:

[0107] If yes, 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 θ of 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 data set D, and θ represents the parameters of the Bayesian model. This function measures the "energy" or "cost" of the observed data point under given parameters θ, and can be used to measure the fit between the parameters and the data. n The smaller the value of (θ), the closer the parameter θ is to the observed data x. n The better the fit. 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 the symmetric function For any θ,θ'∈Θ, the energy difference is bounded by |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 small batches of data points, thus greatly reducing the amount of computation.

[0114] In this embodiment, it is assumed that each energy function U n YesLn -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 may 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 pseudo code shown in Table 2), the probability in line 9 of the pseudo code shown in Table 2 is updated through one 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 that only calculates the probability of a 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 batch size and convergence rate of the algorithm. A larger χ value may increase the batch size and improve the convergence rate, but may also lead to higher computational costs. A smaller χ value may reduce the amount of computation per iteration and reduce the convergence rate. One of the key innovations of the SS-MCMC acceleration method of the present invention is the self-shrinking mechanism of the batch size, as shown in lines 20 to 24 of the pseudocode shown in Table 2. This mechanism can dynamically adjust the required subset size by adjusting the hyperparameter according to the sample acceptance rate, ensuring the best 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 Embodiment 1 on an FPGA. Figure 1As shown, it includes: sample proposal evaluator, small 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 θ' of the i-th iteration, calculate the candidate parameter sample θ' and the parameter sample θ of 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 used to: read the energy difference 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 in the training data set 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 pre-fetching technique that can best utilize on-chip memory is introduced, as shown in Table 2. The SS-MCMC pseudo code is based on 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 it is 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 is 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 As shown. This module is designed with a pipeline method and can effectively calculate the likelihood of each data point. The likelihood probability calculation module completes the probability of data removal and the possibility calculation of line 9 of pseudo code 2. If the data is removed, the likelihood probability will be filled with 0. If the data is selected, its likelihood probability is calculated and output.

[0128] The algorithm iterates i-1 to generate 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 simultaneously process multiple likelihood calculation submodules. The likelihood probability calculation of this embodiment has two likelihood calculation submodules, namely Figure 3 Likelihood module 1 and likelihood module 2 in the example respectively calculate candidate parameter samples θ' and current parameters θ in each iteration. i-1 The likelihood is based on the mini-batch used in the i-th iteration. This is because SS-MCMC requires

[0129] The parameter sample generation module is used to: calculate the acceptance probability α of the candidate parameter sample θ' using the likelihood LH of the acceptance probability finally updated by the small 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, that is, θ i =θ i-1 .

[0130] The acceptance flag accumulator is used to count the number of acceptances of the proposed parameter samples within each interval of a given number of iterations;

[0131] The hyperparameter adjustment module is used to dynamically adjust the hyperparameter χ according to the number of times the candidate parameter samples are accepted within each interval of a given number of iterations.

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

[0133] In this embodiment, a hardware accelerator is implemented on an FPGA to implement the MCMC acceleration method based on adaptive data subsampling described in Embodiment 1. The specific implementation steps are as follows:

[0134] 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, 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, 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 mini-batch size B, waiting for the next operation.

[0136] S3: After S1 is executed, the small batch sampler will start sampling without relying on the execution result of S2. When the small batch sampler runs for the first time, it will read the weight from the head address of the BRAM into the cumulative counter, and store the calculated cumulative weight in the adjacent position of the data; then it will take the weight of the next data and input it into the cumulative counter, and also store the result in the adjacent position of the data; until the cumulative 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 and multiplies it with the local limit parameter C, and then compares it with the cumulative weight of the data one by one until the value of the cumulative weight is compared to be not less than the product value, and finally outputs the data and weight value corresponding to the cumulative weight to the FIFO on the chip. This process will be executed continuously to ensure that the data pipeline meets the design requirements.

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

[0138] S5: When the number of executions of S4 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.

[0139] Example verification:

[0140] Logistic regression is used on the MNIST dataset of handwritten digits to evaluate the performance of the SS-MCMC accelerator of the present invention. It is also compared with several different MCMC methods (other methods are also deployed on FPGA devices). The remaining methods are:

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

[0142] flyMCMC: A commonly used contrasting subsampling MCMC method proposed in the paper “Firefly Monte Carlo: ExactMCMC with Subsets of Data”.

[0143] TunaMH: A state-of-the-art subsampling MCMC method proposed in the “Asymptotically Optimal ExactMinibatch Metropolis-Hastings” paper.

[0144] A standard MH algorithm run was first used to determine the true prediction mean, and then risk was calculated as the mean squared error (MSE) of the estimates 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 fast 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 the present invention and other methods on examples is compared, and the line width in the figure represents the standard deviation observed in these runs (each algorithm is performed 100 iterations). The results confirm that the SS-MCMC method of the present invention is consistently better than the MH and flyMCMC methods in terms of test accuracy. Although the SS-MCMC of the present invention and TunaMH achieve similar high accuracy, the SS-MCMC of the present invention reaches this level faster, highlighting its efficiency in producing accurate results.

[0146] Figure 6 The performance of the present invention, 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 Focus on the first-dimensional probability distribution of parameter samples generated by these methods. The experimental results show that the mean of the SS-MCMC acceleration method of the present invention is closest to the mean of the conventional MH method, which proves the unbiasedness of the SS-MCMC acceleration method of the present invention from an experimental perspective. This finding further verifies 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. An MCMC acceleration method based on adaptive data subsampling, characterized in that: Applications to Bayesian modeling, including: Step 1: Given the initial parameter sample θ0 and initial hyperparameter χ0 of the Bayesian model; Step 2: Generate candidate parameter sample θ' for the i-th iteration, and calculate 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 ,θ'); Step 3: According to the upper bound of energy difference M(θ i-1 ,θ') and the current hyperparameter χ, generate the number of mini-batch sub-sampling B; 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, and update the likelihood LH of the acceptance probability of the candidate parameter sample θ'; 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 θ'; 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 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, that is, θ i =θ i-1 ; Step 7, for each interval with a given number of iterations, the hyperparameter χ is dynamically adjusted according to the number of candidate parameter samples accepted in the interval; Step 8, let i = i + 1, return to step 2, and continue to generate the parameter sample θ of 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 .

2. The MCMC acceleration method based on adaptive data subsampling according to claim 1, characterized in that: Step 2 uses a Gaussian random walk with a step size of s to generate candidate parameter samples θ', expressed as: θ'~θ i-1 +Normal(0,s 2 I D ) In the formula, Normal(0,s 2 I D ) represents a normally distributed random variable, and the mean vector of the normal distribution is 0, and the covariance matrix is ​​s 2 I D ; Among them, I D is a D×D identity matrix.

3. The MCMC acceleration method based on adaptive data subsampling according to claim 1, characterized in that: The upper bound of energy difference M(θ i-1 ,θ') is calculated as: M(θ i-1 ,θ')=||θ'-θ||.

4. The MCMC acceleration method based on adaptive data subsampling according to claim 1, characterized in that: The number of mini-batch sub-samples B is sampled from a Poisson distribution and is expressed as: B~Poisson(χC 2 M 2 (i i-1 ,θ')+CM(θ i-1 ,i')) Where Poisson represents Poisson distribution, and C is the sum of the weights of the training data set.

5. The MCMC acceleration method based on adaptive data subsampling 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 n of the sampled data point b satisfy: P(n b =n)=c n / C; In the formula, P(n b =n) represents the probability of the nth training data being 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; Step 4.2, calculate two intermediate variables φ1 and φ2: Where U(θ) is the component energy function of the Bayesian model, and Indicates that the index is m b The parameters of the Bayesian model for the data points are θ' and θ i-1 The component energy at Indicates index is n b The weight of the data point; Step 4.3, calculate the acceptance probability of each data point: In the formula, Indicates index is n b 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 u1, u1~Uniform(0,1); if Then add the data point to the mini-batch data set and update the likelihood LH of the acceptance probability of the candidate parameter sample θ'; 6. The MCMC acceleration method based on adaptive data subsampling 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 as the parameter sample θ of the i-th iteration of the Bayesian model in step 6 i , then the value of the acceptance flag accumulator is increased by 1; Step 7 is specifically as follows: at each 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.

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 a candidate parameter sample θ' of the i-th iteration, calculate the candidate parameter sample θ' and the parameter sample θ of 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; The sample proposal evaluator is further used to: read the energy difference upper bound M(θ i-1 ,θ') and the current hyperparameter χ, and according to the read M(θ i-1 ,θ') and χ generate the number of mini-batch sub-sampling B; The mini-batch sampler is used to: sample in the training data set according to the number B of mini-batch sub-sampling; 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 θ'; The parameter sample generation module is used to: calculate the acceptance probability α of the candidate parameter sample θ' using the likelihood LH of the acceptance probability finally updated by the small 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, that is, θ i =θ i-1 ; The hyperparameter adjustment module is used to dynamically adjust the hyperparameter X according to the number of times the candidate parameter samples are accepted in each interval of a given number of iterations; The random number generator is used to generate a random number u2 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, wherein 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: An address mapping table stored in the cache is also included, which is used to map the training data to their actual locations in the on-chip memory or the off-chip memory.

10. The hardware accelerator according to claim 7, characterized in that: The likelihood probability calculation module includes two likelihood calculation submodules, which are used to calculate the candidate parameter sample θ' and the parameter sample θ i-1 The likelihood of acceptance probability.

Images