Title: NON-LINEAR DIMENSIONALITY REGULARIZER FOR SOLVING INVERSE PROBLEMS

Abstract: Consider an ill-posed inverse problem of estimating causal factors from observations, one of which is known to lie near some (unknown) low-dimensional, nonlinear manifold expressed by a predefined Mercer-kernel. Solving this problem requires simultaneous estimation of these factors and learning the low-dimensional representation for them. In this work, we introduce a novel non-linear dimensionality regularization technique for solving such problems without pre-training. We re-formulate Kernel-PCA as an energy minimization problem in which low dimensionality constraints are introduced as regularization terms in the energy. To the best of our knowledge, ours is the first attempt to create a dimensionality regularizer in the KPCA framework. Our approach relies on robustly penalizing the rank of the recovered factors directly in the implicit feature space to create their low-dimensional approximations in closed form. Our approach performs robust KPCA in the presence of missing data and noise. We demonstrate state-of-the-art results on predicting missing entries in the standard oil flow dataset. Additionally, we evaluate our method on the challenging problem of Non-Rigid Structure from Motion and our approach delivers promising results on CMU mocap dataset despite the presence of significant occlusions and noise.

1. INTRODUCTION:

Dimensionality reduction techniques are widely used in data modeling, visualization and unsupervised learning. Principal component analysis (PCAJolliffe (2002) ), Kernel PCA (KPCASchölkopf et al. (1998) ) and Latent Variable Models (LVMsLawrence (2005) ) are some of the well known techniques used to create low dimensional representations of the given data while preserving its significant information.

One key deployment of low-dimensional modeling occurs in solving ill-posed inference problems. Assuming the valid solutions to the problem lie near a low-dimensional manifold (i.e. can be parametrized with a reduced set of variables) allows for a tractable inference for otherwise underconstrained problems. After the seminal work of Candès & Recht (2009) ; Recht et al. (2010) on guaranteed rank minimization of the matrix via trace norm heuristics Fazel (2002) , many ill-posed computer vision problems have been tackled by using the trace norm -a convex surrogate of the rank function -as a regularization term in an energy minimization frameworkCandès & Recht (2009) ; Zhou et al. (2014) . The flexible and easy integration of low-rank priors is one of key factors for versatility and success of many algorithms. For example, pre-trained active appearance models Cootes et al. (2001) or 3D morphable models Blanz & Vetter (1999) are converted to robust feature tracking Poling et al. (2014) , dense registration Garg et al. (2013b) and vivid reconstructions of natural videos Garg et al. (2013a) with no a priori knowledge of the scene. Various bilinear factorization problems like background modeling, structure from motion or photometric stereo are also addressed with a variational formulation of the trace norm regularization Cabral et al. (2013) .

On the other hand, although many non-linear dimensionality reduction techniques -in particular KPCA -have been shown to outperform their linear counterparts for many data modeling tasks, they are seldom used to solve inverse problems without using a training phase. A general (discriminative) framework for using non-linear dimensionality reduction is: (i) learn a low-dimensional representation for the data using training examples via the kernel trick (ii) project the test examples on the learned manifold and finally (iii) find a data point (pre-image) corresponding to each projection in the input space. This setup has two major disadvantages. Firstly, many problems of interest come with corrupted observations -noise, missing data and outliers -which violate the low-dimensional modeling assumption.Secondly, computing the pre-image of any point in the low dimensional feature subspace is non-trivial: the pre-image for many points in the low dimensional space might not even exist because the non linear feature mapping function used for mapping the data from input space to the feature space is non-surjective.

Previously, extensions to KPCA like Robust KPCA (RKPCANguyen & De la Torre ( 2009)) and probabilistic KPCA (PKPCASanguinetti & Lawrence (2006) ) with missing data have been proposed to address the first concern, while various additional regularizers have been used to estimate the pre-image robustly Bakir et al. (2004) ; Mika et al. (1998) ; Kwok & Tsang (2004) ; Abrahamsen & Hansen (2009) .

Generative models like LVMs Lawrence (2005) are often used for inference by searching the lowdimensional latent space for a location which maximizes the likelihood of the observations. Problems like segmentation, tracking and semantic 3D reconstruction Prisacariu & Reid (2011); Dame et al. (2013) greatly benefit from using LVM. However, the latent space is learned a priori with clean training data in all these approaches. Almost all non-linear dimensionality reduction techniques are non-trivial to generalize for solving ill-posed problems (See section 4.2) without a pre-training stage. Badly under-constrained problems require the low-dimensional constraints even for finding an initial solution, eliminating applicability of the standard "projection + pre-image estimation" paradigm. This hinders the utility of nonlinear dimensionality reduction and a suitable regularization technique to penalize the non-linear dimensionality is desirable. 

Sum and Substance::

A closer look at most non-linear dimensionality reduction techniques reveals that they rely upon a non-linear mapping function which maps the data from input space to a (usually) higher dimensional feature space. In this feature space the data is assumed to lie on a low-dimensional hyperplane -thus, linear low-rank prior is apt in the feature space. Armed with this simple observation, our aim is to focus on incorporating the advances made in linear dimensionality reduction techniques to their non-linear counterparts, while addressing the problems described above. Figure 1 In particular we propose a novel robust energy minimization algorithm which handles the implicitness of the feature space to directly penalize its rank by iteratively: (i) creating robust low-dimensional representation for the data given the kernel matrix in closed form and (ii) reconstructing the noise-free version of the data (pre-image of the features space projections) using the estimated low-dimensional representations in a unified framework. The proposed algorithm: (i) provides a novel closed form solution to robust KPCA; (ii) yields stateof-the-art results on missing data prediction for the well-known oil flow dataset; (iii) outperforms state-of-the-art linear dimensionality (rank) regularizers to solve NRSfM; and (iv) can be trivially generalized to incorporate other cost functions in an energy minimization framework to solve various ill-posed inference problems.

2. PROBLEM FORMULATION:

This paper focuses on solving a generic inverse problem of recovering causal factor

S = [s 1 , s 2 , • • • s N ] ∈ X × N from N observations W = [w 1 , w 2 , • • • w N ] ∈ Y × N such that f (W, S) = 0.

Here function f (observation,variable), is a generic loss function which aligns the observations W with the variable S (possibly via other causal factors. e.g. R or Z in Section 4.1 and 4.2).

If, f (W, S) = 0 is ill-conditioned (for example when Y X ), we want to recover matrix S under the assumption that the columns of it lie near a low-dimensional non-linear manifold. This can be done by solving a constrained optimization problem of the following form:

EQUATION

where

Φ(S) = [φ(s 1 ), φ(s 2 ), • • • , φ(s N )] ∈ H × N

is the non-linear mapping of matrix S from the input space X to the feature space H (also commonly referred as Reproducing Kernel Hilbert Space), via a non-linear mapping function φ : X → H associated with a Mercer kernel K such that K(S) i,j = φ(s i ) T φ(s j ).

In this paper we present a novel energy minimization framework to solve problems of the general form (1).

As our first contribution, we relax the problem (1) by using the trace norm of Φ(S) -the convex surrogate of rank function -as a penalization function. The trace norm M * =: i λ i (M ) of a matrix M is the sum of its eigenvalues λ i (M ) and was proposed as a tight convex relaxation1 of the rank(M ) and is used in many vision problems as a rank regularizer Fazel (2002) . Although the rank minimization via trace norm relaxation does not lead to a convex problem in presence of a non-linear kernel function, we show in 3.2 that it leads to a closed-form solution to denoising a kernel matrix via penalizing the rank of recovered data (S) directly in the feature space.

With these changes we can rewrite (1) as:

min S f (W, S) + τ Φ(S) * (2)

where τ is a regularization strength.2 

It is important to notice that although the rank of the kernel matrix K(S) is equal to the rank of Φ(S), K(S) * is merely Φ(S) 2 F . Thus, directly penalizing the sum of the singular values of K(S) will not encourage low-rank in the feature space. 3Although we have relaxed the non-convex rank function, (2) is in general difficult to minimize due to the implicitness of the feature space. Most widely used kernel functions like RBF do not have a explicit definition of the function φ. Moreover, the feature space for many kernels is high-(possibly infinite-) dimensional, leading to intractability. These issues are identified as the main barriers to robust KPCA and pre-image estimation Nguyen & De la Torre (2009) . Thus, we have to reformulate (2) by applying kernel trick where the cost function (2) can be expressed in terms of the kernel function alone.

The key insight here is that under the assumption that kernel matrix K(S) is positive semidefinite, we can factorize it as: K(S) = C T C. Although, this factorization is non-unique, it is trivial to show the following:

EQUATION

where λ i (.) is the function mapping the input matrix to its i th largest eigenvalue.

The row space of matrix C in (3) can be seen to span the eigenvectors associated with the kernel matrix K(S) -hence the principal components of the non-linear manifold we want to estimate.

Using (3), problem (2) can finally be written as:

min S,C f (W, S) + τ C * s.t. K(S) = C T C (4)

The above minimization can be solved with a soft relaxation of the manifold constraint by assuming that the columns of S lie near the non-linear manifold.

min S,C f (W, S) + ρ 2 K(S) -C T C 2 F + τ C * (5)

As ρ → ∞, the optimum of ( 5) approaches the optimum of (4) . A local optimum of ( 4) can be achieved using the penalty method of Nocedal & Wright (2006) by optimizing (5) while iteratively increasing ρ as explained in Section 3.

Before moving on, we would like to discuss some alternative interpretations of ( 5) and its relationship to previous work -in particular LVMs. Intuitively, we can also interpret (5) from the probabilistic viewpoint as commonly used in latent variable model based approaches to define kernel function Lawrence (2005) . For example a RBF kernel with additive Gaussian noise and inverse width γ can be defined as: K(S) i,j = e -γ si-sj 2 + , where ∼ N (0, σ). In other words, with a finite ρ, our model allows the data points to lie near a non-linear low-rank manifold instead of on it. Its worth noting here that like LVMs, our energy formulation also attempts to maximize the likelihood of regenerating the training data W , (by choosing f (W, S) to be a simple least squares cost) while doing dimensionality reduction.

Note that in closely related work Geiger et al. (2009) , continuous rank penalization (with a logarithmic prior) has also been used for robust probabilistic non-linear dimensionality reduction and model selection in LVM framework. However, unlike Geiger et al. (2009); Lawrence (2005) where the non-linearities are modeled in latent space (of predefined dimensionality), our approach directly penalizes the non-linear dimensionality of data in a KPCA framework and is applicable to solve inverse problems without pre-training.

3. OPTIMIZATION:

We approach the optimization of (5) by solving the following two sub-problems in alternation:

min S f (W, S) + ρ 2 K(S) -C T C 2 F (6) min C τ C * + ρ 2 K(S) -C T C 2 F (7)

Algorithm 1 outlines the approach and we give a detailed description and interpretations of both sub-problems ( 7) and ( 6) in next two sections of the paper.

3.1. PRE-IMAGE ESTIMATION TO SOLVE INVERSE PROBLEM.:

Subproblem ( 6) can be seen as a generalized pre-image estimation problem: we seek the factor s i , which is the pre-image of the projection of φ(s i ) onto the principle subspace of the RKHS stored in Algorithm 1: Inference with Proposed Regularizer.

Input: Initial estimate S 0 of S.

Output: Low-dimensional S and kernel representation C. Parameters: Initial ρ 0 and maximum ρmax penalty, with scale ρs.

-S = S 0 , ρ = ρ 0 ; while ρ ≤ ρmax do while not converged do -Fix S and estimate C via closed-form solution of (7) using Algorithm 2; -Fix C and minimize (6) to update S using LM algorithm; -ρ = ρρs ;

C T C, which best explains the observation w i . Here ( 6) is generally a non-convex problem, unless the Mercer-kernel is linear, and must therefore be solved using non-linear optimization techniques.

In this work, we use the Levenberg-Marquardt algorithm for optimizing (6).

Notice that (6) only computes the pre-image for the feature space projections of the data points with which the non-linear manifold (matrix C) is learned. An extension to our formulation is desirable if one wants to use the learned non-linear manifold for denoising test data in a classic pre-image estimation framework. Although a valuable direction to pursue, it is out of scope of the present paper.

3.2. ROBUST DIMENSIONALITY REDUCTION:

Algorithm 2: Robust Dimensionality Reduction.

Input: Current estimate of S.

Output: Low-dimensional representation C. Parameters: Current ρ and regularization strength τ .

-[U Λ U T ] = Singular Value Decomposition of K(S); // Λ is a diagonal matrix, storing N singular values λi of K(S). for i = 1 to N do -Find three solutions (lr : r ∈ {1, 2, 3}) of:

l 3 -lλi + τ 2ρ = 0 ; -set l4 = 0; -lr = max(lr, 0) ∀r ∈ {1, 2, 3, 4} ; -r = arg min r { ρ 2 λi -l 2 r 2 + τ lr } ; -λi = lr ; -C = ΛU T ; // Λ is diagonal matrix storing λi.

One can interpret sub-problem (7) as a robust form of KPCA where the kernel matrix has been corrupted with Gaussian noise and we want to generate its low-rank approximation. Although ( 7) is non-convex we can solve it in closed-form via singular value decomposition. This closed-form solution is outlined in Algorithm 2 and is based on the following theorem: Theorem 1. With S n A 0 let A = U ΣU T denote its singular value decomposition. Then

min L ρ 2 ||A -L T L|| 2 F + τ ||L|| * (8) = n i=1 ρ 2 (σ i -γ * 2 i ) 2 + τ γ * i . (9)

A minimizer L * of (8) is given by

L * = Γ * U T (10) with Γ * ∈ D n + , γ * i ∈ {α ∈ R + | p σi,τ /2ρ (α) = 0} {0}, where p a,b denotes the depressed cubic p a,b (x) = x 3 -ax + b. D n

+ is the set of n-by-n diagonal matrices with non-negative entries.

Theorem 1 shows that each eigenvalue of the minimizer C * of (7) can be obtained by solving a depressed cubic whose coefficients are determined by the corresponding eigenvalue of the kernel matrix and the regularization strength τ . The roots of each cubic, together with zero, comprise a set of candidates for the corresponding eigenvalue of C * . The best one from this set is obtained by choosing the value which minimizes (9) (see Algorithm 2).

As elaborated in Section 2, problem (7) can be seen as regularizing sum of square root (L 1/2 norm) of the eigenvalues of the matrix K(S). In a closely related work Zongben et al. ( 2012), authors advocate L 1/2 norm as a better approximation for the cardinality of a vector then the more commonly used L 1 norm. A closed form solution for L 1/2 regularization similar to our work was outlined in Zongben et al. ( 2012) and was shown to outperform the L 1 vector norm regularization for sparse coding. To that end, our Theorem 1 and the proposed closed form solution (Algo 2) for ( 7) can 

4. EXPERIMENTS:

In this section we demonstrate the utility of the proposed algorithm. The aims of our experiments are twofold: (i) to compare our dimensionality reduction technique favorably with KPCA and its robust variants; and (ii) to demonstrate that the proposed non-linear dimensionality regularizer consistently outperforms its linear counterpart (a.k.a. nuclear norm) in solving inverse problems.

4.1. MATRIX COMPLETION:

The nuclear norm has been introduced as a low rank prior originally for solving the matrix completion problem. Thus, it is natural to evaluate its non-linear extensions on the same task. Assuming W ∈ R m×n to be the input matrix and Z a binary matrix specifying the availability of the observations in W , Algorithm 1 can be used for recovering a complete matrix S with the following choice of f (W, Z, S):

f (W, Z, S) = Z • (W -S) 2 F (11) where • represents Hadamard product.

To demonstrate the robustness of our algorithm for matrix completion problem, we choose 100 training samples from the oil flow dataset described in section 3.2 and randomly remove the elements from the data with varying range of probabilities to test the performance of the proposed algorithm against various baselines. Following the experimental setup as specified in Sanguinetti & Lawrence (2006) , we repeat the experiments with 50 different samples of Z. We report the mean and standard deviation of the root mean square reconstruction error for our method with the choice of τ = 0.1, alongside five different methods in Table 1 . Our method significantly improves the performance of missing data completion compared to other robust extensions of KPCA Tipping & Bishop (1999); Sanguinetti & Lawrence (2006) ; Nguyen & De la Torre (2009), for every probability of missing data.

Although we restrict our experiments to least-squares cost functions, it is vital to restate here that our framework could trivially incorporate robust functions like the L 1 norm instead of the Frobenius norm -as a robust data term f (W, Z, S) -to generalize algorithms like Robust PCA Wright et al. (2009) to their non-linear counterparts. Non-rigid structure from motion under orthography is an ill-posed problem where the goal is to estimate the camera locations and 3D structure of a deformable objects from a collection of 2D images which are labeled with landmark correspondences Bregler et al. (2000) . Assuming s i (x j ) ∈ R 3 to be the 3D location of point x j on the deformable object in the i th image, its orthographic projection w i (x j ) ∈ R 2 can be written as w i (x) = R i s i (x j ), where R i ∈ R 2×3 is a orthographic projection matrix Bregler et al. (2000) . Notice that as the object deforms, even with given camera poses, reconstructing the sequence by least-squares reprojection error minimization is an ill-posed problem. In their seminal work, Bregler et al. (2000) proposed to solve this problem with an additional assumption that the reconstructed shapes lie on a low-dimensional linear subspace and can be parameterized as linear combinations of a relatively low number of basis shapes. NRSfM was then cast as the low-rank factorization problem of estimating these basis shapes and corresponding coefficients. 

4.2. KERNEL NON-RIGID STRUCTURE FROM MOTION:

Z i (x j ) w i (x j ) -R i s i (x j ) 2 F (12)

where S is the shape matrix whose columns are 3N dimensional vectors storing the 3D coordinates S i (x j ) of the shapes and Z i (x j ) is a binary variable indicating if projection of point x j is available in the image i.

Assuming the projection matrices to be fixed, this problem is convex and can be exactly solved with standard convex optimization methods. Additionally, if the 2D projections w i (x j ) are noise free, optimizing (12) with very small τ corresponds to selecting the the solution -out of the many solutions -with (almost) zero projection error, which has minimum trace norm Dai et al. (2014) . Thus henceforth, optimization of ( 12) is referred as the trace norm heuristics (TNH). We solve this problem with a first order primal-dual variant of the algorithm given in Garg et al. (2013a) , which can handle missing data. The algorithm is detailed and compared favorably with the state of the art NRSfM approaches (based on linear dimensionality regularization) Appendix C.

A simple kernel extension of the above optimization problem is:

min S,R τ Φ(S) * + F i=1 N j=1 Z i (x j ) w i (x j ) -R i s i (x j ) 2 F f (W,Z,R,S) (13)

where Φ(S) is the non-linear mapping of S to the feature space using an RBF kernel.

With fixed projection matrices R, (13) is of the general form (2), for which the local optima can be found using Algorithm 1. In our experiments we use ground truth camera projection matrices to compare our algorithm against TNH. The advantage of this setup is that with ground-truth rotation and no noise, we can avoid the model selection (finding optimal regularization strength τ ) by setting it low enough. We run the TNH with τ = 10 -7 and use this reconstruction as initialization for Algorithm 1. For the proposed method, we set τ = 10 -4 and use following RBF kernel width selection approach:

• Maximum distance criterion (d max ): we set the maximum distance in the feature space to be 3σ. Thus, the kernel matrix entry corresponding to the shape pairs obtained by TNH with maximum Euclidean distance becomes e -9/2 . • Median distance criterion (d med ): the kernel matrix entry corresponding to the median euclidean distance is set to 0.5.

Following the standard protocol in Dai et al. (2014) ; Akhter et al. (2009) , we quantify the reconstruction results with normalized mean 3D errors e 3D = 1 σF N i j e ij , where e ij is the euclidean distance of a reconstructed point j in frame i from the ground truth, σ is the mean of standard deviation for 3 coordinates for the ground truth 3D structures, and F, N are number of input images and number of points reconstructed.

Table 2 shows the results of the TNH and non-linear dimensionality regularization based methods using the experimental setup explained above, both without missing data and after randomly removing 50% of the image measurements. Our method consistently beats the TNH baseline and improves the mean reconstruction error by ∼ 40% with full data and by ∼ 25% when used with 50% missing data. Figure 2 shows qualitative comparison of the obtained 3D reconstruction using TNH and proposed non-lienar dimensionality regularization technique for some sample frames from various sequences. We refer readers to Appendix B for results with simultaneous reconstruction pose optimization.

5. CONCLUSION:

In this paper we have introduced a novel non-linear dimensionality regularizer which can be incorporated into an energy minimization framework, while solving an inverse problem. The proposed algorithm for penalizing the rank of the data in the feature space has been shown to be robust to noise and missing observations. We have picked NRSfM as an application to substantiate our arguments and have shown that despite missing data and model noise (such as erroneous camera poses) our algorithm significantly outperforms state-of-the-art linear counterparts.

Although our algorithm currently uses slow solvers such as the penalty method and is not directly scalable to very large problems like dense non-rigid reconstruction, we are actively considering alternatives to overcome these limitations. An extension to estimate pre-images with a problem-specific loss function is possible, and this will be useful for online inference with pre-learned lowdimensional manifolds.

Given the success of non-linear dimensionality reduction in modeling real data and overwhelming use of the linear dimensionality regularizers in solving real world problems, we expect that proposed non-linear dimensionality regularizer will be applicable to a wide variety of unsupervised inference problems: recommender systems; 3D reconstruction; denoising; shape prior based object segmentation; and tracking are all possible applications.

A PROOF OF THEOREM 3.1

Proof. We will prove theorem 1 by first establishing a lower bound for (8) and subsequently showing that this lower bound is obtained at L * given by (10). The rotational invariance of the entering norms allows us to write (8) as:

EQUATION

Expanding ( 14) we obtain

EQUATION

The inequality in (17) follows directly by applying Hölder's inequality to ( 16) and using the property that the column vectors wi are unitary.

Next, with L = ΓU T in (8) we have

EQUATION

Finally, since the subproblems in (18) are separable in γi, its minimizer must be KKT-points of the individual subproblems. As the constraints are simple non-negativity constraints, these KKT points are either (positive) stationary points of the objective functions or 0. It is simple to verify that the stationary points are given by the roots of the cubic function p σ i ,τ /2ρ . Hence it follows that there exists a

EQUATION

∀γi ≥ 0, which completes the proof.

A.1 VALIDATING THE CLOSED FORM SOLUTION:

Given the relaxations proposed in Section 2, our assertion that the novel trace regularization based non-linear dimensionality reduction is robust need to be substantiated. To that end, we evaluate our closed-form solution of Algorithm 2 on the standard oil flow dataset introduced in Bishop & James (1993).

This dataset comprises 1000 training and 1000 testing data samples, each of which is of 12 dimensions and categorized into one of three different classes. We add zero mean Gaussian noise with variance σ to the training data5 and recover the low-dimensional manifold for this noisy training data S σ with KPCA and contrast this with the results from Algorithm 2. An inverse width of the Gaussian kernel γ = 0.075 is used for all the experiments on the oil flow dataset.

It is important to note that in this experiment, we only estimate the principal components (and their variances) that explain the estimated non-linear manifold, i.e. matrix C by Algorithm 2, without reconstructing the denoised version of the corrupted data samples.

Both KPCA and our solution require model selection (choice of rank and τ respectively) which is beyond the scope of this paper. Here we resort to evaluate the performance of both methods under different parameters settings. To quantify the accuracy of the recovered manifold (C) we use following criteria: • Manifold Error : A good manifold should preserve maximum variance of the data -i.e. it should be able to generate a denoised version K(S est ) = C T C of the noisy kernel matrix K(S σ ). We define the manifold estimation error as K(S est ) -K(S GT ) 2 F , where K(S GT ) is the kernel matrix derived using noise free data. Figure 3 shows the manifold estimation error for KPCA and our method for different rank and parameter τ respectively.6 

• Classification error: The accuracy of a non-linear manifold is often also tested by the nearest neighbor classification accuracy. We select the estimated manifold which gives minimum Manifold Error for both the methods and report 1NN classification error (percentage of misclassified example) of the 1000 test points by projecting them onto estimated manifolds.

B KERNEL NRSFM WITH CAMERA POSE ESTIMATION:

Extended from section 4.2

Table 4 shows the reconstruction performance on a more realistic experimental setup, with the modification that the camera projection matrices are initialized with rigid factorization and were refined with the shapes by optimizing (2). To solve NRSfM problem with unknown projection matrices, we parameterize each R i with quaternions and alternate between refining the 3D shapes S and projection matrices R using LM. The regularization strength τ was selected for the TNH method by golden section search and parabolic interpolation for every test case independently. This ensures the best possible performance for the baseline. For our proposed approach τ was kept to 10 -4 for all sequences for both missing data and full data NRSfM. This experimental protocol somewhat disadvantages the non-linear method, since its performance can be further improved by a judicious choice of the regularization strength. However our purpose is primarily to show that the non-linear method adds value even without timeconsuming per-sequence tuning. To that end, note that despite large errors in the camera pose estimations by TNH and 50% missing measurements, the proposed method shows significant (∼ 10%) improvements in terms of reconstruction errors, proving our broader claims that non-linear representations are better suited for modeling real data, and that our robust dimensionality regularizer can improve inference for ill-posed problems.

As suggested by Dai et al. (2014) , robust camera pose initialization is beneficial for the structure estimation. We have used rigid factorization for initializing camera poses here but this can be trivially changed. We hope that further improvements can be made by choosing better kernel functions, with cross validation based model selection (value of τ ) and with a more appropriate tuning of kernel width. Selecting a suitable kernel and its parameters is crucial for success of kernelized algorithms. It becomes more challenging when no training data is available. We hope to explore other kernel functions and parameter selection criteria in our future work.

We would also like to contrast our work with Gotardo & Martinez (2011a), which is the only work we are aware of where non-linear dimensionality reduction is attempted for NRSfM. While estimating the shapes lying on a two dimensional non-linear manifold, Gotardo & Martinez (2011a) additionally assumes smooth 3D trajectories (parametrized with a low frequency DCT basis) and a pre-defined hard linear rank constraint on 3D shapes. The method relies on sparse approximation of the kernel matrix as a proxy for dimensionality reduction. The reported results were hard to replicate under our experimental setup for a fair comparison due to non-smooth deformations. However, in contrast to Gotardo & Martinez (2011a), our algorithm is applicable in a more general setup, can be modified to incorporate smoothness priors and robust data terms but more importantly, is flexible to integrate with a wide range of energy minimization formulations leading to a larger applicability beyond NRSfM.

C TNH ALGORITHM FOR NRSFM:

In section 4.2, we have compared the proposed non-linear dimensionality reduction prior against a variant of Garg et al. (2013a) which handles missing data by optimizing:

min S,R τ S * + F i=1 N j=1 Z i (x j ) w i (x j ) -R i s i (x j ) 2 (21)

This problem is convex in S given noise free projection matrix R i 's but non-differentiable. To optimize (21), we first rewrite it in its primal-dual form by dualizing the trace norm7 :

max Q min S,R τ < S, Q > + F i=1 N j=1 Z i (x j ) w i (x j ) -R i s i (x j ) 2 s.t. Q s ≤ 1 (22)

where Q ∈ R X ×N stores the dual variables to S and . s represent spectral norm (highest eigenvalue) of a matrix. We choose quaternions to perametrize the 2 × 3 camera matrices R i to satisfy orthonormality constraints as done in Garg et al. (2013a) and optimize the saddle point problem (22) using alternation.

In particular, for a single iteration: (i) we optimize the camera poses R i 's using LM, (ii) take a steepest descend step for updating S and (ii) a steepest ascend step for updating Q which is followed by projecting its spectral norm to unit ball. Given ground truth camera matrices ( without step (i)), alternation (ii-iii) can be shown to reach global minima of ( 22). Algorithm 3 outlines TNH algorithm.

As the main manuscript uses NRSfM only as a practical application of our non-linear dimensionality reduction prior, we have restricted our NRSfM experiments to only compare the proposed method against its linear counterpart. For the timely evaluation, the reported experiments we conducted on sub-sampled CMU mocap dataset. Here, we supplement the arguments presented in the main manuscript by favorably comparing the linear dimensionality reduction based NRSfM algorithm(TNH) to other NRSfM methods on full length CMU mocap sequences.

:

More precisely, M * was shown to be the tight convex envelope of rank(M )/ M s, where M s represent spectral norm of M .

1/τ can also be viewed as Lagrange multiplier to the constraints in (1).

Although it is clear that relaxing the rank of kernel matrix to K(S) * is suboptimal, works likeHuang et al. (2012);Cabral et al. (2013) with a variational definition of nuclear norm, allude to the possibility of kernelization. Further investigation is required to compare this counterpart to our tighter relaxation.

Since our main goal is to validate the usefulness of the proposed non-linear dimensionality regularizer, we opt for a reduced size dataset for more rapid and flexible evaluation.

Note that our formulation assumes Gaussian noise in K(S) where as for this evaluation we add noise to S directly.

Errors from non-noisy kernel matrix can be replaced by cross validating the entries of the kernel matrix for model selection for more realistic experiment.

For more details on primal dual formulation and dual norm of the trace norm seeRockafellar (1974);Recht et al. (2010);Chambolle & Pock (2011).

Figures and Tables:

Figure 1: Non-linear dimensionality regularizer for NRSfM. The top part of the figure explains the ill-posed inverse problem of recovering the causal factors (1); projection matrices Ri and 3D structures Si, from 2D image observations (2) Wi's, by minimizing the image reprojection error f (W, R, S) = i Wi -RiSi 2 . Assuming that the recovered 3D structures (Si's) lies near an unknown non-linear manifold (represented by the blue curve) in the input space, we propose to regularize the dimensionality of this manifold (3) -span of the non-linearly transformed shape vectors φ(Si)'sby minimizing φ(S) * . The non-linear transformation φ is defined implicitly with a Mercer kernel and maps the non-linear manifold to a linear low rank subspace (shown in blue line) of RKHS.

explains this central idea and proposed dimensionality regularizer in a nutshell with Non Rigid Structure from Motion (NRSfM) as the example application. Our Contribution: In this work we propose a unified for simultaneous robust KPCA and preimage estimation while solving an ill-posed inference problem without a pre-training stage.

Figure 2: Non-linear dimensionality regularisation improves NRSfM performance compared to its linear counterpart.Figure shows the ground truth 3D structures in red wire-frame overlaid with the structures estimated using: (a) proposed non-linear dimensionality regularizer shown in blue dots and (b) corresponding linear dimensionality regularizer (TNH) shown in black crosses, for sample frames of CMU mocap sequence. Red circles represent the 3D points for which the projections were known whereas squares annotated missing 2D observations. See text and Table 2 for details.

Figure 2: Non-linear dimensionality regularisation improves NRSfM performance compared to its linear counterpart.Figure shows the ground truth 3D structures in red wire-frame overlaid with the structures estimated using: (a) proposed non-linear dimensionality regularizer shown in blue dots and (b) corresponding linear dimensionality regularizer (TNH) shown in black crosses, for sample frames of CMU mocap sequence. Red circles represent the 3D points for which the projections were known whereas squares annotated missing 2D observations. See text and Table 2 for details.

Recent work, likeDai et al. (2014);Garg et al.  (2013a)  have shown that the trace norm regularizer can be used as a convex envelope of the low-rank prior to robustly address ill-posed nature of the problem. A good solution to NRSfM can be achieved by optimizing:

Figure 3: Performance comparison between KPCA and our Robust closed-form solution with dimensionality regularization on oil flow dataset with additive Gaussian noise of standard deviation σ. Plots show the normalized kernel matrix errors with different rank of the model. Kernel PCA results are shown in dotted line with diamond while ours are with solid line with a star. Bar-plot show the worst and the best errors obtained by our method for a single rank of recovered kernel matrix.

Trace norm Heuristics.Input: Initial estimates S 0 , R 0 of S and R. Output: Low-dimensional S and camera poses R. Parameters: Regularization strength τ , measurements W and binary mask Z.-S = S 0 , R = R 0 ; // set iteration count η step size σ and duals Q η = 0; -σ = 1/τ ; -Q = 0;while not converged do // projection matrix estimation -Fix S, Q and refine Ri for every image i with LM;// steepest descend update for Sij for each point xj and each frame ifor i = 1 to F do for j = 1 to N do -S η+1 ij = I2×2 + σ(ZijR T i Ri) -1 (S η ij -στ Q η ij + σR T i (Zij • wij));// accelerated steepest ascend update for Q-Q = Q η + τ σ(2S n+1 -S n ); -U DV T = singular value decomposition of Q η ; -D = min(D, 1); -Q η+1 = U DV T ;// Go to next iteration η = η + 1

Performance comparison on missing data completion on Oil Flow Dataset: Row 1 shows the amount of missing data and subsequent rows show the mean and standard deviation of the error in recovered data matrix over 50 runs on 100 samples of oil flow dataset by: (1) The mean method (also the initialization of other methods) where the missing entries are replaced by the mean of the known values of the corresponding attributes, (2) 1-nearest neighbor method in which missing entries are filled by the values of the nearest point, generalization ofZongben et al. (2012) to include the L 1/2 matrix norms for which a simplified proof is included in the Appendix A. It is important to note however, that the motivation and implication of using L 1/2 regularization in the context of non-linear dimensionality reduction are significantly different to that of Zongben et al. (2012) and related workDu et al. (2013);Zhao et al. (2014) which are designed for linear modeling of the causal factors. The core insight of using L 1 regularization in the feature space via the parametrization given in 3 facilitates a natural way for non-linear modeling of causal factors with low dimensionality while solving an inverse problem by making feature space tractable.

Table content: <table><tr><td colspan="5">(3) PPCA Tipping &amp; Bishop (1999), (4) PKPCA of Sanguinetti &amp; Lawrence (2006), (5)RKPCA Nguyen &amp; De</td></tr><tr><td>la Torre (2009) and our method.</td><td/><td/><td/><td/></tr><tr><td>p(del)</td><td>0.05</td><td>0.10</td><td>0.25</td><td>0.50</td></tr><tr><td>mean</td><td>13 ± 4</td><td colspan="2">28 ± 4 70 ± 9</td><td>139 ± 7</td></tr><tr><td>1-NN</td><td>5 ± 3</td><td colspan="2">14 ± 5 90 ± 20</td><td>NA</td></tr><tr><td>PPCA</td><td>3.7 ± .6</td><td colspan="3">9 ± 2 50 ± 10 140 ± 30</td></tr><tr><td>PKPCA</td><td>5 ± 1</td><td colspan="3">12 ± 3 32 ± 6 100 ± 20</td></tr><tr><td colspan="3">RKPCA 3.2 ± 1.9 8 ± 4</td><td>27 ± 8</td><td>83 ± 15</td></tr><tr><td>Ours</td><td>2.3±2</td><td>6±3</td><td>22±7</td><td>70±11</td></tr><tr><td>be seen as</td><td/><td/><td/><td/></tr></table>

3D reconstruction errors for linear and non-linear dimensionality regularization with ground truth camera poses. Column 1 and 4 gives gives error for TNH while column (2-3) and (5-6) gives the corresponding error for proposed method with different width of RBF kernel. Row 5 reports the mean error over 4 sequences.We use a sub-sampled version of CMU mocap dataset by selecting every 10 th frame of the smoothly deforming human body consisting 41 mocap points used inDai et al. (2014). 4

Table content: <table><tr><td/><td/><td colspan="2">No Missing Data</td><td/><td colspan="2">50% Missing Data</td></tr><tr><td>Dataset</td><td>Linear</td><td colspan="2">Non-Linear d max d med</td><td>Linear</td><td colspan="2">Non-Linear d max d med</td></tr><tr><td>Drink</td><td>0.0227</td><td>0.0114</td><td>0.0083</td><td>0.0313</td><td>0.0248</td><td>0.0229</td></tr><tr><td>Pickup</td><td>0.0487</td><td>0.0312</td><td>0.0279</td><td>0.0936</td><td>0.0709</td><td>0.0658</td></tr><tr><td>Yoga</td><td>0.0344</td><td>0.0257</td><td>0.0276</td><td>0.0828</td><td>0.0611</td><td>0.0612</td></tr><tr><td>Stretch</td><td>0.0418</td><td>0.0286</td><td>0.0271</td><td>0.0911</td><td>0.0694</td><td>0.0705</td></tr><tr><td>Mean</td><td>0.0369</td><td>0.0242</td><td>0.0227</td><td>0.0747</td><td>0.0565</td><td>0.0551</td></tr><tr><td colspan="3">4.2.1 RESULTS ON THE CMU DATASET</td><td/><td/><td/><td/></tr></table>

Robust dimensionality reduction accuracy by KPCA versus our closed-form solution on the full oil flow dataset. Columns from left to right represent: (1) standard deviation of the noise in training samples (2-3) Error in the estimated low-dimensional kernel matrix by (2) KPCA and (3) our closed-form solution, (4-5) Nearest neighbor classification error of test data using (4) KPCA and (5) our closed-form solution respectively.

Table content: <table><tr><td/><td/><td colspan="3">Manifold Error</td><td/><td colspan="4">Classification Error</td></tr><tr><td colspan="2">STD</td><td colspan="4">KPCA Our CFS</td><td colspan="4">KPCA Our CFS</td></tr><tr><td colspan="2">.2</td><td>0.1099</td><td/><td>0.1068</td><td/><td colspan="2">9.60%</td><td colspan="2">9.60%</td></tr><tr><td colspan="2">.3</td><td>0.2298</td><td/><td>0.2184</td><td/><td colspan="4">19.90% 15.70%</td></tr><tr><td colspan="2">.4</td><td>0.3522</td><td/><td>0.3339</td><td/><td colspan="4">40.10% 22.20%</td></tr><tr><td/><td>0.4</td><td/><td/><td/><td/><td/><td/><td/></tr><tr><td/><td>0.35</td><td/><td/><td/><td/><td/><td/><td/></tr><tr><td>Manifold error</td><td>0.2 0.25 0.3</td><td/><td/><td/><td/><td/><td/><td>KPCA,σ=.2 Ours,σ=.2</td></tr><tr><td/><td/><td/><td/><td/><td/><td/><td/><td>KPCA,σ=.3</td></tr><tr><td/><td/><td/><td/><td/><td/><td/><td/><td>Ours,σ=.3</td></tr><tr><td/><td>0.15</td><td/><td/><td/><td/><td/><td/><td>KPCA,σ=.4</td></tr><tr><td/><td/><td/><td/><td/><td/><td/><td/><td>Ours,σ=.4</td></tr><tr><td/><td>0.1</td><td/><td/><td/><td/><td/><td/><td/></tr><tr><td/><td>0</td><td>2</td><td>4</td><td>6</td><td>8</td><td>10</td><td>12</td><td>14</td><td>16</td></tr><tr><td/><td/><td/><td/><td colspan="3">Rank of kernel matrix</td><td/><td/></tr></table>

3D reconstruction errors for linear and non-linear dimensionality regularization with noisy camera pose initialization from rigid factorization and refined in alternation with shape. The format is same as Table2.

Table content: <table><tr><td/><td/><td colspan="2">No Missing Data</td><td/><td colspan="2">50% Missing Data</td></tr><tr><td>Dataset</td><td>Linear τ = τ  *</td><td colspan="2">Non-Linear τ = 10 -4</td><td>Linear τ = τ  *</td><td colspan="2">Non-Linear τ = 10 -4</td></tr><tr><td/><td/><td>dmax</td><td>d med</td><td/><td>dmax</td><td>d med</td></tr><tr><td>Drink</td><td>0.0947</td><td>0.0926</td><td>0.0906</td><td>0.0957</td><td>0.0942</td><td>0.0937</td></tr><tr><td>Pickup</td><td>0.1282</td><td>0.1071</td><td>0.1059</td><td>0.1598</td><td>0.1354</td><td>0.1339</td></tr><tr><td>Yoga</td><td>0.2912</td><td>0.2683</td><td>0.2639</td><td>0.2821</td><td>0.2455</td><td>0.2457</td></tr><tr><td>Stretch</td><td>0.1094</td><td>0.1043</td><td>0.1031</td><td>0.1398</td><td>0.1459</td><td>0.1484</td></tr><tr><td>Mean</td><td>0.1559</td><td>0.1430</td><td>0.1409</td><td>0.1694</td><td>0.1552</td><td>0.1554</td></tr></table>

3D reconstruction errors for different NRSfM approaches and our TNH Algorithm given ground truth camera projection matrices. Results for all the methods (except TNH) are taken fromDai et al. (2014).

Table content: <table><tr><td>Dataset</td><td colspan="3">PTAAkhter et al. (2009) CSF2Gotardo &amp; Martinez (2011b) BMMDai et al. (2014)</td><td>TNH</td></tr><tr><td>Drink</td><td>0.0229</td><td>0.0215</td><td>0.0238</td><td>0.0237</td></tr><tr><td>Pick-up</td><td>0.0992</td><td>0.0814</td><td>0.0497</td><td>0.0482</td></tr><tr><td>Yoga</td><td>0.0580</td><td>0.0371</td><td>0.0334</td><td>0.0333</td></tr><tr><td>Stretch</td><td>0.0822</td><td>0.0442</td><td>0.0456</td><td>0.0431</td></tr></table>

Peer review (iclr2017):
The paper proposes a nonlinear regularizer for solving ill-posed inverse problems. The latent variables (or causal factors) corresponding to the observed data are assumed to lie near a low dimensional subspace in an RKHS induced by a predetermined kernel. The proposed regularizer can be seen as an extension of the linear low-rank assumption on the latent factors. A nuclear norm penalty on the Cholesky factor of the kernel matrix is used as a relaxation for the dimensionality of the subspace. Empirical results are reported on two tasks involving linear inverse problems -- missing feature imputation, and estimating non-rigid 3D structures from a sequence of 2D orthographic projections -- and the proposed method is shown to outperform linear low-rank regularizer. The clarity of the paper has scope for improvement (particularly, Introduction) - the back and forth b/w dimensionality reduction techniques and inverse problems is confusing at times. Clearly defining the ill-posed inverse problem first and then motivating the need for a regularizer (which brings dimensionality reduction techniques into the picture) may be a more clear flow in my opinion. The motivation behind relaxation of rank() in Eq 1 to nuclear-norm in Eq 2 is not clear to me in this setting. The relaxation does not yield a convex problem over S,C (Eq 5) and also increases the computations (Algo 2 needs to do full SVD of K(S) every time). The authors should discuss pros/cons over the alternate approach that fixes the rank of C (which can be selected using cross-validation, in the same way as is selected), leaving just the first two terms in Eq 5. For this simpler objective, an interesting question to ask would be -- are there kernel functions for which it can solved in a scalable manner? The proposed alternating optimization approach in the current form is computationally intensive and seems hard to scale to even moderate sized data -- in every iteration one needs to compute the kernel matrix over S and perform full SVD over the kernel matrix (Algo 2). Empirical evaluations are also not extensive -- (i) the dataset used for feature imputation is old and non-standard, (ii) for structure estimation from motion on CMU dataset, the paper only compares with linear low-rank regularization, (iii) there is no comment/study on the convergence of the alternating procedure (Algo 1).