# A Logical Framework of the Evidence Function Approximation Associated with Relevance Vector Machine

Authors: Dawei Zou; Liang Tong; Jicheng Wang; Shoufang Fan; Jinxing Ji
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Relevance vector machine (RVM) is a Bayesian sparse kernel method for regression in statistical learning theory, which can avoid principal limitations of the support vector machine (SVM) and result in faster performance on test data while maintaining comparable generalization error. In this paper, we develop a logic framework of the evidence function approximation associated with RVM based on Taylor expansion instead of traditional technology called “completing the square.” While constructing the term of completing the square, we have to find the term of completing the square by making use of some skill, which in practice increases the difficulty in dealing with the evidence function approximation associated with RVM. The logical framework in this paper based on Taylor expansion shows some advantages compared with the conventional method of completing the square, which is easier to be enforced due to the clear logical framework and avoids the difficulty in looking for the term of completing the square intentionally. From symmetry of covariance in a multivariate Gaussian distribution and algebraic knowledge, we derive approximation and maximization of the evidence function associated with RVM, which is consistent with the previous result using the novel logical framework. Finally, we derive the EM algorithm for RVM, which is also consistent with the previous result except that we use the precision matrix as the covariance.


---

## Body

## 1. Introduction

When fitting polynomial curve to a given training data, we encounter the problem of choosing the order of the polynomial, which make us discuss an important concept called model selection. Model selection is an intractable problem in machine learning field since we have to face all kinds of complexity of the model according to the complexity of the problem being solved. We try our best to look for suitable approaches to fix the problem. For example, to avoid overfitting problem, we add regularization term [1, 2] to the sum of the squares of the errors between the prediction for each data point and the corresponding target values or partition the available data into a training set which determines the coefficient and a separate validation set, used to optimize the model complexity, also called a hold-out set. However, this is too wasteful of valuable training data and we have to look for more sophisticated methods. Model selection from a Bayesian perspective is an alternative technology [3] whose merit is that the overfitting problem with maximum likelihood can be avoided by marginalizing (summing for discrete case or integrating for continuous case) over the model parameters compared with making point estimates and models can be compared directly on the training data without the need for a validation set, which allows all available data to be used to train, avoids a lot of training operations for each model associated with crossvalidation, and multiple complexity parameters can be determined simultaneously as part of the training process.Bayesian model comparison simply involves the use of probabilities to depict uncertainty in the model choice by a consistent application of the sum and product rules of probability [4]. The core of model selection is to evaluate evidence function, assuming that all models are given equal prior probability. Relevance vector machine (RVM) [5, 6] is a Bayesian sparse kernel method for classification and regression that has many characteristics of the support vector machine (SVM) while avoiding its main limitation. Relevance vector machine (RVM) generally leads to much sparser models than support vector machine (SVM), which can result in faster performance, which can be applied in automatic emulation and get better performance [7]. In [5], the evidence function approximation associated with relevance vector machine needs some skill to compute and makes use of the logical framework called “completing the square” [8]. However, the method of completing the square has difficulty in looking for the term of completing the square, especially in complex problems in the machine learning field. The paper will develop a novel logical framework based on Taylor expansion, which is different from the classic technique of completing the square over parameter of the model for evaluating the evidence function associated with relevance vector machine shortly. The method presented in this paper is convenient to evaluate the evidence function with relevance vector machine (RVM).


## 2. Preliminary Mathematical Knowledge

Lemma 1 (Woodbury identity).
A
+
B
D
−
1
C
−
1
=
A
−
1
−
A
−
1
B
D
+
C
A
−
1
B
−
1
C
A
−
1.Lemma 2.
IfA and B are matrices of size N×M, then(1)IN+ABT=IM+ATB.Specially, IN+abT=1+aTb, where a and b are N-dimensional column vectors.Lemma 3.
∂
/
∂
x
ln
A
=
T
r
A
−
1
∂
A
/
∂
x.Lemma 4.
∂
/
∂
x
A
−
1
=
−
A
−
1
∂
A
/
∂
x
A
−
1.Lemma 5.
Ifx and y are n-dimensional column vectors, then TrxyT=xTy, where Tr⋅ denotes the trace of a matrix.The proofs of the abovementioned four lemmas refer to [9–11].Definition 1.
The operatorhT∇ is defined by hT∇0f=f and hT∇k+1f=hT∇hT∇kf for k∈N.
These symbolic powers act as operators on a function ofn variables, where ∇ is the gradient operator.
For example, we consider a two-dimensional case. Leth=h1,h2T and ∇=∂/∂x1,∂/∂x2, we have hT∇=h1∂/∂x1+h2∂/∂x2 and hT∇f=h1∂f/∂x1+h2∂f/∂x2. After the operator hT∇f acts on x, we obtain(2)hT∇fx=h1∂fx∂x1+h2∂fx∂x2=hT∇fx.
We also note that(3)hT∇2f=hT∇h1∂f∂x1+h2∂f∂x2=h1∂∂x1h1∂f∂x1+h2∂f∂x2+h2∂∂x2h1∂f∂x1+h2∂f∂x2=h12∂2f∂x12+2h1h2∂2f∂x1∂x2+h22∂2f∂x22.
The operatorhT∇2f also acts on x, and we obtain(4)hT∇2fx=h12∂2fx∂x12+2h1h2∂2fx∂x1∂x2+h22∂2fx∂x12.
Furthermore, we obtain(5)hT∇2fx=hTHh,where h=h1,h2T and H=∇2fx=∇∇fx=∂2fx/∂xi∂xj, which is called the Hessian matrix.
We usehT∇f=h1∂f/∂x1+h2∂f/∂x2 again and obtain(6)hT∇3f=hT∇h12∂2f∂x12+2h1h2∂2f∂x1∂x2+h22∂2f∂x22.
From (6), we obtain(7)hT∇3f=h13∂3f∂x13+3h12h2∂3f∂x12∂x2+3h1h22∂3f∂x1∂x22+h23∂3f∂x23.
By making use of the exerting operatorhT∇3f on variable x, we obtain(8)hT∇3fx=h13∂3fx∂x13+3h12h2∂3fx∂x12∂x2+3h1h22∂3fx∂x1∂x22+h23∂3fx∂x23.
From (8), we can draw a conclusion that hT∇3fx includes all third-order partial derivative of fx to variables.Lemma 6 (Taylor’s formula).
Letf:Bx,r⟶ℝ be a function that belongs to the class CBx,r, where Bx,r⊆ℝk. If h∈ℝk is such that h<r, then there exists θ∈0,1 such that(9)fx+h=∑m=0n−11m!hT∇mfx+1n!hT∇nfx+θh.The proof of Lemma6 is found in [12].Corollary 1.
Iff:Bw,r⟶ℝ is a quadratic form defined by(10)fw=wTAw+wTb+c,where A is an n×n inverse and symmetric matrix, b and w=w1,w2,⋯,wnT are n-dimensional column vectors, and c is a constant. Then, the Taylor expansion of fw at the stationary point w0=1/2A−1b is given by(11)fw=fw0+12w−w0THw−w0,where the element of Hessian matrix H is given by ∂2f/∂wi∂wj.Proof.
The Taylor formula (9) implies(12)fw=fw0+hT∇fw0+12!hT∇2fw0+13!hT∇3fw0+θh,where h=w−w0.For obtaining the stationary point forfw, setting∇fw=2Aw+b=0, we obtain the stationary point w0=1/2A−1b for fw. Since fw is a quadratic function of w, we can derive hT∇3fx=0 at an arbitrary variable x. From the abovementioned operations, we obtain (11).The equality of (11) plays a central role in our discussion because it is clearly different from conventional technology called “completing the square.” We will build up a novel logical framework and discuss its merit compared with the traditional method when dealing with evaluating the evidence function with relevance vector machine (RVM) shortly.


## 3. Evaluation and Maximization of the Evidence Function with RVM

### 3.1. Linear Basis Function Models and Bayesian Linear Regression

The simplest linear model for regression is the one that includes a linear combination of the input variables:(13)yx,w=w0+w1x1+⋯+wDxD.However, a linear combination of the input variables imposes significant limitations when we use the model. So, we extend the model by using linear combinations of fixed nonlinear functions of the input variables, which take the form(14)yx,w=w0+∑i=1M−1wiϕix,where ϕix is defined as basis function and w0 is called a bias parameter. For convenience, we define an additional dummy “basis function” ϕ0x=1 and rewrite(15)yx,w=∑i=0M−1wiϕix=wTϕx,where w=w0,⋯,wM−1T and ϕx=ϕ0x,⋯,ϕM−1xT [8].The target variablet is given by a deterministic function yx,w with additive Gaussian noise such that(16)t=yx,w+ε,where ε is a zero mean Gaussian random variable with precision (inverse variance) β. So, we obtain(17)ptx,w,β=Ntyx,w,β−1.Now, we consider a data set of inputsX=x1,⋯,xN with corresponding target value t=t1,⋯,tN. From the assumption that these data points are drawn independently from the Gaussian distribution (17), we obtain the following likelihood function:(18)ptX,w,β=∏n=1NNtnwTϕxn,β−1.Since (18) is the exponential of a quadratic function of w, we can define a corresponding conjugate prior given by(19)pw=Nwm0,S0.From (18) and (19), we obtain the posterior distribution in the form(20)pwt,X,β=NwmN,SN,where(21)mN=SNS0−1m0+βΦTt,(22)SN−1=S0−1+βΦTΦ,where Φ is an N×M matrix, called the design matrix, whose elements are given by Φnj=ϕjxn so that(23)Φ=ϕ0x1…ϕM−1x1⋮⋱⋮ϕ0xN⋯ϕM−1xN.The proofs of (21) and (22) in detail refer to Appendix A.


### 3.2. Bayesian Model Comparison and Evidence Approximation

Assume that we face a set ofL models Mi, where i=1,⋯,L. Our task is to compare them and choose an optimal model based on some standard from the Bayesian perspective, which can avoid the overfitting phenomenon associated with maximum likelihood. The uncertainty is expressed through a prior probability distribution pMi. Given a training set, we suppose that all models are given the equal prior probability, which is feasible because there should be no preference for L given models in practice. We then wish to evaluate the posterior distribution:(24)pMiD∝pMipDMi,where pMiD is defined as model evidence, which is also called the marginal likelihood since it can be viewed as a likelihood function over the space of models, in which the parameters have been marginalized out.Now, we introduce prior distribution over the hyperparametersα and β to complete a fully Bayesian treatment of the linear basis function model and we can obtain prediction by marginalizing with respect to the hyperparameters α and β and the parameter w. However, the complete marginalization over all of these variables is analytically intractable. We then discuss an estimation, in which we determine the hyperparameters α and β by maximizing the marginal likelihood function obtained by first integrating over the parameters w. The framework is viewed in the statistics field as empirical Bayes [13, 14], or type-2 maximum likelihood [15], or generalized likelihood [16], and is also called the evidence approximation in the machine learning literature [8, 17, 18].Here, if we introduce hyperprior over theα and β, then we obtain predictive distribution by marginalizing over the hyperparameters α and β and the parameter w as follows:(25)ptt=∭ptx,w,βpwt,α,βpα,βtdwdαdβ,where ptx,w,β is defined by (17) and pwt,α,β is given by (20) with mN and SN defined by (21) and (22).If the posterior distributionpα,βt is sharply peaked around α^ and β^, we can obtain the predictive distribution simply by marginalizing over w in which α and β are set to the value α^ and β^:(26)ptt≃ptt,α^,β^=∫ptx,w,β^pwt,α^,β^dw.The posterior distribution forα and β is given by(27)pα,βt∝ptα,βpα,β.We can obtain the valuesα^ and β^ by maximizing the marginal likelihood function ptα,β in the evidence framework if the prior is relatively flat. We continue evaluating the marginal likelihood for the linear basis function and obtain its maxima which will allow us to determine the values of these hyperparameters from the training data set only.


### 3.3. Evaluation of the Evidence Function with RVM

Relevance vector machine (RVM) is a Bayesian sparse kernel method for regression. The relevance vector machine for regression is a linear model of the form given by (15). The following likelihood function is given by (18).The weight prior takes the form(28)pwα=∏n=1NNwi0,αi−1,where α=α1,⋯,αMT and αi is hyperparameter and represents the precision of the corresponding parameter wi.By making use of (18), we obtain(29)ptX,w,β=∏n=1NNtnwTϕxn,β−1=β2πN/2exp−β2∑i=1Ntn−wTϕxn2=β2πN/2exp−β2t−Φw2.Similarly, by making use of (28), we obtain(30)pwα=∏n=1NNwi0,αi−1=α1α2⋯αM2πM/2exp−12∑i=1Mαiwi2=A1/22πM/2exp−12wTAw.We obtain the posterior distribution for the weights, which is also Gaussian and takes the form(31)pwt,X,α,β=Nwm,Σ,where the covariance and mean are given by(32)m=βΣΦTt,(33)Σ=A+βΦTΦ−1,where A=diagαi.The proofs of (32) and (33) in detail refer to Appendix B.We can use evidence approximation, which is also known as type-2 maximum likelihood, to determine the values of hyperparametersα and β. We first obtain the marginal likelihood function by integrating out the weight parameters:(34)ptX,α,β=∫ptX,w,βpwαdw,and then maximize the marginal likelihood function.We substitute (32) and (33) into (34) and then obtain the following result:(35)ptX,α,β=β2πN/2A1/22πM/2∫exp−Ewdw,where we have defined(36)Ew=β2t−Φw2+12wTAw=β2t−ΦwTt−Φw+12wTAw=β2tTt−2wTΦTt+wTΦTΦw+12wTAw.For the stationary point ofEw, we first derive the gradient of ∇Ew by using Corollary 1:(37)∇Ew=βΦTΦ+Aw−βΦTt,(38)∇∇Ew=A+βΦTΦ=Σ−1.Let(39)∇Ew=βΦTΦ+Aw−βΦTt=0.From (39), we obtain(40)w=ββΦTΦ+A−1ΦTt=βΣΦTt=m.From (11), (38), and (40), it follows that(41)Ew=Em+12w−mTΣ−1w−m,together with(42)Em=β2t−Φm2+12mTAm.We substitute (41) into (35) and then get the following form:(43)ptX,α,β=β2πN/2A1/22πM/2∫exp−Emexp−12w−mTΣ−1w−mdw=β2πN/2A1/22πM/2exp−Em∫exp−12w−mTΣ−1w−mdw=β2πN/2A1/2Σ1/2exp−Em,where we have used 1/2πM/2Σ1/2∫exp−1/2w−mTΣ−1w−mdw=1.From (32) and (33), we obtain(44)βΣΦTt=Σ−1m.We substitute (33) into (44) and obtain(45)βΣΦTt=A+βΦTΦm.We expand (45) and then obtain(46)βΣΦTt=Am+βΦTΦm.We multiply both sides of (46) by mT and derive(47)βmTΣΦTt=mTAm+βmTΦTΦm.From (42), we obtain(48)Em=β2t−Φm2+12mTAm=12βt−ΦmTt−Φm+mTAm=12βtTt−2βmTΦTt+βmTΦTΦm+mTAm=β2tTIN−mTΦTt=12tTβIN−β2ΦΣΦTt,where we make use of (32) and (47).According to Lemma1, we obtain(49)βIN−β2ΦΣΦT−1=βIN−βΦA+βΦTΦ−1βΦT−1=β−1IN+ΦA−1ΦT.By making use of (48) and (49), we obtain(50)ptX,α,β=β2πN/2A1/2Σ1/2exp−Em=β2πN/2A1/2Σ1/2exp−12tTβ−1IN+ΦA−1ΦT−1t=12πN/2β−1IN+ΦA−1ΦT1/2exp−12tTβ−1IN+ΦA−1ΦT−1t,where we have used the following result according to Lemma 2 and (33):(51)β−1IN+ΦA−1ΦT=1βNIN+βΦA−1ΦT=1βNIM+βA−1ΦTΦ=1βNA−1A+βΦTΦ=1βNA−1Σ−1=1βNA−1Σ−1.From (50), we obtain(52)ptX,α,β=Nt0,C,together with C=β−1IN+ΦA−1ΦT.From (43), we can then write the log of the marginal likelihood in the form:(53)lnptX,α,β=N2lnβ−Em+12∑i=1Mlnαi+12lnΣ−N2ln2π.


### 3.4. Maximization of the Evidence Function with RVM

Firstly, we consider the maximizingptX,α,β with respect to αi. By making use of (33) and Lemma 3, we obtain(54)dlnΣdαi=TrΣ−1dΣdαi=TrA+βΦTΦdA+βΦTΦ−1dαi.From (54) and Lemma 4, we obtain(55)dlnΣdαi=−TrdA+βΦTΦdαiΣ,where we have used(56)dA+βΦTΦ−1dαi=−A+βΦTΦ−1dA+βΦTΦdαiA+βΦTΦ−1and (33).From (55), we obtain(57)dlnΣdαi=−Σiiwhere Σii is the ith component of the posterior covariance Σ defined by (33).From (42), we obtain(58)∂Em∂αi=12mi2.From (57) and (58), we easily obtain the stationary point of (53) with respect to αi such that(59)12αi−12Σii−12mi2=0,where mi is the ith component of the posterior mean m defined by (32).After multiplying both sides of (59) by 2αi and arranging, we have the following result:(60)αimi2=1−αiΣii,(61)αi=γimi2,together with γi=1−αiΣii.Similarly, we can maximize the log marginal likelihood (53) with respect to β. Firstly, we obtain the derivative of lnΣ with respect to β using (33) and Lemma 3:(62)dlnΣdβ=TrΣ−1dΣdβ=TrA+βΦTΦdA+βΦTΦ−1dβ.By making use of (33) and Lemma 4, we obtain(63)dlnΣdβ=−TrΦTΦΣ.SinceA+βΦTΦA+βΦTΦ−1=IM, we rearrange and get the following result:(64)ΦTΦΣ=β−1IM−AΣ.By using definition of matrixA, (63) and (64), we obtain(65)dlnΣdβ=−1β∑i1−αiΣii.From (42), we obtain(66)∂Em∂β=12t−Φm2.From (65) and (66), we obtain the stationary point of (53) with respect to β such that(67)Nβ−t−Φm2−1β∑iγi=0.By multiplying both sides of (67) by β and arranging, we obtain(68)1β=t−Φm2N−∑iγi.The result of (61) and (68) is consistent with [5, 8]. We note that (61) and (68) are implicit solutions for α and β since both γi and 1/β depend on α and β. To address the problem, an iterative procedure is devised. We use (32) and (33) to evaluate the mean and covariance of the posterior distribution by choosing initial value α and β, respectively, and then alternately re-evaluate the hyperparameters using (61) and (68). Again, learning continues to use (32) and (33) to re-evaluate the mean and covariance of the posterior distribution, until a suitable convergence standard is satisfied.


## 3.1. Linear Basis Function Models and Bayesian Linear Regression

The simplest linear model for regression is the one that includes a linear combination of the input variables:(13)yx,w=w0+w1x1+⋯+wDxD.However, a linear combination of the input variables imposes significant limitations when we use the model. So, we extend the model by using linear combinations of fixed nonlinear functions of the input variables, which take the form(14)yx,w=w0+∑i=1M−1wiϕix,where ϕix is defined as basis function and w0 is called a bias parameter. For convenience, we define an additional dummy “basis function” ϕ0x=1 and rewrite(15)yx,w=∑i=0M−1wiϕix=wTϕx,where w=w0,⋯,wM−1T and ϕx=ϕ0x,⋯,ϕM−1xT [8].The target variablet is given by a deterministic function yx,w with additive Gaussian noise such that(16)t=yx,w+ε,where ε is a zero mean Gaussian random variable with precision (inverse variance) β. So, we obtain(17)ptx,w,β=Ntyx,w,β−1.Now, we consider a data set of inputsX=x1,⋯,xN with corresponding target value t=t1,⋯,tN. From the assumption that these data points are drawn independently from the Gaussian distribution (17), we obtain the following likelihood function:(18)ptX,w,β=∏n=1NNtnwTϕxn,β−1.Since (18) is the exponential of a quadratic function of w, we can define a corresponding conjugate prior given by(19)pw=Nwm0,S0.From (18) and (19), we obtain the posterior distribution in the form(20)pwt,X,β=NwmN,SN,where(21)mN=SNS0−1m0+βΦTt,(22)SN−1=S0−1+βΦTΦ,where Φ is an N×M matrix, called the design matrix, whose elements are given by Φnj=ϕjxn so that(23)Φ=ϕ0x1…ϕM−1x1⋮⋱⋮ϕ0xN⋯ϕM−1xN.The proofs of (21) and (22) in detail refer to Appendix A.


## 3.2. Bayesian Model Comparison and Evidence Approximation

Assume that we face a set ofL models Mi, where i=1,⋯,L. Our task is to compare them and choose an optimal model based on some standard from the Bayesian perspective, which can avoid the overfitting phenomenon associated with maximum likelihood. The uncertainty is expressed through a prior probability distribution pMi. Given a training set, we suppose that all models are given the equal prior probability, which is feasible because there should be no preference for L given models in practice. We then wish to evaluate the posterior distribution:(24)pMiD∝pMipDMi,where pMiD is defined as model evidence, which is also called the marginal likelihood since it can be viewed as a likelihood function over the space of models, in which the parameters have been marginalized out.Now, we introduce prior distribution over the hyperparametersα and β to complete a fully Bayesian treatment of the linear basis function model and we can obtain prediction by marginalizing with respect to the hyperparameters α and β and the parameter w. However, the complete marginalization over all of these variables is analytically intractable. We then discuss an estimation, in which we determine the hyperparameters α and β by maximizing the marginal likelihood function obtained by first integrating over the parameters w. The framework is viewed in the statistics field as empirical Bayes [13, 14], or type-2 maximum likelihood [15], or generalized likelihood [16], and is also called the evidence approximation in the machine learning literature [8, 17, 18].Here, if we introduce hyperprior over theα and β, then we obtain predictive distribution by marginalizing over the hyperparameters α and β and the parameter w as follows:(25)ptt=∭ptx,w,βpwt,α,βpα,βtdwdαdβ,where ptx,w,β is defined by (17) and pwt,α,β is given by (20) with mN and SN defined by (21) and (22).If the posterior distributionpα,βt is sharply peaked around α^ and β^, we can obtain the predictive distribution simply by marginalizing over w in which α and β are set to the value α^ and β^:(26)ptt≃ptt,α^,β^=∫ptx,w,β^pwt,α^,β^dw.The posterior distribution forα and β is given by(27)pα,βt∝ptα,βpα,β.We can obtain the valuesα^ and β^ by maximizing the marginal likelihood function ptα,β in the evidence framework if the prior is relatively flat. We continue evaluating the marginal likelihood for the linear basis function and obtain its maxima which will allow us to determine the values of these hyperparameters from the training data set only.


## 3.3. Evaluation of the Evidence Function with RVM

Relevance vector machine (RVM) is a Bayesian sparse kernel method for regression. The relevance vector machine for regression is a linear model of the form given by (15). The following likelihood function is given by (18).The weight prior takes the form(28)pwα=∏n=1NNwi0,αi−1,where α=α1,⋯,αMT and αi is hyperparameter and represents the precision of the corresponding parameter wi.By making use of (18), we obtain(29)ptX,w,β=∏n=1NNtnwTϕxn,β−1=β2πN/2exp−β2∑i=1Ntn−wTϕxn2=β2πN/2exp−β2t−Φw2.Similarly, by making use of (28), we obtain(30)pwα=∏n=1NNwi0,αi−1=α1α2⋯αM2πM/2exp−12∑i=1Mαiwi2=A1/22πM/2exp−12wTAw.We obtain the posterior distribution for the weights, which is also Gaussian and takes the form(31)pwt,X,α,β=Nwm,Σ,where the covariance and mean are given by(32)m=βΣΦTt,(33)Σ=A+βΦTΦ−1,where A=diagαi.The proofs of (32) and (33) in detail refer to Appendix B.We can use evidence approximation, which is also known as type-2 maximum likelihood, to determine the values of hyperparametersα and β. We first obtain the marginal likelihood function by integrating out the weight parameters:(34)ptX,α,β=∫ptX,w,βpwαdw,and then maximize the marginal likelihood function.We substitute (32) and (33) into (34) and then obtain the following result:(35)ptX,α,β=β2πN/2A1/22πM/2∫exp−Ewdw,where we have defined(36)Ew=β2t−Φw2+12wTAw=β2t−ΦwTt−Φw+12wTAw=β2tTt−2wTΦTt+wTΦTΦw+12wTAw.For the stationary point ofEw, we first derive the gradient of ∇Ew by using Corollary 1:(37)∇Ew=βΦTΦ+Aw−βΦTt,(38)∇∇Ew=A+βΦTΦ=Σ−1.Let(39)∇Ew=βΦTΦ+Aw−βΦTt=0.From (39), we obtain(40)w=ββΦTΦ+A−1ΦTt=βΣΦTt=m.From (11), (38), and (40), it follows that(41)Ew=Em+12w−mTΣ−1w−m,together with(42)Em=β2t−Φm2+12mTAm.We substitute (41) into (35) and then get the following form:(43)ptX,α,β=β2πN/2A1/22πM/2∫exp−Emexp−12w−mTΣ−1w−mdw=β2πN/2A1/22πM/2exp−Em∫exp−12w−mTΣ−1w−mdw=β2πN/2A1/2Σ1/2exp−Em,where we have used 1/2πM/2Σ1/2∫exp−1/2w−mTΣ−1w−mdw=1.From (32) and (33), we obtain(44)βΣΦTt=Σ−1m.We substitute (33) into (44) and obtain(45)βΣΦTt=A+βΦTΦm.We expand (45) and then obtain(46)βΣΦTt=Am+βΦTΦm.We multiply both sides of (46) by mT and derive(47)βmTΣΦTt=mTAm+βmTΦTΦm.From (42), we obtain(48)Em=β2t−Φm2+12mTAm=12βt−ΦmTt−Φm+mTAm=12βtTt−2βmTΦTt+βmTΦTΦm+mTAm=β2tTIN−mTΦTt=12tTβIN−β2ΦΣΦTt,where we make use of (32) and (47).According to Lemma1, we obtain(49)βIN−β2ΦΣΦT−1=βIN−βΦA+βΦTΦ−1βΦT−1=β−1IN+ΦA−1ΦT.By making use of (48) and (49), we obtain(50)ptX,α,β=β2πN/2A1/2Σ1/2exp−Em=β2πN/2A1/2Σ1/2exp−12tTβ−1IN+ΦA−1ΦT−1t=12πN/2β−1IN+ΦA−1ΦT1/2exp−12tTβ−1IN+ΦA−1ΦT−1t,where we have used the following result according to Lemma 2 and (33):(51)β−1IN+ΦA−1ΦT=1βNIN+βΦA−1ΦT=1βNIM+βA−1ΦTΦ=1βNA−1A+βΦTΦ=1βNA−1Σ−1=1βNA−1Σ−1.From (50), we obtain(52)ptX,α,β=Nt0,C,together with C=β−1IN+ΦA−1ΦT.From (43), we can then write the log of the marginal likelihood in the form:(53)lnptX,α,β=N2lnβ−Em+12∑i=1Mlnαi+12lnΣ−N2ln2π.


## 3.4. Maximization of the Evidence Function with RVM

Firstly, we consider the maximizingptX,α,β with respect to αi. By making use of (33) and Lemma 3, we obtain(54)dlnΣdαi=TrΣ−1dΣdαi=TrA+βΦTΦdA+βΦTΦ−1dαi.From (54) and Lemma 4, we obtain(55)dlnΣdαi=−TrdA+βΦTΦdαiΣ,where we have used(56)dA+βΦTΦ−1dαi=−A+βΦTΦ−1dA+βΦTΦdαiA+βΦTΦ−1and (33).From (55), we obtain(57)dlnΣdαi=−Σiiwhere Σii is the ith component of the posterior covariance Σ defined by (33).From (42), we obtain(58)∂Em∂αi=12mi2.From (57) and (58), we easily obtain the stationary point of (53) with respect to αi such that(59)12αi−12Σii−12mi2=0,where mi is the ith component of the posterior mean m defined by (32).After multiplying both sides of (59) by 2αi and arranging, we have the following result:(60)αimi2=1−αiΣii,(61)αi=γimi2,together with γi=1−αiΣii.Similarly, we can maximize the log marginal likelihood (53) with respect to β. Firstly, we obtain the derivative of lnΣ with respect to β using (33) and Lemma 3:(62)dlnΣdβ=TrΣ−1dΣdβ=TrA+βΦTΦdA+βΦTΦ−1dβ.By making use of (33) and Lemma 4, we obtain(63)dlnΣdβ=−TrΦTΦΣ.SinceA+βΦTΦA+βΦTΦ−1=IM, we rearrange and get the following result:(64)ΦTΦΣ=β−1IM−AΣ.By using definition of matrixA, (63) and (64), we obtain(65)dlnΣdβ=−1β∑i1−αiΣii.From (42), we obtain(66)∂Em∂β=12t−Φm2.From (65) and (66), we obtain the stationary point of (53) with respect to β such that(67)Nβ−t−Φm2−1β∑iγi=0.By multiplying both sides of (67) by β and arranging, we obtain(68)1β=t−Φm2N−∑iγi.The result of (61) and (68) is consistent with [5, 8]. We note that (61) and (68) are implicit solutions for α and β since both γi and 1/β depend on α and β. To address the problem, an iterative procedure is devised. We use (32) and (33) to evaluate the mean and covariance of the posterior distribution by choosing initial value α and β, respectively, and then alternately re-evaluate the hyperparameters using (61) and (68). Again, learning continues to use (32) and (33) to re-evaluate the mean and covariance of the posterior distribution, until a suitable convergence standard is satisfied.


## 4. EM Algorithm for RVM

A powerful and elegant technique for finding likelihood solutions for probabilistic models with latent variables is called the expectation-maximization algorithm or EM algorithm which also forms the basis for derivation of the variational framework [19, 20].The integration with respect to the weights to obtain the marginal likelihood can make use of the EM algorithm for maximization with the weights as latent variables very naturally. Since it is most common to use the EM for maximization,L will now become the (positive) log evidence. To derive the EM algorithm for the RVM, we firstly obtain the log marginal likelihood by marginalizing the weights from the joint distribution over t and w:(69)L=lnptX,α,β=ln∫pt,wX,α,βdw,where the joint distribution over t and w is both equal to the marginal likelihood times the posterior on w and to the likelihood multiplied by the prior on w:(70)pt,wX,α,β=ptX,w,βpwα=pwt,X,α,βptX,α,β.By defining a variational probability distributionqw over the weights and using Jensen’s inequality, we obtain a lower bound on L:(71)L=ln∫qwpt,wX,α,βqwdw≥∫qwlnpt,wX,α,βqwdw=Fq,α,β.The EM algorithm achieves the maximization of the log marginal likelihoodL by iteratively maximizing the lower bound. In the Expectation step (E-step), F is maximized with respect to the variational distribution qw for fixed parameters α and β, and in the Maximization step (M-step), F is maximized with respect to the hyperparameters α and β for fixed qw. Insight into how to perform the E-step can be gained by rewriting the lower bound F as(72)F=L+KLqwpwt,X,α,β,where KLqwpwt,X,α,β is the Kullback–Leibler (KL) divergence between the variational distribution qw and the posterior distribution of the weights. The KL divergence is always a positive number and is only equal to zero if the two distributions are identical. The E-step corresponds thus to making qw equal to the posterior on w, which implies F=L and qw=pw|t,X,α,β. Since the posterior is Gaussian, the E-step reduces to computing its mean m and covariance Σ defined by (32) and (33).To perform the M-step, we rewriteF in a different way:(73)F=∫qwlnpt,wX,α,βdw−Hqw,where Hqw is entropy of qw, which is an irrelevant constant because it does not depend on α or β, which implies ∂Hqw/∂αi=0 and ∂Hqw/∂β=0. Therefore, the M-step is to maximize the average of the log joint distribution of t and w over qw with respect to the parameters α and β:(74)F=∫qwlnpt,wX,α,βdw−Hqw=∫qwlnpwα+lnptX,w,βdw−Hqw=∫qwlnpwαdw+∫qwlnptX,w,βdw−Hqw.From the first term of (74), we obtain(75)∫qwlnpwαdw=∫qw12lnA−M2ln2π−12wTAwdw=∫qw12lnA−M2ln2πdw−12∫qwwTAwdw=12lnA−M2ln2π−12∫qwwTAwdw,where we have used ∫qww=1.From the third term of (75) and Lemma 5, we obtain(76)∫qwwTAwdw=∫qwTrAwwTdw=Tr∫qwAwwTdw=TrA∫qwwwTdw=TrAEqwwT=TrAmmT+Σ=TrAmmT+TrAΣ=mTAm+TrAΣ,where Eq⋅ denotes the expectation of a random variable over the distribution qw and we have used EqwwT=mmT+Σ, where the mean m and covariance Σ is defined by (32) and (33).From (75) and (76), we obtain(77)∫qwlnpwαdw=12lnA−M2ln2π−12∫qwwTAwdw=12lnA−M2ln2π−12mTAm−12TrAΣ.We start to calculate the second term of (74) and obtain(78)∫qwlnptX,w,βdw=∫qwlnβ2πN2exp−β2t−Φw2dw=∫qwN2lnβ−ln2π−β2t−Φw2dw=∫qwN2lnβ−N2ln2πdw−β2∫qwt−Φw2dw=N2lnβ−N2ln2π−β2∫qwt−Φw2dw=N2lnβ−N2ln2π−β2∫qwt−ΦwTt−Φwdw=N2lnβ−N2ln2π−β2∫qwtTt−2wTΦTt+wTΦTΦwdw.From the third term of (78), we obtain(79)∫qwtTt−2wTΦTt+wTΦTΦwdw=∫qwtTtdw−2∫qwwTΦTtdw+∫qwwTΦTΦwdw=tTt−2∫qwwTΦTtdw+∫qwTrΦwwTΦTdw=tTt−2∫qwwTΦTtdw+Tr∫qwΦwwTΦTdw=tTt−2∫qwwTdwΦTt+TrΦ∫qwwwTdwΦT=tTt−2EqwTΦTt+TrΦEqwwTΦT=tTt−2mTΦTt+TrΦmmT+ΣΦT=tTt−2mTΦTt+TrΦmmTΦT+ΦΣΦT=tTt−2mTΦTt+TrΦmmTΦT+TrΦΣΦT=tTt−2mTΦTt+mTΦTΦm+TrΦΣΦT=t−ΦmTt−Φm+TrΦΣΦT=t−Φm2+TrΦΣΦT=t−Φm2+TrΣΦTΦ=t−Φm2+TrΦTΦΣ,where we have used Lemma 5 and EqwwT=mmT+Σ.From (78) and (79), we obtain(80)∫qwlnptX,w,βdw=N2lnβ−N2ln2π−β2t−Φm2+TrΦTΦΣ.From (77) and (80), we obtain(81)F=∫qwlnpt,wX,α,βdw−Hqw=12lnA−M+N2ln2π−12mTAm−12TrAΣ+N2lnβ−β2t−Φm2+TrΦTΦΣ−Hqw.We easily obtain the stationary point of (81) with respect to αi and β such that(82)∂F∂αi=12αi−12mi2+Σii=0,(83)∂F∂β=N2β−12t−Φm2+TrΦTΦΣ=0.So, we obtain the update rules as follows:(84)αinew=1mi2+Σii,(85)βnew=Nt−Φm2+TrΦTΦΣ,where mi is the ith component of the posterior mean m defined by (32) and Σii is the ith component of the posterior covariance Σ defined by (33).By making use of (64) and (65), we further obtain the following result for (85):(86)βnew=Nt−Φm2+βold−1∑iγi,where γi=1−αiΣii.From the abovementioned discussion, we obtain the update rules by using the EM algorithm to maximize the log marginal likelihood, which is also consistent with Quiñonero Candela’s result, except that we use the precision matrix as the covariance [21].


## 5. Conclusions

From the abovementioned discussion, we can draw the conclusion that the novel logical framework presented in this paper based on Taylor expansion shows some advantages compared with the conventional method of completing the square, which is easier to be implemented due to the clear logical framework, and avoids the limitation in seeking the term of completing the square intentionally. The approximation and maximization of the evidence function associated with RVM is consistent with Tipping’s result of (61) and (68) using the novel logical framework. We obtain the update rules of (84) and (85) by using the EM algorithm for maximization of the log marginal likelihood, which is also consistent with Quiñonero Candela’s result, except that we use the precision matrix as the covariance.



---

*Source: 2548310-2020-03-28.xml*