A data augmentation approach for a class of statistical inference problems

We present an algorithm for a class of statistical inference problems. The main idea is to reformulate the inference problem as an optimization procedure, based on the generation of surrogate (auxiliary) functions. This approach is motivated by the MM algorithm, combined with the systematic and iterative structure of the Expectation-Maximization algorithm. The resulting algorithm can deal with hidden variables in Maximum Likelihood and Maximum a Posteriori estimation problems, Instrumental Variables, Regularized Optimization and Constrained Optimization problems. The advantage of the proposed algorithm is to provide a systematic procedure to build surrogate functions for a class of problems where hidden variables are usually involved. Numerical examples show the benefits of the proposed approach.


Introduction
Problems in statistics and system identification often involve variables for which measurements are not available. Among others, real-life examples can be found in communication systems [1,2] and systems with quantized data [3,4]. In Maximum Likelihood (ML) estimation problems, the likelihood function is in general difficult to optimize by using closed-form expressions, and numerical approximations are usually cumbersome. These difficulties are traditionally avoided by the utilization of the Expectation-Maximization (EM) algorithm [5], where a surrogate (auxiliary) function is optimized instead of the main objective function. This surrogate function includes the complete data, i.e. the measurements and the random variables for which there are no measurements. The incorporation of such hidden data or latent variables is usually termed as data augmentation, where the main goal is to obtain, in general, simple and fast algorithms [6].
On the other hand, the MM (MM stands for Maximization-Minorization or Minimization-Majorization, depending on the optimization problem that needs to be solved) algorithm PLOS ONE | https://doi.org/10.1371/journal.pone.0208499 December 10, 2018 1 / 24 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 [7] is generally employed for solving more general optimization problems, not only for ML and Maximum a Posteriori (MAP) estimation problems. In general, the main motivation for using the MM algorithm is the lack of closed-form expressions for the solution of the optimization problem or dealing with objective cost functions that are not convex. Applications where the MM algorithm has been utilized include communication systems problems [8] and image processing [9]. For constrained optimization problems, an elegant solution is presented by Marks and Wright [10], where the constraints are incorporated via the formulation of surrogate functions. Surprisingly, Marks' approach has not received the same attention from the scientific community when it comes to compare it with the EM and the MM algorithms. In fact, these three approaches are contemporary, but the EM algorithm has attracted most of the attention (out of the three methods), and it has been used for solving linear and nonlinear statistical inference problems in biology and engineering, see e.g. [11][12][13][14][15][16], amongst others. On the other hand, as shown in [7], the MM algorithm has obtained much less attention, while Marks' approach is mostly known to a limited audience in the the Communication Systems community. These three approaches have important similarities: i) a surrogate function is defined and optimized in place of the original optimization problem, and ii) the solution is obtained iteratively. In general, these algorithms are "principles and recipes" [17] or a "philosophy" [7] for constructing solutions to a broad variety of optimization problems.
In this paper we adopt the ideas behind [5,7,10] to develop an algorithm for a special class of functions. Our approach generalizes the ones of [5,7,10] by reinterpreting the E-step in the EM algorithm and expressing the cost function in terms of an infinite mixture or kernel. This kind of problems can be interpreted as inverse problems that, in turn, can be posed as integral equations, such as channel modelling in wireless communications [18], estimation of the distribution of stellar rotational velocities [19], and mass estimation in particle physics problems [20,21]. In particular cases, the kernel corresponds to a variance-mean Gaussian mixture (VMGM), see e.g. [22]. VMGMs, also referred to as normal variance-mean mixtures [23] and normal scale mixtures [24], have been considered in the literature for formulating EM-based approaches to solve ML [25] and MAP problems, including regularized sparse estimation problems [22,26,27]. Sparse estimation problems have been widely studied in the last two decades and several techniques have been developed that include different types of penalties or contraint, see e.g [28,29] and different strategies for the formulation of those penalties/constraints [30]. In this paper we show that our proposal can also be considered for sparsity problems, however the analysis of the solution is out of the scope of the paper. Our approach is applicable to a wide class of functions, which allows for defining the likelihood function, the prior density function, and constraints as kernels, extending also the work in [10]. Thus, our work encompasses the following contributions: i) a systematic approach to constructing surrogate functions for a class of cost functions and constraints, ii) a class of kernels where the unknown quantities of the algorithm can be easily computed, and iii) a generalization of [5,7,10] by including the cost function and the constraints in one general expression. Our proposal is based, among other things, on a particular way to apply Jensen's inequality [31]. In addition, we provide the details on how to construct quadratic surrogate functions for cost functions and constraints.
Our algorithm is tested by two examples. In the first one we considered the problem of estimating the rotational velocities of stars. The system model corresponds to the convolution of two probability density functions (pdf's) and thus it is an infinite mixture. We show that our reinterpretation of the EM algorithm allows for the direct application of our proposal for the correct estimation of the parameter of interest. In the second example, we considered the estimation of the channel in wireless communications, where the true distribution can be either Rayleigh or Rice, depending on environment where the electromagnetic waves propagate. The problem is solved considering a sum of a Rayleigh and a Rice term, allowing for a more complex channel distribution. To select the more representative distribution, Akaike's Information Criterion [32] was also considered in order to obtain the least complex model that exhibits the best possible fitting.

The EM algorithm
Let us consider an estimation problem and its corresponding log-likelihood function defined as ℓ(θ) = log p(y|θ), where p(y|θ) is the likelihood function, θ 2 R p , and y 2 R N . Denoting the complete data by z 2 O(y), and using Bayes' theorem, we can obtain: Using Jensen's inequality [31], it is possible to show that for any value of θ, the function Hðθ;θ ðiÞ Þ is decreasing. Hence, the optimization is only carried out on the auxiliary function Qðθ;θ ðiÞ Þ because, by maximizing Qðθ;θ ðiÞ Þ, the new parameterθ ðiþ1Þ is such that the likelihood function increases (see e.g. [5,33]

Data augmentation in inference problems
Data augmentation algorithms are based on the construction of the augmented data and its many-to-one mapping O(y). This augmented data is assumed to describe a model from which the observed data y is obtained via marginalization [36]. That is, a system with a likelihood function p(y|θ) can be understood to arise from where the augmented data corresponds to (y, x) and x is the latent data [6,36]. This idea has been utilized for supervised learning [37] and the development of the Bayesian Lasso [38], to mention a few examples. In those problems, the Laplace distribution is expressed as a twolevel hierarchical-Bayes model. This equivalence is obtained from the representation of the Laplace distribution as a VMGM: In fact, there are several pdf's than can be expressed as VMGMs, as shown in Table 1 [22], where g(θ) is the penalty term that can be expressed as a pdf. In addition, in [26] it was developed an early version of the methodology presented in this paper, exploring the estimation of a sparse parameter vector utilizing the ℓ q -(pseudo)norm, with 0 < q < 1.

A systematic approach to construct surrogate functions for a class of inference problems
Here, we consider a general optimization cost defined as: where θ is a parameter vector, y is a given data (i.e. measurements), z is the complete data (comprised of the observed data y and the hidden variables (unobserved data), O(y) is a mapping from the sample space of z to the sample space of y, K(�, �) is a (positive) kernel function, and μ(�) is a measure, see e.g [31]. The definition in (9) is based on the definition of the auxiliary function Q in the EM algorithm [5], where it is assumed throughout the paper that there is a mapping that relates the not observed data to the observed data, and that the complete data lies in O(y) [5]. Notice that in (9) the kernel function may not be a pdf. However, several functions can be expressed in terms of a pdf. The most common cases are Gaussian kernels (yielding VMGMs) [23] and Laplace kernels (yielding Laplace mixtures) [39]. Remark 2. Notice that, as explained in Section 1, once the hidden data has been selected, the data augmentation procedure comes with the definition of VðθÞ in (9). From here, we follow the systematic nature of the EM and MM algorithms in terms of the iterative nature of the technique.

Constructing the surrogate function
Since we are considering the optimization of the function VðθÞ, we can also consider the optimization of the function Without modifying the cost function in (10), we can multiply and divide by the logarithm of the kernel function, obtaining: Hence, for any value of θ, the function Hðθ;θ ðiÞ Þ in (14) is a decreasing function. Remark 3. The kernel function K(z, θ) satisfies the standing assumption K(z, θ) > 0 since the proposed scheme is built, among others, on the logarithm of the kernel function K(z, θ). The definition of the kernel function in (9) allows for kernels that are not pdf's. On the other hand, some kernels may correspond to a scaled version of a pdf. In that sense, for the cost function in (9) we can define a new kernel and a new measure as The idea behind using a surrogate function is to obtain a simpler algorithm for the optimization of the objective function when compared to the original optimization problem. This can be achieved iteratively if the Fisher Identity for the surrogate function and the objective function is satisfied. That is, Lemma 1. For the class of objective functions in (9), the surrogate function Qðθ;θ ðiÞ Þ in (13) satisfies the Fisher identity defined in (17).
Proof. From (12) we have:  (13) can be utilized to obtain an adequate surrogate function that satisfies the properties in Mark's approach in (38)- (40).
We summarize our proposed algorithm in Table 2.

Surrogate functions for inverse problems
The objective function in (9) can be understood as an integral equation [40], where the unknown function of the integral equation corresponds to the kernel K(z, θ). In this kind of problems, samples from VðθÞ are available, whilst dμ(z) is assumed known. Several problems can be posed as this kind of problems. In the following, we explain how to use the approach presented in this paper to solve different inverse problems that arise in the integral equation form.

Stellar rotational velocity estimation.
One of the many problems in Astronomy deals with is the estimation of rotational velocities of stars. This particular problem is of great importance, since it allows astronomers to describe and model the stars formation, their internal structure and evolution, as well as how they interact with other stars, see e.g. [19,41,42].
Modern telescopes allow for the measurement of the rotational velocities from the telescope point of view, that is, a projection of the true rotational velocity. This is modelled (spatially) as the convolution of the true rotational velocity pdf and a uniform distribution over the sphere (for more details see e.g. [19]): ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where p(y|σ) is the uniform projected rotational velocity pdf and p(x|σ) is the true rotational velocity pdf to be estimated, and σ a hyperparameter. Thus, we can define the kernel function as K(x, σ) = p(x|σ) and dmðxÞ ¼ y x ffi ffi ffi ffi ffi ffi ffi ffi dx. This definition allows for the direct utilization of the expressions in (13) in order to estimate the parameter σ that defines the unknown rotational velocity pdf. We illustrate with one example in Section 6.1.

Channel estimation in wireless communications.
It is well known in the Communications community that the wireless channel corresponds to the superposition of different copies of the transmitted signal that have been reflected, refracted and scattered. Thus, those copies arrive at the receiver with different phase. Those components are referred to as multipath components [43]. On the other hand, it has been shown empirically that a good model for the multipath channel corresponds to either Rayleigh or Rice distributions [44]. However, there are cases when those distributions do not provide a good model for the channel. One of those cases corresponds to the presense of different channel models in the vicinity of the transmitter and the vicinity of the receiver. This situation ocurrs particularly in the so called urban scenarios, where the channel exhibits different behaviouirs in different places. For this scenario, an adequate model that takes into consideration different models and a transition from a local distribution and a global distribution is a continuous mixture of the form [18] pðxÞ ¼ where K(x, σ) is a pdf in a local area, μ(σ) is a distribution of σ which depends on the path Table 2. Proposed algorithm.
Step 7: i = i + 1 and back to Step 4 until convergence. from transmitter to the local cluster, and σ is a vector in the parameter space. Our approach can be directly used in order to estimate the true nature of the channel when expressed as a mixture. This can be done since the four most common chanel distributions are Rayleigh, Rice, log-normal and Nakagami-m [45]. Hence, assuming that the local and global distribution families are known, the attainment of the auxiliary function Qðσ;σ ðiÞ Þ is straightforward and thus the ML estimate of σ.

Neutrino mass search in particle physics.
In the Particle Physics community there is a plethora of works dealing with the estimation of masses of neutrinos, see e.g. [20] and the references therein. Among the methods that are generally used to determine the absolute masses of neutrinos we find the β-decay and the direct determination of neutrino mass, see e.g [21]. In β-decay methods, the neutrinos are analized based on their energy spectrums, where the measurements, corresponding to the observed β-spectra, are associated with an integral equation of the form [20,21] FðUÞ ¼ where b is a constant that represents the measurement noise, R(E) is the emitted β-spectra, and T 0 (E, U) is the impulse response of the equipment. Again, by regarding T 0 (E, U) as the kernel function and R(E)dE as a measure, the attainment of the auxiliary function Q is straightforward.

Estimation of mixture distributions.
Mixture distributions have been widely studied in the literature, particularly finite mixtures, see e.g. [46,47]. Their representation can be expressed in a general fashion by the notation [48] pðyjθÞ ¼ where K(z, θ) is a suitable function that may be either a pdf (for continuous random variables) or a probability function (for discrete random variables). The expression in (22) for an infinite mixture. Our approach can be tailored for the estimation of the parameters of a finite mixture of the form where we have assumed that pðyjβÞ ¼ Q N k¼1 pðy k jβÞ, pðyjβÞ ¼ P M j¼1 l j � j ðy k ; y j Þ, and λ = [λ 1 ,. . ., λ M ] are the mixing weights, ϕ j (y, θ j ) are the components densities parametrised by θ j . The kernel function defined in our approach can be utilized to represent jth component in the discrete mixture in (25) as where the β = [β 1 ,. . ., β M ] is the vector of parameters to be estimated, with β j = [λ j , θ j ]. Notice that the dependence with respect to the variable z is implicit. Utilizing the expression derived in (13) we obtain the following E-step Notice that, as shown in [49], the auxiliary function Qðβ;β ðiÞ Þ in (27) is the same that we obtain when we consider that the data are fully categorized, i.e. each y k , k = 1,. . ., N, is assumed to be drawn from only one distribution of the mixture. This assumption yields a data augmentation problem that is solved using the EM algorithm [46].

Remark 7. If we consider a combination of infinite mixtures and finite mixtures of the form
with pðx k jβÞ ¼ we can utlize the same approach described here to solve the problem of estimating the parameters in (28), β. In this case, the jth kernel is defined as [50] K j ðx k ; β j Þ ¼ l j �ðx k ; y j Þ; ð30Þ and measure Then, the log-likelihood function can be expressed as

Constrained problems in statistical inference
Statistical Inference and System Identification techniques include a variety of methods that can be used in order to obtain a model of a system from data. Classical methods, such as Least Squares, ML, MAP [51], Prediction Error Method, Instrumental Variables [52], and Stochastic Embedding [53] have been considered in the literature for such task. However, the increasing complexity of modern system models has motivated researchers to revisit and reconsider those techniques for some problems. This has resulted in the incorporation of constraints and penalties, yielding an often more complicated optimization problem. For instance, it has been shown that the incorporation of linear equality constraints may improve the accuracy of the parameter estimates, see e.g [54]. On the other hand, the incorporation of regularization terms (or penalties) also improves the accuracy of the estimates, reducing the effect of noise and eliminating spurious local minima [55]. Regularization can be mainly incorporated in two ways: by adding regularizing constraints (a penalty function) or by including a probability density function (pdf) as a prior distribution for the parameters, see e.g. [27]. Another way to improve the estimation is by incorporating inequality constraints, where certain functions of the parameters may be required, for physical reasons amongst others, to lie between certain bounds [56]. From this point of view, it is possible to consider the classical methods with constraints or penalties, as in [53,[55][56][57]. Perhaps one of the most utilized approaches for penalized estimation (with complicated non-linear expressions) is the MM algorithm-for details on the MM algorithm see Section 2.2. This technique allows for the utilization of a surrogate function that is simple to handle, in terms of derivatives and optimization techniques, and that is, in turn, iteratively solved. However, its inequality constraint counterpart, here referred to as Marks' approach [10], is not as well known as the MM algorithm. Moreover, there is no straightforward manner to obtain such surrogate function. In this paper we focus on a systematic way to obtain the corresponding surrogate function using Marks' approach for a class of constraints.

Mark's approach
The approach in [10] deals with inequality constraints by using a similar approach to the EM and MM algorithms. The basic idea is, again, to generate a surrogate function that allows for an iterative procedure whose optimum value is the optimum value of the original optimization problem.
Let us consider the following constrained optimization problem: Provided the above properties are satisfied, then the following approximation of (37): ð41Þ iteratively converges to the solution of the optimization problem (37). As shown in [10], the optimization problem in (41) is equivalent to the original problem (37), since the solution of (41) converges to a point that satisfies the Karush-Kuhn-Tucker conditions of the original optimization problem.

Remark 8. Mark's approach can be considered as a generalization of the MM algorithm, since the latter can be derived (for a broad class of problems) from the former. Let us consider the following problem:
Using the epigraph representation of (42) [58], we obtain the equivalent problem Using Mark's approach (41), we can iteratively find a local optimum of (42) via where Qðθ;θ ðiÞ Þ in (44) is a surrogate function for f(θ) in (42). From the epigraph representation we then obtain which is the definition of the MM algorithm (see 2.2) for more details.

A quadratic surrogate function for a class of kernels
In this section we focus on a special class of the kernel functions K(z, θ). For this particular class, the following is satisfied: where A(z) is a matrix and b is a vector, both of adequate dimensions.  (17) we have that from which we can solve for R in some cases. In other cases, the matrix R can also be computed using Monte Carlo algorithms. In particular, if A(z) is a diagonal matrix, then R is also a diagonal matrix defined by R = diag[r 1 , r 2 ,. . .]. Thus, we have where θ i is the ith component of the parameter vector θ,θ ðiÞ i is the ith component of the estimateθ ðiÞ , r i is the ith element of the diagonal of R, and b i is the ith element of the vector b.
Hence, when optimizing the auxiliary function Qðθ;θ ðiÞ Þ we obtain @ @θ Qðθ;θ ðiÞ Þ ¼ . . Equivalently,θ This implies that in our approach, it is not necessary to obtain the auxiliary function Q and optimize it. Instead, by computing R and b at every iteration, the new estimate can be obtained. [59] or quadrature methods [60].

Remark 10. The computation of the matrix R can be cumbersome when the matrix A(z) is not diagonal. In those cases, the integral that defines R can be computed utilizing Markov Chain Monte Carlo, quasi-Monte Carlo
The class defined in (46) arises naturally when dealing with VMGM, because the corresponding kernel function is normal and, thus, its logarithm is a quadratic function.
In particular, the utilization of VMGM encompasses different expressions commonly used for parameter estimation. Indeed, we have: (i) Lasso: The Lasso [61], expressed as a Laplace pdf, is represented by a VMGM [38], with p(θ|z) a zero-mean Gaussian distribution (of iid terms) as the kernel and p(z) an exponential distribution with parameter γ 2 /2. That is, (ii) Elastic-Net: The Elastic-Net penalty [62] is interpreted as a pdf if it corresponds to the product of two pdf's, a Laplacian (as in the Lasso case) and a Gaussian pdf. In this sense, by grouping those pdf's, we obtain [63] pðθÞ ¼ k EN with pðl j Þ / ffi ffi ffi ffi ffi ffiffi ð1À kÞ 2 k l j .
(iii) Group-Lasso: The Group-Lasso penalty is obtained via VMGM-representation as where w 2 l is the chi-squared distribution with l degrees of freedom. For the class of kernels here described, the proposed method for constructing surrogate functions can also be understood as part of sequential quadratic programming (SQP) methods [64] when, for instance, the above penalties are utilized as constraints in a constrained ML estimation problem. Indeed, the general case of equality and inequality-constrained minimization problems is defined as [65]: which is solved by iteratively defining quadratic functions that approximate the objective function and the inequality constraint around a current iterateθ ðiÞ . In the same way, our proposal generates an algorithm with quadratic surrogate functions, where an auxiliary functioñ Qðθ;θ ðiÞ Þ in (18) must be constructed for f(θ) and/or g(θ) in (53).

Numerical examples
In this section, we illustrate our proposed algorithm with two numerical examples.

Example 1: Estimation of the distribution of stellar rotational velocities
A common model for p(x|σ) found in the Astronomy literature is the Maxwellian distribution (see e.g. [19,66]) In practice, the measurements correspond to realizations of p(y|σ) [19], from which the likelihood function can be defined as: where y = [y 1 ,. . ., y N ] T , ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi x k is Maxwellian distributed, and N is the number of measurement points. Hence, the log-likelihood function can be expressed as: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi If we define the complete data z = (x, y), the kernel function K(�, �) and the measure μ(�) in (9) can be defined as and dmðx k ; y k Þ ¼ y k x k ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Then, the log-likelihood function in (56) can be written as: Thus, the ML estimator is obtained from: Since the parameter that is needed to be estimated is part of the convolution in (19), the optimization problem in (61) cannot be solved in a straightforward manner. Instead, we utilize the re-interpretation of the EM algorithm that we propose for solving (61 In Table 3 we summarized the specialisation of our proposed algorithm for this example. For the numerical simulation, we have considered the problem solved in [19], with the true dispersion parameter σ 0 = 8. The measurement data y = [y 1 ,. . ., y N ] was generated using the Slice Sampler (see e.g. [67]) applied to (19). The simulation setup is as follows: • The data length is given by N = 10000.
• The number of Monte Carlo (MC) simulations is 50.
• The stopping criterion is given by: kŝ ðiÞ Àŝ ðiÀ 1Þ k=kŝ ðiÞ k < 10 À 6 ; or the maximum number of iterations of 100 has been reached.
The results are shown in Fig 1, were the estimated p(x|σ) for each MC simulation is shown. It is clear that the estimated Maxwellian distributions are very similar to the true density distribution. The mean value of the estimated parameter wasŝ ¼ 7:9920. The estimation from each MC simulation is shown in Fig 2. It can be clearly seen that the estimated parameterŝ is close to the true value.

Example 2: Channel estimation in wireless communications
When modelling the wireless channel, a popular technique that is commonly used corresponds to the transmission of a sine tone at a given frequency, see e.g. [68]. The received power is then modelled as a random variable. The corresponding distribution has been widely studied in the literature from measurements, and the empirical data have shown that the two most common distributions are Rayleigh an Rice [44,45]. Hence, using similar ideas as in [50], in this example we formulate the channel distribution as a discrete sum based on a Rayleigh and a Rice component to determine the nature of the wireless channel. Table 3. Proposed algorithm for Maxwellian distribution estimation in Example 1.
Step 4: Computeŝ ðiþ1Þ (65) Step 5: i = i + 1 and back to Step 3 until convergence. https://doi.org/10.1371/journal.pone.0208499.t003 Data augmentation for statistical inference problems  First, the discrete mixture that we want to adjust from data is given by and I 0 (�) is the modified Bessel function of zeroth order. In addition, we must also include the constraint λ 1 + λ 2 = 1 so p(x) is a pdf. Thus, we can directly apply our proposed approach by assuming that each measurement point can be associated with a hidden variable that describes if such data point comes from the Rayleigh component or the Rice component, as it is traditionally formulated when dealing with discrete mixtures [46]. Hence, the auxiliary function Qðθ;θ ðiÞ Þ is given by whereθ ðiÞ is the current estimate, z is the unobserved (hidden) data and z tj is an indicator parameter such that z tj = 1 if the t-th observation comes from component j and is zero otherwise. It is given by: The estimateθ ðiþ1Þ at the next iteration is then given by: We also consider the utilization of Akaike's Information Criterion (AIC) in order to obtain an accurate yet simple model and, thus, discriminating from a Rayleigh channel, a Rice channel, and a mixture of both.
With the above formulation, we consider two cases: a Rayleigh distributed channel and a Rice distributed channel.
6.2.1 Rayleigh distributed channel. In this example, the random variable x is drawn from the Rayleigh distribution with σ 2 = 1, using the Slice Sampler [69]. The best estimated model corresponds to the single Rayleigh component in the mixture. The corresponding estimation of σ 2 yieldŝ s 2 1 ¼ 1:0007 � 9:003 � 10 À 4 . Fig 3 shows the true Rayleigh distribution and the mean estimated pdf from the 50 MC simulations. We observe an important agreement between the true pdf and the estimated model. with v = 4 and σ 2 = 1, using the Slice Sampler. The best model is selected as a single Rician component. The corresponding estimated parameters arev 2 ¼ 4:003 � 1:3 � 10 À 3 and s 2 2 ¼ 0:9858 � 2:2 � 10 À 3 . In Fig 4 we show the true Rice distribution and the mean estimated pdf from 50 MC simulations. We can observe that the estimator exhibits a good performance for the estimation of a Rice distribution.

Conclusions and future work
In this paper we have presented a systematic approach for constructing surrogate functions in a wide range of inference problems. Our approach can be utilized for constructing surrogate functions for both the cost function and the constraints, generalizing the popular EM and MM algorithms. Our approach is based on the utilization of data augmentation and kernel functions, yielding simple optimization algorithms when the kernel can be expressed as VMGM. We have shown how our proposal can be utilized to solve inverse problems that are expressed as integral equations and mixture distributions.
In addition, we have shown that our approach can be utilised for constrained/penalized ML and MAP estimations problems. In particular, common problems in statistical inference can directly be solved using our proposal since they can be posed as Variance Mean Gaussian Mixtures (VMGM), yielding quadratic surrogate functions.
In the last two decades the problem of sparse estimation has attracted a lot of attention. Since our approach can be utilized in those problems, and since it is based on the principles of the MM algorithm, a detailed analysis can be done in terms of accuracy and convergence of our technique, and compared against other techniques, such as the ones in [28,30], and [29], where different Lasso-type problems are compared, the MM algorithm is utilized in constrained problems for ML estimation in generalized linear model regression, and the MM algorithm is used for (unconstrained) sparse estimation under non-convex penalties, respectively.