JP6721535B2 - LLE calculation device, LLE calculation method, and LLE calculation program - Google Patents

Description

The present invention relates to an LLE calculation device, an LLE calculation method, and an LLE calculation program.

LLE (Locally Linear Embedding) is one of the typical methods for reducing nonlinear dimensions. The basic idea of LLE is to approximate each data point by regression with its neighboring data points and calculate the embedding of each data point into the low dimension. LLE solves two optimization problems that minimize reconstruction cost and embedding cost. In these optimizations, generation of initial values by random numbers, setting of learning rate, etc. are not necessary, so that LLE can effectively reduce the dimension.

Sam T. Roweis and Lawrence K. Saul, Nonlinear Dimensionality Reduction by Locally Linear Embedding, Science, 290(5500):2323-2326, 2000.

However, the conventional LLE has a problem that the calculation cost is high. For example, where K is the number of edges from neighboring nodes to each data point, the conventional LLE algorithm calculates a Gram matrix of size K×K for each data point, and the Lagrange undetermined multiplier method. Is used to calculate the edge weight from the inverse matrix of the Gram matrix (for example, see Non-Patent Document 1). Since the cost of calculating the inverse matrix is the cube of the number of edges to each data point, the cost of calculating the edge weight is O(NK ³ ) when N is the number of data points.

Also, the conventional LLE algorithm calculates the eigenvalue vector of the LLE kernel in order to calculate the embedding of each data point in the lower dimension. Since the size of the LLE kernel is N×N, the cost required to calculate the eigenvector is O(N ³ ). Furthermore, a memory cost of O(N ² ) is required to hold the eigenvalue vector. Therefore, applying LLE to large-scale data is unrealistic because it requires a large amount of calculation cost and memory cost.

The LLE calculation device of the present invention includes a singular value decomposition calculation unit that calculates a singular value decomposition of a multidimensional matrix having a plurality of data points, and a data point selection unit that selects a data point to be calculated from the plurality of data points. For each of the plurality of data points, an approximate distance between the data point to be calculated and an approximate distance calculation unit that calculates based on the singular value decomposition, and the approximate distance is less than or equal to a predetermined value. A distance estimation unit that calculates an estimated value of the Euclidean distance between each of the data points and the data point of the calculation target, each of the data points for which the estimated value of the Euclidean distance is a predetermined value or less, and the calculation target A proximity determining unit that determines the Euclidean distance between the data point and the data point whose calculated Euclidean distance is less than or equal to a predetermined value as a data point near the data point to be calculated, and the calculation target. Edge weight calculator for calculating the edge weight between each data point and each of the data points in the vicinity of the data point to be calculated, and an LLE kernel for the plurality of data points based on the edge weight. And a vector calculation unit that calculates an eigenvector of

According to the present invention, the cost required for calculating LLE can be reduced.

FIG. 1 is a diagram illustrating an example of the configuration of the LLE calculation device according to the first embodiment. FIG. 2 is a diagram showing Lemma 1. FIG. 3 is a diagram showing Lemma 2. FIG. 4 is a diagram showing Lemma 3. FIG. 5 is a diagram illustrating an example of an algorithm of the LLE calculation process according to the first embodiment. FIG. 6 is a diagram showing Theorem 1, Theorem 2 and Theorem 3. FIG. 7 is a flowchart showing a processing flow of the LLE calculation device according to the first embodiment. FIG. 8 is a diagram illustrating an example of a computer that executes the LLE calculation program.

Hereinafter, embodiments of an LLE calculation device, an LLE calculation method, and an LLE calculation program according to the present application will be described in detail with reference to the drawings. The present invention is not limited to the embodiments described below.

[Conventional LLE]
First, the conventional LLE will be described. Let XεR ^N×M be a matrix of N data points of M dimensions. Let x _i =(x _i [1], x _i [2],..., x _i [M]) be the i th row vector of matrix X, then x _i corresponds to the i th data point.

LLE computes a low-dimensional dataset YεR ^N×m from a high-dimensional dataset X. Let y _i =(y _i [1], y _i [2],..., y _i [m]) be a row vector of length m, and y _i is the mapping of data point x _{i to} the low dimension. Corresponding to. The basic idea in LLE is that even if a data point has a globally nonlinear structure, each data point and its neighboring data points can be locally linearly approximated. By expressing each data point from the neighborhood, LLE can express high-dimensional data in low dimension.

The conventional LLE algorithm is composed of two optimization problem solving parts that minimize the reconstruction cost and the embedding cost. In the first part, for each data point, the neighboring data points are calculated from the Euclidean distance, and the edge weight is calculated by their regression analysis. The conventional LLE algorithm minimizes the reconstruction cost of equation (1), where x _p is the data point represented by the regression analysis.

|| · || in formula (1) is the L2 norm of the ^vector, the ~ _{x p} is the result of regression analysis of the calculated _{x p} by the equation (2). Here, ^~ -indicates a symbol in which ^~ is written immediately above.

In Equation (2), N[x _p ] is a set of K neighboring data points of the data point x _p , and W[i,j] is an edge weight of the data points x _j to x _i . To calculate the weights, the reconstruction cost is minimized so that the sum of the edge weights for the data points x _p in equation (1) is 1. Apply Lagrange's undetermined multiplier method to the Gram matrix to minimize the reconstruction cost. And in particular with the K near the _{x i} and _{x j} for the data points _{x p,} elements of the Gram matrix _{G p} of K × K for the data point _{x p is,} G p _[i, j] the Gram matrix _{G p} When the [i,j] component is used, it is calculated as in Expression (3).

Here, ···· is the inner product of two vectors. When G _p ⁻¹ is an inverse matrix of the matrix G _p , the edge weight of each data point is calculated as in Expression (4).

In order to calculate the edge weight as in Expression (4), it is necessary to obtain the inverse matrix G _p ⁻¹ . The second part of the LLE algorithm determining the vector y _p for the data points x _p by minimizing the embedded cost of formula (5).

When I is an identity matrix and K=(I−W) ^T (I−W) is an LLE kernel of size N×N, and W is a matrix of N×N edge weights, the embedding cost is It can be rewritten as Y ^T KY. Here, · ^T is a transposed matrix of matrix.

There are two restrictions on the embedding cost. The first constraint is that Y ^T Y=I. This means that the solution is rank m. Further, when 0 is a zero vector, Σ _p y _p =0. Therefore, the second part of the LLE algorithm is to solve the optimization problem of equation (6).

Expression (6) is an eigenvalue problem in which the m+1 eigenvectors at the bottom of the kernel K are solutions. Here, the eigenvector having the smallest eigenvalue is discarded because all the elements are the same.

The conventional LLE algorithm requires high calculation cost. In the first part, the matrix G _p is calculated and its inverse matrix is calculated. Since the size of the matrix G _p is K×K and its elements are calculated from the inner product of the vectors of length M, the calculation cost for obtaining the matrix G _p for all data points is O(NK ² M). Becomes Further, since it is necessary to calculate the inverse matrix of the matrix G _p , the calculation cost of O(NK ³ ) is required to obtain the edge weight. Further, the calculation cost of O(N ² M) is required to obtain the neighbors for all the nodes. In addition, the second part requires a calculation cost of O(N ³ ). This is because it is necessary to perform the eigenvalue decomposition of the N×N kernel K. Therefore, the calculation cost of the conventional LLE is O(N(K ² M+K ³ )+N ² M+N ³ ). The memory cost is O(N ² ) because the eigenvalue decomposition of the kernel K is performed. As a result, it is not realistic to apply LLE to large-scale data.

[Configuration of First Embodiment]
The configuration of the first embodiment will be described with reference to FIG. FIG. 1 is a diagram illustrating an example of the configuration of the LLE calculation device according to the first embodiment. As shown in FIG. 1, the LLE calculation device 10 includes a singular value decomposition calculation unit 101, a data point selection unit 102, an approximate distance calculation unit 103, a distance estimation unit 104, a neighborhood determination unit 105, an edge weight calculation unit 106, and rearrangement. It has a unit 107, an LU decomposition calculation unit 108, and a vector calculation unit 109.

[LLE of the first embodiment]
First, the calculation method of LLE and the mathematical background in this embodiment will be described. Then, the processing performed by each functional unit of the LLE calculation device 10 will be described.

In the present embodiment, if two data points are close to each other in a multidimensional space, the LLE is calculated at high speed based on the knowledge that most of their neighborhoods are common. Here, equation (3) as seen in the definition of Gram matrix of, in the conventional LLE, when the data points x _i and x _j and K near the data points x _p, the data point x _p of the Gram matrix [I,j] becomes G _p [i,j]=<x _p −x _i, x _p −x _j >. _Therefore, x p = _x to the data points _{x p} and _{x q} such that _{q, x p} and _{x q} are [i their Gram matrix even share data points _{x i} and _{x j} as a neighboring, j] components are not the same.

Therefore, in the present embodiment, in order to perform the processing at high speed, the Lagrange's undetermined multiplier method is applied to the data points x _p by a method different from the original algorithm, and the matrix C _p is calculated. In the present embodiment, since the [i,j] component of the matrix C _p is C _p [i,j]=<x _i, x _j >, different data can be obtained if the same data point is shared as a neighborhood. The matrix C _p will hold the same components for the points. Therefore, the matrix C _p can be calculated at high speed.

Furthermore, by using this method, the inverse matrix can be calculated at high speed. Specifically, since the matrix C _p has the same components as the matrix C _q by this method, the matrix C _p ⁻¹ can be gradually calculated from the matrix C _q ⁻¹ by using the Woodbury formula. ..

In order to reduce the embedding cost, the present invention calculates the eigenvalue vector of the kernel K at high speed using the inverse power method (see, for example, Reference 1 (Brian Bradie, A Friendly Introduction to Numerical Analysis, 2007, Pearson.)). ). However, in this embodiment, this method is not used for naiveness. This is because it is necessary to calculate the inverse matrix of the kernel K in order to use this method, which causes a calculation cost of O(N ³ ). There is also a problem that a memory cost of O(N ² ) occurs. As described above, the kernel K is calculated as K=(I−W) ^T (I−W). Here, since W is calculated from edge weights from neighboring data points, the kernel K has a sparse data structure. However, even if the kernel K is sparse, its inverse matrix has a dense structure, which requires high calculation cost and memory cost.

In order to solve these problems, in this embodiment, the LU decomposition of the matrix I-W is calculated, and the eigenvector of the kernel K is calculated. Here, the LU decomposition is a method of decomposing a matrix into a lower triangular matrix and an upper triangular matrix (for example, reference 2 (William H. Press, Saul A. Teukolsky, William T Vetterling, Brian P. Flannery, Numerical Recipes 3rd. Edition, 2007, Cambridge University Press.)). Since the eigenvector can be calculated from the lower triangular matrix and the upper triangular matrix, the calculation of the kernel K can be avoided. Further, since the lower triangular matrix and the upper triangular matrix have a sparse data structure, the calculation cost and the memory cost can be suppressed.

The present embodiment has the same calculation result as the conventional LLE algorithm. This is because the present embodiment performs calculations that minimize the reconstruction cost and the embedding cost. That is, this embodiment can reduce the dimension with a small calculation cost and memory cost without changing the calculation result of the original algorithm.

[Calculation of Edge Weight in First Embodiment]
Next, a method for calculating edge weights at high speed is described. In the present invention, the Lagrange undetermined multiplier method is used to calculate the matrix C _p of the data points x _p . The Lagrangian undetermined multiplier method can obtain a steady state of a function having a plurality of variables under constant constraints (for example, see Reference 3 (Christopher M. Bishop, Pattern Recognition and Machine Learning, 2010, Springer.). reference). In LLE, the variable is the edge weight that minimizes the reconstruction cost, and the constraint is that the sum of the edge weights for each data point is one. Letting γ be a Lagrange undetermined multiplier, the Lagrange function L for each data point is:

Therefore, for the edge weight W[p,i] of the data point x _i such that x _i εN[x _p ], the condition for the function L to be in a steady state for W[p,i] is ∂L /∂W[p,i]=0. Similarly, the condition of the function L for γ is ∂L/∂γ=0. Therefore, the condition that the function L is in the steady state is as shown in Expression (12).

In Equation (12), C _p is a (K+1)×(K+1) matrix, and w and p are column vectors of length K+1 given below.

Here, x _i , x _j εN[x _p ]. The elements of the matrix _{C p} from equation _{(13) C p [i,} j] = <x i, x j> can be calculated with the elements of the matrix _{C p} be can be calculated independently of _{x p} I understand. Therefore, if the data points x _p and x _q share the same neighborhood, the elements of the matrix C _p of the data points x _p can be calculated at high speed by referring to the matrix C _q of the data points x _q . If d is the number of different neighbors at the data points x _p and x _q (that is, d=K−|N[x _p ]∩N[x _q ]|), the matrix C _p is obtained in this embodiment. Requires a calculation cost of O(dKM).

[Calculation of Inverse Matrix of First Embodiment]
Then we describe a method for obtaining the inverse matrix _C ^{p -1} of the matrix _{C p} at high speed using the formula Woodbury. As shown in Expression (14), the element of i in the vector w corresponds to the edge weight W[p,i]. Therefore, the edge weight can be calculated as w=C _p ⁻¹ from the equation (12) in which C _p ⁻¹ is the inverse matrix of the matrix C _p . But since the directly from the matrix C _p Request matrix C _p ^-1 consuming computational cost O (K ^3), calculate the cost of this approach is high.

A method for calculating edge weights at high speed using the Woodbury formula will be described. First, when d=1, that is, when the data point x _p has only one different neighborhood with respect to the data point x _q (K _{- | N [x p] ∩N} [x q] | = 1) is described. Let x _{k be} the only different neighbor data point at data point x _p , and let x _k be the data point neighbor of data point x _q with respect to x _k . Therefore _{x k = N [x p]} -N [x p] is _{∩N [x q] x k =} N [x q] -N [x p] ∩N [x q] become.

Here, it is assumed that the variance of each vector corresponding to a data point is 1, that is <x _i, x _i >=1. This assumption will be removed later. The following matrix ΔC is used when updating the inverse matrix C _p ⁻¹ corresponding to the data points x _k and x _{k′ in} the k th row and column of the matrices C _p and C _q .

The lemma 1 (lemma 1) shown in FIG. 2 holds for the matrix ΔC. FIG. 2 is a diagram showing Lemma 1. Further, Lemma 1 shows that the matrix ΔC is a symmetric matrix calculated from the asymmetric matrices C _p and C _q , and only k rows and columns of the matrix ΔC have nonzero elements. From this lemma, Lemma 2 (lemma 2) for the matrix ΔC shown in FIG. 3 is established. FIG. 3 is a diagram showing Lemma 2.

Since V ^T V=I and ΔCv ₁ =||f||v ₁ and ΔCv ₂ =−||f||v ₂ , the lemma 2 and the formula (19) are the matrix ΔC. It can be seen that it corresponds to the eigenvalue decomposition of rank 2. By using this formula, the matrices V and F can be obtained at a calculation cost of O(K). This is because the length of the vector f is K+1. Therefore, the eigenvalue decomposition of the matrix ΔC can be obtained from the lemma 2 at a calculation cost of O(K).

It can be progressively updated inverse matrix _C ^{p -1} data points _{x p} from the inverse matrix _C ^{q -1} of the data points _{x q} by using the official Woodbury. Since a _C p = _{C q} + _ΔC from equation (16), consisting of Lemma 2 and _C p = _{C q} + _ΔC. Therefore, the inverse matrix C _p ⁻¹ can be calculated as follows by using the Woodbury formula.

In Formula (20), since F ⁻¹ and V ^T C _q ⁻¹ V are 2×2 matrices, the calculation cost for _obtaining the inverse matrix (F ⁻¹ +V ^T C _q ⁻¹ V) ⁻¹ is constant. .. The sizes of the matrices C _q ⁻¹ and V are (K+1)×(K+1) and (K+1)×2, respectively. Therefore, the calculation cost of O(K ² ) is required to obtain C _p ⁻¹ from the equation (20).

Although the inverse matrix C _p ^-1 can be obtained at high speed from the equation (20), this equation assumes that <x _i, x _i >=1 holds in the matrix C _p . Here, the matrix C _p ^-1 is updated using this assumption. And when the corresponding matrix _C'p matrix _{C p,} the matrix ΔC of (K + 1) × (K + 1) for the matrix _C'p is calculated as Equation (21).

It can be seen from the matrix ΔC′ that the differences between the matrices C′ _p and C _p are limited to [k, k] elements, and the matrix ΔC′ drops the assumption that <x _i , x _i >=1 in matrix C _p. It turns out that it corresponds to. Therefore, if the assumption is dropped, it can be seen that C _p =C _q +ΔC+ΔC′ holds. If e _k is a unit vector that is 0 except for the k-th element, the matrix ΔC′ can be obtained as in Expression (22).

Therefore the inverse matrix of the matrix _C'p after formula (20) by using the official Woodbury the ^{_(C'p) -1} can be computed as in Equation (23).

Similar to Expression (20), the inverse matrix can be calculated at high speed from Expression (23) at a calculation cost of O(K ² ). In equations (20) and (23), it was assumed that d=1 holds, that is, for data point x _q , data point x _p has one different neighboring data point. However, even if d>1, iteratively repeats the above process for each data point such that x _i ε{x _j |x _j εN[x _p ] and x x _j /εN[x _q ]}. The inverse matrix can be obtained at high speed by repeatedly applying the two. Note that, here, x/εN indicates that x is not included in the set N.

[Calculation of data points in the vicinity of the first embodiment]
Next, we describe a method to calculate the neighboring data points at high speed for each data point. In the present invention, the approximate distance is calculated using SVD, and the neighborhood is calculated. Since the SVD can find the lower limit value of the distance, it can accurately find the neighborhood. If ^~ _xp =( ^~ _xp [1], ^~ _xp [2], ..., ^~ _xp [s]) is a data point _xp =( _xp [1], _xp [2], ..., x _p [M]) is the dimension reduction by SVD of rank s, and E[x _p , x _q ] is the Euclidean distance between data points x _p and x _q , then SVD defines the lower bound of the distance. Giving means that E[ _xp , _xq ] ≥ E[ ^~ _xp , ^~ _xq ] holds.

In the existing research to calculate the neighborhood, first, the approximate distance is calculated, and then, if it is determined that the data point may be the neighborhood due to the approximate distance, the Euclidean distance is calculated (for example, reference 4 ( Dennis Shasha, High Performance Discovery in Time Series: Techniques and Case Studies, 2004, Springer-Verlag New York Inc..)). Note that the SVD of each data point can be obtained at a calculation cost of O(NMlogs) using an existing method (for example, Reference 5 (Nathan Halko, Per-Gunnar Martinsson, Joel A. Tropp, Finding Structure with Randomness: Probabilistic Algorithms for Constructing Approximate Matrix Decompositions, 2011, SIAM Review.)).

In order to improve the accuracy of the approximate distance, the Euclidean distance is estimated using the structure of the neighborhood graph. The data point x _q is near the data point x _r . In this case ^{_{E [~ x q, ~ x}} r] by using SVD and _{E [x} q, x _r] holds. Further, if if calculating the Euclidean distance of the data point _{x p} data points _{x r} in the process of obtaining the vicinity ^{_{of, E [~ x p, ~}} x r] and E [x p, x r] holds. Assuming that E_[x _p , x _q ] is the estimated distance from the data point x _p to x _q , the estimated distance is calculated as in Expression (24) in this embodiment. Here, ·_ indicates a symbol with _ written directly under ·.

Here _u r _{[x p]} is ^{_{u r [x p] = (}} ~ x p [s + 1] - ~ x r [s + 1], ~ x p [s + 2] - ~ xr [s + 2], ..., ~ x a vector of length Ms such that _p [M]- ^to _xr [M]). Lemma 3 (Lemma 3) shown in FIG. 4 holds for the estimated distance. FIG. 4 is a diagram showing Lemma 3.

Used for Lemma 3 error for _{E [x} p, x _q] of the estimated distance _{E [x} p, x _q] is less than ^{_{E [~ x p, ~ x}} q], E [x p, x q] also shows that can improve the accuracy of the approximation of ^{_{E [~ x p, ~ x}} q] by. Further, E[x _p, x _q ]≧E[x _p, x _q ] of Lemma 3 shows that the property of the lower limit value is satisfied, that is, the neighborhood can be accurately obtained.

In this method, the L ₂ norm of the vectors u _r [x _p ] and u _r [x _q ] used in equation (24) is calculated as in equations (25) and (26).

Here if calculates the Euclidean distance of the data point x _p and x _r to check the x _r as the vicinity of the data point x _p, if the data point x _q is in the vicinity of x _r, in the formula (25) By using (26), the calculation cost for _obtaining E[x _p, x _q ] is O(1). In this case this, _{E [x} p, xr] and ^{_{E [~ x p, ~ x}} r] and _{E [x} q, x _r] and ^{_{E [~ x q, ~ x}} r] the calculation of O (1) This is because it can be calculated at cost.

Next, a method for calculating the eigenvector at high speed will be described. The conventional LLE algorithm calculates the embedding Y in the low dimension of the data set X by the optimal solution to the embedding cost. For the edge weight matrix W, the LLE kernel is calculated as K=(I−W) ^T (I−W), but this optimal solution corresponds to the bottom eigenvector of the kernel K. Therefore, if λ _i is |λ ₁ |≧|λ ₂ |≧. ． ． If the i-th eigenvalue of the kernel K such that ≧|λ _N |, and z _i is the corresponding eigenvector, then the m-dimensional embedding Y for the dataset X is Y=[ from the m+1 bottom eigenvectors. z _N-m, z _N-m+1 ,. ． ． , Z _N−1 ].

In order to calculate the eigenvector at high speed, the inverse power method is used in this embodiment. The inverse power method is a method of calculating the bottom eigenvector from the inverse matrix of the matrix for calculating the eigenvalue, and specifically, the eigenvector z _N can be calculated from the inverse matrix K−1. Therefore, if a ₀ is an arbitrary column vector of length N, the base eigenvalue λ _Na and its eigenvector z _N can be calculated by applying the power method to a _τ =K ⁻¹ a _τ ⁻¹ . However, since the inverse matrix K-1 has a dense structure, the calculation cost required for this method is O(N ³ ) and the memory cost is O(N ² ).

In order to avoid the inverse matrix of the direct kernel K, we compute the LU decomposition of I-W. Specifically, when L and U are a lower triangular matrix and an upper triangular matrix, the matrix L and U are calculated as LU=I-W. As a result, equation (27) is obtained.

In the present embodiment, the powers of the matrices L and U are used to calculate a _τ-1 to a _τ as in Expression (28).

In equation (28), forward substitution and backward substitution are applied to the lower triangular matrix and the upper triangular matrix in order to calculate b, b′, b″, and a _τ (for example, see Reference 2). Equation (28) shows that obtained in the form of ^{_{a τ-1 = U T L}} T LU aτ vector a _tau from the vector a _tau-1 from the lower and upper triangular matrices. Since the equation ^{(27) K = U T L} T LU holds, it is possible to calculate the ^{_{a τ = K -1 a τ-}} 1 without using the inverse matrix of kernel K directly. The base eigenvalue λ _N and the eigenvector z _N can be calculated by applying the power method until the vector a _τ converges as in Expressions (28) to (29).

Since the lower triangular matrix and the upper triangular matrix of the kernel K have a sparse structure, the present invention can quickly calculate the vector a _τ from the equation (28). Next, a method for calculating the other base eigenvalues λ _N-i and eigenvectors z _N-i (i=1, 2,..., _M ) will be described. The Hotelling method is famous as a method for obtaining such an eigenvalue and an eigenvector (for example, see Reference 1). The basic idea of the Hotelling method is to shift the matrix for calculating the eigenvalues so that the maximum eigenvalue is 0 and other eigenvalues remain unchanged. As a result, the second largest eigenvalue becomes the largest eigenvalue in the shifted matrix. Specifically, λ _N-i and z _N-i are calculated by calculating the eigenvalues and eigenvectors of the matrix H _i of Expression (30) using the power method.

Here, when h _i,τ is a column vector of length N, the eigenvalue λ _N-i and the eigenvector z _N-i are hi _i,τ =H _i hi _i,τ-1 You can ask. However, this method has a problem that a matrix for H _i has a dense structure, memory cost of O computational cost (N ³⁾ and O (N ²⁾ is required.

In order to avoid directly calculating the matrix H _i , in the present invention, the vector h _i,τ is calculated from the vector h _i,τ−1 as shown in Expression (31).

Here, similarly to the equation (28), the vector a _τ is calculated from the vector h _{i, τ-1} as shown in the equation (32) using forward substitution and backward substitution as follows.

In equation (31), the vector of length N z _N−j z _N−j ^T hi _{, τ−1} /λ _N−j is calculated from the right to the left at the calculation cost of O(N). Ask. Since K ⁻¹ =(U ^T L ^T LU) ⁻¹ , the formula (33) is established from the formulas (30), (31), and (32).

Expression (33) shows that the vector h _i,τ can be calculated from the expressions (31) and (32) by using the matrix H _{i of the} Hotelling method. In the present embodiment, the eigenvalue λ _N−i of the kernel K and the eigenvector z _N−i (1≦i≦m) are calculated by the power method based on the equation (34) until the vector h _i,τ converges.

Since the matrix H _i is not calculated directly, the present invention can calculate eigenvalues and eigenvectors at high speed.

Next, a method for calculating the LU decomposition of the matrix I-W used in the inverse power method at high speed will be described. The elements of the matrices L and U can be calculated as in equations (35) and (36) by the Kraut's algorithm (see eg reference 2) if the matrix W′ is given as W′=I−W.

It can be seen that in equations (35) and (36), the columns of the matrices L and U can be calculated from left to right, and the elements in each column from top to bottom. Therefore, the elements of the matrices L and U can be calculated from the elements above and to the left of the corresponding matrices W', L, U.

In this embodiment, the number of non-zero elements in the matrices L and U is reduced to calculate the LU decomposition at high speed. This means that if L[i,k]=0 or U[k,j]=0, then the calculation of L[i,k]U[k,j] is effectively performed from equations (35) and (36). This is because it can be omitted. In this embodiment, the upper and left elements of the matrices L and U become zero if the corresponding elements of the matrix W′ are zero, and the upper and left elements of the matrix W′ are W′=I−W. Therefore, if the left and upper elements of the corresponding matrix W become zero, the fact that they become zero is used. Since the matrix W is calculated from edge weights, LU decomposition can be performed at high speed if the upper and left elements become sparse by rearranging the data points of the matrix W.

In this embodiment, the LU decomposition is performed at high speed by rearranging the data points in ascending order of the order. Here, the order is the number of edges generated from the data point. By rearranging the elements, the elements above and to the left of the matrix W can be made sparse, so that the LU decomposition of the matrix IW can be performed at high speed. Sorting the data points can be done at a computational cost of O(N) by using a distributed counting sort (see eg Reference 6 (Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest, Clifford Stein, Introduction to Algorithms, 2009, The MIT Press.)). If n is the number of non-zero elements of each average row vector in the triangular matrix, the memory cost for holding the matrices L and U is O(Nn). Further, from the equations (35) and (36), the LU decomposition of the matrix I-W can be performed at a calculation cost of O(Nn ² ).

[Algorithm of First Embodiment]
The algorithm of this embodiment will be described together with the processing of each functional unit of the LLE calculation device 10 with reference to FIG. FIG. 5 is a diagram illustrating an example of an algorithm of the LLE calculation process according to the first embodiment.

In Algorithm 15, S is the set of data points selected to calculate the edge weights, D is the set of data points in the neighborhood of the selected data points, and C[x _p ] is the data in the neighborhood graph. A set of data points directly connected to the point x _p . In the algorithm 15, the SVD is calculated in order to calculate the neighborhood at a high speed. However, if s is a rank of the SVD such that s=2 ^l , in the present embodiment, the ranks are 2 ¹ , 2 ² ,. ． ． , 2 ^l, and gradually increase to find the neighborhood. From the viewpoint of calculation cost, if the approximation is rough, the approximate distance can be calculated at high speed, but the accuracy of the approximation is not high. On the other hand, if the accuracy is fine, the approximate distance can be calculated with high accuracy, but the calculation cost will be high. Therefore, in the present embodiment, data points that are not close to each other are pruned by coarse approximation, and data points that are likely to be close to each other are calculated by fine approximation to calculate an approximate distance. As a result, the present invention can determine the neighborhood by appropriately setting the approximation accuracy according to the distance of each data point.

As shown in Algorithm 15, first, the matrix X, the low-dimensional number m, and the rank s (=2 ^l ) of the SVD are input. Then, the singular value decomposition calculation unit 101 calculates the singular value decomposition (SVD) of the multidimensional matrix having a plurality of data points. Specifically, the singular value decomposition calculation unit 101 initializes S and calculates the SVD of the matrix X (1-2nd line).

The data point selection unit 102 selects a data point to be calculated from a plurality of data points. Specifically, the data point selection unit 102 selects a data point having many edges to other data points (fourth row). This is because in the present embodiment, the distance to other data points is estimated using the neighborhood graph. Then, the data point selection unit 102 initializes D and θ (5th to 6th lines), and initializes N[x _p ] by adding K data points at which the distance becomes ∞ (7th line). Eye).

The approximate distance calculation unit 103 calculates the approximate distance between each of the plurality of data points and the data point to be calculated based on the singular value decomposition. At this time, the approximate distance calculation unit 103 gradually increases the rank from 2 ¹ to 2 ¹ in order to calculate the neighborhood, but if the approximate distance becomes larger than θ, the data point is pruned (9- (Line 13).

The distance estimation unit 104 also calculates an estimated value of the Euclidean distance between each data point whose approximate distance is equal to or less than a predetermined value and the data point to be calculated. That is, the distance estimation unit 104 estimates the distance using the graph structure when the data point cannot be pruned using the SVD (14th to 18th lines).

The neighborhood determining unit 105 calculates the Euclidean distance between each data point whose estimated value of the Euclidean distance is less than or equal to a predetermined value and the data point to be calculated, and the calculated Euclidean distance is less than or equal to the predetermined value. Establish the point to a data point near the data point to be calculated. When the proximity determining unit 105 determines the data point to be in the vicinity of the selected data point, the proximity determining unit 105 updates N[x _p ] using the distance of the data point (lines 19-25).

Then, the neighborhood determination unit 105 finds a data point having the largest number of neighborhoods (26th line). This is because the edge weight calculation unit 106 calculates the edge weight using the neighborhood of another data point.

The edge weight calculation unit 106 calculates edge weights between a data point to be calculated and each data point near the data point to be calculated. The edge weight calculation unit 106 calculates the edge weight using the Woodbury formula when there is a shared neighborhood (lines 27-28), and otherwise uses the Gram matrix to calculate the edge weight. (Lines 29-30). In this way, the edge weight calculation unit 106 can use the Woodbury formula to calculate the inverse matrix of the data points in the vicinity of the data point to be calculated from the inverse matrix of the data points in the vicinity of the calculated data point. ..

In other words, the edge weight calculation unit 106 calculates if the data points in the vicinity of the calculated data points for which the edge weights have been calculated include at least some of the data points in the vicinity of the data point to be calculated. Based on the calculation result of the edge weight for the data point, the edge weight between the data point to be calculated and each data point in the vicinity of the data point to be calculated can be calculated.

The rearrangement unit 107 rearranges the matrix representing the edge weights of a plurality of data points in ascending order of the order (lines 32-33). The LU decomposition calculation unit 108 also calculates the LU decomposition of the adjacency matrix of a plurality of data points.

The vector calculation unit 109 calculates the eigenvector of the LLE kernel for a plurality of data points based on the edge weight. For example, the vector calculation unit 109 calculates the eigenvector and eigenvalue of the LLE kernel using the inverse power multiplication method based on the LU decomposition calculated by the LU decomposition calculation unit 108 (lines 34-36). ).

The edge weight calculation unit 106 may always execute the 29th line in the 28th to 30th lines. In this case, similarly to the conventional LLE, the edge weight calculation unit 106 calculates the edge weight using the Gram matrix. Further, the processing by the rearrangement unit 107 and the LU decomposition calculation unit 108 may not be executed. In this case, the vector calculation unit 109 calculates the eigenvector and the eigenvalue by the same method as the conventional LLE. Further, for the algorithm 15 shown in FIG. 5, the properties shown in FIG. 6 (Theorems 1 to 3) hold. FIG. 6 is a diagram showing Theorem 1, Theorem 2 and Theorem 3.

[Processing of First Embodiment]
The processing flow of the LLE calculation device 10 will be described with reference to FIG. 7. FIG. 7 is a flowchart showing a processing flow of the LLE calculation device according to the first embodiment. First, the matrix X, the number of dimensions m, and the rank s are input to the LLE calculation device 10 (step S101). Next, the singular value decomposition calculation unit 101 initializes S to an empty set (step S102). Then, the singular value decomposition calculation unit 101 calculates the SVD of the rank s of the matrix X (step S103).

The data point selection unit 102 selects a data point to be calculated from a plurality of data points while increasing i from 1 to N (steps S104, S105, S132). Then, the data point selection unit 102 initializes D and θ (steps S106 and S107), and further initializes N[x _p ] by adding K data points at which the distance is ∞ (step S106). S108).

The approximate distance calculation unit 103 calculates an approximate distance between each of the plurality of data points and the data point to be calculated based on the singular value decomposition (steps S109, S110, S111). At this time, the approximate distance calculation unit 103 gradually increases the rank from 2 ¹ to 2 ¹ in order to calculate the neighborhood, but if the approximate distance becomes larger than θ (step S112, Yes), the data points are calculated. Pruning is performed (step S113). If the approximate distance is not larger than θ (No in step S112), the approximate distance calculation unit 103 proceeds to the next process (step S114).

Here, when there is no data point that has not been pruned (No in step S115), the LLE calculation device 10 proceeds to step S128. On the other hand, when data points that have not been pruned remain (step S115, Yes), the distance estimation unit 104 sets a distance between each data point whose approximate distance is equal to or less than a predetermined value and the data point to be calculated. The estimated value of the Euclidean distance is calculated (steps S116 and S117). If the estimated distance value is larger than θ (step S118, Yes), the distance estimation unit 104 prunes the data point (step S119). If the approximate distance is not larger than θ (No in step S118), the distance estimation unit 104 proceeds to the next process (step S120).

Here, when there is no data point that has not been pruned (No in step S121), the LLE calculation device 10 proceeds to step S128. On the other hand, when data points that have not been pruned remain (step S121, Yes), the neighborhood determining unit 105 determines each of the data points for which the Euclidean distance estimated value is equal to or less than a predetermined value and the data point to be calculated. The Euclidean distance between the calculated data points is calculated (step S122), and the data points whose calculated Euclidean distance is less than or equal to a predetermined value (step S123, Yes) are determined as data points near the data point to be calculated, and N[ x _p ] is updated using the distance of the data point (steps S124, S125, S126). Further, the proximity determining unit 105 updates θ (step S127). Then, the neighborhood determination unit 105 finds a data point having the largest number of neighbors (step S128).

The edge weight calculation unit 106 calculates the edge weight using the Woodbury formula when there is a shared neighborhood (step S129, Yes) (step S130), and when not (step S129, No), the gram matrix is calculated. The edge weight is calculated by using (step S131).

The rearrangement unit 107 rearranges the matrix representing the edge weights of a plurality of data points in ascending order of the order (step S133). Also, the LU decomposition calculation unit 108 calculates the LU decomposition of the adjacency matrix of the plurality of data points (step S134).

The vector calculation unit 109 calculates the eigenvector of the LLE kernel for a plurality of data points based on the edge weight. For example, the vector calculation unit 109 calculates the eigenvector and eigenvalue of the LLE kernel using the inverse power method based on the LU decomposition calculated by the LU decomposition calculation unit 108 (steps S135, S136, S137, S138). Then, the LLE calculation device 10 outputs the dimension-reduced matrix Y (step S139).

[Effects of First Embodiment]
The singular value decomposition calculation unit 101 calculates the singular value decomposition of a multidimensional matrix having a plurality of data points. Further, the data point selection unit 102 selects a data point to be calculated from the plurality of data points. Further, the approximate distance calculation unit 103 calculates an approximate distance between each of the plurality of data points and the data point to be calculated based on the singular value decomposition. The distance estimation unit 104 also calculates an estimated value of the Euclidean distance between each data point whose approximate distance is equal to or less than a predetermined value and the data point to be calculated. Also, the neighborhood determining unit 105 calculates the Euclidean distance between each data point whose estimated value of the Euclidean distance is less than or equal to a predetermined value and the data point to be calculated, and the calculated Euclidean distance is less than or equal to the predetermined value. Establish a data point as a data point near the data point to be calculated. Further, the edge weight calculation unit 106 calculates the weight of the edge between the data point to be calculated and each data point in the vicinity of the data point to be calculated. Further, the vector calculation unit 109 calculates the eigenvectors of the LLE kernel for the plurality of data points based on the edge weights. As described above, in the present embodiment, the target for which the Euclidean distance is actually calculated by the neighborhood determining unit 105 is reduced in advance by using the approximate distance and the estimated value of the Euclidean distance. Therefore, according to this embodiment, it is possible to reduce the calculation cost and the memory cost when the dimension is reduced by LLE.

The edge weight calculation unit 106 calculates the calculated data points when the data points in the vicinity of the calculated data points for which the edge weights have been calculated include at least some of the data points in the vicinity of the data point to be calculated. The edge weight between the data point to be calculated and each of the data points in the vicinity of the data point to be calculated can be calculated based on the calculation result of the edge weight with respect to. In this way, for data points having the same neighborhood, by using the already calculated result, it is possible to reduce the cost required for calculating the edge weight.

The edge weight calculation unit 106 can use the Woodbury formula to calculate the inverse matrix of the data points in the vicinity of the data point to be calculated from the inverse matrix of the data points in the vicinity of the calculated data point. As described above, by using the Woodbury formula, it is possible to reduce the cost particularly required for the calculation of the inverse matrix.

The sorting unit 107 sorts the matrix representing the edge weights of a plurality of data points in ascending order of the order. The LU decomposition calculation unit 108 also calculates the LU decomposition of the adjacency matrix of a plurality of data points. At this time, the vector calculation unit 109 calculates the eigenvector of the LLE kernel based on the LU decomposition calculated by the LU decomposition calculation unit 108. By performing LU decomposition in this way, it is possible to reduce the cost required to calculate the eigenvector.

[System configuration, etc.]
Further, each constituent element of each illustrated device is functionally conceptual, and does not necessarily have to be physically configured as illustrated. That is, the specific form of distribution/integration of each device is not limited to the one shown in the figure, and all or part of the device may be functionally or physically distributed/arranged in arbitrary units according to various loads and usage conditions. It can be integrated and configured. Furthermore, each processing function performed in each device is realized in whole or in part by a CPU (Central Processing Unit) and a program that is analyzed and executed by the CPU, or a hardware by a wired logic. Can be realized as.

Further, of the processes described in the present embodiment, all or part of the processes described as being automatically performed may be manually performed, or the processes described as manually performed may be performed. All or part of the process can be automatically performed by a known method. In addition, the processing procedures, control procedures, specific names, and information including various data and parameters shown in the above-mentioned documents and drawings can be arbitrarily changed unless otherwise specified.

[program]
As one embodiment, the LLE calculation device 10 can be implemented by installing an LLE calculation program that executes the above LLE calculation as package software or online software in a desired computer. For example, by causing the information processing apparatus to execute the above LLE calculation program, the information processing apparatus can be caused to function as the LLE calculation apparatus 10. The information processing device mentioned here includes a desktop or notebook personal computer. In addition, the information processing apparatus also includes a mobile communication terminal such as a smartphone, a mobile phone, a PHS (Personal Handyphone System), and a slate terminal such as a PDA (Personal Digital Assistant) in its category.

The LLE calculation device 10 can also be implemented as a LLE calculation server device that uses a terminal device used by a user as a client and provides the client with a service related to the above LLE calculation. For example, the LLE calculation server device is implemented as a server device that provides an LLE calculation service in which a multidimensional matrix is input and a dimension-reduced matrix is output. In this case, the LLE calculation server device may be implemented as a Web server, or may be implemented as a cloud that provides the above-mentioned LLE calculation service by outsourcing.

FIG. 8 is a diagram illustrating an example of a computer that executes the LLE calculation program. The computer 1000 has, for example, a memory 1010 and a CPU 1020. The computer 1000 also has a hard disk drive interface 1030, a disk drive interface 1040, a serial port interface 1050, a video adapter 1060, and a network interface 1070. These units are connected by a bus 1080.

The memory 1010 includes a ROM (Read Only Memory) 1011 and a RAM (Random access memory) 1012. The ROM 1011 stores, for example, a boot program such as a BIOS (Basic Input Output System). The hard disk drive interface 1030 is connected to the hard disk drive 1090. The disk drive interface 1040 is connected to the disk drive 1100. For example, a removable storage medium such as a magnetic disk or an optical disk is inserted into the disk drive 1100. The serial port interface 1050 is connected to, for example, a mouse 1110 and a keyboard 1120. The video adapter 1060 is connected to the display 1130, for example.

The hard disk drive 1090 stores, for example, an OS (Operating System) 1091, an application program 1092, a program module 1093, and program data 1094. That is, the program that defines each process of the LLE computing device 10 is implemented as a program module 1093 in which code executable by a computer is described. The program module 1093 is stored in the hard disk drive 1090, for example. For example, the hard disk drive 1090 stores a program module 1093 for executing the same processing as the functional configuration of the LLE computing device 10. The hard disk drive 1090 may be replaced by SSD.

Further, the setting data used in the processing of the above-described embodiment is stored as the program data 1094 in the memory 1010 or the hard disk drive 1090, for example. Then, the CPU 1020 reads out the program module 1093 and the program data 1094 stored in the memory 1010 or the hard disk drive 1090 to the RAM 1012 as necessary and executes them.

The program module 1093 and the program data 1094 are not limited to being stored in the hard disk drive 1090, but may be stored in, for example, a removable storage medium and read by the CPU 1020 via the disk drive 1100 or the like. Alternatively, the program module 1093 and the program data 1094 may be stored in another computer connected via a network (LAN (Local Area Network), WAN (Wide Area Network), etc.). Then, the program module 1093 and the program data 1094 may be read by the CPU 1020 from another computer via the network interface 1070.

10 LLE Calculator 101 Singular Value Decomposition Calculator 102 Data Point Selector 103 Approximate Distance Calculator 104 Distance Estimator 105 Neighbor Determinator 106 Edge Weight Calculator 107 Sorting Unit 108 LU Decomposition Calculator 109 Vector Calculator

Claims (6)

A singular value decomposition calculator that calculates the singular value decomposition of a multidimensional matrix having multiple data points,
A data point selection unit that selects a data point to be calculated from the plurality of data points,
For each of the plurality of data points, an approximate distance between the data point to be calculated, an approximate distance calculation unit that calculates based on the singular value decomposition,
A distance estimation unit that calculates an estimated value of the Euclidean distance between each of the data points whose approximate distance is less than or equal to a predetermined value and the data point to be calculated,
The Euclidean distance between each of the data points whose estimated value of the Euclidean distance is less than or equal to a predetermined value and the data point to be calculated is calculated, and the calculated Euclidean distance is equal to or less than a predetermined value. A neighborhood determining unit that determines the data points near the data point to be calculated,
An edge weight calculation unit that calculates a weight of an edge between the data point to be calculated and each of the data points in the vicinity of the data point to be calculated,
A vector calculator that calculates an eigenvector of the LLE kernel for the plurality of data points based on the edge weights;
An LLE calculation device comprising: The edge weight calculator calculates the edge weights when the data points in the vicinity of the calculated data points for which the edge weights have been calculated include at least a part of the data points in the vicinity of the data point to be calculated. Calculating an edge weight between the data point to be calculated and each of the data points in the vicinity of the data point to be calculated based on the calculation result of the weight of the edge for the processed data point. The LLE calculation device according to claim 1. The edge weight calculator may use a Woodbury formula to calculate an inverse matrix of data points in the vicinity of the data point to be calculated from an inverse matrix of data points in the vicinity of the calculated data point. The LLE calculation device according to claim 2. A sorting unit that sorts the matrix representing the edge weights of the plurality of data points in ascending order of order,
An LU decomposition calculator for calculating an LU decomposition of an adjacency matrix of the plurality of data points,
Further has
The LLE calculation device according to claim 1, wherein the vector calculation unit calculates the eigenvector of the LLE kernel based on the LU decomposition calculated by the LU decomposition calculation unit. An LLE calculation method executed by an LLE calculation device, comprising:
A singular value decomposition calculation step for calculating the singular value decomposition of a multidimensional matrix having a plurality of data points,
A data point selecting step of selecting a data point to be calculated from the plurality of data points,
For each of the plurality of data points, an approximate distance between the data points to be calculated, an approximate distance calculation step of calculating based on the singular value decomposition,
A distance estimating step of calculating an estimated value of the Euclidean distance between each of the data points whose approximate distance is equal to or less than a predetermined value and the data point to be calculated,
The Euclidean distance between each of the data points whose estimated value of the Euclidean distance is less than or equal to a predetermined value and the data point to be calculated is calculated, and the calculated Euclidean distance is equal to or less than a predetermined value. A neighborhood confirmation step of confirming data points in the vicinity of the data point to be calculated,
An edge weight calculation step of calculating a weight of an edge between the data point to be calculated and each of the data points in the vicinity of the data point to be calculated,
A vector calculation step of calculating an eigenvector of the LLE kernel for the plurality of data points based on the edge weights;
An LLE calculation method comprising: An LLE calculation program for causing a computer to function as the LLE calculation device according to claim 1.

Images