Scaling Locally Linear Embedding

Yasuhiro Fujiwara†‡*, Naoki Marumo†, Mathieu Blondel†, Koh Takeuchi†,
Hideaki Kim†, Tomoharu Iwata†, Naonori Ueda††NTT Communication Science Laboratories,‡NTT Software Innovation Center,‡*Osaka University

Introduction• Locally Linear Embedding (LLE) is one of the
popular approaches of nonlinear dimensionality
reduction• However, it requires high computational cost

• The goal in this study is to enhance the processing speed of LLE• We propose a novel approach, Ripple, as asolution2

LLE: overview• LLE is an approach that can effectively reducethe dimensionality of multi-dimensional data*1– Effectively represents manifold (clustering structure)of a dataset by embedding• 3 steps:
1. k-NN graph
2. Edge weight by regression
3. Eigen decomposition based on adjacency matrix• Eigen vectors correspond to embedding*1 Roweis et al., “Nonlinear Dimensionality Reduction by Locally Linear Embedding”, Science, 2000
3

LLE: step1 k-NN graph• Use naive approach– Pick up data points (nodes) one by one– Compute top 𝐾𝐾 based on Euclidean distance– Computation cost: 𝑂𝑂 𝑁𝑁 2 𝑀𝑀
• 𝑁𝑁: # data points, 𝑀𝑀: # dimensions

4LLE: step2 edge weight• Compute edge weights by regression (sum ofedge weights is normalized to 1)• Lagrange multiplier method– Regression error of node 𝑥𝑥𝑝𝑝
𝜀𝜀 = 𝑥𝑥𝑝𝑝 − �w[𝑝𝑝, 𝑖𝑖]𝑥𝑥𝑖𝑖

𝑥𝑥𝑖𝑖 ∈ℕ[𝑥𝑥𝑝𝑝 ]2regression=�

𝑥𝑥i ,𝑥𝑥j ∈ℕ[𝑥𝑥𝑝𝑝 ]𝐺𝐺[i, j]ℕ[𝑥𝑥𝑝𝑝 ]: top K data points of 𝑥𝑥𝑝𝑝
G[i, j]: (𝑥𝑥𝑝𝑝 -𝑥𝑥i )(𝑥𝑥𝑝𝑝 -𝑥𝑥𝑗𝑗 )– Edge weights are obtained from inverse of matrix 𝐺𝐺
w p, i = �𝐾𝐾𝐺𝐺

𝑗𝑗=1−1[𝑖𝑖, 𝑗𝑗]��

𝐾𝐾𝑖𝑖=1

�𝐾𝐾

𝐺𝐺 −1 [𝑖𝑖, 𝑗𝑗]

𝑗𝑗=1– Computation cost: 𝑂𝑂 𝑁𝑁(𝑀𝑀𝐾𝐾 2 + 𝐾𝐾 3 )• 𝑂𝑂(𝑀𝑀𝐾𝐾 2 ) for computing matrix G• 𝑂𝑂(𝐾𝐾 3 ) for computing matrix 𝐺𝐺 −1w[p, i]: edge weight from 𝑥𝑥𝑖𝑖 to 𝑥𝑥𝑝𝑝

5

LLE: step 3 Eigen decomposition• Minimize error of dimensionality reduction– Compute Eigen decomposition of matrix K obtained
from adjacency matrix 𝑊𝑊 of step 2
K = 𝐼𝐼 − 𝑊𝑊 𝑇𝑇 (𝐼𝐼 − 𝑊𝑊)– “Bottom (the smallest) ” Eigen vectors of matrix K are
the embedding to low dimensionalities• Not top (the highest) Eigen vectors– Computation cost: 𝑂𝑂 𝑁𝑁 3
• Size of matrix K is 𝑁𝑁 × 𝑁𝑁

6

LLE: problem• Impractical for large size of dataset– It needs high computation cost
– Computation cost: 𝑂𝑂 𝑁𝑁 2 𝑀𝑀 + 𝑁𝑁 𝑀𝑀𝐾𝐾 2 + 𝐾𝐾 3 + 𝑁𝑁 3

• Especially Eigen decomposition needs high computation cost
since it requires 𝑂𝑂 𝑁𝑁 3 time7

Proposed method• We exploit 3 approaches– Incremental weight computation• Efficiently compute weights by using inner product of datapoints and the Woodbury formula– Improve lower bound of Euclidean distance• Enhance effectiveness of pruning by exploiting graphstructure• Efficiently obtain top-k nodes– LU decomposition based Eigen decomposition• Compute eigenvectors “ascending” order of Eigen values– Power method compute eigenvectors in “descending” order• We avoid computing inverse matrix by LU decomposition• Significantly reduce computation cost of Eigen decomposition
8

Incremental weight computation• Basic idea
– Compute k-NN graph by picking up nodes one by one– Since such nodes share nearest neighbors, we usesthe Wood bury formula to compute inverse matrix– But, each node has totally different matrix 𝐺𝐺𝑝𝑝
𝐺𝐺𝑝𝑝 [i, j]: (𝑥𝑥𝑝𝑝 -𝑥𝑥i )(𝑥𝑥𝑝𝑝 -𝑥𝑥𝑗𝑗 )– We use matrix 𝐶𝐶𝑝𝑝 instead of 𝐺𝐺𝑝𝑝 that can shareelement in applying Lagrange multiplier method
Share nearest
neighbors

Processing node forweight computation
Already processed nodefor weight computation

9

Incremental weight computation• Let y = ∑𝑥𝑥𝑖𝑖 ∈ℕ[𝑥𝑥𝑝𝑝 ] 𝑊𝑊[𝑝𝑝, 𝑖𝑖]𝑥𝑥𝑖𝑖 , regression error 𝜀𝜀𝜀𝜀 = 𝑥𝑥𝑝𝑝 − �𝑥𝑥𝑖𝑖 ∈ℕ[𝑥𝑥𝑝𝑝 ]𝑊𝑊[𝑝𝑝, 𝑖𝑖]𝑥𝑥𝑖𝑖

2= 𝑥𝑥𝑝𝑝 (𝑥𝑥𝑝𝑝 )𝑇𝑇 −2𝑥𝑥𝑦𝑦 𝑇𝑇 + 𝑦𝑦(𝑦𝑦)𝑇𝑇• Let 𝐿𝐿 = 𝜀𝜀 + 𝛾𝛾𝑊𝑊 be function of Lagrange multiplier method∂𝐿𝐿𝑇𝑇𝑥𝑥𝑖𝑖 𝑥𝑥𝑗𝑗 𝑊𝑊 𝑝𝑝, 𝑖𝑖 + 𝛾𝛾= −2𝑥𝑥𝑝𝑝 𝑥𝑥𝑖𝑖 𝑇𝑇 + 2 �∂𝑊𝑊[𝑝𝑝, 𝑖𝑖]𝑥𝑥j ∈ℕ 𝑥𝑥𝑝𝑝

𝑇𝑇

∂𝐿𝐿= ∑𝑥𝑥 ∈ℕ 𝑥𝑥𝑝𝑝 𝑊𝑊 𝑝𝑝, 𝑖𝑖 -1∂𝛾𝛾i• If 𝐶𝐶𝑝𝑝 [𝑖𝑖, 𝑗𝑗] = 𝑥𝑥𝑖𝑖 𝑥𝑥𝑗𝑗 and 𝑝𝑝[𝑖𝑖] = 𝑥𝑥𝑝𝑝 𝑥𝑥𝑖𝑖 𝑇𝑇 , we have 𝐶𝐶𝑝𝑝 𝑊𝑊 = 𝑝𝑝
from Lagrange multiplier method, thus 𝑊𝑊 = 𝑝𝑝(𝐶𝐶𝑝𝑝 )−1• (𝐶𝐶𝑝𝑝 )−1 is independent from 𝑥𝑥𝑝𝑝 , so incrementally updated10

Incremental weight computation• Here we assume that the first nearest neighbor node isdifferent from node 𝑥𝑥𝑝𝑝 and 𝑥𝑥𝑞𝑞
– We also assume that norms of 𝑥𝑥𝑝𝑝 and 𝑥𝑥𝑞𝑞 are 1
1𝐶𝐶𝑞𝑞 = 𝑥𝑥1 𝑥𝑥2 𝑇𝑇
𝑥𝑥1 𝑥𝑥3 𝑇𝑇

𝑥𝑥1 𝑥𝑥2 𝑇𝑇
1𝑥𝑥2 𝑥𝑥3 𝑇𝑇

𝑥𝑥1 𝑥𝑥3 𝑇𝑇
𝑥𝑥2 𝑥𝑥3 𝑇𝑇
11𝐶𝐶𝑝𝑝 = 𝑥𝑥1 𝑥𝑥2 𝑇𝑇
𝑥𝑥1 𝑥𝑥3 𝑇𝑇

different

𝑥𝑥1 𝑥𝑥2 𝑇𝑇
1𝑥𝑥2 𝑥𝑥3 𝑇𝑇

𝑥𝑥1 𝑥𝑥3 𝑇𝑇
𝑥𝑥2 𝑥𝑥3 𝑇𝑇
1• If ∆C = 𝐶𝐶𝑝𝑝 − 𝐶𝐶q , ∆C = 𝑉𝑉 𝑇𝑇 𝐷𝐷𝐷𝐷 where 𝑉𝑉 and 𝐷𝐷 are rank-2– Since ∆C has particular form, it need 𝑂𝑂 𝐾𝐾 to compute 𝑉𝑉 𝑇𝑇 𝐷𝐷𝐷𝐷

• We can compute (𝐶𝐶𝑝𝑝 )−1 at 𝑂𝑂 𝐾𝐾 2 by the Woodbury formula– We apply each different nearest neighbor node to compute (𝐶𝐶𝑝𝑝 )−1(𝐶𝐶𝑝𝑝 )−1 = (𝐶𝐶𝑞𝑞 )−1 −(𝐶𝐶𝑞𝑞 )−1 V( 𝐹𝐹 −1 + (𝑉𝑉)𝑇𝑇 (𝐶𝐶𝑞𝑞 )−1 𝑉𝑉)−1 (𝑉𝑉)𝑇𝑇 (𝐶𝐶𝑞𝑞 )−111

Improve lower bound• Prune distance computation in computing k-NN graph byapproximating Euclidean distance E 𝑥𝑥𝑝𝑝 , 𝑥𝑥𝑞𝑞
– SVD is a popular approach for approximation• Improve lower bound by SVD E[�𝑥𝑥𝑝𝑝 , 𝑥𝑥�𝑞𝑞 ] as 𝐸𝐸[𝑥𝑥𝑝𝑝 , 𝑥𝑥𝑞𝑞 ] in thefollowing form:
𝐸𝐸 𝑥𝑥𝑝𝑝 , 𝑥𝑥𝑞𝑞 =(E[�𝑥𝑥𝑝𝑝 , 𝑥𝑥�𝑞𝑞 ])2 +(𝑢𝑢𝑟𝑟 𝑥𝑥𝑝𝑝 − 𝑢𝑢𝑟𝑟 𝑥𝑥𝑞𝑞 )2

where 𝑢𝑢𝑟𝑟 𝑥𝑥𝑝𝑝 =Distanceby SVD

Lower bound by dimensions
not used in SVD(E[𝑥𝑥𝑝𝑝 , 𝑥𝑥𝑟𝑟 ])2 − (E[�𝑥𝑥𝑝𝑝 , 𝑥𝑥�𝑞𝑞 ])2

Norm of dimensions notused in SVD

We use triangular inequality in this approach12

LU decomposition based Eigen decomposition• Power method is the most popular approach in computing
Eigen vector– But, it computes the largest not smallest Eigen values– Embedding is the smallest Eigen vector of K = 𝐼𝐼 − 𝑊𝑊 𝑇𝑇 (𝐼𝐼 − 𝑊𝑊)• Inverse power method computes the smallest Eigen value– It apply power method for the inverse matrix– Its computation cost is 𝑂𝑂 𝑁𝑁 3– Impractical for large-size of dataset• We avoid the inverse matrix by LU decomposition– We have sparse matrices after LU decomposition– We can apply this approach of large graphs13

LU decomposition based Eigen decomposition• Compute LU decomposition for 𝐼𝐼 − 𝑊𝑊 (LU = 𝐼𝐼 − 𝑊𝑊)– Thus we have K = 𝐼𝐼 − 𝑊𝑊 𝑇𝑇 𝐼𝐼 − 𝑊𝑊 = 𝑈𝑈 𝑇𝑇 𝐿𝐿𝑇𝑇 𝐿𝐿𝐿𝐿

• The smallest Eigen value 𝜆𝜆N and its Eigen vector 𝑍𝑍𝑁𝑁 canbe computed as follows similar to power method:𝜆𝜆N = (𝑎𝑎𝜏𝜏 )𝑇𝑇 𝑎𝑎𝜏𝜏 ⁄ (𝑎𝑎𝜏𝜏 )𝑇𝑇 𝑎𝑎𝜏𝜏−1𝑧𝑧N = 𝑎𝑎𝜏𝜏 ⁄𝑎𝑎𝜏𝜏

𝑇𝑇𝑇𝑇 ′′′′ ′′where 𝑎𝑎𝜏𝜏−1 = 𝑈𝑈 𝑏𝑏, 𝑏𝑏 = 𝑈𝑈 𝑏𝑏 , 𝑏𝑏 = 𝐿𝐿𝑏𝑏 , 𝑏𝑏 = 𝑈𝑈𝑎𝑎𝜏𝜏• Since vector 𝑎𝑎𝜏𝜏 is updated as 𝑎𝑎𝜏𝜏−1 = 𝑈𝑈 𝑇𝑇 𝐿𝐿𝑇𝑇 𝐿𝐿𝐿𝐿𝑎𝑎𝜏𝜏 , we cancompute the smallest Eigen value– Note we have 𝑎𝑎𝜏𝜏 = 𝐾𝐾 −1 𝑎𝑎𝜏𝜏−1 since K = 𝑈𝑈 𝑇𝑇 𝐿𝐿𝑇𝑇 𝐿𝐿𝐿𝐿14

Theoretical analysis• Ripple can efficiently obtain the sameembedding results as the original approaches

15Experiment: preliminaries• We used the following five datasets– USPS; 7291 items and 256 features– SensIT; 78,823 items and 100 features– ALOI; 108,000 items and 128 features– MSD; 515,345 items and 90 features– INRIA; 1,000,000 items and 128 features• Comparison methods– CLLE: k-means based approach*3– LLL: Nystrom method based approach*4– VN: Nystrom method based approach*5*3 Hui et al., “Clustering-based Locally Linear Embedding”, ICPR, 2008
*4 Vladymyrov et al., “Locally Linear Landmarks for Large-scale Manifold Learning”, ECML/PKDD, 2013
*5 Vladymyrov et al., “The Variational Nystr¨om Method for Large-Scale Spectral Problems”, ICML, 201616

Experiment: efficiency• Wall clock time
– Ripple is much faster than existing methods
Up to 560,330, and 260 times
faster than previous methods

Ripple is scalable
to large data
Up to 2,300
times faster17

Experiment: exactness (CLLE)• Ripple yields the same result as the original approach• CLLE has trade-off between efficiency and accuracy
against # clusters of k-means method
CLLE increases
error against #clusters

Ripple isexact

Accuracy

CLLE increases
efficiency against# clusters

But, Ripple ismore efficient

Efficiency

18

Conclusions• This study proposed an efficient approach for
Locally linear embedding (LLE)• Our approach, Ripple, (1) incrementally computeedge weights, (2) improve the lower bounds in
obtaining k-NN graph, and (3) exploits LU
decomposition in computing Eigen vectors• Experimental results show that our approach isfaster than the previous approach

Thank you for your attention