# Patent

## Title
Method for mapping high-dimensional samples to reduced-dimensional manifolds

## Abstract

A method determines mappings between high-dimensional measured samples and a reduced-dimensional manifold. Samples are acquired of a physical system. The samples have a high number of dimensions. A low number of dimensions are determined for a manifold embedded in a space of the high-dimensional samples. Local charts having the dimensions of the low number of dimensions of the manifold are fitted to selected high-dimensional samples. The charts are then connected to determine a forward mapping from any high-dimensional sample to a coordinate on the manifold and a reverse mapping from any coordinate on the manifold to a corresponding point in high-dimensional sample space.

## Background
FIELD OF THE INVENTION
This invention relates generally to sample analysis, and more particularly to reducing the dimensionality of high-dimensional sample spaces and mapping the samples to a lower dimensional sample space.
BACKGROUND OF THE INVENTION
It is often the case that a physical system or object having a modest number of degrees of freedom is measured with an apparatus that has many more degrees of freedom than the system or object. For example, a human face has a few hundred (d) degrees of freedom, while a single digital mega-pixel image of the human face has many million (D) degrees of freedom, for example, i.e., at least three color intensities for each pixel, and thus d<<D.
The space of all possible images of varied facial expressions is thus a low-dimensional manifold embedded in the high-dimensional sample space of all possible digital images. Herein, a manifold is exemplified by a smooth 2D surface embedded in a 3D space (coordinate system), and more generally defined as a connected subregion of a high-dimensional space that can be flattened to yield a lower-dimensional space that preserves the connectivity structure of local neighborhoods in the subregion.
In the context of the present invention, a manifold is a representation of all possible configurations of a measured physical system or object. For example, embedded in the high-dimensional space of all possible images of a hand there is a 22-dimensional manifold representing all possible joint angles. The invention shows how to estimate and use this manifold and its embedding from a finite number of measurements, e.g., a set of digital images.
In the field of mathematics known as differential geometry, charting is known as the problem of assigning a local low-dimensional coordinate system to every point on a manifold whose geometry is known analytically, e.g., characterized by mathematical formulae. Of particular interest is the case where these local charts can be connected to give a globally consistent chart, which assigns a unique low-dimensional coordinate.
In the present invention, it is shown how to perform analogous operations on measured data in a high-dimensional sample space of real-world physical systems. The invention uses the chart to edit, visualize or render, compress, de-noise, and otherwise manipulate the measured signal while preserving the fidelity of the results.
With charting, it is presumed that the data lie on or near some low-dimensional manifold embedded in the sample space, and that there exists a one-to-one smooth non-linear transform between the manifold and a low-dimensional vector space.
It is desired to estimate smooth continuous mappings between the sample and coordinate spaces. Often this analysis sheds light on intrinsic variables of the data-generating phenomenon, for example, revealing perceptual or configuration spaces.
Topology-neutral non-linear dimensionality reduction processes can be divided into those that determine mappings, and those that directly determine low-dimensional embeddings.
DeMers et al., in “Nonlinear Dimensionality Reduction,” Advances in Neural Information Processing Systems, Volume 5, 1993, described how auto-encoding neural networks with a hidden layer “bottleneck,” can be used to effectively cast dimensionality reduction as a compression problem.
Hastie et al., defined principal curves as non-parametric ID curves that pass through the center of “nearby” data samples in “Principal Curves,” J. Am. Statistical Assoc, 84(406):502-516, 1989.
A number of known methods properly regularize that approach and extend it to surfaces, see Smola et al., in “Regularized Principal Manifolds,” Machine Learning, 1999. They provided an analysis of such methods in the broader framework of regularized quantization methods, and developed an iterative process as a means of investigating sample complexity.
A number of significant advances have been made recently in embedding processes. Gomes et al., in “A variational approach to recovering a manifold from sample points,” Proc. ECCV, 2002, treated manifold completion as an anisotropic diffusion problem where sample points are iteratively expanded until they connect to their neighbors.
The ISOMAP process represented remote distances as sums of a trusted set of distances between immediate neighbors, and used multi-dimensional scaling to determine a low-dimensional embedding that minimally distorts all distances, see Tenenbaum et al., “A Global Geometric Framework For Nonlinear Dimensionality Reduction,” Science, 290:2319-2323, Dec. 22, 2000.
A locally linear embedding (LLE) method represented each sample point as a weighted combination of its trusted nearest neighbors, then determined a low-dimensional embedding that minimally distorted these barycentric coordinates, see Roweis et al., “Nonlinear Dimensionality Reduction By Locally Linear Embedding,” Science, 290:2323-2326, Dec. 22, 2000.
Those methods have complementary strengths: ISOMAP handles missing data well but fails if the data hull are non-convex; and vice versa for the LLE method. Both offer embeddings without mappings. However, trusted-set methods are vulnerable to noise because they only consider the subset of point-to-point relationships that tend to have the lowest signal-to-noise relationship. It is known that trusted methods can be unstable because small changes to the size of the trusted neighbor set can induce large changes in the set of constraints that the low-dimensional embedding must satisfy.
In a return to mapping, Roweis et al., in “Global coordination of linear models, “Advances in Neural Information Processing Systems,” Volume 13, 2002, described global coordination. In that method, a mixture of locally linear projections from sample to coordinate space are learned. They devised a posterior probability measure that penalizes distortions in the mapping, and gave an expectation-maximization (EM) rule for locally maximizing the posterior probability.
That work heralded an intriguing new use of variational methods to learn representations, but also highlighted the difficulty of even hill-climbing their multi-modal posterior probability.
It is desired to provide a method for charting a manifold that is a closed-form solution and avoids hill-climbing a non-uniform posterior probability, and to find mappings between sample space and charting space.
SUMMARY OF THE INVENTION
The invention provides non-linear forward and reverse mappings between a high-dimensional (D) sample space and a low-dimensional (d) manifold. The mappings effectively recover a Cartesian coordinate system for the manifold from which the data are sampled. The mappings preserve local distances and angles in the manifold, and are pseudo-invertible, smooth, continuous and well-behaved in the presence of sampling bias and noise.
The invention estimates an intrinsic low dimensionality of the manifold, decomposes the data sample into locally linear low-dimensional patches or “charts,” merges these charts into a single low-dimensional coordinate system, and determines forward and reverse mappings between the sample and coordinate spaces. All solutions are unique, optimal, and of a closed form.
First, the intrinsic dimensionality (d) of the measured data manifold is estimated from a growth process of point-to-point distances. Then the invention treats the d-dimensional charting problem from the perspective of minimizing information loss. Then, the problem of mapping data-points to a vector space or manifold is performed by parameterizing a kernel-based mixture of linear projections. The posterior probability combines a standard Gaussian mixture model (GMM) likelihood function with a prior probability that penalizes uncertainty in the locations of a projected point due to inconsistent projections.
For the special case where the posterior probability is unimodal and maximizable in closed form, the invention provides a GMM whose covariances reveal local linear subspaces on the manifold. With the invention method, any GMM and choice of reduced dimension d has a unique, closed-form solution for a minimal-distorting non-linear forward projection of the data from the sample space to the d-dimensional coordinate space, as well as an inverse mapping defining the surface of the manifold in the sample space.
In analogy to differential geometry, the subspaces are called “charts” and their merger the “connection.” The method can be used to untie knots, unroll and untwist sheets, and to visualize, compress and model high-dimensional bio-mechanical data, such as digital images of human body parts.



BRIEF DESCRIPTION OF THE DRAWINGS
FIG. 1 is a flow diagram of a method for representing high-dimensional measured data samples by reduced-dimensional manifolds;
FIG. 2 is a graph of a point-count growth process;
FIG. 3 is a graph comparing slopes of growth curves; and
FIG. 4 is a block diagram of a chart.



DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT
A method 100 according to the invention begins by measuring or otherwise acquiring 110 samples Y 111 of a physical system or object 101. For example, the samples 111 are pixel intensities in a digital image of a face.
The space (coordinate system) of the samples characteristically has a very high number (D) of dimensions. The invention provides a method for constructing a low-dimensional (d) manifold which approximates the distribution of data points in the sample space, and determining the intrinsic d-dimensional coordinate system of the manifold. The coordinate system is assigned so that every point on the manifold has a unique low-dimensional set of ordinates, and nearby points have similar ordinates.
The method also produces functions ƒ and ƒ1 141-142 that respectively map forward and backward between the original high-dimensional sample space and the low-dimensional coordinate system of the manifold. These mapping functions support operations such as signal compression, de-noising, interpolation, extrapolation, and synthesis 151-156, described in further detail below.
If the samples are pixel intensities in a digital image of a human faces, then the images can be represented by the low-dimensional coordinates of their corresponding points on the manifold, allowing high-resolution images to be represented with a just few hundred numbers. In other words, the sample space or images are greatly compressed 151. Furthermore, images of arbitrary facial expressions can be synthesized 152 from arbitrary points on the manifold.
A path on the manifold connecting two points representing two facial expressions give a sequence of facial images that interpolates 153 between the two expressions, producing a physically plausible video. A straight line in the recovered coordinate system gives one such path. However, a straight line connecting two points in the measured image space gives a useless sequence of blends between the two images.
The method can also be used to correct 154 a noisy or corrupted samples, e.g. images. Such images are points in image space that do not lie on the manifold. To de-noise a corrupted image, find a closest point on the manifold to the corrupted data, and synthesize a corresponding image pixel from that point.
The distribution of all samples, e.g., images, can be visualized (rendered) 155 in 2D or 3D by plotting each sample point according to the first two or three ordinates of its manifold coordinates. This can reveal trends and groupings that are impossible to see in the original sample space.
Statistical analyses 156 that are not tractable in the original high-dimensional sample space become feasible in the low-dimensional coordinate system of the manifold.
Estimating Locally Linear Scale and Intrinsic Dimensionality
As stated above, the invention begins with matrix of measured or otherwise acquired sample points Y 111 of a physical system or object 101. For example, the sample points are pixel intensity values of digital images. The matrix is represented by

    
    
        Y{dot over (=)}[y1, . . . , yN], ynεD 
    
    


populating a D-dimensional sample space, where D is very large, 106<<D. It is desired to fit a d-dimensional manifold M to this sample space so that d<<D, e.g. d=102 for a reduction of approximately four orders of magnitude.
The invention finds a non-linear forward projection ƒ 141 to a subspace

    
    
        G(Y)→X{dot over (=)}[x1, . . . , xN], xn=εd,
    
    


and a one-to-one inverse map ƒ1 142

    
    
        G−1(X)→Y
    
    


such that local relations between nearby points on the manifold are preserved. It is desired that the forward map G is non-catastrophic, that is, smooth continuous and without sharp folds. Parallel lines on the manifold in D should map to continuous smooth non-intersecting curves in d.
Smoothness means that at some scale r the mapping from a neighborhood on M to d is effectively linear. As one approaches the limit of infinite noiseless samples, finding r is done as follows.
Consider a sphere of radius r centered on an arbitrary data point. Let Σr be the covariance of the nr points within that sphere. On the true manifold, both nr and the volume of Σr grow as the radius rd of the sphere increases, and therefore at very small scales the rank of Σr is d.
When dealing with finite samples, the scale at which the manifold is locally linear is the interval 0≦r≦R for which rank (Σr)≦d. At larger scales the curvature of the embedded manifold inflates the rank. This is complicated by the presence of sample noise of standard deviation σ; for 0≦r≦σ, rank(Σr)=D>d. In practical settings, linearity is obtained in the scale interval σ≦r≦R.
The immediate problem is to estimate r not knowing σ or R. From the volume growth process, the determinant |Σr|∝rd. Thus it follows that there exists a scale interval in which v(r){dot over (=)}d/drlog|Σr| satisfies v(r)=d/r. For reasons of computational and sample complexity, tracking the growth of the covariance at every point is not a practical procedure. However, it is practical to track the point-count growth process c(r){dot over (=)}dlog r/d lognr which behaves similarly and peaks at c(r)=1/d when r is at the locally linear scale. The point-count growth process on a 1D manifold noisily embedded in a 2D space is shown in FIG. 2 where the samples 111 are dots.
Characteristic curves of the process are shown in FIG. 3 for radial growth and 1D-3D hypothesises. The point-count growth process is used to find the intrinsic dimensionality of a 2D manifold in three -space. Lines of slope 1/3, 1/2, and 1/1 are fitted to sections of the log r/lognr curve. For neighborhoods of radius r ̃1 with roughly n≈10 points, the slope peaks at 1/2 indicating a dimensionality of d=2. Below that, the data appears 3D because it is dominated by noise, except for n≦D points. Above, the data appears as >2D because of manifold curvature. As the r-sphere expands to cover the entire data-set the dimensionality appears to drop to 1 as the process begins to track the edges of the sheet.
At sub-noise scales 211, c(r<σ)≈1/D<1/d, because noise has distributed points in all directions with equal probability. At the locally linear scale 212, c(r)=1/d, because points are distributed only in the directions of the local tangent space of the manifold. At the scale 213 at which curvature becomes significant, c(r>R)<1/d, because the manifold is no longer perpendicular to the surface of the sphere, so the sphere does not have to grow as fast to accommodate new points.
Therefore, the maximum slope of c(r) gives a good estimate of the local dimensionality d of the manifold, provided that the sphere has not expanded to a manifold boundary, i.e., edges and corners, because the boundary has lower dimension than the manifold. Even for low-dimensional manifolds such as sheets, the boundary manifolds are very small relative to the full manifold, so the boundary effect is typically limited to a small rise in c(r) as r approaches the scale of the entire sample space. In practice, the method according to the invention simply expands a sphere of radius r at every point or a random selection of points and looks for the first peak in c(r). The values d and r can be estimated globally, per-point.
Fitting Local Charts
For the charting fitting step 130, the method 100 finds a soft partitioning of the data samples into neighborhoods N 131. As shown in FIG. 4, the charts can be expressed as a Gaussian, or elliptically shaped sphere 401, with the dominant axes 410 lying in the plane of the manifold, and the minor axis 420 perpendicular to the manifold. The charts are locally linear and low-dimensional, in the sense that the data project into a subspace associated with each neighborhood with (1) minimal loss of variance and (2) maximal agreement of the projections of nearby points into nearby neighborhoods. Criterion (1) is expressed by maximizing the likelihood function p of a Gaussian mixture model (GMM) density fitted to the data samples: 



  
    
      
        
          p
          ⁡
          
            (
            
              
                y
                i
              
              |
              
                μ
                ⁢
                
                   
                
                ⁢
                Σ
              
            
            )
          
        
        ⁢
        
          =
          .
        
        ⁢
        
          
            
              ∑
              j
            
            ⁢
            
              
                p
                (
                
                  
                    
                      y
                      i
                    
                    |
                    
                      μ
                      j
                    
                  
                  ,
                  
                    Σ
                    j
                  
                
                )
              
              ⁢
              
                p
                j
              
            
          
          =
          
            
              ∑
              j
            
            ⁢
            
              
                N
                (
                
                  
                    
                      y
                      i
                    
                    ;
                    
                      μ
                      j
                    
                  
                  ,
                  
                    Σ
                    j
                  
                
                )
              
              ⁢
              
                
                  p
                  j
                
                .
              
            
          
        
      
    
    
      
        (
        1
        )
      
    
  




Each Gaussian component defines a local neighborhood N centered around μj with axes defined by the eigenvectors of Σj. The amount of data variance along each axis is indicated by the eigenvalues of Σj. If the manifold is locally linear in the vicinity of the μj, then all but the d dominant eigenvalues are near-zero. This implies that the associated eigenvectors constitute the optimal variance-preserving local coordinate system.
To some degree likelihood maximization naturally realize this property. Likelihood maximization requires that the GMM components shrink in volume to fit the data as tightly as possible. This shrinking is best accommodated by positioning the components so that they project onto locally flat collections of data points. However, this is easily violated by degenerate zero-variance GMM components or components fitted to small enough locales so that data density not on the manifold is comparable to density on the manifold, e.g., at the noise scale. Consequently a prior probability is needed.
Criterion (2) implies that neighboring partitions should have dominant axes that span similar subspaces, because disagreement, i.e., large subspace angles, would lead to inconsistent projections of the same point and therefore uncertainty about its location in low-dimensional coordinate space. The principal insight is that criterion (2) is exactly the cost of coding the location of a point in one neighborhood when it is generated by another neighborhood. The cross-entropy between the Gaussian models defining the two neighborhoods is 



  
    
      
        
          D
          ⁡
          
            (
            
              
                N
                1
              
              ||
              
                N
                2
              
            
            )
          
        
        =
        
        ⁢
        
          ∫
          
            
              d
              y
            
            ⁢
            
               
            
            ⁢
            
              N
              ⁡
              
                (
                
                  
                    y
                    ;
                    
                      μ
                      1
                    
                  
                  ,
                  
                    Σ
                    1
                  
                
                )
              
            
            ⁢
            
               
            
            ⁢
            log
            ⁢
            
              
                N
                ⁡
                
                  (
                  
                    
                      y
                      ;
                      
                        μ
                        1
                      
                    
                    ,
                    
                      Σ
                      1
                    
                  
                  )
                
              
              
                N
                ⁡
                
                  (
                  
                    
                      y
                      ;
                      
                        μ
                        2
                      
                    
                    ,
                    
                      Σ
                      2
                    
                  
                  )
                
              
            
          
        
      
    
    
       
    
    
       
    
    
       
    
    
       
    
    
      
              
        ⁢
        
          (
          2
          )
        
      
    
  
  
    
      
        =
        
        ⁢
        
          
            log
            ⁢
            
               
            
            ⁢
            
              
              
                
                  Σ
                  1
                  
                    -
                    1
                  
                
                ⁢
                
                  Σ
                  2
                
              
              
            
          
          +
          
            trace
            (
            
              
                Σ
                1
                
                  -
                  1
                
              
              ⁢
              
                Σ
                2
              
            
            )
          
          +
        
      
    
    
       
    
    
       
    
    
       
    
    
       
    
    
      
        
        ⁢
        
          (
          3
          )
        
      
    
  
  
    
      
        
        ⁢
        
          
            
              
                (
                
                  
                    μ
                    1
                  
                  -
                  
                    μ
                    2
                  
                
                )
              
              T
            
            ⁢
            
              
                Σ
                1
                
                  -
                  1
                
              
              ⁡
              
                (
                
                  
                    μ
                    1
                  
                  -
                  
                    μ
                    2
                  
                
                )
              
            
          
          +
          
            D
            /
            2.
          
        
      
    
    
       
    
    
       
    
    
       
    
    
       
    
    
      
         
      
    
  




Roughly speaking, the terms on the second line measure differences in size, orientation, and position of the charts 131, respectively, of two coordinate frames located at the means μ1, μ2 with axes specified by the eigenvectors of Σ1, Σ2. All three terms decline to zero as the overlap between the two frames is maximized. To maximize consistency between adjacent neighborhoods, the prior

    
    
        a. p(μ, Σ){dot over (=)}exp−Σl≠jml(yl)D(N1∥N2),
    
    


where ml(yj) is a measure of co-locality.
The dominant axes in the neighboring charts 131 do not need to aligned, they only need to span nearly the same subspace. This is a much easier objective to satisfy, and it contains a useful special case where the posterior p(μ,Σ|Y)∝Σip(yi|μΣ)p(μ;Σ) is unimodal and can be maximized in closed form.
A Gaussian neighborhood chart is associated with each data-point, setting μi=yi. All neighborhoods are taken to be a priori equally probable, setting pi=1/N, and define the co-locality measure mi(y)∝N(y;yi;σ2), with the scale parameter σ specifying the expected size of a neighborhood on the manifold in sample space. A reasonable choice is σ=r/2, so that 2erf(2)>99.5% of the density of mi(y) is contained in the area around yi where the manifold is expected to be locally linear. With μi, pi and mi(y) fixed, the MAP estimates of the GMM covariances are 



  
    
      
        
          
            
              
                Σ
                i
              
              =
              
              ⁢
              
                
                  ∑
                  j
                
                ⁢
                
                  
                    
                      m
                      i
                    
                    ⁡
                    
                      (
                      
                        y
                        j
                      
                      )
                    
                  
                  ⁢
                  
                    (
                    
                      
                        
                          (
                          
                            
                              y
                              j
                            
                            -
                            
                              μ
                              i
                            
                          
                          )
                        
                        ⁢
                        
                          
                            (
                            
                              
                                y
                                j
                              
                              -
                              
                                μ
                                i
                              
                            
                            )
                          
                          T
                        
                      
                      +
                    
                  
                
              
            
          
        
        
          
            
              
                
                ⁢
                
                  
                    
                      (
                      
                        
                          μ
                          j
                        
                        -
                        
                          μ
                          i
                        
                      
                      )
                    
                    ⁢
                    
                      
                        (
                        
                          
                            μ
                            j
                          
                          -
                          
                            μ
                            i
                          
                        
                        )
                      
                      T
                    
                  
                  +
                  
                    Σ
                    j
                  
                
                )
              
              /
              
                
                  ∑
                  j
                
                ⁢
                
                  
                    m
                    i
                  
                  ⁡
                  
                    (
                    
                      y
                      j
                    
                    )
                  
                
              
            
          
        
      
    
    
      
        (
        4
        )
      
    
  




Note that each covariance Σi is dependent on all other Σj. The MAP estimators for all covariances can be arranged into a set of fully constrained linear equations and solved exactly for their mutually optimal values. This key step brings non-local information about the shape of the manifold into the local description of each neighborhood.
This step also ensures that adjoining neighborhoods have similar covariances, and thus small angles between their respective subspaces, particularly subspaces of dimension d or higher. Even if the local set of data points are dense in a direction perpendicular to the manifold, the local charts orient parallel to the manifold as part of a globally optimal solution.
Connecting the Charts
The next step 140 constructs a connection for any set of charts specified as an arbitrary non-degenerate GMM. The GMM yields a soft partitioning of the dataset into neighborhoods of mean μk and covariance μk. The optimal variance-preserving low-dimensional coordinate system for each neighborhood derives from its weighted principal component analysis, which is exactly specified by the eigenvectors of its covariance matrix.
Eigendecompose VkΛkVTk←Σk with eigenvalues in descending order on the diagonal of Λk, and define the chart projection Wk{dot over (=)}[Id,0]VTk, which extracts the eigenvectors in Vk corresponding to the d dominant eigenvalues.
Define per-chart probability py|k(y)VTk{dot over (=)}N(y;μk, Σk), local chart coordinate μki{dot over (=)}Wk(yi−μk), and let Uk{dot over (=)}(Uk1, . . . , UkN) be the local coordinates of all points in chart k.
The goal is to connect all charts into a globally consistent low-dimensional coordinate system. For each chart, there is a low-dimensional affine transform Gk that projects Uk into the global coordinate system. Integrating over all charts, the weighted average of the projections of point yi into the low-dimensional vector space is 



  
    
      
        
          
            x
            |
            y
          
          ^
        
        ⁢
        
          =
          .
        
        ⁢
        
          
            
              
                ∑
                j
              
              ⁢
              
                
                  
                    G
                    j
                  
                  ⁡
                  
                    [
                    
                      
                        
                          
                            
                              W
                              j
                            
                            ⁡
                            
                              (
                              
                                y
                                -
                                
                                  μ
                                  j
                                
                              
                              )
                            
                          
                        
                      
                      
                        
                          1
                        
                      
                    
                    ]
                  
                
                ⁢
                
                   
                
                ⁢
                
                  
                    p
                    
                      y
                      |
                      j
                    
                  
                  ⁡
                  
                    (
                    y
                    )
                  
                
                ⁢
                
                   
                
                ⁢
                
                  p
                  j
                
              
            
            ⇒
            
              
                
                  x
                  i
                
                |
                
                  y
                  i
                
              
              ^
            
          
          ⁢
          
            =
            .
          
          ⁢
          
            
              ∑
              i
            
            ⁢
            
              
                
                  G
                  j
                
                ⁡
                
                  [
                  
                    
                      
                        
                          u
                          ji
                        
                      
                    
                    
                      
                        1
                      
                    
                  
                  ]
                
              
              ⁢
              
                
                  p
                  
                    y
                    |
                    j
                  
                
                ⁡
                
                  (
                  
                    y
                    i
                  
                  )
                
              
              ⁢
              
                
                  p
                  j
                
                .
              
            
          
        
      
    
    
      
        (
        5
        )
      
    
  




If a point has non-zero probabilities in two charts, then the affine 



  
    
      
        
          
            
              G
              ⁢
              
                =
                .
              
              ⁢
              
              ⁢
              
                [
                
                  
                    G
                    1
                  
                  ,
                  ⋯
                  ⁢
                  
                     
                  
                  ,
                  
                    G
                    K
                  
                
                ]
              
            
          
        
        
          
            
              =
              
              ⁢
              
                arg
                ⁢
                
                   
                
                ⁢
                
                  
                    min
                    
                      
                        G
                        k
                      
                      ,
                      
                        G
                        j
                      
                    
                  
                  ⁢
                  
                    
                      ∑
                      i
                    
                    ⁢
                    
                      
                        
                          p
                          
                            y
                            |
                            k
                          
                        
                        ⁡
                        
                          (
                          
                            y
                            i
                          
                          )
                        
                      
                      ⁢
                      
                        
                          p
                          
                            y
                            |
                            j
                          
                        
                        ⁡
                        
                          (
                          
                            y
                            i
                          
                          )
                        
                      
                      ⁢
                      
                        p
                        k
                      
                      ⁢
                      
                        p
                        i
                      
                      ⁢
                      
                        
                          
                          
                            
                              
                                G
                                k
                              
                              ⁡
                              
                                [
                                
                                  
                                    
                                      
                                        u
                                        ki
                                      
                                    
                                  
                                  
                                    
                                      1
                                    
                                  
                                
                                ]
                              
                            
                            -
                            
                              
                                G
                                j
                              
                              ⁡
                              
                                [
                                
                                  
                                    
                                      
                                        u
                                        ji
                                      
                                    
                                  
                                  
                                    
                                      1
                                    
                                  
                                
                                ]
                              
                            
                          
                          
                        
                        F
                        2
                      
                    
                  
                
              
            
          
        
      
    
    
      
        (
        6
        )
      
    
  




transforms of those two charts maps the point to the same place in coordinate space.

This can be set up as a weighted least-squares problem solved in closed form. Equation (6) generates a homogeneous set of equations that uniquely determines a solution up to an affine transform of G.
For the moment, pick one neighborhood and locate it at the center, thereby eliminating the affine indeterminacy. This adds the constraint

    
    
        G1=(I, 0)T.
    
    


To solve, define an indicator matrix Fk{dot over (=)}(0, . . . , 0, I, 0, . . . , 0)T with the identity matrix occupying the kth block, such that Gk=GFk.
Let Pk{dot over (=)}diag([pyjk(y1), . . . , pyjk(yN)]pk) contain the per-chart probability of all points on its diagonal. The squared error is then 



  
    
      
        
          
            
              
                ∑
                k
              
              ⁢
              
                [
                
                  
                    Φ
                    ⁡
                    
                      (
                      
                        
                          G
                          ⁢
                          
                            
                              U
                              _
                            
                            k
                          
                        
                        -
                        
                          [
                          
                            
                              
                                
                                  U
                                  1
                                
                              
                            
                            
                              
                                0
                              
                            
                          
                          ]
                        
                      
                      )
                    
                  
                  ⁢
                  
                    P
                    k
                  
                  ⁢
                  
                    P
                    1
                  
                
                )
              
            
            +
            
              
                ∑
                
                  j
                  ≠
                  k
                
              
              ⁢
              
                
                  Φ
                  ⁡
                  
                    (
                    
                      
                        G
                        ⁢
                        
                          
                            U
                            _
                          
                          j
                        
                      
                      -
                      
                        G
                        ⁢
                        
                          
                            U
                            _
                          
                          k
                        
                      
                    
                    )
                  
                
                ⁢
                
                  P
                  j
                
                ⁢
                
                  P
                  k
                
                ⁢
                
                  )
                
              
            
          
          ]
        
        ,
      
    
    
      
        (
        7
        )
      
    
  




where 



  
    
      U
      _
    
    k
  
  ⁢
  
    =
    .
  
  ⁢
  
    
      F
      k
    
    ⁡
    
      [
      
        
          
            
              U
              k
            
          
        
        
          
            1
          
        
      
      ]
    
  




and Φ(X){dot over (=)}∥X∥F2=trace(XTX). The first Φ(·) term penalizes

inconsistency with the anchor chart and the Φ(.) second term penalizes all other pair-wise patch inconsistencies. Differentiating with respect to G, 



  
    
      
        
          
            
              
                d
                
                  d
                  G
                
              
              =
              
              ⁢
              
                
                  ∑
                  k
                
                ⁢
                
                  [
                  
                    
                      2
                      ⁢
                      
                        
                          U
                          _
                        
                        k
                      
                      ⁢
                      
                        P
                        k
                        2
                      
                      ⁢
                      
                        P
                        1
                        2
                      
                      ⁢
                      
                        
                          U
                          _
                        
                        k
                        T
                      
                      ⁢
                      
                        G
                        T
                      
                    
                    -
                    
                      2
                      ⁢
                      
                        
                          U
                          _
                        
                        k
                      
                      ⁢
                      
                        P
                        k
                        2
                      
                      ⁢
                      
                        
                          
                            P
                            1
                            2
                          
                          ⁡
                          
                            [
                            
                              
                                
                                  
                                    U
                                    1
                                  
                                
                              
                              
                                
                                  0
                                
                              
                            
                            ]
                          
                        
                        T
                      
                    
                    +
                  
                
              
            
          
        
        
          
            
              
              ⁢
              
                
                  ∑
                  
                    j
                    ≠
                    k
                  
                
                ⁢
                
                  
                    (
                    
                      
                        
                          U
                          _
                        
                        k
                      
                      -
                      
                        
                          U
                          _
                        
                        j
                      
                    
                    )
                  
                  ⁢
                  
                    P
                    k
                    2
                  
                  ⁢
                  
                    
                      
                        P
                        j
                        2
                      
                      ⁡
                      
                        (
                        
                          
                            
                              U
                              _
                            
                            k
                          
                          -
                          
                            
                              U
                              _
                            
                            j
                          
                        
                        )
                      
                    
                    T
                  
                  ⁢
                  
                    G
                    T
                  
                
              
              ]
            
          
        
      
    
    
      
        (
        8
        )
      
    
  
  
    
      
        
          G
          T
        
        =
        
          
            
              (
              
                
                  
                    ∑
                    k
                  
                  ⁢
                  
                    
                      
                        U
                        _
                      
                      k
                    
                    ⁢
                    
                      
                        P
                        k
                        2
                      
                      (
                      
                        
                          ∑
                          
                            j
                            ≠
                            k
                          
                        
                        ⁢
                        
                          P
                          j
                          2
                        
                      
                      )
                    
                    ⁢
                    
                      
                        U
                        _
                      
                      k
                      T
                    
                  
                
                -
                
                  
                    ∑
                    
                      j
                      ≠
                      k
                    
                  
                  ⁢
                  
                    
                      
                        U
                        _
                      
                      k
                    
                    ⁢
                    
                      P
                      k
                      2
                    
                    ⁢
                    
                      P
                      j
                      2
                    
                    ⁢
                    
                      
                        U
                        _
                      
                      j
                      T
                    
                  
                
              
              )
            
            
              -
              1
            
          
          ⁢
          
            (
            
              
                ∑
                k
              
              ⁢
              
                
                  
                    U
                    _
                  
                  k
                
                ⁢
                
                  P
                  k
                  2
                
                ⁢
                
                  
                    
                      P
                      1
                      2
                    
                    ⁡
                    
                      [
                      
                        
                          
                            
                              U
                              1
                            
                          
                        
                        
                          
                            0
                          
                        
                      
                      ]
                    
                  
                  T
                
              
            
            )
          
        
      
    
    
      
        (
        9
        )
      
    
  




Now, remove the dependence on an arbitrary reference neighborhood 



  
    
      
        
          
            
              G
              =
              
              ⁢
              
                arg
                ⁢
                
                   
                
                ⁢
                
                  
                    min
                    G
                  
                  ⁢
                  
                    
                      ∑
                      
                        j
                        ≠
                        k
                      
                    
                    ⁢
                    
                      
                        Φ
                        ⁡
                        
                          (
                          
                            
                              G
                              ⁢
                              
                                
                                  U
                                  _
                                
                                j
                              
                            
                            -
                            
                              G
                              ⁢
                              
                                
                                  U
                                  _
                                
                                k
                              
                            
                          
                          )
                        
                      
                      ⁢
                      
                        P
                        j
                      
                      ⁢
                      
                        P
                        k
                      
                      ⁢
                      
                        )
                      
                    
                  
                
              
            
          
        
        
          
            
              =
              
              ⁢
              
                arg
                ⁢
                
                   
                
                ⁢
                
                  
                    min
                    G
                  
                  ⁢
                  
                    Φ
                    ⁢
                    
                       
                    
                    ⁢
                    
                      (
                      GQ
                      )
                    
                  
                
              
            
          
        
        
          
            
              =
              
              ⁢
              
                arg
                ⁢
                
                   
                
                ⁢
                
                  
                    min
                    G
                  
                  ⁢
                  
                     
                  
                  ⁢
                  
                    trace
                    ⁢
                    
                       
                    
                    ⁢
                    
                      GQQ
                      T
                    
                    ⁢
                    
                      G
                      T
                    
                  
                
              
            
          
        
      
    
    
      
        (
        10
        )
      
    
  




G1 by combining the two Φ(.) terms in equation 7, where Q{dot over (=)}Σt≠jUjUkPjPk). If the requirements is that GGT=I to prevent degenerate solutions, then equation 10 is solved, up to an arbitrary rotation in coordinate space, by setting GT to the eigenvectors associated with the smallest eigenvalues of QQT.

The eigenvectors can be determined efficiently without explicitly forming QQT. Other numerical efficiencies derive from the fact that most of the probabilities in each Pk are vanishingly small.
After the connection is solved, equation (5) gives the forward projection ƒ 141 of any point y into coordinate space. There are several numerically distinct candidates for the backprojection ƒ1 142, posterior mean, mode, or exact inverse. Noting that chart-wise projection defines a complementary density in global chart space 



  
    
      
        p
        
          x
          |
          k
        
      
      ⁡
      
        (
        x
        )
      
    
    =
    
      N
      ⁡
      
        (
        
          
            x
            ;
            
              
                G
                k
              
              ⁡
              
                [
                
                  
                    
                      0
                    
                  
                  
                    
                      1
                    
                  
                
                ]
              
            
          
          ,
          
            
              
                G
                k
              
              ⁡
              
                [
                
                  
                    
                      
                        
                          [
                          
                            
                              I
                              d
                            
                            ,
                            0
                          
                          ]
                        
                        ⁢
                        
                           
                        
                        ⁢
                        
                          
                            
                              Λ
                              k
                            
                            ⁢
                            
                               
                            
                            [
                            
                              
                                I
                                d
                              
                              ,
                              0
                            
                            ]
                          
                          T
                        
                      
                    
                    
                      0
                    
                  
                  
                    
                      0
                    
                    
                      0
                    
                  
                
                ]
              
            
            ⁢
            
               
            
            ⁢
            
              G
              k
              T
            
          
        
        )
      
    
  
  ,




and the back-projection of x into neighborhood k on the surface of the manifold, p(y|x,k), is a Dirac delta whose mean is a linear function of x.
Under this assumptions, the posterior mean map from low-dimensional coordinate space to high-dimensional data space is obtained by integrating out uncertainty over which chart generates x: 



  
    
      
        
          
            y
            |
            x
          
          ⋒
        
        =
        
          
            
              ∑
              k
            
            
               
            
          
          ⁢
          
            
              
                p
                
                  x
                  ⁢
                  
                    
                    k
                  
                
              
              ⁡
              
                (
                x
                )
              
            
            ⁢
            
              
                
                  p
                  k
                
                ⁡
                
                  (
                  
                    
                      μ
                      k
                    
                    +
                    
                      
                        
                          
                            W
                            k
                            T
                          
                          ⁡
                          
                            (
                            
                              
                                G
                                k
                              
                              ⁡
                              
                                [
                                
                                  
                                    
                                      I
                                    
                                  
                                  
                                    
                                      0
                                    
                                  
                                
                                ]
                              
                            
                            )
                          
                        
                        +
                      
                      ⁢
                      
                        (
                        
                          x
                          -
                          
                            
                              G
                              k
                            
                            ⁡
                            
                              [
                              
                                
                                  
                                    0
                                  
                                
                                
                                  
                                    1
                                  
                                
                              
                              ]
                            
                          
                        
                        )
                      
                    
                  
                  )
                
              
              .
            
          
        
      
    
    
      
        (
        11
        )
      
    
  




where X+ denotes pseudo-inverse. In general, this back-projecting function does not reconstruct the original points. Instead, it generates a surface that passes through the weighted average of the μi of all the neighborhoods in which yi has non-zero probability, much like a principal curve passes through the center of each local group of points.
EFFECT OF THE INVENTION
The invention breaks the problem of kernel-based dimensionality reduction into three steps: first, estimate a locally linear scale and intrinsic dimensionality, second find a set of data point-centered on locally linear neighborhoods or “charts” such that adjoining charts span maximally similar subspaces, and third determine a minimal distortion merger or “connection” of all charts.
Steps two and three have Bayesian settings, but by offloading the non-linearity onto the kernels, these steps become least-squares problems, solvable in closed form. Charting according to the invention performs well on data samples where the signal-to-noise ratio confounds state-of-the-art methods.
There are two possible reasons for this. First, some non-local information is used to construct the system of charts and their connection, and second the mapping only preserves the component of local point-to-point relationships that project onto the manifold; relationships perpendicular to the manifold are discarded. The first reason can be understood as a smoothness bias that makes the charting more robust to noise, and the second reason has a simple and direct analogy to noise-suppression by principle component analysis.
Although the invention has been described by way of examples of preferred embodiments, it is to be understood that various other adaptations and modifications may be made within the spirit and scope of the invention. Therefore, it is the object of the appended claims to cover all such variations and modifications as come within the true spirit and scope of the invention.

## Claims

1. A method for determining mappings between high-dimensional measured samples and a reduced-dimensional manifold, comprising:
acquiring samples of a physical system, the samples having a high number of dimensions; 
determining a low number of dimensions of a manifold embedded in a space of the high-dimensional samples; 
fitting local charts having the dimensions of the low number of dimensions of the manifold to selected high-dimensional samples, wherein each chart is represented by an elliptically shaped sphere with a set of dominant axes lying in a plane of the manifold, and a set of minor axis perpendicular to the manifold; and 
connecting the charts to determine a forward mapping from any high-dimensional sample to a coordinate on the manifold and a reverse mapping from any coordinate on the manifold to a corresponding point in high-dimensional sample space.

2. The method of claim 1 wherein the number of dimensions of the manifold is substantially less than the number of dimensions of the samples.

3. The method of claim 1 further comprising:
estimating the number of dimensions of the manifold from the samples.

4. The method of claim 1 wherein the mappings between the sample space and the coordinate space are estimated from a Gaussian mixture model.

5. The method of claim 4 wherein the Gaussian mixture model is fitted to maximize an overlap between neighboring mixture components by minimizing cross-entropies.

6. The method of claim 4 wherein the Gaussian mixture model is fitted near each sample.

7. The method of claim 1 further comprising:
acquiring the samples from an image of a face.

8. The method of claim 1 further comprising:
compressing the samples by the forward mapping.

9. The method of claim 1 further comprising:
synthesizing new samples by the forward mapping.

10. The method of claim 1 further comprising:
interpolating between the samples using by the forward mapping.

11. The method of claim 1 further comprising:
removing noise from the samples using the forward and reverse mappings.

12. The method of claim 1 further comprising:
rendering the samples using the forward mapping.

13. The method of claim 1 further comprising:
analyzing the samples using the forward and reverse mappings.

14. The method of claim 7 further comprising:
synthesizing an arbitrary facial expressions using the forward and reverse mappings.

15. The method of claim 1 wherein the number of dimensions of the samples is D and the number of dimensions of the manifold is d, and 106<<D and d ̃102 for a reduction of approximately four orders of magnitude.

16. The method of claim 1 wherein the manifold is smooth, and continuous.