{'fraction_non_alphanumeric': 0.04144197500077928, 'fraction_numerical': 0.013247716717059942, 'mean_word_length': 4.293103448275862, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Model Building for Semiparametric Mixtures

An important and yet difficult problem in fitting multivariate mixture models is determining the mixture complexity. We develop theory and a unified framework for finding the nonparametric maximum likelihood estimator of a multivariate mixing distribution and consequently estimating the mixture complexity. Multivariate mixtures provide a flexible approach to fitting high-dimensional data while offering data reduction through the number, location and shape of the component densities. The central principle of our method is to cast the mixture maximization problem in the concave optimization framework with finitely many linear inequality constraints and turn it into an unconstrained problem using a"penalty function". We establish the existence of parameter estimators and prove the convergence properties of the proposed algorithms. The role of a"sieve parameter'' in reducing the dimensionality of mixture models is demonstrated. We derive analytical machinery for building a collection of semiparametric mixture models, including the multivariate case, via the sieve parameter. The performance of the methods are shown with applications to several data sets including the cdc15 cell-cycle yeast microarray data.


Introduction
Multivariate mixture modeling is a bridge between clustering and nonparametric multivariate density estimation. The estimated multivariate mixture model provides both an estimate of the density for the overall data and partitions the data into several components or clusters. Determining the mixture complexity is challenging even in one dimension (Laird, 1978;Jewell, 1982;Titterington et al., 1985;Lesperance and Kalbfleisch, 1992;Roeder, 1994;Lindsay, 1995;McLachlan and Peel, 2001;Pilla and Loader, 2003;Scott, 2004a;Pilla and Charnigo, 2005). A fundamental problem in high-dimensional modeling is determining the number of components and their centers. One popular model-free approach to high-dimensional modeling is the K-means algorithm (see Hastie et al. (2001) and the references therein). Model-based techniques such as density estimation (Scott, 1992;James et al., 2001;Scott, 2004b) and multivariate mixture models provide a reliable and flexible approach to high-dimensional modeling while providing a data reduction through the number, location and shape of the component densities. In the context of mixture models, the problem becomes determining the number of mixture components and estimating the corresponding location parameter vectors. Furthermore, mixture models provide much of the flexibility of the nonparametric approaches, while retaining many advantages of the parametric approaches (Laird, 1978;Roeder, 1992;Lesperance and Kalbfleisch, 1992;Lindsay, 1995;Charnigo and Pilla, 2005;Scott, 2004b).
One of the main reasons for the popularity of model-free approaches such as the K-means algorithm for high-dimensional modeling is the lack of a unified and powerful technique for fitting multivariate mixtures. The focus of this article is to develop analytical machinery for building a collection of semiparametric mixture models, including the multivariate case. The theory and methods developed in this article have applications to image analysis, high-dimensional clustering and data mining, to name a few. The practical applications of semiparametric mixture models are broad and include case-control studies with errors-in-variables (Roeder et al., 1996), random effects models and empirical Bayes method (Lindsay, 1995). A natural outcome of applying a multivariate mixture model for high-dimensional clustering is that (1) each cluster is statistically represented by a para-metric distribution; for instance, normal in continuous case and Poisson in discrete case and (2) it provides the proportion of observations in each cluster through the estimated mixture probability. Furthermore, statistical tests can be developed easily based on the parameter estimators of the multivariate mixture models to answer various scientific or biological questions. Due to high levels of noise inherent in many of the massive data sets, including the microarray technology, it is highly desirable to carry out the high-dimensional data analysis within a statistical framework.
Let m := |supp(Q)| be the size of the support set of a mixing measure Q; i.e., the mixture complexity. If Q is finitely supported, m is finite and otherwise it is infinite. The current standard approach for finding the maximum likelihood estimator (MLE) of Q, when m is known a priori or fixed, is the well-known EM algorithm developed in the seminal article by Dempster et al. (1977).
In the absence of the knowledge of mixture complexity, it is instructive to start with an overparameterized mixture model and search over the whole continuous parameter space effectively to obtain a parsimonious mixture model. Overparameterization here refers to fitting a model with many components relative to the actual number in the nonparametric MLE (NPMLE) of Q; hence, there is a redundancy of components in the mixture model. Such a scheme would be robust to parameter starting values chosen for fitting the mixture model. To accomplish this, one requires a powerful mixture algorithm that pushes most of the mixture probabilities to zero, which is on the boundary of the parameter space. The focus of this article is to develop theory and create robust (to starting values) as well as powerful algorithms to address this problem.
The popular and widely employed EM-based algorithms are particularly slow for fitting such an overparameterized mixture model and it is very difficult, if not impossible, to remove the unnecessary components; since the algorithm can never reach such a boundary point.
Moreover, the EM algorithm is sensitive to parameter starting values and fails to converge in certain mixture problems. Figure 4 in Section 7.1 demonstrates this aspect of the algorithm.
Other examples where the EM algorithm converges to saddle points or fails to converge are noted by McLachlan and Krishnan (1997).

Statistical Framework
Let F := {f θ (x) : θ ∈ Ω ⊂ ℜ p } be a family of probability density functions with respect to a σ-finite dominating measure µ for a p-dimensional random vector x ∈ X ⊂ ℜ n and a p-dimensional location parameter vector θ ∈ Ω ⊂ ℜ p , a measurable space. Assume that the component density f θ (x) is bounded in θ for each x ∈ X . Let G be the space of all probability measures on Ω with the σ-field generated by its Borel subsets. For a given Q ∈ G, we assume that data vector x arises from the marginal density which is referred to as a mixture density. The mixture model (1) is also applicable to empirical Bayes estimation, where Q is an unknown prior distribution and the objective becomes estimation of the posterior distribution of θ without assuming a functional form for the prior distribution.
The goal is to estimate the mixing measure Q by finding the probability measure Q ∈ G that maximizes the nonparametric mixture loglikelihood log {L(Q)} = n i=1 log {g Q (x i )}. It is well known that finding the NPMLE of Q is computationally intensive (Lesperance and Kalbfleisch, 1992;Roeder, 1992;Lindsay, 1995;Bickel et al., 1998;Susko et al., 1999). Although Q is an arbitrary probability measure, under mild conditions Lindsay (1983a,b) showed that finding the MLE involved a standard problem of convex optimization, that of maximizing a concave function over a convex set. One consequence of this is that, as long as l(Q) is bounded, the MLE of Q is concentrated on a support of cardinality at most that of d-the number of distinct observed data vectors. This is a very useful, albeit surprising, result since a potentially difficult nonparametric estimation problem is reduced to that of a finite dimension; hence algorithms can be constructed to find the solution. Hence, we restrict the attention to discrete probability measures Q having p-dimensional support vectors θ 1 , . . . , θ m collected in a matrix Θ with a corresponding vector of masses or mixing probabilities denoted by π = (π 1 , . . . , π m ) T such that π is in the unit simplex We consider model fitting for both discrete and continuous data. Therefore, it is instructive to define f θ := {f θ (y 1 ), . . . , f θ (y d )} T to be the d-dimensional vector of distinct likelihood terms, where T denotes transpose, (y 1 , . . . , y d ) ∈ Y are the distinct observation vectors arising from the original data vectors (x 1 , . . . , x n ) ∈ X . Let n i be the number of times y i occurs in the sample of x vectors. For continuous data, n = d and for discrete data, often n ≫ d.
The observed data matrix of dimension (n × p) with row vectors y i (i = 1, . . . , n) is assumed to arise from the mixture density g Q (y i ) = m j=1 π j f θ j (y i ), and the discrete mixing measure can be represented as Q = j π j ̺(θ j ), where ̺(θ) is a discrete measure with mass one at θ ∈ Ω. If m is fixed, the model g Q (y i ) will be referred to as the m-component mixture model, and one can always find the NPMLE of Q using m equal to d (Lindsay, 1983a,b). However, the actual number of distinct support vectors with positive mixture probability, referred to as active supports, can be as small as one. The mixture loglikelihood of Q becomes (2) The goal is to find Q ∈ G such that l( Q) = sup Q ∈ G l(Q).
The biggest practical problem one faces in solving the loglikelihood equations in (2) is that the number of inequalities is equal to the number of elements of the parameter space Ω.
There are some important problems where this number is finite, although possibly very large, such as in target recognition, hyperspectral image analysis and positron emission tomography. For these problems discretization of the parameter space is directly relevant.
In other problems, Ω may be a continuous space; hence, one needs a machinery for approximating the parameter space to solve these equations (see Section 2.1).

Main Results
In this article, we develop a unified framework for finding the NPMLE of a multivariate mixing distribution Q and consequently for building a collection of semiparametric mixture models. The key ingredients for building these models are the "sieve parameter" controlling the dimensionality of the mixture problem (as shown in Figures 2 and 3) and the ability to fit overparameterized mixture models. This collection of models enable us to investigate the role of many overlapping densities, thereby creating an ideal situation for solving large-scale practical problems.
We create a powerful technique referred to as the "Penalized Dual method" and an efficient algorithm for fitting overparameterized mixture models. Consequently we have a method for estimating the mixture complexity. This algorithm is a step in the direction of developing a unified framework for building a collection of semiparametric mixture models.
The underlying principle for our method is that the mixture loglikelihood forms a concave functional on the convex set of all probability distributions which implies that there exists a dual optimization problem [Section 5.3, Lindsay (1995)]. Lesperance and Kalbfleisch (1992) and Susko et al. (1999) exploited this to create an elegant algorithm for finding the NPMLE of Q in univariate mixtures. This research is in the same spirit but extends these ideas by introducing a "penalty term". A fundamental feature of our approach is that it eliminates ad hoc procedures to estimate the penalty parameter. The dual problem has a statistical interpretation analogous to the least squares problem and the formulation is strikingly similar to the one that arises in empirical likelihood (Owen, 2001) framework (see Section 3.1).
Our main results are summarized as follows.
1. In Section 2, we cast the mixture problem in the dual optimization framework. We first develop a machinery for approximating the continuous parameter space Ω and next create an algorithm (based on the Penalized Dual method) for finding the maximum of l(Q). Consequently, we propose an algorithm for estimating the mixture complexity.
We establish the convergence of this algorithm to MLE in Section 5.
2. In Section 3, we develop theory for the Penalized Dual method in solving the dual optimization problem while presenting the statistical interpretation for our framework.
By exploiting the inherent advantage of the penalty formulation, we derive a technique for converting parameter estimators from the Penalized Dual problem into the mixture probability parameters. We show that the Penalized Dual estimators converge to the mixture probability estimators as the penalty is increased, with the correct limits.
3. Section 4 establishes the existence of parameter estimators and derives convergence results for the Penalized Dual algorithm, for fitting overparameterized mixture models.  Figure 4, Section 7.1.
7. Section 7 illustrates the power of the proposed methods using several applications.
We compare our method with the EM algorithm, due to lack of a unified and/or stable algorithms for fitting the overparameterized mixture problems and for building semiparametric mixture models. For the univariate mixture case, we compared with Rotated EM algorithm (an accelerated version of the EM only applicable to the univariate mixtures) proposed by Pilla and Lindsay (2001). Our empirical investigation demonstrate the faster rate of convergence of our algorithm, compared with the EM algorithm.
Section 8 presents the conclusions and the Appendix derives technical details.

Relevant Literature
Widely employed model-free methods for high-dimensional modeling include the K-means algorithm, hierarchical clustering and agglomerative and divisive algorithms (Hastie et al., 2001). However, none of these techniques take advantage of the inherent statistical structure of the data.
The existing model-based mixture algorithms include those for finding the NPMLE of Q (Lesperance and Kalbfleisch, 1992;Susko et al., 1999;Connolly et al., 2001). These algorithms are either not fast enough for high-dimensional modeling or not applicable for the following mixture problem: (1) the component densities are poorly separated and/or (2) many of the estimated mixture probabilities are on the boundary of the parameter space.
In analyzing Sloan Digital Sky Survey data by fitting the multivariate normal mixtures, Connolly et al. (2001) noted that many existing techniques are not computationally efficient and their mixture EM algorithm obtains an improvement of only three orders of magnitude. Therefore, developing a powerful method for fitting multivariate mixtures is desirable.
Although there have been some promising developments on accelerating the EM algorithm (see McLachlan and Krishnan (1997) and the references therein), none of these methods address the overparameterized mixture problem described earlier. To overcome the above difficulties, Lindsay (1996, 2001) proposed alternative augmentation schemes based on the principles of the EM that provide a significantly improved convergence rate of the EM algorithm for a class of finite mixture models. At this time, it is not clear how to extend these methods to multivariate mixture models; however, they do provide an important class for comparison with our algorithm in univariate mixture problems (comparisons are made in Section 7.1). Lindsay (1995, Section 6.3) discusses several algorithmic methods based on directional derivatives such as the vertex direction method and vertex exchange method to find the NPMLE of Q. These methods also require searching over a discrete parameter space and have certain computational disadvantages (Lesperance and Kalbfleisch, 1992).

Mixture Maximum Likelihood Problems
In this section we first formulate the mixture problem as a convex optimization problem and next create a framework for approximating the continuous parameter space. Lastly, we develop an algorithm for finding the NPMLE of Q. This algorithm forms the basis for building a collection of semiparametric mixture models developed in Section 6.

Maximizing l(Q) via Approximating G
If the number of components in Q is fixed, but the location parameter vectors are unknown, then l(Q) can have several local maxima (Lesperance and Kalbfleisch, 1992;Lindsay, 1995;Pilla and Lindsay, 2001). Both the EM and the K-means algorithms can get trapped at a local maximum while requiring a priori knowledge of the mixture complexity m. To overcome this problem, researchers often randomly perturb the parameter starting values and recompute the local maxima (Hall and Zhou, 2003;Hunter, 2004 To combat the difficulties with the parameter starting value problem and an a priori knowledge of m, we develop a technique in which we approximate G by the set of all probability measures on a discrete parameter space of Ω. Approximating G: Approximate G by G m , where G m is a set of discrete distributions generated by a finite subset of Ω. We set this finite subset to be Θ m = (θ 1 , . . . , θ m ). As m → ∞ and Θ m becomes dense in Ω, the set G m → G. In practice, a sufficiently large m is chosen such that G m approximates G well. Therefore, the cardinality of Θ m , namely m, determines how close the MLE is to the global MLE over all measures on Ω. In approximating G, it is important to select a suitable Θ m while keeping computations manageable. This will be addressed in Section 6.3.
In what follows, we distinguish between the three mixture problems.
1. The fixed support mixture problem is equivalent to maximizing over the parameter space Π while treating the support set Θ m ⊂ Ω as fixed. This is the primal or mixture problem for which we define a "dual" in Section 3. Note that dim(π) = (m − 1) and grows with the cardinality of Θ m , which is a major obstacle when dim(Θ m ) is large. However, the dimension of our dual optimization problem is fixed at d, the number of distinct observed data vectors.
2. We fix the number of components in the mixing distribution Q to be m but treat the θ parameter vectors as unknown for each component. Therefore, the continuous support mixture model problem becomes simultaneously estimating π and θ parameter vectors by maximizing l(Q) over Π × Θ m for a fixed m.
3. In the absence of knowledge of mixture complexity m, maximizing the mixture loglikelihood in (2) yields an NPMLE that is a discrete distribution on the parameter space with a random number of component densities (Lindsay, 1995;Pilla and Lindsay, 2001). This will be referred to as the nonparametric mixture model. The goal in turn becomes finding the probability measure Q ∈ G that maximizes (2).
For discrete mixture problems, such as binomial or Poisson, often d ≪ n; therefore, the dual methods are able to reduce the dimension of the mixture problem. The effect of this dimensionality on the performance of the algorithms will be demonstrated in Section 7.1.

Characterization of the NPMLE of
Under compactness of Γ, we can define the convex hull of Γ as Conv(Γ) = {g Q : and a corresponding maximizing measure Q can be characterized in terms of the gradient function as shown next.
Definition 1 (Finite identifiability): For a given family F, suppose that Q 1 , Q 2 ∈ G have finite support. Suppose that Q j ∈ G yields the mixture density g Q j (y) for j = 1, 2. If g Q 1 (y) = g Q 2 (y) for all y ∈ Y implies Q 1 = Q 2 , then the corresponding collection of mixture densities is said to have the finite identifiability property.
An important aspect of our technique is based on the following fundamental property.
For the NPMLE Q, the ith fitted model g Q (y i ) is guaranteed to be unique (regardless of identifiability of the mixture density), and that one can determine these fitted values by solving for the residual w i , on a log-scale, defined as In ordinary parametric likelihood problems the solution is characterized by the likelihood equations. We extend these ideas to our problem to show that the fitted values and the corresponding mixing distribution Q ∈ G can be further characterized in terms of a set of gradient equations. That is, Q is an NPMLE if and only if where the gradient function, the directional derivative of the mixture loglikelihood in the direction of a component density, is defined as If a candidate maximizing measure Q violates the gradient inequality in (5) at some θ ∈ Θ m , then one is not at the maximum. In particular, one can increase the loglikelihood by placing some positive probability at θ ∈ Θ m .

Finding the Maximum of l(Q) via Θ m
For a fixed m, mixture estimation is challenging due to the fact that l(Q) is not concave and hence there are several local maxima (Lesperance and Kalbfleisch, 1992;Lindsay, 1995;McLachlan and Peel, 2001;Pilla and Lindsay, 2001). We create an algorithm that is robust to the choice of parameter starting values and reaches closer to the global maximum of l(Q).
We find the NPMLE of Q adaptively as follows.

Algorithm 1 [Finding the Maximum of l(Q)]
1. Consider Θ m ⊂ Ω to be the support set of Q. Solve the fixed support mixture problem by maximizing (3) over the parameter space Π on the support set Θ m , while treating θ ∈ Θ m as fixed. It is worth noting that the larger the cardinality of Θ m , the higher the value of the loglikelihood at convergence.
2. Apply the MLE π (with the corresponding fixed support set Θ m ⊂ Ω) obtained in Step 1, as parameter starting values for the continuous support mixture problem and maximize (2) over the product parameter space Π × Θ m . For Step 1, one requires a stable and powerful mixture algorithm and is derived in the next sections. In particular, Algorithm 2 presented in Section 4.2 can be employed in Step 1. The Step 2 may include estimation of other parameters in the model such as Σ ∈ S, in the multivariate normal mixtures context.
For the continuous support mixture model, it will be shown in Section 7 that Algorithm 1 reaches closer to the global maximum, if not to the global maximum. Our empirical evidence suggests that Algorithm 1 is superior to EM-type algorithms that start with random (or arbitrary) parameter values.

The Dual Optimization Problem: Properties of Estimators
We now present the problem that is dual to the primal problem (3) considered by Lindsay (1983a) and develop theory for effectively solving it. The dual problem is to maximize subject to the constraints w = (w 1 , . . . , Let w ∈ ℜ d + be the solution to the above dual (or concave) optimization problem. The solution satisfies the relationship (4), so that solving the dual problem for w is equivalent to finding the log-scale residuals. Hence, indirectly, via (4), we obtain the model fitted values g Q . A challenging step is that one must solve for the parameter estimates for the model from these fitted values. We create a method that exploits the particular choice of our penalty term. Note that the constraints are linear in the parameter vector w ∈ ℜ d + , and that the number of free parameters equals d, while the number of constraints equals the cardinality of Θ m . The dual optimization is with respect to w whose dimension equals d. This is especially advantageous with large data sets containing, say, thousands of observations (see Section 7.1). In Appendix A.1, we establish the relationship between the primal and dual problems at the solution.

Statistical Interpretation of the Dual Problem
The formulation in (7) and (8) is strikingly similar to the one that arises in empirical likelihood framework (Owen, 2001) in which the function l(w) is maximized over a similar set of linear constraints. The empirical likelihood problem also has a dual problem, although it does not appear to be computationally useful.
There is a natural interpretation of the dual problem that is analogous to the linear model framework. In an application of the least squares problem, one finds the fitted values y directly by projecting the data y ∈ Y onto the model space X (i.e., P X y = y) or solves for the residual e by projecting y onto the orthogonal complement of the model space (i.e., P X ⊥ y = e, where ⊥ denotes the orthogonal projection). In turn, we solve for the fitted values using y = (y − e). The primal and dual problems have the same relationship as the projection and complementary projection of linear models. The parallel with the linear model framework holds if we let the data y i equal log (n i /n), the fitted model y i equal log {g Q (y i )} and the log-scale residual e i equal log ( w i ). This approach again falls very much into the spirit of the empirical likelihood, where (n i /n) is the NPMLE of the probability of observing y i ∈ Y.

The Penalized Dual Method: Theory
The goal in this section is to turn the constrained dual optimization problem defined in (7) and (8) into an unconstrained one using a "penalty function". This is referred to as the Penalized Dual method. Our method is in the spirit of the log-barrier method (Renegar, 2001) for convex programming; however it differs in two important respects as will be shown.

The Penalized Dual method maximizes
over w ∈ ℜ d + , where γ is a tuning parameter and P(w, γ) is a penalty function that ensures that the Penalized Dual solution does not violate the constraints; the dual solution always stays in the interior of the constraint set. One choice for the penalty function is where the penalty parameter γ is some large power and the constraint function is defined as That is, the dual problem constraints have the form p θ j (w) ≤ 1. We first show that by increasing γ, P(w, γ) will eventually create an infinite penalty on any w ∈ ℜ d + that violates the constraints and advances the solution towards the dual problem solution.
The proofs for our technical results are derived in the Appendix.
Proposition 1 For a given θ ∈ Ω and w ∈ ℜ d + , the term in the summand of the penalty function P(w, γ) satisfies: If p θ (w) < 1, the penalty function is decreasing in γ for all γ ∈ ℜ + .
Two main elegant features of our penalty function are the following: (1) We can directly construct an estimator for π parameters from the penalized dual solution.
(2) It is simple to calculate the gradient function to assess the algorithmic convergence using the relation (20), defined in Section 3.4.
It is common in the optimization literature to employ a "barrier function" to build the penalty. For example, the log-barrier function defined as approaches −∞ as w ∈ ℜ d + approaches the boundary of the feasible set from the interior (Roos et al., 1997;Renegar, 2001). The effect of the penalty can be diminished by making γ close to 0. Our focus here is on a soft penalty of the form (10) which is well behaved outside the feasible set; however, as will be shown, it does force the solution into the interior.
The penalized problem is unconstrained; therefore, we can find the "Penalized Dual optimal estimator" denoted w γ = ( w 1,γ , . . . , w d,γ ) T , given by (A.6) in Appendix A.2, by solving For γ = 1, there exists an explicit solution to the above equation as This is an initial interior point solution for the algorithm. On the other hand, the conventional log-barrier methods do not automatically produce a starting value for the "barrier parameter".

Existence of Parameter Estimators
We consider the following method to solve for the π parameters from the dual problem solution by exploiting the penalized structure.
Recovering the Primal Estimators: Using (4), the model fitted values g Q are found from the penalized dual solution w. However, such a solution does not immediately provide an estimator for π ∈ Π and the technique for obtaining it is derived next.
The penalized dual residuals are used to obtain a natural estimator for the mixture or primal problem, denoted by π ⋆ γ = π ⋆ 1,γ , . . . , π ⋆ m,γ T . The statistic, which is referred to as the Penalized Dual estimator is where In Appendix A.2, it is shown that the estimator w i,γ , derived in (A.6), can be approximated in terms of {p θ j ( w γ )} (γ−1) . However, these latter quantities with the power (γ − 1) do not sum to one, and hence are turned into a candidate estimator via normalization: This candidate estimator is used to obtain π ⋆ γ with its elements having the power γ using the following theorem.
Theorem 2 (a) For a given γ ∈ ℜ + , the Penalized Dual estimator is one EM-step from the candidate estimator π † γ ; consequently, π ⋆ γ yields a higher likelihood value. (b) The estimators are in the unit simplex The Penalized Dual solution w γ satisfies p θ j ( w γ ) ≤ 1 (j = 1, . . . , m) and hence the estimator π ⋆ γ remains in the feasible region defined by (8).
All the proofs are relegated to Appendix A.4.
The estimator π ⋆ provides a direct way to obtain the primal estimator π from our penalized dual solution, avoiding the problems of selection and inversion.

Properties of the Penalized Dual Estimators
In this section, we derive several statistical properties of the estimators. First, we establish that π ⋆ γ converges to the MLE π as the penalty parameter γ increases (Theorem 3 below). Along the way, we establish several important properties of the primal-gradient function that are necessary for solving the primal-dual problem.
Although π ⋆ γ represents an EM improvement over π † γ , the candidate estimator, it is easier to establish optimization results for the latter. The following theorem shows that for a sufficiently large penalty, the Penalized Dual estimator will be close to the primal estimator.
Let Q † γ be the mixing distribution at the π † γ solution.
Our goal is to obtain the mixture estimation problem from the penalized dual one using π ⋆ γ . Therefore, it is important to determine directly from the dual problem how accurate is the estimator π ⋆ γ . We derive the gradient function corresponding to the mixing distribution Q † γ to accomplish this. It is easier to calculate this for π † knowing that π ⋆ γ can only be better. From the primal-gradient function in (6), the gradient function for the estimator Theorem 4 The primal-gradient function at the candidate estimator π † can be written as where Theorem 4 expresses the gradient function in terms of the dual solution and leads to a simpler device for checking the accuracy of the estimators.
Corollary 5 At the candidate estimator π † , the primal-gradient function satisfies where the term on the right-hand-side does not depend on j.
We have established that one can refine the NPMLE of Q to the required accuracy by increasing m and γ appropriately.

The Structure of the Penalized Dual Algorithm
In this section, we first investigate the structure of the penalized dual problem viewed as a function of w and γ and next present a strategy for their joint estimation. Next, we present the Penalized Dual algorithm to effectively search over the discretized (but large) parameter space Θ m ⊂ Ω. Lastly, convergence properties of the algorithms are derived.
We let z = log (w) to eliminate the constraint w ∈ ℜ d + .
Theorem 6 (a) The function is strictly concave in (z, γ), where p θ j (z) = i exp(z i ) f θ j (y i ). For any z ∈ ℜ in the feasible region defined by (8), the function K(z, γ) is strictly increasing as a function of γ. (b) The function K(z, γ) is bounded above and achieves its maximum at z = z and γ = ∞.

Automatic Selection of γ
A fundamental aspect of our algorithm is that we can maximize K(z, γ) simultaneously with respect to z ∈ ℜ and γ ∈ ℜ + ; different from the approach employed in the conventional log-barrier methods in which this was not possible. Therefore, we can select the penalty parameter γ automatically. From Theorem 6, the global maximum over z and γ is attained when z = z and γ → ∞.
Remark 1: It may seem paradoxical to treat γ as an unknown parameter even though it has an optimum value of ∞. When γ is large, K(z, γ) has very severe curvature at the constraint boundary. This limits the range of effectiveness of quadratic approximation methods. Therefore, one should start with a small value for γ and increase it as the algorithm progresses through the parameter space. This could possibly be achieved in some other systematic fashion; however, our empirical investigations suggest that systematic methods were not as efficient as our approach. A possible explanation could be that our strategy takes the curvature of the function K(z, γ) into account, in providing the relevant information for determining the increments for γ.

Searching Effectively Over the Discretized Parameter Space
An algorithm for efficiently searching over the large discretized parameter space Θ m ⊂ Ω (required for Step 1 of Algorithm 1 described in Section 2.1) is derived next. In effect, the following algorithm is used for fitting the fixed support mixture model.

Algorithm 2 (The Penalized Dual Algorithm)
1. Consider γ = 1 and its corresponding explicit solution z (1) given in (13) as the starting solution for the algorithm.
2. Maximize the concave function K(z, γ) simultaneously with respect to z ∈ ℜ and γ ∈ ℜ + using a modified Newton-Raphson algorithm [constrained the step size to ensure monotonicity in K(z, γ)] until the following convergence criterion is satisfied; namely, the L 2 -norm of the change in the value of K(z, γ) is less than 10 −6 .
3. Fix γ at γ (k) obtained in Step 2 and find z γ (k) = arg max z∈ℜ K z, γ (k) using the modified Newton-Raphson algorithm. The algorithm is considered to have converged to the maximum at step t, when the inequality based on the primal-gradient function is satisfied since it guarantees convergence to a similar accuracy in the loglikelihood.
As described in Section 2, the supremum of the gradient function sup θ ∈ Θm D Q (θ) provides an assessment of the progression to the maximum and hence the criterion in (22) has a solid theoretical justification.

Remark 2: The
Step 3 of the algorithm is necessary since after Step 2, the Primal Dual estimator π ⋆ γ obtained via (14) are often not sufficiently close to the primal estimator π. This is because the algorithm does not necessarily satisfy the condition {∂K(z, γ)/∂z} = 0 with sufficient accuracy. In our applications, however, the primal-gradient inequality (22) was always achieved at the tolerance of 0.005; in fact, often reached significantly greater accuracy in the Penalized Dual estimators.
The Penalized Dual Algorithm with Inactive Constraints: In Algorithm 2, if an estimated mixture probability π j (j = 1, . . . , m) is zero, then the corresponding constraint in the dual problem is inactive. We can dynamically update the active constraints by removing the inactive ones while adding new ones, whenever the support set violated the gradient inequality. From the Penalized Dual estimator in (14), it follows that if π j → 0, then In the former, one can essentially remove the corresponding density f θ j (y i ). It will be shown in Section 7 that the above algorithm, denoted by PD IC , produced a further reduction in computational time.
As a consequence of the concavity of K(z, γ) established in Theorem 6, the Hessian H for K(z, γ) (derived in equation (A.7) in Appendix A.3) is always non-singular and the sequence obtained from the Penalized Dual algorithm (i.e., Algorithm 2) are well defined. Even for large-scale problems such as the yeast microarray data considered in Section 7.3 in which H is of dimension 697, our modified Newton-Raphson algorithm was stable and efficient.
Theoretically, owing to Theorem 6, the Step 2 of the Penalized Dual algorithm produces a sequence z (k) , γ (k) k≥1 such that the sequence of functions K(z (k) , γ (k) ) k≥1 → K( z, ∞) as k → ∞. This effectively implies that the sequence z (k) k≥1 → z and the sequence γ (k) k≥1 → ∞ as k → ∞. However, in practice, convergence of γ is slow; therefore, we terminate the modified Newton-Raphson algorithm in Step 2 when γ (k) is sufficiently large and maximize K z, γ (k) over z for a fixed γ (k) .
In our experience, a direct maximization of the mixture loglikelihood l(π) over Π using a modified Newton-Raphson algorithm was unstable and failed to converge to the maximum π.

Convergence Properties of the Algorithms
In this section, we establish the convergence properties, including the rate of convergence, of the algorithms.
First, we consider the algorithm for fitting the continuous support mixture model; i.e., an algorithm employed in Step 2 of Algorithm 1. We prove that the sequence of estimates {β (s) } s≥1 obtained from the Step 2 of Algorithm 1 converges to an MLE of β ∈ Ω, namely β, for a given data y ∈ Y as s increases. For instance, in the multivariate normal mixture framework, β becomes (Q, Σ). Assume that the sequence of estimates {β (s) } s≥1 monotonically increases the loglikelihood l(β). An algorithm is said to converge if β ⋆ = lim s β (s) exists, for a parameter vector β ∈ Ω. Wu (1983) established that monotonicity of l(·) does not imply the convergence of the sequence to a stationary point; however, if the sequence {l(β (s) )} s≥1 is bounded above, then it does converge monotonically to a stationary point of l(β). The convergence of β (s) to β implies the convergence of l(β (s) ) to l( β) according to the Theorem 5, under the regularity conditions, derived by Wu (1983).
Owing to Theorem 6, the Algorithm 2 (or Step 1 of Algorithm 1) produces a sequence of estimates {π (t) } t≥1 that is guaranteed to converge to the unique MLE π. Combined this result with Theorem 5 in Wu (1983) establishes the convergence of Algorithm 1 to β.

Convergence Criteria
For the applications and simulation experiment, we used the convergence criterion based on the gradient function for the Penalized Dual (PD) and discrete EM (i.e., for fitting the fixed support mixture model) algorithms. That is, the algorithm has converged to the MLE π if the criterion in (22) is satisfied. For the rest of the article, we denote the discrete EM by D-EM algorithm.
The D-EM algorithm is a sublinearly convergent algorithm (Pilla and Lindsay, 2001); therefore, a conventional convergence criterion based on the loglikelihood change or changes in parameters, such as for a given tolerance τ can be very misleading in the sense that the actual distance to the final loglikelihood can be orders of magnitude different from τ . That is, this criterion may be met even though the parameter values are far from the correct solution (Titterington et al., 1985;Pilla and Lindsay, 2001). However, such rules are widely employed and therefore we conducted an experiment to assess the two criteria on two data sets.
The most important assessment of the convergence of an ML algorithm is the value of the loglikelihood, as it provides information about the accuracy of parameter estimators on a confidence interval scale. Therefore, loglikelihood-based criterion is a useful one to employ in assessing the convergence of an algorithm in finding the MLE of the parameters (Lindsay, 1995;Pilla and Lindsay, 2001).

Simulation Experimental Design:
We consider the simulated data by generating a sample of size n = 270 from N p (Q, I) with p = 3, where N p (Q, I) represents a measure of a p-dimensional normal random variable with mean Q and an identity variance-covariance matrix. The true mixing measure for Q ∈ G, is chosen by selecting the coordinates of θ j (j = 1, . . . , m) from the set {−5, 0, 5} in all possible combinations, with equal mass at each support vector. This resulted in a total of m = 3 3 mixture components.
Fisher Iris Data: We fit a mixture of multivariate normal distributions to Fisher iris data (Fisher, 1936). The data consists of n = 150 observations collected on flowers of three iris species (Setosa, Verginica and Versicolor). Each observation is a vector of p = 4 variables sepal length (y 1 ), sepal width (y 2 ), petal length (y 3 ) and petal width (y 4 ). For each of the data sets, we selected the observed data matrix y for Θ m and also set Σ = S, the sample variance-covariance matrix. In order to assess the accuracy of the algorithms at a given step t, we found the final loglikelihood value l( π) to a high degree of accuracy using the PD algorithm for a sufficiently large t. Next, we fit mixtures of multivariate normal distributions to the simulated and iris data sets via the D-EM algorithm using the convergence criterion (23) with τ = 0.0001. The Λ (t) values, presented in Table 1, demonstrate that the convergence criterion (23) would result in substantially less than four decimal accuracy for the Fisher iris data. On the other hand, the criterion based on Ψ(Q (t) ) in (22) guarantees the final accuracy.

Empirical Assessment of Convergence Rate
We empirically assess the rate of convergence of the PD algorithm relative to the D-EM algorithm by defining Λ (t) in (24) as the residual of the loglikelihood at the tth step.
To be precise, for some π (0) ∈ Π, let π (t) t≥1 be a sequence in Π generated by an algorithm (such as the PD and D-EM algorithms). The algorithm can be expressed as π (t) ∈ M π (t−1) for t ≥ 1, where the map M : Π → 2 Π is a point-to-set mapping. If π (t) converges to π and M(·) is continuous, then π must satisfy π ∈ M( π).
Definition 2 (Asymptotic Convergence Rate): Assume that π = M( π) and that the sequence {π (t) } t≥1 is generated by the map M such that lim t→∞ π (t) = π. Under the regularity conditions given by Wu (1983), this implies that lim t→∞ l π (t) = l( π). The asymptotic convergence rate of the loglikelihood sequence l π (t) t≥1 at l( π) generated by an algorithm is defined as From the following lemma (Pilla and Lindsay, 2001), the smaller the r for any given loglikelihood sequence, the faster it is progressing towards the MLE.
Lemma 7 If the sequence {l(π (t) )} t≥1 is converging linearly, then as t → ∞, the slope of the curve obtained by plotting log {Λ (t) } against t converges to log (r), where r is the asymptotic rate of convergence of the loglikelihood sequence generated by an algorithm.
In order to assess the rate of convergence of the PD, relative to the D-EM, algorithm, we consider the Fisher iris data considered earlier for fitting a collection of semiparametric mixture of multivariate normal distributions N p (Q, δ Σ), where p = 4, Σ ∈ S and δ ∈ ℜ + (details in Section 6.2). As before, we selected the observed data y for Θ m and set Σ = S in fitting the PD and D-EM algorithms. Figure 1 demonstrates the behavior of the algorithms for δ ∈ {5, 2, 1, 0.5}. We used logarithmic scaling of the vertical axis since a linearly convergent algorithm will become linear on this scale as t → ∞. Note that for the fixed support mixture model, the D-EM algorithm is converging sublinearly whereas the PD is converging linearly to the MLE π; a significant improvement in convergence rate. In fact, Pilla and Lindsay (2001) observed a similar behavior of sublinear convergence of the D-EM algorithm for a class of univariate finite mixture problems.

Semiparametric Mixtures of Multivariate Normal Distributions
The methodology developed in this article is applicable to a wide range of problems, including multivariate t mixtures. However, the particular interest here is in difficult problems with multivariate normal mixtures due to its ubiquitous applications. In semiparametric mixture setting, the mixing distribution Q is modeled nonparametrically in the presence of an unknown Σ ∈ S, the variance-covariance matrix common to all m components.

Structural Properties
Let g Q (y i ; Σ) = m j=1 π j f µ j (y i ; Σ) for y i ∈ Y be a finite mixture of p-dimensional normal distributions, where µ j ∈ ℜ p ⊂ Ω is the mean vector of the jth component density f µ j (y i ; Σ) and Σ ∈ S is common to all m components. Note that in the continuous case d = n. The corresponding loglikelihood is expressed as In the univariate case, Charnigo and Pilla (2005) establish that for a general family of mixture models with a structural parameter β (e.g., σ 2 in the normal case), the likelihood framework breaks down when joint estimation of m, Q and β is attempted: at best the joint estimator of m, Q and β is degenerate, and at worst it does not even exist. The ML fails in this setting since taking finite samples from continuous probability distributions yields discrete data sets. When models that closely mimic discrete probability distributions are available, as they are when there are no restrictions on Q and β, the likelihood will favor such models. The NPMLE results of Lindsay (1995, Section 2.6) cannot be applied if Σ is unknown; however, the following result holds.
We define the gradient function for the multivariate normal mixture distributions as Next, we define an NPMLE of Q ∈ G for a fixed Σ ∈ S as Theorem 8 (Unique NPMLE of Q) Assume Σ > 0 is fixed.
(1) Suppose Q Σ satisfies then Q Σ is an NPMLE of Q ∈ G.
For the multivariate normal mixture model with a common Σ, we restrict attention to finite discrete latent distributions Q, then the pair (Q, Σ) is identifiable (Lindsay, 1995).
For a general family of univariate mixtures, Charnigo and Pilla (2005) establish that joint estimation of m, Q and β is a well-defined problem if Q is finitely supported. If Q is not finitely supported, then g Q (y; Σ) need not determine Q and β uniquely. Hence, an ML approach to the joint estimation of m, Q and Σ fails. However, since we fix m and consider Q to be finitely supported, joint estimation of Q and Σ is feasible. Therefore, we can apply Algorithm 1 described in Section 2.1 to jointly estimate Q and Σ.
The following theorem establishes that joint identifiability of (Q, Σ) fails if Q is not finitely supported. The proof follows from the univariate nesting structure result, under mild regularity conditions, given by Charnigo and Pilla (2005). where N p (Q, Σ) represents a measure of a p-dimensional normal random variable with mean Q and a variance-covariance matrix Σ ∈ S.

Models
We investigate building the sieve of models N p (F, δ Σ), where Σ ∈ S and δ ∈ ℜ + is a sieve parameter (similar to the smoothing parameter employed in density estimation). The sieve parameter controls the dimensionality of a mixture model as will be demonstrated later.
We derive theory for building a collection of semiparametric mixture models, including the multivariate case.
In order to create a general family of mixture models, we consider the class {N p (Q, δ Σ) : Q ∈ G, Σ ∈ S, δ ∈ ℜ + }. For δ 1 > δ 0 , Theorem 9 implies that hence, the collection of models becomes richer as δ → 0 (see Figures 2 and 3). Moreover, every p-dimensional distribution F can be obtained as the weak limit of N p (F, δ Σ) as δ → 0.
Therefore, we can approximate any distribution by choosing δ small. As a consequence, the principle of maximum likelihood cannot be applied to select δ in the model N p (Q, δ Σ) since the likelihood becomes unbounded as δ → 0. Charnigo and Pilla (2005) develop theory for the univariate mixtures and demonstrate the effect of small δ for a general family of univariate mixtures which extends to the multivariate case considered here.
We create a strategy for building a collection of models {(Q δ , δ Σ) : Q ∈ G, δ ∈ ℜ + } using the Penalized Dual algorithm. As δ → 0, the NPMLE Q δ converges in distribution to n −1 i ϑ(y i ), where ϑ(y) is a discrete measure concentrated at y.
To demonstrate the effect of δ on the mixture complexity, we create a collection of models for both the univariate and multivariate data. In the univariate case, we simply have a σ parameter. The univariate application considers the galaxy data set [ Table 1 of Roeder (1990)] of 82 observations of relative velocities for galaxies from six well separated conic sections of the Corona Borealis region. Scientific interest lies in identifying substructures in clusters of galaxies. Multimodality is evidence of voids and superclusters in the far universe. Roeder (1990) obtained σ = 0.95 using least squares cross validation. We set µ ∈ Θ m = {9, . . . , 35} with a grid size of 0.02 for building a collection of semiparametric mixture models using Algorithm 2. The plot of log (σ) against the support set µ corresponding to the estimate π obtained using the PD algorithm (namely, Algorithm 2) is shown in Figure   2. That is, at each fixed log (σ), the plot displays µ parameter values that have positive mixture probability. The figure demonstrates the effect of σ on the mixture complexity m.
Next, we consider the Fisher's iris data described earlier. Once again, we selected observed data y for Θ m and set Σ = S. We let δ = {0.1, . . . , 5} for building a collection of semiparametric multivariate mixture models using Algorithm 2. Figure 3 shows the effect of δ on the mixture model complexity when only the two variables, namely the petal length and petal width are considered. The galaxy and Fisher iris data sets demonstrate that the number of components is a consequence of the choice of σ (or δ as the case may be) rather than a pre-selected parameter. 6.3 Selection of the Support Set of Q As described in Section 2.1, in approximating G, the biggest challenge is in selecting a suitable Θ m while keeping computations manageable. This is addressed in this section.
In the absence of a prior knowledge of the mixture complexity, correct specification of Θ m ⊂ Ω, the support set of Q, is very important for the Step 1 of Algorithm 1 (or equivalently for Algorithm 2). As expected, the final loglikelihood depends on this choice. In this section, we illustrate through the simulated data described earlier how the observed data matrix y provides the best choice for approximating the continuous parameter space Ω. In effect, we select {θ 1 = y 1 , θ 2 = y 2 , . . . , θ m = y d }. Note that choosing y for the discrete parameter space Θ m clearly covers the region of likely support vectors for the normal means and has the advantage of adapting naturally in richness to the sample size of the problem.
To assess the effectiveness of using y for Θ m (which is approximating the continuous parameter space Ω) we consider the simulated data described in Section 5.1. The true mixing measure for Q chosen for the simulation experiment is denoted by "True Support" in Table   2. The "Equi-Distant" set for Θ m was constructed on a lattice by choosing the elements in θ j = (θ j1 , θ j2 , θ j3 ) for j = 1, . . . , 8 3 from the set {−7, −5, . . . , 5, 7} resulting in a total of m = 8 3 support vectors; this set also included all the true support vectors. Table 2  results obtained using Algorithm 2 (i.e., fixed support mixture model of estimating π for a given Θ m ) and Step 2 of Algorithm 1 (i.e., continuous support mixture model of estimating Q and Σ for a fixed m).
In general, y should be an effective choice for Θ m given that equi-distant is still a subjective one in the absence of any knowledge about the length of the distance. From the theory presented in Section 2.1, as m → ∞, Θ m → Ω. However, in practice, choosing m = n is effectively creating a dense set for Θ m and in fact approximating Ω very well. Table 2: Effect of Θ m ⊂ Ω, the support set for Q, on the estimated loglikelihood. We set Σ = S in finding l( π) using the Penalized Dual (PD) algorithm. The solution π obtained from the PD algorithm with the corresponding fixed support set and S are employed as parameter starting values for finding l Q δ , δ Σ for a fixed δ = 0.2 using the continuous -2181.9 -1936.9 Equi-Distant -2182. 8 -1901.7 Observed Data -2178.6 -1876.0

Applications and Simulation Experiment
The applications in this section are used to investigate the roles of many overlapping components which create an ideal situation for solving the large-scale practical problems.
We assess the performance of the algorithms in finding the NPMLE of Q and for fitting the collection of semiparametric mixture models with applications to several data sets. The data sets, the parameter estimates and the Matlab software for fitting mixtures are available from the first author. For the Step 1 of Algorithm 1, we set Σ = S; however, we estimate it in Step 2.

Mortality Data
Our first application illustrates the tremendous advantage of our method in reducing the dimension of discrete mixture problems. In these problems the magnitude of d, the number of distinct observed data points, could be much smaller than m, the cardinality of Θ m .
We consider the data on death rates which gives the number of death notices for women aged 80 and over, from the Times newspaper for each day in the three-year period 1910 to 1912 (Titterington et al., 1985). For the later data sets, we chose the observed data matrix y as the support set Θ m . However, for this application, we selected the support set to be We fit a mixture of Poisson distributions to the mortality data using the PD and D-EM algorithms. Table 3  At η = 0.01, the D-EM has retained 199 support points with non zero probability at convergence (obtaining a smaller loglikelihood value) whereas the PD has retained just 26 support points and reached the MLE in a reasonable number of steps; a significant reduction Table 3: Building a collection of finite mixture models Q η : η ∈ ℜ + using Algorithm 1 for the mortality data. First step involved setting Θ m = {0, 0 + η, . . . , 9 − η, 9} with η ∈ {1, 0.5, 0.1, 0.01} and finding l( π) using the PD-based and discrete EM (D-EM) algorithms.
The estimate π obtained from the PD algorithm with the corresponding fixed support set is used as parameter starting values for finding l Q η using the continuous EM (C-EM) algorithm.

Simulation Experiment
We consider the simulated data described in Section 5.1. The data were generated from the multivariate normal mixture densities N p (Q, I) with p = 3 and n = 270 (see Table 2) by selecting the true mixing measure for Q ∈ G as the coordinates of θ j (j = 1, . . . , m) from the set {−5, 0, 5} in all possible combinations, with equal mass at each support vector.
Following Section 6.2, we apply the Penalized Dual algorithm in the context of building a collection of semiparametric mixture models for selected values of the sieve parameter δ.
This will yield estimators with both many and few active support vectors; thereby providing a mechanism to demonstrate the superiority of our method over the D-EM algorithm across a range of applications. Both the PD and PD IC algorithms provide uniformly better performance, producing 6 to 40-fold improvement in CPU factor over the D-EM algorithm.
As illustrated in Section 6.2, in fitting N p (Q, δ Σ), there is a trade-off between decrease in the sieve parameter δ and the increase in mixture complexity m; by increasing δ, we obtain a reduction in the mixture complexity m.
As discussed in Section 5.1, an important attribute of the convergence of an algorithm is the value of loglikelihood, as it indicates accuracy on a confidence interval scale. Therefore,  using Algorithm 1 for the simulated, Fisher iris and Yeast microarray data sets. First step involved choosing the observed data matrix y for Θ m and setting Σ = S in finding l( π) using the PD and discrete EM (D-EM) algorithms. The estimate π obtained from the PD algorithm with the corresponding fixed support set and δ S are employed as parameter starting values for finding l Q δ , δ Σ using the continuous EM (C-EM) algorithm.  Table 4, we present the loglikelihood values obtained using various algorithms. The CPU factor for the PD algorithm over the EM algorithm ranged from ten to over forty-fold for δ ∈ {5, 2, 1, 0.5, 0.2}. As predicted by the theory, for the PD and D-EM algorithms, the final Ψ(Q (t) ) given by (22) does provide a guarantee on the level of algorithmic convergence Λ (t) . Indeed, in some cases the bound Λ (t) ≤ Ψ(Q (t) ) was very conservative. Moreover, the convergence criteria for the PD algorithm described in Section 4.1 achieved the desired accuracy in Λ (t) ; however, typically the PD algorithms terminated at a considerably higher accuracy than the D-EM algorithm. In order to measure this effect, we continued the D-EM algorithm to the same level of accuracy as that of the PD for δ = 1. In this case, for the EM algorithm, N (t) = 9, 987 at convergence, resulting in a CPU factor of 60 instead of 24.

Fisher Iris and Yeast Microarray Data Sets
We fit a mixture of multivariate normal distributions to the Fisher iris data described earlier by finding the NPMLE of Q and by building a collection of semiparametric mixture models for selected values of δ; results are presented in Table 4. The performance of the PD and D-EM algorithms was similar to that of the simulated data.
Instead of choosing the PD solution as the parameter starting values for the C-EM, we consider random values to demonstrate their effect on a given algorithm. It is important to recognize that the C-EM algorithm requires an a priori knowledge of m. We considered the Fisher iris data with δ = 1 and randomly selected m = 15 data vectors from n = 150 as parameter starting values for the algorithm. In the ten runs of the C-EM algorithm with random starting values, l( Q, Σ) ranged from -151.99 to -179.13; all of which are sub-optimal modes due to the "poor choice" of starting values. Similar behavior of the C-EM algorithm, in reaching a sub-optimal solution, was observed by Pilla and Lindsay (2001) for the galaxy data. Without an a priori knowledge of m, choosing m can be quite a challenge for using the C-EM algorithm in large-scale practical problems.
The main technology for conducting high-throughput experiments in functional genomics is the microarray-a technical approach for assaying the abundance of mRNA for several genes simultaneously (see Hastie et al. (2001)  across samples. One way to organize gene expression data is to cluster genes on the basis of their expression patterns. One can think of the genes as points in ℜ 12 , which we want to cluster together in some fashion.
We fit multivariate normal mixtures to the yeast microarray data by finding the NPMLE of Q. This is an example of high-dimensional modeling. We observed similar performance of the algorithms to the previous examples. Table 4 presents the loglikelihood values. Since each observation is a point in ℜ 12 , at δ = 0.2, we obtain the empirical CDF as the MLE.
That is, each observation is its own component; hence the solution is not interesting.

Discussion
In this article we developed a framework for approximating the continuous parameter space and created an algorithm (based on the Penalized Dual method) for finding the maximum of l(Q); consequently an algorithm for estimating the mixture complexity. We established convergence properties of the proposed algorithm. By exploiting the inherent advantage of the penalty formulation, we derived a technique for converting the parameter estimators from the Penalized Dual problem into those for the mixture probability parameters. We established the existence of parameter estimators and derived convergence results for the Penalized Dual algorithm, for fitting overparameterized mixture models. It was shown empirically that the Penalized Dual algorithm has a faster rate of convergence, compared with the discrete EM algorithm for overparameterized mixture problems.
The algorithm based on the Penalized Dual method reaches closer to the global maximum and is robust to the choice of the support set Θ m ⊂ Ω (dimensionality of the problem).
These are desirable features for (1) analyzing high-dimensional data, and (2) for building a collection of semiparametric mixture models. The dimension of the dual optimization problem is fixed at d, the number of distinct observed data vectors; whereas that of the discrete EM grows with the cardinality of Θ m . For discrete mixture problems, such as binomial or Poisson, often d ≪ n; therefore, there is no dimensionality cost with the dual problem. When the cardinality of Θ m is large, the discrete EM algorithm fails to converge to the MLE, for all practical purposes, in certain mixture problems.
We derived several important structural properties of multivariate normal mixtures in which Q is modeled nonparametrically in the presence of an unknown variance-covariance matrix Σ ∈ S common to all m components. The role of the sieve parameter in reducing the dimension of the mixture problem was demonstrated by creating new graphical devices, namely the Mixture Tree plots.
The proposed methods are very powerful in searching over the whole discretized parameter space and in yielding a parsimonious mixture model. The discrete EM algorithm can be very difficult, if not impossible, in yielding a parsimonious model in problems with hundreds or thousands of parameters. Such problems are becoming increasingly common due to the rapid explosion of high-throughput data in microarray data and data mining. The applications for the methods described in this article are rich. Multivariate normal mixtures arise in many different practical scenarios, including data mining, knowledge discovery, data compression, pattern recognition and pattern classification.
for j = 1, . . . , m. From this definition of w i , the dual problem in (A.1) can be expressed as Equivalently, the problem is max w i n i log (w i ) subject to w ∈ ℜ d + which is the dual optimization problem in (7).
A.2 Derivation of the Penalized-Dual Estimator π ⋆ j,γ First, from the primal-dual relationship, it follows that for y i ∈ Y; i = 1, . . . , d. (A.5) Second, the following fixed-point equation is obtained by solving (12), By comparing the right-hand sides of (A.5) and (A.6), it is clear that g Q (y i ) parallels the term j {p θ j ( w γ )} (γ−1) f θ j (y i ) and that the latter expression resembles a mixture density with {p θ j ( w γ )} (γ−1) playing the role of π j .
A.3 Hessian Matrix of the Function K(z, γ) Let F = (f θ 1 , . . . , f θm ) T be an (m × d) matrix where f θ j = {f θ j (y 1 ), . . . , f θ j (y d )} T is the d-dimensional vector. In the sequel, we denote a vector of ones by 1 (with dimension clear from the context) and the diagonal matrix with elements x by diag(x). Therefore, K(z, γ) = 1 n n T · z − 1 γ 1 T · p γ for z ∈ ℜ and γ ∈ ℜ + , where n = (n 1 , . . . , n d ) T and the constraint vector p = p θ 1 (z), . . . , p θm (z) T with p θ j (z) (sometimes written as p j for exposition) is as in (11) expressed in terms of z ∈ ℜ.
The compactness of the parameter space Π can in turn be used to establish the convergence of g Q † γ to the maximizing value g Q . If the vector of masses π for g Q are uniquely determined, then the masses must converge as well. This in turn implies that as γ → ∞, the mixing distribution Q † γ with π † as the vector of masses is the NPMLE. Consequently, π † j,γ → π j as γ → ∞ for j = 1, . . . , m.
Proof of Theorem 6. The proof of part (a) is a consequence of the negative definiteness of H and the properties of P(w, γ) given in Proposition 1. For part (b), let z γ = arg max z∈ℜ K(z, γ) for any fixed γ ∈ ℜ + . That is, z γ is the maximizer of the K(z, γ) for a fixed γ. From equation (21) for a given z γ , it follows that where p j,γ is expressed in terms of z γ . From part (c) of Theorem 2, we have { p j,γ } γ ≤ 1 for j = 1, . . . , m. Therefore, log ( p j,γ ) ≤ 0 and hence ∂K( z γ , γ)/∂γ > 0 for any finite γ ∈ ℜ + and fixed z γ . That is, for a fixed z ∈ ℜ, the only point at which the function K(z, γ) can approach its supremum is at γ = ∞.
Proof of Theorem 8. For exposition we drop the subscript Σ from Q ⋆ and Q; however it is understood that the mixing measures have a dependence on the fixed Σ.
Part (1): We first establish that Q is an NPMLE. We start with creating a path in G from Q to Q ⋆ , by letting Q α = (1 − α) Q + α Q ⋆ for α ∈ [0, 1], Q ⋆ ∈ G and Q satisfying the relation (27). Note that Q α ∈ G; therefore, G is a convex set. The loglikelihood along this path satisfies l (Q α ; Σ) ≥ (1 − α) l Q; Σ + α l (Q ⋆ ; Σ) for α ∈ [0, 1] and for a fixed Σ. Therefore, l(Q α ; Σ) is a concave function for a fixed Σ. The directional derivative of l(Q α ; Σ) at g Q toward g Q ⋆ can be expressed, after simplification, as From (27), it follows that (A.9) is ≤ 0 for all Q ⋆ ∈ G and Q satisfying (27). This result combined with the concavity of l(Q α ; Σ) implies that Q is an NPMLE of Q ∈ G.