#  Random projections in gravitational-wave searches of compact binaries II:
efficient reconstruction of the detection statistic

Amit Reza<EMAIL_ADDRESS>Anirban Dasgupta<EMAIL_ADDRESS>Anand
S. Sengupta<EMAIL_ADDRESS>Indian Institute of Technology Gandhinagar,
Gujarat 382355, India.

###### Abstract

Low-latency gravitational wave search pipelines such as GstLAL take advantage
of low-rank factorization of the template matrix via singular value
decomposition (SVD). The matrix factors can be used to ’reconstruct’ the
detection statistic to the desired precision by linearly combining the vectors
obtained from filtering the data against the top-few basis vectors that span
the template vector space. With unprecedented improvements in detector
bandwidth and sensitivity in advanced-LIGO and Virgo detectors, one expects
the size of template banks to increase by several orders of magnitude in
upcoming searches. Naturally, this poses a formidable computational challenge
in factorizing extremely large template matrices. Previously, [Kulkarni et al.
[6]], we had introduced the idea of random projection-based matrix
factorization as a computationally viable alternative to SVD, applicable for
large template banks. In this follow-up paper, we demonstrate the application
of a block-wise randomized matrix factorization (RMF) algorithm for computing
low-rank factorizations at a preset average fractional loss of SNR. This new
scheme is shown to be more efficient in the context of the LLOID framework of
the GstLAL search pipeline. Further, it is well-known that for very large
template banks, the total computational cost of the search is dominated by the
cost of reconstructing the detection statistic as compared to that of
filtering the data. However, the issue of optimizing the reconstruction cost
has not been addressed satisfactorily so far in the available literature. We
show that it is possible to approximately reconstruct the time-series of the
matched-filter detection statistic at a fraction of the total cost using the
matching pursuit algorithm. The combination of the two algorithms presented in
this paper can handle online searches involving large template banks more
efficiently. We have analyzed the total computational cost in detail and offer
various tips for optimally applying the RMF scheme in different parts of the
parameter space. The algorithms presented in this paper are designed in a
suitable manner that can be efficiently implemented over a distributed
computing architecture. Results from several numerical simulations have been
presented to demonstrate their efficacy.

Keywords \- Gravitational-wave (GW), GstLAL, Singular Value Decomposition
(SVD), Random Projection (RP), Randomized matrix factorization (RMF),
Computational Complexity

††preprint: LIGO DCC number LIGO-P2100004

##  I Introduction

The detection of a short gamma-ray burst (GRB-$170817$A) [44] by the Gamma-ray
Burst Monitor (GBM) instrument aboard the Fermi satellite, coincident with the
discovery of the GW-$170817$ [43] event in data from the advanced-LIGO and
Virgo detectors event marks an important epoch in GW-astronomy: that of an era
of multi-messenger observations. More multi-messenger discoveries can provide
a better understanding of sources of GWs, and the underlying astrophysical
processes leading to their progeny. It is obvious that the prompt detection of
the GWs is very crucial to chase the transient events associated with the GWs
sources. Hence the design of efficient and real-time detection pipelines has
become essential.

While some recently proposed machine-learning based detection pipelines have
been demonstrated to be quite efficient at near real-time detection of
astrophysical GW signals, most of the well-established, traditional data
analysis pipelines are based on the matched-filtering [40; 28] scheme. Under
the matched-filtering analysis framework, one computes the cross-correlation
between a large number of theoretically modelled CBC waveforms (called the
template bank) and the detector output, which may contain astrophysical GW
signals from coalescing compact binaries buried in instrumental and
environment noise that couple to the apparatus. The cross-correlation is
computed over the sensitive bandwidth of the detector and weighted inversely
by the noise power spectrum of the detector. For additive Gaussian noise, the
matched filter technique can be shown to be optimum in that it yields the
maximum signal-to-noise ratio (SNR). While the noise in real LIGO or Virgo
detectors is neither wide-sense stationary nor Gaussian, a phenomenological
final detection statistic using the matched-filtering output is employed to
ensure near optimally.

New GW detectors such as LIGO-India, KAGRA are slated to join the global
network of terrestrial GW detectors in the next few years. The combined data
from such a network can not only help reconstruct the physical properties of
the CBC sources with unprecedented accuracy but also help test new and subtle
phenomena. It is reasonable to expect that a larger number of next-generation
GW detectors will bring in new challenges to the matched-filtering based
detection pipelines, especially in their ability at low-latency detection of
transient signals: (a) the most straightforward challenge will arise from the
large volume of data that will be needed to be analyzed in real-time (b)
improvement in detector sensitivity at lower frequencies, and the deployment
of theoretical template waveforms that include precession effects will
increase the number of templates in the bank by several orders of magnitude,
(c) the improvement of detector sensitivity at low frequencies will lead to
many more cycles of GW signals within the detector’s sensitive bandwidth,
which in turn will increase the time duration of these transient signals. For
example, the signal from a typical coalescing binary neutron-star system will
last several tens of minutes, which will have a precipitous effect on the
amount of data samples that need to be analyzed. The issues outlined above set
the stage for the analysis presented in this paper, where we explore
alternative, more economical implementations of the matched filtering scheme
useful for the real-time GW searches from CBC sources.

The development of a low-latency search pipeline (GstLAL) [10; 8; 5] by the
LSC has significantly reduced the computational cost of implementing a matched
filtering search leading to a reduction in latency of detection. This has also
facilitated rapid EM follow-up observations. The GstLAL pipeline calculates
the approximate matched-filtering SNR based on the truncated SVD [37]
framework [1]. The basic idea is as follows: instead of using the templates
directly for filtering the data, a set of basis vectors spanning the template
bank is pre-computed, and only a fraction of the basis vectors (in descending
order of their corresponding singular values) are used to filter the data. In
other words, the top-few basis vectors are used as surrogate templates to
filter the data, which are then combined to approximately reconstruct the
matched-filter SNR time-series. Such an ad-hoc truncation of the bases
inevitably leads to a loss in SNR. Operationally, one works with enough number
of bases to guarantee an average SNR loss
$\langle{\delta\rho}/\rho\rangle\simeq 10^{-3}$. The power of this method
comes from the rapid decay in the spectrum of singular values, which enables
us to use only a very small fraction of bases vectors to reach the desired
accuracy in SNR, thereby leading to a significant reduction in the filtering
cost.

The time complexity of SVD factorization does not scale well with the increase
in matrix size leading to several practical problems in its implementation for
GW searches with very large template banks. It is a formidable computational
challenge to SVD factorize a very large template bank in its entirety. As a
practical workaround, such large template banks are split into sub-banks and
the SVD-based matched filtering scheme is applied independently to each of the
split banks [8; 5]. However, such a strategy diminishes the linear dependence
of the template waveforms, and as a consequence, the anticipated computational
advantage of the SVD method is substantially compromised. In a recent work
[6], we have shown that the computation of the basis vectors for an entire
template bank is possible by implementing a probabilistic low-rank matrix
approximation framework which can adequately address the scaling issue of SVD
factorization of huge template matrices. We have shown theoretically that the
random projection (RP) [36] based matrix factorization technique can be a
practical alternative for computing the set of basis in the GstLAL search
pipeline. Randomized matrix factorization (RMF) algorithms [11; 13] are
computationally more efficient, numerically stable, highly parallelizable, and
can be fine-tuned to work optimally for the problem identified above. It
incurs fewer floating-point operations than the standard deterministic matrix
factorization as it uses a fixed amount of passes over the data.

The average fractional loss of SNR under the truncated SVD method is directly
proportional to the discarded number of singular values. Depending on the
target average fractional loss of SNR, one can decide the required number of
basis vectors. Hence, a SVD based filtering scheme allows using an acceptable
average fractional loss of SNR as a tunable parameter to determine the search
efficiency using discarded set of basis vectors based on the less weighted
singular values.

However, the truncated SVD technique can measure the average fractional loss
of SNR corresponding to the number of discarded basis vectors. Hence, the
whole spectrum of singular values and basis vectors must be used to fix the
cut-off on the average fractional loss of SNR. For a large template bank,
where the number of required basis vectors is significantly less in comparison
to the total number of templates corresponding to a fixed average fractional
SNR loss, the number of singular values and basis vectors that need to be
discarded are is the large fraction of the number of templates. Hence
computation and storing of the whole set of singular values and basis vectors
is a waste of computational resources, which is a computationally expensive
step and slows down the entire analysis, and needs to be optimized to make the
whole process faster. It is more effective to reverse the problem where one
can find a fixed number of basis vectors based on predefined average
fractional SNR threshold, and as a result, there is no need to obtain the less
important singular and basis vectors. But using the truncated SVD framework,
it is hard to address this issue for a large template bank. Hence, in this
work, we introduce block-wise RP based matrix factorization of a template
matrix that can easily handle this kind of issue and is able to find an
optimal number of bases based on a predefined average fractional SNR loss. The
basic idea of block-wise RMF is to obtain a few essential singular values and
basis vectors in each iteration depending on the size of the block of the RP
matrix and check for the optimal criteria on the average fractional SNR loss.
Once it reaches the predefined fractional SNR loss, there is no need to
project further. Mathematically, the block-wise projection helps obtain a
fixed importance sub-space of the template matrix in each iteration. The basis
vectors can be computed from that sub-space only. If the block size is small,
one has to use a small sub-space to calculate the basis vectors in each
iteration.

The SVD based matched filtering scheme’s primary focus has been to reduce the
filtering cost. The set of basis vectors are firstly computed, and filtering
is performed against those bases. Since the number of bases is significantly
less than the number of templates, the computational cost is saved. But, this
is not the complete scenario as here we are discarding the cost of
reconstruction of the SNR time-series. After getting the filter output from
the correlation between the data and basis, we need to multiply it with the
coefficient vectors to obtain the SNR time-series for each template. This
reconstruction is an issue for a large bank as the cost of reconstructing the
coefficient is high. Hence, just reducing the filtering cost does not suffice
as the ultimate solution to this problem. We need to handle the issue of the
reduction of the reconstruction cost too, which is challenging as well as an
open research problem. To the best of our knowledge, there are no research
articles available to address this issue. In this work, we are providing an
approach to reduce the reconstruction cost for the first time. The method
obtains a sparse coefficient corresponding to a specific template waveform
based on the matching pursuit (MP) algorithm. This sparse coefficient can be
used to reconstruct the SNR time-series for each of the templates. The MP
algorithm transforms the original dense coefficient vectors into sparse ones.
Therefore, these transforms can reduce the reconstruction cost as the number
of floating-point operations will be less than the use of the dense
coefficient matrix. Thus, a randomized matrix factorization scheme and
matching pursuit algorithms together can increase the power of this kind of
basis-based filtering scheme, and in a distributed set-up, it can provide an
optimal framework compared to the SVD based filtering scheme.

The rest of the paper is organized as follows. In section-II, we have
summarized the matched filtering scheme for the CBC searches, the SVD based
matched filtering scheme, and the total computational cost analysis of this
scheme. In section-III, the conceptual idea behind the RMF has been described.
The RMF method based on a predefined average SNR loss has been described in
sub-section-III.1 which is followed by a section ( section-IV) in which
reconstruction of the SNRtime-series using MP has been proposed. The optimal
RMF schemes for GstLAL pipeline have been designed in section-V. In the
following sections, details of the advanced RMF method and a new approach for
reducing the reconstruction cost combining RMF with MP have been described.

##  II Compact Binaries searches

The coalescence of compact binaries e.g. neutron stars, black holes, is a
promising source of GW signal. The standard matched filtering technique is
used as an initial step to identify a signal’s presence in a detector output.
The matched filtering operation has to be performed between the detector
output and theoretical waveform. A large number of theoretically well-modelled
waveforms are needed to probe the component mass parameter space for this
matched filtering operation. This set of theoretical waveforms are called the
template waveforms $\bm{h}^{\alpha}(t):\alpha=1,2,\cdots N_{T}$ and can be
thought of as vectors lying in a vector space. The number of time-series
samples determines the dimension of each vector $(N_{s}=f_{s}\times T)$, where
$f_{s}$ and $T$ are the sampling frequency and time, respectively. Suppose the
number of template waveforms for a fixed parameter space is $N_{T}$. Each
template $\bm{h}^{\alpha}(t)$ consists of two orthogonal template
$\bm{h}^{\alpha}_{0}$ and $\bm{h}^{\alpha}_{\pi/2}$ such that

$\bm{h}^{\alpha}(t)=\mathcal{A}_{0}\,\bm{h}^{\alpha}_{0}\big{(}t,\theta\big{)}+\mathcal{A}_{\pi/2}\,\bm{h}^{\alpha}_{\pi/2}\big{(}t,\theta\big{)}\,,$
(1)

where $\theta$ represents the set of intrinsic and extrinsic parameters, hence
a parameter space consisting of $N_{T}$ number of templates can be considered
as a template matrix $\mathbf{H}$ of size $2N_{T}\times N_{s}$, where
normalized whitened template waveforms are stacked row-wise. The template
waveforms are whitened by the power spectral density of the noise. Therefore,
we can define the matched-filter output for a $\alpha^{\texttt{th}}$ template
at a specific time against detector output $\bm{s}$ as follows:

$\rho_{\alpha}(\Delta
t)=\sqrt{\big{(}H_{2\alpha-1}\,\cdot\bm{s}^{T}\big{)}^{2}\,+\,\big{(}H_{2\alpha}\,\cdot\bm{s}^{T}\big{)}^{2}}\,,$
(2)

where, $H_{\alpha}$ denotes the $\alpha^{\texttt{th}}$ row of template matrix
$\mathbf{H}$ and $\bm{s}^{T}$ implies the transpose of the detector output
$\bm{s}$. Note that the Eq.2 shows the matched filter output for a fixed time-
stamp. It is always preferable to compute the correlation between the whitened
template and detector output in a frequency domain wherein the computational
cost of matched filtering with a specific template against data considering
all time-stamps is $N_{s}\,\log N_{s}$. Whereas, to obtain the matched filter
output for all time-stamps in a time domain, the template waveform has to be
shifted by $\Delta t$ every time to obtain the output defined in Eq.2. This
shifted version of each of the rows of $\mathbf{H}$ can be done by
constructing a circulant matrix structure for each of the row vectors. Hence,
in the time domain, the matched filter output for a fixed template with all
the time stamps can be redefined as follows:

$\rho_{\alpha}(\Delta
t)=\sqrt{\big{(}\mathcal{C}(H_{2\alpha-1})\,\cdot\bm{s}^{T}\big{)}^{2}\,+\,\big{(}\mathcal{C}(H_{2\alpha})\,\cdot\bm{s}^{T}\big{)}^{2}}\,,$
(3)

where $\mathcal{C}$ represents a circulant matrix for a specific row of
$\mathbf{H}$. The dimension of $\mathcal{C}$ is $N_{s}\times N_{s}$. Hence,
the computational cost for matched filter output for all time-stamps defined
in Eq.3 is ${N_{s}}^{2}$. But this cost can be reduced to $N_{s}\,\log N_{s}$
as circulant matrix can be factored into discrete Fourier components. Hence,
the computational cost to obtain the matched output for all time-stamps using
frequency domain or time domain representation is the same, and it requires
$N_{s}\,\log N_{s}$ floating-point operations. So, the total computational
time of matched filtering scheme considering all the templates is
$N_{T}\,N_{s}\,\log N_{s}$. In such a search, for a large template matrix
$\mathbf{H}$, for example, one having the number of rows (templates) and the
number of columns (time samples) of $\mathcal{O}(10^{5})$,
$\mathcal{O}(10^{6})$ respectively, the matched filtering task becomes very
expensive as the number of floating-point operations will be of order
$10^{12}$. For the upcoming searches, it is expected that the number of
template waveforms and the number of time samples will increase significantly,
e.g., it can be of order $\mathcal{O}(10^{5})$, $\mathcal{O}(10^{6})$
respectively. Therefore the computational cost of performing matched filter
operation will also be increased, which eventually makes the search process
slow and tedious. The SVD based matched filtering scheme [1] can reduce the
computational cost of matched filtering. In the SVD based matched filtering
scheme, a set of essential basis vectors have been computed first from a group
of template waveforms. Further, data (i.e., detector output) is filtered
against these important bases only. Since the number of essential bases is
less than the total number of template waveforms, thus it can reduce the
computational cost of performing matched filters significantly. The basis
representation of the template waveforms and matched filter operations in
terms of basis and data is shown briefly in the next section.

###  II.1 Singular Value decomposition based Matched Filtering

SVD can be used to compute the matched filtering output more efficiently.
Instead of calculating the cross-correlation between the data vector $\bm{s}$
and template vector $\bm{h}^{\alpha}$ directly, the idea is to compute a set
of orthogonal basis from template vectors using SVD and perform cross-
correlation of the data stream with those set of basis vectors.

Using SVD, one can numerically compute a set of orthogonal basis
($\bm{v}_{\nu}$) in such a way that a specific linear combination of these
basis vectors can represent the template vectors. After applying SVD, the
template matrix $\mathbf{H}_{2N_{T}\times N_{s}}$ can be decomposed into three
special matrices $\mathbf{U}$, $\mathbf{\Sigma}$, and $\mathbf{V}$, where
$\mathbf{V}$ and $\mathbf{U}$ consist of a set of orthonormal vectors and
$\mathbf{\Sigma}$, the diagonal matrix contains the singular values in
descending order of magnitude. $\mathbf{V}$ is a matrix of orthonormal bases
whose each column are basis vectors.

$H_{\mu j}=\sum_{\nu=1}^{2N_{T}}{u_{\mu\nu}\sigma_{\nu\nu}\,{v_{\nu
j}}^{T}}=\sum_{\nu=1}^{2N_{T}}{A_{\alpha\nu}\bm{v}_{\nu}}\,,$ (4)

where $\bm{v}_{\nu}=\\{\bm{v}_{1},\bm{v}_{2},...,\bm{v}_{2N_{T}}\\}$ is a set
of basis vector. $A_{\alpha\nu}$ is the reconstruction matrix and each row
defines the corresponding weights of the basis vectors for representing the
row vectors of template matrix $\mathbf{H}$. From Eq.(4), it is clear that the
required number of basis vectors for representing any row vector of the
template matrix $\mathbf{H}$ is equal to the full-rank of the template matrix
i.e.
$2N_{T}=\textrm{min}\big{(}2N_{T},N_{s}\big{)},\,\text{if}\,2N_{T}<N_{s}$.

For the non-precession and aligned spin waveform, the neighbouring template
vectors are differed by phase only. Hence all of them are almost similar.
Therefore, only a few (let us suppose $\ell$, where $\ell\ll 2N_{T}$) number
of basis vectors are sufficient to reconstruct the template vectors, and
indeed it will approximate the matched filter output with high precision.
These $\ell$ number of basis vectors are the top most important singular
directions, ranked based on the top-$\ell$ singular values. Hence it is
possible to approximately reconstruct the template matrix
$\mathbf{H}^{\prime}$ based on this set of bases. The reconstructed matrix
$\mathbf{H}^{\prime}$ can be written as follows:

$H^{\prime}_{\mu j}=\sum_{\nu=1}^{\ell}{A_{\alpha\nu}\bm{v}_{\nu}}$ (5)

Exclusion of $(2N_{T}-\ell)$ number of less important basis vectors culminates
the approximation in the length of each template vector. Hence, the overall
effect will be observed on the preservation of the matrix norm (energy) of
$\mathbf{H}$, which indeed turns out to find out the optimal value of $\ell$
for which the approximated template matrix ($\mathbf{H}^{\prime}$) captures
the maximum energy, i.e. $\|\mathbf{H}-\mathbf{H}^{\prime}\|_{F}<\epsilon$.
The spectra of the singular values (diagonal elements of the matrix
$\mathbf{\Sigma}$) can help dictate the right choice of $\ell$. If singular
values fall sharply, then it is enough to take
$\ell\ll\textrm{min}\big{(}2N_{T},N_{s}\big{)}$, and based on the cut-off on
the spectra, one can compute the value of $\ell$. This reduction of the number
of basis vectors directly affects the matched filtering output $\rho(\Delta
t)$. Thus the reconstruction of $\rho(\Delta t)$ based on topmost singular
values can be defined as follows:

$\rho^{\prime}_{\alpha}(\Delta
t)=\sqrt{\sum_{\nu=1}^{\ell}{A_{(2\alpha-1)\nu}\big{(}\bm{v}_{\nu}\cdot\bm{s}^{T}\big{)}}+\sum_{\nu=1}^{\ell}{A_{2\alpha\nu}(\bm{v}_{\nu}\cdot\bm{s}^{T})}}$
(6)

The expected average fractional SNR loss ($\langle{\delta\rho}/\rho\rangle$)
can be defined as a function of truncated singular values as defined in [1].

$\langle{\delta\rho}/\rho\rangle=\frac{1}{4N_{T}}\,\Bigg{(}\frac{\sum_{\nu=\ell+1}^{2N_{T}}{\sigma_{\nu}^{2}}}{\sum_{\nu=1}^{2N_{T}}{\sigma_{\nu}^{2}}}\Bigg{)}$
(7)

Therefore average expected fractional loss of SNR could be used as a tuning
parameter to fix a certain threshold for deciding the number of reduced basis
vectors to reconstruct the original template waveforms. From Eq.(6), it is
clear that we can approximate the SNR time-series calculation corresponding to
each template waveform using only $\ell$ number of filtering operations
between the top-$\ell$ basis vectors and the detector output. The time
complexity of the filtering is $\mathcal{O}(\ell\,N_{s}\,\log N_{s})$, is very
less in comparison to the required time complexity of the direct correlation
between template vectors and data vectors because of
$\ell\ll\textrm{min}\big{(}2N_{T},N_{s}\big{)}$. In this way, reducing the set
of basis vectors helps to reduce the total number of the matched filtering
operations and, indeed, the time complexity of filtering. Note that this
computational cost excluded the cost of reconstructing the SNR time-series
after performing the filtering operation of the bases against data. From Eq.6,
it is clear that to obtain every SNR time-series corresponding to each
template waveform; one has to multiply with the corresponding coefficient
vectors $A_{2\alpha-1}$ and $A_{2\alpha}$. For a small template matrix
$\mathbf{H}$, the reconstruction cost is negligible, but for a large template
matrix, it is not suggested to ignore the computational cost of reconstruction
of the SNR time-series. Hence, we have demonstrated the complete cost analysis
in the next section, including reconstruction cost and matrix factorization
cost for the SVD based matched filtering scheme.

###  II.2 Computational cost analysis of SVD-based matched filtering approach

As mentioned above, the total computational cost of matched filtering using
SVD can be decoupled as a sum up of two separate costs.

  1. 1.

Filtering cost is associated with calculating the cross-correlation values of
the set of basis vectors with the data vectors. The total cost amounts to
$\ell\,N_{s}\,\log N_{s}$.

  2. 2.

Reconstruction cost is associated with the reconstruction of SNR time-series
for each template. This involves multiplying each filter output with the
corresponding coefficient vector. The total reconstruction cost is
$2N_{T}\,\ell\,N_{s}$, as the coefficient matrix
$A_{\alpha\nu}\in\mathbb{R}^{2N_{T}\times\ell}$, and the filter output matrix
$\\{\mathbf{V}\cdot\bm{s}^{T}\\}\in\mathbb{R}^{\ell\times N_{s}}$.

Apart from the above two costs, there is another cost for performing matrix
factorization of the template matrix $\mathbf{H}$ using SVD. This
computational cost is $\mathcal{O}(N_{s}\,{N_{T}}^{2})\,:\,2N_{T}<N_{s}$. Note
that matrix factorization cost can be thought of as off-line cost, which is a
part of data pre-processing, whereas the filtering and reconstruction cost is
an online cost. But, it is crucial that this cost is included because for a
large template matrix, if the factorization cost is huge, then computation of
the basis vectors beforehand becomes an impossible task. For example, for a
template matrix $\mathbf{H}$ of size $10^{5}\times 10^{6}$, factorization
using SVD is also unbearable as the time complexity will increase drastically
as $\mathcal{O}(10^{16})$. Also, applying SVD to such a huge matrix requires a
large run-time memory space which is also difficult to obtain. Therefore, it
is clearly an issue because of the required time complexity and large memory
space.

It is clear from the above discussion that the SVD based matched filtering
approach adds an extra computational burden in terms of the reconstruction
cost of the SNR time-series, whereas the direct matched filtering between
template and data has only the filtering cost. So, comparing only the
filtering cost for both methods is inequitable. The comparison of filtering
cost using these two approaches is only reasonable if the reconstruction cost
in case of SVD based matched filtering scheme is negligible, which is only
possible if the number of templates in a template bank is few, which is an
improbable scenario for the upcoming CBC searches. Therefore, it is crucial to
investigate the possible ways to reduce the reconstruction cost.

If the template matrix contains a large number of template waveforms, then the
size of the coefficient matrix is large; therefore, the reconstruction cost of
the SNR will be very high compared to the filtering cost, and hence
reconstruction cost will dominate over the filtering cost. For example,
consider size of $\mathbf{H}$ as $10^{5}\times 10^{6}$. Now if $\ell=10^{4}$,
then the filtering cost becomes $\mathcal{O}(10^{10})$, whereas the
reconstruction cost of SNR time-series becomes $\mathcal{O}(10^{15})$. Hence,
the reconstruction cost is $10^{5}$ times the filtering cost. In fact, for any
large template matrix, the reconstruction cost is always dominant no matter
how many basis vectors are considered. In that scenario, only reducing the
filtering cost will not be sufficient. Hence, it is a constraint on the matrix
factorization-based matched filtering scheme. Currently, no method is
available to reduce this reconstruction cost. In this work, for the first
time, we prescribe an efficient solution procedure to address this issue, as
is presented in detail in section IV.

Additionally, to calculate the $\langle{\delta\rho}/\rho\rangle$ after
considering a fixed number ($\ell$) of importance basis vectors for the
approximation of SNR time-series, firstly, it is required to compute all the
basis vectors ($2N_{T}$) and their corresponding singular values. It is then
required to fix a threshold for $\langle{\delta\rho}/\rho\rangle$ and choose
those many singular values for which the threshold can be achieved. Only those
sets of basis are considered as an important basis to compute the SNR time-
series. In practice, generally $\ell\ll 2N_{T}$. Rest of $2N_{T}-\ell$ bases
are considered as less important, and hence it is discarded from the final set
of important bases. As these sets of bases are not used for further
computation of SNR time-series calculation, computation of those bases is
wastage of computational resources and wastage of large memory space required
to store them. For a large template bank, storing those unimportant bases
takes ample memory space even when they are not of any specific use. Only the
corresponding singular values are used to evaluate the corresponding
$\langle{\delta\rho}/\rho\rangle$. For illustrating the wastage of memory
space to store unimportant bases, let us consider a template bank with a
$10^{5}$ number of template waveforms, such that each waveform has $10^{6}$
time-stamps. Then the required space to store the whole template matrix is
$200\,\texttt{GB}$. Now for such a template bank, if hypothetically, one
considers that the required number of essential basis vectors $\ell=10^{4}$,
then, the memory needed for storing the unimportant ($2N_{T}-\ell=19\times
10^{4}$) basis vectors is $190\,\texttt{GB}$. For this specific example,
$95\%$ of the memory is unnecessarily occupied by those bases which are not
used in the SVD based matched filtering scheme. Also, if only $5\%$ bases are
required for the matched filtering calculation, then using SVD decomposition,
one should compute those many bases. In this way, computation costs can be
further reduced. However, in this method, the whole set of bases needs to be
computed to obtain a $\langle{\delta\rho}/\rho\rangle$. Therefore, although it
is not optimal to compute all the bases computationally, it is essential to
calculate $\langle{\delta\rho}/\rho\rangle$ and obtain a corresponding $\ell$.
Thus in the case of the SVD based matched filtering strategy, one has to
compute all the basis vectors knowing the facts that $80-95\,\%$ of the bases
are useless and need to be discarded from the final list.

It is a real challenging problem constraining the limited computational
resources. Ideally, one has to pre-defined a $\langle{\delta\rho}/\rho\rangle$
and compute a fixed set of bases that can satisfy it. In this way, the
utilization of computational resources will be optimized. However, this kind
of reverse calculation is not feasible using SVD factorization. Addressing
this problem using truncated-SVD set-up is hard. This crucial problem has to
be sorted out unless the computational resources will be misapplied for a
large template bank. As a result, designing the real-time detection pipeline
can not be succeeded. Hence, though the SVD based matched-filtering approach
is the best-known approach for low-latency CBC searches, it has some specific
limitations for a large template bank, as mentioned above. One can summarize
the limitations as follows:

  1. 1.

Large template matrix can not be decomposed using SVD due to the high
computational cost of decomposition.

  2. 2.

For a large template matrix, the reconstruction cost is the dominant cost.
Hence, only reducing the filtering cost is not sufficient.

  3. 3.

Performing entire SVD and then choosing $\ell$ bases to attain predefined
$\langle{\delta\rho}/\rho\rangle$ is inefficient.

To handle the first two problems, the current low-latency CBC search pipeline
(GstLAL) divides the full bank into a group of sub-banks, and for each sub-
bank, the SVD based matched filtering is performing independently. However,
this is a sub-optimal solution as an optimal way of dividing the whole bank
into sub-banks is an open problem. Secondly, dividing into sub-banks may cause
for losing the linear dependency between the nearby waveforms. As a result, it
increases the number of essential basis vectors for a fixed
$\langle{\delta\rho}/\rho\rangle$, which indeed increases the filtering cost.
In our previous work [6], we have computed $\beta$, the ratio of the number of
bases summed across all the sub-banks to the number of basis from the SVD
factorization of the full bank. In Figure-$1$ [6], we have plotted $\beta$
against $\langle{\delta\rho}/\rho\rangle$ for six different template bank. The
figure shows that the number of important bases for a whole template bank is
less compared to the sum of the bases considering all the sub-banks for a
fixed $\langle{\delta\rho}/\rho\rangle$. The value of $\beta$ increases with
the increasing size of the template bank. This is because the computation of a
global set of the basis for a whole template bank is more relevant than
computing a local set of bases based on the different sub-banks of the entire
template bank as the number of the basis for the previous case is less (see
Figure-$1$ of [6] for more details). Hence, to reduce the filtering cost
optimally, one has to compute the global set of bases by directly applying SVD
to the entire template bank. But that is not possible because the
computational cost for the factorization for an entire template bank is high.
Hence, each sub-bank matrix factorization is the best-known solution for the
SVD based matched-filtering scheme. If SVD factorization for a large template
bank is feasible, we may resolve the matched filtering cost optimally compared
to the overall filtering cost combining all the sub-banks factorization
schemes. However, for a large template bank, matched filtering cost is not
only the cost. The reconstruction of the SNR time-series is also
computationally expensive. In fact, for a large template matrix,
reconstruction cost is the dominant cost. Therefore, in the current scheme, to
reduce the reconstruction cost, one must divide the full bank into sub-banks.
As for each sub-banks, the reconstruction cost is minimal. However, the
overall filtering cost is increased due to make a large number of sub-banks.
Hence, there is a clear trade-off between the reconstruction cost and
filtering cost for a fixed bank. Since the problem of reduction of the
filtering cost and the reconstruction cost is complementary. It is not
possible to obtain the optimal solution for both issues simultaneously using
any matrix factorization (e.g. SVD) based matched-filtering scheme. The
division into sub-banks can reduce the reconstruction cost, whereas it can
increase the filtering cost as the number of bases increases by a few factors.
The number of basis vectors is directly proportional to the size of the bank.
Therefore it can reduce the filtering cost if one computes the global set of
bases by factorizing the entire bank instead of obtaining the basis vectors
locally from the sub-banks. The optimal time complexity can be related to the
division of the bank. Not only that, the whole comparison depends on the
number of sub-banks, the size of the template bank, and the number of desired
top basis vectors. Therefore, it is crucial to analyze the time complexity
calculation in detail for both cases by properly considering all these
factors. In section F, we have shown a numerical analysis of the optimized way
of division into sub-banks.

In our previous work [6], we have shown that the RP based matrix factorization
technique can efficiently compute the set of basis and the coefficient matrix
for a large template bank, thus reducing the filtering cost. In this work, we
have used a further improved version of the RP based matrix factorization by
which it is possible to identify the value of $\ell$ corresponding to the pre-
defined $\langle{\delta\rho}/\rho\rangle$. Further, we demonstrate the
reduction of the reconstruction cost using the Matching Pursuit (MP) algorithm
[46]. Finally, we prescribe an algorithm that combines the matrix
factorization of a template matrix with a fixed
$\langle{\delta\rho}/\rho\rangle$ and reduces the reconstruction cost of SNR
time-series of each template waveforms.

##  III Random projection-based template matrix factorization

In our previous work [6], we have demonstrated the utility of the RP-based
matrix factorization for the GW searches from the compact binary coalescence.
Consider a template matrix $\mathbf{H}_{2N_{T}\times N_{s}}$, $\exists$ a
$\ell$-dimensional dominant subspace
$\bar{\mathbf{H}}_{2N_{T}\times\ell}=\mathbf{H}\,\mathbf{\Omega}\in\mathbb{R}^{\ell}$.
All the row vectors of the template matrix $\mathbf{H}$ have been projected
into a $\ell$-dimensional space using an RP matrix $\mathbf{\Omega}$. Hence
the dimension of each of the row vectors reduced from $N_{s}\rightarrow\ell$,
where $\ell\ll N_{s}$. Thus, for a large template matrix $\mathbf{H}$, to
minimize the factorization cost, it is optimal to use $\bar{\mathbf{H}}$ to
obtain the basis of the original matrix as it almost preserves the geometrical
structure of the row-space. Hence $\ell$-dimensional representation of the row
vectors of $\mathbf{H}$ can be used to compute the basis vectors $\mathbf{Q}$,
which are nearly in the same direction as the top-$\ell$ eigenvector
directions of the original eigenvector of the row space of $\mathbf{H}$.
Concisely, RP involves taking the projection of a high-dimensional vector to
map it into a lower-dimensional space while providing some guarantees on the
approximate preservation of pair-wise distance between the vectors. The mere
fact that the idea of RP came directly from Johnson-Lindenstrauss lemma [33]
guarantees the preservation of the pair-wise distance with a certain accuracy
between a set of points which are projected to a lower-dimensional space from
a higher dimensional space using RP operator. The theoretical and practical
bound of the distortion factor with an example is shown in the subsection A of
section VI. This RP based matrix factorization scheme also provides similar
factors to obtain rank-$\ell$ approximation ($\mathbf{H}^{(\ell)}$) as using
SVD except for the fact that it factorizes the template matrix into two
factors $\mathbf{Q}$ and $\mathbf{B}$ whereas SVD factorizes it into three
different factors. We can use $\mathbf{B}$ as a surrogate template to filter
against the data vector $\bm{s}$. Note that the low-rank value ($\ell$) is
defined as the dimension of the projected lower-dimensional space; hence, all
the template waveforms are projected from $N_{s}\rightarrow\ell$ using RP
operator $\mathbf{\Omega}$ and hence the set of basis $\mathbf{Q}$ has been
computed using those template waveforms embedded in a lower-dimensional space.
Therefore, the obtained set of basis $\mathbf{Q}$ can be approximated using
the top-$\ell$ basis vectors of $\texttt{range}\,(\mathbf{H})$. The
coefficient matrix ($\mathbf{B}$) can be obtained by projecting all the row
vectors of $\mathbf{H}$ onto the set of basis vectors $\mathbf{Q}$, i.e.
$\mathbf{B}=\mathbf{Q}^{T}\,\mathbf{H}$. Hence, we can approximate
$\mathbf{H}_{2N_{T}\times
N_{s}}\approx\mathbf{Q}_{2N_{T}\times\ell}\,\mathbf{B}_{\ell\times N_{s}}$.
Therefore approximated the SNR time-series
($\rho^{\prime\prime}_{\alpha}(\Delta t)$) for each of the templates can be
computed as follows:

$\rho^{\prime\prime}_{\alpha}(\Delta
t)=\sqrt{\sum_{\nu=1}^{\ell}{Q_{(2\alpha-1)\nu}\,\big{(}\bm{B}_{\nu}\cdot\,\bm{s}^{T}\big{)}(\Delta
t)\,+\,\sum_{\nu=1}^{\ell}{Q_{(2\alpha)\nu}\,\big{(}\bm{B}_{\nu}\,\cdot\bm{s}^{T}\big{)}(\Delta
t)}}}$ (8)

Based on the projected dimension $\ell$, one can compute the corresponding
$\langle{\delta\rho}/\rho\rangle$ (See [6]). Hence, $\mathbf{Q}\,\mathbf{B}$
factorization guarantees
$\|\mathbf{H}-\mathbf{Q\,B}\|_{F}\leq\langle{\delta\rho}/\rho\rangle$. It is
notable that in this case the obtained $\langle{\delta\rho}/\rho\rangle$ can
not be predefined before factorization. However, one can also project the
column vectors of $\mathbf{H}$ onto a $\ell$ dimensional space by using a RP
matrix $\mathbf{\Omega}$ of dimension $2N_{T}\times\ell$. In that case, the
factorization can be defined as $\mathbf{H}_{2N_{T}\times
N_{s}}\approx\mathbf{B}_{2N_{T}\times\ell}\,{\mathbf{Q}^{T}}_{\ell\times
N_{s}}$. The details of the algorithm is described in appendix C. This kind of
factorization is required for the construction of the sparse coefficient
matrix using matching pursuit algorithm, described in section IV.
Additionally, the $\ell$ dimensional sub-space formation is possible by
projecting both the row and column vectors of $\mathbf{H}$ using two RP
matrices. This kind of compression of rows and columns is useful if both the
number of rows and columns of the template matrix are large.

Figure 1: The figure shows the pictorial description of the RMF scheme with a
fixed rank $\ell$ for a template matrix $\mathbf{H}_{2N_{T}\times N_{s}}$ (See
Alg-3). In this scheme, the computation of the orthonormal basis for a
template matrix has been done from a lower-dimensional representation of the
template matrix, denoted as $\bar{\mathbf{H}}$, which is computed by
projecting all the template waveforms to the column vectors of a RP matrix
$\Omega$. The steps are shown in the stage-$1$ of this figure. Since, $\ell\ll
N_{s}$, hence the lower-dimensional representational of the template matrix is
a much smaller matrix as compare to the original template matrix. Hence,
computational expenses can be reduced using this intermediate representation
of the template waveforms for a large template matrix. In Stage-$2$, the
computation of the coefficient matrix has been shown after obtaining the basis
vectors. The Stage-$3$ is an optional step, as it is only useful to obtain the
singular value spectra of the coefficient matrix $\mathbf{B}$, which
approximates the top-$\ell$ singular values of the template matrix
$\mathbf{H}$. The rank-$\ell$ truncation of $\mathbf{H}$ based on SVD can be
equivalently written based on the factors obtained from this RMF method as
$\mathbf{H}_{2N_{T}\times
N_{s}}\approx\mathbf{U}_{2N_{T}\times\ell}\,\Sigma_{\ell\times\ell}\,{\mathbf{V}^{T}}_{\ell\times
N_{s}}$. Similarly, it can be decomposed as $\mathbf{H}_{2N_{T}\times
N_{s}}\approx\mathbf{Q}_{2N_{T}\times\ell}\,\mathbf{B}_{\ell\times N_{s}}$
using RMF. Further, the factors obtained from SVD and RMF can be comparable if
we decompose the coefficient matrix $\mathbf{B}$ as shown in the Stage-$3$ of
this figure. Therefore, one can approximate the left singular basis
$\mathbf{U}_{2N_{T}\times\ell}\approx\mathbf{Q}_{2N_{T}\times\ell}\,{\mathbf{U}_{\mathbf{B}}}_{\ell\times\ell}$.
Similarly, the top-$\ell$ singular values of $\mathbf{H}$ can be approximated
using the singular values of $\mathbf{B}$, i.e.
$\Sigma_{\ell\times\ell}\approx{\Sigma_{\mathbf{B}}}_{\ell\times\ell}$. The
right singular basis can also be approximated similar way:
$\mathbf{V}_{N_{s}\times\ell}\approx{\mathbf{V}_{\mathbf{B}}}_{N_{s}\times\ell}$.

The RP based matrix factorization has the following two crucial benefits:

  * •

It reduces the time and space complexity of decomposition into factors since
the number of dimensions is now quite manageable.

  * •

Due to the involvement of the simple computational steps, it is easy to make a
simple workflow over a high-throughput computing (HTC) environment as well as
over distributed-memory architecture, e.g., High-performance computing (HPC).

Notably, the matrix factorization scheme described in Alg-3 computes
top-$\ell$ basis vectors based on some fixed value of $\ell$. Hence,
computation of the corresponding $\langle{\delta\rho}/\rho\rangle$ is only
possible after obtaining the top-$\ell$ singular values. As
$\langle{\delta\rho}/\rho\rangle$ is the only tuning parameter for adjusting
the approximated SNR time-series obtained by bases-based matched filtering
scheme. Therefore it is crucial to design an algorithm that can provide a set
of basis vectors based on a pre-defined $\langle{\delta\rho}/\rho\rangle$. SVD
or Alg-3 based matched filtering scheme is inefficient to do that. Therefore
in the next section, we describe an algorithm used to obtain the matrix
factors corresponding to a pre-defined $\langle{\delta\rho}/\rho\rangle$.
Mainly the algorithm is an extended version of Alg-3 in which some optimal
rank-$\ell$ approximation of the template matrix has been computed based on
fixed $\langle{\delta\rho}/\rho\rangle$ iteratively. Identification of optimal
$\ell$ based on a pre-defined $\langle{\delta\rho}/\rho\rangle$ is also
possible iteratively using SVD factorization by obtaining SVD factors in each
iteration. But for a large template matrix, where applying SVD for once is
computationally expensive, using SVD over many iterations is not advisable due
to the limitation of the computational resources. However, the randomized
Alg-1 is automated and also easily operative on distributed memory
architectures.

###  III.1 RP-based template matrix factorization with a fixed average SNR-
loss

From the previous section, it is clear that the standard-RMF scheme (
described in Figure-1, and also in Alg-3) only provides a predefined number of
basis vectors, and hence there is no such amenability on the approximated
error. Therefore, to control the corresponding error due to the factorization,
Alg-3 needs to be modified. Due to this different way of defining the problem
statement, RMF algorithms [13; 11; 14] can be classified into the following
two classes.

  1. 1.

Fixed-rank based factorization: The user provides the predefined rank $\ell$
to obtain the matrix factors. It means that the number of required basis
vectors is already fixed before performing the matrix factorization scheme.
Hence, the dimension of the projected space is already predefined, and one can
do the matrix factorization in that specific dimensional space only.
Therefore, the tolerance error due to the factorization will depend on the
user-defined rank $\ell$. Alg-3 employs this kind of factorization, where the
user-defined rank of the matrix is $\ell$, and all the higher dimensional
vectors are projected onto $\ell$ dimensional space using the RP-operator.

  2. 2.

Auto-rank based factorization: For this class of algorithms, the numerical
rank of a data matrix’s approximated factors is not decided beforehand. The
algorithm is designed in such a way that the rank corresponding to a specific
preset precision error can be revealed automatically. These kinds of
algorithms are greedy by nature, wherein the initial guess about the rank has
to be made. It is better to start with a small value as a guess of the rank
and compute the factorization based on that rank by projecting all the vectors
onto that projected space. One can continue the process until the precision
error matches the predefined error. This is a kind of block-wise iterative
matrix factorization scheme in which, in every iteration, we are obtaining
some essential basis and computing the precision error to decide how many
basis vectors are sufficient to get the desired precision error. Auto-rank-
based factorization can automatically determine the desired rank based on the
constraint on the error tolerance. Here, for a template matrix factorization,
we have set $\langle{\delta\rho}/\rho\rangle$ as the error-tolerance.

Input: Template Matrix {$\mathbf{H}\in\mathbb{R}^{2N_{T}\times
N_{s}}:2N_{T}\leq N_{s}$}.

Output: $\mathbf{Q}_{2N_{T}\times\ell}$, $B_{\ell\times N_{s}}$

1 for  _i = 1,2,3, $\cdots$_ do

2 $\mathbf{\Omega}_{i}\in\mathbb{R}^{N_{s}\times
b}\;:\mathbf{\Omega}_{i}\;\in\mathcal{N}\big{(}0,1\big{)}\;$

$\bar{\mathbf{H}}=\mathbf{H}\,\mathbf{\Omega}_{i}$

// random sampling of column space of $\mathbf{H}$

3

$\mathbf{Q}_{i}=\textrm{qr}\,\big{(}\mathbf{Q}_{i}-\sum_{j=1}^{i-1}{\mathbf{Q}_{j}\,{\mathbf{Q}_{j}}^{T}\,\mathbf{Q}_{i}}\big{)}$

// Re-orthogonalization

4

5 $\mathbf{B}_{i}={\mathbf{Q}_{i}}^{T}\,\mathbf{H}$

$\mathbf{H}=\mathbf{H}-\mathbf{Q}_{i}\,\mathbf{B}_{i}$

// Null space projection

6

7 if  _$\|\mathbf{H}\|_{F}\geq\langle\frac{\delta\rho}{\rho}\rangle$_ then

$\mathbf{Q}=\big{[}\mathbf{Q}_{1}|\mathbf{Q}_{2}|\mathbf{Q}_{3}|\cdots|\mathbf{Q}_{i}\big{]}$

// Column-wise stacking

8

$\mathbf{B}=\big{[}{\mathbf{B}_{1}}^{T}|{\mathbf{B}_{2}}^{T}|\cdots|{\mathbf{B}_{i}}^{T}\big{]}$

// Row-wise stacking

9

10

11 end if

12

13 end for

Algorithm 1 RMF with a fixed error-$\epsilon$( Block-wise RMF)

Alg-1 [14] employs the RP of the row vectors iteratively to obtain the factors
within a fixed error bound to determine the approximated numerical rank of a
data matrix automatically. It is inspired by the column pivoting Gram-Schmidt
scheme and combines this scheme along with the random sampling (projection),
the blocking to obtain the $\mathbf{Q}\,\mathbf{B}$ factorization. Instead of
taking the random projection of the row vectors into $\ell$-dimensional
subspace directly, a few numbers of blocks $b<\ell$ of RP matrix have been
generated. Firstly, the rows of the template matrix have been projected onto a
$b$ dimensional space, and then partial QR decomposition has been used to
obtain the orthogonal basis of that specific ($b$ dimensional) subspace. It is
clear that for the first iteration the number of basis vectors will be $b$ and
the corresponding $\mathbf{Q}^{(i)}\,\mathbf{B}^{(i)}:i=1,2,3,\cdots,$
represents a $b$ rank approximation of $\mathbf{H}$. After obtaining
$\mathbf{Q}^{(i)}$ and $\mathbf{B}^{(i)}$, it is easy to verify the
corresponding $\langle{\delta\rho}/\rho\rangle$. The process continues until
it converges to a pre-defined $\langle{\delta\rho}/\rho\rangle$. Iterative
improvement in the set of orthogonal basis vectors
($\mathbf{Q}=\big{[}\mathbf{Q}^{(1)}|\mathbf{Q}^{(2)}|\cdots|\mathbf{Q}^{(i)}\big{]}$)
and dense skeleton matrix
($\mathbf{B}=\big{[}\mathbf{B}^{(1)}|\mathbf{B}^{(2)}|\cdots|\mathbf{B}^{(i)}$)
has been taken into account until it reaches the desired accuracy. By choosing
a suitable value of $\langle{\delta\rho}/\rho\rangle$, the factorization can
be utilized to reveal the value of $\ell$ which provides an optimal
rank-$\ell$ approximation of $\mathbf{H}$.

The RP matrix
$\mathbf{\Omega}_{N_{s}\times\ell}=\Big{[}\mathbf{\Omega^{(1)}_{N_{s}\times
b}}|\mathbf{\Omega^{(2)}_{N_{s}\times
b}}|\cdots|\mathbf{\Omega^{(p)}_{N_{s}\times b}}\Big{]}$ can be thought of as
a collection of $p=\lceil\ell/b\rceil$ disjoint sub-random matrices, where the
dimension of each sub-random matrices is $N_{s}\times b$. Therefore it turns
out to be optimally finding out of the value of $p$ corresponding to a fixed
$\langle{\delta\rho}/\rho\rangle$. Each of these small blocks can be used to
project the template matrix $\mathbf{H}$ into $b$ dimensional space
($\bar{\mathbf{H}}_{2N_{T}\times b}=\mathbf{H}\,\mathbf{\Omega}$). As the
dimension of the projected matrix $\bar{\mathbf{H}}$ is small, thus the
computation cost for obtaining $b$ number of orthonormal column of
$\mathbf{Q}$ using a standard QR decomposition will be inexpensive. The
corresponding $b$ rows of the coefficient matrix $\mathbf{B}$ can also be
computed using $\mathbf{Q}$. After getting the factors, the template matrix
$\mathbf{H}$ is updated by projecting it out perpendicular to the $b$ basis
vectors, described in step-$7$ of Alg-1.

This step is computationally expensive for a large template matrix, as in each
iteration, it is required to access the whole template matrix. Thus, we can do
this step optimally in a distributed architecture. It is an essential step as
it updates the relative error in the matrix approximation in terms of the
Frobenius norm. The latter can be easily related to
$\langle{\delta\rho}/\rho\rangle$. At this point, we need first to check
either the target accuracy is reached or not. If not, another pass is made
through the updated $\mathbf{H}$ using subsequent blocks of the RP matrix. It
is expected, after a few iterations, one may achieve the optimal rank $\ell$.
While the block-projection scheme described above can lead to computational
advantages, it can also lead to the aggregation of round-off errors to compute
the basis. Also, a block of basis vectors from a specific iteration is
generally not orthogonal to the block of basis vectors obtained from another
iteration. Therefore, to ensure the orthogonality of the block of the basis
vectors with previously obtained basis vectors, we need to incorporate a Gram-
Schmidt like re-orthogonalization procedure to construct the final set of
orthonormal basis vectors $\mathbf{Q}$. This operation is shown in step $5$ of
Alg-1. This re-orthogonalization step is again computationally expensive.
Hence for a large template matrix, this step needs to be performed in a
distributed way.

Figure 2: The figure shows the estimated number of basis corresponding to a
fixed value of $\langle{\delta\rho}/\rho\rangle$. The estimation of basis
vectors has been computed using the block-wise RMF Alg-1. The profile shows
that number of required basis increases as the value of
$\langle{\delta\rho}/\rho\rangle$ gets smaller, which is expected. A template
matrix $\mathbf{H}$ of size $(2N_{T}×N_{s})\equiv 9130\times 65536$ is used in
this analysis. This template bank has been constructed using non-spinning
template waveform for component mass parameters $\big{(}m_{1,2}\big{)}$ in the
range $2.5M_{\odot}\leq m_{1},m_{2}\leq 17.5M_{\odot}$. In each iteration, a
block of size $b=10$ is used to project the template waveform into the lower
dimensional space. In every iteration, the average SNR loss is calculated as
$\langle{\delta\rho}/\rho\rangle\leq\frac{\|\mathbf{H}-\mathbf{H}^{(\ell)}\|_{F}}{\|\mathbf{H}\|_{F}}$.
The computation of the Frobenius norm of the original template matrix
($\|\mathbf{H}\|_{F}$) has been pre-computed, and Frobenius norm of
$\|\mathbf{H}-\mathbf{H}^{\ell}\|_{F}$ is computed in each iteration. The
computation of $\mathbf{H}-\mathbf{H}^{(\ell)}$ is done in step-$8$ of
Algorithm-1.

##  IV Construction of sparse coefficient matrix using Matching pursuit
algorithm

We already have seen that the rows of the template matrix $\mathbf{H}$ can be
represented as the linear combination of the set of basis vectors, which can
be obtained from SVD or RMF. If we consider top-$\ell$ basis vectors, then the
SVD decomposition of the template matrix can be shown as follows:

$\mathbf{H}_{2N_{T}\times
N_{s}}\approx\mathbf{U}_{2N_{T}\times\ell}\,\mathbf{\Sigma}_{\ell\times
N_{s}}\,{\mathbf{V}^{T}}_{\ell\times N_{s}}$ (9)

It is clear that for the case of SVD decomposition, $\mathbf{U}$ and
$\mathbf{\Sigma}$ together represent a coefficient matrix of dimension
$2N_{T}\times\ell$. Similarly, one can obtain RMF of $\mathbf{H}$ as follows:

$\mathbf{H}_{2N_{T}\times
N_{s}}\approx\mathbf{B}_{2N_{T}\times\ell}\,{\mathbf{Q}^{T}}_{\ell\times
N_{s}}$ (10)

Therefore in generic notation, $H^{\alpha}$ can be written as a combination of
the $\alpha^{\texttt{th}}$ coefficient vector and $\ell$ number of basis
vectors. Generally, these coefficient vectors (computed using SVD or RMF) and
indeed the coefficient matrix is a dense matrix, which implies that all the
$\mathcal{O}(N_{T}\,\ell)$ numbers of elements of coefficient matrix are non-
zero. This implies that all these weights are important for the construction
of the SNR time-series. Therefore, further reduction of the reconstruction
cost is not possible. The only possible way to reduce the reconstruction cost
is to obtain a sparse coefficient matrix in which every row has a, let say,
$k$ number of non-zero weights, where $k<\ell$. For dense coefficient matrix
obtained from the matrix factorization scheme, the reconstruction cost is
$2N_{T}\,\ell\,N_{s}$ floating-point operations. But if we consider a sparse
coefficient matrix, then the reconstruction cost becomes $2N_{T}\,k\,N_{s}$.
Hence, the reconstruction cost can be reduced by order of $\ell/k$. Note that
we need to choose such a value of $k$ for which the error in the
reconstruction SNR time-series is negligible. One can think of making it
sparse by converting some of these non-zero coefficients into zero directly.
One possible way is to make the last $\ell-k$ number of weights equal to zero
directly. Since all $\ell$ basis vectors are important; hence, all the
corresponding weights of the basis have the same importance. Therefore
directly converting the last $\ell-k$ number of weights into zero increases
the waveform reconstruction error, and indeed it has an immense effect on the
reconstruction of the SNR time-series. Alternatively, one can sort out the
coefficients and then convert the last $\ell-k$ coefficients into zero. But
this way of transforming from dense to sparse does not work as the effect
reflects the reconstruction accuracy of the waveforms. Hence the
transformation from the dense to a sparse coefficient matrix by making some
non-zero weights into zero-weights is directly is not suggestive. It is
essential to assign a new value for the rest of the non-zero weights in such a
way that these updated non-zero weights can preserve the length of each
template waveforms (i.e., the rows of $\mathbf{H}$) with high accuracy. This
kind of transformation from dense to sparse coefficients by updating the
weights corresponding to each basis vector can be done using the MP algorithm
[46]. This section investigated the possibility of constructing the sparse
coefficient matrix using the MP algorithm. For this purpose, we need a set of
basis vectors that can be obtained by factorized template matrix
$\mathbf{H}_{2N_{T}\times N_{s}}$ using RMF or SVD. After factorization, we
only choose top-$\ell$ basis vectors for the representation of the waveforms.
Hence, the dense coefficient matrix is of dimension $2N_{T}\times\ell$. if we
want to replace the dense coefficient matrix with a sparse one, that implies
we want to fix a non-zero number of weights in each row vector as $k$, where
$k<\ell$. Then we want to find a new coefficient matrix ($\bar{\mathbf{B}}$,
in case of RMF) which has the same dimension ($2N_{T}\times\ell$) as the
previous dense coefficient matrix, but the number of non-zero elements becomes
$k$ for each row. If $k<\ell$ and the corresponding reconstructed row vectors
of the template matrix can be approximated with high accuracy, we can replace
the dense coefficient matrix with the sparse one. The finding of the new
coefficients using the MP algorithm is an optimization problem that can be
defined as follows:

$\displaystyle\underset{k}{\text{minimize}}$
$\displaystyle\|\bm{x}^{\alpha}\|_{0},\;\alpha=1,\ldots,2N_{T}$ (11) subject
to $\displaystyle\tilde{\mathbf{H}}^{\alpha}=\bm{x}^{\alpha}\,\mathbf{Q}^{T}.$

The Eq.11 represents a sparsity constraint-based optimization problem, where
we want to approximate the row vectors of $\mathbf{H}$ with a specific
predefined number of sparse coefficients $k$. In this optimization problem, we
have assumed that every row vector can be approximated with the same number of
non-zero components $k$. One can also choose a different value of $k$ for
approximating the different row vectors. Note that, for any condition, the
optimal value of $k$ should be computed in a such a way that the approximated
row vector of $\mathbf{H}$ can follow
$\|\mathbf{H}^{\alpha}-\tilde{\mathbf{H}}^{\alpha}\|_{2}<\delta,\;\alpha=1,\ldots,2N_{T}$,
where $\delta$ is the $\texttt{L}_{2}$ error for the approximation of the row
vectors and it should be small.

The above optimization problem can also be redefined as a
$\texttt{L}_{2}$-error ($\epsilon$) constraint-based optimization problem
where the number of sparse coefficients is evaluated based on some predefined
error bound. Hence, it can be considered as follows:

$\displaystyle\epsilon=\underset{\bm{x}^{\alpha}}{\textrm{min}}$
$\displaystyle\|\mathbf{H}^{\alpha}-\bm{x}^{\alpha}\,\mathbf{Q}^{T}\|_{2},$
(12) subject to $\displaystyle\|\bm{x}^{\alpha}\|_{0}\leq
k,\;\alpha=1,\ldots,2N_{T}.$

Alternately, we can specify $\epsilon$, the upper bound on the desired target
error, and try to minimize $k$, the sparsity, subject to this constraint. Both
the optimization problems defined in Eq.12, or the alternative, are non-convex
optimization problems. In fact, it is a NP-hard problem. It can be sub-
optimally solved using the MP algorithm [46], which is an iterative procedure
to obtain the re-weighted non-zero coefficients. It finds each element of a
coefficient vector in the step-by-step iterative process. Given a basis
$\mathbf{Q}$ and a row vector $\mathbf{H}^{\alpha}$, first fixed initial
residual ($\bm{r}$) as $\bm{r}=\mathbf{H}^{\alpha}$ and choose an unselected
basis vector $\mathbf{Q}^{(i)}$ from the set of basis $\mathbf{Q}$ and
recalculate the residual as
$\bm{r}=\mathbf{H}^{\alpha}-\bm{x}\,\mathbf{Q}^{(i)}$. The procedure needs to
be repeated until either $\|\bm{x}\|_{0}\leq k$ or
$\|\bm{r}\|_{2}\leq\epsilon$ is reached.

In general, for the MP algorithm, the Fourier basis, the Haar basis can be
used to obtain the sparse coefficient matrix. However, as we already have a
set of basis vectors obtained from SVD or RMF, we can use these set of basis
vectors directly for the computation of the sparse coefficient vectors. Hence,
this sparse coefficient vectors construction method using MP is easily fitted
with the SVD or RMF based match filtering scheme. Here we have used a set of
basis vectors obtained from RMF as an input to the MP algorithm. For optimal
computational cost reduction using the RMF based matched filtering scheme, we
have combined the RMF method with the MP algorithm. For a fixed $\ell$, we
have demonstrated the proposed RMF algorithm with sparse coefficients obtained
from MP in Alg-2.

Input: Template Matrix {$\mathbf{H}\in\mathbb{R}^{2N_{T}\times
N_{s}}:2N_{T}\leq N_{s}$}, $\ell$, $k$

Output: $\mathbf{Q}_{2N_{T}\times k}$, $\mathbf{\tilde{B}}_{k\times N_{s}}$

1 $\mathbf{\Omega}\in\mathbb{R}^{\ell\times
2N_{T}}\;:\mathbf{\Omega}_{ij}\;\in\mathcal{N}\big{(}0,1\big{)}\;$

2 ${\bar{\mathbf{H}}}_{\ell\times N_{s}}=\mathbf{\Omega}_{\ell\times
2N_{T}}\,\mathbf{H}_{2N_{T}\times N_{s}}$

3 $\mathbf{Q}_{\ell\times N_{s}}=\textrm{qr}(\bar{\mathbf{H}})$

4 for  _$\alpha$ = 1,2,3, $\cdots$, $2N_{T}$_ do

$\bm{\tilde{b}}^{\alpha}=\texttt{MP}\big{(}H^{\alpha},\mathbf{Q}^{T}\big{)}$

// Using Matching Pursuit algorithm

${\mathbf{\tilde{B}}}_{2N_{T}\times\ell}=[\bm{\tilde{b}}^{1}|\bm{\tilde{b}}^{2}|\cdots|\bm{\tilde{b}}^{\alpha}]$

// Row-wise stacking

5

6

7 end for

Algorithm 2 RMF with sparse coefficient Figure 3: The left panel of the figure
shows the distribution of the norm ($\texttt{L}_{2}$) of the reconstructed
waveform using a sparse coefficient vector obtained from the MP algorithm, and
the bases are computed using RMF. The size of the template matrix is
$1162\times 16384$. $\ell=200$ number of basis chosen based on the
$\langle{\delta\rho}/\rho\rangle=2\times 10^{-4}$ and out of $200$ non-zero
coefficients, $k=150$ coefficients are consider as non-zero for the formation
of the sparse coefficient vectors. The figure clearly shows that using sparse
coefficient vectors; we can also preserve the norm of the template waveform
within an error of $10^{-2}$.

We now present the efficiency of the MP algorithm to obtain the sparse
coefficient vectors for each template waveform (i.e., the row vectors of
$\mathbf{H}$). As we used normalized template waveforms (i.e., the
$\texttt{L}_{2}$-norm is unit), it is expected, after reconstruction of the
template waveforms using sparse coefficients obtained from the MP algorithm
and a set of top-$\ell$ basis vectors obtained from RMF can also approximate
the norm of each template waveform with high accuracy nearly to the unit.
Therefore, we have performed a simulation and computed the
$\texttt{L}_{2}$-norm of reconstructed rows of $\mathbf{H}$ using $k$-sparse
coefficient vector obtained using MP algorithm and a set of top-$\ell$ basis
vectors obtained from RMF. Figure-3 shows the histogram of the
$\texttt{L}_{2}$-norm of each row vector after reconstruction using Alg-2. For
this example, we consider a template bank $\mathbf{H}$ containing $N_{T}=581$
templates covering the component mass space: $5\leq m_{1,2}/M_{\odot}\leq 15$
using non-spinning TaylorT4 waveforms. Each waveform was taken to be $8$
seconds long, sampled at $2048$ $\mathrm{Hz}$, thereby setting $N_{s}=16384$.
The required number of basis $\ell=200$ is fixed based on the
$\langle{\delta\rho}/\rho\rangle=2\times 10^{-4}$. The sparsity for each of
the coefficient vectors is considered as $k=150$. It is expected that due to
the conversion from the dense coefficient vector to the sparse coefficient
vector, the $\texttt{L}_{2}$-norm of the reconstructed rows can not be of
exact unit magnitude. But from Figure-3, it is clear that the maximum error in
the norm is $10^{-2}$. Hence the sparsification of the coefficient vector can
preserve the norm of each of the row vectors with high accuracy. Therefore, we
can use the sparse coefficient matrix for the computation of the SNR time-
series. Hence, for this specific example, using the sparse coefficient matrix,
we reduced the reconstruction cost by $25\%$ of the previous reconstruction
cost considering the dense coefficient matrix. This is because we choose a
value of $k=\frac{3}{4}\,\ell$.

Figure 4: The figure shows the reconstructed waveform corresponding to a mass
component $m_{1},m_{2}=36\,,40\,M_{\odot}$. We have chosen a set of $2000$
template waveforms from the range of mass components $30$-$50\,M_{\odot}$. Out
of $2000$ basis vectors, top-$100$ ($\ell=100$) basis vectors have been used
to reconstruct the waveform. However, for the reconstruction of the waveform,
we have not used the coefficient vector obtained from the SVD. We have used
the MP scheme to make the more sparse coefficients. We have used half of the
coefficient vector’s length (i.e. out of $100$, i.e., $k=50$) to construct
sparse coefficients. The reconstruction is done using top-$100$ SVD-basis and
top-$100$ RMF basis (shown in orange and blue colours). The reconstruction
accuracy for both cases is almost the same.

##  V Implementation of RMF for LLOID-type framework

In the GstLAL pipeline, a full template bank is divided into several sub-
banks. The sub-banks division is done based on the threshold on the chirp-mass
and the duration of the template waveforms. The template waveforms from each
sub-bank can further be split into several time slices based on different
sampling frequency rates [10]. The early inspiral signal can be sampled using
a low sampling rate as it enters the detection bandwidth with low frequency.
At the same time, merger and ringdown (up-to-last stable orbit) signals
contain high-frequency components. Thus required sampling rate for the merger
and ringdown will be high as compared to the inspiral part. Also, the time
duration for inspiral, merger, and ringdown is different. Due to the low
frequency, the inspiral waveform duration is longer compared to the merger,
and similarly, the merger waveform has a long duration compared to the
ringdown segment. Adapting this idea, any waveforms from a sub-bank can be
sliced into early inspiral, late inspiral, ringdown, merger, and sampled with
different sampling rates using the down-sample method. For example, The early
and late inspiral parts of waveforms can be sampled with a low sampling rate
e.g. $32$-$256$ $\mathrm{Hz}$. The merger and ringdown parts can be sampled
with $512$-$2048$ $\mathrm{Hz}$ depending on the duration of the template
waveform. Figure-19 demonstrates an example of the time-frequency evolution of
a waveform from a NSBH system and how several sampling frequencies are used to
represent the early inspiral, late inspiral, merger, and ringdown parts of the
waveform. A split bank is defined by combining a specific time slice of all
template waveforms in a given sub-bank. Each sub-bank has several split banks.
Therefore, instead of computing the basis vectors for a sub-bank, it computes
several independent SVD decompositions based on split-banks. Each split bank
decomposition only contributes to calculating a specific fraction of the SNR
time-series corresponding to a particular sampling rate. Further, to obtain
the complete SNR time-series based on a fixed (uniform) sampling rate, one
needs to use the sinc interpolation scheme to upsampled the matched filter
output for those time-slices for which the sampling rate is low. Combining
them with the SNR time-series of the other slices of the high sampling rate
will provide the full approximated SNR time-series based on the top few basis
vectors. The mathematical construction of the LLOID scheme is described in the
Ref. [8]. The approximation of SNR time-series using LLOID strategy plays a
vital role on a region of CBC parameter space centered on BNS masses, for BNS
system, the waveforms are very long, hence required to apply LLOID scheme for
the reduction of the filtering cost. However, the total number of individual
SVD decomposition operations has been increased in that process. Also, an
extra computation cost is incurred due to the use of the sinc interpolation
scheme for the upsampling of the partial SNR time-series. This section
presents RMF schemes that can be useful to obtain basis vectors for a split
bank. Our definition of a ’split bank’ is similar to the description used in a
LLOID framework. However, the main difference in the definition is that we
used the fixed (uniform) sampling rate for all the phases (i.e., early and
late inspiral, merger, and ringdown) of the template waveform. That implies
the sampling rate to generate the merger and ringdown phase is the same as the
inspiral one. Since our main objective is to obtain common basis vectors
combining all the split banks, we need to use the uniform sampling rate for
the generation of each time-slices. We have used the finally obtained common
set of basis vectors for the computation of the SNR time-series. As we used a
fixed sampling rate, downsampling for each phase and the upsampling step via
sinc interpolation in the LLOID scheme are easily excluded in our approach to
obtain the SNR time-series. Thus, our proposed process is simple and
computationally less expensive. For each split bank, the computation of
independent basis vectors via SVD is not required. Hence, we not only reduce
the overall cost of performing several SVD considering all split banks but,
further, the cost of downsampling and upsampling have vanished. Figure-5
demonstrates a split bank for our purpose. The blue box defines the first
split bank, which contains the time-slices of the template waveforms for the
early inspiral part. Similarly, the green and the red boxes represent the
second and $p^{\texttt{th}}$ split bank. It is clear from the figure that the
red box contains the merger time-slice part only. The proposed algorithm is
mainly designed for the handling of a set of long waveforms (i.e., for BNS
system). However, it can also be successfully implemented for a group of short
(i.e., for high mass BBH system) or medium-duration (i.e., for low mass BBH or
NSBH system) waveforms. In Figure-17, we have pictorially shown the
applicability of our various RMF algorithms to the different regimes of the
parameter space for a CBC sources.

Figure 5: The figure demonstrates the construction of the split banks. For a
specific sub-bank, if the waveforms are long (i.e. the case for BNS, NSBH, and
low-mass BBH systems, we can split every waveform into several parts and
stacked them together to construct a split bank. The blue, green, and red
boxes describe such a split bank. In this figure, the red box represents a
split template matrix containing the merger and ringdown part, whereas the
blue and green boxes contain the early inspiral part. Alg-7, and Figure-6
describe the strategy to combine these split banks to obtain the set of basis
vectors for the full one. Alg-8 shows the same scheme but with a fixed error
instead of a predefined rank $\ell$.

###  V.1 Proposed method

For the description of the scheme shown in Alg-6, we assumed that Figure-5
depicts a bank with longer waveforms. That implies that the column dimension
of the constructed template matrix is greater than the row dimension i.e.
$N_{s}\gg 2N_{T}$. Let us consider
$\mathbf{H}=\\{\mathbf{H}_{1}|\mathbf{H}_{2}|\cdots|\mathbf{H}_{i}\\}:i=1,2,\cdots
p$ is a template matrix and $\mathbf{H}_{i}$ represents a split-bank
corresponding to the $i^{\texttt{th}}$ time-slices of the template waveform.
The dimension of each $\mathbf{H}_{i}:i=1,2,\cdots,p$ becomes
$2N_{T}\times\frac{N_{s}}{p}$. For simplicity, we have assumed that each
split-bank has the same dimension. But, one can also choose a different size
for different time-slices. In that case, the column-dimension of the template
matrix corresponding to each individual split-bank will be different. However,
our approach is independent of the split-banks’ column dimension as our final
objective is to obtain the rank-$\ell$ factorization of a template matrix by
accessing the individual split-bank only. Suppose we have only access to each
of these $\mathbf{H}_{i}$ independently that implies these template matrices
are stored in the independent machines. Thus we want to propose a RMF scheme
which can be useful to obtain a $\ell$-rank factorization of the whole
template matrix $\mathbf{H}$ as
$\mathbf{Q}_{2N_{T}\times\ell}\,\mathbf{B}_{\ell\times N_{s}}$ using these
individual template matrices $\mathbf{H}_{i}$. The designed scheme is shown in
Alg-7. Similar scheme using block-wise RMF is also shown in Alg-8. The steps
involved in this algorithm are also discussed in the Figure-6. In this figure,
the stage-$0$ represents a distributed architecture, in which independent
machine has individual fraction of the template waveform i.e.
$\mathbf{H}_{2N_{T}\times\frac{N_{s}}{p}}$. Further each machine can generate
a specific fraction of the RP matrix i.e.
${\Omega_{i}}_{\frac{N_{s}}{p}\times\ell}$ to project a fraction of template
waveforms in a lower (say $b$) dimensional space. Hence, the dimension of each
row vectors of $\bar{\mathbf{H}}_{i}$ becomes $2N_{T}\times\ell$. This step is
shown in the Figure-6, Stage-$1(a)$. The benefit of this approach is that no
need to generate a large RP matrix $\mathbf{\Omega}_{N_{s}\times\ell}$ to
obtain a lower dimensional represent of long template waveforms. Further,
${\mathbf{\Omega}_{i}}_{\frac{N_{s}}{p}\times\ell}$ can be generated
independently in each of machines, hence there are zero communication cost is
involved. In this way, reduction of the space as well as computational cost
for lower-dimensional representation of the waveform can be possible. However,
individual lower-dimensional representation of the split-bank can not be
suitable to obtain the top-$\ell$ basis vectors of the compressed row vector
of the full template matrix. Thus to compute the top-$\ell$ basis vectors from
the lower dimension representation of the whole template matrix, we need to
add all these individual projected components
$\bar{\mathbf{H}}_{i}:i=1,2,\cdots,p$ together. The addition of these
projected components, mathematically shown as
$\bar{\mathbf{H}}=\sum_{i=1}^{p}{\bar{\mathbf{H}}_{i}}$ (Shown in
Stage-$1(b)$. In the next stage (Stage-$2(a)$), the top-$\ell$ basis vectors
can be computed using QR decomposition. The next stage is optional and only
required if we want to compare the singular values of the original data matrix
along with the compressed matrix $\mathbf{B}$. To compute the individual
coefficient matrix $\mathbf{B}_{i}$, we need to re-distribute the
corresponding fraction of $\mathbf{Q}_{i}$ to each of the machines which
already has stored the fraction of the original template matrix. Finally, in
Stage-$3(b)$, we need to stack them together to get the full coefficient
matrix corresponding to the original template matrix. Alternatively,
computation of the coefficient matrix can be done using the MP algorithm; in
that case, no need to generate the coefficient matrix $\mathbf{B}$ using the
projection operation shown by stage-$3$. Figure-7 and Figure-8 described the
efficiency of the proposed scheme shown in Alg-7.

Figure 6: The figure demonstrates the steps involved in the RMF scheme (See
Alg-7 of the appendix) similar to the LLOID-type framework.  Figure 7: The
figure shows the comparison of singular values ($\sigma$) for a template
matrix $\mathbf{H}$ of size $(2N_{T}\times N_{s})\equiv 1324\times 225280$,
normalized by the maximum singular value ($\sigma_{max}$) as obtained from SVD
and RMF algorithm shown in Alg-7. The LLOID-type framework has been
constructed and prescribed RMF (Alg-7) is performed in the target dimensions
$\mathbb{R}^{\ell}$, where $\ell=600$. In this example, template bank was
constructed using a non-spinning waveform model for component mass parameters
($m_{1,2}$) in the range $1.0M_{\odot}\leq m_{1},m_{2}\leq 1.1M_{\odot}$. The
waveforms are generated using the IMRPhenomD waveform model. The lower cut-off
frequency for the waveform generation is chosen as $10$ $\mathrm{Hz}$.
Therefore the duration of the longest waveform of this bank becomes $110$ sec.
Further, the $2048$ $\mathrm{Hz}$ sampling rate is used for the generation of
the waveforms. The waveforms are split into $11$ split-banks, and hence the
duration of each time-slices becomes $10$ sec. Therefore, each template
matrix’s size corresponding to the split bank becomes $1324\times 22528$. As
seen here, the top-$\ell$ singular-values obtained by RMF agree very well with
the spectrum obtained by traditional SVD factorization. Figure 8: The figure
shows the distribution of $\texttt{L}_{2}$-norm of the reconstructed template
waveform using $\ell=100,300,600$ for the template matrix $\mathbf{H}$ used in
Figure-7. It is clear from the figure that the mean of the reconstructed
waveforms’ norm improves with the increasing value of the rank $\ell$.
Similarly, the variance of the norm decreases with a rising value of $\ell$.
This plot shows the efficacy of the proposed RMF scheme shown in Alg-7 as the
obtained basis vectors along with the coefficient matrix can reconstruct the
waveform with very high accuracy.

####  V.1.1 RMF: combining all sub-banks

It is difficult for a huge bank to factorize the full template matrix together
in SVD set-up even in a distributed architecture. Hence, we aim to investigate
the whole template matrix’s factorization by exploring the RMF set-up in a
distributed manner. Firstly consider that the template matrix $\mathbf{H}$ is
too large to store in a single machine, and each machine can store a fixed
number of templates. That implies that each machine can contain a specific
sub-bank. In Figure-9, we assumed that a template matrix corresponding to a
whole template bank is divided into $q$ number of sub-matrices that are
equivalent to the fact that the entire bank is divided into $q$ number of sub-
banks. We can apply RMF or SVD independently to obtain the corresponding basis
vectors for each sub-bank. However, we aim to design a scheme using RMF for
getting a standard set of bases by combining all the sub-banks, which is
beyond the scope of the current framework of the GstLAL pipeline.

Figure 9: The figure demonstrates the division of a template bank into several
sub-banks. The sub-bank division for non-precession waveforms is generally
done based on the chirp-mass threshold ($\mathcal{M}$) of the template
waveforms. This figure assumes that the waveforms in a template bank are
sorted based on the chirp-mass value, and the different color represents
different sub-banks. For each sub-bank, we can generate a template matrix. The
separate colored box of waveforms represents the sub-template matrix
corresponding to each sub-bank. Using SVD or RMF, one can compute a set of the
basis for each sub-banks. However, it is hard to combine all the sub-banks and
obtain a global collection of the basis for the whole bank if the bank size is
large. However, using our proposed RMF algorithm, we can combine them and
obtain a global set of bases, which can be useful for representing the
template waveforms for the whole template bank. The proposed algorithm is
shown in Alg-5 and also described using simple steps in Figure-10.  Figure 10:
The figure demonstrates the steps involved in the RMF algorithm (See Alg-5)
designed to compute a global collection of bases, combining all the sub-banks.

The problem is similar to the problem defined in Alg-7. The only difference is
that in this set-up, we have a fixed number of time-samples $N_{s}$ for each
template waveforms and stored a fraction of the total number of waveforms in a
fixed machine, i.e., the template matrix corresponding to a sub-bank. The
whole scheme can be designed in a distributed set-up where each of these
split-banks $\mathbf{H}_{i}$ can be stored in an independent machine. As the
entire template matrix is distributed over different devices, it is
impractical to perform the SVD or fixed rank RMF scheme onto the whole
template matrix $\mathbf{H}$ directly. Thus, we have prescribed a scheme using
block-wise RMF by which it is possible to get a rank-$\ell$ approximation of
the full template matrix. The details of the algorithm defined in Alg-5,
Alg-6.

The steps involved in this scheme are shown in Alg-5. This scheme is a
distributed version of Alg-3, where the steps are designed as the same spirit
of Alg-3, and also in each step, simple matrix algebraic operations are
adapted. Suppose we have $p$ number of machines and one central node that can
communicate between other machines. In Alg- 5, the step-$1$ shows the
generation of RP matrix $\mathbf{\Omega}$ of size $N_{s}\times\ell$ in each of
the machine. Each machine contains $2{N_{T}}_{b}$ number of rows of the
template matrix $\mathbf{H}$. Therefore, the lower-dimensional representation
using RP of these $2{N_{T}}_{b}$ number of rows are shown in step-$2$ by
constructing $\bar{\mathbf{H}}_{i}$. If $b$ is small, then the column
dimension of these $\bar{\mathbf{H}}_{i}$ matrices are small enough to fit
them into local memory (RAM) each of the machines easily. Next, we should
transfer all this lower-dimensional representation (i.e., in a $\ell$
dimensional space) of the row vectors to the central machine for obtaining the
$\ell$ number of basis vectors. Hence, step-$3$ shows the communication
between all the individual machines with the central machine to collect all
$\ell$ dimensional representation of the rows of $\mathbf{H}$ defining by
$\bar{\mathbf{H}}$. After collecting all lower-dimensional representations of
the row vectors in the central node, one can perform a QR decomposition on
$\bar{\mathbf{H}}$ to obtain $\ell$ number of basis vectors $\mathbf{Q}$. The
step-$4$ illustrates this operation. To compute the coefficient matrix
$\mathbf{B}$ in a distributed set-up, again, we need to pass $2{N_{T}}_{b}$
numbers of rows of $\mathbf{Q}$ to each of the machines, which is shown in
step-$5$. The step-$6$ showed the computation of $\mathbf{B}_{i}$ of dimension
$\ell\times N_{s}$ by projecting the $2{N_{T}}_{b}$ basis vectors onto
$2{N_{T}}_{b}$ rows of $\mathbf{H}$ independently in each machine. But the
final $\mathbf{B}$ should be a projection of $\ell$ basis vectors onto all the
rows of $\mathbf{H}$. Hence after getting $\mathbf{B}_{i}$, one has to send
all this $\mathbf{B}_{i}$ to the central node to compute $\mathbf{B}$, which
is shown in step-$7$. Note that a total of three communications between the
central node and the local machines are required. These are reflected by the
steps $3$, $5$, and $7$, respectively. The communication cost for the
steps-$3$ and $5$ are $2{N_{T}}_{b}\times\ell$ floating-point operations for
each of the local machines. Similarly, for step-$7$, the communication cost is
$\ell\times N_{s}$.

Figure 11: The figure shows the comparison of singular values ($\sigma$) for a
template matrix $\mathbf{H}$ of size $(2N_{T}\times N_{s})\equiv 10888\times
18432$, normalized by the maximum singular value ($\sigma_{max}$) as obtained
from SVD and RMF algorithm shown in Alg-5. The prescribed RMF (Alg-5) scheme
is performed in the target dimensions $\mathbb{R}^{\ell}$, where $\ell=2000$.
In this example, the template bank was constructed using aligned-spin signal
model for component mass parameters ($m_{1,2}$) in the range $5.0M_{\odot}\leq
m_{1},m_{2}\leq 35.0M_{\odot}$. The range of dimensionless spin magnitude was
taken to upto $0.1$. Stochastic template bank placement algorithm is used to
generate the template bank. Further, IMRPhenomPv2 signal model is used for the
generation of the waveforms. The lower cut-off frequency for the waveform
generation is chosen as $30$ $\mathrm{Hz}$. Further, the $2048$ $\mathrm{Hz}$
sampling rate is used for the generation of the waveforms. The waveforms are
split into $11$ sub-banks i.e. $\mathbf{H}=|{\mathbf{H}_{1}}_{1000\times
18432}|{\mathbf{H}_{2}}_{1000\times 18432}|\cdots|{\mathbf{H}_{11}}_{888\times
18432}|$. As seen here, the top-$\ell$ singular-values obtained by RMF agree
very well with the spectrum obtained by traditional SVD factorization. Figure
12: The figure shows the distribution of $\texttt{L}_{2}$-norm of the
reconstructed template waveform using $\ell=500,10^{3},2\times 10^{3}$ for the
template matrix $\mathbf{H}$ used in Figure-11. The mean of the reconstructed
waveforms’ norm improves with the increasing value of the rank $\ell$.
Similarly, the variance of the norm decreases with a rising value of $\ell$.
This plot shows the efficacy of the proposed RMF scheme shown in Alg-5 as the
obtained basis vectors along with the coefficient matrix can reconstruct the
waveform with high accuracy.

###  V.2 Numerical Simulation

Here, we have done a Monte-Carlo simulation to reconstruct the SNR time-series
for a $500$ randomly sampled injection from a specific parameter space. We
consider a template bank $\mathbf{H}$ which contains $N_{T}=581$ templates
covering the component mass space: $5\leq m_{1,2}/M_{\odot}\leq 15$ using non-
spinning $\texttt{TaylorT4}4$ waveforms. Each waveform was taken to be $8$
seconds long, sampled at $2048$ $\mathrm{Hz}$, thereby setting $N_{s}=16384$.
Hence the size of the template matrix is $1162\times 16384$. The required
number of basis $\ell=200$ is fixed based on the
$\langle{\delta\rho}/\rho\rangle=2\times 10^{-4}$. First, we have computed the
top-$200$ basis vectors, and we have used these bases for obtaining a sparse
coefficient for each of the template waveforms using the MP algorithm. For
sparse representation, we have chosen $k=150$ non-zero elements out of $200$,
for each coefficient vector. The simulation result is shown in Figure-13.

Figure 13: It shows the distribution of $\langle{\delta\rho}/\rho\rangle$
obtained from a Monte-Carlo injection study. This experiment is similar to the
numerical study carried out in [1]. Here, we have applied RMF-based
factorization of the template banks. Further, the SNR time-series
reconstruction is done based on the sparse-coefficient obtained using the MP
algorithm. This study is reflecting the performance of Alg-2, in which
$\ell=200$ and $k=150$ have been chosen.

##  VI Discussion and Conclusion

This work demonstrated the practical implementation strategy for computation
of matched filter output between data and template waveforms combining block-
wise RMF and MP algorithm. Block-wise RMF can be used to compute an essential
set of basis vectors, whereas MP can be used for obtaining sparse coefficients
of the template waveforms. Further, this work prescribed some advanced RMF
algorithms (block-wise) for computing a basis for a large bank by combining
all the sub-banks and similarly for long waveforms incorporating all split
banks. In addition, we have designed an alternative and computationally
efficient framework similar to the LLOID based on prescribed algorithms and
explored several regions (BNS, BBH, NSBH) of the parameter space to
demonstrate the efficiency of our framework. The work presented here has
mainly three aspects.

  1. 1.

The designed framework is efficient to obtain the required number of basis
vectors at a fixed average SNR loss.

  2. 2.

Secondly, block-wise RMF algorithms with a combination of MP-based sparse
coefficient construction can also address the issue of reducing the
reconstruction cost of computing matched-filter output. It is notable that the
MP-based sparse coefficient construction scheme is entirely independent of the
RMF algorithm and can be easily adaptable along with other basis finding
methods, e.g., SVD. Hence, the sparse coefficient computation using MP can be
directly applicable to the GstLAL pipeline in the current framework if someone
wants to reduce the cost of reconstruction. The highlights of the work are as
follows:

  3. 3.

Here, we have demonstrated the advanced block-wise RMF Algorithms as an
alternative to the LLOID framework. This scheme can split the waveform into
several parts, similar to the waveform splitting in the LLOID method. However,
the fundamental difference between the two approaches is that in the LLOID
scheme, downsampled is used for the different parts ( i.e., early inspiral,
late inspiral, merger, and ring-down) of the waveform to reduce the size of
the split-bank. In our approach, we have used a uniform sampling frequency
throughout every part of the waveform. Due to the use of downsampling, the
process of obtaining the whole SNR time series in the LLOID framework is
complicated, computationally expensive, as after receiving the SNR from each
split part, upsampling is also required before combining them to get the full-
one. Our scheme is simple, as we use the same sampling frequency, so no need
to upsample the obtained SNR from each split part. The use of downsampling-
based representation can be possible for our proposed setup. In that case, we
need to compute the basis vectors of this split-part separately, similar to
the current framework. Since our objective is to obtain one set of bases,
combining all the split parts. Hence, our designed algorithms are solely made
to get a standard set basis using all the divided parts. Incorporating RMF in
a LLOID framework is a straightforward task in which we need to replace the
SVD method with RMF. We will explore the current LLOID framework of GstLAL
using the RMF setup in the upcoming work.

We believe that the prescribed algorithms can be useful for the fast time-
domain matched filtering calculation in the GstLAL pipeline. However, it is
also necessary to investigate the practical challenges of implementing these
algorithms in the current pipeline. In the current GstLAL pipeline, streamer
multi-media framework has been used and HTC based distributed setup is
developed to run all the steps involved in this pipeline [7]. In recent work,
Gittens et al. [47] showed the adaptability of RMF algorithm in a Apache-SPARK
set-up. SPARK-optimized code with enhanced computation power, making our
algorithms much faster. Such fast computation of basis is unlikely to be
possible for the current SVD-based GstLAL framework, which is computationally
expensive. In the future, we plan to explore further the possibility of
integrating SPARK based RMF with the Gstreamer based matched-filtering setup,
which can reduce the computational complexity of the matched filtering cost
for the CBC search.

###### Acknowledgements.

AR would like to thank Chad Hanna and Sarah Caudill for useful discussions on
the GstLAL pipeline. The authors would like to thank Dilip Krishswamy for
useful discussions and feedback. This work was carried out with the generous
funding available from DST’s grant no. T-150 under the ICPS special call.

##  Appendix A Implementation of Random Projection

Let us consider a template matrix $\mathbf{H}$ has $2N_{T}$ number of rows
$\\{\bm{h}_{1},\cdots,\bm{h}_{i}\,:i=1,2,\cdots,2N_{T}\\}\in\mathbb{R}^{N_{s}}$.
and $\exists$ a linear function
$\Phi:\mathbb{R}^{N_{s}}\rightarrow\mathbb{R}^{\ell}$ which transformed
$\bm{h}_{i}\rightarrow\bm{h}^{\prime}_{i}=\Phi(\bm{h}_{i})$ such a way that
for any pair of row vectors $\bm{h}_{i}$, $\bm{h}_{j}$ follows the following
inequality.

$(1-\epsilon)\,d(\bm{h}_{i},\bm{h}_{j})\leqslant
d(\Phi(\bm{h}_{i})-\Phi(\bm{h}_{j}))\leqslant(1+\epsilon)\,d(\bm{h}_{i},\bm{h}_{j})\,,$
(13)

where $d(\bm{h}_{i},\bm{h}_{j})$ represents the the Euclidean distance between
two row vectors $\bm{h}_{i}$ and $\bm{h}_{j}$ of the template matrix
$\mathbf{H}$. Also,
$\Phi(\bm{h}_{i})=\bm{h}_{i}\,\mathbf{\Omega}\in\mathbb{R}^{\ell}$ represents
the RP of each of the row vectors in $\ell$-dimensional space. The inequality
shown in Eq.13 known as Johnson-Lindenstrauss (JL)lemma. The JL lemma
describes that by projecting any set of vectors from their original
dimensional space to a lower-dimensional space, preserving the pairwise
distances between the vectors can be possible within a distortion factor. The
distortion factor measures the difference between the distance in the original
and projected feature space. That implies, if the distortion factor is
minimal, then the projected space dimension can preserve the pairwise distance
between the vectors optimally. Therefore, the number of projected (reduced)
dimensions is directly proportional to the distortion factor. This
transformation from original dimensional to lower dimension space can be done
easily by taking the linear combination of the rows of the template matrix
($\mathbf{H}_{2N_{T}\times N_{s}}$) with the columns of the RP matrix
$\mathbf{\Omega}_{N_{s}\times\ell}$ i.e.,
$\bar{\mathbf{H}}_{2N_{T}\times\ell}=\sum_{i=1}^{2N_{T}}\sum_{j=1}^{\ell}{\bm{h}_{i}\bm{\Omega}_{j}}=\mathbf{H}\mathbf{\Omega}$
such that
$\mathbb{E}\big{[}\mathbf{\Omega}\,\mathbf{\Omega}^{T}\big{]}=\mathbb{I}$.
This kind of matrix can be generated using several way [34], [36], [35]. One
of the well-known approaches is to use standard Gaussian distribution to
generate the rows of the matrix
$\mathbf{\Omega}:\Omega_{ij}\in\mathcal{N}(0,1)$.

|
---|---
(a) | (b)
Figure 14: (a) The figure shows the distribution of the error (percentile) of
computing the pair-wise distance between a set of whitened template waveforms
in the original dimensional space and projected space (with the dimension half
of the original space) using the two RP operators, Gaussian and sparse. Here,
the error is defined as
$\texttt{Error}=\frac{|d_{ij}-d^{\prime}_{ij}|}{d_{ij}}$, where $d_{ij}$ is
the distance between $i^{\texttt{th}}$ and $j^{\texttt{th}}$ rows of the
template matrix $\mathbf{H}$ in the original space. Similarly,
$d^{\prime}_{ij}$ is the distance between $i^{\texttt{th}}$ and
$j^{\texttt{th}}$ rows of the template matrix $\mathbf{H}$ in the projected
space. The template matrix has been constructed by stacking a set of whitened
template waveforms. (b) The difference between the percentile error
distribution occurred due to two different RP operators. $\Delta Error$
defines the difference between the error obtained from two operators.

We have carried an experiment to demonstrate the inequality described by the
JL lemma. The Figure-14 shows a histogram plot for the error of approximating
the pair-wise distance for a set of whitened template waveform in the original
as well in the projected (lower) dimensional space.

##  Appendix B Subspace Iteration

Alg-3 describes the basic RMF of a template matrix $\mathbf{H}$. The accuracy
of the basic RMF scheme depends on the characteristic of the profile of the
singular values. If the singular values fall sharply, that implies that the
number of important bases is less, whereas the profile is flat, implying that
all bases have equal importance. Hence the error incurred due to the low-rank
representation of the template matrix can be defined as a function of the
singular values as follows:

$\epsilon=\|\mathbf{H}-\mathbf{H}^{(\ell)}\|_{F}=\Big{(}\sum_{i=\ell+1}^{{\textrm{min}\\{2N_{T},N_{s}\\}}}{\sigma_{i}^{2}}\Big{)}^{\frac{1}{2}}$
(14)

From the above relation (Eq.14), it is clear that if the singular spectrum of
$\mathbf{H}$ decays slowly, RMF can have a high reconstruction error. In this
scenario, Halko et. al. [12] proposed to use power iteration scheme (For
details see Alg-$4.4$ of [11] ). The main benefit of applying the power scheme
is that the bases $\mathbf{U}$ and $\mathbf{V}$ remain the same after
operating the power scheme on the template matrix $\mathbf{H}$. That implies,
after applying the power method, the template matrix looks different, but the
computed bases are the same as the original template matrix. Therefore the
vectors of the range of the original template matrix $\mathbf{H}$ are still
written as the linear combination of those basis vectors. However, the act of
the power scheme on $\mathbf{H}$ changes the relative weights of the singular
values and forcibly changes the spectrum of the singular values from slow
decay to rapid decay. Formally, if the singular values of $\mathbf{H}$ are
$\Sigma$, the singular values of
$\big{(}\mathbf{H}\mathbf{H}^{T}\big{)}^{P}\mathbf{H}$ are $\Sigma^{2P+1}$.
Here, we have outlined this relation mathematically. Let us consider, the
template matrix $\mathbf{H}$ decomposes as follows:
$\mathbf{H}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{T}$. Then the representation
of the template matrix after applying the power method can be decomposed as
follows:

$\begin{split}\Big{(}\mathbf{H}\mathbf{H}^{T}\Big{)}^{P}\,\mathbf{H}&=\Big{(}\big{(}\mathbf{U}\,\mathbf{\Sigma}\,\mathbf{V}^{T}\big{)}\,\big{(}\mathbf{U}\,\mathbf{\Sigma}\,\mathbf{V}^{T}\big{)}^{T}\Big{)}^{P}\mathbf{H}\\\
&=\Big{(}\mathbf{U\Sigma(V^{T}V)\Sigma U^{T}}\Big{)}^{P}\mathbf{H}\\\
&=\Big{(}\mathbf{U}\mathbf{\Sigma}^{2}\mathbf{U}^{T}\Big{)}^{P}\mathbf{H}\
\big{[}\text{Since},\,\mathbf{V}^{T}\mathbf{V}=\mathbb{I}_{2N_{T}\times
2N_{T}}\big{]}\\\
&=\Big{(}\mathbf{U}\mathbf{\Sigma}^{2}\mathbf{U}^{T}\Big{)}^{P}\,\Big{(}\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{T}\Big{)}\\\
&=\mathbf{U}\,\mathbf{\Sigma}^{2P+1}\,\mathbf{V}^{T}\ \Big{[}\text{Using
mathematical induction:}\,P=1,2,\cdots\Big{]}\end{split}$ (15)

The above relation shown in Eq.15 holds due to the fact that the matrices
$\mathbf{H}$ and $\big{(}\mathbf{H}\,\mathbf{H}^{T}\big{)}^{P}\,\mathbf{H}$
have the same left and right singular vectors. Since the SVD of
$\Big{(}\mathbf{H}\mathbf{H}^{T}\Big{)}^{P}\mathbf{H}$ is
$\mathbf{U}\mathbf{\Sigma}^{2P+1}\mathbf{V}^{T}$ and the obtained singular
vectors obtained from the SVD is unique [41], hence, the singular vectors for
$\mathbf{H}$ and $\big{(}\mathbf{H}\,\mathbf{H}^{T}\big{)}^{P}$ are the same.
Hence it is clear from the above relation that the singular values spectrum
decays exponentially with the number of power iterations, but the singular
vectors remain same for both the matrices. From Eq.14, it is clear that if the
singular values of the template matrix fall rapidly, then the reconstruction
error is less. However, the error is large for a slowly decaying singular
values spectrum. Hence, in practice, the minimal approximation error for
$\\{\sigma_{\ell+1},\cdots,\sigma_{\textrm{min}\\{2N_{T},N_{s}\\}}\\}$ will be
large. Therefore randomized scheme fails to provide the best $\ell$-dominant
sub-space, and consequently, the scheme fails to provide a best rank-$\ell$
approximation of the template matrix. This problem can frequently occur for
any large template matrix. To overcome this problem, one can sample
$\big{(}\mathbf{H}\mathbf{H}^{T}\big{)}^{P}\mathbf{H}$ instead of sampling
$\mathbf{H}$. One can show that
$\big{(}\mathbf{H}\mathbf{H}^{T}\big{)}^{P}\mathbf{H}$ has the same left and
right singular vectors as $\mathbf{H}$. However, the singular values of
$\big{(}\mathbf{H}\mathbf{H}^{T}\big{)}^{P}\mathbf{H}$ are the
$(2P+1)^{\texttt{th}}$ power of the singular values of $\mathbf{H}$ , i.e.,
$\sigma{\big{(}\mathbf{H}\mathbf{H}^{T}\big{)}^{P}\mathbf{H}}=\sigma_{i}^{2P+1}\big{(}\mathbf{H}\big{)}$
[ As shown in Eq.15]. This relation shows, even though the singular values
spectrum of $\mathbf{H}$ falls slowly, then also the singular values of
$\big{(}\mathbf{H}\mathbf{H}^{T}\big{)}^{P}\mathbf{H}\big{)}$ fall sharply.
Hence it is always better to operate the RP on
$\big{(}\mathbf{H}\mathbf{H}^{T}\big{)}^{P}\mathbf{H}\big{)}$ instead of
$\mathbf{H}$. Finally, one has to compute the orthogonal basis vectors
$\mathbf{Q}$ of
$\big{(}\mathbf{H}\mathbf{H}^{T}\big{)}^{P}\mathbf{H}\mathbf{\Omega}$.
Step-$2$ $\&$ $3$ of Alg-(3) can be replaced as follows:

  * •

$\bar{\mathbf{H}}=\big{(}\mathbf{H}\mathbf{H}^{T}\big{)}^{P}\mathbf{H}\mathbf{\Omega}$

  * •

$\mathbf{Q}=qr\big{(}\bar{\mathbf{H}}\big{)}=qr\big{(}\big{(}\mathbf{H}\mathbf{H}^{T}\big{)}^{P}\mathbf{H}\mathbf{\Omega}\big{)}$

The optimal value of $P$ entirely depends on the low-rank structure of the
template matrix. For a large template matrix, it is better to apply the power
scheme to reduce the approximation error—however, the computational complexity
increases due to the implementation of the power scheme. Hence, using the
power scheme for a large template matrix in a distributed set-up is always
recommended. Also, the power scheme can be easily adaptable in the distributed
system architecture. In general, a small value of $P\,(\,=1,2,3)$ can also
provide sufficient decay on the singular value spectrum; therefore, highly
accurate results can be attained using more computational resources. The
Figure-15 demonstrates an example in which singular values spectrum using SVD,
RMF (Alg-3) and RMF with power iteration. The last few singular values
obtained using RMF are not matched exactly with the original one. After
applying the same algorithm with power iteration, we obtained the same
top-$600$ singular values obtained from SVD.

Figure 15: The figure shows a singular value spectrum of a template matrix
$\mathbf{H}\equiv 1324\times 225280$. The singular values are computed using
SVD, RMF with fixed rank-$\ell$, and RMF with power scheme. Here $\ell=600$,
and $P=3$ is chosen.

##  Appendix C Details of the random matrix factorization (RMF)

###  C.1 Fixed-rank RMF (Category-I)

In this section, we have described two categories of fixed-rank RMF schemes.
The fixed-rank RMF (Category-I) describes by Alg-3, whereas Alg-4 represents
the fixed-rank RMF (Category-II). The only difference between these two
categories is the RP on the template matrix. In the first category, the row
vectors of the template matrix are projected in the lower-dimensional space.
However, in the latter category, the column vectors of the template matrix are
cast in the lower dimensional space. That implies that the obtained set of
basis vectors from the first category span the column space, and the obtained
basis vectors from the second category span the row-space considering all the
template waveforms represent a vector space. The mathematical concept is
similar to the SVD decomposition, where $\mathbf{U}$ and $\mathbf{V}$ are the
set of basis vectors of the row space and column space, respectively. The
pictorial description of Alg-3 is shown in Figure-1.

Input: Template Matrix {$\mathbf{H}_{2N_{T}\times
N_{s}}\in\mathbb{R}^{N_{s}}:2N_{T}\geq N_{s}$}.

Output: $\mathbf{Q}_{2N_{T}\times\ell},\mathbf{B}_{\ell\times N_{s}}$

$\mathbf{\Omega}_{N_{s}\times\ell}:\mathbf{\Omega}_{ij}\in\mathcal{N}(0,1)\\}\;\&$

// Generate a random matrix

$\bar{\mathbf{H}}=\mathbf{H\Omega}$

// Compression of the column space

$\mathbf{Q}=qr(\bar{\mathbf{H}})$

//  Set of basis vectors of the projected space

$\mathbf{B}=\mathbf{Q}^{T}\mathbf{H}$

// Surrogate matrix s.t.
${\|{\mathbf{H}}\|_{F}}^{2}-{\|{\bf{B}}\|_{F}}^{2}\leq\delta$

Algorithm 3 RMF with a fixed rank ($\ell$)

###  C.2 Fixed-rank RMF (Category-II)

Input: Template Matrix {$\mathbf{H}\in\mathbb{R}^{2N_{T}\times
N_{s}}:2N_{T}\leq N_{s}$}, $\ell$

Output: $\mathbf{Q}_{2N_{T}\times\ell}$, $\mathbf{B}_{\ell\times N_{s}}$

1
$\mathbf{\Omega}_{i}\in\mathbb{R}^{N_{s}\times\ell}\;:\mathbf{\Omega}_{ij}\;\in\mathcal{N}(0,1)\;$

2 for  _P = 1, 2_ do

3 ${\bar{\mathbf{H}}}_{N_{s}\times\ell}={\mathbf{H}^{T}}_{N_{s}\times
2N_{T}}\,\big{(}\mathbf{H}_{2N_{T}\times
N_{s}}\,{\mathbf{\Omega}_{i}}_{N_{s}\times\ell}\big{)}$

4 ${\mathbf{\Omega_{i}}}_{N_{s}\times\ell}=\textrm{qr}(\bar{\mathbf{H}})$

5

6 end for

$\mathbf{Q_{i}}_{N_{s}\times\ell}=\mathbf{\Omega_{i}}$

// Consider top-$b$ basis vectors out of $\ell$ basis vectors.

7 $\mathbf{B_{i}}_{2N_{T}\times\ell}=\mathbf{H}\,\mathbf{Q_{i}}$

8

Algorithm 4 RMF ($\mathbf{H}\approx\mathbf{B}\,\mathbf{Q}^{T}$) with power
scheme

###  C.3 Mathematical preliminaries of block-wise RMF

Mathematical explanation of Null-space projection:

In this sub-section, we described the mathematical understanding of the null-
space projection of the row vectors of the template matrix $\mathbf{H}$ which
is an essential step for the RMF algorithm with a fixed error. Let us
consider, a template matrix $\mathbf{H}_{2N_{T}\times N_{s}}$ decomposed as
$\mathbf{H}\approx\mathbf{Q}_{2N_{T}\times\textrm{min}{(2N_{T},N_{s})}}\,\mathbf{B}_{\textrm{min}(2N_{T},N_{s})\times
N_{s}}$. Where
$\mathbf{Q}=\\{\bm{q}_{i}\\}_{i=1}^{\textrm{min}\big{(}2N_{T},N_{s}\big{)}}$
as a set of basis and
$\mathbf{B}=\\{\bm{b}_{i}\\}_{i=1}^{\textrm{min}\big{(}2N_{T},N_{s}\big{)}}$
as a set of coefficient vectors.

  1. 1.

The template matrix $\mathbf{H}_{2N_{T}\times N_{s}}$ can be decomposed as a
sum of rank-one matrices as follows:

$\mathbf{H}=\mathbf{H}_{1}+\mathbf{H}_{2}+\cdots+\mathbf{H}_{\textrm{min}\big{(}2N_{T},N_{s}\big{)}}\,,$
(16)

where
$\mathbf{H}_{i}\Big{(}=\bm{q}_{i}\bm{b}_{i},i=1,2,\cdots,\textrm{min}(2N_{T},N_{s})\Big{)}$
represents the rank-one matrices.

  2. 2.

Using the concept of rank-one matrix decomposition, one can similarly
decompose a matrix $\mathbf{H}$ into a sum of rank-$b$ matrix as follows:

$\mathbf{H}=\mathbf{Q}_{1}\mathbf{B}_{1}+\mathbf{Q}_{2}\mathbf{B}_{2}+\cdots+\mathbf{Q}_{\textrm{min}\big{(}2N_{T},N_{s}\big{)}}\mathbf{B}_{\textrm{min}\big{(}2N_{T},N_{s}\big{)}}\,,$
(17)

where the dimension of each basis matrix
$\mathbf{Q}_{i}:i=1,2,\cdots,\textrm{min}(2N_{T},N_{s})$ is
$b\times\textrm{min}(2N_{T},N_{s})$, and dimension for the coefficient matrix
$B_{i}:i=1,2,\cdots,\textrm{min}(2N_{T},N_{s})$ is
$\textrm{min}(2N_{T},N_{s})\times N_{s}$. For simplicity, we considered that
all matrices $\mathbf{Q}_{i}$ contains same number of basis i.e., $b$.

  3. 3.

Let us define a residual matrix $\mathbf{R}$ as follows
$\mathbf{R}=\mathbf{H}-\mathbf{Q}\mathbf{B}$.
After first iteration, it can be computed as:
$\mathbf{R}^{(1)}=\mathbf{H}-\mathbf{Q}_{1}\mathbf{B}_{1}$. If the
$\mathbf{Q}_{1}\mathbf{B}_{1}$ represents the $b$-rank approximation of the
template matrix $\mathbf{H}$, then the residual matrix (after the first
iteration) implies that it only contains the contribution of the
$\Big{(}\textrm{min}(2N_{T},N_{s})-b\Big{)}$ number of the basis vectors as
the contribution for the top-$b$ basis vectors already subtracted out from it.
In the second iteration, we want to compute the second set of important basis
vectors $b+1,b+2,\cdots,b+b$. Therefore, to compute this set, we can use the
residual matrix $\mathbf{R}^{(1)}$ as it contains a contribution of the rest
of the basis vectors. Thus, we need to apply the QR decomposition of
$b$-dimensional representation of the residual matrix $\mathbf{R}^{(1)}$ to
obtain the set of basis $\mathbf{Q}_{2}$ i.e.
$\mathbf{Q}_{2}=\texttt{qr}\big{(}{\mathbf{R}^{(1)}\mathbf{\Omega_{2}}}\big{)}$.
Further, the coefficient matrix $\mathbf{B}_{2}$ can be computed by projecting
all the basis vectors on $\mathbf{R}^{(1)}$, i.e.
$\mathbf{B}_{2}=\mathbf{Q}_{2}^{T}\mathbf{R}^{(1)}$. This process is called
the null-space projection of the basis vectors. In every iteration, we
subtract the energy (contribution) of a specific set of bases, and then the
next set of importance bases is computed from the residual matrix. In this
way, we reach the null-space after a $i^{\texttt{th}}$ iteration, and the
total number of basis vectors up to that iteration is called the rank of the
matrix, e.g., if in each iteration, we get top-$b$ basis vectors and the rank
of the matrix is $\ell$, then after $i^{\texttt{th}}$ iteration we should
obtain top $\ell$ basis vectors, that implies $\ell=ib$. The rest of
$(\ell-\textrm{min}(2N_{T},N_{s})$ basis vectors are less important and after
$i^{\texttt{th}}$ iteration the energy of the matrix is also minimal,
therefore, these $(\ell-\textrm{min}(2N_{T},N_{s})$ are considered the basis
vectors for the null-space. We can compute the residual matrix after the
second iteration as follows:
$\mathbf{R}^{(2)}=\mathbf{H}-\mathbf{Q}_{2}\mathbf{B}_{2}=\mathbf{H}-\mathbf{Q}_{2}\Big{(}\mathbf{Q}_{2}^{T}\mathbf{R}^{(1)}\Big{)}$.
This implies, after second iteration, we already got the $2b$ of basis vectors
combining $\big{[}\mathbf{Q}_{1},\mathbf{Q}_{2}\big{]}$, and the residual
matrix now contains the energy of the space span by the maximum
$(\textrm{min}(2N_{T},N_{s})-2b)$ basis vectors. Similarly, we can generalize
the calculation of the residual part after $i^{\texttt{th}}$ iteration as
follows:
$\mathbf{R}^{(i)}=\mathbf{H}-\mathbf{Q}_{i}\mathbf{B}_{i}=\mathbf{H}-\mathbf{Q}_{i}\Big{(}\mathbf{Q}_{i}^{T}\mathbf{R}^{(i-1)}\Big{)}$.

We used the null-space projection for the computation of the basis and
coefficient matrices $\\{\mathbf{Q_{i}}\\}_{i=1}^{r}$ and
$\\{\mathbf{B_{i}}\\}_{i=1}^{r}$ in each iteration in the block-wise RMF
algorithm. We first initiate the algorithm by setting
$\mathbf{H}^{(0)}=\mathbf{H}$. Then we follow following step for each
iteration $i=1,2,\cdots\ell$.

  1. 1.

Computation of the basis vectors from a set of vectors from a lower-
dimensional row-space.
$\mathbf{Q}_{i}=\texttt{orth}\big{(}\mathbf{H}^{(i-1)}\mathbf{\Omega}_{i}\big{)}$
In the first iteration, the row vectors of the data matrix are projected to a
$b$-dimensional space using RP matrix ${\mathbf{\Omega}_{i}}_{b\times N_{s}}$.
However, from the first iteration onwards, we computed the residual energy (in
terms of Frobenius norm) of the data matrix, and from the second iteration
onwards, the set of basis, e.g., $\\{\mathbf{Q}_{2},\mathbf{Q}_{3},\cdots$ has
been computed using the null-space projection of the row-space obtained from
the previous iteration.

$\displaystyle\mathbf{Q}_{2}$
$\displaystyle=\texttt{qr}\big{(}\mathbf{R}^{(1)}\mathbf{\Omega}_{2}\big{)}$
$\displaystyle=\texttt{qr}\big{(}\mathbf{H}\mathbf{\Omega}_{2}-\mathbf{Q}_{1}\mathbf{B}_{1}\mathbf{\Omega}_{2}\big{)}$

  2. 2.

In each iteration, it is required to obtain the coefficient matrix to measure
the corresponding error.
$\mathbf{B}_{i}={\mathbf{Q}_{i}}^{T}\mathbf{H}^{(i-1)}$. From the second
iteration onwards, $\mathbf{H}^{(i-1)}$ represents the null-space projection.
Therefore, the coefficient matrix is computed using the residual matrix of the
previous iteration. That implies that the coefficient vectors are computed by
projecting the residual matrix onto the block of basis vectors. However, using
the following steps, one can show mathematically that this is exactly
equivalent to the projection of the orthogonal basis onto the original data
matrix if the computed basis vectors in the current stage are orthogonal to
the previous stage, i.e., $\mathbf{Q}_{i-1}^{T}\,\mathbf{Q}_{i}=0$.

$\displaystyle\mathbf{B}_{2}$
$\displaystyle=\mathbf{Q}_{2}^{T}\mathbf{R}^{(1)}$
$\displaystyle=\mathbf{Q}_{2}^{T}\mathbf{H}-\mathbf{Q}_{2}^{T}\mathbf{Q}_{1}\mathbf{B}_{1}$
$\displaystyle=\mathbf{Q}_{2}^{T}\mathbf{H}\,\big{[}\mathbf{Q}_{j}^{T}\mathbf{Q}_{i}=0,i\neq
j\big{]}$

  3. 3.

In each iteration, it is required to trace the residual energy of the data
matrix, and it can be done by subtracting out the $b$-rank approximation of
the matrix from the data matrix in each iteration.
$\mathbf{H}^{(i)}=\mathbf{H}^{(i-1)}-\mathbf{Q}_{i}\mathbf{B}_{i}$

$\mathbf{H}^{(i)}$ holds precisely the residual remaining after
$i^{\texttt{th}}$ step. This means that incorporating adaptive rank
determining can be possible by simply computing the residual energy in each
iteration i.e. $\|\mathbf{H}^{(i)}\|_{F}\leq\epsilon$. The pictorial
representation of this block-wise RMF has been demonstrated in Fig.16.

Figure 16: This figure demonstrates an iterative method to compute basis
vectors at a predefined average SNR loss. The first step is to project all the
template waveforms in a small ($b$) dimensional space and then compute the
top-$b$ basis vectors. If these set of basis vectors are sufficient to
approximate the $\langle{\delta\rho}/\rho\rangle$ with a predefined error
$\epsilon$, then no need to process further. The process has to be continued
unless the error reaches the desired accuracy. The corresponding set of basis
vectors are the final optimal set of bases.

###  C.4 Type-I RMF (row-wise)

Alg-5 shows the computation of a standard set of basis vectors by combining
all the sub-banks of a specific template bank. We assumed that the number of
essential basis vectors by combining all the sub-banks is $\ell$. Hence, the
algorithm provides top-$\ell$ basis vectors and their corresponding
coefficients. Further, we assumed that the whole template matrix is not be
stored in a single machine due to the size of the template matrix. Hence, the
entire template matrix is divided into a sub-template matrix. These matrices
are placed over a set of local machines. These sub-template matrices are
similar to the sub-banks defined in the GstLAL pipeline. Currently, for the
non-precession waveforms, the sub-bank division in the GstLAL has been done
based on chirp-mass ($\mathcal{M}$) and duration of the waveforms. We have
also followed the same way of obtaining sub-banks. In GstLAL, each sub-banks
are treated separately, and computed basis vectors span the waveforms from
those sub-banks only. Thus the basis vectors obtained from each sub-banks are
independent of each other. Here we present a scheme ( described by Alg-5) to
get a basis combining all the sub-banks. That implies that the obtained basis
vectors are common to all the sub-banks.

In this section, we have described Alg-5 and estimates the required storage
memory in each step. Let us consider, a template bank has $N_{T}$ number of
template waveform and each of the waveform has $N_{s}$ time samples. Hence,
the size of the template matrix is $H_{2N_{T}\times N_{s}}$. Also, consider
that we have $n_{b}$ number of individual machines and each machine has
$16\texttt{GB}$ RAM. Each machine can contains $N_{T_{b}}$ number of
templates. Hence, the total number of templates in a template bank is
$N_{T}=n_{b}\,N_{T_{b}}$. Let $\ell$ is the required number of important bases
corresponding to a specific $\langle{\delta\rho}/\rho\rangle$. In RMF, the
template waveforms are projected in a lower-dimensional space using a RP
operator. In Alg-5, the compressed version (or the projected version) of the
waveforms for each sub-banks is formulated using a RP matrix
$\mathbf{\Omega}$; the step is shown in the step-$2$. All the sub-banks are
distributed over different machines; hence, we have two options to perform the
projection operation. The first option is that we can provide the same
$\mathbf{\Omega}$ to all machines and then complete the projection operation
on each sub-bank separately. The alternative option is that we can generate
the same $\Omega$ for all machines. The first option is computationally
expensive, as there is a communication cost is involved for the passing of the
same $\mathbf{\Omega}$ to all machines. However, the second option does not
require any communication cost. Hence it is optimal. But the problem is that,
in general, each machine can generate different $\mathbf{\Omega}$, which
potentially affects the projection procedure as the waveforms from all the
sub-banks should be projected using the same RP vectors to the lower-
dimensional subspace. It can be resolved if we use a fixed seed for the
generation of all random vectors of $\mathbf{\Omega}$ for all the machines.
The step $3,4$ of Alg-5 needs to be done in a central node connected with all
the other local machines. This node can be called a ”master node.” Step- $3$
shows the collection of all lower-dimensional represents of all template
waveforms from all the sub-banks. In step-$4$, the basis is computed from that
$\ell$ dimensional sub-space. To estimate each step’s required storage memory,
let us consider a template matrix containing $10^{5}$ number of template
waveform, and each waveform has $10^{6}$ of sample points. The storing of such
a large waveform required $160\texttt{GB}$ disk space. Let us consider; each
machine has $16\texttt{GB}$ RAM memory space. Therefore, each machine can
contain $10^{2}$ number of template waveform that means the size of the each
sub-matrix $\mathbf{H_{i}}$ is $2N_{T_{b}}\times N_{s}\equiv 200\times
10^{6}$. The required memory to store such a sub-bank is $1.6\texttt{GB}$. Let
us consider $\ell=10^{4}$, therefore to generate the RP matrix required memory
space is $80\texttt{GB}$. Since the machine does not have that much memory, it
is impossible to generate the whole RP matrix in one machine. We can generate
a block of the RP matrix at a time. Therefore $\mathbf{\Omega}$ defined in
step-$1$ can be generated as follows:
$\mathbf{\Omega}=\Big{[}{\mathbf{\Omega}_{1}}_{2N_{T_{b}}\times
b}|{\mathbf{\Omega}_{2}}_{2N_{T_{b}}\times
b}|\cdots|{\mathbf{\Omega}_{j}}_{2N_{T_{b}}\times b}\Big{]}$. If we consider
$b=10^{3}$ and $j=10$, then $\mathbf{\Omega}_{j}$ can be of size $10^{6}\times
10^{3}$. The required memory to generate each of RP block is $1$ GB. Hence,
the step-$2$ can be computed using RP block as follows:
$\bar{\mathbf{H}}_{i}=\mathbf{H}_{i}\,\mathbf{\Omega}_{1}+\mathbf{H}_{i}\,\mathbf{\Omega}_{2}+\cdots+\mathbf{H}_{i}\,\mathbf{\Omega}_{j}$.
To store each of these sub-matrices $\mathbf{H}_{i}\,\mathbf{\Omega}_{j}$
required memory is $16\texttt{MB}$. Thus, the required memory to store all
these sub-matrices and the final output matrix $\bar{\mathbf{H}}_{i}$ shown in
step-$2$ is $320\,(=16\times 10+160)\texttt{MB}$. The required memory to store
randomly projected template matrix $\bar{\mathbf{H}}$ shown in step-$3$ is
$16\texttt{GB}$. Note that this matrix is stored in the master node. The
master node should have a large memory space compared to the local distributed
nodes (machines). There is a communication cost also involved in this step.
$2\times 10^{6}$ number of floating-point operations are required to send the
matrix $\bar{\mathbf{H}}_{i}$ from the individual machine to the master node.
Overall, the master node required $32\texttt{GB}$ of memory space to save
$\bar{\mathbf{H}}$ and $\mathbf{Q}$ in the memory. Again there is also a
$10^{6}$ number of floating-point operations are needed to send back again the
corresponding bases $\mathbf{Q}_{i}$ to the individual machine from the master
node. To perform the QR decomposition in step-$4$ required $2\times 10^{13}$
floating-point operations. The required memory to store $\mathbf{Q}_{i}$ in
each machine is $16\texttt{MB}$. Also, to store $B_{i}$ it is required
$10\texttt{GB}$. Another communication involves from local machines to the
central node to compute the computation of the coefficient matrix $\mathbf{B}$
for which the communication cost is $10^{10}$ number of floating-point
operations. Alg-6 describes the block-RMF version of Alg-5 at a fixed error.
However, if we combine RMF with MP, then the computation of the coefficient
matrix $\mathbf{B}$ is not required, and in that case, we should not count the
estimated cost for the construction of the coefficient matrix.

Input: Sub-bank {$(\mathbf{H}_{i})_{2{N_{T}}_{b}\times N_{s}}$}

Output: $\mathbf{Q}_{2N_{T}\times\ell}$, $\mathbf{B}_{\ell\times N_{s}}$

1
$\mathbf{\Omega}\in\mathbb{R}^{N_{s}\times\ell}:\Omega_{ij}\in\mathcal{N}(0,1)$

2
$({\bar{\mathbf{H}}_{i}})_{2{N_{T}}_{b}\times\ell}=\mathbf{H}_{i}\,\mathbf{\Omega}$

$(\bar{\mathbf{H}})_{2N_{T}\times\ell}=[\mathbf{\bar{H}}_{1}|\mathbf{\bar{H}}_{2}|\cdots|\mathbf{\bar{H}}_{i}]$

// Row-wise stacking

3
$\mathbf{Q}_{2N_{T}\times\ell},\mathbf{R}_{\ell\times\ell}=\text{qr}(\bar{\mathbf{H}})$

4
$(\mathbf{Q}_{i})_{2{N_{T}}_{b}\times\ell}=\mathbf{Q}[(i-1){N_{T}}_{b}\,:i{N_{T}}_{b},:];i=1,2,\cdots,n_{b}$

5 ${\mathbf{B}_{i}}_{{}_{\ell\times
N_{s}}}={\mathbf{Q}_{i}}^{T}\,\mathbf{H}_{i}$

6 $\mathbf{B}_{\ell\times N_{s}}=\sum_{i=1}^{n_{b}}{\mathbf{B}_{i}}$

7

Algorithm 5 Type-I RMF (row-wise)

Input: Sub-bank {$\big{(}\mathbf{H}_{i}\big{)}_{2{{N_{T}}_{b}}\times N_{s}}$}

Output: $\mathbf{Q}_{2N_{T}\times\ell}$, $\mathbf{B}_{\ell\times N_{s}}$

1 for  _i = 1,2,3, $\cdots$_ do

2 $\mathbf{\Omega}_{i}\in\mathbb{R}^{N_{s}\times
b}:\Omega_{ij}\in\mathcal{N}\big{(}0,1\big{)}$

3 $({\bar{\mathbf{H}}_{i}})_{2{{N_{T}}_{b}}\times
b}=\mathbf{H}_{i}\,\mathbf{\Omega}_{i}$

$(\bar{\mathbf{H}})_{2N_{T}\times
b}=\big{[}\mathbf{\bar{H}}_{1}|\mathbf{\bar{H}}_{2}|\cdots|\mathbf{\bar{H}}_{i}\big{]}$

// Row-wise stacking

4 $\mathbf{Q}_{2N_{T}\times b},\mathbf{R}_{b\times
b}=\text{qr}(\bar{\mathbf{H}})$

5
$\mathbf{Q}_{i}=\textrm{qr}\,\big{(}\mathbf{Q}_{i}-\sum_{j=1}^{i-1}{\mathbf{Q}_{j}{\mathbf{Q}_{j}}^{T}\mathbf{Q}_{i}}\big{)}$

6 $\big{(}\mathbf{Q}_{i}\big{)}_{2{{N_{T}}_{b}}\times
b}=\mathbf{Q}\big{[}(i-1){N_{T}}_{b}\,:i{N_{T}}_{b},:\big{]};i=1,2,\cdots,n_{b}$

7 ${\mathbf{B}_{i}}_{{b\times N_{s}}}={\mathbf{Q}_{i}}^{T}\,\mathbf{H}_{i}$

8 $\mathbf{B}_{b\times N_{s}}=\sum_{i=1}^{n_{b}}{\mathbf{B}_{i}}$

9 $\mathbf{H}=\mathbf{H}-\mathbf{Q}_{i}\,\mathbf{B}_{i}$

10 if  _$\|\mathbf{H}\|_{F}\geq\langle\frac{\delta\rho}{\rho}\rangle$_ then

11
$\mathbf{Q}=\big{[}\mathbf{Q}_{1}|\mathbf{Q}_{2}|\mathbf{Q}_{3}|\cdots|\mathbf{Q}_{i}\big{]}$

12
$\mathbf{B}=\big{[}{\mathbf{B}_{1}}^{T}|{\mathbf{B}_{2}}^{T}|\cdots|{\mathbf{B}_{i}}^{T}\big{]}$

13 end if

14

15 end for

Algorithm 6 Type-I RMF (row-wise) at a fixed error)

###  C.5 Type-II RMF Algorithm (column-wise division)

The Alg-7 shows the scheme of obtaining basis at a fixed rank-$\ell$ over a
distributed architecture framework. In this case, we assumed that the number
of rows of a template matrix is less than the number of columns. Hence, the
columns of the template waveforms are distributed over a set of local
machines.

The prescribed framework is similar to the split-bank framework of the GstLAL
pipeline. In a split-bank (or time-slices) framework, the template waveforms
are sliced over time for a specific sub-bank. In GstLAL, split banks are made
by considering different sampling frequencies for the different parts of the
waveform. Individual SVD has performed to a specific region to obtain the
basis vectors. This approach helps handle a long waveform. We split a template
into several regions for our purpose, but the sampling rate for generating
these parts is the same. Also, we have combined all the parts to obtain a set
of bases by combining all the split-banks. Here we assumed that each local
machine could store only a fraction of the sampling points.

Input: Split bank {$\mathbf{H}_{2N_{T}\times
N_{s}}=\big{[}\mathbf{H}_{1}|\mathbf{H}_{2}|\cdots|\mathbf{H}_{i}\big{]}:\mathbf{H}_{i}\in\mathbb{R}^{2N_{T}\times(N_{s}/p)}:2N_{T}\leq(N_{s}/p)$}.

Output: $\mathbf{Q}_{2N_{T}\times\ell}$, $\mathbf{B}_{\ell\times N_{s}}$

1
$\mathbf{\Omega}\in\mathbb{R}^{(N_{s})\times\ell}:\Omega_{ij}\in\mathcal{N}(0,1)$

2
$\mathbf{\Omega}_{i}=\mathbf{\Omega}[(i-1)(N_{s}/p):i(N_{s}/p),:],i=1,2,\cdots,p$

3
$({\bar{\mathbf{H}}_{i}})_{2{N_{T}}\times\ell}=\mathbf{H}_{i}\,\mathbf{\Omega}_{i}$

4
$(\bar{\mathbf{H}})_{2N_{T}\times\ell}=\sum_{i=1}^{n_{b}}{\bar{\mathbf{H}}_{i}}$

5
$\mathbf{Q}_{2N_{T}\times\ell},\mathbf{R}_{\ell\times\ell}=\text{qr}(\bar{\mathbf{H}})$

6 ${\mathbf{B}_{i}}_{{}_{\ell\times(N_{s}/p)}}=\mathbf{Q}^{T}\,\mathbf{H}_{i}$

$\mathbf{B}_{\ell\times
N_{s}}=[\mathbf{B}_{1}|\mathbf{B}_{2}|\cdots|\mathbf{B}_{i}]$

// Column-wise stacking

Algorithm 7 Type-II RMF (column-wise)

Input: Split bank {$\mathbf{H}_{2N_{T}\times
N_{s}}=\big{[}\mathbf{H}_{1}|\mathbf{H}_{2}|\cdots|\mathbf{H}_{i}\big{]}:\mathbf{H}_{i}\in\mathbb{R}^{2N_{T}\times\big{(}N_{s}/p\big{)}}:2N_{T}\leq\big{(}N_{s}/p\big{)}$}.

Output: $\mathbf{Q}_{2N_{T}\times\ell}$, $\mathbf{B}_{\ell\times N_{s}}$

1 for  _i = 1,2,3, $\cdots$_ do

2 $\mathbf{\Omega}_{i}\in\mathbb{R}^{(N_{s}\times
b)}:\Omega_{ij}\in\mathcal{N}\big{(}0,1\big{)}$

3
${\mathbf{\Omega}_{i}}^{(i)}=\mathbf{\Omega}\big{[}(i-1)(N_{s}/p):i(N_{s}/p),:\big{]},i=1,2,\cdots,p$

4 $\big{(}{\bar{\mathbf{H}}_{i}}\big{)}_{2N_{T}\times
b}=\mathbf{H}_{i}\,{\mathbf{\Omega}_{i}}^{(i)}$

5 $\big{(}\bar{\mathbf{H}}\big{)}_{2N_{T}\times
b}=\sum_{i=1}^{n_{b}}{\mathbf{H}_{i}}$

6 $\mathbf{Q}_{2N_{T}\times
b},\mathbf{R}_{b\times\ell}=\textrm{qr}(\bar{\mathbf{H}})$

7 ${\mathbf{B}_{i}}_{{b\times(N_{s}/p)}}=\mathbf{Q}^{T}\,\mathbf{H}_{i}$

$\mathbf{B}_{b\times
N_{s}}=\big{[}\mathbf{B}_{1}|\mathbf{B}_{2}|\cdots|\mathbf{B}_{i}\big{]}$

// Column-wise stacking

8 $\mathbf{H}=\mathbf{H}-\mathbf{Q}_{i}\,\mathbf{B}_{i}$

9 if  _$\|\mathbf{H}\|_{F}\geq\langle\frac{\delta\rho}{\rho}\rangle$_ then

10
$\mathbf{Q}=\big{[}\mathbf{Q}_{1}|\mathbf{Q}_{2}|\mathbf{Q}_{3}|\cdots|\mathbf{Q}_{i}\big{]}$

11
$\mathbf{B}=\big{[}{\mathbf{B}_{1}}^{T}|{\mathbf{B}_{2}}^{T}|\cdots|{\mathbf{B}_{i}}^{T}\big{]}$

12 end if

13

14 end for

Algorithm 8 Type-II RMF (column-wise) with a fixed error

Input: Template Matrix {$\mathbf{H}\in\mathbb{R}^{2N_{T}\times
N_{s}}:2N_{T}\leq N_{s}$}, $\langle{\delta\rho}/\rho\rangle$

Output: $\mathbf{Q}_{2N_{T}\times k}$, $\mathbf{\tilde{B}}_{k\times N_{s}}$

1 for  _$i=1,2,\cdots$_ do

2 $\mathbf{\Omega}_{i}\in\mathbb{R}^{N_{s}\times
b}\;:\mathbf{\Omega}_{ij}\;\in\mathcal{N}\big{(}0,1\big{)}\;$

3 for  _P = 1, 2_ do

4 ${\bar{\mathbf{H}}}_{N_{s}\times b}={\mathbf{H}^{T}}_{N_{s}\times
2N_{T}}\,\big{(}\mathbf{H}_{2N_{T}\times
N_{s}}\,{\mathbf{\Omega}_{i}}_{N_{s}\times b}\big{)}$

5 ${\mathbf{\Omega_{i}}}_{N_{s}\times b}=\textrm{qr}(\bar{\mathbf{H}})$

6

7 end for

$\mathbf{Q}_{i}=\mathbf{\Omega}_{i}$

// Consider top-$b$ basis vectors out of $\ell$ basis vectors.

8 $\mathbf{B}_{i}=\mathbf{H}\,\mathbf{Q}_{i}$

9 $\mathbf{H}=\mathbf{H}-\bf{Q}_{i}\,\mathbf{B}_{i}$

10 if  _$\|\mathbf{H}\|_{F}\geq\langle{\delta\rho}/\rho\rangle$_ then

11
$\mathbf{Q}=\big{[}\mathbf{Q}_{1}|\mathbf{Q}_{2}|\mathbf{Q}_{3}|\cdots|\mathbf{Q}_{i}\big{]}$

12 end if

13

14 end for

15$\text{Choose}\,k\,\text{based on}\,\ell.$

16 $\mathbf{Q}=\mathbf{Q}\big{[}:,1:k\big{]}$

17 for  _$\alpha$ = 1,2,3, $\cdots$, $2N_{T}$_ do

$\text{Find}\,\vec{\tilde{b}}^{\alpha}:H^{\alpha}\approx\vec{\tilde{b}}^{\alpha}_{1\times
k}\,\mathbf{Q}_{k\times N_{s}}$

// Using Matching Pursuit algorithm

${\mathbf{\tilde{B}}}_{2N_{T}\times
k}=[\vec{\tilde{b}}^{1}|\vec{\tilde{b}}^{2}|\cdots|\vec{\tilde{b}}^{\alpha}]$

// Row-wise stacking

18

19

20 end for

Algorithm 9 RMFwith fixed error + sparse coefficient using MP.

##  Appendix D Construction of a sparse reconstruction matrix (RMF\+ MP)

Alg-9 describes the hybrid scheme for obtaining basis vectors along with the
sparse coefficient by combining block-wise RMF and MP algorithm. A set of