JP5186322B2 - Time series data analysis system, method and program - Google Patents

Description

  The present invention relates to a system, method, and program for analyzing time series data acquired from a sensor or the like.

  Conventionally, a process of monitoring presence / absence of an abnormality by monitoring time-series data acquired from a sensor attached to a manufacturing apparatus or a monitoring apparatus has been performed.

  On the other hand, in recent years, an attempt is made to run a vehicle equipped with a sensor and analyze time-series data obtained from such a sensor. For automobiles, engine cooling water temperature sensor, engine intake air temperature sensor, oil temperature sensor, intake pipe pressure sensor for fuel injection device, turbocharger supercharging pressure sensor, throttle position sensor, steering rudder angle sensor , Vehicle height sensor, liquid level sensor, electromagnetic rotational speed sensor, knock sensor, acceleration sensor, flow rate sensor, oxygen sensor, lean air-fuel ratio sensor, etc. It reaches 100.

One of the things that we would like to do using the time series output data of a plurality of sensors is abnormality diagnosis. That is, FIG. 1B shows test data for which such an abnormality diagnosis is desired. Such test data
It is assumed that D = {x ^(t) | x ^(t) ∈ R ^p , t = 1, 2,..., N}. Here, R ^p represents a p-dimensional real number space, and therefore x ^(t) is a p-dimensional vector having a real number component. The one component corresponds to one sensor. t is an index of time and corresponds to a time for each increment. Therefore, x ⁽¹⁾ , x ⁽²⁾ ,..., X ^(N) are N vectors along the time series.

On the other hand, it is assumed that normal data for reference is prepared in the same format as test data. This is, for example, a sensor value measured during normal driving of the automobile. Normal data for reference is
It is assumed that D˜ = {x˜ ^(t) | x˜ ^(t) ∈R ^p , t = 1, 2,..., N}.

  Then, performing abnormality diagnosis means comparing D to D to determine which sensor indicates an abnormality.

As an approach for performing such a comparison, methods conventionally known include the following.
(1) A technique using time series modeling such as an autoregressive model as an approach from the limit having a strong time correlation, for example, a technique disclosed in Japanese Patent Laid-Open No. 6-109498.
(2) A technique that uses a statistical analysis method based on the assumption of independent identical distribution.
(3) A technique using chaos time series analysis based on the law of motion of the system, for example, a technique as disclosed in JP-A-8-278815.

  However, especially in the case of automobile sensor data, (a) it is very dynamic, (b) there is often a strong correlation between time series signals, and (c) no prior knowledge is given for each sensor. There is. That is, since the signal is very dynamic and complex, the time-series aspect changes greatly if the data is different. Therefore, even if the signals are from the same sensor, it does not make sense to “compare with each other”. In other words, stationarity from the time series cannot be assumed. For this reason, it is difficult to create a time series prediction model using a traditional time series prediction method and to solve the problem using a technique of detecting a deviation from the prediction.

  However, the sensor data cannot be considered to follow an independent and identical distribution, and the probability distribution is generally complicated. Furthermore, there is generally a dependency relationship between sensors, and each sensor must be treated as a multivariate system. Here, the dependency relationship is, for example, the relationship of “acceleration when the accelerator is depressed”, which is related to the output from the throttle position sensor and the output from the rotation speed sensor or the acceleration sensor. It is like that. However, in multivariate systems, in general, the number of data required to establish statistical asymptotic theory depends exponentially on p, which is the number of sensors. Therefore, N → ∞ asymptotics cannot be used in many realistic situations.

  On the other hand, in the case of automobiles, power generation facilities, etc., the target to which the sensor is installed is a very complex physical system, and it is almost hopeless to establish the equation of motion of the entire system.

  For these reasons, the above prior art cannot be applied to data regarded as a multivariate system, such as automobile sensor data.

  Accordingly, the inventors, including one of the inventors of the present application, in the US patent application Ser. No. 11/668745, filed on Jan. 30, 2007, assigned to the present application, elucidates the correlation abnormality of sensor signals. Based on the correlation matrix, each sensor is embedded in the Euclidean space using the multidimensional scaling method, and a method is proposed in which the repulsive potential in the embedded space is linked to the degree of abnormality. However, this technique is not sufficient to apply the above problem in that the noise removal function is not sufficient, the embedding process tends to be unstable, and the calculated degree of abnormality is not standardized. .

  In addition, the inventors, including one of the inventors of the present application, based on the idea of proximity conservation principle in US patent application Ser. No. 11/959073 filed on Dec. 18, 2007, assigned to the present application. Therefore, we proposed an anomaly detection method intended to solve the above difficulties. The feature of this method is that it finds anomalies from broken similarity between sensor signals, does not consider each sensor other than nearby sensors, and calculates the degree of abnormality from changes in the tightness of the neighborhood graph That is.

  This technique has been used successfully for analyzing automobile data and detecting abnormalities in medium-scale power generation facilities. However, the limitation of this technique is that it does not have a means to find the vicinity of the sensor well. For example, in a variable that is independent of any other, there are no friends that can be called neighbors. On the other hand, in a sensor in which the measurement system is multiplexed, it should be possible to find friends according to the multiplicity. However, means for analyzing such differences between sensors in advance is unknown, and there is no choice but to determine the number k of neighbors ignoring the individuality of each variable. For this reason, there is a problem that false positives (not regarded as a problem but regarded as a problem) frequently occur particularly in sensors that are highly independent from others.

  Furthermore, since the definition of the degree of abnormality is based on heuristics, there is a problem that an error based on positive reversal of a single sensor cannot be detected or the detection capability is not sufficient in a noisy situation.

JP-A-6-109498 JP-A-8-278815 US patent application Ser. No. 11 / 668,745 filed Jan. 30, 2007 US patent application Ser. No. 11/959073 filed Dec. 18, 2007 J. Friedman, T. Hastie, and R. Tibahirani, Sparse inverse covariance estimation with the graphical lasso. Biostatics, 9 (3): 432-441, 2008 O. Banarjee, L.E.Ghaoui, and G. Natsoulis, Convex optimization techniques for fitting sparse gaussian graphical models.In in Proceedings of ICML, pages 89-96.Press, 2006. J. Friedman, T. Hastie, H. Hofling and Tibshinari, Pathwise coordinate optimization. Annals of Applied Statistics. I (2): 302-332

  An object of the present invention is to solve the above-mentioned problems of the prior art by providing a function for automatically inferring a proximity graph structure around each sensor from data.

In the present invention, the following are given as inputs.
(1) Test data D = {x ^(t) | x ^(t) ∈ R ^p , t = 1,2, ..., N}
(2) Normalization constant ρ Here, 0 <ρ <1, and the value is appropriately selected by the user.
(3) Reference normal data D˜ = {x ~ ^(t) | x ~ ^(t) ∈ R ^p , t = 1,2,..., N} correlation coefficient matrix S˜

As a result, the abnormalities e ₁ , e ₂ ,..., E _p of each sensor are output.

但し、l_i(A|B)というのは、negentropy(A,B)のことである。
４．次のステップは、下記の式によって、i = 1,2,...,pに対する異常度を計算することである。
e_i = max{l_i(Λ|Ｓ) - l_i(Λ〜|Ｓ),l_i(Λ|Ｓ〜) - l_i(Λ〜|Ｓ〜)}
これは、相対エントロピー、またはカルバック・ライブラー距離として知られているものである。 The outline of the algorithm according to the present invention is as follows. Note that all of the following processing is automatically executed by a computer program. First, a procedure for obtaining a correlation coefficient matrix is represented by a function display called corr (). Then, the first step is to obtain the correlation coefficient matrix S of the test data D as follows.
1. S ← corr (D)
2. The next step is the calculation of the accuracy matrix for D and D˜. The accuracy matrix is basically an inverse matrix of the correlation coefficient matrix. In particular, in the present invention, the accuracy matrix is sparse by solving an expression to which an L ₁ norm normalization term is added. Get a solution that becomes a matrix. One suitable algorithm for solving such an optimization equation is graphical lasso. If it is expressed by a function display called glasso (),
With Λ ← glasso (ρ, S), Λ˜ ← glasso (ρ, S˜), the accuracy matrix Λ of the test data D and the accuracy matrix Λ˜ of the reference normal data D˜ are obtained.
3. The next step is the calculation of the feature matrix for i = 1,2, ..., p. The feature quantity matrix is expressed as follows.

However, l _i (A | B) means negentropy (A, B).
4). The next step is to calculate the anomaly for i = 1,2, ..., p by the following formula:
e _i = max {l _i (Λ | S) −l _i (Λ˜ | S), l _i (Λ | S˜) −l _i (Λ˜ | S˜)}
This is what is known as relative entropy, or Cullback Ribler distance.

  As described above, according to the present invention, the processing is performed so that the accuracy matrix becomes a sparse matrix, and the optimization equation is solved preferably using graphic lasso, and the Cullback Live is obtained for each sensor. Since the error distance represents the degree of anomaly, the function for automatically inferring the proximity graph structure around each sensor from the data is incorporated into the anomaly detection method of time series data based on the proximity preservation principle. Is given a theoretically consistent degree of anomaly based on a direct comparison of local statistical models.

  This provides a more appropriate sensor time-series data abnormality determination technique with fewer false positives.

  Embodiments of the present invention will be described below with reference to the drawings. It should be understood that these examples are for the purpose of illustrating preferred embodiments of the invention and are not intended to limit the scope of the invention to what is shown here. Further, throughout the following drawings, the same reference numerals denote the same objects unless otherwise specified.

    Referring to FIG. 2, a block diagram of computer hardware for realizing a system configuration and processing according to an embodiment of the present invention is shown. In FIG. 2, a CPU 204, a main memory (RAM) 206, a hard disk drive (HDD) 208, a keyboard 210, a mouse 212, and a display 214 are connected to the system bus 202. The CPU 104 is preferably based on a 32-bit or 64-bit architecture, such as Intel Pentium (trademark) 4, Intel Core (trademark) 2 DUO, AMD Athlon (trademark), etc. Can be used. The main memory 204 preferably has a capacity of 512 KB or more, more preferably a capacity of 1 GB or more.

  The hard disk drive 208 stores in advance an operating system, a processing program according to the present invention, and the like, which are not shown individually. The operating system may be any compatible with the CPU 204, such as Linux (trademark), Microsoft Windows Vista, Windows XP (trademark), Windows (trademark) 2000, Apple Mac OS (trademark).

  The hard disk drive 208 further stores normal time series data for sensor reference and time series data for sensor test as separate files. A plurality of data recorded in different situations may be stored as individual files. Then, the user can detect anomalies in desired time-series data by selecting one of the stored files during the test.

  Similarly, a plurality of reference normal time-series data files may be prepared and stored in the hard disk drive 208.

  The keyboard 210 and the mouse 212 are used to operate graphic objects such as icons, taskbars, and windows displayed on the display 214 according to a graphic user interface provided by the operating system. The keyboard 210 and the mouse 212 are also used for designating a file in which test sensor time-series data for abnormality detection is stored.

  The display 214 is preferably, but not limited to, a 32-bit true color LCD monitor with a resolution of 1024 × 768 or higher. The display 214 is used for displaying a waveform of time series data and a result of abnormality detection.

  The hard disk drive 208 further stores a program for detecting an abnormality according to the present invention. This program can be written in any existing programming language such as C ++, C #, Java ™, Perl, or Ruby. When using Windows Vista, Windows XP (trademark), Windows (trademark) 2000, or the like as an operating system, it can be implemented as an application program including a GUI by using the function of the Win32 API. However, the program for performing abnormality detection according to the present invention can also be implemented as a CUI.

FIG. 3 is a diagram illustrating the relationship between data from sensors. As shown, the mechanism of anomaly detection of the present invention, when an arbitrary integer between the i from 1 to p, per data x _i from the sensor, the data x ₁ except data x _i,. .., x _j , ..., x _i-1 , x _{i + 1} , ..., x _p An algorithm is applied to make the relationship more prominent. That is, the relevance of data between weak relationships is ignored, and only the relevance between data having a considerably large relationship is left. When such a relationship is displayed in a matrix, a sparse matrix in which most components are 0 is obtained. In other words, it is an object of the present invention to obtain a solution in which the accuracy matrix, which is a matrix representing such a relationship, is a sparse matrix. FIG. 3 shows that the data x _i has a certain relationship with x _j and a strong relationship with x _{i + 1,} and other relationships represented by dotted lines are regarded as irrelevant.

  Refer to section 2.3.1 of C.M. Bishop, “Pattern Matching and Machine Learning”. Springer Verlag for details on the accuracy matrix.

  FIG. 4 is a flowchart showing an outline of processing according to the present invention. A program for executing this flowchart is stored in the hard disk drive 208, and can be loaded into the main memory 206 and executed by the operation of the operating system by a user operation.

In step 402, the correlation coefficient matrix of the test data D is calculated by the procedure corr ().
When the definition of the test data D is reproduced again, {x ^(t) | x ^(t) ∈ R ^p , t = 1, 2,..., N}. Here, the dimension number p of the time series data is about 50, for example, and may be several hundreds. N is the number of time divisions, and is generally about 100 to 1000. The appropriate time division interval is 0.1 seconds here, although it depends on the application. However, it is the number of dimensions p rather than N that affects the subsequent calculation amount. This is because p is directly the number of rows and columns of the correlation coefficient matrix. This data is stored in the hard disk drive 208 and is read into the main memory 206 by the processing program of the present invention for processing.

The correlation coefficient matrix of the test data D is calculated as follows. First, the i-th component of x ^(t), to be expressed as x _i ^(t). Here, i∈ {1, .., p}.

Therefore, the processing program calculates the average m _i and standard deviation σ _{i of} x _i ^(t) for t = 1, 2,. Then, the processing program standardizes x _i ^(t) by the following equation using the average m _i and the standard deviation σ _i .

In this way, the processing program obtains the correlation coefficient matrix S by the following equation.

Here, the normal reference data ^{D~ = {x ~ (t)} | x ~ (t) ∈R p, t = 1,2, ..., N} correlation matrix S~ of a given If it is not calculated yet, calculate it here as well.

In step 404, the processing program calculates the accuracy matrix Λ of the test data D and the accuracy matrix Λ˜ of the reference normal data D˜ using the function glasso () as follows. “glasso” means the graphical lasso algorithm.
Λ = glasso (ρ, D)
Λ ~ = glasso (ρ, D ~)

By definition, the accuracy matrix Λ is an inverse matrix of S of the correlation coefficient matrix. However, it is generally impossible to calculate the inverse matrix of S. This is because S is generally not regular. This is the case, for example, when there are variable pairs with some degree of strong correlation. It must be considered that such a situation always occurs in the time series data of the sensor system. Fortunately, even if S is regular, it is not necessary to simply calculate the inverse matrix to obtain the accuracy matrix Λ. This is because the accuracy matrix Λ is not a sparse matrix and does not satisfy the aim of the present invention of automatically selecting a neighborhood. Let's assume that the precision matrix Λ, which is a sparse matrix, has been found, as will be explained later. Then, according to the theory of graphical modeling, if the (i, j) component (i ≠ j) of the accuracy matrix Λ is 0, x _i and x _j are the conditions when all other variables are given. It is independent. That is,

The present invention makes use of this conditional independence for neighborhood search. That is, if it is determined that it is not statistically independent, it is regarded as a neighborhood, otherwise it is not a neighborhood.
For example, for x ₂ , if j is only 4 and 9 for which the accuracy matrix Λ _{2, j} ≠ 0, it can be seen that the neighborhood of x ₂ is only x ₄ and x ₉ . Therefore, optimally determining the precision matrix Λ as a sparse matrix with many matrix elements being zero is equivalent to the automatic selection of neighborhoods.

ここで、detは行列式をあらわす。 As already mentioned, there is no guarantee that the sensor data itself follows a multivariate normal distribution. Instead, the multivariate Gaussian model is used as a criterion for determining whether or not it is a neighborhood. Therefore, since it is assumed that the mean is zero, the model to be applied to the data is a p-dimensional Gaussian distribution with the mean zero and the precision matrix Λ. It is expressed as the following equation.

Here, det represents a determinant.

Here, attention is paid to a specific (i, j), for example, (i, j) = (1,2). And, as appropriate numbers other than x ₁ and x _2, for example, try to put a 0. Then, the shoulder of the exponential function is as follows.

Therefore, if Λ _1,2 = 0, there is no cross term connecting x ₁ and x ₂ and it can be seen that they become independent when all other data is stopped.

Therefore, the discussion moves on the sparsity of the accuracy matrix Λ. We leave the sparseness and try to estimate the accuracy matrix from the data with maximum likelihood. Then, the log likelihood of the test data D is as follows.

  Where tr is the diagonal sum of the matrix. Therefore, in order to obtain an accuracy matrix that maximizes the likelihood, it is sufficient to maximize {} in the above equation.

However, the accuracy matrix obtained in this way is generally not sparse. Therefore, the present invention utilizes the property that a sparse matrix can be obtained when maximum likelihood estimation is performed by adding a normalization term based on the L ₁ norm. The reason why the method of obtaining a sparse solution with L ₁ normalization term is called lasso (Least absolute shrinkage and selection operator) in the field of statistics is the reason for the name of graphical lasso used here.

ここで、|Λ_i,j|は、行列Λの(i,j)要素の絶対値をあらわす。 Therefore, the function f (Λ; S, ρ) to which the L ₁ norm is assigned is defined as follows.

Here, | Λ _{i, j} | represents the absolute value of the (i, j) element of the matrix Λ.

Then, obtaining Λ using this function f (Λ; S, ρ) results in solving the following optimization problem.

  A preferred method for solving this optimization problem is an algorithm called graphical lasso. A detailed explanation of this can be found in J. Friedman, T. Hastie, and R. Tibahirani, Sparse inverse covariance estimation with the graphical lasso. Biostatics, 9 (3): 432-441, 2008.

の存在によって、Λの多くの行列要素は、正定値性を失わない範囲で厳密にゼロになる。定数ρは、直感的に言うと、どの程度相関係数を無視するかの尺度であり、0 < ρ < 1の範囲でユーザが与える。graphical lassoの場合の好適な１つの値は、ρ = 0.3である。ρが1に近いと、得られる近接グラフは非常に小さくなり、ρ = 0だと、基本的に全てのセンサが結合した完全グラフが得られる。 Normalization term

, Many matrix elements of Λ are exactly zero within a range that does not lose positive definiteness. Intuitively speaking, the constant ρ is a measure of how much the correlation coefficient is ignored, and is given by the user in the range of 0 <ρ <1. One preferred value for graphical lasso is ρ = 0.3. When ρ is close to 1, the obtained proximity graph is very small, and when ρ = 0, a complete graph is basically obtained in which all sensors are connected.

  The graphical lasso algorithm will be described in more detail later.

  In this way, when the accuracy matrix Λ of the test data D and the accuracy matrix Λ˜ of the reference normal data D˜ are obtained in step 404, the matrix of the negative entropy (negentropy) in step 406 is calculated. Then, the process proceeds to calculation of the abnormality score in step 408. Steps 406 and 408 are sequentially performed for i = 1, 2,.

By the way, if there is no essential difference between the test data D and the reference normal data D˜, there should be no big difference in the state of the proximity graph represented by the accuracy matrices Λ and Λ˜ obtained from each. is there. That is, in view of the algorithm for obtaining the accuracy matrix, the proximity graph itself should be obtained regardless of the apparent difference of the time series data as shown in FIG. Conversely, when focusing on a variable x _i , if there is a large change in the adjacent graph around it, there is a high probability that something abnormal has occurred in that variable.

ここで、< .. >_Dという表示は、データＤでの経験分布をあらわす。
これらの式が、前述のガウス分布の式から見導かれることは、この分野の当業者に明らかであろう。これらの式に基づき、ステップ４０６で、個々のi=1,...,pにつき、負のエントロピーが計算される。 Now, since the accuracy matrix Λ of the test data D and the accuracy matrix Λ˜ of the reference normal data D˜ are obtained,
Probability model p (x) = N (x | 0, Λ ⁻¹ ) of test data D,
The probability model p ~ (x) = N (x | 0, Λ ~ ^-1 ) of the test data D ~ is obtained. Therefore, in consideration of these, in calculating the degree of abnormality of the variable x _i, all variables x _1, x ₂ other than the _{x i, ..., x i-} 1, x i + 1, ... , the conditional distribution p of x _i of when given the x _p | a _{(x i rest), p ~} | it is natural to evaluate the difference between (x _i rest). Here, because rest, all variables x ₁ except _{x i, x 2, ...,} x i-1, x i + 1, ..., is that the x _p. Therefore, for convenience of notation, if all variables x ₁ , x ₂ , ..., x _i-1 , x _{i + 1} , ..., x _p other than x _i are set to z _i ∈R ^p−1 , The function of negentropy is as follows.

Here, the display <..> _D represents the empirical distribution in the data D.
It will be apparent to those skilled in the art that these equations are derived from the aforementioned Gaussian equation. Based on these equations, at step 406, negative entropy is calculated for each i = 1,.

この結果は、見て取れるように、相対エントロピー、すなわちカルバック・ライブラーの距離である。
そして、e_i ≡ max(d₁,d₂)により、i番目のデータの異常度が求められる。なお、前述の２００７年１２月１８日出願の米国特許出願第１１／９５９０７３号明細書に記載の技法では、相違度行列Ｄの要素d_i,jが、- ln|S_i,j|で定義されていたので、相関係数の符合反転を検知することができなかったが、本発明の技法には、そのような制約はない。 Based on these results, the processing program calculates the degree of abnormality e _i for i = 1, 2,.

The result, as can be seen, is the relative entropy, i.e. the distance of the Cullback liber.
Then, the degree of abnormality of the i-th data is obtained by e _i ≡ max (d ₁ , d ₂ ). In the technique described in the aforementioned US patent application Ser. No. 11/959073 filed on Dec. 18, 2007, the element d _{i, j} of the dissimilarity matrix D is defined by −ln | S _{i, j} |. However, the sign inversion of the correlation coefficient could not be detected, but the technique of the present invention has no such limitation.

  In some experiments, it has been found that the technique of the present invention can detect the degree of abnormality with high accuracy even in a situation where the noise is large and some of the sensors are disconnected. It is considered that this is because the effects of noise and disconnection are suppressed to a considerable extent when the test data D and the reference normal data D˜ are compared as a sparse accuracy matrix.

そこで、この式をΛで偏微分する。このとき、下記の公式に留意する。

<Supplementary explanation about calculation of grahical lasso>
As described above, the problem to be solved by grahical lasso is to obtain a matrix Λ that maximizes the following equation.

Therefore, this equation is partially differentiated by Λ. At this time, pay attention to the following formula.

ここで、sign(Λ)とは、符号関数sign()を、行列の要素に適用した結果を要素とする行列である。すなわち、行列sign(Λ)の(i,j)要素は、Λ_i,j≠0のときsign(Λ_i,j)であり、Λ_i,j=0であるなら、[-1,1]の何らかの値を返すように定義される。 Then, the equation of the gradient resulting from the partial differentiation at Λ is as follows.

Here, sign (Λ) is a matrix whose elements are the results of applying the sign function sign () to the elements of the matrix. That is, the (i, j) element of the matrix sign (Λ) is sign (Λ _{i, j} ) when Λ _{i, j} ≠ 0, and [-1,1] if Λ _{i, j} = 0. Defined to return some value of.

  Λ is obtained by making the above gradient equation equivalent to the zero matrix, but because of its nature, Λ is all diagonal elements, so all diagonal elements of sign (Λ) are 1. .

From this, Σ _{i, j} = S _{i, j} + ρ and the diagonal elements can be determined exactly. Here, Σ≡Λ ⁻¹ was defined.

As described above, diagonal elements can be obtained strictly, but non-diagonal elements are not easily obtained. Therefore, a technique called block gradient method is used. For the block gradient method, see "O. Banarjee, LEGhaoui, and G. Natsoulis, Convex optimization techniques for fitting sparse gaussian graphical models. In in Proceedings of ICML, pages 89-96. Press, 2006." and "J. Friedman, T. Hastie, H. Hofling and Tibshinari, Pathwise coordinate optimization. Annals of Applied Statistics. I (2): 302-332 ".
In the block gradient method, the matrix is decomposed as follows.

  Here, L and W are (p−1) × (p−1) matrices, and λ and σ are scalars. Therefore, w and l are (p-1) -dimensional vectors.

  Here, the block gradient method solves the optimization problem by assuming that one column vector w and l are unknown and the other λ, σ, L, and W are known.

但し、

である。また、sは、相関係数行列Ｓにおける、wに対応する位置のベクトルである。この式から見て取れるように、βは、L₁正規化項付きの２次計画問題を解くことで求まる。よってβは、一般的な２次計画ソルバーを使っても解けるが、好適な解法があり、それについては、後述する。 The condition that Equation 15 becomes zero with respect to w determined in this way is equivalent to the following equation. β is a vector parameter variable.

However,

It is. Further, s is a vector at a position corresponding to w in the correlation coefficient matrix S. As can be seen from this equation, β can be obtained by solving a quadratic programming problem with an L ₁ normalization term. Therefore, although β can be solved using a general quadratic programming solver, there is a suitable solution, which will be described later.

The expression of the part of the vector at the position corresponding to w in Σ≡Λ ⁻¹ under the condition that the number 15 is equal to zero can be written as w − s − sign (l) = 0.

この式の右上のブロックから得られる l = -λＷ^-1wを使い、また、β≡Ｗ^-1wとおけば、
Ｗβ - s +ρsign(β) = 0
数１７の偏微分を実行することによって、この式に等しくなることが確認できる。 Therefore, when ΣΛ is calculated, it becomes as follows.

Using l = -λW ^-1 w obtained from the upper right block of this equation, and if β≡W ^-1 w,
Wβ-s + ρsign (β) = 0
By performing the partial differentiation of Equation 17, it can be confirmed that it is equal to this equation.

In the above explanation, the block gradient method is applied to the p-th variable, but it can be said that the steps are more generally as follows.
(1) First, for example, Σ = S + ρ1 _p is set.
(2) Starting from i = 1, the rows and columns are rearranged for the i-th variable so that the i-th row is in the p-th row and the i-th column is in the p-th column.
(3) After optimally determining l and λ, return the array of rows and columns.
(4) Finish if the difference is small from the original matrix. Otherwise, increase i by one. If i exceeds p, set i to 1.
(5) Return to step (2).

At this time, it is a step for optimally determining l and λ, and w ^T l + σλ = 1 is obtained from the lower right part of Equation 19, and from this and l = −λW ⁻¹ w, the following equation is obtained: Is obtained.

Note that W ⁻¹ does not appear explicitly in the resulting formula. In other words, when the correlation coefficient matrix has dropped in rank, generally W ⁻¹ cannot be obtained. In this sense, this calculation method that does not use an inverse matrix is advantageous in that it can cope with a correlation coefficient matrix whose rank is lowered.

この目的関数をg(β;Ｗ,ρ)と書き、β_iで偏微分すると、

が成り立つ。これに基づく最適化条件を次の２つの場合に分けて考える。 Now, the optimization problem corresponding to Equation 17 is written down.

When this objective function is written as g (β; W, ρ) and partial differentiation is performed with β _i ,

Holds. The optimization conditions based on this are considered in the following two cases.

と解ける。但し、

と置いた。数２３で、A_i - ρ < 0 であれば、もともとの条件β_i > 0と矛盾する。このときは、

のように単調であることから、定義域の左端、すなわち、β_i → 0が解となる。 (1) When β _i > 0, the condition that number 22 is 0 is that β _m (m ≠ i)

It can be solved. However,

I put it. If A _i −ρ <0 in Equation 23, it contradicts the original condition β _i > 0. At this time,

Therefore, the solution is the left end of the domain, that is, β _i → 0.

と、解ける。但し、A_i + ρ > 0であれば、もともとの条件β_i < 0と矛盾する。このときは、

のように単調であることから、定義域の左端、すなわち、β_i → 0が解となる。 (2) When β _i <0, the condition that number 22 is 0 is

It can be solved. However, if A _i + ρ> 0, it is inconsistent with the original condition β _i <0. At this time,

Therefore, the solution is the left end of the domain, that is, β _i → 0.

という置換を繰り返せばよい。添え字iが最後の変数まで行ったら、また最初に戻って、収束するまで繰り返す。 As described above, when the two cases are summarized, the optimization problem of Equation 21 is as follows.

You can repeat this substitution. When subscript i goes to the last variable, go back to the beginning and repeat until convergence.

As can be seen from the equation Σ _{i, j} = S _{i, j} + ρ, W _{i, i} is always positive, and this solution has no specificity. Further, the two graphical lasso description methods described above can be described by a simple iterative formula, and the calculation efficiency is high. In addition, as described above, even if Σ is dropped in rank, it is not necessary to use the inverse matrix for the calculation, so that the calculation can be executed stably.

  The above example was a case of detecting anomalies in data acquired from an automobile sensor, but is not limited to this, and any measurement data such as a power plant or a plant can be expressed as multidimensional time series data. Applicable to the case.

It is a figure which shows the example of the waveform of reference normal time series data and the time series data for a test. It is a block diagram of the hardware for implementing this invention. It is a figure which shows the neighborhood relationship between variables. It is a figure which shows the flowchart of the process of this invention.

Claims (18)

A system for calculating the degree of abnormality of multidimensional time series data by a computer,
Means for inputting first multidimensional time-series data and second multi-dimensional time-series data having the same dimension as the first multi-dimensional time-series data;
A first accuracy matrix that is an inverse matrix of the covariance matrix of the first multidimensional time series data, and a second accuracy matrix that is an inverse matrix of the covariance matrix of the second multidimensional time series data. Means for calculating each precision matrix to be a sparse matrix;
Means for calculating a negative entropy for each dimension of the multidimensional time-series data by a multivariate probability distribution model based on the first accuracy matrix and a multivariate probability distribution model based on a second accuracy matrix;
Means for calculating the degree of anomaly for each dimension based on the relative entropy distance based on the calculated negative entropy;
Anomaly calculation system for multidimensional time series data.   The system according to claim 1, wherein the first multi-dimensional time series data is reference normal data, and the second multi-dimensional time series data is test data for detecting an abnormality.   The system of claim 1, wherein the multivariate probability distribution model is a multivariate Gaussian model. The system according to claim 1, wherein the means for calculating the accuracy matrix to be a sparse matrix performs calculation by optimizing an expression to which a normalization term of L ₁ norm is added.   The system according to claim 4, wherein the means for calculating the precision matrix to be a sparse matrix uses a graphical lasso algorithm.   The system of claim 1, wherein the multi-dimensional time series data is data obtained from a plurality of sensors of an automobile, and the multi-dimensional individual elements correspond to individual sensors. A method for calculating the degree of abnormality of multidimensional time series data by a computer,
Inputting first multidimensional time-series data and second multi-dimensional time-series data having the same dimension as the first multi-dimensional time-series data;
A first accuracy matrix that is an inverse matrix of the covariance matrix of the first multidimensional time series data, and a second accuracy matrix that is an inverse matrix of the covariance matrix of the second multidimensional time series data. Calculating each precision matrix to be a sparse matrix;
Calculating a negative entropy for each dimension of the multidimensional time-series data by a multivariate probability distribution model based on the first accuracy matrix and a multivariate probability distribution model based on a second accuracy matrix;
Calculating the degree of abnormality for each dimension based on the relative entropy distance based on the calculated negative entropy.
Abnormality calculation method for multidimensional time series data.   8. The method according to claim 7, wherein the first multi-dimensional time series data is reference normal data, and the second multi-dimensional time series data is test data for detecting an abnormality.   The method of claim 7, wherein the multivariate probability distribution model is a multivariate Gaussian model. 8. The method according to claim 7, wherein the step of calculating the accuracy matrix so as to be a sparse matrix is performed by optimizing an expression to which an L ₁ norm normalization term is added.   The system of claim 10, wherein the step of calculating the precision matrix to be a sparse matrix uses a graphical lasso algorithm.   The method of claim 7, wherein the multi-dimensional time series data is data obtained from a plurality of sensors of an automobile, and the multi-dimensional individual elements correspond to individual sensors. A program for calculating the degree of abnormality of multidimensional time series data by a computer,
In the computer,
Inputting first multidimensional time-series data and second multi-dimensional time-series data having the same dimension as the first multi-dimensional time-series data;
A first accuracy matrix that is an inverse matrix of the covariance matrix of the first multidimensional time series data, and a second accuracy matrix that is an inverse matrix of the covariance matrix of the second multidimensional time series data. Calculating each precision matrix to be a sparse matrix;
Calculating a negative entropy for each dimension of the multidimensional time-series data by a multivariate probability distribution model based on the first accuracy matrix and a multivariate probability distribution model based on a second accuracy matrix;
Calculating the degree of abnormality for each dimension based on the relative entropy distance based on the calculated negative entropy;
program.   14. The program according to claim 13, wherein the first multi-dimensional time series data is reference normal data, and the second multi-dimensional time series data is test data whose abnormality is to be detected.   The program according to claim 13, wherein the multivariate probability distribution model is a multivariate Gaussian model. 14. The program according to claim 13, wherein the step of calculating the precision matrix to be a sparse matrix is performed by optimizing an expression to which an L ₁ norm normalization term is added.   The program according to claim 16, wherein the step of calculating the precision matrix to be a sparse matrix uses a graphical lasso algorithm.   14. The program according to claim 13, wherein the multidimensional time series data is data acquired from a plurality of sensors of an automobile, and the individual elements of the multidimensional correspond to individual sensors.

Images