Convolutional non-negative pattern learning method based on maximum boundary likelihood estimation under generalized Kullback-Leibler divergence

Abstract

The invention discloses a convolutional non-negative pattern learning method based on maximum boundary likelihood estimation under generalized Kullback-Leibler (KL) divergence, and the method comprises the steps: researching a likelihood function of a convolutional non-negative pattern under the generalized Kullback-Leibler divergence, and carrying out the product division of an activation function from a model in an integral mode, according to the method provided by the invention, the local correlation of the columns in the observation data V can be modeled, and the used maximum likelihood estimation method has the characteristic of automatic model order selection. When an integral without an analytic solution is processed, the lower boundary of the model is continuously optimized and improved through a variational inference method, the target of maximum marginal likelihood is achieved, a final result not only keeps the simplicity of an observation data form, but also has the computational complexity obviously lower than that of a sampling algorithm represented by a Markov chain Monte Carlo algorithm, and the method has the advantages that the method is simple and easy to implement. And the calculation complexity is similar to that of a convolution non-negative matrix factorization algorithm based on multiplication iteration.

Description

Technical Field

[0001] The present invention belongs to the technical field of convolution learning, and in particular relates to a convolution non-negative pattern learning method based on maximum boundary likelihood estimation under generalized Kullback-Leibler divergence. Background Art

[0002] The goal of non-negative matrix factorization (NMF) is to decompose non-negative data V of size F×N into non-negative matrices W and H of size F×K and K×N. When describing the difference between V and WH, scholars have proposed a lot of objective functions. Among them, the generalized (Kullback-Leibler, KL) divergence is a commonly used measurement method. The generalized KL divergence is widely used in non-negative matrix factorization because it can better handle the relationship between time-frequency points with small values ​​and time-frequency points with large values. In addition, assuming that the observed data V follows the generative model of the Poisson distribution, minimizing the generalized KL divergence between V and WH is exactly equivalent to the maximum likelihood estimation. However, the non-negative pattern W learned by non-negative matrix factorization has the problem of asymptotic optimality. This is because, in the model, the number of parameters that need to be estimated by maximum likelihood is F×K+K×N, that is, the number of parameters will continue to grow with the increase of the number of observations N. In order to solve this problem, it is necessary to consider the marginal likelihood function about W. In order to deal with the integral in the maximum marginal likelihood estimation (MMLE), variational Bayes and Monte-Carlo methods are often used. Studies have shown that the automatic model order selection contained in the maximum marginal likelihood estimation method can avoid the overfitting problem of the model. However, when using NMF to process time series signals, the learned non-negative pattern W still encounters some problems. Usually, time correlation can be added to the activation function, such as introducing sparsity, time continuity or Markov model-based methods for the activation function. However, these methods are not suitable for solving the marginal likelihood function because the activation function in the model has been divided by integral. As another alternative, convolutional non-negative matrix factorization (CNMF) is used to capture repetitive patterns in time series data. Convolutional non-negative matrix factorization introduces an additional time dimension for each basis function on the basis of the original basis function, so that the basis function is transformed from a basic basis vector to a continuous frame form; however, the current maximum likelihood estimation based on CNMF has the technical problem of being unable to directly extract data basis functions with short-term correlation from the observed data, and needs to be optimized by alternately optimizing the basis function and the activation function, resulting in the problem that the total number of parameters that need to be estimated will increase with the increase in the number of observed samples.On the other hand, many of the growing parameters to be estimated are useless, and storing them will waste computer computing resources. The present invention solves the problem of how to directly extract convolution non-negative basis functions with short-time time series structure from observed data, avoiding the problem of asymptotic optimality in traditional optimization algorithms. Summary of the invention

[0003] The purpose of the present invention is to overcome the deficiencies of the prior art and to provide a convolutional non-negative pattern learning method based on maximum marginal likelihood estimation under generalized Kullback-Leibler divergence.

[0004] To achieve the above object, the technical solution adopted by the present invention is:

[0005] A convolutional non-negative pattern learning method based on maximum marginal likelihood estimation under generalized Kullback-Leibler divergence includes the following steps:

[0006] S1. Based on the generative model, the generative model of the observed target is assumed to obtain the generative model of the observed target. The expression is as follows:

[0007]

[0008] In the formula, v f,n is an element in V, w(t) f,k represents the elements in W, represents a Poisson distribution, and Γ(·) represents the Gamma distribution, h k,n-t represents the elements in the matrix obtained by shifting the matrix H right by t columns;

[0009] S2. Assuming that the observed data are independently distributed and conditionally dependent on W and H, the generative model in step S1 is written as the following combination:

[0010]

[0011]

[0012] In the formula, c t,k,f,n is the implicit source, v f,n is an element in V, w(t) f,k represents the elements in W, h k,n-t represents the elements in the matrix obtained by shifting the matrix H right by t columns, represents Poisson distribution;

[0013] S3. According to the generative model obtained in step S2, the log-likelihood of the generative model is calculated. The expression of the log-likelihood is:

[0014] CML (V|W)=logp(V|W)

[0015] =log∫ H p(V|W,H)p(H)dH

[0016] In the formula, C ML (V|W) represents the log-likelihood of V given the parameter W;

[0017] S4. Use Jason's inequality to derive the lower bound of the log-likelihood expression obtained in step S3, and regard the obtained lower bound as a function of W. The expression of the lower bound is:

[0018] Q(W)∝∫ C,H q(C,H)logp(V,C,H|W)dCdH

[0019] Where the symbol ∝ indicates equality under the condition of ignoring the constant term, q(C,H) is the distribution of the tools used, and p(V,C,H|W) represents the joint probability distribution under the given W condition;

[0020] S5. Based on the lower boundary expression of step S4, in order to solve the integral of the lower boundary expression in step S4, a tool distribution q(C, H) in a factorized form is selected. The expression of the tool distribution q(C, H) is:

[0021]

[0022] In the formula, q(C f,n ) is the probability p t,k,f,n is the multinomial distribution of probability, q(h k,n ) is Gamma distribution, C f,n Represents a matrix of size T×K;

[0023] S6. According to the principle of variational Bayesian EM algorithm, according to the expression of the tool distribution q(C,H) selected in step S5, q(h k,n ) and q(C f,n ) and obtain the update formula of the optimal tool distribution of q(C f,n ) and q(h k,n ), the sufficient statistic is expressed as:

[0024]

[0025] In the formula, and is a sufficient statistic for the instrument distribution, and Where ψ represents the digamma function, which is expressed as w(t) f,k represents the elements in W, c t,k,f,n+t Represents the elements in the tensor C, α k is the shape parameter of the Gamma distribution, β k is the scale parameter of the Gamma distribution;

[0026] S7. Based on Bayes' rule, the lower bound is regarded as an objective function about W, and the constant term irrelevant to W is discarded. Then the lower bound expression of step S4 can be rewritten as follows:

[0027]

[0028] In the formula, the symbol ∝ indicates equality under the condition of ignoring the constant term, <·> q represents the expectation under a given q distribution, c t,k,f,n represents the elements in the tensor C, w(t) f,n represents the elements in W, h k,n-t represents the elements in the matrix obtained by shifting the matrix H right by t columns;

[0029] S8, combined with step S6 to obtain the optimal tool distribution q (C f,n ) and q(h k,n ) and the lower boundary expression obtained in step S7, the update criterion is obtained, and the formula of the update criterion is as follows:

[0030]

[0031] Where W(t) represents the matrix representation of the t-time slice in W. It represents the result after shifting all the columns in the matrix H to the right by t columns, and 1 represents an all-1 matrix of the same size as V.

[0032] Preferably, step S6 comprises the following steps:

[0033] S61, the parameter factor q(h) in the expression of the tool distribution q(C,H) selected in step S5 k,n ) and q(C f,n ) has the following form:

[0034]

[0035] In the formula, <·> π It represents the expectation under the probability distribution π, symbol A -i,j It means that the elements in A are removed after a i,j The remaining elements, logp(V,C,H|W) represents the log-likelihood of the joint probability density distribution with W as the parameter;

[0036] S62, expand logp(V,C,H|W) in step S61 to obtain the following form:

[0037]

[0038] In the formula, δ(·) represents the Kronecker delta function, v f,n is an element in V, w(t) f,k represents the elements in W, h k,n-t represents the elements in the matrix obtained by shifting the matrix H right by t columns, c t,k,f,n Represents the elements in the tensor C, h k,n represents the elements in the matrix H, α k is the shape parameter of the Gamma distribution, β k is the scale parameter of the Gamma distribution;

[0039] S63, according to the expected expression of step S61 and the expansion obtained in step S62, derive q(h k,n ) and q(C f,n )The update formula of the optimal tool distribution is obtained, and q(C f,n ) and q(h k,n ), the sufficient statistic is expressed as:

[0040]

[0041] In the formula, and is a sufficient statistic for the instrumental distribution, and Where ψ represents the digamma function, which is expressed as w(t) f,k represents the elements in W, c t,k,f,n+t Represents the elements in the tensor C, α k is the shape parameter of the Gamma distribution, β k is the scale parameter of the Gamma distribution;

[0042] Preferably, in step S6, the q(h k,n ) and q(C f,n )The update formulas of the optimal tool distribution are as follows:

[0043]

[0044] In the formula, <·> represents the expectation under a given distribution, h k,n represents the elements in the matrix H, C is a four-dimensional tensor representing the implicit source, w(t)f,k represents the elements in W, c t,k,f,n is the implicit source, c t,k,f,n+t Represents the elements in the tensor C, α k is the shape parameter of the Gamma distribution, β k is the scale parameter of the Gamma distribution; h k,n-t represents the elements in the matrix obtained by shifting the matrix H right by t columns, C f,n represents a matrix of size T×K, q(C f,n ) is the probability p t,k,f,n is the multinomial distribution of probability,

[0045] Preferably, step S8 comprises the following steps:

[0046] S81, optimize the lower boundary obtained in step S7 to maximize it, obtain the partial derivative, and set it equal to 0, and obtain the following relationship:

[0047]

[0048] In the formula, Q(W) represents the part of the lower boundary to be optimized that is related to W, and h k,n represents the elements in the matrix H, C represents the four-dimensional tensor of the implicit source, whose elements are c t,k,f,n ,w(t) f,k represents the elements in W, h k,n-t represents the elements in the matrix obtained by shifting the matrix H right by t columns;

[0049] S82. According to the relationship obtained in step S81, an update formula of the element W is derived. The update formula is as follows:

[0050]

[0051] In the formula, h k,n represents the elements in the matrix H, C is a four-dimensional tensor representing the implicit source, and its elements are c t,k,f,n ,w(t) f,k represents the elements in W, h k,n-t represents the elements in the matrix obtained by shifting the matrix H right by t columns;

[0052] S83, combining q(C f,n ) and q(h k,n ) and the update formula obtained in step S82 to obtain the update criterion, and the formula of the update criterion is as follows:

[0053]

[0054] Where W(t) represents the matrix representation of the t-time slice in W. It represents the result of shifting all the columns in the matrix H to the right by t columns, and 1 represents an all-1 matrix of the same size as V.

[0055] Compared with the prior art, the present invention has the following beneficial effects:

[0056] The method proposed in the present invention can directly learn the convolution non-negative basis function with local correlation from the observed data V, avoiding the suboptimal solution method of alternating optimization required by the traditional convolution non-negative matrix decomposition, and maintaining the inherent model order selection characteristics of the maximum marginal likelihood method. Through the variational inference method, the huge amount of calculation of the sampling method is avoided, and the marginal likelihood estimation can be maximized with a calculation amount comparable to that of the conventional CNMF algorithm. Experimental results show that when K takes a sufficiently large value, the proposed algorithm can obtain more accurate and more interpretable results than the traditional CNMF method. BRIEF DESCRIPTION OF THE DRAWINGS

[0057] Figure 1 This is a graph of the synthetic dataset;

[0058] Figure 2 The convolution non-negative patterns learned from the synthetic data sets by the method provided in the embodiment of the present invention and the existing MM1 algorithm respectively;

[0059] in, Figure 2 (a) is the method of the present invention; Figure 2 (b) is the existing MM1 algorithm;

[0060] Figure 3 The vectors learned by the method provided in the embodiment of the present invention and the existing MM1 algorithm respectively;

[0061] in, Figure 3 (a) is the method of the present invention; Figure 3 (b) is the existing MM1 algorithm. DETAILED DESCRIPTION

[0062] The technical solutions in the embodiments of the present invention are described clearly and completely below. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0063] 1. Introduction to the basic model related to the present invention

[0064] 1. Basic assumptions

[0065] Based on the generative model, the generative model of the observed target is assumed to obtain the generative model of the observed target. The expression is as follows:

[0066]

[0067] In the formula, v f,n is an element in V, w(t) f,k represents the elements in W, h k,n-t represents the elements in the matrix obtained by shifting the matrix H right by t columns, represents the Poisson distribution, which is defined as follows:

[0068]

[0069] Where Γ(·) is the Gamma function.

[0070] Assume that the observed data are independently distributed and conditionally dependent on W and H. Since Poisson distribution has superposition, the hidden source also obeys Poisson distribution. Therefore, the present invention uses a combined model to represent the generative model of the above observed target, and the expression is as follows:

[0071]

[0072]

[0073] In the formula, c t,k,f,n is the implicit source, v f,n is an element in V, w(t) f,k represents the elements in W, h k,n-t represents the elements in the matrix obtained by shifting the matrix H right by t columns, represents Poisson distribution;

[0074] Activation function coefficient h k,n is a random variable that follows a Gamma distribution, and its expression is as follows:

[0075]

[0076]

[0077] In the formula, is the Gamma distribution;

[0078] Because of the conjugate relationship between the Gamma distribution and the Poisson distribution, the Gamma distribution is preferred.

[0079] 2. Maximum Margin Likelihood

[0080] According to the obtained generative model, the marginal log-likelihood of the generative model with respect to W is calculated. The log-likelihood can be written as follows:

[0081]

[0082] In the formula, C ML (V|W) represents the log-likelihood of V given the parameter W;

[0083] The integral in the formula is difficult to solve, that is, the analytical solution of the marginal model cannot be obtained. It should be pointed out that in the current setting, the problem studied is not the full Bayesian interference of integrating all variables (including those in W), but treating W as a deterministic parameter and focusing on the marginal likelihood of integrating H.

[0084] 3. Lower Bound of Marginal Likelihood

[0085] Because the integral is difficult to solve and the analytical solution after the integral calculation cannot be obtained, we can resort to Jason's inequality and define the lower bound of the log-likelihood expression obtained above. Using latent variables C and H, the data V is expanded to obtain the following relationship:

[0086]

[0087] Therefore, the lower bound of the objective function can be obtained, and the lower bound is expressed as follows:

[0088] Q(W)∝∫ C,H q(C,H)logp(V,C,H|W)dCdH(9)

[0089] Where ∝ represents equality when the constant term is ignored, p(V,C,H|W) represents the joint probability distribution under the given W condition, and q(C,H) is the instrumental distribution. When q(C,H) is equal to the posterior distribution p(C,H|V,W), the inequality becomes an equality. The maximum marginal likelihood can be achieved by optimizing this lower bound. In order to maximize the integral in formula (9), its partial derivative is calculated and set to 0, and the exact posterior p(C,H|V,W) is approximated by the variational approximation q(C,H), and the expectation maximization algorithm is used to maximize this lower bound.

[0090] 2. Specific construction process of the method provided in the embodiment of the present invention

[0091] The embodiment of the present invention provides a convolutional non-negative pattern learning method based on maximum marginal likelihood estimation under generalized Kullback-Leibler divergence, comprising the following steps:

[0092] 1. Mean Field Approximation

[0093] In the variational method, mean field approximation is a powerful technical means. Its core idea is to use the complete factorization of all latent variables to approximate the posterior probability. According to the mean field approximation theory and based on the lower bound expression, the instrument distribution q(C,H) in factorized form can be used. The expression of the instrument distribution q(C,H) is:

[0094]

[0095] In the formula, q(C f,n ) is p t,k,f,n is a multinomial distribution of probability, and q(h k,n ) are the shape parameter and scale parameter respectively and Gamma distribution, C f,n represents a matrix of size K×T, and C f,n The expression is as follows:

[0096]

[0097] Where C is a tensor with four dimensions.

[0098] Factor q(h k,n ) and q(C f,n ) should meet the following expectations:

[0099]

[0100]

[0101] In the formula, <·> π represents the expectation under the probability distribution π, A -i,j Indicates that except for element a i,j The set consisting of all remaining elements in A, logp(V,C,H|W) represents the log-likelihood of the joint probability density distribution with W as the parameter.

[0102] 2. Step E

[0103] According to the principle of variational Bayesian EM algorithm, in order to calculate the expectation of equations (12) and (13), logp(V,C,H|W) is expanded into the following form:

[0104]

[0105] Where δ(·) is the Kronecker function, v f,n is an element in V, w(t) f,k represents the elements in W, h k,n-t represents the elements in the matrix obtained by shifting the matrix H right by t columns, c t,k,f,n Represents the elements in the tensor C, h k,n represents the elements in the matrix H, α k is the shape parameter of the Gamma distribution, β k is the scale parameter of the Gamma distribution;

[0106] According to formulas (12)-(14), we can get q(h k,n ) and q(C f,n )The update formula of the optimal tool distribution is as follows:

[0107]

[0108]

[0109] In the formula, <·> represents the expectation under a given distribution, h k,n represents the elements in the matrix H, C is a four-dimensional tensor representing the implicit source, w(t) f,k represents the elements in W, c t,k,f,n is the element in the implicit source, α k is the shape parameter of the Gamma distribution, h k,n-t represents the elements in the matrix obtained by shifting the matrix H right by t columns, C f,n represents a matrix of size T×K, q(C f,n ) is the probability p t,k,f,n is the multinomial distribution of probability,

[0110] The optimal tool distribution q(C f,n ) and q(h k,n ) has the following sufficient statistics:

[0111]

[0112]

[0113] In the formula, and is a sufficient statistic for the instrumental distribution, and Where ψ represents the digamma function, which is expressed as w(t)f,k represents the elements in W, c t,k,f,n represents the elements in the tensor C, α k is the shape parameter of the Gamma distribution, β k is the scale parameter of the Gamma distribution, q(C f,n ) obeys with probability p t,k,f,n The multinomial distribution is:

[0114]

[0115] Where w(t) f,k and w(τ) f,l represents the elements in W, h k,n-t represents the elements in the matrix obtained by shifting the matrix H right by t columns, h l,n-τ Represents the elements in the matrix obtained by shifting the matrix H right by τ columns.

[0116] The expectation in formula (17) <c t,k,f,n > and the expectation in (18) <logh k,n >Can be obtained by:

[0117]

[0118]

[0119] In the formula, ψ is defined as The digamma function of ; and ρ(n) can be defined by equation (22):

[0120]

[0121] 3. M-step

[0122] Using Bayes' rule, the lower bound is regarded as an objective function about W, and the constant term that is irrelevant to W is discarded. The lower bound expression (9) can be rewritten as follows:

[0123]

[0124] In the formula, the symbol ∝ indicates equality under the condition of ignoring the constant term, <·> q represents the expectation under a given q distribution, c t,k,f,n represents the elements in the tensor C, w(t) f,n represents the elements in W, h k,n-t is the element in the matrix H shifted right by t columns;

[0125] The lower bound of the objective function can be derived using the maximum marginal a posteriori (MMAP) estimation method of formula (23). The partial derivative of the lower bound with respect to a specific convolution non-negative mode is calculated and set equal to 0, and the following relationship can be obtained:

[0126]

[0127] In the formula, Q(W) represents the part of the lower boundary to be optimized that is related to W, and h k,n represents the elements in the matrix H, C represents the four-dimensional tensor of the implicit source, whose elements are c t,k,f,n ,w(t) f,k represents the elements in W, h k,n-t represents the elements in the matrix obtained by shifting the matrix H right by t columns;

[0128] Divide the negative and positive parts obtained by the partial derivative calculation above and multiply them by the dictionary element to derive the update formula for the element W. The update formula is as follows:

[0129]

[0130] In the formula, h k,n represents the elements in the matrix H, C is a four-dimensional tensor representing the implicit source, and its elements are c t,k,f,n ,w(t) f,k represents the elements in W, h k,n-t Represents the elements in the matrix obtained by shifting the matrix H right by t columns.

[0131] Combining formulas (20) and (25), the update criterion becomes:

[0132]

[0133]

[0134] Where W(t) represents the matrix representation of the t-time slice in W. It represents the result of shifting all the columns in the matrix H to the right by t columns, and 1 represents an all-1 matrix of the same size as V.

[0135] 3. The performance and effect of the method provided by the embodiment of the present invention are evaluated below

[0136] 1. Mean Field Approximation

[0137] based on Figure 1 The performance and effect of the method proposed in the embodiment of the present invention are evaluated by using synthetic data. Figure 1As can be seen, the synthetic data contains three repeating structural patterns. Each spans 20 columns in the synthetic data. The values ​​in the repeating pattern vary from 7 to 13 (shown in black or gray), and the values ​​of other white areas are all 0. The roof-like structure is called type I, while the monotonically rising and monotonically falling structures are type II and type III, respectively. It can be seen from the synthetic data that type I and type II are repeated 4 times each, while type III is only repeated 2 times.

[0138] The hyperparameter is fixed to α k =1 (i.e., the Gamma distribution degenerates into an exponential distribution), and let β k =1, that is, according to the literature (O.Dikmen and C.Févotte, "Maximum marginal likelihood estimation for nonnegative dictionary learning," in 2011 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Prague, 2011, pp.1992-1995), h k,n ~exp(-h k,n ). As a baseline method, the MM1 algorithm in the literature (D. Fagot, H. Wendt, C. Févotte and P. Smaragdis, "Majorization-minimization algorithms for convolutive NMF with the beta-divergence," In Proc. IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 8202-8206, Brighton, United Kingdom, 2019) is used for comparison. Let K = 10, that is, 10 convolution basis functions are used to learn repeated patterns. Each algorithm is repeated 1000 times, and the learned convolution non-negative patterns are sorted in descending order of the elements in the vector κ. The elements in the vector κ represent the importance of the convolution non-negative pattern, and the expression is:

[0139]

[0140] 2. B. Experimental Results

[0141] The convolution non-negative model learned by the marginal maximum likelihood estimation (MMLE) method proposed in the present invention and the MM1 method in the literature is Figure 2 Through Figure 2 It is obvious that the maximum marginal likelihood estimation can learn a very simple representation of the observed data, which only contains a few convolution basis functions. In addition, by Figure 2 It can also be seen that the MMLE method proposed in this invention can correctly learn Type I and Type III, while Type II is represented by two basis functions. Figure 2 As can be seen in the right sub-figure, the MM1 algorithm makes full use of all basis functions and learns trivial patterns that may not reflect the inherent structure of the data. Therefore, in order to avoid overfitting, when K is significantly higher than the actual model order, MMLE shows an advantage over MM1.

[0142] Further detect the elements in the vector κ and the corresponding Figure 2 The convolution non-negative mode in the result is as follows Figure 3 As shown, through Figure 3 It can be found that the MMLE method proposed in the present invention contains the ability of automatic model order selection, and the basis functions that are irrelevant to the representation of the observed data (the basis functions with smaller corresponding values ​​in the vector) can be easily removed. However, under the MM1 condition, the learned basis functions all have the same contribution in representing the observed data, making it impossible to remove any of these basis functions to achieve a simple representation of the observed data.

[0143] In summary, the present invention studies the marginal likelihood estimation problem of convolution non-negative matrix decomposition by integral product division of activation function under generalized KL divergence, and proposes a method based on maximum marginal likelihood estimation. The algorithm proposed in the present invention can model the local correlation of columns in the observed data V and maintain the inherent model order selection characteristics of the MMLE method. Through variational inference, the marginal likelihood estimation can be maximized, and its computational complexity is comparable to that of the conventional CNMF algorithm. Experimental results show that when K is a sufficiently large value, the proposed algorithm can obtain more accurate and more interpretable results than the traditional CNMF method.

[0144] Although the embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions and variations may be made to the embodiments without departing from the principles and spirit of the present invention, and that the scope of the present invention is defined by the claims and their equivalents.

Claims

1. A convolutional non-negative pattern learning method based on maximum marginal likelihood estimation under generalized Kullback-Leibler divergence, characterized in that: The following steps are involved: S1. Based on the generative model, the generative model of the observed target is assumed to obtain the generative model of the observed target. The expression is as follows: In the formula, v f,n is an element in V, w(t) f,k represents the elements in W, represents a Poisson distribution, and Γ(·) represents the Gamma distribution, h k,n-t represents the elements in the matrix obtained by shifting the matrix H right by t columns; S2. Assuming that the observed data are independently distributed and conditionally dependent on W and H, the generative model in step S1 is written as the following combination: In the formula, c t,k,f,n is the implicit source, v f,n is an element in V, w(t) f,k represents the elements in W, h k,n-t represents the elements in the matrix obtained by shifting the matrix H right by t columns, represents Poisson distribution; S3. According to the generative model obtained in step S2, the log-likelihood of the generative model is calculated. The expression of the log-likelihood is: C ML (V|W)=logp(V|W) =log∫ H p(V|W,H)p(H)dH In the formula, C ML (V|W) represents the log-likelihood of V given the parameter W; S4. Use Jason's inequality to derive the lower bound of the log-likelihood expression obtained in step S3, and regard the obtained lower bound as a function of W. The expression of the lower bound is: Q(W)∝∫ C,H q(C,H)logp(V,C,H|W)dCdH Where the symbol ∝ indicates equality under the condition of ignoring the constant term, q(C,H) is the distribution of the tools used, and p(V,C,H|W) represents the joint probability distribution under the given W condition; S5. Based on the lower boundary expression of step S4, in order to solve the integral of the lower boundary expression in step S4, a tool distribution q(C, H) in a factorized form is selected. The expression of the tool distribution q(C, H) is: In the formula, q(C f,n ) is p t,k,f,n is the multinomial distribution of probability, q(h k,n ) is the Gamma distribution, C f,n Represents a matrix of size T×K; S6. According to the principle of variational Bayesian EM algorithm, according to the expression of the tool distribution q(C,H) selected in step S5, q(h k,n ) and q(C f,n ) and obtain the update formula of the optimal tool distribution of q(C f,n ) and q(h k,n ), the sufficient statistic is expressed as: In the formula, and is a sufficient statistic for the instrument distribution, and Where ψ represents the digamma function, which is expressed as w(t) f,k represents the elements in W, c t,k,f,n+t is an element in the four-dimensional tensor C, α k is the shape parameter of the Gamma distribution, β k is the scale parameter of the Gamma distribution; S7. Based on Bayes' rule, the lower bound is regarded as an objective function about W, and the constant term irrelevant to W is discarded. Then the lower bound expression of step S4 can be rewritten as follows: In the formula, the symbol ∝ indicates equality under the condition of ignoring the constant term, <·> q represents the expectation under a given q distribution, c t,k,f,n Represents the elements in the tensor C, w(t) f,n represents the elements in W, h k,n-t represents the elements in the matrix obtained by shifting the matrix H right by t columns; S8, combined with step S6 to obtain the optimal tool distribution q (C f,n ) and q(h k,n ) and the lower boundary expression obtained in step S7, the update criterion is obtained, and the formula of the update criterion is as follows: Where W(t) represents the matrix representation of the t-time slice in W. It represents the result after shifting all the columns in the matrix H to the right by t columns, and 1 represents an all-1 matrix of the same size as V.

2. The convolutional non-negative pattern learning method based on maximum marginal likelihood estimation under generalized Kullback-Leibler divergence according to claim 1, characterized in that: Step S6 comprises the following steps: S61, the parameter factor q(h) in the expression of the tool distribution q(C,H) selected in step S5 k,n ) and q(C f,n ) has the following form: In the formula, <·> π It represents the expectation under the probability distribution π, symbol A -i,j It means that the elements in A are removed after a i,j The remaining elements, logp(V,C,H|W) represents the log-likelihood of the joint probability density distribution with W as the parameter; S62, expand logp(V,C,H|W) in step S61 to obtain the following form: In the formula, δ(·) represents the Kronecker delta function, v f,n is an element in V, w(t) f,k represents the elements in W, h k,n-t represents the elements in the matrix obtained by shifting the matrix H right by t columns, c t,k,f,n represents the elements in the tensor C, h k,n represents the elements in the matrix H, α k is the shape parameter of the Gamma distribution, β k is the scale parameter of the Gamma distribution; S63, according to the expected expression of step S61 and the expansion obtained in step S62, derive q(h k,n ) and q(C f,n )The update formula of the optimal tool distribution is obtained, and q(C f,n ) and q(h k,n ), the sufficient statistic is expressed as: In the formula, and is a sufficient statistic for the instrumental distribution, and Where ψ represents the digamma function, which is expressed as w(t) f,k represents the elements in W, c t,k,f,n+t Represents the elements in the tensor C, α k is the shape parameter of the Gamma distribution, β k is the scale parameter of the Gamma distribution.

3. The convolutional non-negative pattern learning method based on maximum marginal likelihood estimation under generalized Kullback-Leibler divergence according to claim 1, characterized in that: In step S6, the q(h k,n ) and q(C f,n )The update formulas of the optimal tool distribution are as follows: In the formula, <·> represents the expectation under a given distribution, h k,n represents the elements in the matrix H, C is a four-dimensional tensor representing the implicit source, w(t) f,k represents the elements in W, c t,k,f,n is the implicit source, c t,k,f,n+t Represents the elements in the tensor C, α k is the shape parameter of the Gamma distribution, h k,n-t represents the elements in the matrix obtained by shifting the matrix H right by t columns, C f,n represents a matrix of size T×K, q(C f,n ) is the probability p t,k,f,n is the multinomial distribution of probability, 4. The convolutional non-negative pattern learning method based on maximum marginal likelihood estimation under generalized Kullback-Leibler divergence according to claim 1, characterized in that: Step S8 comprises the following steps: S81, optimize the lower boundary obtained in step S7 to maximize it, obtain the partial derivative, and set it equal to 0, and obtain the following relationship: In the formula, Q(W) represents the part of the lower boundary to be optimized that is related to W, and h k,n represents the elements in the matrix H, C represents the four-dimensional tensor of the implicit source, whose elements are c t,k,f,n ,w(t) f,k represents the elements in W, h k,n-t represents the elements in the matrix obtained by shifting the matrix H right by t columns; S82. According to the relationship obtained in step S81, an update formula of the element W is derived. The update formula is as follows: In the formula, h k,n represents the elements in the matrix H, C is a four-dimensional tensor representing the implicit source, and its elements are c t,k,f,n ,w(t) f,k represents the elements in W, h k,n-t represents the elements in the matrix obtained by shifting the matrix H right by t columns; S83, combining q(C f,n ) and q(h k,n ) and the update formula obtained in step S82 to obtain the update criterion, and the formula of the update criterion is as follows: Where W(t) represents the matrix representation of the t-time slice in W. It represents the result of shifting all the columns in the matrix H to the right by t columns, and 1 represents an all-1 matrix of the same size as V.

Images