Covariance Matrix

Calculation and analysis of the network covariance matrix.

class covseisnet.covariance.CovarianceMatrix(input_array: ndarray | CovarianceMatrix | list, stats: list[Stats] | list[dict] | None = None, stft: ShortTimeFourierTransform | None = None)

Bases: ndarray

This class inherits all the methods of standard numpy arrays, along with extra methods tailored for array processing. Any numpy method or function applied to an instance of CovarianceMatrix will return an instance of CovarianceMatrix.

Mathematical definition

We consider the set of network seismograms \(\{v_i(t)\}_{i=1}^N\) with \(N\) traces, all recorded on the same time samples (synchronized). For the sake of simplicity, we use the continuous time notation. We first calculate the short-time Fourier transform of each trace within a sliding window of duration \(T\), and center time \(t_m\) (with the index \(m\) the time window index), defined on the time interval \(\tau_m = [t_m - T/2, t_m + T/2]\). The short-time Fourier transform of the traces is defined as

\[u_{i,m}(f) = \mathcal{F} v_i(t \in \tau_m)\]

The spectral covariance matrix is defined as the average of the outer product of the short-time Fourier transform of the traces over a set of time windows \(\{\tau_m\}_{m=1}^M\), with \(M\) the number of windows used to estimate the covariance. This lead to the definition of the spectral covariance matrix \(C_{ij}(f)\) as

\[C_{ij}(f) = \sum_{m=1}^M u_{i,m}(f) u_{j,m}^*(f)\]

where \(M\) is the number of windows used to estimate the covariance, \(^*\) is the complex conjugate, and \(\tau_m\) is the time window of index \(m\). The covariance matrix is a complex-valued matrix of shape \(N \times N\).

Practical implementation

The short-time Fourier transform is calculated with the help of the ShortTimeFourierTransform class. Both the calculation of the short-time Fourier transform and the covariance matrix can run in parallel depending on the available resources. Depending on the averaging size and frequency content, the covariance matrix can have a different shape:

• Shape is (n_traces, n_traces) if a single frequency and time window is given. This configuration is also obtained when slicing or reducing a covariance matrix object.

• Shape is (n_frequencies, n_traces, n_traces) if only one time frame is given, for n_frequencies frequencies, which depends on the window size and sampling rate.

• Shape is (n_times, n_frequencies, n_traces, n_traces) if multiple time frames are given.

All the methods defined in the the CovarianceMatrix class are performed on the flattened array with the flat() method, which allow to obtain as many \(N \times N\) covariance matrices as time and frequency samples.

Tip

The CovarianceMatrix class is rarely meant to be instantiated directly. It should most often be obtained from the output of the calculate_covariance_matrix() function as in the following example:

>>> import covseisnet as csn
>>> import numpy as np
>>> stream = csn.NetworkStream()
>>> time, frequency, covariance = csn.calculate_covariance_matrix(
...     stream, window_duration=10.0, average=10
... )
>>> print(type(covariance))
<class 'covseisnet.covariance.CovarianceMatrix'>

Note

If for some reason you need to create a CovarianceMatrix object from a numpy array, you have the several solutions provided by the numpy subclassing mechanism. We recommend using the contructor of the class, which allows to pass the stats and stft attributes, as in the following example:

>>> import covseisnet as csn
>>> import numpy as np
>>> covariance = csn.CovarianceMatrix(
...     np.zeros((4, 4)), stats=None, stft=None
... )
>>> covariance
CovarianceMatrix([[ 0.,  0.,  0.,  0.],
                  [ 0.,  0.,  0.,  0.],
                  [ 0.,  0.,  0.,  0.],
                  [ 0.,  0.,  0.,  0.]])

In that case, the covariance matrix object passes through several tests such as the Hermitian property of the two last dimensions, and the number of dimensions of the input array. If the input array does not pass the tests, a ValueError is raised. The input array can be any array-like object, as it is passed to the numpy.asarray() function, which also allow to turn the input array into a complex-valued numpy array. The input array is then cast into a CovarianceMatrix object with the view casting mechanism of numpy, and the attributes stats and stft are added to the object. You can also directly employ the view method of the numpy array, but you will have to set the stats and stft attributes manually.

coherence(kind='spectral_width', epsilon=1e-10)

Covariance-based coherence estimation.

The coherence is obtained from each covariance matrix eigenvalue distribution, calculated with the eigenvalues() method. For a given matrix \(\mathbf{C} \in \mathbb{C}^{N \times N}\), we obtain the eigenvalues \(\boldsymbol{\lambda} = \{\lambda_i\}_{i=1\ldots N}\), with \(\lambda_i \in [0, 1]\) normalized by the sum of the eigenvalues. The coherence is obtained from \(c = f(\boldsymbol{\lambda})\), with \(f\) being defined by the kind parameter:

• The spectral width is obtained with setting kind="spectral_width" , and returns the width \(\sigma\) of the eigenvalue distribution (obtained at each time and frequency) such as

  \[\sigma = \sum_{i=0}^n i \lambda_i\]

• The entropy is obtained with setting kind="entropy", and returns the entropy \(h\) of the eigenvalue distribution (obtained at each time and frequency) such as

  \[h = - \sum_{i=0}^n i \lambda_i \log(\lambda_i + \epsilon)\]

• The Shanon diversity index is obtained with setting kind="diversity", and returns the diversity index \(D\) of the eigenvalue distribution (obtained at each time and frequency) such as the exponential of the entropy:

  \[D = \exp(h + \epsilon)\]

In each case, \(\epsilon\) is a regularization parameter to avoid the logarithm of zero. The coherence is calculated for each time and frequency sample (if any), and the result is a coherence matrix of maximal shape (n_times, n_frequencies). The \(\log\) and \(\exp\) functions are the natural logarithm and exponential.

Parameters:

• kind (str, optional) -- The type of coherence, may be "spectral_width" (default), "entropy", or "diversity".

• epsilon (float, optional) -- The regularization parameter for the logarithm.

Returns:

numpy.ndarray -- The coherence of maximal shape (n_times, n_frequencies), depending on the input covariance matrix shape.

Raises:

ValueError -- If the covariance matrix is not Hermitian.

See also

eigenvalues()

eigenvalues(norm: ~typing.Callable = <function max>) → ndarray

Eigenvalue decomposition.

Given and Hermitian matrix \(\mathbf{C} \in \mathbb{C}^{N \times N}\), the eigenvalue decomposition is defined as

\[\mathbf{C} = \mathbf{U D U}^\dagger\]

where \(\mathbf{U} \in \mathbb{C}^{N \times N}\) is the unitary matrix of eigenvectors (which are not calculated here, see eigenvectors() for this purpose), \(\mathbf{U}^\dagger\) is the conjugate transpose of \(\mathbf{U}\), and \(\mathbf{D} \in \mathbb{R}^{N \times N}\) is a diagonal matrix of eigenvalues, as

\[\begin{split}\mathbf{D} = \pmatrix{\lambda_1 & 0 & \cdots & 0 \\ 0 & \lambda_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \lambda_N}\end{split}\]

with \(\lambda_i\) the eigenvalues. For convenience, the eigenvalues are returned as a vector such as

\[\boldsymbol\lambda = \pmatrix{\lambda_1, \lambda_2, \ldots, \lambda_N}\]

The eigenvalues are sorted in decreasing order and normalized by the maximum eigenvalue by default, but can be normalized by any function provided by numpy. Since the matrix \(\mathbf{C}\) is Hermitian by definition, the eigenvalues are real- and positive-valued. Also, the eigenvectors are orthogonal and normalized.

The eigenvalue decomposition is performed onto the two last dimensions of the CovarianceMatrix object. The function used for eigenvalue decomposition is numpy.linalg.eigvalsh(). It sassumes that the input matrix is 2D and hermitian, so the decomposition is performed onto the lower triangular part to save time.

Parameters:

norm (function, optional) -- The function used to normalize the eigenvalues. Can be numpy.max(), (default), any other functions provided by numpy, or a custom function. Note that the function must accept the axis keyword argument.

Returns:

numpy.ndarray -- The eigenvalues of maximal shape (n_times, n_frequencies, n_stations).

Notes

The eigenvalue decomposition is performed onto the two last dimensions of the CovarianceMatrix object. The matrices are first flattened with the flat() method, so the eigenvalues are calculated for each time and frequency sample. The eigenvalues are sorted in decreasing order, and normalized by the maximum eigenvalue by default, before being reshaped to the original shape of the covariance matrix. This maximizes the performance of the eigenvalue decomposition.

See also

eigenvectors(), numpy.linalg.eigvalsh()

Examples

Calculate the eigenvalues of the example covariance matrix:

>>> import covseisnet as cn
>>> import numpy as np
>>> cov = np.random.randn(3, 3, 3) + 1j * np.random.randn(3, 3, 3)
>>> cov = np.array([cov @ cov.T.conj() for cov in cov])
>>> cov = csn.CovarianceMatrix(cov)
>>> cov.eigenvalues()
    array([[1.        , 0.14577012, 0.02521345],
           [1.        , 0.13510247, 0.00369051],
           [1.        , 0.22129766, 0.0148769 ]])

eigenvectors(rank: int | tuple | slice | None = None, return_covariance: bool = False, weights: ndarray | None = None) → ndarray | CovarianceMatrix

Extract eigenvectors of given rank.

Given a Hermitian matrix \(\mathbf{C} \in \mathbb{C}^{N \times N}\), the eigenvector decomposition is defined as

\[\mathbf{C} = \mathbf{U D U}^\dagger\]

where \(\mathbf{U} \in \mathbb{C}^{N \times N}\) is the unitary matrix of eigenvectors, \(\mathbf{U}^\dagger\) is the conjugate transpose of \(\mathbf{U}\), and \(\mathbf{D} \in \mathbb{R}^{N \times N}\) is a diagonal matrix of eigenvalues. The eigenvectors are normalized and orthogonal, and the eigenvalues are real- and positive-valued. The eigenvectors are returned as a matrix

\[\mathbf{U} = \pmatrix{\mathbf{u}_1, \mathbf{u}_2, \ldots, \mathbf{u}_R}\]

with \(\mathbf{v}_i\) the eigenvectors, and \(R\) the maximum rank of the eigenvectors. The eigenvectors are sorted in decreasing order of the eigenvalues, and the eigenvectors are normalized. The function used for extracting eigenvectors is scipy.linalg.eigh(). It assumes that the input matrix is 2D and hermitian, so the decomposition is performed onto the lower triangular part.

If the covariance parameter is set to False (default), the eigenvectors are returned as a matrix of shape (n_times, n_frequencies, n_stations, rank) if the parameter covariance is False, else (n_times, n_frequencies, n_stations, n_stations), resulting from the outer product of the eigenvectors.

If the covariance parameter is set to True, the return value is a covariance matrix object with the outer product of the eigenvectors multiplied by the given weights \(\mathbf{w} \in \mathbb{R}^{R}\):

\[\tilde{\mathbf{C}} = \sum_{i=1}^R w_i \mathbf{u}_i \mathbf{u}_i^\dagger\]

The weights are the eigenvalues of the covariance matrix if weights is None (default), else the weights are the input weights. The weights are used to scale the eigenvectors before the outer product. In particular, if the weights are zeros and ones, this function can be used to apply spatial whitening to the covariance matrix.

Parameters:

• rank (int, tuple, slice, or None, optional) -- The rank of the eigenvectors. If None, all eigenvectors are returned. If a tuple or slice, the eigenvectors are sliced according to the tuple or slice. If an integer is given, the eigenvectors are sliced according to the integer. Default is None.

• return_covariance (bool, optional) -- If True, the outer product of the eigenvectors is returned as a covariance matrix object. Default is False.

• weights (numpy.ndarray, optional) -- The weights used to scale the eigenvectors before the outer product. If None, the eigenvalues are used as weights. Default is None.

Returns:

covseisnet.covariance.CovarianceMatrix -- The complex-valued eigenvector array of shape (n_times, n_frequencies, n_stations, rank) if covariance is False, else (n_times, n_frequencies, n_stations, n_stations).

Raises:

ValueError -- If the covariance matrix is not Hermitian.

See also

eigenvalues(), scipy.linalg.eigh()

property flat

Covariance matrices with flatten first dimensions.

The shape of the covariance matrix depend on the number of time windows and frequencies. The property flat allows to obtain as many \(N \times N\) covariance matrices as time and frequency samples.

Returns:

np.ndarray -- The covariance matrices in a shape (n_times * n_frequencies, n_traces, n_traces).

Example

>>> import covseisnet as csn
>>> import numpy as np
>>> cov = np.zeros((5, 4, 2, 2))
>>> cov = csn.CovarianceMatrix(cov)
>>> cov.shape
    (5, 4, 2, 2)
>>> c.flat.shape
    (20, 2, 2)

triu(**kwargs)

Extract upper triangular matrices.

This method is useful when calculating the cross-correlation matrix associated with the covariance matrix. Indeed, since the covariance matrix is Hermitian, the cross-correlation matrix is symmetric, so there is no need to calculate the lower triangular part.

The method triu() is applied to the flattened array, then reshaped to the original shape of the covariance matrix. The last dimension of the returned matrix is the number of upper triangular elements of the covariance matrix.

Parameters:

**kwargs (dict, optional) -- Keyword arguments passed to the numpy.triu_indices() function.

Returns:

CovarianceMatrix -- The upper triangular part of the covariance matrix, with a maximum shape of (n_times, n_frequencies, n_traces * (n_traces + 1) // 2).

Example

>>> import covseisnet as csn
>>> import numpy as np
>>> cov = np.zeros((5, 4, 2, 2))
>>> cov = csn.CovarianceMatrix(cov)
>>> cov.triu().shape
    (5, 4, 3)

twosided(axis: int = -3) → CovarianceMatrix

Get the full covariance spectrum.

Given that the covariance matrix is Hermitian, the full covariance matrix can be obtained by filling the negative frequencies with the complex conjugate of the positive frequencies. The method twosided() performs this operation.

The frequency axis is assumed to be the second axis of the covariance matrix. The function returns a new covariance matrix with the negative frequencies filled with the complex conjugate of the positive frequencies.

Parameters:

axis (int, optional) -- The frequency axis of the covariance matrix. Default is -3.

Returns:

CovarianceMatrix -- The full covariance matrix.

property is_hermitian: bool

Check if the covariance matrix is Hermitian.

Given a covariance matrix \(\mathbf{C} \in \mathbb{C}^{N \times N}\), the matrix is Hermitian if the matrix is equal to its conjugate transpose, that is

\[\mathbf{C} = \mathbf{C}^\dagger\]

to some extent provided by the tol parameter. This check is performed on all the covariance matrices of the object, e.g., for each time and frequency sample if any.

Parameters:

tol (float, optional) -- The tolerance for the comparison.

Returns:

bool -- True if the covariance matrix is Hermitian, False otherwise.

covseisnet.covariance.calculate_covariance_matrix(stream: NetworkStream, window_duration: float, average: int, window_step: float | None = None, average_step: int | None = None, whiten: str = 'none', water_level: float = 1e-10, **kwargs) → tuple[ndarray, ndarray, CovarianceMatrix]

Calculate covariance matrix.

The covariance matrix is calculated from the Fourier transform of the input stream. The covariance matrix is calculated for each time window and frequency, and averaged over a given number of windows. Please refer to the CovarianceMatrix class, rubric Mathematical definition for a detailed explanation of the covariance matrix. You may also refer to the paper [1] for deeper insights.

Covariance-level whitening

The whitening parameter can be used to normalize the covariance matrix directly, in addition to the trace-level whitening implementation in the covseisnet.stream.NetworkStream.whiten() method. The parameter can be set to "none" (default), "slice", or "window". The "none" option does not apply any whitening to the covariance matrix. The "slice" option normalizes the spectra \(u_{i,m}(f)\) by the mean of the absolute value of the spectra within the same group of time windows \(\{\tau_m\}_{m=1}^M\), so that

\[u_{i,m}(f) = \frac{u_{i,m}(f)}{\sum_{i=1}^M |u_{i,m}(f)|}\]

The "window" option normalizes the spectra \(u_{i,m}(f)\) by the absolute value of the spectra within the same time window \(\tau_m\) so that

\[u_{i,m}(f) = \frac{u_{i,m}(f)}{|u_{i,m}(f)|}\]

These additional whitening methods can be used in addition to the trave-level whiten() method to further whiten the covariance matrix.

Parameters:

• stream (NetworkStream) -- The input data stream.

• window_duration (float) -- The duration of the window used to calculate the short-time Fourier transform, in seconds.

• window_step (float) -- The step of the window used to calculate the short-time Fourier, in seconds. Default is None, which is equivalent to half of the window duration.

• average (int) -- The number of window used to estimate the sample covariance.

• average_step (int, optional) -- The sliding window step for covariance matrix calculation (in number of windows).

• whiten (str, optional) -- The type of whitening applied to the covariance matrix. Can be "none" (default), "slice", or "window". This parameter can be used in addition to the whiten() method to further whiten the covariance matrix.

• water_level (float, optional) -- The water level used to avoid division by zero in the "window" whitening method.

• **kwargs (dict, optional) -- Additional keyword arguments passed to the ShortTimeFourierTransform class.

Returns:

• numpy.ndarray -- The time vector of the beginning of each covariance window, expressed in matplotlib time.

• numpy.ndarray -- The frequency vector.

• CovarianceMatrix -- The complex covariance matrix, with a shape depending on the number of time windows and frequencies, maximum shape (n_times, n_frequencies, n_traces, n_traces).

Example

Calculate the covariance matrix of the example stream with 1 second windows averaged over 5 windows:

>>> import covseisnet as csn
>>> stream = csn.read()
>>> t, f, c = csn.calculate_covariance_matrix(
...     stream, window_duration_sec=1., average=5
... )
>>> print(c.shape)
    (27, 51, 3, 3)

References

[1]

Léonard Seydoux, Nikolai M. Shapiro, Julien De Rosny, Florent Brenguier, and Matthieu Landès. Detecting seismic activity with a covariance matrix analysis of data recorded on seismic arrays. Geophysical Journal International, 204(3):1430–1442, March 2016. URL: https://academic.oup.com/gji/article-lookup/doi/10.1093/gji/ggv531, doi:10.1093/gji/ggv531.