Heavy-Tailed Symmetric Stochastic Neigbor Embedding

by Zhirong Yang, Irwin King, Zenglin Xu and Erkki Oja

t-SNE is a very useful visualization algorithm, which inherits the idea from SNE but modifies the neighborhood probabilities to t-distributions instead of the original. This heavy-tailed distribution works pretty well on many data. This paper discusses a more general case.

After transforming the original problem into its Lagrange dual, the authors get a unified algorithm (fixed-point iteration) for a family of heavy tailed distribution (including t-distribution). This is somewhat not as interesting as I have expected.

The authors also discussed how to integrate supervised information into the SNE algorithm. But their idea is kind of rough: inserting a similarity computed from the supervised information into the similarity in the high-dimensional space. I don't know whether this would truly work...

K-means Clustering via Principal Component Analysis

by Cris Ding and Xiaofeng He

Th most important result of this paper is that principal components are the continuous solutions to the discrete cluster membership indicators for K-means clustering and the subspace spanned by the cluster centroids are given by the spectral expansion of the data covariance matrix truncated at K-1 terms.

I'd like to say their discovery is very similar to that of normalized cut and the spectral expansion of Laplacian matrix. Let

J_K = \sum_{k = 1}^K \sum_{i \in C_k} \|x_i - m_k \|^2

be the objective of K-means clustering,

d(C_k, C_l) = \sum_{i \in C_k} \sum_{j \in C_l} \| x_i - x_j \|^2

be the distance between two clusters. It can be shown that

J_K = N \bar{y}^2 - \frac{1}{2} J_D,

where in 2-way clustering case

J_D = \frac{N_1 N_2}{N} \left( 2 \frac{d(C_1, C_2)}{N_1 N_2} - \frac{d(C_1, C_1)}{N_1^2} - \frac{d(C_2, C_2)}{N_2^2}\right).

Further we have

\frac{d(C_1, C_2)}{N_1 N_2} = \frac{d(C_1, C_1)}{N_1^2} + \frac{d(C_2, C_2)}{N_2^2} + (m_1 - m_2)^2.

The important theorem they proved says in 2-way clustering, the continuous solution is given by the leading principal components, those positive in one cluster and negative in the other. The optimal K-means solution in this case can be bounded with the leading spectrum,

n \bar{y}^2 - \lamnda_1 \leq J_{2} \leq n \bar{y}^2.

If we create a fully connected graph with each weight \frac{1}{N}, the linearization of the N-cut algorithm is PCA and the N-cut is simply K-means, is it right? The difference is PCA is done in the feature space and spectral clustering employs the embedding space Laplacian eigenmaps learns.

For K-way clustering, we have similar results. But unlike 2-way clustering, we can not interpret the PC in a similar way. Maybe the matting paper has better interpretation. Their second result is about the centroids.

This result has several implications. Spectral clustering is done in the embedding space learned by Laplacian eigenmap; PCA/K-means is done in the feature space; KPCA/Kernel K-means is done in the kernel-induced space. So if we considering the unsupervised KDR idea, can we find a similar interpretation such as clustering?

Visualizaing Data using t-SNE

by Laurens van der Maaten and Geoffrey Hinton

This paper proposes a variant of SNE. The SNE was previously proposed by Hinton and Roweis. Its idea is to find a projection such that the neighboring probability is preserved: given the coordinates x_i, the neighboring probability is

p_{i \mid j} = \frac{g(x_i, x_j)}{\sum_k g_j(x_j, x_k)}.

In the projection, if we define a similar probability q_{i \mid j}, we wish to have a most similar distribution of these samples, therefore we intend to minimize

\sum_j \sum_i p_{i \mid j} \log \frac{p_{i \mid j}}{q_{i \mid j}}.

The original setting is taking g_j(x_i, x_j) with Gaussian kernel (with different variances). By setting a so-called perplexity, we then binary-search a proper variance for each sample. This perplexity corresponds to the continuous version of effective neighbors. And we compute the gradient of y_i:

\frac{\partial \mathrm{KL}}{\partial y_j} = 2 \sum_i ( p_{i \mid j} - q_{i \mid j} + p_{j \mid i} - q_{j \mid i})(y_j - y_i).

The gradient decent is suppllemented with a momentum.

The so-called t-SNE has two major modification to the SNE model. For one thing, it uses a symmetric p_{i, j} instead, which is simply (p_{i \md j} + p_{j \mid i}) / 2, this will simplify the derivatives of the objective function,

\frac{\partial \mathrm{KL}}{\partial y_j} = 4 \sum_i (p_{i, j} - q_{i, j})(y_j - y_i).

For another thing, we use a different pdf for the projection instead of the Gaussian kernel. The "t" comes from here, since they propose using t-distribution with 1 degree of freedom (therefore a Cauchy distribution). Then the derivatives become

\frac{\partial \mathrm{KL}}{\partial y_j} = 4 \sum_i (p_{i, j} - q_{i,j})(y_j - y_i)(1 + \| y_j - y_i\|^2)^{-1}.

The great thing about t-SNE is its capability to yield a good embedding. But somehow it is a little difficult to increase the embedding dimension. Not sure about it though, doing some experiments with it.

High-dimensional Data Analysis: The Curses and Blessings of Dimensionality

by David L. Donoho

In this non-technical paper, Donoho prospected the development of data analysis and suggested the difficulty in dealing with high dimensional data (the curse of dimensionality) and the possible advantages we may take (the blessing of dimensionality). His advisor John Tuckey coined the word data analysis, which differs from the traditional mathematics or statistics in the way of rigorous derivation: more heuristic ideas are applied. He argued that nowadays people choose quite a different method of exploration in science: only a few variables are selected to make a model which can be analytically formulated and studied in the traditional way while now it is more widely employed to collect a large amount of data or the information about one sample is enriched via all kinds of ways.

Curse of Dimensionality:

• For a naive search algorithm in a high dimensional discrete space, the number of candidates is exponential of the dimensionality.
• For optimization, we need (1/\epsilon)^d evaluations for a Lipschitz continuous objective function for an approximation of error less than \epsilon.
• In function approximation to a Lipschitz continous function, we also need O\big((1/\epsilon)^d \big) evaluations for a uniform approximation of error \epsilon.
• In numeric integration of a Lipschitz continous function, again we need O\big( (1/\epsilon)^d\big) evaluations on a grid in order to obtain an approximate integration of error \epsilon.
• Nonparametric extimation will be difficult. The much slower convegence rate is the ugly head of the curse of dimensionality.
• Regularization is usually employed to tackle the curse of dimensionality.
  

Blessing of Dimensionality:

• The concentration of measure phenomenon is for a lot of distribution in the high dimensional space. Let (\Omega, \mathcal{F}, \mu) be a metric measurable space, d(\cdot, \cdot) is the distance and \mu is the probability measure of the \sigma-field \mathcal{F}. Let \alpha(\epsilon) = \sup_A \{ \mu( \Omega \backslash A_\epsilon) \mid \mu(A) = 1/2\} where A_\epsilon = \{ x \in \Omega \mid d(x, A) < \epsilon\} is the \epsilon-extension of the measurable set A. The so-called Levy family (\Omega_n, \mathcal{F}_n, \mu_n), which satisfies \forall \epsilon > 0, \lim_{n \to \infty} \alpha_n(\epsilon) = 0, means that the probability measure will concentrate in the center part of the support when n increase. If the convergence to 0 is C \exp(-c n \epsilon^2), then it is also called normal Levy family. In many cases, the index n can be the dimensionality, which means as the dimensionality goes up, we may observe the concentration of measure. For a normal Levy family, the tail's behavior is like a Gaussian. The concentration of measure means in high dimensional space it is more likely that the data you sample contains much information as you don't need to go further from the data, you have already covered the support. The probability decays as you leave a set with a fixed nonzero measure set.
• Dimension asymptotics. This is the refinement of the previous argument.
• Approach to continuum: there is an underlying compactness to the space of observed data which will be reflected by approximate finite-dimensionality and an increasing simplicity of analysis for large d. (sounds like the reason for compressive sensing).
• Approximation through randomness. Random projections, CS are in this range.
• Approximation through computational harmonic analysis, which is a way to the continuum.
  

I am not familiar with harmonic analysis. Many points in this paper are still vague to me. I must reread it when I gain more knowledge.

Random Projection for Manifold Learning

by Chinmay Hedge and Michael B. Wakin and Richard G. Baraniuk

This paper tells us that we may project the data from a higher dimensional space (n dimensions) into a random subspace (m dimensions) and we will not lose much information of the intrinsic structural information. Therefore we may apply one of the manifold learning algorithms to the randomly-projected data and still recover their structure in a proper space (d dimensions). And d < m << n and m = O(d log n).

The first theorem comes from RIP satisfaction, when

m \geq O\left( \frac{d \log (n V \tau^{-1}) \log(\rho^{-1})}{\epsilon^2}\right)

with probability no less than 1-ρ, for any x, z, we have

(1 - \epsilon) \sqrt{\frac{m}{n}} \leq \frac{d_i( \Phi x, \Phi z)}{d_i( x, z )} \leq (1 - \epsilon) \sqrt{\frac{m}{n}},

where Φ is the random orthogonal projection, V is the volume and 1/τ is the condition number of the latent manifold. To estimate the intrinsic dimensions, the author propose to employ Grassberger-Procaccia algorithm (based on correlation). They found bounds for the estimated dimensions: when

m \geq O\left( \frac{d \log( n V \tau^{-1}) \log (\rho^{-1})}{\beta^2 \delta^2} \right),

then the estmated correlation dimension satisfies

(1 - \delta) \hat{d} \leq \hat{d}_\Phi \leq (1 + \delta) \hat{d},

with probability exceeding 1-ρ. And the residual variance of the ISOMAP embedding is also bounded

R_\Phi \leq R + C \Gamma^2 \epsilon,

where Γ is the diameter of the point cloud.

Since in practice we are not sure about the intrinsic dimensionality, we usually try to add more orthogonal vectors and see whether the residual variance is small enough.

Well it is atypical way of applying CS in machine learning. Put all data in a random projection subspace and prove there are bounds for later learning algorithms. It is not the way I'd like to use CS.

Sufficient Dimensionality Reduction

by Amir Globerson and Naftali Tishby

The authors proposed a nonparametric way for sufficient dimensionality reduction. Their idea is very similar to CME. We don't know Pr( x, y ) but we may estimated Pr(x), Pr(y) and several measurement φ(x) (the reduced features). Therefore, given a φ(x), the minimum mutual information I(q) of q(x, y) (the estimated joint density) can be seen as the goodness of φ(x) and the maximization of this minimum leads to the best φ(x) available to us.

The minimum mutual information I(q) can be obtained with variational methods. We are actually finding the best q(x, y) such that q(x) = Pr(x), q(y) = Pr(y) and the expectation of φ(x) on the conditional q(x | y) and Pr(x | y) is the same. The solution has a CME flavor (the solution is a so-called I-projection of q(x, y) on all the distributions that satisfy the constraints): convex and unique. Let the corresponding Lagrange multipliers for φ(x) be ψ(y). Then in a sense, the functional I(q) can be regarded as the Lagrange function with the parametrization of q(x, y).

The maximization step can be converted into a minimization step, which can still be solved by I-projection. This is because I(q) = I(Pr(x, y)) - KL( Pr(x, y) | q ) where q(x) in the ``dual feasible set.'' And therefore we need two I-projections.

The algorithm then works in this way: fix the φ(x) part, project the conditional q(y | x) onto the subspace for tunable ψ(y) spaces; update the joint density with q(y | x) Pr(x); fix ψ(y) part, project the conditional q(x | y) onto the subspace for tunable φ(x) spaces; update the joint density with q(x | y) Pr(y).

I am not quite sure about my own interpretation though but the convergence (to a local minimum) of the algorithm is not difficult to prove. The I-projections can be done partially with a few GIS iterations. The most important idea in this paper is from information geometry: the Pythagoras theorem can be interpreted as a geodesic distance on the manifold and the dual interpretation of φ(x) and ψ(y) is very intuitive.

A related stat article:
Francesca Chiaromnte and R. Dennis Cook: Sufficient dimension reduction and graphics in regression, AISM.

Dimensionality Reduction for Supervised Learning with Reproducing Kernel Hilbert Spaces

by Kenji Fukumizu, Francis R. Bach and Michael I. Jordan

This paper tells us one way of SDR (sufficient dimensionality reduction) using KDR (kernel dimensionality reduction), though it has a somewhat different path to reach the conclusion. I will scan the earlier nonparametric version.

This model proposed by the author is kind of semi-parametric. The sufficiency of a dimensionality reduction algorithm for a regression problem is that y is independent of x conditioned on the projection of x in a subspace (let the projection be Px = u and v is the residual). Therefore either we maximize the mutual information of y and u or we minimize the mutual information of y and v. Consider about the KDR technique: the mutual information can be approximated with Kernel CCA or KGV (kernel generalized variance).

The different derivation starts with two definitions: one for covariance operator and the other conditional covariance operator. In a way, these operators formulates the idea of ``two r.v.s are independent when arbitrary functions of them are independent,'' but in a more abstract and rigorous way. With these concept we know that the conditional covariance operator Σ_YY|U is no less than Σ_YY|X and the equality happens when Y is independent of V. And Σ_YX|U = 0 when Y is independent of V.

When we concretize these operators with matrices, we have the same formulation as in KDR. Therefore the optimization is done in the same way. This procedure is analogous to feature selection and may be applied to variable selection.

Localized Sliced Inverse Regression

by Qiang Wu, Sanyan Mukherjee and Feng Liang

Sliced inverse regression (SIR) for classification is equivalent to FDC. That is to solve a generalized eigenvalue problem Γβ=λΣβ, where Γ is the between-class covariance matrix and Σ is the total covariance matrix. The largest eigenvalues of this problem correspond to the disired directions. Well, in some cases we use the in-class covariance matrix B for Σ.

This has been discussed in Jieping Ye's paper extensively. In a way Σ is more likely to be non-singular than B. Due to the fact Σ = Γ + B, there exists a non-singular matrix A such that they can be diagonalized simultaneously with a transform Az = x. If they are non-singular, it matters not whether we use B or Σ. When B or Σ is singular, we might use its principal subspace or its null space for approximation (well it's not approximation I guess). I wonder what we shall do for SIR. Maybe we do the same rubbish too.

Now we come to the LSIR. The authors proposes a localized Γ by choosing the kNN samples of the same class. The only difference between LSIR and LFDC, I think, is the denominator part (i.e. a global covariance class v.s. a localized in-class covariance). Maybe I have to explore the corresponding graph to see their difference more clearly.

Their extension for SSL is to put all unlabelled samples into to all slices. That might be natural in the graph embedding framework. Why should I read this paper then?

Probabilistic Non-linear Principal Component Analysis with Gaussian Process Latent Variable Models

by Neil Lawrence

This might be the first article that covers most of GP-LVM, which is of great interest to me since last time Fei Sha mentioned this model related to a former research of mine. The idea is not very new but to me it is something fresh.

The basic model comes from Probabilistic PCA (PPCA). The PPCA model is a latent variable model that interpret PCA in a probabilistic framework. The latent variable z generate our observations with a linear transformation x = Wz + b + e, where e is an isotropis noise. The bias b could be estimated with sample mean independently. In PPCA, z has a prior of standard Gaussian. The estimation of W and the variance for the noise could be achieved with EM algorithm. So you might understand, the EM algorithm actually maximizes the marginal distibution of x and we marginalize z.

In a linear model, it's easy to get the reduction function and the reconstruction function easily. But as for nonlinear model, usually we can only have one-way mapping. The inverse mapping is usually difficult to get. So does GP-LVM. We might think a GP-LVM like this. The first part comes if we marginalize W instead of z in the PPCA model. We endow W with a similar Gaussian prior and z could be obtained with eigen decomposition of X X^T. Actually the solution resembles PCA and classic MDS. For a GP-LVM, unlike kernel PCA putting a kernelized Gram matrix for observed data (so the solution is still obtained via eigen decomposition), we have to find the best configurationsfor all z, whose kernelized matrix best approximated X X^T. Therefore, although we have a quite similar optimization problem, it is much harder to solve now due to the nonlinearity of the kernel function. Therefore a gradient based search such as CG must be applied.

But as we might foresee, the computational cost is intolerable for practical problems. The paper then focuses on designing an efficient algorithm for GP-LVM. The most important technique they employ is informative vector machine, which selects a comparatively small active set for the kernel methods (as Nystrom method).

Regression on Manifold Using Kernel Dimension Reduction

by Jens Nilsson, Fei Sha and Michael I. Jordan

This paper describes a regression problem: the samples come from a manifold embedded in a high dimensional space and the response is function defined on the manifold. There are 3 techniques applied in this problem: sufficient dimension reduction (SDR), kernel dimension reduction (KDR) and manifold learning.

The basic idea is from KDR, in which we seek a subspace (to be the central subspace) B such that the projection onto the subspace are sufficient to regress the responses. This idea (in linear case, namely SDR) requires solving the following optimization problem,

where K_Z^c is projected Gram matrix. Here with manifold learning, the projected sample Bx_i in the RKHS is approximated with the affine transformation of the coordinates obtained by the manifold learning algorithm (e.g. Laplacian Eigenmap). Therefore the optimization problem could be formulated as

After solving the problem, the regression could be done in the corresponding subspace. The difficult part is how to solve this nonlinear and nonconvex problem. It is solved with projected gradient descent method.

But it seems hard to extend to new data. It is argued that by computing the latent B, it is possible to extend the regressor on new data. But it doesn't say how. I guess it must be much tougher :-(

Trace Ratio vc. Ratio Trace for Dimensionality Reduction

by Huan Wang, Shuicheng Yan, Dong Xu, Xiaoou Tang, Thomas Huang

This is a CVPR paper in 2007. I read Yan's CVPR paper in 2005. So I am quite familiar with this direction. Basically I agree with their argument on using trace ratio, esp. for those graph embedding algorithms. It is more sounding by directly maximizing the ratio of between-class scatterness to with-in class compactness. the optimization becomes difficult when trace ratio is the objective instead of the originally used formulations, e.g. iteratively defined optimization, trace formulation and determinant formulation.

The constraints of their proposed optimization ensure orthogonality of the projection. The algorithm, also an approximation in the principal subspace of S_t, iteratively updates the projection V, such that

V^* = arg max tr V^T (S^u_p - λⁿ S^u_t ) V s.t. V^TV = I

The iteration will increase the objective function monotonously and ensures a global optimum (main contribution). For other constraints, such as one in kernelization, an eigen decomposition is first applied.

Hmm... another work as an extension to this (vector -> tensor) says a global optimum might not be obtainable, which I might check later.

Uncovering Shared Structures in Multiclass Classification

by Yonatan Amit, Michael Fink, Nathan Srebro and Shimon Ullman

In this paper, we have Srebro, one of the authors of previous papers on MMMF. Now we are finding linear projection, a matrix W. To get a good projection, they define a loss function:

which is regularized by the norm of W. Thus we have,

But, they argue the Frobenius norm can be substituted for the trace norm, resulted in the following problem which is formally the same as MMMF:

After analyzing the dual problem of the equation above, it is found out the problem can be kernelized as

Then to solve these problem, SDP or PR-GC can be used.

Null Space versus Orthogonal Linear Discriminant Analysis

by Jieping Ye and Tao Xiong

This article tells me one important thing. LDA or Fisher discriminant criterion is still being studied by lots of guy. Last time I read one by Manli Zhu, in which they propose another different criterion than the traditional ones when the co-variance matrix is rank-deficient.

One of the simplest idea is applying a ridge purturbation to the singular with-in class covariance matrix, which is usu. called regularization. The other is finding an approximation in the principal subspace of S_w. Zhu's idea roots from the latter one, in which all principal directions are selected for later general eigenvalue problem. Zhu rejects those principal directions that are almost perpendicular to the subspace spanned by S_b.

In this paper, the author mentions several other strategies for dealing with the singularity of S_w, two of which they put attention to Orthognal LDA (OLDA) and Null Space LDA (NLDA). Their main result is they yield the same result under certain condtions. The condition is mild when dealing with high dimensional data (the number of dimension is much higher than the number of samples).

I guess I will review this article after I check those LDAs, which personally I suspect their usability. Here are those literatures:

• A new LDA-based face recognition system which can solve the small size problem, Pattern Recognition 33 (2000) by Chen and et al.
• Solving the small sample size problem of LDA, ICPR 2002, by Huang and et al.
• Characterization of a family of algorithms for generalized discriminant analysis on undersampled problems, JMLR 2005.
• Penalized discriminant analysis, Anals of Statistics 23 (1995) by Hastie and et al.