PaperlayerView exact revision on arXiv

Source-linked AI summary

Signal Processing on Graphs: Causal Modeling of Unstructured Data

Jonathan Mei, José M. F. Moura

arXiv:1503.00173v6cs.ITstat.ML

TL;DR

The paper addresses how to recover graph structure from large collections of temporally dependent, otherwise unstructured time series. It introduces a tractable causal graph-process framework and estimation algorithms for a directed, weighted adjacency matrix, with analyses and evaluations showing close recovery in simulations and consistency with physical knowledge in real data.

  • Problem

    Unknown relationships among large collections of time series must be represented despite data diversity, while static dimensionality-reduction methods omit temporal structure.

  • Method

    The paper models temporal dependencies with an autoregressive process whose coefficients are graph-polynomial filters of a possibly directed, weighted adjacency matrix, then estimates the graph and coefficients in separate tractable steps.

  • Results

    The estimated adjacency matrices are close to the true graphs in simulated data and consistent with prior physical knowledge in the tested real dataset.

  • Takeaways & Limitations

    The methodology provides a computationally tractable way to structure multiple time series through spatial and temporal dependencies using a graph that may capture causal relations.

Abstract

from arXiv · show

Many applications collect a large number of time series, for example, the financial data of companies quoted in a stock exchange, the health care data of all patients that visit the emergency room of a hospital, or the temperature sequences continuously measured by weather stations across the US. These data are often referred to as unstructured. A first task in its analytics is to derive a low dimensional representation, a graph or discrete manifold, that describes well the interrelations among the time series and their intrarelations across time. This paper presents a computationally tractable algorithm for estimating this graph that structures the data. The resulting graph is directed and weighted, possibly capturing causal relations, not just reciprocal correlations as in many existing approaches in the literature. A convergence analysis is carried out. The algorithm is demonstrated on random graph datasets and real network time series datasets, and its performance is compared to that of related methods. The adjacency matrices estimated with the new method are close to the true graph in the simulated data and consistent with prior physical knowledge in the real dataset tested.

I. INTRODUCTION

The paper addresses graph-structure estimation for large collections of time series whose relationships are unknown and temporally dependent. It models these dependencies with a possibly directed, weighted adjacency matrix associated with causal network effects.

  • Motivation: Large, simultaneous time-series collections arise across disciplines and are commonly treated as unstructured data.Examples include finance, health care, utilities, biology, and social interactions.
  • Motivation: Graph representations provide low-dimensional structure for statistics, inference, signal processing, and analysis of network diffusion.These methods extract knowledge from observed data supported on graphs.
  • Problem: Unknown graph structure must be inferred because data diversity and unstructuredness challenge first-principles modeling.Earlier dimensionality-reduction methods use static snapshots and omit time.
  • Problem: The target is a possibly directed, weighted adjacency matrix A capturing dependencies among time series.The paper associates the graph with causal network effects rather than only Markov relations or Granger-style causality.
  • Relation to Prior Work: Unlike inverse-covariance and SVAR approaches, the proposed model uses one weighted adjacency matrix with graph-filter coefficients.This can be more parsimonious than estimating multiple weighted graphs and supports analysis through the DSPG framework.

C. Graph Signal Processing using the Laplacian

The paper uses the adjacency matrix as the fundamental graph-shift operator and models time-series evolution through causal graph processes. Its formulation incorporates directed, weighted structure, temporal propagation constraints, identifiability conditions, and sparse estimation.

  • C. Graph Signal Processing using the Laplacian: Laplacian-based methods often assume symmetric matrices, which restricts them to undirected relations and enables convex spectral machinery.Some asymmetric Laplacian approaches additionally impose zero row or column sums.
  • C. Graph Signal Processing using the Laplacian: The proposed adjacency matrix A may contain positive and negative weights and complex eigenvalues, while still permitting Laplacian computation.The adjacency matrix is adopted as the basic building block.
  • III. CAUSAL GRAPH PROCESSES: A causal graph process is autoregressive in time, with coefficients represented as matrix polynomials of the adjacency matrix.The model is not Markovian in graph nodes and edges.
  • III. CAUSAL GRAPH PROCESSES: The polynomial order is limited by both the time lag i and the degree NA of A's minimum polynomial.The paper typically assumes model order M ≪ NA, so effects from i steps ago have at most i-th order network reach.
  • III. CAUSAL GRAPH PROCESSES: Identifiability is enforced by setting c10 = 0 and c11 = 1 so that P1(A, c) = A.The reduced parameterization uniquely specifies A and c under the stated nondegeneracy assumption.
  • Estimation: Estimation is decomposed into recovering filter matrices Ri, recovering A, and estimating polynomial coefficients cij.The decomposition is intended to make the estimation more tractable than the joint formulation.

A. Solving for Pi(A, c)

The first estimation stage solves for graph-filter matrices using a multi-convex formulation and block coordinate descent. Sparsity is exploited to reduce computational cost, after which A and the polynomial coefficients are recovered.

  • A. Solving for Pi(A, c): The optimization over filter matrices Ri is non-convex but multi-convex, becoming convex in each Ri when the others are fixed.This structure motivates block coordinate descent.
  • A. Solving for Pi(A, c): Each Ri subproblem is an ℓ1-regularized least-squares problem solvable with standard gradient-based methods.The optimization separates estimation across M blocks.
  • A. Solving for Pi(A, c): The naive worst-case complexity is O(M^2N^3 + KMN^2), dominated by matrix operations across the M blocks.Sparse matrix multiplication reduces the matrix-multiplication contribution when the Ri are sparse.
  • A. Solving for Pi(A, c): After estimating Ri, the adjacency matrix can be obtained from R1 or by jointly using the estimated filter matrices.The second approach has worst-case complexity O(MN^3), reduced for sparse matrices.
  • A. Solving for Pi(A, c): The coefficients cij can be estimated from the recovered A and Ri or from A together with the observed data.The resulting optimization can use standard ℓ1-regularized least-squares methods.

D. Base Estimation Algorithm

The paper develops a three-step procedure to estimate the adjacency matrix and filter coefficients, offering a more efficient alternative to direct nonconvex optimization. A simplified variant removes the commutativity term for easier computation and analysis, while an extended variant uses the estimates to initialize the direct approach.

  • Generalization: The modeling framework can accommodate non-Gaussian noise and certain nonlinear dependencies through a likelihood-based loss and regularization functions.The resulting optimization remains nonconvex because of the matrix-polynomial terms.
  • Base algorithm: The base algorithm alternates estimation of matrix coefficients, adjacency structure, and filter coefficients in three steps.It estimates the matrix coefficients, sets or derives the adjacency estimate, and then solves for the filter coefficients.
  • Base algorithm: The three-step procedure is described as more efficient and well behaved than directly optimizing the joint nonconvex formulation.
  • Simplified algorithm: The simplified algorithm removes the commutativity term, making the joint optimization over all matrix coefficients convex.Its worst-case cost is O(M(K−M)N^2), reducible to O(M(K−M)SMN) with sparse matrix-vector products in the best case.
  • Extended algorithm: The extended algorithm initializes the direct joint optimization with estimates from the base algorithm to seek better local minima.

V. CONVERGENCE OF ESTIMATION

The convergence analysis establishes conditions under which the proposed algorithms converge or achieve theoretical guarantees. These conditions include standard time-series assumptions, boundedness and sparsity, sufficient samples, and stronger-than-stability restrictions.

  • Convergence results: Convergence of the base and extended algorithms follows from existing results, while the simplified algorithm receives a dedicated performance proof.
  • Base algorithm: Block coordinate descent is not generally guaranteed to reach a global or local minimum, but converges for the matrix-coefficient updates under continuity and compact convex sublevel-set assumptions.
  • Extended algorithm: The extended algorithm uses a bounded-below objective with compact sublevel sets and appropriately chosen subgradient step sizes to produce nonincreasing objective values.
  • Assumptions: The theoretical analysis assumes a stationary CGP with uncorrelated i.i.d. Gaussian noise, sparse bounded parameters, and sufficient samples relative to process stability, network size, and sparsity.
  • Limitations of assumptions: The analysis identifies explicit sparsity as its most restrictive assumption and notes that the stronger-than-stability conditions are not implied by stationarity.The sparsity restriction corresponds roughly to graphs whose higher-order adjacency powers remain sufficiently sparse.

2) Theoretical Performance:

The paper proves high-probability estimation and prediction guarantees for the simplified algorithm under assumptions (A1)–(A8). The guarantees improve with sample size while allowing network size and autoregressive order to grow slowly relative to the observations.

  • Main result: Theorem 1 shows that assumptions (A1)–(A8) suffice for small average excess prediction risk with high probability.The result applies as the number of nodes N and time observations K grow, provided growth remains sufficiently controlled.
  • Scaling behavior: The guarantees permit the autoregressive order M and network size N to grow, as long as they grow slowly enough relative to the sample size K.For large K, the theorem states that the resulting error δ is much smaller than the bound L.
  • Proof strategy: The proof first bounds adjacency estimation error and filter-coefficient estimation error, then translates those bounds into a prediction-error guarantee.
  • Adjacency estimation: With an appropriate regularization parameter λ1, the adjacency estimate is consistent with high probability as the sample size increases.Lemma 2 provides an exponentially decreasing failure probability εA = d1 exp{−d2K/ω^2}.
  • Filter estimation: Under the same assumptions and an appropriate λ2, the filter-coefficient estimate is also good with high probability.

VI. EXPERIMENTS

The experiments evaluate CGP graph estimation on temperature-network and synthetic time-series data, comparing directed learned graphs with MRF and undirected distance-graph models. On temperature data, CGP achieves lower test error than MRF at comparable sparse parameterizations and captures geographically meaningful weather patterns.

  • A. Temperature Data: The evaluation uses a 150-city, 365-day temperature dataset and compares CGP Algorithms 1–3 with sparse vector autoregressive MRF and distance-graph models.The temperature series are detrended before training and split into even and odd time indices for training and testing.
  • A. Temperature Data: Directed graphs estimated from data outperform undirected distance graphs in testing, including at high sparsity where pnnz < 0.3.The comparison uses prediction error on held-out time-series data.
  • A. Temperature Data: Prediction performance is also evaluated against the total number of nonzero model parameters, with the same trends reported for the sparsity comparison.For CGP, the parameter count includes nonzero adjacency entries, diagonal entries, and polynomial coefficients.
  • A. Temperature Data: CGP has lower test error than MRF with fewer nonzero parameters when Nnnz < 2000.At the same sparsity level, MRF uses approximately M times more nonzero parameters than CGP.
  • A. Temperature Data: The estimated CGP temperature graph identifies west-to-east wind influence and the north-northwest-to-south-southeast influence of the Rocky Mountains.The graph is estimated with order M = 2 and sparsity pnnz = 0.05, using longitude and latitude for visualization.

B. Synthetic Data

Synthetic experiments test graph estimation across K-regular, stochastic block-model, Erdős–Rényi, and power-law topologies at varying graph and time-series sizes. Estimation errors generally improve with more samples for several topologies, while behavior varies with topology-specific agreement with the theoretical assumptions.

  • B. Synthetic Data: Synthetic data use K-regular, stochastic block-model, Erdős–Rényi, and power-law graphs with N ∈ {1000, 1500}, fixed M = 3, and varying sample sizes.The generated processes include sparse polynomial coefficients and Gaussian noise, with 20 independent Monte Carlo datasets used for estimation.
  • B. Synthetic Data: Figure 3 compares true and estimated absolute edge weights for the four topologies, using display parameters that favor fewer false alarms over missed edges.Higher absolute weights appear darker blue; layouts are selected for visualization because the graphs are large.
  • B. Synthetic Data: For representative parameter settings, the miss-rate and false-alarm tradeoff varies across topologies, with PM values from 0.07 for KR to 0.52 for PL.The listed pairs are KR: PFA = 0.01, PM = 0.07; SBM: 0.01, 0.40; ER: 0.03, 0.25; PL: 0.03, 0.52.
  • B. Synthetic Data: Average adjacency and prediction errors decrease with N and K for KR and SBM, while ER shows no clear N trend and PL prediction error fluctuates around zero.More samples improve ER adjacency estimation and prediction performance even without a clear N trend.
  • B. Synthetic Data: The KR results empirically conform to the theory, whereas other topologies may require loosened assumptions or additional network statistics such as diameter or maximum degree.The difference is attributed to topology-dependent network properties and their relationship to the theoretical sparsity assumptions.
  • B. Synthetic Data: The experiments illustrate both advantages and limitations of the proposed approach across real and simulated network time series.The paper’s overall experimental program combines temperature data with four synthetic graph families.

APPENDIX A PROOF OF LEMMA 2

The appendix proves a recovery result for the graph adjacency matrix under stability, sample-size, sparsity, and regularity assumptions. The proof bounds system quantities and invokes restricted-eigenvalue and deviation conditions to establish the lemma.

  • APPENDIX A PROOF OF LEMMA 2: With sufficiently large sample size, the lemma gives a high-probability ℓ1 error bound for the estimated adjacency matrix bA.The bound uses εA = d1 exp{−d2K/ω^2} and holds with probability at least 1 − εA under the stated assumptions.
  • APPENDIX A PROOF OF LEMMA 2: The proof begins by defining quantities that represent system stability and are used throughout the argument.The stacked state is represented in companion form through a system matrix.
  • APPENDIX A PROOF OF LEMMA 2: The stacked state obeys ex[k] = eAex[k−1] + ew[k], where the state contains current and lagged time-series vectors and the noise is stacked correspondingly.This companion-form representation converts the autoregressive process into a first-order state-space relation.
  • APPENDIX A PROOF OF LEMMA 2: Under assumption (A5), the proof bounds µmin(A) and µmax(A) through B(z) = I − A(z), producing the stability quantity Q.The bounds are then inserted into results from prior propositions to obtain Lemma 2.
  • APPENDIX A PROOF OF LEMMA 2: The proof of the extended result relies on restricted-eigenvalue and deviation conditions that characterize the geometry of the sparse optimization objective.Together, these conditions imply the desired recovery result under assumptions (A1)–(A8).

A. Eigenvalue and Deviation Conditions

The analysis establishes restricted eigenvalue and deviation conditions first for matrices built from the true parameter, then transfers them to matrices built from the estimated parameter. These conditions yield a high-probability ℓ1 error bound for estimating c.

  • Condition meanings: The Re condition ensures sufficient objective curvature in sparse directions, linking large objective perturbations to small parameter-estimation perturbations.It is defined through the matrix Γ and vectors θ, and is used in sparse estimation.
  • Condition meanings: The De condition requires small gradients near the true c, so the objective minimizer remains near the true parameter.It is expressed as ∥γ −Γc∥∞≤q.
  • Guarantee: ∥bc −c∥1 ≤64sMλ2/α1 for λ2 ≥q1 when the transferred Re and De conditions hold with high probability.This is the resulting sparse-estimation error guarantee for c.
  • Proof strategy: The estimated matrices are difficult to analyze because B(bA) and Y(bA) are complicated random functions of the true parameters and data rather than multivariate Gaussian variables.For the true-parameter matrices B(A) and Y(A), the corresponding conditions are comparatively straightforward to show.
  • Proof strategy: The proof first establishes Re(α0,τ) and De(q0) for B(A) and Y(A), then transfers them to B(bA) and Y(bA) available in the second stage.This transfer addresses the fact that the second-stage estimator cannot use the true A directly.

D. Proof of Supplementary Proposition 2

Supplementary Proposition 2 establishes the deviation condition for estimating c by bounding the relevant noise-gradient term. The proof uses concentration and union-bound arguments under sparsity and time-series assumptions.

  • Proposition statement: P(∥B(A)⊤e∥∞≥q0) ≤εDe, with q0 = 2LT 1−β√NtNπσug(Q) and εDe = 6M exp{−d4T 1−2β}.Here T = K −M and 0 < β < 1/2.
  • Proof: The proof defines e = Y(A)−B(A)c and applies Proposition 2.4 together with a union bound to control the maximum coordinate of B(A)⊤e.The argument uses the maximum sparsity tN of the evolution matrices.
  • Proof: Choosing η = T −β yields the stated deviation bound because η2 < η for 0 < β < 1/2.The resulting failure probability is bounded by 6M exp{−d4T 1−2β}.
  • Proof: The subsequent argument combines the deviation probability with an additional exponential term bounded using T < K.The combined bound is εDe +(6M +1) exp{−d4T}.

APPENDIX C PROOF OF THEOREM 1

Theorem 1 combines the preceding estimation bounds under assumptions (A1)–(A8). Its proof controls the estimated parameters and then bounds the resulting terms in the theorem’s error expression.

  • Theorem statement: For any 0 < β < ν < 1/2, assumptions (A1)–(A8) are sufficient for the error ϵ in (16) to satisfy the theorem’s stated bound.The theorem also invokes a universal constant d1.
  • Proof: The proof applies a union bound to the results from Lemmas 2 and 3 before proceeding under ∥bA −A∥2 ≤δA and ∥bc −c∥1 ≤δc.These inequalities provide the working parameter-error bounds used in the remainder of the proof.
  • Proof: The proof bounds ∥bA∥2 by L+δA and ∥bc∥1 by ρ + δc, then controls the remaining terms in the theorem’s error expression.The later bounds use the parameter-error inequalities established earlier.
Loading 1503.00173v6…