LEARNING SPARSITY-CONSTRAINED GAUSSIAN GRAPHICAL MODELS IN ANOMALY DETECTION

Abstract

A first dependency graph is constructed based on a first data set by solving an objective function constrained with a maximum number of non-zeros and formulated with a regularization term comprising a quadratic penalty to control sparsity. The quadratic penalty in constructing the second dependency graph is determined as a function of the first data set. A second dependency graph is constructed based on a second data set by solving the objective function constrained with the maximum number of non-zeros and formulated with the regularization term comprising a quadratic penalty. The quadratic penalty in constructing the second dependency graph is determined as a function of the first data set and the second data set. An anomaly score is determined for each of a plurality of sensors based on comparing the first dependency graph and the second dependency graph, nodes of which represent sensors.

Description

STATEMENT REGARDING PRIOR DISCLOSURES BY THE INVENTOR OR A JOINT INVENTOR

The following disclosure is submitted under 35 U.S.C. 102(b)(1)(A): D. T. Phan, T. Ide, J. Kalagnanam, M. Menickelly and K. Scheinberg, “A Novel l0-Constrained Gaussian Graphical Model for Anomaly Localization,” 2017 IEEE International Conference on Data Mining Workshops (ICDMW), New Orleans, La., USA, November 2017, pp. 830-833.

BACKGROUND

The present disclosure relates to computer-implemented anomaly detection and learning of dependency graphs via which anomalies are detected, for example, in industries such as manufacturing and heavy industries.

In industrial or manufacturing processes, or other like processes, a network of sensors is deployed to enable automatic monitoring of various operating conditions. Data mining techniques, for example, can convert low-level sensor data collected from Internet-of-Things (IoT) devices into actionable insights, for example, for condition-based maintenance. For instance, identifying an unusual occurrence or change in data can signal anomaly.

In traditional problem settings, a goal of anomaly detection is to compute or determine the degree of anomalousness for a multivariate measurement, giving an overall anomaly score. In sensor networks, for example, wireless sensor networks, it may be inadequate simply to indicate whether or not a network is behaving anomalously, for example, to perform anomaly detection. For example, if a sensor network takes measurements from different car parts of a prototype automobile on a test track, then it may be also valuable to the field engineers to know which car parts are contributing to an anomalous behavior, for example, in addition to knowing that the car as a whole is behaving anomalously. Similar problem settings can be found across different application domains, e.g., monitoring mechanical stresses in buildings and bridges, identifying the sources of flooding, and finding sources of threats in computer networks.

BRIEF SUMMARY

A method and system, for example, which detects anomaly in sensor networks, may be provided. The method, in one aspect, may include receiving a first data set comprising sensor data detected by a plurality of sensors coupled to equipments in operation during a first period of time. The method may also include constructing a first dependency graph based on the first data set by solving a first objective function constrained with a maximum number of non-zeros and formulated with a regularization term comprising a quadratic penalty to control sparsity, wherein the quadratic penalty is determined as a function of the first data set in constructing the first dependency graph. The method may further include receiving a second data set comprising sensor data detected by the plurality of sensors coupled to the equipments in operation during a second period of time. The method may also include constructing a second dependency graph based on the second data set by solving the first objective function constrained with the maximum number of non-zeros and formulated with the regularization term comprising a quadratic penalty, wherein the quadratic penalty is determined as a function of the first data set and the second data set in constructing the second dependency graph. The method may also include determining an anomaly score for each of the plurality of sensors based on comparing the first dependency graph and the second dependency graph. In a further aspect, the method may also include, responsive to determining that an anomaly score associated with a sensor in the plurality of sensors meets a threshold value indicative of abnormal functioning, automatically performing a corrective action to correct equipment coupled to the sensor.

A system, in one aspect, may include a hardware processor coupled with a communication interface and a memory device coupled to the hardware processor. The hardware processor may be operable to receive via the communication interface, a first data set comprising sensor data detected by a plurality of sensors coupled to equipments in operation during a first period of time. The hardware processor may be further operable to construct a first dependency graph based on the first data set by solving a first objective function constrained with a maximum number of non-zeros and formulated with a regularization term comprising a quadratic penalty to control sparsity, wherein the quadratic penalty is determined as a function of the first data set in constructing the first dependency graph. The hardware processor may be further operable to receive a second data set comprising sensor data detected by the plurality of sensors coupled to the equipments in operation during a second period of time. The hardware processor may be further operable to construct a second dependency graph based on the second data set by solving the first objective function constrained with the maximum number of non-zeros and formulated with the regularization term comprising a quadratic penalty, wherein the quadratic penalty is determined as a function of the first data set and the second data set in constructing the second dependency graph. The hardware processor may be further operable to determine an anomaly score for each of the plurality of sensors based on comparing the first dependency graph and the second dependency graph. In a further aspect, the hardware processor may be further operable to, responsive to determining that an anomaly score associated with a sensor in the plurality of sensors meets a threshold value indicative of abnormal functioning, automatically perform a corrective action to correct equipment coupled to the sensor.

A computer readable storage medium storing a program of instructions executable by a machine to perform one or more methods described herein also may be provided.

Further features as well as the structure and operation of various embodiments are described in detail below with reference to the accompanying drawings. In the drawings, like reference numbers indicate identical or functionally similar elements.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a diagram illustrating anomaly detection for multivariate noisy sensor data in one embodiment.

FIG. 2 is a diagram illustrating computing of a change score of a sensor in a dependency graph in one embodiment.

FIG. 3 is a diagram showing components of a system in one embodiment, which performs localization for multivariate time series datasets detects anomalies in sensor data.

FIG. 4 is a diagram illustrating a screen display of an analytics tool in an embodiment, which integrate anomaly detection technique in one embodiment.

FIG. 5 illustrates a schematic of an example computer or processing system that may implement a system in one embodiment of the present disclosure.

DETAILED DESCRIPTION

In embodiments, methods, systems and techniques are provided, which detect anomalies in sensor networks. In embodiments, methods, systems and techniques may anomaly localization, where a variable-wise anomaly score is desired. In embodiments, methods, systems and techniques include learning a dependency graph of multivariate data, which may be noisy. The methods, system and techniques automatically identify the number of dependencies and perform transfer learning from a training phase to a testing phase. Learned dependency graphs are used to detect anomalies in sensor networks. In one aspect, a learning performed according to an embodiment in the present disclosure may be useful in tracking inter-variable dependency, reducing the number of runs to obtain a graph that is sparse enough for interpretability (e.g., in detecting anomalies), while capturing real dependencies in the multivariate data.

In an embodiment, methods, systems and techniques may be implemented as a computer-implemented solution library, for instance, which can be integrated or plugged into a computer-implemented platform or server which may function as analytics tool, for example, in performing operations managements in an industry.

In an embodiment, an optimization is generated, which allows to automatically identify the number of conditional dependencies for a pair of variables and transfer a learning sparsity pattern and magnitude of correlations from a normal state to a testing phase. In an embodiment, a two-step algorithm based on projections and active set identification is presented. In an embodiment, a full gradient may be used for solution quality and a reduced subspace Newton direction to accelerate it; Feasibility with respect to membership in sparsity constraint is handled via projection; Symmetric positive-definiteness is ensured through a line-search; a continuation step improves solution quality and running time.

Anomaly detection, for example, via learned sparse dependency graph, in one embodiment identifies asset behavior that is outside of a defined norm that constitutes normal behavior. In an embodiment, anomaly detection algorithm is based on a probabilistic graphical model, and may detect anomalies with higher accuracy and fewer false alarms. In an embodiment, early detection of anomalies or faults in a system can lead to considerable savings of time and cost. A system and/or method of the present disclosure, for example, reduces maintenance cost for example, incurred in vessel monitoring, increases asset utilization, prevents unexpected equipment breakdowns by identifying anomalous conditions that can lead to accidents or equipment damage, and reduces false positive and false negative alerts. In industrial processing, for example, productivity of assets may be enhanced, and downtime and maintenance costs may be reduced.

A system and/or method, in an embodiment, solved a sparsity-constrained model, allows transfer learning via an algorithm. The algorithm, in one aspect, learns and transfers connectivity pattern and dependency magnitude from the normal training phase to the testing phase for the dependency matrix. In an embodiment, an additional regularization term may be included, which keeps the magnitudes of precision matrix entries to manageable amount. In an embodiment, a system and/or method selects the sparsity level of the dependency graphs a priori.

FIG. 1 is a diagram illustrating anomaly detection for multivariate noisy sensor data in one embodiment. A model building phase 102 in an embodiment includes using past dataset associated with sensors, for example, operating in a domain (e.g., industry domain), for example, in normal condition (e.g., without failure). Such dataset is considered to be normal data. For instance, dataset 104, for example, of past period (e.g., selected historical time period, e.g., a number of time series sample data points (e.g., 100 sample data points) associated with each sensor), may be obtained and a dependency graph model 106 is built based on the normal data. The dependency graph model includes nodes which represent sensors and links between the sensors during operating conditions.

During a runtime phase 108, a current dependency graph 112 from real-time sensor data 110 (data points associated with current period of time, or window of time) is built or computed. The dependency graph 106 associated with the normal state is compared with the current dependency graph 112. A change in the dependency among or between the nodes indicates anomaly in the current processing in the domain. In this way, anomalies among sensors are detected in real-world situations to predict when and where maintenance is required.

For example, at 114, a dependency graph 106 may be discovered or built based on the past sensor data 104 taken during normal operations of a process (e.g., manufacturing or industrial process). At 116, a current dependency graph 112 is computed from real-time sensor data 108. The dependency graphs 106 and 112 may be compared to identify or determine any anomaly in the current operations of the process.

As an example application, nodes in a dependency graph represent sensors, for example, monitoring different components of a process such as a manufacturing or industrial process. Dependency arcs or edges between nodes represent that the data of the sensors move together, or have a pattern of moving together. If, for example, two nodes are conditionally independent, then the value of dependency edge is zero. In one embodiment, detecting anomaly is carried out by comparing the change of dependency graphs between the training phase and testing phase.

FIG. 2 is a diagram illustrating computing of a change score of a sensor in a dependency graph in one embodiment. A preparation phase 202, in an embodiment, includes obtaining or receiving data matrix X, which includes data such as sensor data output during an operation of a manufacturing or industrial process, for example, by a system coupled with sensors for operating the manufacturing or industrial process. For instance, for each sensor, a time series of data points (e.g., of a selected time window) associated with a system operating under normal operating conditions (e.g., no errors, faults or failures of equipment or process), may be obtained.

At 206, the data of data matrix X is input to an optimization model, Equation (a), and Equation (a) is solved subject to a specified constraint.

In an embodiment, the constraint specifies a maximum number of non-zero entries. Non-zero entries indicate a relationship between sensors. The partial correlation between two sensors, for example, is non-zero when the corresponding element of dependency matrix is non-zero. At 208, a solution to Equation (a) creates a dependency graph G (e.g., a Gaussian graphical model).

An operation phase 210, for example, at time t, in an embodiment, includes obtaining or receiving data matrix X^(t) 212, for example, sensor data output during a manufacturing or industrial process at time window (t), which may be the current time or near current time. At 214, the data of data matrix X^(t) 212 is input to the optimization model, Equation (a), and Equation (a) is solved subject to a specified constraint. At 216, the solution of the optimization model (Equation (a)) creates a dependency graph G_t, (e.g., a Gaussian graphical model) associated with the current data (operations data of time t window).

In another aspect, at 214, Equation (c) may be solved using both data matrix X 210 and data matrix X^(t) 212 to create the dependency graph G_t, for example, in lieu of Equation (a).

The dependency graph G 208 and the dependency graph G_t 216 are input to another model, Equation (b), and Equation (b) is solved at 218. Solving Equation (b) produces a change score a_i (anomaly score for sensor i) at time t 220, for each sensor considered. For example, an anomaly score for each sensor is determined, for example, by solving Equation (b).

$\begin{array}{cc}\underset{X\succ 0}{\min }〈S,X〉-\log \det \left(X\right)+\frac{\rho }{2}{X}_Fs.t.{X}_0\le \kappa ,s.t.X_{i,j}=0,\mathrm{if}x_i\mathrm{and}x_j\mathrm{are}\mathrm{conditionally}\mathrm{independent}& \mathrm{Equation}\left(a\right)\\ a_i=\int {\mathrm{dx}}_{-i}p\left(x_{-i}D\right)\int {\mathrm{dx}}_ip\left(x_ix_{-i},D\right)\ln \left(\frac{p\left(x_ix_{-i},D\right)}{p\left(x_ix_{-i},D^{\prime }\right)}\right)& \mathrm{Equation}\left(b\right)\end{array}$

where D is the training data set and D′ is the testing data set. A detailed formula is given in Equation (8) below.

$\begin{array}{cc}\underset{X\succ 0}{\min }〈S,X〉-\log \det \left(X\right)+\frac{\rho }{2}{X-X_{\mathrm{train}}}_Fs.t.{X}_0\le \kappa ,s.t.X_{i,j}=0,\mathrm{if}x_i\mathrm{and}x_j\mathrm{are}\mathrm{conditionally}\mathrm{independent}& \mathrm{Equation}\left(c\right)\end{array}$

In the above formulations, S represents the sample covariance of the user-inputted data, X represents the dependency graph to be learned at time t (also referred to as X^(t) above), and X_train represents the dependency graph from the normal operation phase (i.e., the training phase), ρ represents a regularization parameter, ∥⋅∥_F represents a norm term. Details of the formulations are described further below.

In an embodiment, the L0 (l₀)-constraint (∥X∥₀≤K) guarantees that the solution admits a level of sparsity. In an embodiment, the quadratic penalty

$\left(\frac{\rho }{2}{X-X_{\mathrm{train}}}_F\right)$

encourages the capacity of selecting groups in the presence of highly correlated variables and transfer learning. The models of Equations (a) and (b) are non-convex models, which may be challenging to solve, and in embodiments, optimization algorithms are provided to solve such non-convex models.

In an embodiment, the problem of anomaly localization in a sensor network for multivariate time-series data is considered by computing anomaly scores for each variable separately. A system and/or method in an embodiment evaluate a difference measure using the conditional expected Kullback-Liebler divergence between two sparse Gaussian graphical models (GGMs) learned from different sliding windows of the dataset. To estimate the sparse GGMs, an optimization model is provided and the system and/or method in an embodiment constrain the number of nonzeros (l₀ norm) in the learned inverse covariance matrix and apply an additional l₂-regularization in the objective. In an embodiment, a projected gradient algorithm efficiently solves this nonconvex problem. Numerical evidence supports the benefits of using the model and method over the usual convex relaxations (i.e., l₁-regularized models) for learning sparse GGMs in this anomaly localization setting.

Anomaly localization of the present disclosure in an embodiment determines or computes a variable-wise anomaly score. In wireless sensor networks, e.g., in addition to determining whether or not a network is behaving anomalously, it is useful to determine or know what aspect of the network is anomalous. For example, if a sensor network takes measurements from different car parts of a prototype automobile on a test track, then it is valuable to the field engineers to know which car parts are contributing to an anomalous behavior, e.g., more than simply knowing that the car as a whole is behaving anomalously. Similar problem settings can be found across different application domains: monitoring mechanical stresses in buildings and bridges, identifying the sources of flooding or forest fires, pinpointing the change points in stock price time series, finding sources of serious threats in computer networks, and identifying car loss.

In an embodiment, sparse dependency graphs (e.g., GGMs) are learned in the context of anomaly localization. In an embodiment, a model is presented to learn a sparse dependecy graph where both sparsity pattern and numeric values are taken into account. GGM-learning model employing both an l₀ constraint and an l₂-regularization, along with an optimization algorithm to solve it, outperforms other models including typical l₁-regularized models. In one aspect, the system and/or method leverage the conditional expected Kullback-Liebler (KL) divergence method in computing anomaly scores for each variable. A comparative study conducted on various GGM-learning models and scoring methods for assigning anomaly scores demonstrates efficiency of a system and method in an embodiment. Using l₀ norm-based methods can present computationally challenging issues. The system and/or method in an embodiment provide l₀-constrained GGM learning algorithms in the anomaly detection setting. In an embodiment, an l₀-constrained l₂-regularized model learns a sparse dependency graph for anomaly localization, and a proximal gradient algorithm handles both an l₀ constraint and a positive-definite constraint.

In an embodiment, the following notations are used. For any X∈R^n×nm ∇f(X) denotes the gradient of f, which is a matrix in R^n×n. Notation tr(X) denotes the trace operator, i.e., Σ_i=1ⁿX_i,i. The vectorization operator vec(X) transforms X into a vector in Rⁿ² by stacking the columns from left to right. σ₁(X) and σ_n (X) denote the minimum and maximum eigenvalues of X, respectively. Given an index set F⊂{1, . . . , n}×{1, . . . , n}, [X]_F denotes the matrix obtained from the identity matrix I_n by replacing the (i, j)-th entry by X_i,j for all (i, j)∈F. Let ⊗ denote the Kronecker product of matrices and let e_i denote a unit vector with a 1 in the i th coordinate. For any finite set S, |S| denotes the number of elements of S.

Learning Sparse Graphical Models

The following describes optimization models that the system and/or method use in an embodiment to learn sparse precision matrices. Given m data samples {y₁, . . . y_m}⊆Rⁿ, the system in an embodiment considers fitting an n-variate Gaussian distribution N(μ, X⁻¹) to the data, where μ denotes the mean and X is the precision matrix. In an embodiment, it is assumed the data is centered so that μ=0. Zeros in X indicate conditional independence between two variables. To estimate the unknown parameter X, a maximum likelihood problem as follows may be solved,

$\begin{array}{cc}\underset{X\succ 0}{\min }\mathrm{tr}\left(SX\right)-\log \det \left(X\right),& \left(1\right)\end{array}$

where X0 denotes a positive definite constraint. In an embodiment, the empirical covariance matrix is defined as

$S=\frac{1}{m}\sum _{i=1}^my_iy_i^T.$

The solution X∈R^n×n is the estimated inverse covariance, which describes a Gaussian graphical model (GGM). However, in general, due to noise in the data, X obtained from Equation (1) is a dense matrix and thus the conditional dependencies suggested by X may reveal very little. For reasons of interpretability, sparsity in the matrix X is desired.

The sparse graphs for anomaly localization have been computed by solving the well-studied l₁-regularized maximum likelihood problem, i.e.,

$\begin{array}{cc}\underset{X\succ 0}{\min }\mathrm{tr}\left(SX\right)-\log \det \left(X\right)+\lambda {X}_1,& \left(2\right)\end{array}$

where ∥X∥₁=Σ_i,j=1ⁿ|X| is the element-wise l₁ norm of the matrix X, and λ>0 is a regularization parameter to control sparsity. Equation (2) is a convex problem, and generally admits a unique global minimizer. Thus, given S, selecting a sparse graph is a matter of computing a regularization path by varying λ and cross-validating the solutions along the path. Since sensor data of industrial applications are very noisy in general, some useful interaction information among the variables in the graph can be omitted.

Enforcing sparsity in the constraint set instead of in the objective function may avoid bias. In an embodiment, the system and/or method of the present present disclose and solve a model that estimates sparse graphs. In an embodiment, the system and/or method directly constrains sparsity by specifying a maximally allowable number of nonzeros κ in the optimal solution and adds a quadratic penalty

$\begin{array}{cc}\underset{X\succ 0}{\min }\mathrm{tr}\left(SX\right)-\log \det \left(X\right)+\frac{\lambda }{2}{X}_F^2s.t.{X}_0\le \kappa ,& \left(3\right)\end{array}$

where λ>0 is a regularization parameter and the Frobenius norm term ∥X∥_F²=Σ_i,j=1ⁿ|X_ij|² is also referred to as l₂-regularization. The present disclosure denotes by ∥⋅∥₀ the number of nonzeros of the argument, the so-called l₀-norm. The l₀-constraint guarantees that the solution admits a level of sparsity, while the quadratic penalty encourages the capacity of selecting groups in the presence of highly correlated variables. Without the regularization, some entries of the purely l₀-constrained model can have relatively large magnitudes. The associated anomaly scores significantly dominate the others, and thus some faulty variables can be overlooked. Hence, in an embodiment, the l₂-regularization term keeps the magnitude of all entries uniformly similar. From a computational point of view, the l₂-term makes the problem of Equation (3) better-behaved. The objective function is strongly convex and the optimal solution set is bounded. Thus, it is expected to be solved more efficiently. In an embodiment, the l₀-constraint and the l₂-regularization are combined. Because of the presence of the l₀ constraint, it may be challenging to solve problem of Equation (3); in an embodiment, an efficient proximal gradient algorithm is presented to handle the problem.

In an embodiment, the model of Equation (3) can be used to learn dependency graphs in both the train phase and the testing phase, for example, learning a sparse GGM for anomaly detection. In an embodiment, the system and/or method may also use the following for the testing phase

$\begin{array}{cc}\underset{X\succ 0}{\min }\mathrm{tr}\left(SX\right)-\log \det \left(X\right)+\frac{\lambda }{2}{X-X_{\mathrm{normal}}}_F^2s.t.{X}_0\le \kappa ,& \left(4\right)\end{array}$

where X_normal is a precision matrix from a normal state. The system and/or method in one embodiment perform transfer learning for sparsity pattern and magnitude of correlations from the normal state to the testing phase.

$\begin{array}{cc}\begin{array}{cc}\underset{X\succ 0}{\min }& \mathrm{tr}\left(\mathrm{SX}\right)-\log \det \left(X\right)\\ s.t.& {X}_0\le \kappa .\end{array}& \left(5\right)\\ \underset{X>0}{\min }\mathrm{tr}\left(\mathrm{SX}\right)-\log \det \left(X\right)+{\lambda }_1{X}_1+{\lambda }_2{X}_F^2.& \left(6\right)\end{array}$

Equations (5) and (6) represent additional models, which may be applied in learning a sparse GGM for anomaly detection algorithms, for instance, in the context of anomaly localization. Equation (6) is an an l₁-based “elastic net” version for GGM model.

Methods for Computing Anomaly Score of Each Variable

The following describes embodiments of techniques that perform a change analysis, e.g., assigning anomaly scores to individual variables to measure the magnitude of their contributions to the observed differences between two GGMs from training and testing datasets. The techniques, for example, in embodiments, implement the processing shown at 120 in FIG. 1.

Conditional Expected KL-Divergence

The system and/or method in an embodiment may consider the conditional expected Kullback-Liebler divergence between two learned GGMs to compute anomaly scores. Given two probabilistic models A and B with respective distributions p_A and p_B, the system and/or method consider the expected KL divergence d_i^AB between the marginal distributions p_A(x_i|z_i) and p_B(x_i|z_i)

$\begin{array}{cc}{d}_i^{\mathrm{AB}}=\int {\mathrm{dz}}_ip_A\left(z_i\right)\int {\mathrm{dx}}_ip_A\left(x_iz_i\right)\ln \frac{p_A\left(x_iz_i\right)}{p_B\left(x_iz_i\right)},& \left(7\right)\end{array}$

for each variable i. By swapping A and B in Equation (7), the system and/or method obtain a symmetric formula for d_i^BA. Because the system and/or method are learning GGMs, the formula of Equation (7) takes a special form that reduces to matrix-vector multiplications. For the i-th variable, and for a learned GGM X, let

$\begin{array}{cc}X=\left[\begin{array}{cc}L& I\\ I^T& \alpha \end{array}\right],X^{-1}=\left[\begin{array}{cc}W& w\\ w^T& \beta \end{array}\right]& \left(7b\right)\end{array}$

be permuted such that the last row and column of X, X⁻¹ correspond to the i-th variable. Then, letting X^A and X^B denote the GGMs learned from datasets A and B respectively, it can be shown

$\begin{array}{cc}{d}_i^{\mathrm{AB}}=w^{A^T}\left(I^B-I^A\right)+\frac{1}{2}\left[\frac{I^{B^T}W^AI^B}{{\alpha }^B}-\frac{I^{A^T}W^AI^A}{{\alpha }^A}\right]+\frac{1}{2}\left[\ln \frac{{\alpha }^A}{{\alpha }^B}+{\beta }^A\left({\alpha }^B-{\alpha }^A\right)\right].& \left(8\right)\end{array}$

The system and/or method then define the anomaly score of the i th variable as

d _i=max{d_i^AB,d_i^BA}.  (9)

When a threshold t>0 is selected, d_i>t is an indication that the i-th variable is anomalous, while d_i≤t is an indication that the i th variable is healthy. Alpha and beta, L and W notations are components of matrices X and X⁻¹ when they are are broken up as in equation (7b).

Stochastic Nearest Neighbors

The system and/or method in an embodiment may consider using a stochastic nearest neighbors approach to anomaly scoring. Let S^A and S^B denote the sample covariance matrices of datasets A and B respectively, and let S_i^A, S_i^B denote the i th columns of S^A and S^B respectively. For a sparse GGM X^A learned for dataset A, define an indicator vector associated with the i th variable l_A,i coordinate-wise by

$\begin{array}{cc}{\left[l_{A,i}\right]}_j=\{\begin{array}{cc}1& \mathrm{if}i\mathrm{and}j\mathrm{are}\mathrm{adjacent}\mathrm{in}X^A\\ 0& \mathrm{otherwise}.\end{array}& \left(10\right)\end{array}$

Then, a measure of the dissimilarity between the neighborhoods of the i th variable between the GGMs for A and B, weighted by the sample covariances of the two datasets, is given by

$\begin{array}{cc}{d}_i^{\mathrm{AB}}=\frac{l_{A,i}^T\left(S_i^A-S_i^B\right)}{\left(1+l_{A,i}^TS_i^A\right)\left(1+l_{A,i}^TS_i^B\right)}& \left(11\right)\end{array}$

Symmetrically, the system and/or method define d_i^BA and then compute an anomaly score d_i as in Equation (9). “T” in the equations stands for the transpose of a vector or matrix.

Sparsest Subgraph Approximation

The system and/or method in an embodiment may consider using sparsest subgraph approximation technique to score anomaly. Given two (sparse) dependency graphs, e.g. X^A and X^B consider a graph given by an (also sparse) adjacency matrix A defined entrywise by Λ_ij=|X_ij^A−X_ij^B|. The problem of finding an independent set (a set of mutually disconnected nodes) of size k in Λ is generalized to finding a “sparsest k-subgraph” in Λ:

$\begin{array}{cc}\underset{d\in {\left\{0,1\right\}}^n}{\min }d^T\Lambda d\mathrm{subject}\mathrm{to}1_n^Td=k,& \left(12\right)\end{array}$

where 1_n denotes a vector of ones. Intuitively, the solution to Equation (12) identifies a set of k nodes that change the least between X^A and X^B, suggesting a healthy (non-anomalous) set of nodes. However, (12) is a mixed-integer quadratic programming (MIQP) formulation, the solution of which is generally NP-hard, and moreover requires foreknowledge of the magnitude of k, the size of the supposedly healthy set of nodes. Therefore, the system and/or method may use a convex QP relaxation (solvable in polynomial time)

$\begin{array}{cc}\underset{d\in R^n}{\min }d^T{\Lambda }_{\mu }ds.\mathrm{to}1_n^Td=1,d\ge 0_n,& \left(13\right)\end{array}$

where 0_n is a vector of zeros and Λ_μ=Λ+μI_n is the original matrix Λ with an added scaled identity with μ≥0 sufficiently bounded away from zero to enforce the positive definiteness of Λ_μ. With this relaxation, the solution d* to Equation (13) can be interpreted as anomaly scores.

In an embodiment, the system and/or method uses the l₀-l₂ based model of Equations (3) and (4) to learn a sparse GGM in the training phase and the testing phase respectively, and computes anomaly scores for each variable by applying the conditional expected KL divergence of Equation (9).

Optimization Algorithm for l₀ Sparse Models

In embodiments, the system and/or method implements a projected gradient algorithm for solving the learning problems or models of Equation (3) and Equation (4).

Projected Gradient Newton Algorithm

In an embodiment, a projected gradient Newton algorithm with warm starts may solve the l₀-constrained inverse covariance learning problems of Equations (3) and (4). The following describes a method for solving Equation (3), and the method also can be applied to Equation (4). It is noted when A is small, Equation (3) is difficult to be solved. In an embodiment, instead of directly solving (3), the system and/or method solves a sequence of subproblems parameterized by a bigger λ>0,

$\begin{array}{cc}\begin{array}{cc}\underset{X\succ 0}{\min }& f\left(X\right)=\mathrm{tr}\left(\mathrm{SX}\right)-\mathrm{logdet}\left(X\right)+\frac{\lambda }{2}{X}_F^2\\ s.t.& {X}_0\le \kappa \\ & X_{i,j}=0\forall \left(i,j\right)\in I,\end{array}& \left(14\right)\end{array}$

where the system and/or method denote the Frobenius norm by ∥X∥_F²=Σ_i,j=1ⁿ|X_ij|². The objective of Equation (14) is strongly convex and the optimal solution set of Equation (14) is bounded. The l₂ regularizer makes the subproblems better-conditioned, possibly allowing for finding higher quality local minima of Equation (14) faster. By driving λ to zero in Equation (14), an optimal solution for Equation (3) can be obtained. Thus, to yield a practical algorithm, the system and/or method in an embodiment consider a sequence of strictly decreasing values {λ_k}; in the k-th iteration of an algorithm, the system and/or method seek an approximate local minimizer to Equation (14) with λ=λ_k. On the (k+1)-th iteration, the system and/or method warm start the solution of the (k+1)-th subproblem with the solution to the k-th subproblem as an initial point.

Before analyzing optimality conditions of (3) or (14), a hardness result may be proved, for example, demonstrating that for any κ>0 in Equation (3), there exists no λ>0 such that an optimal solution to the l₁ regularized model of Equation (2) with regularization parameter λ recovers the global optimum of Equation (3). This suggests solving Equation (3) without resorting to convex formulations.

Theorem 1. For any κ≥n, there exists no value of λ>0 such that an optimal solution of (2) globally minimizes (3).

Proof 1. Towards a contradiction, suppose there exists κ≥n and λ>0 such that X is a global minimizer of both (2) and (3). Define Y=X⁻¹ and Ψ(X)=−log det(X)+tr(SX). Since S is positive semi-definite and nonzero, there exists at least one diagonal entry S_ii>0. Without loss of generality, suppose S₁₁>0. For any α such that |α|<min(X₁₁, σ₁(X)), both ∥X+αe₁e₁^T∥₀=∥X∥₀≤κ and X+αe₁e₁^T are positive definite. Additionally, since (X+αe₁e₁^T)_kl=0 for any (k,l)∈| it may be concluded that X+αe₁e₁^T is feasible in Equation (3) for every α such that |α|<min(X₁₁, σ₁(X)).

Moreover, using a property of determinants,

$\begin{array}{cc}\psi \left(X+\alpha e_1e_1^T\right)-\psi \left(X\right)=& -\log \left(\det \left(X\right)\left(1+\alpha Y_{11}\right)\right)+\\ & \mathrm{tr}\left(S\left(X+\alpha e_1e_1^T\right)\right)+\\ & \mathrm{logdet}\left(X\right)-\mathrm{tr}\left(\mathrm{SX}\right)\\ =& -\log \left(1+\alpha Y_{11}\right)+\alpha S_{11}\underset{\_}{\underset{\_}{\Delta }}g\left(\alpha \right).\end{array}$

Note that if X is a minimizer of Equation (3), then g(α)≥0 for every |α|<min(X₁₁,σ₁(X)). Since g is a convex one-dimensional function of α in the neighborhood of zero and g(0)=0, it may be concluded that g′(0)=0. This implies that S₁₁=Y₁₁. From a characterization of optimal solutions of Equation (2), if X is optimal for Equation (2), then S₁₁=Y₁₁−λ. Since λ>0, this is a contradiction.

First-Order Stationarity Conditions

The following characterizes properties of Equation (14), for example, concerning the boundedness of the set of optimal solutions and appropriate optimality conditions for Equation (3). In an embodiment, it may be assumed that the l₂-regularization parameter λ>0. For ease of notation, the sparsity constraint set is defined as

ωΔ{X∈R^n×n:∥X∥₀≤κ,X_i,j=0∀(i,j)∈|},

and the projection operator

$P_{\Omega }\left(X\right)\underset{\_}{\underset{\_}{\Delta }}\arg \underset{Y\in \Omega }{\min }{X-Y}_F.$

It may be noted that the gradient of the objective function in Equation (14) is

∇f(X)=S−X⁻¹+λX,

and the Hessian is

∇²f(X)=X⁻¹⊗X⁻¹+λI_n².

It is shown in the following theorem that Equation (14) can be safely reformulated in such a way that the feasible set is compact.

Theorem 2. Assume that problem of Equation (14) is feasible and λ>0. Consider the problem

$\begin{array}{cc}\begin{array}{cc}\underset{X\succcurlyeq 0}{\min }& \mathrm{tr}\left(\mathrm{SX}\right)-\mathrm{logdet}\left(X\right)+\frac{\lambda }{2}{X}_F^2\\ s.t.& {X}_0\le \kappa \\ & X_{i,j}=0,\forall \left(i,j\right)\in I\\ & \underset{\_}{\sigma }I_nX\stackrel{\_}{\sigma }I_n,\end{array}& \left(15\right)\end{array}$

where there are fixed constants

$\stackrel{\_}{\sigma }=\frac{n^2+\sqrt{n^4+2\lambda \mathrm{nf}\left({\mathrm{tI}}_n\right)-2\lambda n^2}}{\lambda },\underset{\_}{\sigma }={\exp }^{-f\left({\mathrm{tI}}_n\right)+\frac{n\lambda }{2}{\exp }^{-2f\left({\mathrm{tI}}_n\right)}{\stackrel{\_}{\sigma }}^{2-2n}}{\stackrel{\_}{\sigma }}^{1-n},$

and

$t=\frac{-\mathrm{tr}\left(S\right)+\sqrt{{\mathrm{tr}\left(S\right)}^2+4n^2\lambda }}{2n\lambda }.$

Then, the sets of global minima of Equation (14) and Equation (15) are equivalent.

Proof 2. Suppose X*∈R^n×n is an optimal solution of Equation (14) and 0<σ₁≤ . . . ≤σ_n are the eigenvalues of X*. The system and/or method in an embodiment seek lower and upper bounds for σ₁ and σ_n, respectively.

Consider X=tI_n for an arbitrary fixed scalar t>0. X is feasible with respect to the constraints of Equation (14), and

$\begin{array}{cc}F\left(\stackrel{\_}{X}\right)\ge f\left(X^{*}\right)=\mathrm{tr}\left({\mathrm{SX}}^{*}\right)-\mathrm{logdet}\left(X^{*}\right)+\frac{\lambda }{2}{X}_F^2\ge -\mathrm{logdet}\left(X^{*}\right)+\frac{\lambda }{2}\sum _{i=1}^nX_{\mathrm{ii}}^{{*}^2}\ge -n\log \left({\sigma }_n\right)+\frac{\lambda }{2n}{\left(\sum _{i=1}^nX_{\mathrm{ii}}^{*}\right)}^2\ge -n\log \left({\sigma }_n\right)+\frac{\lambda }{2n}{\sigma }_n^2.& \left(16\right)\end{array}$

The second inequality is due to the fact that S±0 and X*0, and thus tr(SX*)≥0. The third inequality uses the Cauchy-Schwarz inequality, and the last inequality holds because Σ_i=1ⁿX_ii*=Σ_i=1ⁿσ_i and σ_i>0. Combining Equation (16) with the fact that log(σ_n)≤σ_n−1, the system and/or method get

$\begin{array}{cc}f\left(\stackrel{\_}{X}\right)\ge -n\left({\sigma }_n-1\right)+\frac{\lambda }{2n}{\sigma }_n^2.& \left(17\right)\end{array}$

The right-hand side of Equation (17) tends to infinity as σ_n→∞; since f(X) is finite for the fixed value of t, σ_n is bounded from above. Specifically, σ_n cannot exceed the larger of the two solutions to the quadratic equation

$\frac{\lambda }{2n}{\sigma }_n^2-n\left({\sigma }_n-1\right)-f\left(\stackrel{\_}{X}\right)=0.$

That is,

$\begin{array}{cc}{\sigma }_n\le \stackrel{\_}{\sigma }\underset{\_}{\underset{\_}{\Delta }}\frac{n^2+\sqrt{n^4+2\lambda \mathrm{nf}\left(\stackrel{\_}{X}\right)-2\lambda n^2}}{\lambda }<\infty .& \left(18\right)\end{array}$

Likewise, a lower bound for σ₁ is estimated as

f(X)≥f(X*)≥−log det(X*)

≥−log(σ₁)−(n−1)log(σ),

from which it may be concluded

σ₁≥exp(−f(X))σ^1-n.  (19)

In an embodiment, this bound can be tightened further. Starting from an intermediate step in the derivation of (16), there is

$f\left(\stackrel{\_}{X}\right)\ge -\mathrm{logdet}\left(X^{*}\right)+\frac{\lambda }{2n}{\left(\sum _{i=1}^n{\sigma }_i\right)}^2\ge -\log \left({\sigma }_1\right)-\left(n-1\right)\log \left(\stackrel{\_}{\sigma }\right)+\frac{n\lambda }{2}{\sigma }_1^2\ge -\log \left({\sigma }_1\right)-\left(n-1\right)\log \left(\stackrel{\_}{\sigma }\right)+\frac{n\lambda }{2}\exp \left(-2f\left(\stackrel{\_}{X}\right)\right){\stackrel{\_}{\sigma }}^{2-2n}.$

Thus,

$\begin{array}{cc}{\sigma }_1\ge \underset{\_}{\sigma }\underset{\_}{\underset{\_}{\Delta }}{\exp }^{-f\left(\stackrel{\_}{X}\right)+\frac{n\lambda }{2}{\exp }^{-2f\left(\stackrel{\_}{X}\right)}{\stackrel{\_}{\sigma }}^{2-2n}}{\stackrel{\_}{\sigma }}^{1-n}>0.& \left(20\right)\end{array}$

From Equations (18) and (20), in order to tighten the bounds σ and σ, the system and/or method may choose t so that f(X)=f(tI_n) is minimized. Observe that

$f\left({\mathrm{tI}}_n\right)=\mathrm{tr}\left(S\right)t-n\log \left(t\right)+\frac{n\lambda }{2}t^2$

is a strongly convex one-dimensional function in t for t>0. Thus f (tI_n) is minimized when the necessary optimality condition

$\frac{d}{\mathrm{dt}}f\left({\mathrm{tI}}_n\right)=0$

is satisfied; thus the following is chosen

$t=\frac{-\mathrm{tr}\left(S\right)+\sqrt{{\mathrm{tr}\left(S\right)}^2+4n^2\lambda }}{2n\lambda },$

which completes the proof.

The following introduces in Theorem 3 optimality conditions for problem (14). For any positive definite matrix H and a scalar t>0, define

$\begin{array}{cc}{h}_{t,H,X}\left(Z\right)=\mathrm{tr}\left(\nabla {f\left(X\right)}^T\left(Z-X\right)\right)+\frac{1}{2t}{Z-X}_H^2,& \left(21\right)\\ S_{t,H}\left(X\right)=\arg \underset{Z\in \Omega }{\min }h_{t,H,X}\left(Z\right)=\arg \underset{Z\in \Omega }{\min }P{Z-\left(X-t〈H,\nabla f\left(X\right)〉\right)}_{H^{-1}}^2,\mathrm{where}{X}_H=\sqrt{{\mathrm{vec}\left(X\right)}^T\mathrm{Hvec}\left(X\right)}.& \left(22\right)\end{array}$

Theorem 3. Suppose that X* is an optimal solution of (14) and H is positive definite. Then there exists T>0 such that

X*=S _t,H(X*) for all t∈(0,T).  (23)

Proof 3. From Theorem 2, there is, any optimal solution X* of (14) satisfies σI_nX*, where σ>0 is independent of X*. Since ∇²f(X*)=X*⁻¹⊗X*⁻¹+λI_n2, it follows that ∇²f(X*)(σ⁻²+λ)I_n2. Hence, there is, from strong convexity that

$\begin{array}{cc}f\left(Y\right)\le f\left(X^{*}\right)+\mathrm{tr}\left(〈\nabla f\left(X^{*}\right),\left(Y-X^{*}\right)〉\right)+\frac{1}{2L_f}{Y-X^{*}}_F^2,& \left(24\right)\end{array}$

for all Y∈Ω, where L_f=(σ⁻²+λ)⁻¹. It is shown that (23) holds for T=σ₁(H)L_f. First, it is demonstrated that X*∈S_t,H (X*) for any t∈[0, T). This first claim is a statement concerning existence of X*. The following subsequently proves uniqueness of X*, which proves the theorem. From (24) and since

$\frac{1}{t}H\succ \frac{1}{L_f}I_n$

for all t∈[0, T), there is, for all Y∈Ω

$f\left(X^{*}\right)\le f\left(Y\right)\le f\left(X^{*}\right)+\mathrm{tr}\left(\nabla f\left(X^{*}\right)\left(Y-X^{*}\right)\right)+\frac{1}{2t}{Y-X^{*}}_H^2.$

It follows that

$h_{t,H,X^{*}}\left(Y\right)=\mathrm{tr}\left(\nabla f\left(X^{*}\right)\left(Y-X^{*}\right)\right)+\frac{1}{2t}{Y-X^{*}}_H^2\ge 0$

for all Y∈Ω. Thus X*∈S_t, (X*) since h_t,H,X*(X*)=0, proving the claim.

To prove uniqueness, suppose towards contradiction that there exists {circumflex over (X)}∈S_t,H(X*) satisfying (23), but {circumflex over (X)}≠X*. Noting that by definition h_t,H,X*(X*)=0, we must have h_t,H,X*({circumflex over (X)})=0, and so

$\begin{array}{cc}\mathrm{tr}\left(\nabla f\left(X^{*}\right)\left(\hat{X}-X^{*}\right)\right)=-\frac{1}{2t}{\hat{X}-X^{*}}_H^2.& \left(25\right)\end{array}$

Combining (24) and (25), there is,

$\begin{array}{cc}f\left(\hat{X}\right)\le f\left(X^{*}\right)+\mathrm{tr}\left(\nabla f\left(X^{*}\right)\left(\hat{X}-X^{*}\right)\right)+\frac{1}{2L_f}{\hat{X}-X^{*}}_F^2=f\left(X^{*}\right)-\frac{1}{2}{\left(\hat{X}-X^{*}\right)}^T\left(\frac{H}{t}-\frac{I_n}{L_f}\right)\left(\hat{X}-X^{*}\right)<f\left(X^{*}\right),& \left(26\right)\end{array}$

which contradicts the global optimality of X* for (14).

For a general positive definite matrix H, solving the problem defining S_t,H (X) in (22) is generally intractable due to the cardinality constraint. As a result, verifying the optimality condition (23) is computationally expensive. The following lemma offers a simpler condition, which is an immediate consequence of Theorem 3 when H=I_n. A similar result was established, but for a slightly different constraint set Ω.

Lemma 1. Suppose that X* is an optimal solution of (14). Then there exists T>0 such that

X*=P _Ω(X*−t∇f(X*)) for all t∈(0,T).  (27)

In one aspect, a subsequent convergence analysis relies on Lemma 1. The following states and proves a technical proposition regarding the operator S_t,In, concerning how it acts on convergent sequences.

Proposition 1. Let t>0. Suppose {X^k} is a sequence in R^n×n converging to X*, and suppose there exists a corresponding sequence {Y^k} such that Y^k∈S_t,In (X^k) for each k. If there exists Y* such that Y^k→Y*, then Y*∈S_t,In(X*).

Proof 4. Since Y^k∈Ω for each k and Ω is closed, it is concluded that Y*∈Ω. Hence, it is obtained from the definition of S_t,In(X*) that

$\begin{array}{cc}{Y^{*}-\left(X^{*}-t\nabla f\left(X^{*}\right)\right)}_F\ge \underset{Y\in \Omega }{\min }{Y-\left(X^{*}-t\nabla f\left(X^{*}\right)\right)}_F& \left(28\right)\end{array}$

If inequality (28) holds as an equality, then the proposition is proven. Towards a contradiction, suppose (28) holds as a strict inequality. Then, for all k,

$\begin{array}{cc}\begin{array}{cc}{Y^{*}-\left(X^{*}-t\nabla f\left(X^{*}\right)\right)}_F>& \mathrm{ }\underset{Y\in \Omega }{\min }{Y-\left(X^{*}-t\nabla f\left(X^{*}\right)\right)}_F\\ \ge & \mathrm{ }{\underset{Y\in \Omega }{\min }Y-\left(X^k-t\nabla f\left(X^k\right)\right)}_F-\\ & \left(X^{*}-t\nabla f\left(X^{*}\right)\right)-\left(X^k-t\nabla f\left(X^k\right)\right){}_F\\ =& \mathrm{ }{Y^k-\left(X^k-t\nabla f\left(X^k\right)\right)}_F-\\ & \left(X^{*}-t\nabla f\left(X^{*}\right)\right)-\left(X^k-t\nabla f\left(X^k\right)\right){}_F,\end{array}& \left(29\right)\end{array}$

where the second inequality follows from the triangle inequality, and the equality follows from the definition of Y^k. As k→∞, (29) implies

∥Y*−(X*−t∇f(X*))∥_F>∥Y*−(X*−t∇f(X*))∥_F,

a contradiction.

Algorithm Description

An algorithm is presented to solve a general problem of the form

min{f(X):X∈R^n×n,∥X∥₀≤κ,X_i,j=0,∀(i,j)∈|,X0,X=X^T},  (30)

where f is twice differentiable on its domain. In an embodiment, a strategy for solving the constrained optimization problem is to embed approximate Newton method into a projected gradient method. In this algorithm, feasibility with respect to sparsity constraint Ω is handled via projection, while the symmetric positive-definiteness of the iterates (i.e. feasibility with respect to {X0}) in ensured through a line-search procedure. Convergence analysis can be applied, for example, when the objective function has the special structure as in (14):

$f\left(X\right)=\mathrm{tr}\left(\mathrm{SX}\right)-\mathrm{logdet}\left(X\right)+\frac{\lambda }{2}{X}_F^2.$

A particular issue in the case of applying a projected gradient method to (30) is that the constraint set of (30) is nonconvex and not closed. This is a different situation from problems in the conventional scope of projected gradient methods, which require the feasible set to be convex and closed. In one aspect, due to the openness of the positive-definite (PD) constraint, directly applying a projected gradient method to Equation (30) is difficult. In an embodiment, a methodology of the present disclosure avoids dealing with the PD constraint in the projection step.

In an algorithm of the present disclosure in an embodiment, the projected gradient direction and the approximate Newton step are alternatively employed. On each gradient search step, the algorithm begins from a feasible iterate X^k and then backtrack along the projection are defined by P_Ω(X^k−α_k∇f(X^k)), initializing the stepsize α_k>0 with a Barzilai-Borwein (BB) stepsize [1]:

${\alpha }_k^{\mathrm{BB}}=\frac{\mathrm{tr}\left({\left(\nabla f\left(X^k\right)-\nabla f\left(X^{k-1}\right)\right)}^T\left(X^k-X^{k-1}\right)\right)}{{X^k-X^{k-1}}_F^2}.$

BB-stepsizes perform better when they are embedded in a nonmonotone line search. In an embodiment, a scheme based on the Grippo-Lampariello-Lucidi (GLL) stepsize rule is considered. Let f_max^k denote the largest of the M most recent function values:

f _max ^k=max{f(X^k−j):0≤j<min(k,M)}.

This quantity f_max^k is used to define a sufficient decrease condition. For instance, the algorithm terminates the backtracking procedure for computing X^k+1=P_Ω(X^k−α_k∇f(X^k)) by gradually reducing α_k once the following conditions (C1) and (C2) are both satisfied:

(C1) A sufficient decrease condition for the objective function value is attained, i.e.

$f\left(X^{k+1}\right)\le f_{\max }^k-\frac{\delta }{2}{X^{k+1}-X^k}_F^2,$

for some algorithmic parameter δ>0.

(C2) The next iterate X^k+1 is feasible with respect to the positive definite constraint, i.e. X^k+10.

In an embodiment, to accelerate convergence and make use of the fact that the optimal solution is often very sparse in many applications, i.e., κ=n², the algorithm in the present disclosure may additionally incorporate projected approximate Newton steps, by using second-order information in a reduced subspace. At each iterate X^k, a local quadratic model of the objective function can be minimized subject to the constraint that all but a working set W_k of the matrix entries will not be fixed to zero in the minimization. In an embodiment of the algorithm, W_k may be chosen as

W _k={(i,j):X_i,j^k≠0}∪O_k  (31)

or

W _k={(i,j):Y_i,j≠0 for some Y∈P_Ω(X^k−α_k∇f(X^k))}∪O_k,  (32)

where α_k>0 is a step size and O_k is a set of indices corresponding to the ┌ηκ┐ largest absolute values of entries in X^k−X^k−1 for some (small) η∈(0,1). The inclusion of O_k helps to close the optimality gap faster when the condition (27) in Lemma 1 does not hold. The determination for updating the working set W_k from (31) or (32) is based on how X^k is computed, which is described in detail. In one aspect, the particular selection of W_k is not critical to proving theoretical convergence results, and the prescriptions in (31), (32) are intended as practical guidance.

For any working set W_k⊆{1, . . . , n}×{1, . . . , n} and a matrix Y, define Y_wk coordinate-wise via

${\left(Y_{W_k}\right)}_{\left(i,j\right)}=\{\begin{array}{cc}{Y}_{\left(i,j\right)}& \mathrm{if}\left(i,j\right)\in W_k\\ 0& \mathrm{otherwise}.\end{array}$

Then, given a current iterate X^k, gradient ∇f(X^k), and Hessian ∇²f(X^k), solve the following reduced quadratic problem to obtain a Newton direction D^k:

$\begin{array}{cc}\underset{D\in R^{n\times n}}{\min }q\left(D\right)\equiv \mathrm{tr}\left(\nabla {f\left(X^k\right)}_{W_k}D_{W_k}\right)+\frac{1}{2}{\mathrm{vec}\left(D_{W_k}\right)}^T{\nabla }^2{f\left(X^k\right)}_{W_k}\mathrm{vec}\left(D_{W_k}\right).& \left(33\right)\end{array}$

To efficiently obtain an approximate solution of (33), an embodiment of a methodology of the present disclosure makes use of a conjugate gradient method. This conjugate gradient method takes advantage of the particular Kronecker product structure of ∇²f(X^k), as well as the observation that |W_k|<<n², for instance, since the methodolgoy searches for a very sparse solution in applications.

A detailed description of our algorithm for solving (30) is given in Algorithm 1.

Algorithm 1: Reduced-space Projected Gradient Algorithm -
RPGA(X0, Ω, λ)
1  Choose parameters
   σϵ (0,1),δ > 0,ηϵ (0,1),αminNewton > 0,[αmin,αmax] ⊂ (0,∞) ,
integers p,q such that 0 < p < q , integerM > 0, grad_step = true .
2  Set k = 0.
3  while some stopping criteria not satisfied do
4  Step 1: Active-set Newton step
5  if (mod(k,q) ≠ p) then
6  if (grad_step = true) then
7  Wk ←Eq. (31)
8  else
9  Wk ←Eq. (32)
   /* end of if-then-else of line 6 */
10 Approximately solve for Dk from (33)
11 Choose αk = αk,0 = 1
12 j ← 0, ls ← true
13 while (ls = true) & (αk ≥ αminNewton) do
14 αk ← σjαk,0
15 Xtrialk+1 ← PΩ (Xk + αkDk)
16 if(f(Xk+1)≤fmaxk-δ2Xk+1-XkF2andXk+1≻0)
   then
17 Xk+1 ← Xtrialk+1
18 ls ← false
19 else
20 j ← j + 1
   /* end of if-then-else of line 16 */
   /* end of while loop of line 13 */
   /* end of if-then of line 5 */
21 Step 2: Gradient projection step
22 grad_step ← false
23 if (mod (k,q) = p) or (αk < αminNewton) then
24 Choose αk,0 = min(αmax,max(αmin,αkBB))
25 j ← 0, ls ← true
26 while (ls = true) do
27 αk ← σjαk,0
28 Xtrialk+1 ← PΩ(Xk − αk∇f(Xk)
29 if(f(Xk+1)≤fmaxk-δ2Xk+1-XkF2andXk+1≻0)
   then
30 Xk+1 ← Xtrialk+1
31 ls ← false
32 else
33 j ← j + 1 /* end of if-then-else of line 29 */
   /* end of while loop of line 26 */
34 grad_step ← true
   /* end of if statement of line 23 */
35 Step 3: Iterate
36 k ← k + 1
   /* end of whiile loop of line 3 */

Remark 1. As mentioned, Ω is a non-convex set, and as such, the operator P_Ω is generally set-valued, i.e. projections are non-unique. A point in P_Ω(X) can be quickly obtained. If X∈R^n×n, set X_i,j=0 for every (i, j)∈. Then P_δ(X) is computed by sorting the n² many entries of |X|, where |⋅| denotes an entry-wise absolute value, and then replacing all but the κ largest values of |X| with 0, where ties could be broken arbitrarily. In order to attain a unique projection, when we compute P_Ω, ties are not broken arbitrarily, but by a dictionary order on the matrix entries. The dictionary order is chosen in such a way to ensure that the projected matrix is symmetric. As a consequence of sorting n² matrix entries, the projection operation can be performed in O(n² log(n)) time.

On each iteration of the backtracking procedure, one would compute an eigenvalue decomposition of the trial step X^k+1 both in order to compute the value of det(X^k+1) in the objective function and to determine the satisfaction of (C2). Computing a determinant is generally an expensive operation. However, since a method in the present disclosure is explicitly searching for a positive definite matrix, it is enough to attempt a Cholesky factorization of X^k+1, which can be done more efficiently. If the factorization fails (a complex number appears on the diagonal of the lower triangular matrix), then the current X^k+1 is not positive definite and the method can immediately retract along the projection arc. If the factorization is successful, then the determinant of the matrix X^k+1 can be computed from the diagonal values of lower triangular matrix. This sequence of attempted Cholesky factorizations is repeated until the stopping criteria is satisfied. The finite termination in the backtracking step is proved below.

In an embodiment, performing a line search in the projected gradient direction in the full space may help a projected approximate Newton algorithm identify a good working set subspace, and enables proving of global convergence. Thus, in an embodiment, the general logic of Algorithm 1 is that, for fixed q∈N, after every q−1 iterations of Step 1 (the active-set Newton steps), the algorithm performs a single iteration of Step 2, the gradient projection step.

As in the case in nonconvex optimization, the iterates of Algorithm 1 need not converge to a global minimizer of (14), but only to a point satisfying necessary conditions for optimality. The difference between (14), which Algorithm 1 is to solve, and (3), the real problem of interest, lies in the l₂ regularization term. In an embodiment, a method of the present disclosure uses an outer homotopy method for solving (3). Algorithm 2 solves a sequence of subproblems of the form (14) by calling Algorithm 1, attempting to find increasingly better-quality minimizers. Two strictly decreasing sets of integer values κ₀>κ₁> . . . >κ_K=κ and real values λ₀>λ₁> . . . >λ_K=A are selected before running the algorithm. A method of the present disclosure may first (approximately) solve (14) using a sparsity parameter κ₀, regularization parameter λ₀, and some initial point X⁰. Then, the method may sequentially solve (14) to some tolerance for each pair (κ_j,λ_j): j=1, . . . , K, initializing each subproblem solve with the previous solution X^j-1. The final regularization parameter λ_K should be quite small in practice.

Algorithm 2: Projected Gradient Newton Algorithm with Warm
Starts - PPNA
1 Choose homotopy sequence {(κj, λj): j = 0,1,2,..., K}, initial point X0.
2 for j = 0,1,...,K do
3   Ωj ← {X:||X||0 ≤ κj, Xi,j = 0 ∀(i, j) ϵ ℑ}
4   Xj+1 ← RPGA(Xj ,Ωj,λj)
/* end of for loop of line 2 */

Convergence Analysis

This section provides a proof that the iterates of Algorithm 1 for solving problem (14) accumulate at a point satisfying the necessary optimality conditions (27). It may be assumed that the initial point X⁰ is feasible with respect to the constraints of (14).

The following states an assumption on κ:

Assumption 1. Assume that κ in (14) (and (15)) is an integer satisfying κ∈[n, n²). Observe that this assumption is reasonable, since as demonstrated earlier, n is a sufficient condition for (15) to have a non-empty feasible set, and if κ=n², then the problem reduces to a convex one. It is demonstrated that the line search in either Step 1 (the active-set Newton step) or Step 2 (the gradient projection step) terminates in a finite number of iterations.

Lemma 2. Suppose that problem (14) satisfies Assumption 1. Let {X^k} denote the sequence of iterates generated by Algorithm 1 for solving problem (14). Suppose the line search starting in either Line 13 or Line 26 of Algorithm 1 is reached in the k-th iteration and suppose X^k is feasible with respect to the constraints of (14). Then, the line search in steps 1 and 2 terminates in a finite number of iterations.

Proof 5. The line search in Step 1 terminates in a finite number of iterations because the termination condition α_k<α_min^Newton is enforced. The following proves the finite termination when the line search occurs in Step 2 (the gradient projection step). Assume ∇f(X^k)≠0, since otherwise the algorithm stops a local solution X^k.

There exists T_p^k>0 such that every matrix in P_Ω(X^k−t∇f(X^k)) is positive definite for all t∈(0, T_p^k). For a nonzero symmetric matrix U, denote the spectral norm of U by ∥U∥₂Δ max(|σ_min(U)|, |σ_max(U)|)>0, where σ_min(U),σ_max(U) denote the smallest and largest eigenvalues of U, respectively. Thus, for all t∈(0,L₁), there is X^k+tU0, where there is defined

$L_1=\frac{{\sigma }_1\left(X^k\right)}{{U}_2}.$

Given X^k, let I and A denote respectively the inactive and active matrix indices of X^k with respect to the cardinality constraint, i.e.:

I Δ{(i,j): X_i,j^k≠0}

A Δ{(i,j): X_i,j^k=0}.

Additionally define

$\begin{array}{ccc}\left(i_x,j_x\right)& \underset{\_}{\underset{\_}{\Delta }}& \arg \underset{\left(i,j\right)\in I}{\min }X_{i,j}^k\\ \left(i_{g,n},j_{g,n}\right)& \underset{\_}{\underset{\_}{\Delta }}& \arg \underset{\left(i,j\right)\in I}{\max }{\left(\nabla f\left(X^k\right)\right)}_{i,j}\\ \left(i_{g,a},j_{g,a}\right)& \underset{\_}{\underset{\_}{\Delta }}& \arg \underset{\left(i,j\right)\in A}{\max }{\left(\nabla f\left(X^k\right)\right)}_{i,j}.\end{array}$

The following demonstrates the existence of the necessary T_p^k>0 by analyzing two separate cases concerning the cardinality of the set of inactive indices, |I|.

Case 1.): For every (i, j)∈I there is

|X_i,j^k−t(∇f(X^k))_i,j|≥|X_i,j^k|−t|(∇f(X^k))_i,j|≥|X_ix,jx^k|−t|(∇f(X^k))_ig,n,jg,n|.

Defining

$T_1^k=\frac{X_{i_x,j_x}^k}{{\left(\nabla f\left(X^k\right)\right)}_{i_{g,n},j_{g,n}}+{\left(\nabla f\left(X^k\right)\right)}_{i_{g,a},j_{g,a}}},$

it is concluded from (34) that for all (i, j)∈I,

|X_i,j^k−t(∇f(X^k))_i,j|>t|(∇f(X^k))_ig,a,jg,a|≥t|(∇f(X^k))_h,i|,  (34)

for all (h,l)∈A and t∈(0, T₁^k). Note that T₁^k is well-defined since |X_ix,jx^k| and |(∇f(X^k))_ig,n,jg,n|+|(∇f(X^k))_ig,a,jg,a| are strictly positive. It is possible to decompose

X ^k −t∇f(X^k)=[X^k−t∇f(X^k)]_I+[X^k−t∇f(X^k)]_A=[X^k−t∇f(X^k)]_I−t[∇f(X^k)]_A.

Combining this decomposition with Remark 1 and (34), it follows that

[X^k−t∇f(X^k)]_I=P_Ω(X^kt∇f(X^k)) for t∈(0,T₁^k).

Now consider U=[∇f(X^k)]_I, and define

$T_2^k\underset{\_}{\underset{\_}{\Delta }}\{\begin{array}{cc}\frac{{\sigma }_1\left(X^k\right)}{{U}_2}& \mathrm{if}U\ne 0\\ 1& \mathrm{otherwise}\end{array}$

Then, letting

T _p ^k=min(T₁^k,T₂^k),

not only is X^k−t[∇f(X^k)]_I=P_Ω(X^k−t∇f(X^k)), but also X^k−t[∇f(X^k)]_I is positive definite for any t∈(0,T_p^k).

Case 2. (|I|<κ): Denote by Â⊆A the set of indices (i, j)∈A corresponding to the κ−|I| largest absolute magnitude entries of [∇f(X^k)]_A.

Denote

$\left(i_{g,a^{\prime }},j_{g,a^{\prime }}\right)\underset{\_}{\underset{\_}{\Delta }}\arg \underset{\left(i,j\right)\in A\backslash \hat{A}}{\max }{\left(\nabla f\left(X^k\right)\right)}_{i,j}.$

Define

$T_3^k=\frac{X_{i_x,j_x}^k}{{\left(\nabla f\left(X^k\right)\right)}_{i_{g,n},j_{g,n}}+{\left(\nabla f\left(X^k\right)\right)}_{i_{g,a^{\prime }},j_{g,a^{\prime }}}}.$

From (34), it follows that

|X_i,j^k+t(∇f(X^k))_i,j|>t|(∇f(X^k))_ig,a,jg,a|≥t|(∇f(X^k))_k,h|,  (35)

for all (i, j)∈I,(k, h)∈A\Â and t∈(0, T₃^k). Similarly to Case 1, it is possible to decompose

$\begin{array}{cc}{X}^k-t\nabla f\left(X^k\right)=& {\left[X^k-t\nabla f\left(X^k\right)\right]}_I+{\left[X^k-t\nabla f\left(X^k\right)\right]}_{\hat{A}}+\\ & {\left[X^k-t\nabla f\left(X^k\right)\right]}_{A\backslash \hat{A}}\\ =& {\left[X^k-t\nabla {f\left(X^k\right)}^k\right]}_I-{t\left[\nabla f\left(X^k\right)\right]}_{\hat{A}}-{t\left[\nabla f\left(X^k\right)\right]}_{A\backslash \hat{A}}.\end{array}$

From this decomposition, Remark 1, and (35), there is [X^k−t∇f(X^k)]_I∪Â∈P_Ω(X^k−t∇f(X^k)) for t∈(0, T₃^k). Define U=[∇f(X^k)]_I∪Â, and analogous to the definition of T₂^k in Case 1, define

$T_4^k\underset{\_}{\underset{\_}{\Delta }}\{\begin{array}{cc}\frac{{\sigma }_1\left(X^k\right)}{{U}_2}& \mathrm{if}U\ne 0\\ 1& \mathrm{otherwise}\end{array}$

Choose

T _p ^k=min(T₃^k,T₄^k).

It is concluded that for any t∈(0,T_p^k), there is both X^k−t[∇f(X^k)]_I∪Â∈P_Ω(X^k−t∇f(X^k)) and X^k−t[∇f(X^k)]_I∪Â is positive definite. Since every matrix in P_Ω(X^k−t∇f(X^k)^k) can be represented as X^k−t[∇f(X^k)]_I∪Â, the existence of the needed T_p^k is proven.

The following proves that there exists a positive number T_δ>0, independent of k, such that

$\begin{array}{cc}f\left(X^{k+1}\right)\le f_{\max }^k-\frac{\delta }{2}{X^{k+1}-X^k}_F^2\mathrm{provided}t\in \left(0,T_{\delta }\right),& \left(36\right)\end{array}$

where X^k+1=P_Ω(X^k−t∇f(X^k)) and suppose that

$f\left(X^k\right)\le f_{\max }^{k-1}-\frac{\delta }{2}{X^k-X^{k-1}}^2.$

Here, the notion X⁻¹=X⁰ and f_max⁻¹=f(X⁰) is adopted. Define

L={X∈R ^n×n :f(X)≤f(X⁰)}.

Note that f(X^k)≤f_max^k−1≤f(X⁰), thus X^k∈L. Replacing X* and X by X and X⁰ respectively in (18) and (20) yields L∪{X∈R^n×n:σI_nXσI_n}, where σ and σ are some constants satisfying 0<σ<σ<∞. Define

K={X∈R ^n×n:0.5σI_nX(σ+0.5σ)I_n}.  (37)

There is L⊂K. Since X^k+1 minimizes h_t,In,xk(Z), it follows

$\begin{array}{cc}{h}_{t,I_n,X^k}\left(X^{k+1}\right)=& \mathrm{tr}\left(\nabla {f\left(X^k\right)}^T\left(X^{k+1}-X^k\right)\right)+\frac{1}{2t}{X^{k+1}-X^k}_F^2\\ \le & \mathrm{ }h_{t,I_nX^k}\left(X^k\right)=0.\end{array}$

Obtain

∥X^k+1−X^k∥_F²≤2ttr(∇f(X^k)^T(X^k−X^k+1))

Taking norms yields

∥X^k+1−X^k∥_F≤2t∥∇f(X^k)∥_F=2t∥S−(X^k)⁻¹+λX^k∥_F

Since X^k lies in the compact set K, there exists a constant c, independent of X^k∈K, such that

∥∇f(X^k)∥_F=∥S−(X^k)⁻¹+λX^k∥_F≤c.

Choose β∈(0,∞) so that 2δc≤0.5σ. X^k+1∈K when t≤β and X^k∈L . Indeed, there is

$\begin{array}{cc}{\sigma }_{\min }\left(X^{k+1}\right)\ge & \mathrm{ }{\sigma }_{\min }\left(X^k\right)-{\sigma }_{\max }\left(X^{k+1}-X_k\right)\\ =& \mathrm{ }{\sigma }_{\min }\left(X^k\right)-{X^{k+1}-X^k}_2\\ \ge & \mathrm{ }\underset{\_}{\sigma }-\frac{{X^{k+1}-X^k}_F}{\sqrt{n}}\\ \ge & \underset{\_}{\sigma }-\frac{\underset{\_}{\sigma }}{2\sqrt{n}}>0.5\underset{\_}{\sigma }\left(\mathrm{since}{X^{k+1}-X^k}_F\le 0.5\underset{\_}{\sigma }\right)\end{array}$

Furthermore, there is

$\begin{array}{cc}{\sigma }_{\max }\left(X^{k+1}\right)\le & {\sigma }_{\max }\left(X^k\right)+{\sigma }_{\max }\left(X^{k+1}-X^k\right)\\ \le & \mathrm{ }\stackrel{\_}{\sigma }+{X^{k+1}-X^k}_F\\ \le & \mathrm{ }\stackrel{\_}{\sigma }+0.5\underset{\_}{\sigma }.\end{array}$

Hence it is concluded that X^k+1∈K. Because the eigenvalues of any matrix in K are bounded, f is Lipschitz continuously differentiable on K. Denote L_f by the Lipschitz constant for f on K. In view of the convexity of K, the line segment connecting X^k and X^k+1 is contained in K. There is

$\begin{array}{cc}f\left(X^{k+1}\right)\le & f\left(X^k\right)+\mathrm{tr}\left(\nabla {f\left(X^k\right)}^T\left(X^{k+1}-X^k\right)\right)+\frac{1}{2L_f}{X^{k+1}-X^k}_F^2\\ \le & \mathrm{ }f\left(X^k\right)-\left(\frac{1}{2t}-\frac{1}{2L_f}\right){X^{k+1}-X^k}_F^2\\ \le & \mathrm{ }f_{\max }^k-\left(\frac{1}{2t}-\frac{1}{2L_f}\right){X^{k+1}-X^k}_F^2\\ \le & \mathrm{ }f_{\max }^k-\frac{\delta }{2}{X^{k+1}-X^k}_F^2\mathrm{if}\delta \le \frac{1}{t}-\frac{1}{L_f}.\end{array}$

Define

$T_{\sigma }=\min \left(\frac{\underset{\_}{\sigma }}{4c},\frac{L_f}{\delta L_f+1}\right),$

then there is

$\begin{array}{cc}f\left(X^{k+1}\right)\le f_{\max }^{k-1}-\frac{\delta }{2}{X^{k+1}-X^k}_F^2\mathrm{provided}t\in \left[0,T_{\delta }\right),& \left(38\right)\end{array}$

i.e., the sufficient decrease condition in Line 29 holds for all t∈[0, T_δ^k). Thus, letting T_p^k be the appropriate T_p^k from either Case 1 or Case 2, and letting σ and σ_k,0 be taken from Algorithm 1, it is concluded that within

$\lceil {\log }_{\sigma }\left(\frac{\min \left\{T_p^k,T_{\delta }\right\}}{{\alpha }_{k,0}}\right)\rceil$

iterations, X^k+1 will satisfy both of the stopping criteria of the line search.

Lemma 3. Suppose that problem (14) satisfies Assumption 1. Let {X^k} denote the sequence of iterates generated by Algorithm 1 for (14). Then, ∥X^k+1−X^k∥→0. Moreover, the sequence {f(X^k)} admits a limit point f*.

Proof 6. If ∇f(X^k)=0 for any k, then the lemma is immediate. So, assume ∇f(X^k)≠0 for every k. Due to the stopping criteria for the line-search in Lines 16 and 28 of Algorithm 1, there is f(X^k+1)≤f_max^k for all k. Thus, f_max^k+1≤f_max^k for all k. The sequence {f_max^k} is monotone decreasing and bounded from below, hence there exists a limit point {f_max^k}→f*. Since f is uniformly continuous on the level set L⊂{X:σI_nXσI_n}, an induction argument can be used to conclude that ∥X^k+1−X^k∥→0.

Theorem 4. Suppose that problem (14) satisfies Assumption 1. Let {X^k} denote the sequence of iterates generated by Algorithm 1 for (14). If X* is an accumulation point of {X^k}, then X* satisfies the necessary optimality conditions stated in Lemma 1, i.e. there exists T>0 such that (27) holds.

Proof 7. Again, assume ∇f(X^k)≠0 for every k, since otherwise the theorem is immediate. Since X* is an accumulation point of {X^k}, there exists a subsequence of indices {k_i}_{i=0, 1, 2, . . .} such that

$\underset{i\to \infty }{\lim }X^{k_i}=X^{*}.$

By the definition or limit, for a sufficiently large k, the positions of the nonzero entries of X^ki are the same for k_i>k, and their magnitudes are uniformly bounded away from zero. In Lemma 2, it was proven that X^k∈L∈{σI_nXσI_n} for all k, the largest eigenvalue of ∇f(X^k)=S−(X^k)⁻¹+λX^k is bounded from above for all k. Thus, reusing the definition of T_p^k from the proof of Lemma 2 for either Case 1 or Case 2, there exists a constant T_p>0 such that T_p^ki≥T_p for all k_i>k. So, the step-size α_ki taken at the end of the k_i-th iteration is lower-bounded as

${\alpha }_{k_i}\ge {\alpha }_{k_i,0}{\sigma }^{\lceil {\log }_{\sigma }\left(\frac{\min \{{\stackrel{\_}{T}}_{p,T_{\delta }\}}}{{\alpha }_{k_i,0}}\right)\rceil }$

and so there exists α>0 such that

$\underset{i\to \infty }{\lim }{\alpha }^{k_i}\ge \stackrel{\_}{\alpha },$

since there are bounds α_min≤α_ki,0≤α_max as algorithmic parameters. Take any T such that 0<T<min(α,T_p). Consider two possible cases for the subsequence {X^ki+1}:

Case 1: Suppose that infinitely many iterates in {X^ki+1}_{i=0, 1, . . .} are generated in Step 2 (the gradient projection step). Then, without loss of generality, it can be assumed that the entire subsequence {X^ki+1}_{i=0, 1, . . .} is generated in Step 2. By triangle inequality, there is

$\underset{i\to \infty }{\lim }{X^{k_i+1}-X^{*}}_F\le \underset{i\to \infty }{\lim }{X^{k_i+1}-X^{k_i}}_F+\underset{i\to \infty }{\lim }{X^{k_i}-X^{*}}_F,$

and so it is concluded from Lemma 3 and the choice of subsequence {k_i} that

$\underset{i\to \infty }{\lim }X^{k_i+1}=X^{*}.$

For any t>0, let Y^ki∈S_t,In(X^ki). It follows from the definition of S_t,In that

$\begin{array}{cc}\mathrm{tr}\left(\nabla {f\left(X^{k_i}\right)}^T\left(X^{k_i+1}-X^{k_i}\right)\right)+\frac{1}{2t}{X^{k_i+1}-X^{k_i}}_F^2\ge \mathrm{tr}\left(\nabla {f\left(X^{k_i}\right)}^T\left(Y^{k_i}-X^{k_i}\right)\right)+\frac{1}{2t}{Y^{k_i}-X^{k_i}}_F^2.& \left(39\right)\end{array}$

Since

$X^{k_i+1}\in S_{{\alpha }_{k_i}.I_n}\left(X^{k_i}\right)$

for all k_i, there also is from the definition of

$S_{{\alpha }_{k_i}.I_n}$

that

$\begin{array}{cc}\mathrm{tr}\left(\nabla {f\left(X^{k_i}\right)}^T\left(Y^{k_i}-X^{k_i}\right)\right)+\frac{1}{2{\alpha }_{k_i}}{Y^{k_i}-X^{k_i}}_F^2\ge \mathrm{tr}\left(\nabla {f\left(X^{k_i}\right)}^T\left(X^{k_i+1}-X^{k_i}\right)\right)+\frac{1}{2{\alpha }_{k_i}}{X^{k_i+1}-X^{k_i}}_F^2& \left(40\right)\end{array}$

Combining (39) and (40), then, it is obtained for any t>0

$\begin{array}{cc}\mathrm{tr}\left(\nabla {f\left(X^{k_i}\right)}^T\left(X^{k_i+1}-X^{k_i}\right)\right)+\frac{1}{2t}{X^{k_i+1}-X^{k_i}}_F^2\ge \mathrm{tr}\left(\nabla {f\left(X^{k_i}\right)}^T\left(X^{k_i+1}-X^{k_i}\right)\right)+\frac{1}{2{\alpha }_{k_i}}{X^{k_i+1}-X^{k_i}}_F^2+\left(\frac{1}{2t}-\frac{1}{2{\alpha }_{k_i}}\right){Y^{k_i}-X^{k_i}}_F^2.& \left(41\right)\end{array}$

If t∈(0, T), then it may be concluded from (42) that

$\underset{i\to \infty }{\lim }{Y^{k_i}-X^{k_i}}_F=0,$

since

$\underset{i\to \infty }{\lim }{X^{k_i+1}-X^{k_i}}_F=0$

by Lemma 3 and since

$\frac{1}{2t}-\frac{1}{2{\alpha }_{k_i}}>0$

by the choice of t. By triangle inequality,

$\underset{i\to \infty }{\lim }{Y^{k_i}-X^{*}}_F\le \underset{i\to \infty }{\lim }{Y^{k_i}-X^{k_i}}_F+\underset{i\to \infty }{\lim }{X^{k_i}-X^{*}}_F,$

and hence

$\underset{i\to \infty }{\lim }Y^{k_i}\to X^{*}.$

So, using Proposition 1, it is concluded that X*∈P_Ω(X*−t∇f(X*)) for all t∈(0, T).

Case 2: Suppose that only finitely many elements of the subsequence {X^ki+1} are generated in Step 2. Then, without loss of generality, it can be assumed that every iterate of the subsequence {X^ki+1}_{i=0, 1 . . .} is generated in Step 1. By the pigeonhole principle, there exists an integer m∈[0, p) such that mod(k_i+1,q)=m for infinite many k_i. Hence one has mod(k_i+1+p−m,q)=p for infinite many k_i; that is, every iterate of the subsequence {X^ki+1+p−m} is generated in Step 2. From the triangle inequality,

$\underset{i\to \infty }{\lim }{X^{k_i+1+p-m}-X^{*}}_F\le \underset{i\to \infty }{\lim }\sum _{j=0}^{p-m}{X^{k_i+j+1}-X^{k_i+j}}_F+\underset{i\to \infty }{\lim }{X^{k_i}-X^{*}}_F,$

and so by Lemma 3 and the choice of the subsequence {X^ki}, X^ki+1+p−m→X*. From here, the same argument as in Case 1 can be applied to the subsequence {X^ki+1+p−m} to get X*∈P_Ω(X*−t∇f(X*)) for all t∈(0, T).It now only remains to prove the uniqueness of P_Ω(X*−t∇f(X*)), show that X*=P_Ω(X*−t∇f(X*)) as a set equality. Two separate cases are analyzed:Case A. (|I(X*)|=κ): Assume that Z*∈P_Ω(X*−t∇f(X)) for all t∈(0,T). Towards a contradiction, suppose that Z*≠X* for some t′∈(0, T). By the definition of P_Ω, it is concluded that because Z*≠X*, |X_i′,j′*|≤|Z_l′,m′*(t′)|, where there is defined

|X^i′,j′*|Δargmin{|X_u,v*|:(u,v)∈I(X*)}

and

|Z_l′*,m′*(t)|Δargmax{|(X*−t∇f(X*))_u,v|:(u,v)∈A(X*)}.

This implies that (∇f(X*))_l′,m′=0, since otherwise there exists {tilde over (t)}∈(0,T}) such that |(X*−{tilde over (t)}∇f(X*))_l′,m′|>|X_i′,j′*|, which would contradict X*∈P_Ω(X*−t∇f(X*)), since (i′, j′)∈I(X*) and (l′, m′)∈A(X*). But, if (∇f(X*))_l′,m′=0, then X*∉(X*−t∇f(X*)) for any t∈(0,T), since (X*−t∇f(X*))_l′,m′=(X*)_l′,m′, and (l′,m′)∈A(X*). This is a contradiction.

Case B. (|I(X*)|<κ): If (∇f(X*))_i,j≠0 for some (i, j)∈A(X*), then there exists {tilde over (t)}∈(0,T) such that |(X*−{tilde over (t)}∇f(X*))_i,|≠0; hence X*∉P_Ω(X*−{tilde over (t)}∇f(X*)), which would be a contradiction. Thus, it is the case that [∇f(X*)]_A(X*)*=0. Consider the reduced projection problem over the inactive indices in I(X*),

$\begin{array}{cc}\underset{Y_{i,j}}{\min }\sum _{\left(i,j\right)\in I\left(X^{*}\right)}{\left(Y_{i,j}-{\left(X^{*}-t\nabla f\left(X^{*}\right)\right)}_{i,j}\right)}^2.& \left(42\right)\end{array}$

Since H_i,j* is a solution to (42), it is derived that [∇f(X*)]_I(X*)=0; thus it is concluded that ∇f(X*)=0, and so it follows that X*=P_Ω(X*−t∇f(X*)) for any t∈(0,T).The following result further shows that any accumulation point produced by Algorithm 1 is a strict local minimizer of (14).

Theorem 5. Any accumulation point X* of {X^k}_{k=0, 1, . . .} generated by Algorithm 1 is a strict local minimizer of (14).

Proof 8. It may be proven that there exists ε>0 such that for all Y so that X*+Y is feasible with respect to the constraints of (14) and ∥Y∥_F∈(0,ε), it holds that f(X*+Y)>f(X*). From Theorem 4, there exists T>0 such that X*=P_Ω(X*−t∇f(X*)) for all t∈(0,T). Equivalently, for all such t, the matrix 0 is the unique global minimizer of

$\arg \underset{Y\text{:}X^{*}+Y\in \Omega }{\min }g_t\left(Y\right)\underset{\_}{\underset{\_}{\Delta }}f\left(X^{*}\right)+\mathrm{tr}\left(\nabla f\left(X^{*}\right)Y\right)+\frac{1}{2t}{Y}_F^2.$

Moreover, f is strongly convex on K={x:0.5σI_nºXº(σ+0.5σ)I_n} defined in (37). Denote ρ by the strong convexity constant. Note that X^k∈L⊂{σI_nºXºσI_n}, hence the limit X*∈{σI_nXσI_n}. In Lemma 2, it was proven that if ∥Y∥_F≤0.5σ then X*+Y∈K. For any ∥Y∥_F≤0.5σ, there is

$\begin{array}{cc}\begin{array}{cc}f\left(X^{*}+Y\right)\ge & \mathrm{ }f\left(X^{*}\right)+\mathrm{tr}\left(〈\nabla f\left(X^{*}\right),Y〉\right)+\frac{1}{2\rho }{Y}_F^2\\ =& f\left(X^{*}\right)+\frac{1}{2\rho }{Y+\rho \nabla f\left(X^{*}\right)}_F^2-\frac{\rho }{2}{\nabla f\left(X^{*}\right)}_F^2,\end{array}& \left(43\right)\end{array}$

where the equality is by completing the square. As in the proof of Theorem 4, the analysis may be split into two cases.

Case A (|I(X*)|=κ): Observe that there is Y_i,j=0 for every (i, j)∈A(X*), in order for X*+Y∈Ω to hold in an arbitrarily small neighborhood of X*. Also observe that (∇f(X*))_i,j=0 for every (i, j)∈I(X*); otherwise, X*∉P_Ω(X*−t∇f(X*)) for sufficiently small t, which would yield a contradiction. Taking

ε<min{0.5σ,min{|X_i,j*|:(i,j)∈I(X*)}},  (44)

there is X*+Y0 since 0.5σ<σ₁(X*). From (43), it is obtained

$\begin{array}{cc}f\left(X^{*}+Y\right)\ge & f\left(X^{*}\right)+\frac{1}{2\rho }\sum _{\left(i,j\right)\in I\left(X^{*}\right)}{\left(Y_{i,j}+{\rho \left(\nabla f\left(X^{*}\right)\right)}_{i,j}\right)}^2+\\ & \frac{1}{2\rho }\sum _{\left(i,j\right)\in A\left(X^{*}\right)}{\left(Y_{i,j}+{\rho \left(\nabla f\left(X^{*}\right)\right)}_{i,j}\right)}^2-\\ & \frac{\rho }{2}\sum _{\left(i,j\right)\in I\left(X^{*}\right)}{\left({\left(\nabla f\left(X^{*}\right)\right)}_{i,j}\right)}^2-\\ & \frac{\rho }{2}\sum _{\left(i,j\right)\in A\left(X^{*}\right)}{\left({\left(\nabla f\left(X^{*}\right)\right)}_{i,j}\right)}^2\\ =& f\left(X^{*}\right)+\frac{1}{2\rho }\sum _{\left(i,j\right)\in I\left(X^{*}\right)}Y_{i,j}^2+\\ & \frac{1}{2\rho }\sum _{\left(i,j\right)\in A\left(X^{*}\right)}{\left({\rho \left(\nabla f\left(X^{*}\right)\right)}_{i,j}\right)}^2-\\ & \frac{\rho }{2}\sum _{\left(i,j\right)\in A\left(X^{*}\right)}{\left({\left(\nabla f\left(X^{*}\right)\right)}_{i,j}\right)}^2\\ =f\left(X^{*}\right)& +\frac{1}{2\rho }\sum _{\left(i,j\right)\in I\left(X^{*}\right)}Y_{i,j}^2\\ >f\left(X^{*}\right),& \end{array}$

where the last strict inequality comes from the fact that not all the entries of Y are 0.

Case B. (|I(X*)|<κ): As in the proof of Case B in Theorem 4, there is ∇f(X*)=0. Hence from (43), there is

$f\left(X^{*}+Y\right)\ge f\left(X^{*}\right)+\frac{1}{2\rho }{Y}_F^2>f\left(X^{*}\right),$

for any ∥Y∥_F≤ε (ε defined in (44)), which completes the proof.

In one aspect, the number of nonzeros κ and the coefficient λ that penalizes the Frobenius norm may be tuned. In one aspect, it is noted that the performance of the model of Equation (3), for example, evaluated based on real data, is stable for a wide range of parameters K and A.

A system and/or method in embodiments provide for anomaly localization for multivariate time series datasets. An l₀-constrained model with an l₂-regularization in an objective function in one embodiment learns sparse dependency graphs, e.g., creates a depedency graph structure. In an embodiment, a projected gradient algorithm for the l₀-constrained problem is developed and proved that it converges to a stationary point. In embodiments, metrics are applied for scoring anomalies. An l₂-regularized l₀-constrained models of Equation (3) and Equation (4), for example, combined with a metric, e.g., the conditional expected Kullback-Liebler divergence anomaly scoring of Equation (8), improves methods for detecting anomalies.

FIG. 3 is a diagram showing components of a system in one embodiment, which performs localization for multivariate time series datasets detects anomalies in sensor data. One or more hardware processors 302 such as a central processing unit (CPU), a graphic process unit (GPU), and/or a Field Programmable Gate Array (FPGA), an application specific integrated circuit (ASIC), and/or another processor, may be coupled with a memory device 304, and constructs a first dependency graph based on a first data set by solving a first objective function constrained with a maximum number of non-zeros and formulated with a regularization term, which includes a quadratic penalty to control sparsity. In an embodiment, the quadratic penalty is determined as a function of the first data set in constructing the first dependency graph. One or more hardware processors 302 also constructs a second dependency graph based on the second data set by solving the first objective function constrained with the maximum number of non-zeros and formulated with the regularization term, which includes a quadratic penalty. In constructing the second dependency graph, the quadratic penalty is determined as a function of the first data set and the second data set.

In an embodiment, one or more hardware processors 302 receive via a communication or network interface 308, a first data set, which includes sensor data detected by a plurality of sensors coupled to equipments in operation during a first period of time, for example, shown at 312. During this first period of time, the equipments are considered to have functioned in normal manner (e.g., without particular errors or malfunctioning), and hence the first data set includes data associated with normal state of operations of the equipments. One or more hardware processors 302 also receive a second data set comprising sensor data detected by the plurality of sensors coupled to the equipments 312 in operation during a second period of time. The second period of time may be the current time, for example, a window of current time. One or more hardware processors 302 determines an anomaly score for each of the plurality of sensors based on comparing the first dependency graph and the second dependency graph.

One or more hardware processors 302, in an embodiment, responsive to determining that an anomaly score associated with a sensor in the plurality of sensors meets a threshold value indicative of abnormal functioning, automatically triggers a corrective action to correct an equipment coupled to the sensor. For instance, one or more hardware processors 302 may send a signal to adjust one or more parameters of one or more equipments 312, via a network interface 308.

In an embodiment, one or more hardware processors 302 solve the first object function implementing a projected gradient algorithm. The projected gradient algorithm embeds an approximate Newton method, wherein feasibility with respect to sparsity constraint is handled via projection, and a symmetric positive-definiteness of iterates is ensured through a line-search procedure.

The memory device 304 may include random access memory (RAM), read-only memory (ROM) or another memory device, and may store data and/or processor instructions for implementing various functionalities associated with the methods and/or systems described herein. A processor 302 may execute computer instructions stored in the memory 304 or received from another computer device or medium. The memory device 304 may, for example, store instructions and/or data for functioning of the one or more hardware processors 302, and may include an operating system and other program of instructions and/or data. In one aspect, input data such as the first data set and the second data set, and program instructions which one or more of the hardware processors 302 may execute, can be stored on a storage device 306 or received via a network interface 308 from a remote device, and may be temporarily loaded into the memory device 304 to construct dependency graphs and detect anomaly. One or more hardware processors 302 may be coupled with interface devices such as a network interface 308 for communicating with remote systems, for example, via a network, and an input/output interface 310 for communicating with input and/or output devices such as a keyboard, mouse, display, and/or others.

FIG. 4 is a diagram illustrating a screen display of an analytics tool in an embodiment, which integrate anomaly detection technique in one embodiment. Learning dependency models and detecting anomalies according to embodiments in the present disclosure can be integrated as an analytics function, for example, shown at 402, which can be selectable via a user interface 400 of an analytics tool. Selecting the function 402, for example, navigates a user through building of dependency graphs and detecting anomalies in connected sensor network, for example, of equipments in operation, for example, in a manufacturing or another industrial process.

FIG. 5 illustrates a schematic of an example computer or processing system that may implement a system in one embodiment of the present disclosure. The computer system is only one example of a suitable processing system and is not intended to suggest any limitation as to the scope of use or functionality of embodiments of the methodology described herein. The processing system shown may be operational with numerous other general purpose or special purpose computing system environments or configurations. Examples of well-known computing systems, environments, and/or configurations that may be suitable for use with the processing system shown in FIG. 5 may include, but are not limited to, personal computer systems, server computer systems, thin clients, thick clients, handheld or laptop devices, multiprocessor systems, microprocessor-based systems, set top boxes, programmable consumer electronics, network PCs, minicomputer systems, mainframe computer systems, and distributed cloud computing environments that include any of the above systems or devices, and the like.

The computer system may be described in the general context of computer system executable instructions, such as program modules, being executed by a computer system. Generally, program modules may include routines, programs, objects, components, logic, data structures, and so on that perform particular tasks or implement particular abstract data types. The computer system may be practiced in distributed cloud computing environments where tasks are performed by remote processing devices that are linked through a communications network. In a distributed cloud computing environment, program modules may be located in both local and remote computer system storage media including memory storage devices.

The components of computer system may include, but are not limited to, one or more processors or processing units 12, a system memory 16, and a bus 14 that couples various system components including system memory 16 to processor 12. The processor 12 may include a module 30 that performs the methods described herein. The module 30 may be programmed into the integrated circuits of the processor 12, or loaded from memory 16, storage device 18, or network 24 or combinations thereof.

Bus 14 may represent one or more of any of several types of bus structures, including a memory bus or memory controller, a peripheral bus, an accelerated graphics port, and a processor or local bus using any of a variety of bus architectures. By way of example, and not limitation, such architectures include Industry Standard Architecture (ISA) bus, Micro Channel Architecture (MCA) bus, Enhanced ISA (EISA) bus, Video Electronics Standards Association (VESA) local bus, and Peripheral Component Interconnects (PCI) bus.

Computer system may include a variety of computer system readable media. Such media may be any available media that is accessible by computer system, and it may include both volatile and non-volatile media, removable and non-removable media.

System memory 16 can include computer system readable media in the form of volatile memory, such as random access memory (RAM) and/or cache memory or others. Computer system may further include other removable/non-removable, volatile/non-volatile computer system storage media. By way of example only, storage system 18 can be provided for reading from and writing to a non-removable, non-volatile magnetic media (e.g., a “hard drive”). Although not shown, a magnetic disk drive for reading from and writing to a removable, non-volatile magnetic disk (e.g., a “floppy disk”), and an optical disk drive for reading from or writing to a removable, non-volatile optical disk such as a CD-ROM, DVD-ROM or other optical media can be provided. In such instances, each can be connected to bus 14 by one or more data media interfaces.

Computer system may also communicate with one or more external devices 26 such as a keyboard, a pointing device, a display 28, etc.; one or more devices that enable a user to interact with computer system; and/or any devices (e.g., network card, modem, etc.) that enable computer system to communicate with one or more other computing devices. Such communication can occur via Input/Output (I/O) interfaces 20.

Still yet, computer system can communicate with one or more networks 24 such as a local area network (LAN), a general wide area network (WAN), and/or a public network (e.g., the Internet) via network adapter 22. As depicted, network adapter 22 communicates with the other components of computer system via bus 14. It should be understood that although not shown, other hardware and/or software components could be used in conjunction with computer system. Examples include, but are not limited to: microcode, device drivers, redundant processing units, external disk drive arrays, RAID systems, tape drives, and data archival storage systems, etc.

The present invention may be a system, a method, and/or a computer program product. The computer program product may include a computer readable storage medium (or media) having computer readable program instructions thereon for causing a processor to carry out aspects of the present invention.

The computer readable storage medium can be a tangible device that can retain and store instructions for use by an instruction execution device. The computer readable storage medium may be, for example, but is not limited to, an electronic storage device, a magnetic storage device, an optical storage device, an electromagnetic storage device, a semiconductor storage device, or any suitable combination of the foregoing. A non-exhaustive list of more specific examples of the computer readable storage medium includes the following: a portable computer diskette, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or Flash memory), a static random access memory (SRAM), a portable compact disc read-only memory (CD-ROM), a digital versatile disk (DVD), a memory stick, a floppy disk, a mechanically encoded device such as punch-cards or raised structures in a groove having instructions recorded thereon, and any suitable combination of the foregoing. A computer readable storage medium, as used herein, is not to be construed as being transitory signals per se, such as radio waves or other freely propagating electromagnetic waves, electromagnetic waves propagating through a waveguide or other transmission media (e.g., light pulses passing through a fiber-optic cable), or electrical signals transmitted through a wire.

Computer readable program instructions described herein can be downloaded to respective computing/processing devices from a computer readable storage medium or to an external computer or external storage device via a network, for example, the Internet, a local area network, a wide area network and/or a wireless network. The network may comprise copper transmission cables, optical transmission fibers, wireless transmission, routers, firewalls, switches, gateway computers and/or edge servers. A network adapter card or network interface in each computing/processing device receives computer readable program instructions from the network and forwards the computer readable program instructions for storage in a computer readable storage medium within the respective computing/processing device.

Computer readable program instructions for carrying out operations of the present invention may be assembler instructions, instruction-set-architecture (ISA) instructions, machine instructions, machine dependent instructions, microcode, firmware instructions, state-setting data, or either source code or object code written in any combination of one or more programming languages, including an object oriented programming language such as Smalltalk, C++ or the like, and conventional procedural programming languages, such as the “C” programming language or similar programming languages. The computer readable program instructions may execute entirely on the user's computer, partly on the user's computer, as a stand-alone software package, partly on the user's computer and partly on a remote computer or entirely on the remote computer or server. In the latter scenario, the remote computer may be connected to the user's computer through any type of network, including a local area network (LAN) or a wide area network (WAN), or the connection may be made to an external computer (for example, through the Internet using an Internet Service Provider). In some embodiments, electronic circuitry including, for example, programmable logic circuitry, field-programmable gate arrays (FPGA), or programmable logic arrays (PLA) may execute the computer readable program instructions by utilizing state information of the computer readable program instructions to personalize the electronic circuitry, in order to perform aspects of the present invention.

Aspects of the present invention are described herein with reference to flowchart illustrations and/or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and/or block diagrams, and combinations of blocks in the flowchart illustrations and/or block diagrams, can be implemented by computer readable program instructions.

These computer readable program instructions may be provided to a processor of a general purpose computer, special purpose computer, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create means for implementing the functions/acts specified in the flowchart and/or block diagram block or blocks. These computer readable program instructions may also be stored in a computer readable storage medium that can direct a computer, a programmable data processing apparatus, and/or other devices to function in a particular manner, such that the computer readable storage medium having instructions stored therein comprises an article of manufacture including instructions which implement aspects of the function/act specified in the flowchart and/or block diagram block or blocks.

The computer readable program instructions may also be loaded onto a computer, other programmable data processing apparatus, or other device to cause a series of operational steps to be performed on the computer, other programmable apparatus or other device to produce a computer implemented process, such that the instructions which execute on the computer, other programmable apparatus, or other device implement the functions/acts specified in the flowchart and/or block diagram block or blocks.

The flowchart and block diagrams in the Figures illustrate the architecture, functionality, and operation of possible implementations of systems, methods, and computer program products according to various embodiments of the present invention. In this regard, each block in the flowchart or block diagrams may represent a module, segment, or portion of instructions, which comprises one or more executable instructions for implementing the specified logical function(s). In some alternative implementations, the functions noted in the block may occur out of the order noted in the figures. For example, two blocks shown in succession may, in fact, be executed substantially concurrently, or the blocks may sometimes be executed in the reverse order, depending upon the functionality involved. It will also be noted that each block of the block diagrams and/or flowchart illustration, and combinations of blocks in the block diagrams and/or flowchart illustration, can be implemented by special purpose hardware-based systems that perform the specified functions or acts or carry out combinations of special purpose hardware and computer instructions.

The terminology used herein is for the purpose of describing particular embodiments only and is not intended to be limiting of the invention. As used herein, the singular forms “a”, “an” and “the” are intended to include the plural forms as well, unless the context clearly indicates otherwise. It will be further understood that the terms “comprises” and/or “comprising,” when used in this specification, specify the presence of stated features, integers, steps, operations, elements, and/or components, but do not preclude the presence or addition of one or more other features, integers, steps, operations, elements, components, and/or groups thereof.

The corresponding structures, materials, acts, and equivalents of all means or step plus function elements, if any, in the claims below are intended to include any structure, material, or act for performing the function in combination with other claimed elements as specifically claimed. The description of the present invention has been presented for purposes of illustration and description, but is not intended to be exhaustive or limited to the invention in the form disclosed. Many modifications and variations will be apparent to those of ordinary skill in the art without departing from the scope and spirit of the invention. The embodiment was chosen and described in order to best explain the principles of the invention and the practical application, and to enable others of ordinary skill in the art to understand the invention for various embodiments with various modifications as are suited to the particular use contemplated.

Claims

1. A method, comprising: receiving a first data set comprising sensor data detected by a plurality of sensors coupled to equipments in operation during a first period of time;constructing a first dependency graph based on the first data set by solving a first objective function constrained with a maximum number of non-zeros and formulated with a regularization term comprising a quadratic penalty to control sparsity, wherein the quadratic penalty is determined as a function of the first data set in constructing the first dependency graph;receiving a second data set comprising sensor data detected by the plurality of sensors coupled to the equipments in operation during a second period of time;constructing a second dependency graph based on the second data set by solving the first objective function constrained with the maximum number of non-zeros and formulated with the regularization term comprising a quadratic penalty, wherein the quadratic penalty is determined as a function of the first data set and the second data set in constructing the second dependency graph; anddetermining an anomaly score for each of the plurality of sensors based on comparing the first dependency graph and the second dependency graph.

2. The method of claim 1, further comprising: responsive to determining that an anomaly score associated with a sensor in the plurality of sensors meets a threshold value indicative of abnormal functioning, automatically performing a corrective action to correct equipment coupled to the sensor.

3. The method of claim 1, wherein the first object function is solved by a projected gradient algorithm.

4. The method of claim 3, wherein the projected gradient algorithm embeds an approximate Newton method, wherein feasibility with respect to sparsity constraint is handled via projection, and a symmetric positive-definiteness of iterates is ensured through a line-search procedure.

5. The method of claim 1, wherein the determining of the anomaly score for each of the plurality of sensors comprises implementing a conditional expected Kullback-Liebler divergence between the first dependency graph and the second dependency graph.

6. The method of claim 1, wherein the determining of the anomaly score for each of the plurality of sensors comprises implementing a stochastic nearest neighbors algorithm measuring dissimilarity between neighborhood of i-th variable representing a sensor between the first dependency graph and the second dependency graph.

7. The method of claim 1, wherein the determining of the anomaly score for each of the plurality of sensors comprises implementing a sparsest subgraph approximation based on the first dependency graph and the second dependency graph.

8. The method of claim 1, wherein the equipments comprise equipments of an industrial process.

9. The method of claim 1, wherein the equipments comprise equipments of a manufacturing process manufacturing a product.

10. A computer readable storage device storing a program of instructions executable by a machine to perform a method comprising: receiving a first data set comprising sensor data detected by a plurality of sensors coupled to equipments in operation during a first period of time;constructing a first dependency graph based on the first data set by solving a first objective function constrained with a maximum number of non-zeros and formulated with a regularization term comprising a quadratic penalty to control sparsity, wherein the quadratic penalty is determined as a function of the first data set in constructing the first dependency graph;receiving a second data set comprising sensor data detected by the plurality of sensors coupled to the equipments in operation during a second period of time;constructing a second dependency graph based on the second data set by solving the first objective function constrained with the maximum number of non-zeros and formulated with the regularization term comprising a quadratic penalty, wherein the quadratic penalty is determined as a function of the first data set and the second data set in constructing the second dependency graph; anddetermining an anomaly score for each of the plurality of sensors based on comparing the first dependency graph and the second dependency graph.

11. The computer readable storage device of claim 10, further comprising: responsive to determining that an anomaly score associated with a sensor in the plurality of sensors meets a threshold value indicative of abnormal functioning, automatically performing a corrective action to correct equipment coupled to the sensor.

12. The computer readable storage device of claim 10, wherein the first object function is solved by a projected gradient algorithm.

13. The computer readable storage device of claim 12, wherein the projected gradient algorithm embeds an approximate Newton method, wherein feasibility with respect to sparsity constraint is handled via projection, and a symmetric positive-definiteness of iterates is ensured through a line-search procedure.

14. The computer readable storage device of claim 10, wherein the determining of the anomaly score for each of the plurality of sensors comprises implementing a conditional expected Kullback-Liebler divergence between the first dependency graph and the second dependency graph.

15. The computer readable storage device of claim 10, wherein the determining of the anomaly score for each of the plurality of sensors comprises implementing a stochastic nearest neighbors algorithm measuring dissimilarity between neighborhood of i-th variable representing a sensor between the first dependency graph and the second dependency graph.

16. The computer readable storage device of claim 10, wherein the determining of the anomaly score for each of the plurality of sensors comprises implementing a sparsest subgraph approximation based on the first dependency graph and the second dependency graph.

17. A system, comprising: a hardware processor coupled with a communication interface;a memory device coupled to the hardware processor;the hardware processor operable to receive via the communication interface, a first data set comprising sensor data detected by a plurality of sensors coupled to equipments in operation during a first period of time;the hardware processor operable to construct a first dependency graph based on the first data set by solving a first objective function constrained with a maximum number of non-zeros and formulated with a regularization term comprising a quadratic penalty to control sparsity, wherein the quadratic penalty is determined as a function of the first data set in constructing the first dependency graph;the hardware processor operable to receive a second data set comprising sensor data detected by the plurality of sensors coupled to the equipments in operation during a second period of time;the hardware processor operable to construct a second dependency graph based on the second data set by solving the first objective function constrained with the maximum number of non-zeros and formulated with the regularization term comprising a quadratic penalty, wherein the quadratic penalty is determined as a function of the first data set and the second data set in constructing the second dependency graph; andthe hardware processor operable to determine an anomaly score for each of the plurality of sensors based on comparing the first dependency graph and the second dependency graph.

18. The system of claim 17, wherein the hardware processor is further operable to, responsive to determining that an anomaly score associated with a sensor in the plurality of sensors meets a threshold value indicative of abnormal functioning, automatically trigger a corrective action to correct an equipment coupled to the sensor.

19. The system of claim 17, wherein the first object function is solved by a projected gradient algorithm.

20. The system of claim 19, wherein the projected gradient algorithm embeds an approximate Newton method, wherein feasibility with respect to sparsity constraint is handled via projection, and a symmetric positive-definiteness of iterates is ensured through a line-search procedure.