Application of the independent component analysis to the iKAGRA data

TL;DR

This paper demonstrates the first application of independent component analysis (ICA) to real gravitational wave detector data, showing it can effectively identify and enhance signals amidst environmental noise.

Contribution

The study introduces the use of ICA on actual iKAGRA data and compares two implementations, highlighting its potential for noise reduction in gravitational wave detection.

Findings

ICA successfully recovers injected signals with improved SNR.

Correlation-based ICA performs best when environmental noise is linear.

First application of ICA to real gravitational wave data.

Abstract

We apply the independent component analysis (ICA) to the real data from a gravitational wave detector for the first time. Specifically we use the iKAGRA data taken in April 2016, and calculate the correlations between the gravitational wave strain channel and 35 physical environmental channels. Using a couple of seismic channels which are found to be strongly correlated with the strain, we perform ICA. Injecting a sinusoidal continuous signal in the strain channel, we find that ICA recovers correct parameters with enhanced signal-to-noise ratio, which demonstrates usefulness of this method. Among the two implementations of ICA used here, we find the correlation method yields the optimal result for the case environmental noises act on the strain channel linearly.

Tables

Table 1. Table 1: Correlation between PEM channels and strain.

dataset	channel	correlation coefficient
strongly	PEM-EX_SEIS_Z_SENSINF_OUT16 (4724ch)	$-0.6409$
correlated	PEM-EY_SEIS_WE_SENSINF_OUT16 (4823ch)	$0.5892$
weakly	PEM-EX_SEIS_Z_SENSINF_OUT16 (4724ch)	$0.3078$
correlated	PEM-EY_SEIS_NS_SENSINF_OUT16 (4774ch)	$-0.2312$

Equations

$${\mbox{\boldmath${x}$}}(t)=\begin{pmatrix}x_{1}(t)\\ x_{2}(t)\end{pmatrix}=A{\mbox{\boldmath${s}$}}(t),~{}~{}~{}{\mbox{\boldmath${s}$}}(t)=\begin{pmatrix}s_{1}(t)\\ s_{2}(t)\end{pmatrix}=\begin{pmatrix}h(t)+n(t)\\ k(t)\end{pmatrix},$$

$$A=\left(\begin{array}{cc}{a}_{11}& a_{12}\\ 0& a_{22}\end{array}\right).$$

$${\mbox{\boldmath${y}$}}=W{\mbox{\boldmath${x}$}},$$

$$p_{y}({\mbox{\boldmath${y}$}})\equiv||W^{-1}||p_{x}({\mbox{\boldmath${x}$}}),$$

$$W=\left(\begin{array}{cc}{w}_{11}& w_{12}\\ 0& w_{22}\end{array}\right).$$

$$D[p_{y}({\mbox{\boldmath${y}$}});q({\mbox{\boldmath${y}$}})]=\int p_{y}({\mbox{\boldmath${y}$}})\ln\frac{p_{y}({\mbox{\boldmath${y}$}})}{q({\mbox{\boldmath${y}$}})}dy=E_{p_{y}}\left[\ln\frac{p_{y}({\mbox{\boldmath${y}$}})}{q({\mbox{\boldmath${y}$}})}\right].$$

$$r_1[s_1(t)]=\frac{1}{\sqrt{2\pi }\sigma }\exp \left[-\frac{1}{2{\sigma }^2}{\left(s_1(t)-h(t,\theta )\right)}^2\right],{\sigma }^2=K(0),$$

$$\left(\begin{array}{c}{y}_1(t)\\ y_2(t)\end{array}\right)=\left(\begin{array}{cc}{w}_{11}& w_{12}\\ 0& w_{22}\end{array}\right)\left(\begin{array}{c}{x}_1(t)\\ x_2(t)\end{array}\right)=\left(\begin{array}{c}{w}_{11}{x}_1(t)+w_{12}{x}_2(t)\\ w_{22}{x}_2(t)\end{array}\right),$$

$$⟨y_1(t)x_2(t)⟩=w_{11}⟨x_1(t)x_2(t)⟩+w_{12}⟨x_2^2(t)⟩=0.$$

$$w_{12}=-\frac{⟨x_1{x}_2⟩}{⟨x_2^2⟩}{w}_{11}.$$

$$\left(\begin{array}{c}{y}_1(t)\\ y_2(t)\\ y_3(t)\end{array}\right)=\left(\begin{array}{ccc}{w}_{11}& w_{12}& w_{13}\\ 0& w_{22}& w_{23}\\ 0& w_{32}& w_{33}\end{array}\right)\left(\begin{array}{c}{x}_1(t)\\ x_2(t)\\ x_3(t)\end{array}\right)=\left(\begin{array}{c}{w}_{11}{x}_1(t)+w_{12}{x}_2(t)+w_{13}{x}_3(t)\\ w_{22}{x}_2(t)+w_{23}{x}_3(t)\\ w_{32}{x}_2(t)+w_{33}{x}_3(t)\end{array}\right),$$

$$\left(\begin{array}{c}{x}_1(t)\\ x_2(t)\\ x_3(t)\end{array}\right)=\left(\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ 0& a_{22}& a_{23}\\ 0& a_{32}& a_{33}\end{array}\right)\left(\begin{array}{c}{s}_1(t)\\ s_2(t)\\ s_3(t)\end{array}\right)=\left(\begin{array}{c}{a}_{11}{s}_1(t)+a_{12}{s}_2(t)+a_{13}{s}_3(t)\\ a_{22}{s}_2(t)+a_{23}{s}_3(t)\\ a_{32}{s}_2(t)+a_{33}{s}_3(t)\end{array}\right).$$

$$y_3(t)={\tilde{x}}_3(t)\equiv \frac{x_3(t)}{\sqrt{⟨x_3^2⟩}}.$$

$$w_{23}=-\frac{⟨x_2{x}_3⟩}{⟨x_3^2⟩}{w}_{22}.$$

$$y_2(t)={\tilde{x}}_2(t)\equiv \frac{x_2^{\prime }(t)}{\sqrt{⟨x_2^{{}^{\prime }2}⟩}},x_2^{\prime }(t)\equiv x_2(t)-\frac{⟨x_2{x}_3⟩}{⟨x_3^2⟩}{x}_3(t).$$

$$y_1(t)=x_1(t)-⟨x_1{\tilde{x}}_2⟩{\tilde{x}}_2(t)-⟨x_1{\tilde{x}}_3⟩{\tilde{x}}_3(t),$$

$$\tilde{{\mbox{\boldmath${x}$}}}(t)=\Lambda^{-1/2}\ {}^{t}\Gamma{\mbox{\boldmath${x}$}}=\Lambda^{-1/2}\ {}^{t}\Gamma A{\mbox{\boldmath${s}$}}\equiv\tilde{A}{\mbox{\boldmath${s}$}},$$

$$\langle\tilde{{\mbox{\boldmath${x}$}}}(t)\,^{t}\tilde{{\mbox{\boldmath${x}$}}}(t)\rangle=\langle\tilde{A}{\mbox{\boldmath${s}$}}\,^{t}(\tilde{A}{\mbox{\boldmath${s}$}})\rangle=\tilde{A}\langle{\mbox{\boldmath${s}$}}\,^{t}{\mbox{\boldmath${s}$}}\rangle^{t}\tilde{A}=\tilde{A}\ ^{t}\tilde{A}=E.$$

$$q({\mbox{\boldmath${y}$}})=\tilde{p}_{y}({\mbox{\boldmath${y}$}})\equiv\prod_{i}\tilde{p}_{i}(y_{i}),~{}~{}~{}\tilde{p}_{i}(y_{i})=\int p_{y}({\mbox{\boldmath${y}$}})dy_{1}...dy_{i-1}dy_{i+1}...,$$

$$L_{\tilde{p}}(W)=D[p_{y}({\mbox{\boldmath${y}$}});\tilde{p}_{y}({\mbox{\boldmath${y}$}})]=-H[{\mbox{\boldmath${x}$}}]-\ln||W||+\sum_{i}H_{i}[y_{i}],$$

$$J[y_i]\equiv H[\nu ]-H[y_i],$$

$$L({\mbox{\boldmath${w}$}}_{i})=\left\{E[G(y_{i})]-E[G(\nu)]\right\}^{2}-\beta\left[|{\mbox{\boldmath${w}$}}_{i}|^{2}-1\right],$$

$$E[\tilde{{\mbox{\boldmath${x}$}}}g(^{t}{\mbox{\boldmath${w}$}}_{i}\tilde{{\mbox{\boldmath${x}$}}})]-\beta{\mbox{\boldmath${w}$}}_{i}=0,$$

$$s(t)=A\sin (2\pi ft),$$

$$\begin{array}{rl}{x}_1(t)& =h(t)+n(t)+\sum _{j=2}^n{a}_{1j}{s}_j(t),\\ x_i(t)& =\sum _{j=2}^n{a}_{ij}{s}_j(t).\end{array}$$

$${\tilde{x}}_1=x_1+\sum _{j=2}^n{b}_{1j}{x}_j,$$

$$\begin{array}{rl}⟨{\tilde{x}}_1(t)x_i(t)⟩& =⟨\left(h(t)+\sum _{j=2}^n{a}_{1j}{s}_j(t)+\sum _{j=2}^n{b}_{1j}\sum _{l=2}^n{a}_{jl}{s}_l(t)\right)\sum _{k=2}^n{a}_{ik}{s}_k(t)⟩\\ & =\sum _{j=2}^n{a}_{1j}{a}_{ij}+\sum _{j=2}^n{b}_{1j}\sum _{k=2}^n{a}_{jk}{a}_{ik}=0,\end{array}$$

$$b_{1j}=\sum _{i=2}^n{a}_{1i}{a}_{ij}^{-1}$$

$$\begin{array}{rl}{\tilde{x}}_1& =x_1+\sum _{i=2}^n{b}_{1i}{x}_i\\ & =h(t)+n(t)+\sum _{i=2}^n{a}_{1i}{s}_i(t)+\sum _{i=2}^n\sum _{j=2}^n\sum _{l=2}^n{a}_{1j}{a}_{ji}^{-1}{a}_{il}{s}_l(t)\\ & =h(t)+n(t)+\sum _{i=2}^n{a}_{1i}{s}_i(t)-\sum _{j=2}^n\sum _{l=2}^n{a}_{1j}{\delta }_{jl}{s}_l(t)\\ & =h(t)+n(t)+\sum _{i=2}^n{a}_{1i}{s}_i(t)-\sum _{j=2}^n{a}_{1j}{s}_j(t)=h(t)+n(t),\end{array}$$

Full text

1 Introduction

Ever since Einstein found the existence of a gravitational wave solution in his theory of general relativity in 1916, it took exactly a century for mankind to succeed in its direct detection. This delay is primarily due to the fact that the gravitational force is an exceedingly weak force compared with other interactions.

The first detection of a gravitational wave by the advanced Laser Interferometer Gravitational wave Observatory (aLIGO) [1] brought a great impact on science and told the beginning of gravitational wave astronomy. Following aLIGO and advanced Virgo, the large-scale cryogenic gravitational wave telescope (LCGT) now known as KAGRA, has been constructed in Kamioka, Japan [2]. KAGRA will play very important roles in the international network of gravitational wave detection by measuring the number of polarization property, which is indispensable to prove the general relativity [3], and by improving the sky localization of each event significantly [4]. As the first underground and cryogenic detector, it will also provide important information to the third-generation detectors.

Because gravity is the weakest force among the four elementary interactions, gravitational waves have high penetrating power. Therefore, unlike electromagnetic waves, they can propagate without being influenced by interstellar medium. In the same way, they enable us to see deep inside dense matter, such as the core of neutron stars, and bring information that electromagnetic waves cannot. On the other hand, due to this property, gravitational wave signals tend to be quite small and its detection becomes very difficult. Thus, it is important to develop methods for extraction of these tiny signals. There are a number of methods which extract signal out of large noises such as matched filtering [5], which yields an optimal result if (and only if) underlying noise is Gaussian distributed. However, the problem is not so simple, as it is known that there exist non-Gaussian noises in real data. They decrease the performance of the analysis methods assuming Gaussianity of the noises. What is worse, these noises may be mistaken for true signals and increase the false alarm probability. Thus, it is necessary to deal with non-Gaussianity properly as stressed in [6]. Characterization, mitigation, and even subtraction of these noises in gravitational wave detector outputs have been extensively studied in the literature. The standard way including pre-data conditioning (whitening, band-passing), line-removal, and $\chi^{2}$ veto are well overviewed in [7]. Many of recent works demonstrate performance of Deep Neural Networks [8, 9, 10, 11, 12, 13], but see also [14].

In this situation, independent component analysis (ICA) [15, 16, 17] , which seperates mixed signal components, occupies a unique position among methods of signal processing because it makes use of non-Gaussianity of signals and noises instead of treating it as an obstacle. ICA has been used in various fields in astronomy, e.g., [18, 19, 20, 21, 22, 23, 24, 25, 26]. For example, reference [18] demonstrated ICA (EFICA and WASOBI) performance on simulated data mimicking two gravitational wave interferometer outputs. The current paper, on the other hand, demonstrates it using real gravitational wave strain data from the iKAGRA detector and multiple real auxiliary channels that recorded status of the detector. ICA can separate various components obeying non-Gaussian distributions, so that it can remove (part of) non-Gaussian noises from strain data that records gravitational wave signals. Then the strain channel would consist of the real signal and (nearly) Gaussian noises. Therefore ICA can support conventional matched filter technique as a non-Gaussian noise subtraction scheme. In addition, ICA can be used even in the cases noises are nonlinearly coupled to the strain channel as demonstrated in [27].

In this paper, we use the correlation method [27] (or the Gram-Schmidt orthogonalization method in the case of multiple channels in general) and FastICA [31]. The former method, although conceptually different in derivation, has practically the same expression as the Wiener filtering [33, 34, 35, 36, 37, 38] which has been used to analyze Caltech 40m and LIGO data, and quite remarkable success was recently reported [38] by using witness sensors including voltage monitors of analogue electronics for power main and photodiodes that monitor beam motion and its size for beam jitter. We report the results of application of these two different ICA methods to the iKAGRA data and discuss its usefulness in gravitational wave data analysis. The paper is organized as follows. In § 2, we introduce ICA in the simplest case where only one environmental channel is incorporated to the strain channel and review analytic formulas of correlation method obtained in our previous paper [27]. Then we extend this method to the case where two different environmental channels are concerned. We also introduce FastICA which is formulated in a different way. In § 3, we present our application of ICA to the iKAGRA data with injected artificial continuous signal. Then we discuss the result focusing on the difference of the two methods in § 4. We argue that for the current setup where noises measured by the environmental channels affect the strain linearly and additively, the what we call correlation method yields the optimal result. The final section § 5 is devoted to conclusion.

2 Independent Component Analysis (ICA)

As is seen in our previous paper [6], signal detection under non-Gaussian noises is much more involved than the case with Gaussian noises since the optimal statistic has much complicated forms. ICA is an attractive method of signal processing because it makes use of non-Gaussian nature of the signals [15, 16, 17] (see [28, 29] for textbooks). We here introduce two methods of ICA as a method of non-Gasussian noise subtraction.

Basically, this method assumes only statistical independence between the signal and noises, and does not impose any other conditions on their distributions. However, a simpler formulation can be achieved by using physical information of gravitational wave detection as expressed in [27]. Following [27], we first formulate the subtraction of non-Gaussian noise in the gravitational wave detection for the case where noise is linearly coupled to the strain. Then we introduce analytic formulas of ICA for this case, which we call the correlation method. Previously, this was two component analysis in [27], but we here developed a multiple component version for combining different environmental channels.

On the other hand, there is a robust formulation which does not incorporate any information of the concerned system, which is called FastICA [31]. We also introduce this method in this section and apply it in our analysis as a comparison.

2.1 Removing non-Gaussian noises

In this paper, we consider the following simple problem as a first step to test applicability of ICA for detection of GWs. Let us consider the case where we have two detector outputs, $x_{1}(t)$ and $x_{2}(t)$ ($t$ stands for time). The former is the output from the laser interferometer, namely, the strain channel, and the latter is an environmental channel such as an output of a seismograph. We wish to separate gravitational wave signal $h(t)$ and non-Gaussian noise $k(t)$ using the data of {}^{t}{\mbox{\boldmath{x}}}(t)=(x_{1}(t),x_{2}(t)).

As the simplest case we assume that there is a linear relation between the outputs and the sources:

[TABLE]

where $A$ is assumed to be a time independent matrix. Since the output of a laser interferometer, of course, suffers from Gaussian noise $n(t)$, we can regard $s_{1}(t)=h(t)+n(t)$ as an original signal. Note that non-Gaussian noise $k(t)$ can contain any Gaussian noise as a part of it. Thus, we have not added any Gaussian noise to $s_{2}(t)$ explicitly.

Since the gravitational wave is so weak that it will not affect any environmental meters such as a seismograph, one may set $A$ as

[TABLE]

The aim of ICA is to find a linear transformation

[TABLE]

such that two components of the transformed variables ${y}$ are statistically independent of each other. Here the distribution of ${y}$, p_{y}({\mbox{\boldmath{y}}}), is constructed from the observed distribution function of ${x}$, p_{x}({\mbox{\boldmath{x}}}), through the transformation (3) as

[TABLE]

where $||X||$ denotes the determinant of matrix $X$. Thanks to the assumption (2), the matrix $W$ also takes a form

[TABLE]

However, since we do not know all the component of $A$, we attempt to determine $W$ to be $A^{-1}$ in such a way that the components of ${y}$, $y_{1}(t)$ and $y_{2}(t)$ to be statistically independent as much as possible. In [27] this was achieved by using the the Kullback-Leibler divergence [30], which represents a distance in the space of statistical distribution functionals. It is defined between two arbitrary PDFs, e.g., p_{y}({\mbox{\boldmath{y}}}) and q({\mbox{\boldmath{y}}}) as

[TABLE]

Here $E_{p_{y}}[\cdot]$ denotes an expectation value with respect to a PDF $p_{y}$. Then we can obtain mutually independent variables ${y}$ by minimizing a cost function L_{q}(W)\equiv D[p_{y}({\mbox{\boldmath{y}}});q({\mbox{\boldmath{y}}})], where q({\mbox{\boldmath{y}}})=q(y_{1})q(y_{2}) is an appropriately chosen distribution function.

The most proper choice of q({\mbox{\boldmath{y}}}) is obviously the true distribution function of independent source variables ${s}$, r({\mbox{\boldmath{s}}})=r_{1}[s_{1}(t)]r_{2}[s_{2}(t)], which is not known a priori. Because $n(t)$ is a Gaussian with vanishing mean in this simple setup, its statistical property is entirely characterized by the two-point correlation function $K(t-t^{\prime})=\langle n(t)n(t^{\prime})\rangle$. Then the marginal distribution function of $s_{1}(t)$ is given by

[TABLE]

where $h(t,\theta)$ is the actual waveform of gravitational radiation emitted from some source, where $\theta$ collectively denotes parameters of the source. On the other hand, we do not specify the PDF of $k(t)$, $r_{2}(s_{2})$, except that it is a super-Gaussian distribution such as a Poisson distribution with a larger tail than Gaussian. We show, however, that we can obtain the matrix $W$ easily for our particular problem with $a_{21}=w_{21}=0$ as we see below.

2.2 Correlation method

From now on we replace the ensemble average $E[\cdot]$ by temporal average of observed values of ${x}$ which we denote by brackets. For the true distribution r({\mbox{\boldmath{y}}}), minimization of the cost function $L_{r}(W)$ results in decorrelating $y_{1}$ and $y_{2}$ [27], i.e. $\langle y_{1}(t)y_{2}(t)\rangle=0$.

From

[TABLE]

it is equivalent to requiring

[TABLE]

We therefore obtain

[TABLE]

Since ICA does not uniquely determine the overall factor of ${y}$ by nature, this relation suffices for our purpose to determine $y_{1}$. These are what we calculated in our previous paper [27] using the Kullback-Leibler divergence.

Here we develop a multiple component method for further analysis and we apply it in § 3.2.2. For three components, {\mbox{\boldmath{y}}}(t) and {\mbox{\boldmath{x}}}(t) become

[TABLE]

and

[TABLE]

In this case also, the minimization of cost function results in decorrelating ${y}$, $\langle y_{1}y_{2}\rangle=\langle y_{2}y_{3}\rangle=\langle y_{3}y_{1}\rangle=0$. This is achieved by the analogy of the Gram–Schmidt process which is a method for orthonormalising a set of vectors, and it can be extended to the case where there are more than three components.

Because of the gauge degree of freedom, we can take $w_{32}=0$ without loss of generality and choose

[TABLE]

We first require $\langle y_{2}(t)y_{3}(t)\rangle=\langle y_{2}(t)x_{3}(t)\rangle=0$. This gives following relation,

[TABLE]

Based on this, we can choose

[TABLE]

If we take

[TABLE]

$\langle y_{1}(t)y_{2}(t)\rangle=\langle y_{2}(t)y_{3}(t)\rangle=\langle y_{3}(t)y_{1}(t)\rangle=0$ is satisfied. Note that Eq. (16) is symmetrical with respect to the permutation of $x_{2}(t)$ and $x_{3}(t)$.

Thus we can observe that the correlation method of ICA shown here is equivalent to the instantaneous Wiener filtering111The Wiener filtering adopted in [33, 34, 35, 36, 37] takes into account the time delay in transfer functions. We can easily incorporate it to our analysis, too, as already demonstrated in [27]., and this is due to the particular character of our problem that only the strain channel is sensitive to the gravitational wave signal with $a_{i1}=0\ (i\neq 1)$ in our linear model.

2.3 FastICA method

Next, we introduce another method to obtain a matrix $W$ called FastICA [31] which can be easily implemented even when {\mbox{\boldmath{x}}}(t)=A{\mbox{\boldmath{s}}}(t) has more than two components. Note that this method can be applied to various cases of signal separation other than the case formulated in §2.1.

In this method, assuming that each component, $s_{i}(t)$, of source vector {\mbox{\boldmath{s}}}(t) is properly normalized with vanishing mean, we first apply whitening to the detector outputs {\mbox{\boldmath{x}}}(t) and take the dispersion of each source $s_{i}(t)$ to be unity. This is achieved in the following way. First let the normalized eigenvector and corresponding eigenvalue of a matrix \langle{\mbox{\boldmath{x}}}\,^{t}{\mbox{\boldmath{x}}}\rangle be {\mbox{\boldmath{c}}}_{i} and $\lambda_{i}$, respectively $(i=1,2,...)$, and define a matrix $\Gamma$ by \Gamma=({\mbox{\boldmath{c}}}_{1},{\mbox{\boldmath{c}}}_{2},{\mbox{\boldmath{c}}}_{3},...), and $\Lambda^{-1/2}$ by $\Lambda^{-1/2}=\mathrm{diag}(\lambda_{1}^{-1/2},\lambda_{2}^{-1/2},...)$. Then the whitened variable \tilde{{\mbox{\boldmath{x}}}}(t) is defined by

[TABLE]

which satisfies

[TABLE]

Here we have used the statistical independence of each component of the normalized source term $s_{i}$. This means that the matrix $\tilde{A}$ is an orthogonal matrix and that $W$ may be identified with ${}^{t}\tilde{A}$ for whitened output data \tilde{{\mbox{\boldmath{x}}}}. Thus we may restrict $W$ to be an orthogonal matrix, too, after appropriate whitening222Note that this procedure is also called as sphering and has nothing to do with the whitening of strain data in frequency domain..

We here choose q({\mbox{\boldmath{y}}}) as a product of marginal distributions,

[TABLE]

since p_{y}({\mbox{\boldmath{y}}})=\tilde{p}_{y}({\mbox{\boldmath{y}}}) is the condition for statistical independence of the variables ${y}$. Then, the cost function defined in terms of the Kullback-Leibler divergence reads

[TABLE]

where H[{\mbox{\boldmath{x}}}]\equiv-\int d{\mbox{\boldmath{x}}}p_{x}({\mbox{\boldmath{x}}})\ln p_{x}({\mbox{\boldmath{x}}}) is the entropy of the distribution of ${x}$, and $H_{i}[y_{i}]\equiv-\int dy_{i}\tilde{p}_{i}(y_{i})\ln\tilde{p}_{i}(y_{i})$ is the entropy of the marginal distribution of $y_{i}$. When $W$ is an orthonormal matrix, only the last term matters to determine $W$. Hence minimization of the cost function for $W$ is achieved by minimizing entropy of the marginal distribution of each variable. This is the spirit of the FastICA method. It has been proposed to maximize the negentropy defined by

[TABLE]

which is a positive semi-definite quantity, instead of the entropy itself. Here $\nu$ is a random Gaussian variable with vanishing mean and unit variance.

In order to achieve easier implementation of the method, however, we minimize a simpler cost function L({\mbox{\boldmath{w}}}_{i}) for each row vector {\mbox{\boldmath{w}}}_{i} constituting the matrix $W$ as W\equiv({\mbox{\boldmath{w}}}_{1},{\mbox{\boldmath{w}}}_{2},...). Since $W$ is an orthogonal matrix now, we find |{\mbox{\boldmath{w}}}_{i}|^{2}=1, so the cost function may be defined as

[TABLE]

where $G$ is an appropriate nonquadratic function and $\beta$ is a Lagrange multiplier. Minimization of Eq. (22) corresponds to solving the following equation:

[TABLE]

where $g(y)=G^{\prime}{}(y)$. FastICA solves for this equation starting from an arbitrary initial choice of {\mbox{\boldmath{w}}}_{i} in terms of the Newton method.

3 Analysis of iKAGRA data

The initial engineering run of KAGRA without the cryogenic system was done in March and April, 2016 [32]. From the results of many time series data that we analyzed, we report those of two datasets of 224 second long. One starts from 20:15:11 UTC on April 14, 2016. The other starts from 01:01:35 UTC on April 17, 2016. For each dataset, we calculated Pearson’s correlation between the strain channel and each of 35 physical environmental monitor (PEM) channels. We found that almost all these channels in the latter (former) data set strongly (weakly) correlated with the strain channel. We call the latter (former) the strongly (weakly) correlated data. The amplitude spectrum density (ASD) of the strain channel for each data set is depicted in Fig. 1.

We chose two channels which showed large correlation with the strain channel for each dataset. Those channels are listed in Table 1.

For both datasets, 4724ch had the largest correlation with the strain. This channel is the output of the seismograph that observes vertical vibration installed at the end of the X arm. Both 4774ch and 4823ch are the outputs of the seismographs installed at the end of the Y arm, and they observe horizontal vibration orthogonal to each other.

We have made mock strain data injecting sinusoidal continuous waves

[TABLE]

to the strain channel and applied two methods of ICA, which were introduced in the previous section, to this mock data and those environmental channels.

We utilized the python implementation of FastICA from scikit-learn333https://scikit-learn.org/stable/. We found that results often depend on initial conditions where the Newton method is started. To mitigate this, we parallelly generated at most thirty realizations and chose one which gives the highest SNR.

3.1 Global performance

First, we analyze how the signal-to-noise ratio (SNR) changes before and after noise separation by ICA for mock data with varying frequencies $f$. We performed matched filter (MF) analysis to both the raw mock strain data and the noise-removed data in terms of the two methods of ICA using 4724ch as an environmental channel. For various $f$ of injected signal (24), we calculated SNR by applying MF with the same frequency as the injected signal. We simultaneously plot the results against the data before and after ICA to assess the global performance of ICA. For strongly correlated data, the results are shown in Fig. 2.

In this dataset, strain had larger amplitude than the other dataset, and we set $A=9\times 10^{-10}$. As one can see from Fig. 2, SNRs are homogeneously enhanced by ICA for $f\gtrsim 0.1$Hz. The correlation method enhances the SNR more than FastICA. However, there are anomalous peaks at frequencies 0.01Hz and 0.04Hz. As shown in Fig. 1(a), even in the absence of injection the strain channel has large amplitude at these frequencies, which is predominantly contributed by seismic noises. We also found that their oscillation phases are more or less stable during the time period we analyzed. Such noises are difficult to be distinguished from our sinusoidal signal waveform and hence yield large SNR of mock strain as shown in Fig. 2. This, however, indicates that by removing contribution of the noises, the SNR can possibly be reduced rather than enhanced provided the injected signal is moderate. This is actually realized in the analysis based on the correlation method as seen in Fig. 2.

On the other hand, in the case of FastICA, the reduction of SNR is not seen. This is solely due to our implementation, which tries to increase the SNR as much as possible as mentioned before. In that sense, around the 0.01Hz and 0.04Hz peaks, blue line in Fig. 2 corresponds to the SNR of the separated noise.

Apart from these low frequencies contaminated by seismic noises, we find that ICA improves SNR significantly throughout the entire frequency range with $f\gtrsim 0.1$Hz. However, based on these considerations, it is deduced that ICA works even near the peak due to seismic noise.

For the weakly correlated data, the results are shown in Fig. 3. The amplitude of strain at this time period is moderate, and we set $A=3\times 10^{-11}$. As is seen in Fig. 3, the SNR of the data with ICA is higher than the mock data in several frequency ranges. Comparing FastICA with the correlation method, the correlation method has fewer frequencies where the SNR falls below that of mock data.

As for the weakly correlated data, 4774ch had the second highest correlation with strain. If we use 4774ch instead of 4724ch as the environmental data, the result changes as shown in Fig. 4. Compared with the case 4724ch is used (Fig. 3), the frequency region where the SNR rises is different. As a whole the improvement of SNR is less significant, which is a natural result considering that the correlation coefficient of 4774ch is smaller than that of 4724ch.

3.2 Parameter estimation for strongly correlated data

3.2.1 Two channels ICA

Next, we perform parameter estimation using the strongly correlated data to examine whether ICA can recover correct parameters of injected signals. We injected the sinusoidal waveform in Eq. (24) with $f=0.125$Hz and $A=1.3\times 10^{-9}$. We applied MF analysis to search for the frequency with the highest SNR which corresponds to the maximum likelihood estimation of the parameter. We compare how the result of parameter estimation changes before and after ICA and how much the SNR changes.

Figure 5 depicts the SNR before and after applying ICA.

In this case, we can see the effect of seismic noise directly. By ICA with 4724ch, SNR at $f\sim 0.01$Hz is reduced and that at the injected frequency $f=0.125$Hz is successfully enhanced. From this result, we deduce that 4724ch is highly correlated to the 0.01Hz peak. On the other hand, the peak of $0.04$Hz is still higher, which turned out to be correlated to 4823ch which had the second largest correlation with the strain, as we will see below.

3.2.2 Multiple channels ICA

・**Correlation method

**As is seen § 3.1, the correlation method shows more stable performance than FastICA, although it is much simpler. This method can be generalized to multi-channel analysis. As a first step to multi-channel analysis, here we investigate the effectiveness of three components analysis, which we developed in §2.2, including two PEM channels which strongly correlated to the strain. For this purpose we have used the mock data including the same signal waveform as in the previous subsection, and applied the three components correlation method to this mock data, 4724ch and 4823ch. The result is depicted in Fig. 6. We simultaneously plotted the results of two-component analysis in which we used 4724ch and 4823ch respectively.

The green and black lines correspond to the cases where noises are removed using one PEM channel. While the 0.01Hz peak was reduced by using 4724ch, the 0.04Hz peak was reduced by using 4823ch. However, both peaks cannot be reduced when we use only one PEM channel. The data with ICA using two PEM channels (cyan line) has much higher SNR than the data with ICA using only one PEM channel. In addition, we successfully reduced both 0.01Hz peak and 0.04Hz peak. This result suggests that by combining many environmental channels we can effectively remove noises with various characteristic frequencies.

　

・**FastICA

**As explained in § 2.3, FastICA can be easily implemented even when there are more than two components. We applied FastICA to the mock data, 4724ch and 4823ch simultaneously. Here, mock data included the same sinusoidal signal as in the previous section. The result is shown in Fig. 7.

As compared to Fig. 5, SNR at fiducial frequency is much higher than the case where only 4724ch was used. In addition, its value is close to the that for three component correlation method (10.60 for FastICA, 10.89 for the correlation method). This result suggests that the use of multiple environmental channels can also enhance the effect of FastICA noise separation. However, compared to the case of the three components correlation method, we may have to make several trials of 3ch FastICA to obtain the best result. This indicates that correlation method is more effective than FastICA for this dataset.

3.3 Parameter estimation for weakly correlated data

We also perform parameter estimation for weakly correlated data. Here, we used 4724ch as an environmental channel. From Fig.3, ICA using 4724ch is most effective for $f=0.227$Hz with this dataset. We injected sinusoidal wave signal with $f=0.227$Hz and $A=3\times 10^{-11}$. Again, we applied MF to search for the frequency with the highest SNR. The result is depicted in Fig. 8.

　The red line represents SNR calculated with the raw mock strain. The green and blue lines correspond to the noise-removed strain by the correlation method and FastICA, respectively. An enlarged figure of the fiducial ($f=0.227$Hz) area is shown in Fig. 8 (b).

As one can see, in the case of the raw mock strain, the position of the SNR peak deviates from the fiducial one. On the other hand, after applying ICA, the SNR is increased and the peak is found at the correct frequency.

　Next, we applied multiple channels ICA for this data using 4724ch and 4774ch. Here, we used the correlation method. Figure 9 depicts the results of analysis. The enlarged figure of the fiducial area is shown in Fig. 9 (b).

The green and black lines correspond to the data with ICA using one PEM channel. When using only 4774ch, enhancement of SNR is small and still the SNR peak deviates from the fiducial frequency. However, the data with ICA using two channels has slightly higher SNR at the correct frequency than the data with ICA using only 4724ch. This result for the weakly correlated data also supports our expectation that the effect of ICA can be enhanced by combining many environmental channels.

4 Discussion

In § 3, we show the performance of ICA as a method of non-Gaussian noise subtraction in GW data. Both ICA methods, namely the correlation method and FastICA subtract the portion of seismic noise. However, the correlation method shows better performance than that of FastICA in most cases. In this section, we consider the reason why this difference appears.

We here use the same notation as in § 2, $x_{1}(t)$ to be the strain channel, $x_{i}(t)\ (i=2,...,n)$ to be the other environmental channels. As we already discussed in § 2.1, the data of these channels can be written in the following form

[TABLE]

Here $s_{i}(t)\ (i=2,...,n)$ are environmental noises that can be measured by the PEM channels $x_{i}(t)$, and $n(t)$ collectively represents noises of the strain channel to which these PEM channels are insensitive. Let us transform $x_{1}$ as

[TABLE]

in order to satisfy $\langle\tilde{x}_{1}x_{i}\rangle=0\ (i=2,...,n)$. This condition can be expanded as

[TABLE]

where we have used $\langle s_{i}s_{j}\rangle=\delta_{ij}$ as in § 2.3. From this equation, we obtain

[TABLE]

with $j=2,...,n$. Note that $a_{ij}^{-1}$ is the inverse matrix of $a_{ij}\ (i,j=2,...,n)$, which is an $(n-1)\times(n-1)$ partial matrix of the mixing matrix $A=(a_{ij})_{1\leq i,j\leq n}$. By substituting this into eq.(26), we obtain

[TABLE]

which shows that all environmental noises $s_{i}$, which are measurable by PEM channels, are removed from $\tilde{x}_{1}$ just by imposing $\langle\tilde{x}_{1}x_{i}\rangle=0\ (i=2,...,n)$. In other words, when we consider auxilialy channels which are not sensitive to gravitational waves, and their target noises affects the strain linearly and additively, we can obtain independent component $h(t)$ by the transformation (26) which eliminates the two-point correlation between the strain and those channels. Although we have not given concrete expression of the transformation (26), $\langle\tilde{x}_{1}x_{i}\rangle=0$ is naturally achieved by the correlation method which is analogous to the Gram-Schmidt orthogonalization. Thus, we find that the correlation method is the optimal filter for linearly coupled noise with $a_{i1}=0$.

On the other hand, FastICA maximizes negentropy after the whitening which makes $\langle\tilde{x}_{1}x_{i}\rangle=0$ without using the condition $a_{i1}=0\ (i\neq 1)$. Since we do not recover this property in general even after maximizing the negentropy, FastICA tends to show less enhancement of SNR than that of the optimal correlation method for most cases in our analysis. This illustrates the importance of incorporating characteristic features of the system as much as possible before applying ICA.

However, the above discussion is only the case where linearly coupled noise with $a_{i1}=0$ is concerned. As for real observational data, there should be much more complicated mixing such as nonlinear coupling of the noise, and we might need a formulation of ICA which treats general mixing of signals. In that sense, it is noteworthy that FastICA, which is formulated without any assumption like $a_{i1}=0$, also shows enhancement of the SNR and improvement of the performance with multi-environmental channels to some extent.

5 Conclusion

In the present paper, we have demonstrated usefulness of ICA in gravitational wave data analysis in application to the iKAGRA strain and environmental channels. Assuming continuous waves as input signals, we have shown that ICA can enhance SNR in particular when the strain channel has large correlation with environmental ones. Moreover, we have shown that ICA can correctly recover input frequencies in parameter estimation. We have also found that combining multiple environmental channels can enhance the effect of ICA to improve SNR.

There are, however, a number of limitations in the analysis presented here because the iKAGRA data contains more low frequency mode than wanted due to the simplified vibration isolation system compared with the full designed specification which will be realized with bKAGRA [39], and iKAGRA configuration was not equipped with environmental monitors that measure hecto-Hertz frequencies. Hence we had to concentrate on relatively lower frequency components as the first step of application of ICA to real data analysis of laser interferometers.

Another limitation is that we have restricted to the case all the environmental noises that can be measured by the PEM channels under consideration act on the strain channel linearly and additively, without incorporating nonlinear couplings. In this particular situation, we have shown that the Gram-Schmidt decorrelation approach, or the instantaneous Wiener filtering which we dubbed the correlation method, gives the optimal result of environmental noise removal as an implementation of ICA. However, ICA can be used even in the cases noises of different origin are nonlinearly coupled to affect the strain channel as demonstrated in [27]. This is one of the merits of ICA absent in other methods. We could not perform such an analysis here due to the limitation of available PEM channels. We plan to return to this issue when the full cryogenic configuration of bKAGRA starts operation with more PEM channels.

**Acknowledgements

**This work was supported by MEXT, JSPS Leading-edge Research Infrastructure Program, JSPS Grant-in-Aid for Specially Promoted Research 26000005, JSPS Grant-in-Aid for Scientific Research on Innovative Areas 2905: JP17H06358, JP17H06361 and JP17H06364, JSPS Core-to-Core Program A. Advanced Research Networks, JSPS Grant-in-Aid for Scientific Research (S) 17H06133, the joint research program of the Institute for Cosmic Ray Research, University of Tokyo in Japan, National Research Foundation (NRF) and Computing Infrastructure Project of KISTI-GSDC in Korea, Academia Sinica (AS), AS Grid Center (ASGC) and the Ministry of Science and Technology (MoST) in Taiwan under grants including AS-CDA-105-M06, the LIGO project, and the Virgo project. This paper carries JGW Document Number JGW-P1910218.

　Jun’ya Kume is supported by a research program of the Leading Graduate Course for Frontiers of Mathematical Sciences and Physics (FMSP). This work was partially supported by JSPS KAKENHI, Grant-in-Aid for Scientific Research No. 15H02082 (Jun’ichi Yokoyama, Yosuke Itoh, Toyokazu Sekiguchi).