A Unified Framework for Compositional Fitting of Active Appearance Models

Active appearance models (AAMs) are one of the most popular and well-established techniques for modeling deformable objects in computer vision. In this paper, we study the problem of fitting AAMs using compositional gradient descent (CGD) algorithms. We present a unified and complete view of these algorithms and classify them with respect to three main characteristics: (i) cost function; (ii) type of composition; and (iii) optimization method. Furthermore, we extend the previous view by: (a) proposing a novel Bayesian cost function that can be interpreted as a general probabilistic formulation of the well-known project-out loss; (b) introducing two new types of composition, asymmetric and bidirectional, that combine the gradients of both image and appearance model to derive better convergent and more robust CGD algorithms; and (c) providing new valuable insights into existent CGD algorithms by reinterpreting them as direct applications of the Schur complement and the Wiberg method. Finally, in order to encourage open research and facilitate future comparisons with our work, we make the implementation of the algorithms studied in this paper publicly available as part of the Menpo Project (http://www.menpo.org).


Introduction
Active Appearance Models (AAMs) [15,29] are one of the most popular and well-established techniques for modeling and segmenting deformable objects in computer vision.AAMs are generative parametric models of shape and appearance that can be fitted to images to recover the set of model parameters that best describe a particular instance of the object being modeled.
Fitting AAMs is a non-linear optimization problem that requires the minimization (maximization) of a global error (similarity) measure between the input image and the appearance model.Several approaches [15,20,29,8,19,18,38,23,43,2,47,28,44,49,21,3] have been proposed to define and solve the previous optimization problem.Broadly speaking, they can be divided into two different groups: -Regression based [15, 20,8,18,43,47,44] -Optimization based [29,19,38,2,28,49,21] Regression based techniques attempt to solve the problem by learning a direct function mapping between the error measure and the optimal values of the parameters.Most notable approaches include variations on the original [15] fixed linear regression approach of [20,18], the adaptive linear regression approach of [8], and the works of [43] and [47] which considerably improved upon previous techniques by using boosted regression.Also, Cootes and Taylor [13] and Tresadern et al. [47] showed that the use of non-linear gradient-based and Haar-like appearance representations, respectively, lead to better fitting accuracy in regression based AAMs.
Optimization based methods for fitting AAMs were proposed by Matthews and Baker in [29].These techniques are known as Compositional Gradient Decent (CGD) algorithms and are based on direct analytical optimization of the error measure.Popular CGD algorithms include the very efficient project-out Inverse Compositional (PIC) algorithm [29], the accurate but costly Simultaneous Inverse Compositional (SIC) algorithm [19], and the more efficient versions of SIC presented in [38] and [49].Lucey et al. [25] extended these algorithms to the Fourier domain to efficiently enable convolution with Gabor filters, increasing their robustness; and the authors of [3] showed that optimization based AAMs using non-linear feature based (e.g.SIFT [24] and HOG [16]) appearance models were competitive with modern state-of-the-art techniques in non-rigid face alignment [55,4] in terms of fitting accuracy.
AAMs have often been criticized for several reasons: i) the limited representational power of their linear appearance model; ii) the difficulty of optimizing shape and appearance parameters simultaneously; and iii) the complexity involved in handling occlusions.However, recent works in this area [38,43,47,25,49,3] suggest that these limitations might have been over-stressed in the literature and that AAMs can produce highly accurate results if appropriate training data [49], appearance representations [47,25,3] and fitting strategies [38,43,47,49] are employed.
In this paper, we study the problem of fitting AAMs using CGD algorithms thoroughly.Summarizing, our main contributions are: -To present a unified and complete overview of the most relevant and recently published CGD algorithms for fitting AAMs [29,19,38,2,28,51,49,21].
To this end, we classify CGD algorithms with respect to three main characteristics: i) the cost function defining the fitting problem; ii) the type of composition used; and iii) the optimization method employed to solve the non-linear optimization problem.-To review the probabilistic interpretation of AAMs and propose a novel Bayesian formulation2 of the fitting problem.We assume a probabilistic model for appearance generation with both Gaussian noise and a Gaussian prior over a latent appearance space.Marginalizing out the latent appearance space, we derive a novel cost function that only depends on shape parameters and that can be interpreted as a valid and more general probabilistic formulation of the well-known project-out cost function [29].Our Bayesian formulation is motivated by seminal works on probabilistic component analysis and object tracking [34,40,46].-To propose the use of two novel types of composition for AAMs: i) asymmetric; and ii) bidirectional.
These types of composition have been widely used in the related field of parametric image alignment [27,32,5,33] and use the gradients of both image and appearance model to derive better convergent and more robust CGD algorithms.-To provide valuable insights into existent strategies used to derive fast and exact simultaneous algorithms for fitting AAMs by reinterpreting them as direct applications of the Schur complement [11] and the Wiberg method [37,45].
The remainder of the paper is structured as follows.Section 2 introduces AAMs and reviews their probabilistic interpretation.Section 3 constitutes the main section of the paper and contains the discussion and derivations related to the cost functions 3.1; composition types 3.2; and optimization methods 3.3.Implementation details and experimental results are reported in Section 5. Finally, conclusions are drawn in Section 6.

Active Appearance Models
AAMs [15,29] are generative parametric models that explain visual variations, in terms of shape and appearance, within a particular object class.AAMs are built from a collection of images for which the spatial position of a sparse set of v landmark points x i = (x i , y i ) T ∈ R 2 representing the shape s = (x 1 , y 1 , . . ., x v , y v ) T ∈ R 2v×1 of the object being modeled have been manually defined a priori.
AAMs are themselves composed of three different models: (i) shape model; (ii) appearance model; and (iii) motion model.
The shape model, which is also referred to as Point Distribution Model (PDM), is obtained by typically applying Principal Component Analysis (PCA) to the set of object's shapes.The resulting shape model is mathematically expressed as: where s ∈ R 2v×1 is the mean shape, and S ∈ R 2v×n and p ∈ R n×1 denote the shape bases and shape parameters, respectively.In order to allow a particular shape instance s to be arbitrarily positioned in space, the previous model can be augmented with a global similarity transform.Note that this normally requires the initial shapes to be normalized with respect to the same type of transform (typically using Procrustes Analysis (PA)) Fig. 1: Exemplar images from the Labeled Faces in-the-Wild (LFPW) dataset [9] for which a consistent set of sparse landmarks representing the shape of the object being model (human face) has been manually defined [41,42].before PCA is applied.This results in the following expression for each landmark point of the shape model: where s, R ∈ R 2×2 and t ∈ R 2 denote the scale, rotation and translation applied by the global similarity transform, respectively.Using the orthonormalization procedure described in [29] the final expression for the shape model can be compactly written as the linear combination of a set of bases: where S = (s * 1 , . . ., s * 4 , s 1 , • • • , s n ) ∈ R 2v×(n+4) and p = (p * 1 , . . ., p * 4 , p 1 , . . ., p n ) T ∈ R (n+4)×1 are redefined as the concatenation of the similarity bases s * i and similarity parameters p * i with the original S and p, respectively.The appearance model is obtained by warping the original images onto a common reference frame (typically defined in terms of the mean shape s) and applying PCA to the obtained warped images.Mathematically, the appearance model is defined by the following expression: where x ∈ Ω denote all pixel positions on the reference frame, and Ā(x), A i (x) and c i denote the mean texture, the appearance bases and appearance parameters, respectively.Denoting a = vec(A(x)) as the vectorized version of the previous appearance instance, Equation 4 can be concisely written in vector form as: where a ∈ R F ×1 is the mean appearance, and A ∈ R F ×m and c ∈ R m×1 denote the appearance bases and appearance parameters, respectively.The role of the motion model, denoted by W(x; p), is to extrapolate the position of all pixel positions x ∈ Ω from the reference frame to a particular shape instance s (and vice-versa) based on their relative position with respect to the sparse set of landmarks defining the shape model (for which direct correspondences are always known).Classic motion models for AAMs are PieceWise Affine (PWA) [14,29] and Thin Plate Splines (TPS) [14,38] warps.
Given an image I containing the object of interest, its manually annotated ground truth shape s, and a particular motion model W(x, p); the two main assumptions behind AAMs are: 1.The ground truth shape of the object can be well approximated by the shape model 2. The object's appearance can be well approximated by the appearance model after the image is warped, using the motion model and the previous shape approximation, onto the reference frame: where i[p] = vec(I(W(x; p))) denotes the vectorized version of the warped image.Note that, the warp W(x; p) which explicitly depends on the shape parameters p, relates the shape and appearance models and is a central part of the AAMs formulation.
Because of the explicit use of the motion model, the two previous assumptions provide a concise definition of AAMs.At this point, it is worth mentioning that the vector notation of Equations 6 and 7 will be, in general, the preferred notation in this paper.

Probabilistic Formulation
A probabilistic interpretation of AAMs can be obtained by rewriting Equations 6 and 7 assuming probabilistic models for shape and appearance generation.In this paper, motivated by seminal works on Probabilistic Component Analysis (PPCA) and object tracking [46,40,34], we will assume probabilistic models for shape and appearance generation with both Gaussian noise and Gaussian priors over the latent shape and appearance spaces3 : where the diagonal matrices Λ = diag(λ s1 , • • • , λ sm ) and Σ = diag(λ a1 , • • • , λ am ) contain the eigenvalues associated to shape and appearance eigenvectors respectively and where ς 2 and σ 2 denote the estimated shape and image noise4 respectively.
This probabilistic formulation will be used to derive Maximum-Likelihood (ML), Maximum A Posteriori (MAP) and Bayesian cost functions for fitting AAMs in Sections 3.1.1and 3.1.2.
The following subsections present a unified and complete view of CGD algorithms by classifying them with respect to their three main characteristics: a) cost function (Section 3.1); b) type of composition (Section 3.2); and c) optimization method (Section 3.3).

Cost Function
AAM fitting is typically formulated as the (regularized) search over the shape and appearance parameters that minimize a global error measure between the vectorized warped image and the appearance model: where D is a data term that quantifies the global error measure between the vectorized warped image and the appearance model and R is an optional regularization term that penalizes complex shape and appearance deformations.

Sum of Squared Differences
Arguably, the most natural choice for the previous data term is the Sum of Squared Differences (SSD) between the vectorized warped image and the linear appearance model5 .Consequently, the classic AAM fitting problem is defined by the following non-linear optimization problem6 : On the other hand, considering regularization, the most natural choice for R is the sum of 2 2 -norms over the shape and appearance parameters.In this case, the regularized AAM fitting problem is defined as follows:

Probabilistic Formulation
A probabilistic formulation of the previous cost function can be naturally derived using the probabilistic generative models introduced in Section 2.1.Denoting the models' parameters as Θ = {s, S, Λ, ā, A, Σ, σ 2 } a ML formulation can be derived as follows: and a MAP formulation can be similarly derived by taking into account the prior distributions over the shape and appearance parameters: where we have assumed the shape and appearance parameters to be independent7 .The previous ML and MAP formulations are weighted version of the optimization problem defined by Equation 11 and 12.In both cases, the maximization of the conditional probability of the vectorized warped image given the shape, appearance and model parameters leads to the minimization of the data term D and, in the MAP case, the maximization of the prior probability over the shape and appearance parameters leads to the minimization of the regularization term R.

Project-Out
Matthews and Baker showed in [29] that one could express the SSD between the vectorized warped image and the linear PCA-based8 appearance model as the sum of two different terms: The first term defines the distance within the appearance subspace and it is always 0 regardless of the value of the shape parameters p: The second term measures the distance to the appearance subspace i.e. the distance within its orthogonal complement.After some algebraic manipulation, one can show that this term reduces to a function that only depends on the shape parameters p: where, for convenience, we have defined the orthogonal complement to the appearance subspace as Ā = I − AA T .Note that, as mentioned above, the previous term does not depend on the appearance parameters c: Therefore, using the previous project-out trick, the minimization problems defined by Equations 11 and 12 reduce to: and respectively.

Probabilistic Formulation
Assuming the probabilistic models defined in Section 2.1, a Bayesian formulation of the previous project-out data term can be naturally derived by marginalizing over the appearance parameters to obtain the following marginalized density: and applying the Woodbury formula 9 [54] to decompose the natural logarithm of the previous density into the sum of two different terms: where As depicted by Figure 2, the previous two terms define respectively: i) the Mahalanobis distance within the linear appearance subspace; and ii) the Euclidean distance to the linear appearance subspace (i.e. the Euclidean distance within its orthogonal complement) weighted by the inverse of the estimated image noise.Note that when the variance Σ of the prior distribution over the latent appearance space increases (and especially as Σ → ∞) c becomes uniformly distributed and the contribution of the first term Fig. 2: The Bayesian project-out formulation fits AAMs by minimizing two different distances: i) the Mahalanobis distance within the linear appearance subspace; and ii) the Euclidean distance to the linear appearance subspace (i.e. the Euclidean distance within its orthogonal complement) weighted by the inverse of the estimated image noise.
loss arises naturally by assuming a uniform prior over the latent appearance space.The probabilistic formulations of the minimization problems defined by Equations 19 and 20 can be derived, from the previous Bayesian Project-Out (BPO) cost function, as and respectively.Where we have defined the BPO operator as

Type of Composition
Assuming, for the time being, that the true appearance parameters c * are known, the problem defined by Equation 11 reduces to a non-rigid image alignment problem [7,35] between the particular instance of the object present in the image and its optimal appearance reconstruction by the appearance model: where a = ā + Ac * is obtained by directly evaluating Equation 4 given the true appearance parameters c * .CGD algorithms iteratively solve the previous nonlinear optimization problem with respect to the shape parameters p by: 1. Introducing an incremental warp W(x; ∆p) according to the particular composition scheme being used.
2. Linearizing the previous incremental warp around the identity warp W(x; ∆p) = W(x; 0) = x.3. Solving for the parameters ∆p of the incremental warp.4. Updating the current warp estimate by using an appropriate compositional update rule. 5. Going back to Step 1 until a particular convergence criterion is met.
Existent CGD algorithms for fitting AAMs have introduced the incremental warp either on the image or the model sides in what are known as forward and inverse compositional frameworks [29,19,38,2,28,49] respectively.Inspired by related works in field of image alignment [27,32,5,33], we notice that novel CGD algorithms can be derived by introducing incremental warps on both image and model sides simultaneously.Depending on the exact relationship between these incremental warps we define two novel types of composition: asymmetric and bidirectional.
The following subsections explain how to introduce the incremental warp into the cost function and how to update the current warp estimate for the four types of composition considered in this paper: i) forward; ii) inverse; iii) asymmetric; and iv) bidirectional.For convenience, in these subsections we will use the simplified cost function defined by Equation 25.Furthermore, to maintain consistency with the vector notation used through out the paper, we will abuse the notation and write the operations of warp composition10 W(x; p) • W(x; ∆p) and inversion 10 W(x; p) −1 as simply p • ∆p and p −1 respectively.

Forward
In the forward compositional framework the incremental warp ∆p is introduced on the image side at each iteration by composing it with the current warp estimate p: Once the optimal values for the parameters of the incremental warp are obtained, the current warp estimate is updated according to the following compositional update rule:

Inverse
On the other hand, the inverse compositional framework inverts the roles of the image and the model by introducing the incremental warp on the model side: Note that, in this case, the model is the one we seek to deform using the incremental warp.
Because the incremental warp is introduced on the model side, the solution ∆p needs to be inverted before it is composed with the current warp estimate:

Asymmetric
Asymmetric composition introduces two related incremental warps onto the cost function; one on the image side (forward) and the other on the model side (inverse): Note that the previous two incremental warps are defined to be each others inverse.Consequently, using the first order approximation to warp inversion for typical AAMs warps ∆p −1 = −∆p defined in [29], we can rewrite the previous asymmetric cost function as: Although this cost function will need to be linearized around both incremental warps, the parameters ∆p controlling both warps are the same.Also, note that the parameters α ∈ [0, 1] and β = (1 − α) control the relative contribution of both incremental warps in the computation of the optimal value for ∆p.
In this case, the update rule for the current warp estimate is obtained by combining the previous forward and inverse compositional update rules into a single compositional update rule: Note that, the special case in which α = β = 0.5 is also referred to as symmetric composition [32,5,33] and that the previous forward and inverse compositions can also be obtained from asymmetric composition by setting α = 1 , β = 0 and α = 0 , β = 1 respectively.

Bidirectional
Similar to the previous asymmetric composition, bidirectional composition also introduces incremental warps on both image and model sides.However, in this case, the two incremental warps are assumed to be independent from each other: Consequently, in Step 4, the cost function needs to be linearized around both incremental warps and solved with respect to the parameters controlling both warps, ∆p and ∆q.
Once the optimal value for both sets of parameters is recovered, the current estimate of the warp is updated using:

Optimization Method
Step 2 and 3 in CGD algorithms, i.e. linearizing the cost and solving for the incremental warp respectively, depend on the specific optimization method used by the algorithm.In this paper, we distinguish between three main optimization methods 11 : i) Gauss-Newton [11, 29,19,28,38,49]; ii) Newton [11,21]; and iii) Wiberg [37,45,38,49].These methods can be used to iteratively solve the non-linear optimization problems defined by Equations 14 and 22.The main differences between them are: 1.The term being linearized.Gauss-Newton and Wiberg linearize the residual r while Newton linearizes the whole data term D.
11 Amberg et al. proposed the use of the Steepest Descent method [11] in [2].However, their approach requires a special formulation of the motion model and it performs poorly using the standard independent AAM formulation [29] used in this work.
2. The way in which each method solves for the incremental parameters ∆c, ∆p and ∆q.Gauss-Newton and Newton can either solve for them simultaneously or in an alternated fashion while Wiberg defines its own procedure to solve for different sets of parameters12 .
The following subsections thoroughly explain how the previous optimization methods are used in CGD algorithms.In order to simplify their comprehension full derivations will be given for all methods using the SSD data term (Equation 11) with both asymmetric (Section 3.2.3)and bidirectional (Section 3.2.4)compositions13 while only direct solutions will be given for the Project-Out data term (Equation 19).

Gauss-Newton
When asymmetric composition is used, the optimization problem defined by the SSD data term is: with the asymmetric residual r a defined as: and where we have introduced the incremental appearance parameters ∆c14 .The Gauss-Newton method solves the previous optimization problem by performing a first order Taylor expansion of the residual: and solving the following approximation of the original problem: where, in order to unclutter the notation, we have defined ∆ = (∆c T , ∆p T ) T and the partial derivative of the residual with respect to the previous parameters, i.e. the Jacobian of the residual, is defined as: where ∇t = (α∇i[p] + β∇(a + Ac)).
When bidirectional composition is used, the optimization problem is defined as: where the bidirectional residual r b reduces to: The Gauss-Newton method proceeds in exactly the same manner as before, i.e. performing a first order Taylor expansion: and solving the approximated problem: where, in this case, ∆ = (∆c T , ∆p T , ∆q T ) T and the Jacobian of the residual is defined as: where

Simultaneous
The optimization problem defined by Equations 38 and 43 can be solved with respect to all parameters simultaneously by simply equating their derivative to 0: The solution is given by: where ∂r

∂∆
T ∂r ∂∆ is known as the Gauss-Newton approximation to the Hessian matrix.
Directly inverting ∂r

∂∆
T ∂r ∂∆ has complexity 15 O((n+ m) 3 ) for asymmetric composition and O((2n + m) 3 ) for bidirectional composition.However, one can take advantage of the problem structure and derive an algorithm with smaller complexity by using the Schur complement 16 [11].
For asymmetric composition we have: Applying the Schur complement, the solution for ∆p is given by: and plugging the solution for ∆p into equation 47 the optimal value for ∆c is obtained by: Using the above procedure the complexity 15 of solving each Gauss-Newton step is reduced to: O(nmF Using bidirectional composition, we can apply the Schur complement either one or two times in order to 15 m and n denote the number of shape and appearance parameters respectively while F denotes the number of pixels on the reference frame. 16Applying the Schur complement to the following system of equations: the solution for x is given by: and the solution for y is obtained by substituting the value of x into the original system.take advantage of the 3×3 block structure of the matrix Applying the Schur complement once, the combined solution for (∆p T , ∆q T ) T is given by: Note that the complexity of inverting this new approximation to the Hessian matrix is O((2n) 3 )17 .Similar to before, plugging the solutions for ∆p and ∆q into Equation 51 we can infer the optimal value for ∆c using: The total complexity per iteration of the previous approach is: The Schur complement can be re-applied to Equation 52 to derive a solution for ∆q that only requires inverting a Hessian approximation matrix of size n × n: where we have defined the projection matrix P as: and the solutions for ∆p and ∆c can be obtained by plugging the solutions for ∆q into Equation 52 and the solutions for ∆q and ∆p into Equation 51 respectively: The total complexity per iteration of the previous approach reduces to: Note that because of their reduced complexity, the solutions defined by Equations 55 and 57 are preferred over the ones defined by Equations 52 and 53.Finally, the solutions using the Project-Out cost function are: -For asymmetric composition: with complexity 18 given by Equation 50.-For bidirectional composition: with complexity 18 given by Equation 58.

Alternated
Another way of solving optimization problems with two or more sets of variables is to use alternated optimization [17].Hence, instead of solving the previous problem simultaneously with respect to all parameters, we can update one set of parameters at a time while keeping the other sets fixed.More specifically, using asymmetric composition we can alternate between updating ∆c given the previous ∆p and then update ∆p given the updated ∆c in an alternate manner.Taking advantage of the structure of the problem defined by Equation 47, we can obtain the following system of equations: 18 In practice, the solutions for the Project-Out cost function can be computed slightly faster than those for the SSD because they do not need to explicitly solve for ∆c.This is specially important in the inverse compositional case because expressions of the form (J T UJ) −1 J T U can be completely precomputed and the computational cost per iteration reduces to O(nF ).
which we can rewrite as: in order to obtain the analytical expression for the previous alternated update rules.The complexity at each iteration is dominated by: In the case of bidirectional composition we can proceed in two different ways: a) update ∆c given the previous ∆p and ∆q and then update (∆p T , ∆q T ) T from the updated ∆c, or b) update ∆c given the previous ∆p and ∆q, then ∆p given the updated ∆c and the previous ∆q and, finally, ∆q given the updated ∆c and ∆p.
From Equation 51, we can derive the following system of equations: from which we can define the alternated update rules for the first of the previous two options: with complexity: The rules for the second options are: and their complexity is dominated by: On the other hand, the alternated update rules using the Project-Out cost function are: -For asymmetric composition: There is no proper alternated rule because the Project-Out cost function only depends on one set of parameters, ∆p.-For bidirectional composition: with equivalent complexity to the one given by Equation 50 because, in this case, the term J T ā ĀJ ā −1 J T ā Ā can be completely precomputed.
Note that all previous alternated update rules, Equations 62, 65, 67 and 98, are similar but slightly different from their simultaneous counterparts, Equations 48 and 49, 52 and 53, 55 and 57, and 60.

Newton
The Newton method performs a second order Taylor expansion of the entire data term D: and solves the approximate problem: Assuming asymmetric composition, the previous data term is defined as: and the matrix containing the first order partial derivatives with respect to the parameters, i.e. the data term's Jacobian, can be written as: On the other hand, the matrix ∂ 2 Da ∂ 2 ∆ of the second order partial derivatives, i.e. the Hessian of the data term, takes the following form: Note that the Hessian matrix is, by definition, symmetric.The definition of its individual terms is provided in Appendix A.1.
A similar derivation can be obtained for bidirectional composition where, as expected, the data term is defined as: In this case, the Jacobian matrix becomes: and the Hessian matrix takes the following form: Notice that the previous matrix is again symmetric.The definition of its individual terms is provided in Appendix A.2.

Simultaneous
Using the Newton method we can solve for all parameters simultaneously by equating the partial derivative of Equation 71 to 0: with the solution given by: Note that, similar to the Gauss-Newton method, the complexity of inverting the Hessian matrix ) for asymmetric composition and O((2n + m) 3 ) for bidirectional composition.As shown by Kossaifi et al. [21] 19 , we can take advantage of the structure of the Hessian in Equations 74 and 77 and apply the Schur complement to obtain more efficient solutions.
The solutions for ∆p and ∆c using asymmetric composition are given by the following expressions: with complexity: where we have defined in order to unclutter the notation.
On the other hand, the solutions for bidirectional composition are given either by: or where we have defined the following auxiliary matrices The complexity of the previous solutions is of: and respectively.
The solutions using the Project-Out cost function are: -For asymmetric composition: with complexity 20 given by Equation 81.-For bidirectional composition: where the projection operator P is defined as: and where we have defined: to unclutter the notation.The complexity per iteration 20 is given by Equation 87.
Note that, the derivations of the previous solutions, for both types of composition, are analogous to the ones shown in Section 3.3.1 for the Gauss-Newton method and, consequently, have been omitted here. 20In practice, the solutions for the project-out cost function can also be computed slightly faster because they do not need to explicitly solve for ∆c.However, in this case, using inverse composition we can only precompute terms of the form J T U and J T UJ but not the entire H −1 J T U because of the explicit dependence between H and the current residual r.

Alternated
Alternated optimization rules can also be derived for the Newton method following the strategy shown in Section 3.3.1 for the Gauss-Newton case.Again, we will simply provide update rules and computational complexity for both types of composition and will omit the details of their full derivation.
For asymmetric composition the alternated rules are defined as: with complexity: The alternated rules for bidirectional composition case are given either by: with complexity: or: with complexity: On the other hand, the alternated update rules for the Newton method using the project-out cost function are: -For asymmetric composition: Again, there is no proper alternated rule because the project-out cost function only depends on one set of parameters, ∆p.-For bidirectional composition: where we have defined: and the complexity at every iteration is given by the following expression complexity: O(nmF Efficient Second-order Minimization (ESM) Notice that, the second order Taylor expanison used by the Newton method means that Newton algorithms are second order optimizations algorithms with respect to the incremental warps.However, as shown in the previous section, this property comes at expenses of a significant increase in computational complexity with respect to (first order)Gauss-Newton algorithms.In this section, we show that the Asymmetric Gauss-Newton algorithms derived in Section 3.3.1 are, in fact, also true second order optimization algorithms with respect to the incremental warp ∆p.The use of asymmetric composition together with the Gauss-Newton method has been proven to naturally lead to Efficient Second order Minimization (ESM) algorithms in the related field of parametric image alignment [27,10,32,33].Following a similar line of reasoning, we will show that Asymmetric Gauss-Newton algorithms for fitting AAMs can also be also interpreted as ESM algorithms.
In order to show the previous relationship we will make use of the simplified data term 21 introduced by 21 Notice that similar derivations can also be obtained using the SSD and Project-Out data terms, but we use the simplified one here for clarity.
Equation 25.Using forward composition, the optimization problem defined by: where the forward residual r f is defined as: As seen before, Gauss-Newton solves the previous optimization problem by performing a first order Taylor expansion of the residual around ∆p: and solving the following approximation of the original problem: However, note that, instead of performing a first order Taylor expansion, we can also perform a second order Taylor expansion of the residual: Then, given the second main assumption behind AAMs (Equation 7) the following approximation must hold: and, because the previous J i and J a are functions of ∆p, we can perform a first order Taylor expansion of J i to obtain: Finally, substituting the previous approximation for ∆p T H i into Equation 105 we arrive at: where the total remainder is cubic with respect to ∆p: The previous expression constitutes a true second order approximation of the forward residual r f where the term 1 2 (J i + J a ) is equivalent to the asymmetric Jacobian in Equation 38 when α = β = 0.5: and, consequently, Asymmetric Gauss-Newton algorithms for fitting AAMs can be viewed as ESM algorithms that only require first order partial derivatives of the residual and that have the same computational complexity as first order algorithms.

Wiberg
The idea behind the Wiberg method is similar to the one used by the alternated Gauss-Newton method in Section 3.3.1,i.e. solving for one set of parameters at a time while keeping the other sets fixed.However, Wiberg does so by rewriting the asymmetric r a (∆c, ∆p) and bidirectional r b (∆c, ∆p, ∆q) residuals as functions that only depend on ∆p and ∆q respectively.
For asymmetric composition, the residual ra (∆p) is defined as follows: where the function ∆c a (∆p) is obtained by solving for ∆c while keeping ∆p fixed: Given the previous residual, the Wiberg method proceeds to define the following optimization problem with respect to ∆p: In this case, the Jacobian ∂r ∂∆p can be obtain by direct application of the chain rule and it is defined as follows: The solution for ∆p is obtained as usual by equating the derivative of 113 with respect to ∆p to 0: where we have used the fact that the matrix Ā is idempotent 22 .Therefore, the Wiberg method solves explicitly, at each iteration, for ∆p using the previous expression and implicitly for ∆c (through ∆c a (∆p)) using Equation 112.The complexity per iteration of the Wiberg method is the same as the one of the Gauss-Newton method after applying the Schur complement, Equation 50.In fact, note that the Wiberg solution for ∆p (Equation 116) is the same as the one of the Gauss-Newton method after applying the Schur complement, Equation 48; and also note the similarity between the 22 Ā is idempotent: solutions for ∆c of both methods, Equations 112 and 49.Finally, note that, due to the close relation between the Wiberg and Gauss-Newton methods, Asymmetric Wiberg algorithms are also ESM algorithms for fitting AAMs.
On the other hand, for bidirectional composition, the residual rb (∆p) is defined as: where, similarly as before, the function ∆c b (∆p, ∆q) is obtained solving for ∆c while keeping both ∆p and ∆q fixed: and the function ∆p b ( ∆c b , ∆q) is obtained by solving for ∆p using the Wiberg method while keeping ∆q fixed: At this point, the Wiberg method proceeds to define the following optimization problem with respect to ∆q: In this case, the Jacobian of the residual can also be obtained by direct application of the chain rule and takes the following form: And, again, the solution for ∆q is obtained as usual by equating the derivative of 121 with respect to ∆q to 0: In this case, the Wiberg method solves explicitly, at each iteration, for ∆p using the previous expression and implicitly for ∆p and ∆c (through ∆p b ( ∆c b , ∆q) and ∆c b (∆p, ∆q)) using Equations 119 and 118 respectively.Again, the complexity per iteration is the same as the one of the Gauss-Newton method after applying the Schur complement, Equation 58; and the solutions for both methods are almost identical, Equations 123, 119 and 118 and Equations 52, 53 and 55.
On the other hand, the Wiberg solutions for the project-out cost function are: -For asymmetric composition: Because the projectout cost function only depends on one set of parameters, ∆p, in this case Wiberg reduces to Gauss-Newton.
-For bidirectional composition: Again, in this case, the solutions obtained with the Wiberg method are almost identical to the ones obtained using Gauss-Newton after applying the Schur complement, Equation 60.

Relation to Prior Work
In this section we relate relevant prior work on CGD algorithms for fitting AAMs [29,19,38,2,28,49,21] to the unified and complete view introduced in the previous Section.

Project-Out algorithms
In their seminal work [29], Matthews and Baker proposed the first CGD algorithm for fitting AAMs, the so-called Project-out Inverse Compositional (PIC) algorithm.This algorithm uses Gauss-Newton to solve the optimization problem posed by the project-out cost function using inverse composition.The use of the projectout norm removes the need to solve for the appearance parameters and the use of inverse composition allows for the precomputation of the pseudo-inverse of the Jacobian with respect to ∆p, i.e.J T ā ĀJ ā −1 J ā Ā.The PIC algorithm is very efficient (O(nF )) but it has been shown to perform poorly in generic and unconstrained scenarios [19,38].In this paper, we refer to this algorithm as the Project-Out Inverse Gauss-Newton algorithm.
The forward version of the previous algorithm, i.e. the Project-Out Forward Gauss-Newton algorithm, was proposed by Amberg et al. in [2].In this case, the use of forward composition prevents the precomputation of the Jacobian pseudo-inverse and its complexity increases to O(nmF +n 2 F +n 3 ).However, this algorithm has been shown to largely outperform its inverse counterpart, and obtains good performance under generic and unconstrained conditions [2, 49] 23  To the best of our knowledge, the rest of Project-Out algorithms derived in Section 3, i.e.: -Project-Out Forward Newton -Project-Out Inverse Newton -Project-Out Asymmetric Gauss-Newton -Project-Out Asymmetric Newton -Project-Out Bidirectional Gauss-Newton Schur -Project-Out Bidirectional Gauss-Newton Alternated -Project-Out Bidirectional Newton Schur -Project-Out Bidirectional Newton Alternated -Project-Out Bidirectional Wiberg have never been published before and are a significant contribution of this work.

SSD algorithms
In [19] Gross et al. presented the Simultaneous Inverse Compositional (SIC) algorithm and show that it largely outperforms the Project-Out Inverse Gauss-Newton algorithm in terms of fitting accuracy.This algorithm uses Gauss-Newton to solve the optimization problem posed by the SSD cost function using inverse composition.In this case, the Jacobian with respect to ∆p, depends on the current value of the appearance parameters and needs to be recomputed at every iteration.Moreover, the inclusion of the Jacobian with respect to the appearance increments δc, increases the size of the simultaneous Jacobian to ∂r ∂∆ = (−A, −J a ) ∈ R F ×(m+n) and, consequently, the computational cost per iteration of the algorithm is As we shown in Sections 3.3.1,3.3.1 and 3.3.3 the previous complexity can be dramatically reduced by taking advantage of the problem structure in order to derive more efficient and exact algorithm by: a) applying the Schur complement; b) adopting an alternated optimization approach; or c) or using the Wiberg 23 Notice that, in [2], Amberg et al. also introduced a hybrid forward/inverse algorithm, coined CoLiNe.This algorithm is a compromise between the previous two algorithms in terms of both complexity and accuracy.Due to its rather ad-hoc derivation, this algorithm was not considered in this paper.method.Papandreou and Maragos [38] proposed an algorithm that is equivalent to the solution obtained by applying the Schur complement to the problem, as described in Section 3.3.1.The same algorithm was reintroduced in [49] using a somehow ad-hoc derivation (reminiscent of the Wiberg method) under the name Fast-SIC.This algorithm has a computational cost per iteration of O(nmF + n 2 F + n 3 ).In this paper, following our unified view on CGD algorithm, we refer to the previous algorithm as the SSD Inverse Gauss-Newton Schur algorithm.The alternated optimization approach was used in [51] and [3] with complexity O(n 2 F + n 3 ) per iteration.We refer to it as the SSD Inverse Gauss-Newton Alternated algorithm.
On the other hand, the forward version of the previous algorithm was first proposed by Martins et al. in [28] 24 .In this case, the Jacobian with respect to ∆p depends on the current value of the shape parameters p through the warped image i[p] and also needs to be recomputed at every iteration.Consequently, the complexity if the algorithm is the same as in the naive inverse approach of Gross et al.In this paper, we refer to this algorithm as the SSD Forward Gauss-Newton algorithm.It is important to notice that Tzimiropoulos and Pantic [49] derived a more efficient version of this algorithm (O(nmF + n 2 F + n 3 )), coined Fast-Forward, by applying the same derivation used to obtain their Fast-SIC algorithm.They showed that in the forward case their derivation removed the need to explicitly solve for the appearance parameters.Their algorithm is equivalent to the previous Project-Out Forward Gauss-Newton.
Finally, Kossaifi et al. derived the SSD Inverse Newton Schur algorithm in [21].This algorithm has a total complexity per iteration of O(nmF + n 2 m + 2n 2 F + n 3 ) and was shown to slightly underperform its equivalent Gauss-Newton counterpart.
The remaining SSD algorithms derived in Section 3, i.e.: -SSD Inverse Wiberg -SSD Forward Gauss-Newton Alternated -SSD Forward Newton Schur -SSD Forward Newton Alternated -SSD Forward Wiberg -SSD Asymmetric Gauss-Newton Schur -SSD Asymmetric Gauss-Newton Alternated -SSD Asymmetric Newton Schur -SSD Asymmetric Newton Alternated -SSD Asymmetric Wiberg -SSD Bidirectional Gauss-Newton Schur -SSD Bidirectional Gauss-Newton Alternated -SSD Bidirectional Newton Schur -SSD Bidirectional Newton Alternated -SSD Bidirectional Wiberg have never been published before and are also a key contribution of the presented work.Notice that, the iterative solutions of all CGD algorithms studied in this paper are given in Appendix B.

Experiments
In this section, we analyze the performance of the CGD algorithms derived in Section 3 on the specific problems of non-rigid face alignment in-the-wild.Results for five experiments are reported.The first experiment compares the fitting accuracy and convergence properties of all algorithms on the test set of the popular Labeled Faces in-the-Wild (LFPW) [9] database.The second experiment quantifies the importance of the two terms in the Bayesian project-out cost function in relation to the fitting accuracy obtained by Project-Out algorithms.In the third experiment, we study the effect that varying the value of the parameters α and β has on the performance of Asymmetric algorithms.The fourth experiment explores the effect of optimizing the cost functions using reduced subsets of the total number of pixels and quantifies the impact that this has on the accuracy and computational efficiency of CGD algorithms.Finally, in the fifth experiment, we report the performance of the most accurate CGD algorithms on the test set of the Helen [22] database and on the entire Annotated Faces in-the-Wild (AFW) [56] database.
Throughout this section, we abbreviate CGD algorithms using the following convention: CF TC OM( OS) where: a) CF stands for Cost Function and can be either SSD or PO depending on whether the algorithm uses the Sum of Squared Differences or the Project Out cost function; b) TC stands for Type of Composition and can be For, Inv, Asy or Bid depending on whether the algorithm uses Forward, Inverse, Asymmetric or Bidirectional compositions; c) OM stands for Optimization Method and can be GN, N or W depending on whether the algorithm uses the Gauss-Newton, Newton or Wiberg optimization methods; and, finally, d) if Gauss-Newton or Newton methods are used, the optional field OS, which stands for Optimization Strategy, can be Sch or Alt depending on whether the algorithm solves for the parameters simultaneously using the Schur complement or using Alternated optimization.For example, following the previous convention the Project Out Bidirectional Gauss-Newton Schur algorithm is denoted by PO Bid GN Sch.
Landmark annotations for all databases are provided by the iBUG group25 [41,42] and fitting accuracy is reported using the point-to-point error measure normalized by the face size 26 proposed in [56] over the 49 interior points of the iBug annotation scheme.
In all face alignment experiments, we use a single AAM, trained using the ∼ 800 and ∼ 2000 training images of the LFPW and Helen databases.Similar to [50], we use a modified version of the Dense Scale Invariant Feature Transform (DSIFT) [24,16] to define the appearance representation of the previous AAM.In particular, we describe each pixel with a reduced SIFT descriptor of length 8 using the public implementation provided by the authors of [52].All algorithms are implemented in a coarse to fine manner using a Gaussian pyramid with 2 levels (face images are normalized to a face size 26 of roughly 150 pixels at the top level).In all experiments, we optimize over 7 shape parameters (4 similarity transform and 3 non-rigid shape parameters) at the first pyramid level and over 16 shape parameters (4 similarity transform and 12 non-rigid shape parameters) at the second one.The dimensionality of the appearance models is kept to represent 75% of the total variance in both levels.This results in 225 and 280 appearance parameters at the first and second pyramid levels respectively.The previous choices were determined by testing on a small hold out set of the training data.
In all experiments, algorithms are initialized by perturbing the similarity transform that perfectly aligns the model's mean shape (a frontal pose and neutral expression looking shape) with the ground truth shape of each image.These transforms are perturbed by adding uniformly distributed random noise to their scale, rotation and translation parameters.Exemplar initializations obtained by this procedure for different amounts of noise are shown in Figure 3. Notice that, we found that initializing using 5% uniform noise is (statistically) equivalent to initializing with the popular OpenCV [12] implementation of the well-known Viola and Jones face detector [53] on the test images of the LFPW database.
Unless stated otherwise: i) algorithms are initialized with 5% uniform noise ii) test images are fitted three times using different random initializations (the same exact random initializations are used for all algorithms); iii) algorithms are left to run for 40 iterations (24 iterations at the first pyramid level and 16 at the second); iv) results for Project-Out algorithms are obtained using the Bayesian project-out cost function defined by Equation 22; and v) results for Asymmetric algorithms are reported for the special case of symmetric composition i.e. α = β = 0.5 in Equation 30.
In order to showcase the broader applicability of AAMs, we complete the previous performance analysis by performing a sixth and last experiment on the problem of non-rigid car alignment in-the-wild.To this end, we report the fitting accuracy of the best performing CGD algorithms on the MIT StreetScene 27 database.
Finally, in order to encourage open research and facilitate future comparisons with the results presented in this section, we make the implementation of all algorithms publicly available as part of the Menpo Project 1 [1].

Comparison on LFPW
In this experiment, we report the fitting accuracy and convergence properties of all CGD algorithms studied in this paper.Results are reported on the ∼ 220 test images of the LFPW database.In order to keep the information easily readable and interpretable, we group algorithms by cost function (i.e.SSD or Project-Out), and optimization method (i.e.Gauss-Newton, Newton or Wiberg).Results for this experiment are reported in Figures 5, 6, 7, 8, 9 and 10.These figures have all the same structure and are composed of four figures and a table.Figures 5a, 6a, 7a, 8a, 9a and 10a report the Cumulative Error Distribution (CED), i.e the proportion of images vs normalized point-to-point error for each of the algorithms' groups.Tables 5e, 6e, 7e, 8e, 9e, and 10e summarize and complete the information on the previous CEDs by stating the proportion of images fitted with a normalized point-to-point error smaller than 0.02, 0.03 and 0.04; and by stating the mean, std and median of the final normalized point-to-point error.The aim of the previous figures and tables is to help us compare the final fitting accuracy obtained by each algorithm.On the other hand, Figures 5b, 6b, 7b, 8b, 9b and 10b report the mean normalized point-to-point error at each iteration while Figures 5c, 5d, 6c, 6d, 7c, 7d, 8c, 8d, 9c, 9d and 10c, 10d report the mean normalized cost at each iteration 28 .The aim of these figures is to help us compare the convergence properties of every algorithm.

SSD Gauss-Newton algorithms
Results for SSD Gauss-Newton algorithms are reported in Figure 5.We can observe that Inverse, Asymmetric and Bidirectional algorithms obtain a similar performance and significantly outperform Forward algorithms in terms of fitting accuracy, Figure 5a and Table 5e.In absolute terms, Bidirectional algorithms slightly outperform Inverse and Asymmetric algorithms.On the other hand, the difference in performance between the Simultaneous Schur and Alternated optimizations strategies are minimal for all algorithms and they were found to have no statistical significance.
Looking at Figures 5b, 5c and 5d there seems to be a clear (and obviously expected) correlation between the normalized point-to-point error and the normalized value of the cost function at each iteration.In terms of convergence, it can be seen that Forward algorithms converge slower than Inverse, Asymmetric and Bidirectional.Bidirectional algorithms converge slightly faster than Inverse algorithms and these slightly faster than Asymmetric algorithms.In this case, the Simultaneous Schur optimization strategy seems to converge slightly faster than the Alternated one for all SSD Gauss-Newton algorithms.

SSD Newton algorithms
Results for SSD Newton algorithms are reported on Figure 6.In this case, we can observe that the fitting performance of all algorithms decreases with respect to their Gauss-Newton counterparts Figure 6a and Table 6e.This is most noticeable in the case of Forward algorithms for which there is ∼ 20% drop in the proportion of images fitted below 0.02, 0.03 and 0.04 with respect to its Gauss-Newton equivalents.For these algorithms there is also a significant increase in the mean and median of the normalized point-to-point error.Asymmetric Newton algorithms also perform considerably worse, between 5% and 10%, than their Gauss-Newton versions.The drop in performance is reduced for Inverse and Bidirectional Newton algorithms for which accuracy is only reduced by around 3% with respect their Gauss-Newton equivalent.
Within Newton algorithms, there are clear differences in terms of speed of convergence 6b, 6c and 6d.Bidirectional algorithms are the fastest to converge followed by Inverse and Asymmetric algorithms, in this order, and lastly Forward algorithms.In this case, the Simultaneous Schur optimization strategy seems to converge again slightly faster than the Alternated one for all algorithms but Bidirectional algorithms, for which the Alternated strategy converges slightly faster.Overall, SSD Newton algorithms converge slower than SSD Gauss-Newton algorithms.

SSD Wiberg algorithms
Results for SSD Wiberg algorithms are reported on Figure 7. Figure 7a and Table 7e and Figures 7b, 7c and  7d show that these results are (as one would expect) virtually equivalent to those obtained by their Gauss-Newton counterparts.

Project-Out Gauss-Newton algorithms
Results for Project-Out Gauss-Newton algorithms are reported on Figure 8.We can observe that, there is significant drop in terms of fitting accuracy for Inverse and Bidirectional algorithms with respect to their SSD versions, 8a and Table 8e.As expected, the Forward algorithm achieves virtually the same results as its SSD counterpart.The Asymmetric algorithm obtains similar accuracy to that of the best performing SSD algorithms.
Looking at Figures 8b, 8c and 8d we can see that Inverse and Bidirectional algorithms converge slightly faster than the Asymmetric algorithm.However, the Asymmetric algorithm ends up descending to a significant lower value of the mean normalized cost which also translates to a lower value for the final mean normalized point-to-point error.Similar to SSD algorithms, the Forward algorithm is the worst convergent algorithm.
Finally notice that, in this case, there is virtually no difference, in terms of both final fitting accuracy and speed of convergence, between the Simultaneous Schur and Alternated optimizations strategies used by the Bidirectional algorithm.

Project-Out Newton algorithms
Results for Project-Out Newton algorithms are reported on Figure 9.It can be clearly seen that Project-Out Newton algorithms perform much worse than their Gauss-Newton and SSD counterparts.The final fitting accuracy obtained by these algorithms is very poor compared to the one obtained by the best SSD and Project-Out Gauss-Newton algorithms, Figures 9a and Table  9e.In fact, by looking at Figures 9b, 9c and 9d only the Forward and Asymmetric algorithms seem to be stable at the second level of the Gaussian pyramid with Inverse and Bidirectional algorithms completely diverging for some of the images as shown by the large mean and std of their final normalized point-to-point errors.

Project-Out Wiberg algorithms
Results for the Project-Out Bidirectional Wiberg algorithm are reported on Figure 9.As expected, the results are virtually identical to those of the obtained by Project-Out Bidirectional Gauss-Newton algorithms.

Weighted Bayesian project-out
In this experiment, we quantify the importance of each of the two terms in our Bayesian project-out cost function, Equation 22.To this end, we introduce the parameters, ρ ∈ [0, 1] and γ = 1 − ρ, to weight up the relative contribution of both terms: Setting ρ = 0, γ = 1 reduces the previous cost function to the original project-out loss proposed in [29]; completely disregarding the contribution of the prior distribution over the appearance parameters i.e the Mahalanobis distance within the appearance subspace.On the contrary, setting ρ = 1, γ = 0 reduces the cost function to the first term; completely disregarding the contribution of the project-out term i.e. the distance to the appearance subspace.Finally setting ρ = γ = 0.5 leads to the standard Bayesian project-out cost function proposed in Section 3.1.2.In order to assess the impact that each term has on the fitting accuracy obtained by the previous Project-Out algorithm we repeat the experimental set up of the first experiment and test all Project-Out Gauss-Newton algorithms for different values of the parameters ρ = 1 − γ.Notice that, in this case, we only report the performance of Gauss-Newton algorithms because they were shown to vastly outperform Newton algorithms and to be virtually equivalent to Wiberg algorithms in the first experiment.
Results for this experiment are reported by Figure 11.We can see that, regardless of the type of composition, a weighted combination of the two previous terms always leads to a smaller mean normalized point-topoint error compared to either term on its own.Note that the final fitting accuracy obtained with the standard Bayesian project-out cost function is substantially better than the one obtained by the original projectout loss (this is specially noticeable for the Inverse and Bidirectional algorithms); fully justifying the inclusion of the first term, i.e the Mahalanobis distance within the appearance subspace, into the cost function.Finally, in this particular experiment, the final fitting accuracy of all algorithms is maximized by setting ρ = 0.1, γ = 0.9, further highlighting the importance of the first term in the Bayesian formulation.

Optimal asymmetric composition
This experiment quantifies the effect that varying the value of the parameters α ∈ [0, 1] and β = 1 − α in Equation 30 has in the fitting accuracy obtained by the Asymmetric algorithms.Note that for α = 1, β = 0 and α = 0, β = 1 these algorithms reduce to their Forward and Inverse versions respectively.Recall that, in previous experiments, we used the Symmetric case α = β = 0.5 to generate the results reported for Asymmetric algorithms.Again, we only report performance for Gauss-Newton algorithms.
We again repeat the experimental set up described in the first experiments and report the fitting accu-racy obtained by the Project Out and SSD Asymmetric Gauss-Newton algorithms for different values of the parameters α = 1−β.Results are shown in Figure 13.For the BPO Asymmetric algorithm, the best results are obtain by setting α = 0.4, β = 0.6, Figures 13a (top)  and 13b.These results slightly outperform those obtain by the default Symmetric algorithm and this particular configuration of the BPO Asymmetric algorithm is the best performing one on the LFPW test dataset.For the SSD Asymmetric Gauss-Newton algorithm the best results are obtained by setting α = 0.2, β = 0.8, Figures 13a (bottom) and 13c.In this case, the boost in performance with respect to the default Symmetric algorithm is significant and, with this particular configuration, the SSD Asymmetric Gauss-Newton algorithm is the best performing SSD algorithm on the LFPW test dataset, outperforming Inverse and Bidirectional algorithms.

Sampling and Number of Iterations
In this experiment, we explore two different strategies to reduce the running time of the previous CGD algorithms.
The first one consists of optimizing the SSD and Project-Out cost functions using only a subset of all pixels in the reference frame.In AAMs the total number of pixels on the reference frame, F , is typically several orders of magnitude bigger than the number of shape, n, and appearance, m, components i.e.F >> m >> n.Therefore, a significant reduction in the complexity (and running time) of CGD algorithms can be obtained by decreasing the number of pixels that are used to optimize the previous cost functions.To this end, we compare the accuracy obtained by using 100%, 50%, 25% and 12% of the total number of pixels on the reference frame.Note that, pixels are (approximately) evenly sampled across the reference frame in all cases, Figure 4.
The second strategy consists of simply reducing the number of iterations that each algorithm is run.Based on the figures used to assess the convergence properties of CGD algorithms in previous experiments, we compare the accuracy obtained by running the algorithms for 40 (24 + 16) and 20 (12 + 8) iterations.
Note that, in order to further highlight the advantages and disadvantages of using the previous strategies we report the fitting accuracy obtained by initializing the algorithms using different amounts of uniform noise.
Once more we repeat the experimental set up of the first experiment and report the fitting accuracy obtained by the Project Out and SSD Asymmetric Gauss-Newton algorithms.Results for this experiment are shown in Figure 12.It can be seen that reducing the number of pixels up to 25% while maintaining the original number of iterations to 40 (24 + 16) has little impact on the fitting accuracy achieved by both algorithms while reducing them to 12% has a clear negative impact, Figures 12a and 12b.Also, performance seems to be consistent along the amount of noise.In terms of run time, Table 12c, reducing the number of pixels to 50%, 25% and 12% offers speed ups of ∼ 2.0x, ∼ 2.9x and ∼ 3.7x for the BPO algorithm and of ∼ 1.8x, ∼ 2.6x and ∼ 2.8x for the SSD algorithm respectively.
On the other hand, reducing the number of iterations from 40 (24 + 16) to 20 (12 + 8) has no negative impact in performance for levels of noise smaller than 2% but has a noticeable negative impact for levels of noise bigger than 5%.Notice that remarkable speed ups, Table 12f, can be obtain for both algorithms by combining the previous two strategies at the expenses of small but noticeable decreases in fitting accuracy.

Comparison on Helen and AFW
In order to facilitate comparisons with recent prior work on AAMs [49, 3, 21] and with other state-of-the-art approaches in face alignment [55,4], in this experiment, we report the fitting accuracy of the SSD and Project-Out Asymmetric Gauss-Newton algorithms on the widely used test set of the Helen database and on the entire AFW database.Furthermore we compare the performance of the previous two algorithms with the one obtained by the recently proposed Gauss-Newton Deformable Part Models (GN-DPMs) proposed by Tzimiropoulos and Pantic in [50]; which was shown to achieve state-of-the-art results in the problem of face alignment in-the-wild.
For both our algorithms, we report two different types of results: i) sampling rate of 25% and 20 (12 + 8) iterations; and ii) sampling rate of 50% and 40 (24+16) iterations, .For GN-DPMs we use the authors public implementation to generate the results.In this case, we report, again, two different types of results by letting the algorithm run for 20 and 40 iterations.
Result for this experiment are shown in Figure 14.Looking at Figure 14a we can see that both, SSD and Project-Out Asymmetric Gauss-Newton algorithms, obtain similar fitting accuracy on the Helen test dataset.Note that, in all cases, their accuracy is comparable to the one achieved by GN-DPMs for normalized pointto-point errors < 0.2 and significantly better for < 0.3, < 0.4.As expected, the best results for both our algorithms are obtained using 50% of the total amount of pixels and 40 (24 + 16) iterations.However, the results obtained by using only 25% of the total amount of pixels and 20 (12 + 8) iterations are comparable to the previous ones; specially for the Project-Out Asymmetric Gauss-Newton.In general, these results are consistent with the ones obtained on the LFPW test dataset, Experiments 5.1 and 5.3.
On the other hand, the performance of both algorithms drops significantly on the AFW database, Figure 14b .In this case, GN-DPMs achieves slightly better results than the SSD and Project-Out Asymmetric Gauss-Newton algorithms for normalized point-topoint errors < 0.2 and slightly worst for < 0.3, < 0.4.Again, both our algorithms obtain better results by using 50% sampling rate and 40 (24 + 16) iterations and the difference in accuracy with respect to the versions using 25% sampling rate and 20 (12 + 8) iterations slightly widens when compared to the results obtained on the Helen test dataset.This drop in performance is consistent with other recent works on AAMs [50, 30,3,31] and it is attributed to large difference in terms of shape and appearance statistics between the images of the AFW dataset and the ones of the training sets of the LFPW and Helen datasets where the AAM model was trained on.
Exemplar results for this experiment are shown in Figures 16 and 17.

Comparison on MIT StreetScene
In this final experiment, we present results for a different type of object: cars.To this end, we use the first view of the MIT StreetScene 27 dataset containing a wide variety of frontal car images obtained in the wild.We use 10-fold cross-validation on the ∼ 500 images of the previous dataset to train and test our algorithms.We report results for the two versions of the SSD Asymmetric Gauss-Newton and the Project-Out Asymmetric Gauss-Newton algorithms used in Experiment 5.5.
Result for this experiment are shown in Figure 15.We can observe that all algorithms obtain similar performance and that they vastly improve upon the original initialization.
Exemplar results for this experiment are shown in Figure 18.

Analysis
Given the results reported by the previous six experiments we conclude that: 1. Overall, Gauss-Newton and Wiberg algorithms vastly outperform Newton algorithms for fitting AAMs.Experiment 5.1 clearly shows that the former algorithms provide significantly higher levels of fitting accuracy at considerably lower computational complexities and run times.These findings are consistent with existent literature in the related field of parametric image alignment [29] and also, to certain extend, with prior work on Newton algorithms for AAM fitting [21].We attribute the bad performance of Newton algorithms to the difficulty of accurately computing a (noiseless) estimate of the full Hessian matrix using finite differences.2. Gauss-Newton and Wiberg algorithms are virtually equivalent in performance.The results in Experiment 5.1 show that the difference in accuracy between both types of algorithms is minimal and the small differences in their respective solutions are, in practice, insignificant.3. Our Bayesian project-out formulation leads to significant improvements in fitting accuracy without adding extra computational cost.Experiment 5.2 shows that a weighted combination of the two terms forming Bayesian project-out loss always outperforms the classic project out formulation.4. The Asymmetric composition proposed in this work leads to CGD algorithms that are more accurate and that converge faster.In particular, the SSD and Project-Out Asymmetric Gauss-Newton algorithms are shown to achieve significantly better performance than their Forward and Inverse counterparts in Experiments 5.1 and 5.3. 5. Finally, a significant reduction in the computational complexity and runtime of CDG algorithms can be obtained by limiting the number of pixels considered during optimization of the loss function and by adjusting the number of iterations that the algorithms are run for, Experiment 5.4.

Conclusion
In this paper we have thoroughly studied the problem of fitting AAMs using CGD algorithms.We have presented a unified and complete framework for these algorithms and classified them with respect to three of their main characteristics: i) cost function; ii) type of composition; and iii) optimization method.Furthermore, we have extended the previous framework by: -Proposing a novel Bayesian cost function for fitting AAMs that can be interpreted as a more general formulation of the well-known project-out loss.We have assumed a probabilistic model for appearance generation with both Gaussian noise and a Gaussian prior over a latent appearance space.Marginalizing out the latent appearance space, we have derived a novel cost function that only depends on shape parameters and that can be interpreted as a valid and more general probabilistic formulation of the wellknown project-out cost function [29].In the experiments, we have showed that our Bayesian formulation considerably outperforms the original projectout cost function.-Proposing asymmetric and bidirectional compositions for CGD algorithms.We have shown the connection between Gauss-Newton Asymmetric algorithms and ESM algorithms and experimentally proved that these two novel types of composition lead to better convergent and more robust CGD algorithm for fitting AAMs.-Providing new valuable insights into existent CGD algorithms by reinterpreting them as direct applications of the Schur complement and the Wiberg method.
Finally, in terms of future work, we plan to: -Adapt existent Supervised Descent (SD) algorithms for face alignment [55,48] to AAMs and investigate their relationship with the CGD algorithms studied in this paper.-Investigate if our Bayesian cost function and the proposed asymmetric and bidirectional compositions can also be successfully applied to similar generative parametric models, such as the Gauss-Newton Parts-Based Deformable Model (GN-DPM) proposed in [50].A Terms in SSD Newton Hessians In this section we define the individual terms of the Hessian matrices used by the SSD Asymmetric and Bidirectional Newton optimization algorithms derived in Section 3.3.2.

A.1 Asymmetric
The individual terms forming the Hessian matrix of the SSD Asymmetric Newton algorithm defined by Equation 74 are defined as follows: where we have defined ∂

A.2 Bidirectional
The individual terms forming the Hessian matrix of the SSD Bidirectional Newton algorithm defined by Equation 77 are defined as follows:

B Iterative solutions of all algorithms
In this section we report the iterative solutions of all CGD algorithms studied in this paper.In order to keep the information structured algorithms are grouped by their cost function.Consequently, iterative solutions for all SSD and Project-Out algorithms are stated in Tables 1 and 2.
Table 2: Iterative solutions of all Project-Out algorithms studied in this paper.
before, then solves approximately by performing a first order Taylor expansion around ∆q: ∆q * = arg min ∆q rb (∆q) + ∂r b ∂∆q ∆q 2

Fig. 3 :
Fig. 3: Exemplar initializations obtained by varying the percentage of uniform noise added to the similarity parameters.Note that, increasing the percentage of noise produces more challenging initialization.

Fig. 4 :
Fig. 4: Subset of pixels on the reference frame used to optimize the SSD and Project-Out cost functions for different sampling rates.
(a) CED on the LFPW test dataset for all SSD Gauss-Newton algorithms initialized with 5% uniform noise.(b) Mean normalized point-to-point error vs number of iterations on the LFPW test dataset for all SSD Gauss-Newton algorithms initialized with 5% uniform noise.(c) Mean normalized cost vs number of first scale iterations on the LFPW test dataset for all SSD Gauss-Newton algorithms initialized with 5% uniform noise.(d) Mean normalized cost vs number of second scale iterations on the LFPW test dataset for all SSD Gauss-Newton algorithms initialized with 5% uniform noise.

Fig. 5 :
Fig. 5: Results showing the fitting accuracy and convergence properties of the SSD Gauss-Newton algorithms on the LFPW test dataset initialized with 5% uniform noise.

( a )
Cumulative error distribution on the LFPW test dataset for all SSD Newton algorithms initialized with 5% uniform noise.(b) Mean normalized point-to-point error vs number of iterations on the LFPW test dataset for all SSD Newton algorithms initialized with 5% uniform noise.(c) Mean normalized cost vs number of first scale iterations on the LFPW test dataset for all SSD Newton algorithms initialized with 5% uniform noise.(d) Mean normalized cost vs number of second scale iterations on the LFPW test dataset for all SSD Newton algorithms initialized with 5% uniform noise.

Fig. 6 :
Fig. 6: Results showing the fitting accuracy and convergence properties of the SSD Newton algorithms on the LFPW test dataset initialized with 5% uniform noise.
(a) CED on the LFPW test dataset for all SSD Wiberg algorithms initialized with 5% uniform noise.(b) Mean normalized point-to-point error vs number of iterations on the LFPW test dataset for all SSD Wiberg algorithms initialized with 5% uniform noise.(c) Mean normalized cost vs number of first scale iterations on the LFPW test dataset for all SSD Wiberg algorithms initialized with 5% uniform noise.(d) Mean normalized cost vs number of second scale iterations on the LFPW test dataset for all SSD Wiberg algorithms initialized with 5% uniform noise.

Fig. 7 :
Fig. 7: Results showing the fitting accuracy and convergence properties of the SSD Wiberg algorithms on the LFPW test dataset.
(a) CED graph on the LFPW test dataset for all Project-Out Gauss-Newton algorithms initialized with 5% uniform noise.(b) Mean normalized point-to-point error vs number of iterations on the LFPW test dataset for all Project-Out Gauss-Newton algorithms initialized with 5% uniform noise.(c) Mean normalized cost vs number of first scale iterations on the LFPW test dataset for all Project-Out Gauss-Newton algorithms initialized with 5% uniform noise.(d) Mean normalized cost vs number of second scale iterations on the LFPW test dataset for all Project-Out Gauss-Newton algorithms initialized with 5% uniform noise.

Fig. 8 :
Fig. 8: Results showing the fitting accuracy and convergence properties of the Project-Out Gauss-Newton algorithms on the LFPW test dataset.
(a) CED graph on the LFPW test dataset for all Project-Out Newton algorithms initialized with 5% uniform noise.(b) Mean normalized point-to-point error vs number of iterations on the LFPW test dataset for all Project-Out Newton algorithms initialized with 5% uniform noise.(c) Mean normalized cost vs number of first scale iterations on the LFPW test dataset for all Project-Out Newton algorithms initialized with 5% uniform noise.(d) Mean normalized cost vs number of second scale iterations on the LFPW test dataset for all Project-Out Newton algorithms initialized with 5% uniform noise.

Fig. 9 :
Fig. 9: Results showing the fitting accuracy and convergence properties of the Project-Out Newton algorithms on the LFPW test dataset.
(a) Cumulative Error Distribution graph on the LFPW test dataset for all Project-Out Wiberg algorithms initialized with 5% uniform noise.(b)

Fig. 10 :
Fig. 10: Results showing the fitting accuracy and convergence properties of the Project-Out Wiberg algorithms on the LFPW test dataset.
(a) Proportion of images with normalized point-to-point errors smaller than 0.02, 0.03 and 0.04 for the Project-Out and SSD Asymmetric Gauss-Newton algorithms for different values of ρ = 1 − γ and initialized with 5% noise.Colors encode overall fitting accuracy, from highest to lowest: red, orange, yellow, green, blue and purple.(b) CED on the LFPW test dataset for Project-Out Forward Gauss-Newton algorithms for different values of ρ = 1 − γ and initialized with 5% noise.(c) CED on the LFPW test dataset for Project-Out Inverse Gauss-Newton algorithms for different values of ρ = 1 − γ and initialized with 5% noise.(d) CED on the LFPW test dataset for Project-Out Asymmetric Gauss-Newton algorithms for different values of ρ = 1 − γ and initialized with 5% noise.(e) CED on the LFPW test dataset for Project-Out Bidirectional Gauss-Newton algorithms for different values of ρ = 1 − γ and initialized with 5% noise.

Fig. 11 :
Fig. 11: Results quantifying the effect of varying the value of the parameters ρ = 1−γ in Project-Out Gauss-Newton algorithms.

( a )
Proportion of images with normalized point-to-point errors smaller than 0.02, 0.03 and 0.04 for the SSD Asymmetric Gauss-Newton algorithm using different sampling rates, 40 (24 + 16) iterations, and initialized with different amounts of noise.Colors encode overall fitting accuracy, from highest to lowest: red, orange, yellow, green, blue and purple.(b) Proportion of images with normalized point-to-point errors smaller than 0.02, 0.03 and 0.04 for the Project-Out Asymmetric Gauss-Newton algorithm using different sampling rates, 40 (24 + 16) iterations, and initialized with different amounts of noise.Colors encode overall fitting accuracy, from highest to lowest: red, orange, yellow, green, blue and purple.Sch ∼ 1680 ms ∼ 930 ms ∼ 650 ms ∼ 590 ms PO Asy GN ∼ 1400 ms ∼ 680 ms ∼ 480 ms ∼ 380 ms (c) Table showing run time of each algorithm for different amounts of sampling and 40 (24 + 16) iterations.(d) Proportion of images with normalized point-to-point errors smaller than 0.02, 0.03 and 0.04 for the Project-Out Asymmetric Gauss-Newton algorithm using different sampling rates, 20 (12+8) iterations, and initialized with different amounts of noise.Colors encode overall fitting accuracy, from highest to lowest: red, orange, yellow, green, blue and purple.(e) Proportion of images with normalized point-to-point errors smaller than 0.02, 0.03 and 0.04 for the SSD Asymmetric Gauss-Newton algorithm using different sampling rates, 20 (12 + 8) iterations, and initialized with different amounts of noise.Colors encode overall fitting accuracy, from highest to lowest: red, orange, yellow, green, blue and purple.Sch ∼ 892 ms ∼ 519 ms ∼ 369 ms ∼ 331 ms PO Asy GN ∼ 707 ms ∼ 365 ms ∼ 235 ms ∼ 211 ms (f)Table showing run time of each algorithm for different amounts of sampling and 20 (12 + 8) iterations.

Fig. 12 :
Fig.12: Results assessing the effectiveness of sampling for the best performing Project-Out and SSD algorithms on the LFPW database.

( a )
Proportion of images with normalized point-to-point errors smaller than 0.02, 0.03 and 0.04 for the Project-Out and SSD Asymmetric Gauss-Newton algorithms for different values of α = 1 − β and initialized with 5% noise.Colors encode overall fitting accuracy, from highest to lowest: red, orange, yellow, green, blue and purple.(b) CED on the LFPW test dataset for Project-Out Asymmetric Gauss-Newton algorithm for different values of α = 1 − β and initialized with 5% noise.(c) CED on the LFPW test dataset for the the SSD Asymmetric Gauss-Newton algorithm for different values of α = 1 − β and initialized with 5% noise.

Fig. 13 :
Fig. 13: Results quantifying the effect of varying the value of the parameters α = 1 − β in Asymmetric algorithms.
(a) CED on the Helen test dataset for the Project-Out and SSD Asymmetric Gauss-Newton algorithms initialized with 5% noise.(b) CED on the AFW database for the Project-Out and SSD Asymmetric Gauss-Newton algorithm initialized with 5% noise.

Fig. 14 :
Fig. 14: Results showing the fitting accuracy of the SSD and Project-Out Asymmetric Gauss-Newton algorithms on the Helen and AFW databases.

Fig. 15 :
Fig. 15: CED on the first view of the MIT StreetScene test dataset for the Project-Out and SSD Asymmetric Gauss-Newton algorithms initialized with 5% noise.

( a )
Exemplar results from the Helen test dataset obtained by the Project-Out Asymmetric Gauss-Newton Schur algorithm.(b) Exemplar results from the Helen test dataset obtained by the SSD Asymmetric Gauss-Newton Schur algorithm.

( a )
Exemplar results from the Helen test dataset obtained by the Project-Out Asymmetric Gauss-Newton Schur algorithm.(b) Exemplar results from the AFW dataset obtained by the SSD Asymmetric Gauss-Newton Schur algorithm.
r b ∂∆c = −A T ∂r b ∂∆c = A T A I (129) ∂ 2 D b ∂∆c∂∆p = ∂ − A T r b ∂∆p = −A T ∂r b ∂∆p = −A T J i (130) ∂ 2 D b ∂∆c∂∆q = ∂ − A T r b ∂∆q = ∂ − A T ∂∆q r b − A T ∂r b ∂∆q = −J T A r b + A T J a

1 J 1 J 1 J− 1
Sch [2, 49] ∆p = − Ĥ−1 i J T i Ār ∆c = A (r + J i ∆p) Ĥi = J T i ĀJ i SSD For GN Alt ∆p = −H −1 i J T i (r − A∆c) ∆c = A (r + J i ∆p) H i = J T i J i SSD For N Sch ∆p = − ĤN i −T i Ār ∆c = A (r + J i ∆p) ĤN i = ∂W ∆p T ∇ 2 i ∂W ∆p r + Ĥi SSD For N Alt ∆p = − H N i −1 J T i Ā (r − A∆c) ∆c = A (r + J i ∆p) H N i = ∂W ∆p T ∇ 2 i ∂W ∆p r + H i SSD For W ∆p = − Ĥ−1 i J T i Ār ∆c = Ar SSD Inv GN Sch [38, 49] ∆p = Ĥ−1 a J T a Ār ∆c = A (r − J a ∆p) Ĥa = J T a ĀJ a SSD Inv GN Alt [51, 3] ∆p = H −1 a J T a (r − A∆c) ∆c = A (r − J a ∆p) H a = J T a J a SSD Inv N Sch ∆p = ĤN a −T a Ār ∆c = A (r − J a ∆p) ĤN a = ∂W ∆p T ∇ 2 a ∂W ∆p r + Ĥa SSD Inv N Alt ∆p = H N a −1 J T a Ā (r − A∆c) ∆c = A (r − J a ∆p) H N a = ∂W ∆p T ∇ 2 i ∂W ∆p r + H a SSD Inv W ∆p = Ĥ−1 a J T a Ār ∆c = Ar SSD Asy GN Sch ∆p = − Ĥ−1 t J T t Ār ∆c = A (r + J t ∆p) Ĥt = J T t ĀJ t SSD Asy GN Alt ∆p = −H −1 t J T t (r − A∆c) ∆c = A (r + J t ∆p) H t = J T t J t SSD Asy N Sch ∆p = − ĤN t −T t Ār ∆c = A (r + J t ∆p) ĤN t = ∂W ∆p T ∇ 2 t ∂W ∆p r + Ĥt SSD Asy N Alt ∆p = − H N t −1 J T t Ā (r − A∆c) ∆c = A (r + J t ∆p) H N t = ∂W ∆p T ∇ 2 t ∂W ∆p r + H t SSD Asy W ∆p = − Ĥ−1 t J T t Ār ∆c = Ar SSD Bid GN Sch ∆p = − Ĥ−1 i J T i Ār 1 ∆q = Ȟ−1 a J T a Pr ∆c = Ar 2 r 1 = (r − J a ∆q) Ȟa = J T a PJ a r 2 = (r + J i ∆p − J a ∆q) P = Ā − ĀJ i Ĥ−1 i J T i Ā SSD Bid GN Alt ∆p = −H −1 i J T i r 3 ∆q = H −1 a J T a r 4 ∆c = Ar 2 r 3 = (r − A∆c − J a ∆q) r 4 = (r − A∆c + J i ∆p) J T a P N r ∆c = Ar 2 ȞN a = ∂W ∆p T ∇ 2 t ∂W ∆p r + Ȟa P N = Ā − ĀJ i ĤN i arg min Table showing the proportion of images fitted with a normalized point-to-point error below 0.02, 0.03 and 0.04 together with the normalized point-to-point error mean, std and median for all SSD Wiberg algorithms initialized with 5% uniform noise.
Mean normalized point-to-point error vs number of iterations graph on the LFPW test dataset for all Project-Out Wiberg algorithms initialized with 5% uniform noise.Table showing the proportion of images fitted with a normalized point-to-point error below 0.02, 0.03 and 0.04 together with the normalized point-to-point error mean, std and median for all Project-Out Wiberg algorithms initialized with 5% uniform noise.
(c) Mean normalized cost vs number of first scale iterations graph on the LFPW test dataset for all Project-Out Wiberg algorithms initialized with 5% uniform noise.(d) Mean normalized cost vs number of second scale iterations graph on the LFPW test dataset for all Project-Out Wiberg algorithms initialized with 5% uniform noise.Algorithm < 0.02 < 0.03 < 0.04 Mean

Table 1 :
Iterative solutions of all SSD algorithms studied in this paper.
i J T i Ā (r − Jā∆q)