{'fraction_non_alphanumeric': 0.02942028519007307, 'fraction_numerical': 0.01721011577129001, 'mean_word_length': 4.649975751697381, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Reduction of transient noise artifacts in gravitational-wave data using deep learning

Excess transient noise artifacts, or glitches impact the data quality of ground-based gravitational-wave (GW) detectors and impair the detection of signals produced by astrophysical sources. Mitigation of glitches is crucial for improving GW signal detectability. However, glitches are the product of short-lived linear and non-linear couplings among the interrelated detector-control systems that include optic alignment systems and mitigation of seismic disturbances, generally making it difficult to model noise couplings. Hence, typically time periods containing glitches are vetoed to mitigate the effect of glitches in GW searches at the cost of reduction of the analyzable data. To increase the available data period and improve the detectability for both model and unmodeled GW signals, we present a new machine learning based method which uses on-site sensors/system-controls monitoring the instrumental and environmental conditions to model noise couplings and subtract glitches in GW detector data. We find that our method reduces 20-70% of the excess power in the data due to glitches. By injecting software simulated signals into the data and recovering them with one of the unmodeled GW detection pipelines, we address the robustness of the glitch-reduction technique to efficiently remove glitches with no unintended effects on GW signals.


Introduction
The first detection of a GW signal from a binary black hole (BBH) merger [1] on September 14th, 2019 was the breakthrough of GW astronomy. Since then, about 50 high confident detections [2][3][4][5][6][7] have been made by the network of GW detectors arXiv:2105.10522v2 [gr-qc] 8 Jul 2021 currently consisting of the two LIGO detectors in the United States [8], the Virgo detector in Italy [9], the GEO-HF detector in Germany [10], the KAGRA detector in Japan [11], and eventually the LIGO-India in India [12]. One of the dominating factors for the breakthough in 2015 was the improvement of the detector sensitivity by a factor of ∼ 10 at the detector's most sensitive frequency compared to the initial LIGO [8]. Reducing transient and periodic noise sources in instrumental and environmental origin may enhance the identification of GW signals by increasing the significance of candidate events of astrophysical signals and avoiding classifying them as sub-threshold triggers [2,13].
The first detection of a GW signal from a BBH merger [1] on September 14th, 2019 was the breakthrough of GW astronomy. One of the dominating factors for the breakthrough in 2015 was the improvement of the detector sensitivity by a factor of ∼ 10 at the detector's most sensitive frequency compared to the initial LIGO [8].
Reducing transient and periodic noise sources in instrumental and environmental origin may enhance the identification of GW signals by increasing the significance of candidate events of astrophysical signals and avoiding classifying them as subthreshold triggers [2,13].
The physical couplings due to the detector design are noise sources such as the fluctuating amplitude of laser light in the arm cavities, fluctuations of the photon arrival time at the output-port photodiode, thermal fluctuations of mirror coatings, and optic suspension [14]. Besides, during observation runs, environmental and instrumental noise sources including wind or ground motions as well as opticcontrolling systems limit the detector sensitivity [15]. Reducing the effect of these noise sources on the detector's output, strain channel, is crucial to improve the detectability of any GW signals and better understand physics in the universe.
The Advanced LIGO (aLIGO) detector use approximately twenty hundred thousand auxiliary channels, or sensors/control-subsystems that monitor different aspects of environmental and instrumental conditions inside and around the detectors in parallel to the strain channel in the time domain. These auxiliary channels can be potential witnesses of couplings of noise sources and can be used to subtract the noise in the strain channel.
Long-lasting noise sources with a duration longer than ∼ 4 seconds have linear and/or non-linear non-stationary coupling mechanisms ( [16] set the analysis window to be 8 seconds for the long-lasting noise). Among others, the main technique to subtract linearly coupled long-lasting noise sources is to calculate coherence between the witness channels and the strain channel. For example, the main source in ∼ 50−1000 Hz frequency band coupled linearly with the jitter of the pre-stabilization laser beam in angle and size was subtracted [17][18][19][20]. For non-linearly and nonstationary coupled noise sources in environmental and instrumental origin, techniques with machine learning algorithms show successful subtractions and improve the detector sensitivity [16,[21][22][23][24].
Other than long-lasting noise, transient noise artifacts, or glitches significantly impair the quality of GW detector data and reduce confidence in the significance of candidate GW events because glitches may resemble astrophysical signals. Glitch removal is crucial for GW signal searches. However, glitches are the product of short-lived linear and non-linear couplings among interrelated detector subsystems, generally making it difficult to find their coupling mechanisms. To mitigate the effect of glitches on GW searches, LIGO-Virgo collaboration vetoes time periods where glitches are present [3,25] or reduces the significance of GW candidates based on the probability of the glitch presence [26]. Software engines such as usepercentage veto (UPV) [27], hierarchical Veto (hVETO) [28], Pointy Poisson [29], PyChChoo [30], and iDQ [31,32] find the statistical correlation of glitches between the detector's output and numerous auxiliary channels to identify the presence of glitches while avoiding the mitigation of astrophysical signals.
The above mitigation techniques have disadvantages: reduction of the analyzable data which might contain GW signals or remaining the effect of the glitch contamination on the data. Therefore, it is desirable to subtract glitches for overcoming the above disadvantages. The functional forms to model the coupling mechanisms of glitches can require detailed knowledge of instrumental subsystems and a large number of parameters [33]. Also, there exist situations where functional forms can not be obtained because of unknown physical mechanisms about glitches. In addition to subtraction techniques for long-lasting noise sources noted above, other machine learning techniques have shown promising applications to GW astronomy [34]. For instance, [35] present a method using various techniques including the wavelet decomposition to classify transient excess-power waveforms injected into the simulated Gaussian data and real data during the LIGO's sixth science run. [36] show updated results of a glitch-classification software Gravity Spy [37] to the third observation run (O3) data with the help of citizen scientists. [31] show the comparison of various machine-learning algorithms to predict the presence of glitches based on auxiliary channels [31]. [38] introduce a method to estimate the data containing a compact binary coalescence (CBC) signal lost due to the presence of an overlapping glitch. These successes motivate us to develop a new machine learningbased method to subtract glitches using auxiliary channels with no dependency from any astrophysical-signal waveforms and no precise prior knowledge of all the system configurations, allowing the method to be easily adaptive to changes in the detector settings.
In this paper, we present a machine learning-based algorithm that subtracts glitches in the detector's output using auxiliary witness channels. Using two classes of glitches that are adversely affecting unmodeled GW detection searches, we characterize the performance of our subtraction technique. Adding simulated GW signals to the data before subtraction, we validate the algorithm not to manipulate or introduce biases to the resultant estimates of the astrophysical parameters.

Glitch subtraction pipeline
We introduce the analysis pipeline applied to the ground-based GW data for glitch subtraction. The pipeline processes the spectrograms magnitude of time series recorded in the strain channel, h(t), and a set of witness auxiliary channels for a class of glitches, a i (t). Those witness auxiliary channels do not causally follow from the detector's output where astrophysical signals are present. Choosing such channels allows the pipeline to only subtract glitches while preserving astrophysical signals. The algorithm uses a 2-dimensional Convolutional Neural Network (CNN) which uses a user-chosen set of witness channels as the input and then outputs the predicted spectrogram magnitude of glitches in h(t). The output of the CNN is then converted to time series using the Fast Griffin-Lim (FGL) transformation [39,40] and conditioned before being subtracted from h(t).

FORMALISM AND LOSS FUNCTION
The time series h(t) in the strain channel can be formulated as where h(t) is the astrophysical signal that may be present in the data, g(t) is a usertargeted glitch waveform that is coupled with witness channels, and n(t) represents the sum of untargeted glitch waveforms and the noise that are not wanted to be subtracted. We design the pipeline to produce an estimate of g(t) from a set of witness channels a i (t) and subtract it from s(t).
Because the amplitude of a glitch waveform varies rapidly and differently in a given frequency bin over its duration, it is more efficient to build the frequencydependent features before the data being fed into the neural network to learn the glitch-couplings. We create the discrete short-time Fourier transform (STFT) from the time series. The STFT divides the time series into small segments and calculates a discrete Fourier transform of each of divided segments. The STFT comprises the real and imaginary parts, which can be inverted to the corresponding time series. However, we find that the neural network more efficiently learns the wanted output with the magnitude of STFT (mSTFT) than using both real and imaginary parts together because the mSTFT has a simpler pattern than the real and imaginary parts individually have. The limitation of using the mSTFT is that it is not invertible. To compensate for this limitation, we use the FGL transformation [39,40] to invert the mSTFT to the corresponding time series by estimating the phase evolution. Because the phase is not fully recovered in this transformation, the phase of the estimated glitch waveform might slightly shift from the phase of the true waveform and the estimated waveform might have the opposite sign of the amplitude. Using the least square fitting method, we allow the amplitude of the estimated glitch waveform to change up to a factor of ±3 and allow the phase to shift up to ±0.02 seconds. We find that correcting the amplitude of the estimated waveform produces a better subtraction while preserving potential astrophysical signals. Also, changing the magnitude of the estimated waveform helps us avoid unwanted data manipulation because the amplitude correction factor tends to be zero when the estimated waveform mismatches the true glitch waveform.
Because of the complexity in estimating g(t) directly, we design the neural network to uses the witness channels and produce an estimate of the mSTFT G(t, f ) of a glitch waveform, where f denotes the frequency. The neural network can be represented as a function F(A i (t, f ); θ) which maps the magnitude mSTFTs A i (f, t) of the witness channels to G(t, f ) given a set of parameters θ. The parameters are obtained by minimizing a loss function J which denotes the difference between the real glitch mSTFT and the predicted counterpart. The operation in the network can be formulated as θ = argmin θ J G(t, f ), F(A i (t, f ); θ ) . ( In the analysis, we choose the loss function to be the mean squared error (MSE) across each pixel of mSTFT as where k denotes each pixel, N is the total number of pixels, andĜ[k; θ] := is an estimate of mSTFT obtained with the neural network. In practice, the dimension of mSTFT depends on the duration of g(t), its sampling rate, and a user-chosen frequency resolution for mSTFT. Using an appropriate combination of the neural internal window function called kernel and the mSTFT dimension, the neural network can efficiently minimize the loss function.

DATA PRE-PROCESSING
The target class of glitches and its witness channels are chosen to create the data set for the neural network. We select a target glitch class that is labeled by a machinelearning classification tool called Gravity Spy [37] which classifies glitches based on their time-frequency representation (TFR) morphology using the Q-transform [41].
Gravity Spy might misclassify some glitches by assigning inappropriate names. Also, some glitches are not loud enough to be worth to be subtracted. Because of the above two reasons, we select a set of glitches in the class with a relatively high signal-to-noise ratio (SNR) threshold (e.g., 10) and classification confidence level (e.g., 0.9) for the glitches.
To find the witness channels for these glitches, we use the software called PyChChoo which allows us to analyze safe auxiliary channels in the coincident windows for the glitches and discover witness channels statistically. Because of the computational efficiency for training the neural network and the achievement of an efficient prediction made by the network, using only witness channels without nonwitness channels is sufficient. We use the auxiliary channels having the probability that the glitch set is louder than the quiet set greater than 0.9 as the witness channels. Because different noise couplings might produce a similar TFR morphology, not all of the glitches in the class have a strong correlation with excess power recorded in the witness channels. Therefore, we select the subset of glitches in a class with their top-ranked witness channel has the probability of an excess-power measure belonging to the glitch set, above 0.9 [30] to make sure that the subset of glitches has a strong correlation with the witness channel. Around the time of glitches, we use a set of time series of the strain and the witness channels with a duration of 36 seconds. Glitches are excess power transients that are distinctively different from long-term varying noise sources. To filter out the long-term varying noise and magnify the characteristic of glitches, we whiten the time series by calculating the convolution between the time series and the time-domain finite-impulse-response filter created with the median value of amplitude spectrum densitys (ASDs), where each ASD is the square root of power spectral density (PSD) obtained by calculating the ratio of the square of the fast Fourier transform (FFT) amplitude with the Hanning window [42] of the divided time series with a duration of 2 seconds, to a given frequency-bin width (see details in [43]). The choice of 36-second duration for the time series is motivated such that the median value of ASD appropriately captures the long-term varying noise characteristic.
To have the same dimension for the mSTFTs of the glitch waveform and a i (t) and save the computational cost, we re-sample the whitened time series with the same sampling rate. For the glitch class used in the results in Sec. 3, we choose the sampling rate to be the lowest sampling rate of witness channels. This choice of sample rate has the Nyquist frequency ∼ 2 times higher than the highest frequency of the glitch so that all the characteristics of glitches are captured. To subtracting a glitch waveform from s(t), this choice also makes the resolution of the predicted glitch waveform small enough to apply a small time shift (∼ 0.02 seconds) for the phase correction.
To extract a glitch waveform g(t), we consider two distinct time-frequency regions. Glitches are expected to be present in one of the regions and not present in the other region. For example, for one of the glitch classes called Scattered light glitches in Sec.3, we consider the two frequency regions above or below the highest frequency of the glitch (e.g., 100 Hz). We assume that the upper-frequency region represents the STFTs of n(t) in Eq. (1) and the lower-frequency region represents the STFT containing glitches (see Appendix A for the verification of this assumption). We keep pixels in both real and imaginary parts in the lower frequency region of the STFT with their magnitude values above 99 percentile (see the study regarding this threshold choice in Appendix C) of the mSTFT value in the upper-frequency region, otherwise, set the pixel values to be zero. Subsequently, we set the pixel values in the upper region to be zero as well. After extracting the STFT of a glitch waveform, we invert the STFT to the time series to obtain the extracted glitch waveform. For Scattered light glitches, we find that the median value of the overlap O (defined in Eq. (4)) between the predicted mSTFTs from the network and the mSTFTs of the extracted glitch waveform with the cutoff frequency of 100 Hz is 1% greater than the overlap with the cutoff frequency of 200 Hz. As Appendix A shows that the mSTFT does not contain excess power above the Gaussian fluctuations in the frequency region above 100 Hz in the data of Scattered light glitches, choosing the cutoff frequency of 200 Hz lets the data keep a few pixels of Gaussian fluctuations above 200 Hz. Therefore, the overlap with the cutoff frequency of 100 Hz is larger than that with the cutoff frequency of 200 Hz. We choose a different choice in splitting the timefrequency region in STFTs for the other class of glitches (see details in Sec. 3). To determine choices of splitting the time-frequency region in STFTs, one can use the method shown in Appendix A and/or use the peak times and peak frequencies of Omicron triggers [41] to find if glitches are isolated in the time domain (see details in Appendix B).
To help the network learn the excess-power couplings more efficiently, we then divide each time series into smaller overlapping segments. Each training sample comprises segments from multiple witness channels. Larger overlap durations increase the data-set size so that the network can have more learnable resources and hence predict the output efficiently at the cost of computational time and memory.
Some of the segments might contain no glitches or only glitches that are not coupled with a set of chosen witness channels. Removing such segments from the data set helps the network learn the couplings more robustly. Using constant Qtransform (CQT) in nnAudio [44], we only select segments with the peak pixel in the CQT of the strain channel being loud enough and being near the peak pixels of at least one of the witness channels within a coincidence window (e.g., less than 1 second). During the above process, we further select the subset of segments with the peak frequency within an expected frequency range. As a criterion of the loudness of the pixel, we consider the peak pixel to be loud enough compared to be the Gaussian fluctuations when the peak pixel value is louder than the 90 percentile of the pixel values by a factor of 3.
After selecting the subset of segments, we finally create the mSTFTs of training sample segments and the corresponding mSTFT of a glitch waveform. To let the network learn efficiently, we normalize the mSTFTs of the strain and witness channels with their mean µ and standard deviation σ, then use these normalized mSTFTs as the input and the true output in Eq. (2) to build the network model parameters.

NEURAL NETWORK ARCHITECTURE
As mentioned in the earlier section, the pipeline uses a 2-dimensional CNN that uses the witness channels and predicts the mSTFT of a glitch waveform in the detector's output. The input for the network is a multi-dimensional image with the width (height, depth) corresponding to time-bins (frequency-bins, channels). The CNN typically consists of a series of convolutional layers, where each layer uses discrete window functions, or kernel with trainable weight. After taking the input image, the layer slides its kernel through the input image and then computes the dot products between the kernels and the portion of the image inside of the kernel. Typically, the kernel's dimension is smaller than that of the image so that the CNN learns local features in the image, making the network suitable to process locally outstanding characteristics of excess-power transients in the image. The output of each layer is passed to a non-linear activation function and becomes the input for the subsequent layer. Because the output of each layer is down-sampled, the subsequent layer represents the input image with a fewer number of features. This sequence of layers is known as an encoder that can extract the glitch-coupled excess-power characteristics and suppress the slowly-varying noise in the image [45]. After the convolutional layer, the network consists of transposed convolutional layers known as a decoder [46][47][48][49]. Each transposed convolution layer inserts pixels of zero between the pixels in the input image and then slides the kernel to computes the dot products of kernels and the modified input image within the kernel. As a consequence, transposed convolutional layers recast the encoded image to have a higher number of pixels by up-sampling. Comprising the encoding and decoding layers (so-called autoencoder ) [46][47][48][49], the output of the network will have the same dimension as the glitch image with extracted glitch-coupled features.
The above considerations motivate us to use the fully connected convolutional autoencoder in our network that provides an estimate of the glitch image from the witness channels' counterpart. In addition to the autoencoder, we employ a convolutional layer before the encoder to normalize the input images, and use a convolution layer after the decoder to make the output-image dimension to be the same as that of the glitch image. More specifically, the input images of the witness channels are first passed to the input convolutional layer and then normalized with Batch Normalization [50]. To make the dimension of the input and output images for the network, the input layer uses a stride of 1 and an appropriate zero padding arrangement. In the encoder, the width and height of the image are reduced by a factor of 2 while the depth (or the number of channels) is increased by a factor of 2. Instead of using pooling layers (e.g., Max Pooling [51]), each layer makes use of a stride of 2 with an appropriate zero padding arrangement. The output of the encoding layer is passed to the decoding layer. In the decoding layer, the width and height of the image are increased by a factor of 2 while the depth is reduced by a factor of 2 by using a stride of 2 with an appropriate zero padding arrangement. The output of the decoding layer is fed into the output convolutional layer, in which an appropriate kernel with a stride of 1 is chosen to make the final output image has the same dimension as that of the glitch image. Except for the output layer, the output of each layer is passed to an activation function before being fed into the subsequent layer.
We adopt the symmetrical structure for the autoencoder similar to [16] because each convolutional layer is commonly known to learn a different level of characteristics of the input images. An earlier (later) layer in the encoder tends to learn a lower (higher) level of characteristics. Hence, the first (last) layer in the encoder extracts a lower (higher) level of characteristics that are then recast by the last (first) layer in the decoder. Likewise, the intermediate levels of characteristics are also extracted and recast by a pair of layers in the encoder and decoder.
In our analysis, the network comprises four convolutional layers for both the encoder and decoder. We choose different sizes (∼ 10) of a kernel for different classes of glitches. For the activation function, we the ReLU [52,53] in the input and the autoencoder layers. We do not use any activation function after the output layer. Each encoding layer increases the number of channels from the value in the input image to 8, 16, 32, and 64, respectively. Each decoding layer decreases the number of channels in inverse order.

TRAINING AND VALIDATION
The analysis of the network can be divided to be two parts; training and validation. During training, the data set are divided into smaller chunks of data, or so-called mini-batches to reduce the computational memory. Data in each mini-batch are fed into the network and the loss function in Eq. (3) is computed by averaging over each mini-batch. The network parameters θ are updated according to the gradient of the loss function with respect to the parameters. For calculating the gradient, we use one of the first-order stochastic gradient descent methods, ADAM [54]. We iteratively update θ by repeating the above calculations over a number of cycles, or epochs. After each epoch, the loss value and the coefficient of determination are calculated using the validation set to prevent over-fitting; the network parameters are tuned with the training set so that the network tends to represent the glitch couplings contained in the training set instead of representing the couplings in a broader data set such as a validation set which is not used to turn the network parameters. To avoid over-fitting, the validation set is chosen not to be overlapping with any data in the training set. We stop the iteration if the network shows over-fitting according to the coefficient of determination of the validation set.

OUTPUT DATA POST-PROCESSING
The output of the network is conditioned before subtracting glitches in the detector's output data. Because the glitch mSTFT is normalized with the mean µ and the standard deviation σ of the pixel values across the training set before being fed into the network, we invert the normalization by multiplying the predicted mSTFT by σ and add µ. While the true mSTFT has no negative values by definition of mSTFT, the predicted mSTFT could have negative values due to the imperfection of the network prediction. We find that the negative values of the predicted mSTFT are typically anti-correlated with the values of the true mSTFT. The median values of the overlap O in 4 between negative pixels of the predicted mSTFTs and the corresponding pixels in the true mSTFTs are approximately −0.1 and −0.35 for Scattered light and Extremely loud glitches, respectively. Taking absolute values of the predicted mSTFTs makes these negative pixels have a positive effect on the overlap of all pixels though the effect is small because the average of the absolute values of the negative pixels is only ∼ 5% and ∼ 9% of the average value of the positive pixels for Scattered light and Extremely loud glitches, respectively. Using the FGL transformation [39,40], we estimate the glitch waveform from the absolute value of mSTFT predicted by the network.
Because the estimated waveform has a slightly larger or smaller amplitude compared to the extracted glitch waveform and the phase of the estimated waveform is slightly shifted due to the network-prediction imperfection and the phaseestimation uncertainty in the FGL transformation, we let the estimated waveform change only its amplitude by a factor up to ±3 and change the phase up to ±0.04 seconds to subtract glitches efficiently and avoid introducing unintended effects on astrophysical signals potentially present. To determine values of the amplitude and phase correction factors, we use the least square fitting method within the time periods where the glitch is present. To determine the portion of the glitch presence, we (1) calculate the absolute values of the estimated glitch waveform, (2) then smooth the curve with the convolution with the rectangular function, and (3) finally determine the time window where the absolute values are above a threshold.
Because the glitch waveforms (including both the estimate and extracted) are fluctuated instead of that their values are smoothly increased towards the peak of amplitude, taking the time window where the absolute values are above a threshold without smoothing the curves makes the time windows to be divided and to cover only sub-portions of glitches (corresponding to high thresholds) or the time window to cover wide portions including the region with no glitches (corresponding to low thresholds). Therefore, we employ the convolution with a rectangular function as one of the smoothing methods. We typically set this threshold to be the ∼ 90 percentile of the absolutes values of a set of the estimated glitch waveforms (see details for each glitch class in Appendix E.2 and Sec. 3). To subtract glitches even more efficiently, we divide this time portion and apply the least square fitting against the strain data within the divided portions (see details for each glitch class in Sec. 3). Smaller lengths of the divided portions allow us to subtract glitches more efficiently but less robustly preserve the waveform of astrophysical signals. To balance the above two factors, we divide the time portion finer around the center time of the glitch waveform because typically glitches have higher frequencies and larger amplitudes around the center time of the waveform. In the least square fitting, we find that the subtraction is better by separating the fitting with the bounds of the amplitude correction factor either in the range (0, 3) or (−3, 0) and then choose the better fitting result based on its coefficient of determination R 2 .

Pipeline performance on LIGO data
We apply our glitch subtraction pipeline to the data of the LIGO-Livingston (L1) detector from January 1st, 2020 to February 3rd, 2020. We choose two district classes of glitches comprising different types of noise couplings that are dominant events for creating background coherent WaveBurst (cWB) [55,56] triggers with high-ranking statistics to study and quantify the performance of our pipeline.
To quantify the network-prediction accuracy, we use the overlap of the true mSTFT and the predicted mSTFT of the glitch given as where G[i] andĜ[i] are the true and estimated mSTFTs of the glitch, respectively, and (i, j, k) run over pixel indices, and N is the total number of pixels. The overlap indicates the direct measure of the network prediction accuracy, ranging from 0 (mismatched) to 1 (perfect matched). SNR is the dominant factor for the detection of GW signal searches. Lowering SNR of noise artifacts improves the detectability of GW signals. . We quantify the performance of the glitch reduction with our pipeline by calculating the fraction of glitch SNR reduction (FNR) after the subtraction as where SNR b and SNR a are the matched-filter SNRs [57] obtained using the extracted glitch waveform as a template and the data before and after the subtraction, respectively (see how to obtain the extracted glitch waveform in Sec 2.2). Values of FNR close to 1 indicate efficient glitch reductions while negative values of FNR imply the increase of the glitch energy in the data.

SCATTERED LIGHT GLITCHES
We first apply our subtraction technique to a class of glitches called Scattered light glitches. Typically, winds and/or earthquakes shake the detector and move the test mass mirrors in the detector arms. The displacements of the mirrors in the longitudinal or rotational directions differ the relative arm lengths to produce archlike glitches in the mSTFT in 10 − 90 Hz, with a duration of ∼ 2 seconds [25,36,58].
To identify the witness channels for this glitch class, we use PyChChoo [30] on the list of Scattered light glitches with SNR above 10 between January 1st, 2020 and February 3rd, 2020 in the L1 detector from the Gravity Spy catalog. We use two witness channels with high confidence: L1:ASC-CSOFT P OUT DQ and L1:ASC-X TR B NSUM OUT DQ which monitor the common length of the two arms and the transmitted light from the mirror at the end of the x-arm, respectively.
To check that the top two witness channels are sufficient, we train the network with various sets of channels with the confidence p g of being witnesses of this glitch class from p g = 0.93 up to p g = 0.71 [30]. We consider 11 different channel sets, where i th set contains up to i th ranked channels. We train the network with a learning rate of 10 −3 (10 −4 /10 −5 ) for the 1-10 (11-30/31-60) epochs, where learning rates determine the gradient to update the network parameters and smaller learning rates correspond to smaller gradients. We terminate the training process if a value of R 2 in the validation set plateaus, i.e., a value of R 2 in the current epoch does not differ from the value in the previous epoch greater than ±0.001%. Figure 1 shows losses in Eq. (3) and overlaps in the validation sets, the validation sample size, and the GPU memories used to train the network for various channel sets. Overlaps at the termination for all channel sets range in ∼ 0.7 ∼ 0.8, where the GPU memory for the 11 th channel set is greater than the memory for the 1 st channel set by a factor of 2.5. Using a higher number of witness channels with high confidence provides a larger amount of glitch-coupling information to the network and let the performance better whereas adding low confident witness channels only provides nonglitch-coupling information to the network and does not improve the performance because the network seems not to use those low confident channels. Furthermore, using low confident channels might add data samples that have chance coincident excess power between these channels and the strain channels and/or subdominant glitch-coupling information (see details for the coincident selection in Sec. 2.2), where the size of the validation set for the 11 th channel set is greater than the size for the 1 st channel set by a factor of 1.5. Therefore, the termination overlap tends to decrease with the use of redundant GPU memories by adding channels above the 5 th rank as shown in the top-right panel in Fig. 1. Using the first two ranked channels is sufficient because the termination overlap for the 2 nd channel set is only less than 3% smaller than the largest value of the termination overlap (obtained by the 4 th channel set) and saves 18% GPU memory. In the following, we use the top-tow ranked channels for Scattered light glitches.
We pre-process the data as described in the previous section with the time series of the strain and the witness channels being re-sampled to a sampling rate of 512 Hz. During the pre-processing, we consider the frequency range above 100 Hz to be the background-noise region and use the 99 percentile pixel-value to extract the glitch waveform below 100 Hz. We also apply a high-pass filter to the strain channel at 10 Hz. We set the sample dimensions of the (training/validation/testing) 10 20 sets to be (9131/2233/678) with the segment overlaps of (93.7/93.7/75)%, where the segment overlap is the fraction of the time-window overlap between segments created by sliding a time window to divide a larger segment into smaller segments (so-called data augmentation [59][60][61]). We set that there is no overlapping time between the three sets and the testing set is later than the other sets. We create the mSTFT with a duration of 8 seconds, a frequency range up to 256 Hz, and (time/frequency) resolutions of 0.0625 seconds and 2 Hz, respectively. We use a square kernel with a size of (8,8) in the autoencoder in the network. During the post-processing, we consider the region where the glitches are present to be the time when the absolute value of the estimated glitch waveform is above the 75 percentile of the corresponding values across the testing set. We choose the divided window length for the least square fitting to be as small as 0.1 seconds and expand the window length from the center of the glitch by a factor of 1.1. Figure 2 shows the distribution of the overlap of the mSTFT and the FNR of the testing set of Scattered light glitches. More similar mSTFTs correspond to more similar waveforms after the FGL transformation, resulting in more efficient reductions of glitch SNRs. We find the overlap ranging from ∼0.6-0.9 and the FNR ranging ∼ 0.4 − 0.7 in 1-σ percentiles and no negative FNR, indicating no additional glitch energy added to the strain data.   by the network. However, only the fractions of this glitch are subtracted due to our selection criterion about the glitch presence mentioned above. In this case, the small value of FNR = 0.02 is also due to the presence of a non-Scattered light glitch at ∼ 0.5 seconds because this glitch contributes to SNR b dominantly. We note that we build the network for a particular class of glitches so that the glitch at ∼ 0.5 seconds in the least case is consistent with the performance of our built network.

EXTREMELY LOUD GLITCHES
Unlike Scattered light glitches where the waveforms in the strain can be analytically modeled with monitored mirror motions and the suspension systems [33], many other glitches are so far not modeled because of an incomplete understanding of their physical non-linear noise-coupling mechanisms. The non-linear activation function used in the network allows us to model non-linear noise couplings and subtract glitches. We apply our method to a class of Extremely loud glitches with SNR above 7.5 between January 1st, 2020 and February 3rd, 2020 in the L1 detector from the Gravity Spy catalog. We use 4 witness channels with high confidence: L1:LSC-POP A LF OUT DQ, L1:LSC-REFL A LF OUT DQ, L1:ASC-X TR A NSUM OUT DQ, and L1:ASC-Y TR B NSUM OUT DQ, identified by PyChChoo [30]. This glitch class is expected to be produced by the laser intensity dips and have extremely high excess power in 10-4096 Hz, lasting ∼ 0.2 seconds.
To check that the choice of witness channels noted above is sufficient, we train the network with various sets containing channels with p g = 0.99 up to p g = 0.55. We consider 11 different channel sets, where i th set contains up to i th ranked channels. Also, we consider the 12 th channel set containing the 1-3 th ranked channels and 5 th ranked channels because the 3 rd ranked channel (L1:ASC-Y TR B NSUM OUT DQ) with p g = 0.98 and the 4 th ranked channel (L1:ASC-Y TR A NSUM OUT DQ) with p g = 0.97 both record the transmitted light in yaw-direction in the alignment length control sub-system and have almost the same glitch-coupling information. The 5 th ranked channel (L1:LSC-POP A LF OUT DQ) records the transmitted light in low frequencies from the power recycling cavity. With the same procedure in Sec. 3.1, we compare losses, overlaps, GPU memories, and size of the validation sets for all channel sets as shown in Fig. 6. The termination overlaps range from 0.77 obtained with the 10 th channel set to 0.84 obtained with the 3 rd channel set. The overlap decreases by adding channels 4-11 th ranked channels. In particular, 6-11 channels have values of p g < 0.57, indicating no evidence of being witnesses so that the network obtain no significant glitch-coupling information from these low confidence channels with the use of redundant GPU memories up to 8.4 GB. In the following, we choose the 12 th channel set containing witness channels noted in the previous paragraph because its termination overlap is only less than 2% compared to the largest value obtained with the 3 rd channel set, and 0.5% increase of the data set.
During the data pre-processing, we re-sample the time series of the strain and the witness channels to a sampling rate of 2048 Hz and apply a high-pass filter at 10 Hz. We consider the time range outside of the 5-second window around the glitch time from the Gravity Spy catalog to be the background-noise region and use the 99 percentile pixel-value to extract the glitch waveform within the 5-second window because these glitches are isolated and not repeating, unlike Scattered light glitches. We set the sample dimensions of the (training/validation/testing) sets to be (3879/940/1233) with the segment overlaps of (96.8/96.8/87.5)%, where there is no overlapping time between the three sets and the testing set is later than the other sets. We create the mSTFT with a duration of 2 seconds, a frequency range up to 1024 Hz, and (time/frequency) resolutions of 0.0156 seconds and 8 Hz, respectively. We use a  rectangular kernel with a size of (13,4) in the autoencoder in the network. During the post-processing, we consider the region of the glitch presence to be the time when the absolute value of the estimated glitch waveform is above the 90 percentile of the corresponding values across the testing set. We choose the divided window length for the least square fitting to be as small as 0.02 seconds and expand the window length from the center of the glitch by a factor of 1.1. Figure 7 shows the distribution of the overlap of the mSTFT and the FNR of the testing set of Extremely loud glitches. We find the overlap ranging from ∼0.7-0.9 and the FNR ranging ∼ 0.1 − 0.6 with 1-σ percentiles and no negative FNR indicating no additional glitch added to the strain data.

INJECTION RECOVERY WITH COHERENT WAVEBURST
Subtracting glitches results in a new strain data which is expected to contain smaller energy due to the presence of glitches, leading to better detectability of astrophysical signals. One way of examining the robustness of our glitch-subtraction method is to add software-simulated signals with known astrophysical parameters into the strain   data before subtraction and use GW detection pipelines to recover the injected signals. In this process, we can assess whether the glitch subtraction technique reduces only the targeted glitches without manipulating the measured astrophysical signals. We use our glitch-subtraction method after injecting a signal in coincidence with a glitch. The presence of glitches adversely affects the unmodeled GW signal searches that do not rely on known waveforms in priori. In first half of the third observation run (O3a), the percentages of the single-detector observing time removed by the data-quality vetoes for the unmodeled searches are greater than the percentages for modeled searches by a factor of ∼ 2.7 and ∼ 7.9 for the LIGO-Hanford (H1)-detector and L1-detector, respectively [3]. Therefore, it is more beneficial to apply glitchsubtraction techniques for unmodeled searches. We use cWB to recover injections before and after subtraction and compare these recovered signals as well as the recovered signal injected in a simulated colored Gaussian noise with the PSD of the L1 data when the glitch subtraction is applied.
To account for the performance of our glitch subtraction operated on only the L1 data, we create a simulated colored Gaussian noise for the H1 data with the same sensitivity as the L1 data and inject signals both in the L1 and H1 detectors where the signal coincides with a glitch in the L1 data. The ranking statistic ρ of cWB accounts for the correlation of a signal injected on the two-detector data so that higher-ranking statistics for a given injection imply that the recovered signal in the L1 data is more similar to the signal in the H1 data, indicating successful glitch subtraction and a better detectability.

Gaussian-modulated Sinusoid Injections
Following studies of unmodeled GW signal searches [62][63][64][65][66][67], we inject a circularly polarized Gaussian-modulated sinusoid signal: where f c is the central frequency, t 0 is the center time, d is the distance to the source, A is an arbitrary amplitude scaling factor, Q determines the length of the signal. Motivated by studies [62][63][64][65][66][67] that use f c = 150 Hz, we also choose f c = 150 Hz. To validate the glitch subtraction technique successfully subtract glitches in the presence of signals with duration compatible glitches, we choose Q = 30, where the injected signal lasts ∼ 1.5 seconds which is compatible with the duration of Scattered light glitches. In addition to the above motivation, we consider signals similar to the first detection of intermediate binary black hole (IMBBH) [5], where the detected signal has a peak frequency of ∼ 50 Hz and ∼ 3 wave cycles. The improvement in the IMBBH detection by subtracting glitches might be useful to understand the mechanism of astrophysical populations [68]. Therefore, we choose f c = 50 Hz and Q = 5 for the second choice.
Using the two representative signal waveforms with different sets of parameters: f c = 50 Hz, Q = 5 and f c = 150 Hz, Q = 30, and choosing the injected SNRs to be uniformly sampled from a set of SNRs, the source direction to be isotropically sampled in the sky, the injected time to be uniformly sampled in a given time window, we examine the pipeline performance on the testing samples with an optimal set of values of FNR = 0.84, O = 0.92 (shown in Fig. 3) and a median value of FNR = 0.58, O = 0.65 (shown in Fig. 4) for Scattered light glitches as well as an optimal set of values of FNR = 0.87, O = 0.93 (shown in Fig. 8) and a median value of FNR = 0.33, O = 0.86 (shown in Fig. 9) for Extremely loud glitches. To assess the effect of the injection time on the pipeline performance, we consider two different injection-time windows: the subtracted portion in the testing-sample data or the full length of the testing-sample data. Because we apply the glitch subtraction in the partial data with excess power detected from the estimated glitch waveform and keep the original data for the rest of the data portion, we inject signals in the subtracted portion to study if our technique can subtract glitches that are overlapping with signals. We also consider the full length of the testingsample data as the injection-time window because subtracting glitches may affect detections of signals near to glitches but not overlapping with them. for Extremely-loud glitch set to assess the cWB detection performance for injections overlapping with glitches with high excess power. We inject 500 (250) waveforms with either high or low f c in the full (partial) injection-time window for each testing sample in each glitch class such that we have 16 injection-test sets. Figure 11 shows the enhancement of ρ after glitch subtraction for Gaussianmodulated sinusoidal injections. With the typical setting in cWB, only a signal with a ranking statistic greater than 6 is reported. We set the statistic for those missed signals to be 6 to quantify the enhancement due to subtraction. The percentages of injections with values of ρ after glitch subtraction greater or equal to the corresponding values before glitch subtraction ranges from 67% (obtained from the set with high-frequency signals injected in the full window of the optimal testing samples of Scattered light) to 100% (obtained from the set with high-frequency signals injected in the partial window of the optimal testing sample of Extremely loud glitches). Similarly, values the enhancement ρ a /ρ b , where ρ a and ρ b are ρ obtained from the data after and before the glitch subtraction, respectively, averaged over injections range from 1.2 (obtained from the set with high-frequency signals injected in the full window of the optimal testing samples of Scattered light) to 3.5 (obtained from the set with high-frequency signals injected in the partial window of the optimal testing sample of Extremely loud glitches).
Removing glitches with their characteristic frequencies close to that of signals typically improves values of ρ effectively because cWB reconstructs signals more effectively. Because Scattered light glitches have the largest power at a frequency of ∼ 30 Hz, the ρ enhancements for the low-frequency injection sets are larger than the enhancement for the high-frequency injection sets by a factor of up to ∼ 1.3. Similarly, Extremely loud glitches have a peak frequency of ∼ 110 Hz on average so that values of the enhancement for the high-frequency injection sets are greater than values for the low-frequency injection sets by a factor of up to ∼ 1.8.
Higher values FNR indicate larger reductions of excess power due to glitches. Hence, values of ρ obtained from the optimal testing samples are larger than values obtained from the median samples. The ρ enhancements in optimal-sample sets are greater than the median-sample sets by a factor ∼ 1.15 ∼ 1.23 (∼ 1.5 ∼ 1.7) and ∼ 1.8 ∼ 2.2 (∼ 1.7 ∼ 2.3) for the full (partial) injected window for Scattered light and Extremely loud glitches. Because Extremely loud glitches typically have extremely loud SNR ∼ 1500 while Scattered glitch glitches have SNR ∼ 17, subtracting Extremely loud glitches improves ρ more than subtracting Scattered light glitches. Also, the partial injection-window sets, where signals are overlapping with glitches tend to correspond to larger ρ enhancements than the full injection-window sets. Values of ρ are typically improved after glitch subtraction for the majority of injections near to glitches but not overlapping with them because the incoherent energy between detectors is reduced and the cWB obtains higher correlations of signals between detectors. Table. 1 shows values of the ρ enhancement and percentages of injection with non-reduced ρ after glitch subtraction for all sets.
Injections with reduced ρ after glitch subtraction are mainly due to 1) the least square fitting process operated between estimated glitch waveforms and the data, or 2) the cWB reconstruct process. The first reason is typically observed when the amplitude of signals is significantly large so that the least-square fitting method dominantly reduces the difference between a signal and an estimated glitch waveform in these cases. Hence, the signal energy is reduced. for example, these cases are observed when signals with high values of ρ obtained from the no-glitch data are injected at the center of glitches. Figure 12 shows an example failure case due to this reason: when an injected signal has larger or comparable to the amplitude of the overlapping glitch, the least square fitting method dominantly minimizes the difference between the estimated glitch waveform and the injection. The second reason is observed when the amplitude of the remaining glitches after subtraction is comparable to the amplitude of the nearby non-overlapping injected signals so that cWB reconstructs the sum of the remaining glitch and the true signal as a signal and the correlation of signals between detectors becomes smaller. The second reason can be seen in 0-1 seconds in the panels in the 3 rd -1,3 th columns in Fig. 11, where the original data without subtraction is used (see Fig. 9 and 8 for the subtracted portions). Figure 13 shows an example of unsuccessful cWB reconstruction for an injection nearby the remaining glitch after subtraction. In this case, the original data with injections is used around the injection time because no excess power is detected at the time of injections from the estimated glitch waveform. However, cWB reconstructs the injection differently before and after glitch subtraction. Before glitch subtraction, the cWB reconstruction process does not use the time portion containing the glitch because the amplitude of the glitch is not compatible with the signal amplitude. After glitch subtraction, the cWB used the data portion containing the remaining glitch whose amplitude is compatible with the amplitude of the injection. As a result, the cWB network correlation coefficient [55,56] is reduced to 0.65 from 0.99 after glitch subtraction, leading to ρ a /ρ b = 0.36.  More accurate signal reconstructions or higher values of the network correlation coefficient produce better estimates of the source direction. To assess the accuracy of the cWB source-direction estimates, we calculate the overlap of sky maps obtained with the no-glitch data and sky maps obtained with pre-subtracted data, where the sky map is provided as the probability distribution over pixelized solid angles and the sky-map overlap is calculated by taking the inner product between two sky maps similar to Eq. (4). Also, we calculate the overlap of sky maps obtained with the no-glitch data and sky maps obtained with the post-subtracted data. For injection missed by cWB with no reported sky maps, we set sky maps to be uniform probability distributions over solid angles according to the maximum entropy principle [69,70] Figure 13. Successful cWB reconstruction for a high-frequency (f c = 150 Hz and Q = 30) Gaussian modulated sinusoidal injection before glitch subtraction (left) and unsuccessful cWB reconstruction after glitch subtraction due to the nearby remaining glitch (right) in the optimal testing sample of Extremely loud glitches. The cWB correlation coefficients are 0.99 and 0.65 for the reconstructed signals before and after glitch subtraction, respectively. The ratio of ρ after and before subtraction is ρ a /ρ b = 0.36. Note that the time scales in the left and right panels are different due to the cWB automated reconstruction process, where the glitch is outside of the reconstruction time window before glitch subtraction.
for the least amount of knowledge about the source direction. We count percentages P sky of injections with the sky-map overlap of the no-glitch-post-subtraction data greater or equal to the sky-map overlap of the no-glitch-pre-subtracted data. Values of P sky greater 50% imply that estimates of the source direction become more accurate after glitch subtraction and P sky 50% indicates that source-direction estimates are compatible before and after glitch subtraction. Figure 14 shows ratios of sky-map overlaps between the no-glitch data and the post-subtracted data to sky-map overlaps with the former and the post-subtracted data for Gaussian modulated sinusoidal injections. Values of P sky range from 60% (obtained with the set with high-frequency injections in the full window of the optimal testing sample of Scattered light glitches) to 94% (obtained with the set with high-frequency injections in the partial window of the optimal testing sample of Extremely loud glitches). Because better signal reconstructions correspond to more accurate source-direction estimates, Values of P sky with optimal-testing-sample sets are greater than values with median-testing-sample sets by a factor of ∼ 0.86 ∼ 1.2 (∼ 1.01 ∼ 1.1) for Scattered light and Extremely loud glitches. The exceptional sets with high-frequency injections in the full window for Scattered light glitches have P sky = 60% for the optimal-testing-sample set and P sky = 69% for the mediantesting-sample set, respectively. However, they are compatible. 90% of injections in the above two exceptional sets have the ratio of the sky-map overlap of the no-glitchpost-subtracted data to the sky-map overlap of the no-glitch-pre-subtracted data in 0.93-2.3 (0.94-1.2) for the optimal (median)-testing sample set because the central frequency f c = 150 Hz is distinctively different from the peak frequency (∼ 30 Hz) of Scattered light glitches. We find that the maximum value of the ratio of the sky-map overlaps to be 150 and 4.5 for the above optimal and median-testing-sample sets, respectively. Values of P sky obtained with the low (high)-frequency sets are greater than values obtained with high (low)-frequency sets for Scattered light (Extremely loud) glitches by a factor of ∼ 1.08 ∼ 1.35 (∼ 1.01 ∼ 1.1) because removing glitches with their characteristic frequencies compatible with central frequencies of injections improve the cWB reconstructions more effectively. Table. 1 shows percentages of injection with the non-reduced ratio of sky-map overlaps after glitch subtraction for all sets.
To assess the accuracy of the cWB estimated central frequencyf c across all injections, we calculate the normalized residual: where f c is the injected central frequency. To quantify the similarity between two distributions, we calculate the two-sided Kolmogorov-Smirnov (KS) statistic S nb [71] between values of ∆f c obtained with the no-glitch data and the data before glitch subtraction as well as the KS statistic S na between values of ∆f c obtained with the former and the data after glitch subtraction. KS statistics are bounded between 0 and 1 and smaller values indicate two distributions are more similar. Values of the ratio R nb na := S nb /S na greater 1 imply that the cWB estimated values off c in the post-subtracted data are more similar to corresponding values in the no-glitch data than the pre-subtracted data while smaller values indicate the cWB estimates in the post-subtracted data are less accurate than the pre-subtracted data. R nb na ∼ 1 implies that the glitch-subtraction technique does not produce unintended effects on the data for the estimates of the central frequency. Figure 15 shows distributions of ∆f c obtained with the no-glitch, the pre-and post-subtracted data. Values of R nb na range from 0.41 (obtained from the set with lowfrequency injections in the full window of the optimal testing sample of Extremely  Fig. 15). Values of R nb na for the low (high)-frequency injection sets are greater than values for the high (low)-frequency injection sets by a factor of ∼ 1.04 ∼ 1.68 (∼ 0.95 ∼ 2.4) because subtracting glitches with their characteristic frequency compatible with signal frequency improves the cWB reconstruction accuracy. For high-frequency injection sets, values of R nb na obtained with optimal-testing-sample sets are greater than values obtained with corresponding median-testing-sample set by a factor ∼ 1.04 − 2.8 across the two glitch classes. For low-frequency injection sets, values of R nb na obtained with median-testing-sample sets are greater than values obtained with optimal-testing-sample sets by a factor of ∼ 1.05 ∼ 2.43 because of the contribution of high-frequency nearby remaining glitches to the cWB signal reconstruction, mentioned above for the set with the lowest value R nb na = 0.41. Table  2  Residual central frequency Counts Figure 15. Distributions of the residual central frequency ∆f c for high-frequency (1-2 th columns) and low-frequency (3-4 th columns) Gaussian-modulated sinusoidal waveforms injected in the full (1,3 th columns) and partial (2,4 th columns) windows of the optimal (1,3 th rows) and median (2,4 th rows) testing samples of Scattered light (1-2 th rows) and Extremely loud (2-4 th rows) glitches, respectively.

Binary Black Hole Injections
In addition to tests with Gaussian-modulated sinusoid signals in Sec. 3.3.1, we also assess the performance of the cWB-signal recovery by injecting non-spinning IMRphenomD BBH merger waveforms [72]. Following the choice of injection parameters used in [16], we choose the component masses to be uniformly distributed in [26,64] M with a constraint of the primarymass m 1 to the secondary-mass m 2 ratio in [0.125, 1], the source direction and binary orientation to be isotropically distributed, and the coalescence phase and the polarization angle to be uniformly distributed in [0, 2π] and [0, π], respectively. We choose the injected SNR sampled from a set of SNRs used in Sec. 3.3.1 and the injection time sampled in the full length of the testing sample data. We have 500 BBH injections for each set so that we have 4 BBH sets. Figure 16 shows the ρ enhancement for BBH injections. Percentages P ρ of injections with non-reduced ρ after glitch subtraction range from 76% (obtained with the median-testing-sample of Scattered light glitches) to 91% (obtained with the optimal-testing-sample of Scattered light glitches). Values of the enhancement ρ a /ρ b averaged over injections range from 1.2 (obtained with the median testing sample of Scattered light glitches) to 2.7 (obtained with the optimal testing sample of Extremely loud glitches). Subtracting significant energy due to glitches improves the cWB reconstruction so that values of the enhancement in Extremely-loud sets are greater than values in Scattered-light sets by a factor of 1.25 and 1.8 for the optimal and median testing sample, respectively. In Scattered-light (Extremely-loud) sets, the value of the enhancement for the optimal test set is greater than the value of the median testing set by a factor of 1.25 (1.8). Table. 1 shows values the ρ enhancement and P ρ for BBH injection sets. Similar to P sky for Gaussian modulated sinusoidal injections, we compute P sky for BBH injection set. Figure 17 shows ratios of sky-map overlaps between the noglitch data and the post-subtracted data to sky-map overlaps with the former and the post-subtracted data for BBH injections. Values of P sky range from 73% (obtained with the median testing sample of Scattered light glitches) to 80% (obtained with the optimal testing sample of Extremely loud glitches). Values of P sky in the optimal testing set are greater than values in the median testing sample set for Scattered light and Extremely-loud glitches by a factor of 1.06 and 1.14, respectively. Table. 1 shows values of P sky for BBH injection sets.
To assess the accuracy of the cWB estimated chirp mass across all injections, we calculate the normalized residual: whereM and M = (m 1 m 2 ) 3/5 /(m 1 + m 2 ) 1/5 are the cWB estimated and injected chirp mass, respectively. Similar to the procedure in the previous section, we calculate the two-sided KS statistic S nb [71] between ∆M obtained with the noglitch data and the pre-subtracted data as well as the KS statistic S na obtained with the former and the post-subtracted data. We calculate values of the ratio R nb na := S nb /S na . Figure 18 shows distributions of ∆M obtained with the no-glitch, the pre-and Extremely loud glitches, respectively. The first and second texts in labels in the legend denote that injections in the pre-subtracted and post-subtracted data are found or missed, respectively.
post-subtracted data in the optimal testing sample sets. The distributions in the optimal testing sample sets are comparable with distributions in the median testing sample sets. As shown in Table 2, values of R nb na range from 0.96 (obtained from the median testing sample of Extremely loud glitches) to 3.26 (obtained from the median testing sample of Scattered light glitches) as shown in Table 2. Values of R nb na are close to or greater than 1, indicating the glitch subtraction technique produces no unintended effect on cWB estimates for the chirp mass or improves the estimates. As shown in Fig. 19, S na and S nb are compatible for BBH sets so that the distribution of ∆M in the post and pre-subtracted data are similar. Values of S na and S nb for BBH sets are typically smaller than values for Gaussian modulated sinusoidal sets because BBH waveforms distinctively differ from glitch waveforms and cWB reconstructs BBH injections more effectively than the Gaussian modulated sinusoidal injections.

False Alarm Rate
The confidence of a GW signal candidate is quantified by the false alarm rate (FAR), or the rate of terrestrial noise events with their Table 1. Enhancements of the cWB ranking statistic ρ a /ρ b , i.e., ratios of the cWB ranking statistic ρ a after and ρ b before glitch subtraction averaged over injections, as well as percentages P ρ of injections with non-reduced ranking statistics and percentages P sky of non-reduced sky-map overlaps after glitch subtraction. The sky-map overlaps are calculated between sky maps obtained with the simulated colored Gaussian data and either of the data before or after glitch subtraction. Values of ρ are set to be 6 for injections missed. The sky-localization estimate is set to be uniformly distributed in solid angles.
Glitch class Testing sample Injection Full window Partial window ρ a /ρ b P ρ P sky ρ a /ρ b P ρ P sky  KS statistic of no-glitch and pre-subtracted data Figure 19. KS statistic obtained with the no-glitch data and the data before glitch subtraction VS. KS statistic obtained with the former and the data after glitch subtraction for residual peak frequency or chirp mass estimated by cWB in the cases of high-frequency (circle marker), low-frequency (cross marker) Gaussianmodulated sinusoid waveforms and BBH waveforms (star marker) injected in the full (slashed hatch) and partial window (dotted hatch) in the optimal (cyan face color) and median (white face color) testing sample of Scattered light (black edge color) and Extremely loud (red edge color) glitches, respectively.

Scattered light
used in the previous section, we use isolated samples in the testing sets with sample sizes of 237 and 156 for Scattered light and Extremely loud glitches from January 7 th , 2020 3:00 UTC to February 2 rd , 2020 03:55 UTC and January 15 th , 2020 17:44 UTC to February 3 rd , 2020 23:55 UTC, respectively. Because the testing sets used in the previous section have segment overlaps to account for statistical errors with larger sample sizes, we use the isolated testing samples from the sets to avoid double subtraction. We use the data containing the above testing samples during the detectors are observing, corresponding to the 20.4-and 20.5-day data from the L1 and H1 detectors, respectively. The percentages of the total duration of the isolated testing samples are 0.07% and 0.026% of 20.4 days for Scattered light and Extremely Table 2. Percentages (P n /P b /P a) of found injections in the (no-glitch data/data before subtraction/data after subtraction), as well as ratios R nb na := S nb /S na of the KS statistic S nb obtained with the no-glitch data and the data before glitch subtraction to the statistic S na obtained with the former and the data after glitch subtraction for residual peak frequency and chirp mass.

Glitch class
Testing sample Injection Full window Partial window Using the L1 data before glitch subtraction with the original H1 data and applying time shifts to the L1 data, we get the background trigger set, where time shits are applied to get triggers representing the noise events coincident between detectors by chance and enlarge the analysis time. Similarly, we also use the L1 data after glitch subtraction with the original H1 data to get another background trigger set. With time shift applied to the L1 data, we obtain 21.2-year equivalent background triggers both before and after glitch-subtracted data. Both trigger sets have the maximum values of ρ = 53.8 and the lowest FAR of 1.5 × 10 −9 Hz (corresponding to iFAR of 21.2 years). Figure 20 shows the FAR of background triggers before and after glitch subtraction. We find that the FAR is typically reduced in the interval from ρ ∼ 7 to ρ ∼ 12. The reduced FAR in this interval can be explained by the reduction of ρ in the subtracted part of the data. Figure 21 shows background triggers within the interval of the subtracted data portions. The distribution of these triggers is due to the quality of the L1 data. The average values of ρ are reduced by 13.2% and 1.9% for Scattered light and Extremely loud glitches. Because Scattered light glitches are typically subtracted with our technique more effectively than Extremely loud glitches, the former has higher percentages than the latter. The FAR is increased in the interval from ρ ∼ 12 to ρ ∼ 23 in the glitch subtracted data because of two triggers with ρ = 18.2 and ρ = 17.5. However, the L1 times of these two triggers are not within the subtracted data portions so that it seems to be due to a realization of cWB trigger-generation process.
Using these two background sets, we first evaluate injections that are not nearby and overlap with glitches. Because we have created the no-glitch data sets in the previous section with the colored Gaussian noise using the PSD of the real L1 data at the time of injections, we can consider injections in the no-glitch data to be those not nearby and overlap with glitches. As a lower limit, we set iFARs of injections with ρ greater than 53.8 to be iFAR = 21.2 years, which is the maximum length of our background.
Across high-and low frequency, and BBH injections, we find that ∼ 90% of injections in the no-glitch data have non-reduced iFARs after glitch subtraction. The iFARs after glitch subtraction is greater than the iFARs before glitch subtraction by a factor of ∼ 1.02 on average over injections because of the reduction of ρ in the background.
As mentioned above, the glitch subtraction technique may increase the ρ of injections near to or overlapping with glitches, causing higher iFARs. Figures 22 and  23 show iFAR distributions before and after glitch subtraction using corresponding backgrounds for Gaussian modulated sinusoidal and BBH injections, respectively. Percentages P i ab of injections with non-reduced iFAR after glitch subtraction range from 88% (obtained from the set with the high-frequency injections in the full window of the median testing sample of Scattered light glitches) to 100% (obtained from the set with high-frequency injections in the partial window of the optimal testing sample of Extremely loud glitches). The sets with the lowest and highest values of P i ab respectively correspond to the sets with the lowest and highest values P ρ because the increases in ρ for injections correspond to the increases in iFAR. The ratio R ab of iFAR values after glitch subtraction to iFAR values before glitch subtraction averaged over injections range from 1.03 (obtained from the set with high-frequency injections in the full window of the optimal testing sample of Scattered light glitches) to 1400 (obtained with high-frequency injections in the partial window of the optimal testing sample of Extremely loud glitches). Because high increases in ρ of injections correspond to higher increases in iFAR, values of R ab for sets with the optimal testing sample typically are greater than values for sets with the median testing sample by a factor of ∼ 0.7 ∼ 150 and ∼ 2.1 ∼ 4.3 for Scattered light and Extremely loud glitches. The sets with high-frequency injections in the full window for Scattered light glitches corresponding to the lowest factor of ∼ 0.7 have comparable values of R ab = 1.3 and R ab = 1.9 for the optimal-and median-testing-sample sets. Subtracting glitches with their peak frequencies close to the characteristic frequencies of injections lets cWB obtain higher values of ρ. Therefore, values of R ab for sets with (low/high)-frequency injections are greater than values of R ab for sets with (high/low) frequency injections by a factor of (∼ 3.1 ∼ 74 /∼ 1.2 ∼ 2.5) for (Scattered light/Extremely loud) glitches.
Weak signals (so-called sub-threshold triggers) near or overlapping with glitches that are missed by cWB or are not confident enough to be classified as astrophysical signals can gain sufficient confidence after glitch subtraction. If we assume an iFAR threshold for weak signals to be a month, percentages P i w of injections with iFAR above a month after glitch subtraction out of injections with iFAR below a month before glitch subtraction range from 1% (obtained from the set with high-frequency injections in the full window of the optimal testing sample of Scattered light glitches) to 57% (obtained with the set with low-frequency injections in the partial window of the optimal testing sample of Scattered light glitches). For Scattered light glitches, sets with low-frequency injections have values of P i w ∼ 40 ∼ 57% and sets with highfrequency and BBH injections have values of P i w ∼ 1 ∼ 4%, where values for the optimal and median-testing-sample sets are compatible. For Extremely loud glitches, sets with the optima-testing sample have values of P i w ∼ 20 − 44% and sets with the median-testing-sample have values of P i w ∼ 5 − 18%, where values for high-frequency sets are greater than values for the low-frequency sets by a factor of ∼ 0.7 ∼ 1.8. Table 3 shows values of (P i ab /P i n /P i w ) and ( R ab / R n ).  Figure 21. Background triggers within the intervals of the subtracted data portions for the data before (gray) and after (red) glitch subtraction for Scattered light (cross) and Extremely loud (circle) glitches, respectively. The dashed and solid lines denote the average values of ρ for these triggers for Scattered light and Extremely loud glitches, respectively. iFAR in the pre-subtracted data [year] iFAR in the post-subtracted data [year] iFAR of injeciotns in the clean time improved by glitch subtraction Figure 22. Distributions of the iFAR for high-frequency (1-2 th columns) and low-frequency (3-4 th columns) Gaussian-modulated sinusoidal waveforms injected in the full (1,3 th columns) and partial (2,4 th columns) windows of the optimal (1,3 th rows) and median (2,4 th rows) testing samples of Scattered light (1-2 th rows) and Extremely loud (2-4 th rows) glitches, respectively. The color scale denotes the ratio of the iFAR of injections in the no-glitch data evaluated with the background ρ distribution obtained with the post-subtracted data to that obtained with the pre-subtracted data. The shaded area denotes the iFAR less than 1 month.

Conclusion
In this paper, we have presented a new machine learning-based algorithm to subtract glitches using a set of auxiliary channels. Glitches are the product of short-live linear and non-linear couplings due to interrelated sub-systems in the detector including the optic alignment systems and mitigation systems of ground motions. Because of the characteristic of glitches, modeling coupling mechanisms is typically challenging. Without prior knowledge of the physical coupling mechanisms, our algorithm takes the data from the sensors monitoring the instrumental and environmental noise transients and then estimates the glitch waveform in the detector's output, providing the glitch-subtracted data stream. Subtracting glitches improves the quality of the iFAR in the pre-subtracted data [year] iFAR in the post-subtracted data [year] iFAR of injeciotns in the clean time improved by glitch subtraction Figure 23. Distributions of iFAR of BBH injections in the full window of the optimal (1/3 th columns) and median (2/4 th columns) testing sample of (Scattered light/Extremely loud) glitches, respectively. The color scale denotes the ratio of the iFAR of injections in the no-glitch data evaluated with the background ρ distribution obtained with the post-subtracted data to that obtained with the pre-subtracted data. The shaded area denotes the iFAR less than 1 month. data and will enhance the detectability of astrophysical GW signals.
Using two classes of glitches with distinct noise couplings in the aLIGO data, we find that our algorithm successfully reduces the SNR of the data due to the presence of glitches by 10 − 70%. Subtracting glitches from the data enhances the cWB ranking statistic by a factor of ∼ 1.03 ∼ 3.5 and ∼ 1.2 ∼ 2.7 averaged over Gaussian modulated sinusoidal injections and BBH injections, respectively. We find that the source-direction, central frequency and chirp mass estimated by cWB after glitch subtraction are comparable or more accurate than that before glitch subtraction. The iFAR of injections in the data portion in the absence of glitches is increased by ∼ 1.02 by subtracting glitches in ∼ 0.1% of the 20.4-day data from the L1 detector. We find that injections near to or overlapping with glitches typically have significant enhancements with glitch subtraction. The iFAR of those injections is increased by a factor ∼ 1.3 ∼ 1400.
In this paper, we focus on the two classes of glitches and apply the glitch subtraction technique to only ∼ 0.1% of the L1 data so that we find no significant reduction of ρ in the background. Creating the CNN network models for other glitch classes and subtract a higher number of glitches both in the L1 and H1 could provide the statistically robust measure of the effect of the glitch subtraction technique on the data.
Currently, the LIGO-Virgo collaboration vetoes glitch classes focused on this Table 3. Percentages (P i ab /P i n ) of injections (near to or overlapping with glitches/in the absence of glitches) with non-reduced iFAR after glitch subtraction and percentages P i w of injections with iFAR above a month after glitch subtraction out of injections with iFAR below a month before glitch subtraction, and ratios ( R ab / R n ) of iFAR after glitch subtraction to iFAR before glitch subtraction averaged over injections (near to or overlapping with glitches/in the absence of glitches).

Glitch class
Testing sample Injection Full window Partial window paper and other glitch classes with witness channels. For example, over the course of the 20.4-day data from the L1 data, ∼ 15000 glitches with SNR above 7.5 in the two classes are present and have a total period of ∼ 1.8% (so-called deadtime), which would be vetoed. By accounting for the deadtime and injections removed by the veto, the comparison of the volume-time integrals [20] between the vetoing method and the glitch subtraction technique allows us to find a better approach. We find that using the spectrograms of the data as the input for the network is more successful than using time series as the input. However, it might improve the glitch subtraction efficiency by using the FGL transformation as well as the amplitude and phase corrections within the loss function to train the network. Improved glitch subtraction would allow us to detect astrophysical signals with higher confidence and brings us a better understanding of the physics in the universe.

Appendix Appendix A. Comparison between Scattered light glitches and quiet times
To show that the frequency region above 100 Hz in time periods containing Scattered light glitches in the strain channel has no excess power and are compatible with the corresponding frequency region of the Gaussian noise, we compare 693 spectrograms containing Scattered light glitches with 306 spectrograms when the strain channel is quiet, statistically evaluate them using the KS test [71].
To create the data set of Scattered light glitches, we whiten the time series of the strain channel with a software called GWpy [43] and then apply a low-pass filter at 512 Hz as used in Sec. 2. We have the Scattered-light set with a sample size of 693 by selecting time periods with a duration of 8 seconds that contains Scattered light glitches. To have a set of quiet data, we use the observing-mode strain channel data with a duration of 4096 seconds beginning from April 2nd, 2019 at 5:04 UTC, without data quality issues such as the corrupting data, the presence of glitches, and hardware injections of simulated signals. We whiten and apply the high-pass filter to the quiet time series and then cut the edge of the whitened time series to remove artifacts of the Fourier transform. By diving the whitened time series into 8-second segments, we have the quiet-data set with a sample size of 301. We create mSTFTs of the Scattered-light set and the quiet set. Figure A1 shows the mSTFT of a Scattered light glitch and a quiet time. Figure A2 shows distributions of the mSTFTs in the frequency region above or below 100 Hz in the Scattered-light set and the quiet set. The Scattered-light (quiet) set has 1.6% and 9.5% (1.6% and 2.1%) of pixels with values above 10 for the frequency region above and below 100 Hz across the set, respectively. To verify the upper-frequency region of the Scattered light set has no excess power above Gaussian fluctuations, we calculate one-sided KS-test statistics for randomly selected 500 pairs of a mSTFT from the Scattered-light set and a mSTFT from the quiet set by taking the mSTFT variations in both sets into account. As the null hypothesis in the one-sided KS test, we consider mSTFT-pixel values of a Scattered-light-glitch data is lower than that of the quiet data because we want to verify if the hypothesis that the upper-frequency region of the Scattered light set has no excess power above Gaussian fluctuations can be rejected. Because the KS-test statistics calculated above contain the variation of mSTFT in both sets, we also calculate one-sided KS-test statistics for randomly selected 500 pairs of two mSTFTs from the quiet set. The right panel in Fig. A3 shows distributions of KS-statistics of these pairs. We then perform a one-sided KS test for the above two distributions of KS statistics. We find that the p-value of the test to be 0.099, which is not confident enough to reject the hypothesis that the upper-frequency region of the data containing Scattered light glitches has no excess power.
As a supplementary study, we perform the same procedure for the frequency region below 100 Hz. The left panel in Fig. A3 shows distributions of KS statistics calculated with mSTFT pairs in the frequency region below 100 Hz. We find the p-value of the test to be 2.4×10 −141 so that Scattered light glitches have excess power below 100 Hz.
To robustly verify that the null hypothesis that the mSTFT of the data containing Scattered light glitches has no excess power in the frequency region above 100 Hz, can not be rejected, we vary the number of pairs and random pair-selection  Figure A3. Distributions of KS statistics calculated with 500 random pairs of mSTFTs in the frequency region below (left) and above (right) 100 Hz. QQ and SQ denote pairs of two mSTFTs from the quiet set and pairs of a mSTFT from the Scattered-light set and a mSTFT from the quiet set, respectively.
realizations for p-values. Similar to the process mentioned above, we randomly select N pairs of a mSTFT from the Scattered-light set and a mSTFT from the quiet set as well as N pairs of two mSTFTs from the quiet set, subsequently calculate one-sided KS statistics for paired mSTFTs in the frequency region above x Hz for both Scattered-light-quiet and quiet-quiet pair sets. Then, we perform the one-sided KS test for the two KS-statistic distributions from two pair-sets. We vary values of N between 100 and 2000 and x between 30 Hz and 210 Hz to show the variation of p-values due to selected numbers of pairs N and cutoff frequencies x. Also, to see the effect of the pair-selection realizations on p-values, we repeat 10 times each test for a given value of N and x. Figure A4 shows the variations of p-values and the corresponding KS statistics. The p-values with cutoff frequencies below 75 Hz significantly smaller, indicating high confidence in rejecting the null hypothesis. For cutoff frequencies below 75 Hz, larger values of N corresponds to smaller p-values because the sample errors are smaller for a set with larger samples even if statistics with larger N are slightly smaller than values with smaller N . The p-values with a cutoff frequency above 100 Hz are comparable, meaning that no excess power is observed in the frequency region above 100 Hz irrespective of values of N and the pair-selection realizations.

Appendix B. Peak time of glitches
To determine choices of splitting time-frequency regions in mSTFTs of the data containing glitches, one can use the peak time of glitches to identify if glitches are isolated in the time domain. Figure B1 shows peak times of Scattered light glitches and Extremely loud glitches in 36-seconds time periods. Extremely loud glitches are isolated in the time domain and their peak times are generally within ±1 around the Omicron-trigger [41] times. Scattered light glitches are repeatedly present because this glitch class is generated due to swinging mirror motions caused by seismic activities. Therefore, to estimate the mSTFT of the Gaussian noise, we use the frequency region above 100 Hz for Scattered light glitches and the time region outside of ±2.5 seconds around the trigger time for Extremely loud glitches.

Appendix C. Pixel threshold to extract glitch waveforms
After determining the time-frequency region of the data to estimate the STFT of the noise, we extract a glitch waveform from the STFT. Because glitches generally can not be modeled, we choose a threshold for pixel values of the STFT to extract glitch waveforms. We keep pixels of the STFT (hereafter called excess pixels) with their magnitude values above a threshold estimated from the STFT representing the noise. The left panel of Fig. C1   simulated Gaussian data. Its right panel shows the histograms of mSTFTs of the injected data, the noise only, and the glitch only. Smaller thresholds let excess pixels have a larger number of noise and glitch pixels while larger thresholds let excess pixels have a smaller number of the noise pixels but smaller glitch pixels.
To determine the best pixel threshold, we take mSTFT pixels with their magnitude values below a threshold (hereafter called un-excess pixels) and then use the two-sided KS test between the un-excess pixels and mSTFT pixels of quiet times. If un-excess pixels are similar to quiet pixels, excess pixels tend to have the majority of glitch pixels and smaller numbers of noise pixels. To quantitatively determine the best threshold, we randomly select 200 pairs of un-excess-pixel mSTFTs with quiet mSTFTs and calculate two-sided KS statistics for each pair. We use the twosided KS statistic because we want to see the similarity of two mSTFT in a pair. Likewise, we randomly select 200 pairs of two quiet mSTFTs and calculate two-sided KS statistics. We take the ratio of the KS statistic averaged over the un-excess-quiet pairs to the KS statistic averaged over the quiet-quiet pairs. Figure C2 shows histograms of KS statistics of 200 quiet-quiet pairs, 200 un-excess-quiet pairs, and 200 excess-quiet pairs with different pixel thresholds. The ratio is closed to one when the un-excess pixels are similar to the quiet pixels, i.e., excess pixels contain the majority of glitch pixels and fewer noise-pixel. Whereas, larger values of the ratio imply that 1) only smaller amplitude noise pixels are contained in un-excess pixels, i.e., excess pixels contain a larger amount of noise pixels (corresponding smaller pixel thresholds), or 2) un-excess pixels have a higher number of glitch pixels, i.e., excess pixels have only a few glitch pixels (corresponding to larger pixel thresholds). We vary the pixel threshold from 50, 90, 99, 99.9, 100 percentiles of the mSTFT in the time-frequency region that is expected to contain no glitches (see the above sections). Also, we consider values of 100-percentile multiplied by 2, 3 and, 5 as the threshold. Therefore, minimizing the ratio over pixel thresholds allows us to find the best pixel threshold. For Scattered light glitches, the ratio reaches the lowest values ∼ 2.1 ∼ 2.2 with thresholds of 99, 99.9, and 100 percentiles. For Extremely loud glitches, the ratio has the minimum value of ∼ 3.9 with 99 percentile. Extremely loud glitches have extremely high SNR so that the ratio is more sensitive to the pixel threshold. According to the above study, we set the best threshold to be 99 percentile for both classes of glitches.

100
200 300 400 500 Pixel threshold 10 1 Ratio of two-sided KS statistics Scattered light Extremely loud Figure C3. The ratio of the unexcess-quiet pair KS statistic to the quiet-quiet pair KS statistic. Pixels thresholds below 100 denote percentiles of mSTFT in the time-frequency region that is expected to contain no glitches. Pixel thresholds above 100 denote 100 percentile multiplied by 2, 3, or 5.

Appendix D. Fast Griffin-Lim transform
Appendix E. Iteration in Fast Griffin-Lim Transformation Figure E1 shows the FNR as a function of FGL iteration with α = 0.99 (recommended value in [40]). Values FNR are comparable after 20 iterations. We choose 32 iterations (default value in [44]) and α = 0.99 in this paper. The data used in this paper does not guarantee the theoretical convergence inherited from the algorithm of the FGL transformation and the convergent point to be the global minimum inherited from the Griffin-Lim (GL) transformation. We address the accuracy of the FGL transformation on the data used in this paper in the next section. where (g i |g) is the noise-weighted inner product [74,75] defined as where S n (|f |) is the one-sided noise PSD,ã andb are the Fourier transform of the time series a(t) and b(t), respectively. The FGL transformed time series has a phase error. Also, we apply the phase correction before subtracting glitches from the data before subtracting glitches. we also calculate the match maximized over phase and time, defined as where the complex inner product (ã,be 2πif t 0 ) complex is defined as Values of M range from -1 (fully anti-correlated) to 1 (perfect match). We use extracted glitch waveforms in the testing set of Scattered light and Extremely loud glitches with sample sizes of 678 and 1233, respectively. Figure E2 shows distributions of the match and maximized match.
Samples with values of M ∼ ±1 and of M max ∼ 1 indicate that the FGL transformed time series are similar to the original extracted-glitch waveforms with some degree of phase shifts. Samples with values M ∼ 0 and M max ∼ 1 indicate that the FGL transformed time series have phase shifts such that they mismatch with the extracted glitch waveforms. Samples with values of M ∼ 0 and M max ∼ 0 indicate that the FGL transformed time series mismatch with the original extractedglitch waveforms with both phases and amplitudes. The median and 90 percentile of absolute values of the time shift t 0 which maximizes the match are ∼ 0.02 (∼ 0.003) and ∼ 0.06 (∼ 0.0034) seconds for Scattered light (Extremely loud) glitches, respectively. Figure E3 shows the optimal and least accurate FGL transformed time series. The optimal and least values of M max are ∼ 0.27 (∼ 0.27) and ∼ 0.99 (∼ 0.97) for Scattered light (Extremely loud) glitches, respectively. The FGL transform seems to produce no significant deviations on the amplitude in the portion where the extracted glitch waveforms have amplitudes close to zero. Therefore, the amount of the mismatch seems to be due to the amplitude difference in the portion where extracted glitch waveforms have large amplitudes.
To understand the meaning of values of M max between the FGL transformed time series and the extracted glitch waveforms in terms of amplitude uncertainty of the FGL transform, we calculate M max between the extracted waveform and themselves injected into the Gaussian noise with zero mean and a given standard deviation. Assuming the FGL transformation provides the amplitude error following a Gaussian distribution with zero mean throughout time series, we can find the standard deviation of the Gaussian noise corresponding to the amplitude uncertainty of the FGL transform. We find that Scattered light and Extremely loud glitches have the amplitude uncertainties equivalent to the ∼ 0.4σ and ∼ 8σ noise, respectively, as shown in Fig. E4. As discussed above, the FGL transformation produces no significant amplitude errors in the portion where the extracted glitch waveforms have amplitudes close to zero. Therefore, the above estimates using the noise with a given σ has a bias for waveforms with larger amplitudes in limited time portions such as Extremely loud glitches. One of the possible approaches to correct the phase and amplitudes of the FGL transformed time series is to split the time series into small segments and change the phase and amplitude in each segment (see details in Sec. 2.5).

Appendix E.2. Time window to subtract glitches
Applying the FGL transform to the estimated mSTFT from the trained network, we obtain the estimated glitch waveforms. Using the estimated glitch waveforms, we subtract glitches from the data. In the glitch-subtraction process, we determine portions of the data containing g glitches based on the estimated glitch waveforms. We only subtract the data in portions containing glitches and use the original data in the rest portions without the glitch subtraction. To determine the data portions, we use the estimated glitch waveforms in the testing set with sample sizes of 678 and 1233 for Scattered light and Extremely loud glitches. We first calculate the absolute values of the estimated glitch waveforms and then smooth the curve. We consider data portions to be the regions where the smoothed curves are above thresholds. Smaller thresholds make the data portions to be larger so that larger data are used in the subtraction process, where larger fractions of the data portions have no glitches and no need to be subtracted. Larger thresholds make the data portions to be smaller so that only small fractions of glitches are subtracted. We consider various percentiles of the absolute value of the estimated glitch waveforms as thresholds. Figure E5 shows the variation of FNRs due to the choices of percentiles. We find that the peak of median FNRs are obtained with and 60 and 55 percentiles for Scattered light and Pandas [78], Matplotlib [79], nds2utils [80], and scikit-learn [81]. This manuscript has been assigned LIGO Document Control Center number LIGO-P2100159.