# Comparative analysis of local angular rotation between the Ring Laser
Gyroscope GINGERINO and GNSS stations

Giuseppe Di Somma<EMAIL_ADDRESS>Dipartimento di Fisica,
Università di Pisa, largo B. Pontecorvo 3, I-56127 Pisa, Italy,  Istituto
Nazionale di Fisica Nucleare (INFN), sez. di Pisa, largo B. Pontecorvo 3,
I-56127 Pisa, Italy.  Nicolò Beverini  Dipartimento di Fisica, Università di
Pisa, largo B. Pontecorvo 3, I-56127 Pisa, Italy,  Giorgio Carelli
Dipartimento di Fisica, Università di Pisa, largo B. Pontecorvo 3, I-56127
Pisa, Italy,  Simone Castellano  Istituto Nazionale di Fisica Nucleare (INFN),
sez. di Pisa, largo B. Pontecorvo 3, I-56127 Pisa, Italy.  Gran Sasso Science
Institute, Viale Francesco Crispi 7, 67100 L’Aquila AQ, Italy  Roberto Devoti
Istituto Nazionale di Geofisica e Vulcanologia, Sez. ONT, sede di Roma, Italy
Enrico Maccioni  Dipartimento di Fisica, Università di Pisa, largo B.
Pontecorvo 3, I-56127 Pisa, Italy,  Istituto Nazionale di Fisica Nucleare
(INFN), sez. di Pisa, largo B. Pontecorvo 3, I-56127 Pisa, Italy.  Paolo
Marsili  Dipartimento di Fisica, Università di Pisa, largo B. Pontecorvo 3,
I-56127 Pisa, Italy,  Angela D.V. Di Virgilio  Istituto Nazionale di Fisica
Nucleare (INFN), sez. di Pisa, largo B. Pontecorvo 3, I-56127 Pisa, Italy.

###### Abstract

The study of local deformations is a hot topic in geodesy. Local rotations of
the crust around the vertical axis can be caused by deformations. In the Gran
Sasso area the ring laser gyroscope GINGERINO and the GNSS array are
operative. One year of data of GINGERINO is compared with the ones from the
GNSS stations, homogeneously selected around the position of GINGERINO, aiming
at looking for rotational signals with period of days common to both systems.
At that purpose the rotational component of the area circumscribed by the GNSS
stations has been evaluated and compared with the GINGERINO data. The
coherences between the signals show structures that even exceed 60$\%$
coherence over the 6-60 days period; this unprecedented analysis is validated
by two different methods that evaluate the local rotation using the GNSS
stations. The analysis reveals that the shared rotational signal’s amplitude
in both instruments is approximately $10^{-13}rad/s$, an order of magnitude
lower than the amplitudes of the signals examined using the coherence method.
The comparison of the ring laser data with GNSS antennas provides evidence of
the validity of the ring laser data for very low frequency investigation,
essential for fundamental physics test.

General Relativity, GNSS, Ring Laser Gyroscopes, Ground deformation.

††preprint: APS/123-QED

##  I INTRODUCTION

The GINGER (Gyroscopes IN GEneral Relativity) project is aiming at building an
array of ring laser gyroscopes (RLGs) in the underground Gran Sasso laboratory
(”Laboratori Nazionali del Gran Sasso”, LNGS). RLGs are top sensitivity
devices able to measure angular rotation rate exploiting the Sagnac effect.
GINGER will provide unique data for geophysical and fundamental physics
investigation Angela D.V. Di Virgilio, et al. (2023), especially for very low
frequency investigation. In this paper the data recorded by the RLG prototype
GINGERINO Belfi et al. (2017), presently active at LNGS, are compared with the
GNSS data in order to provide a test of the validity of the signal
reconstruction of the RLG signals. The GNSS technique is normally used for
studies of the 3D kinematic deformation patterns on the surface of the Earth
in local or global areas. In particular, GNSS time series allow to reconstruct
plate kinematics at global scale (e.g. Kreemer et al. (2014); Altamimi et al.
(2023)) and local strain rate tensors that identify deformation zones at any
scale on the Earth Okazaki et al. (2021); Serpelloni et al. (2022); Palano et
al. (2023). Here we exploit the possibility that a GNSS network and the Ring
Laser Gyroscope (RLG) facility may both sense a regional rotational signal of
the ground motion. Recent publications demonstrate that ground deformations in
seismically active areas are not only dominated by tectonic motions but are
also influenced by important hydrological processes Leone et al. (2023);
Pintori et al. (2021); Michel et al. (2021); Devoti et al. (2015), non-tidal
and atmospheric loading processes Mémin et al. (2020); Martens et al. (2020).
All these processes may, in principle, induce regional deformations with
periodicity ranging from a few days up to multi-annual signals, mainly driven
by climatic variations of the water cycle.
In this low frequency range the use of RLG data may be a promising tool to
monitor local geodetic signals showing their rotational nature. GINGERINO is
sensitive to rotations around the vertical axis of global and local origin
Tercjak and Brzeziński (2017). It is the first RLG based on a simple
mechanical structure which has been able to run in a continuous basis with
duty cycle larger than $90\%$.
The aim in this work is to define a relevant set of observable that can be
used to correlate GNSS and RLG observations, in order to isolate possible
regional geophysical processes that induce significant rotation rates. To our
knowledge this is the first attempt to observe common signals in these
instruments.
In the following the first section is dedicated to the description of the
GINGERINO apparatus and the analysis applied to the data to be compared to the
GNSS array. In the second, we illustrate the methods used to derive the
rotational component of GNSS stations. Next, we show some tests performed on
the analysis tool used to derive coherences and in the last section some
comments on the results and how to extend the analysis in the future.

##  II GINGERINO, the RLG prototype of the Gran Sasso

The Sagnac effect states that two counter-propagating light beams traveling in
a closed ring path complete the round trip with a difference in time
proportional to the absolute (or inertial) angular rotation rate of the rigid
closed path itself.
Based on this principle several devices have been developed for inertial
navigation, using optical fibers, optical cavity and cold atoms. A RLG is
composed of a high finesse polygonal optical cavity, where an active medium is
contained. Exciting the active medium the laser process is generated in the
two directions. If the optical cavity is rotating, the frequency of the two
counter-rotating beams is different, due to the different round trip time.
Then, combining outside the cavity the two beams on a photodetector, a beat
signal is generated at the so called Sagnac frequency$f_{s}$:

$f_{s}=4\frac{A}{P\lambda}\cdot\Omega\cos{\theta}$ (1)

where A is the area delimited by the light path, P the perimeter, $\Omega$ the
absolute value of the angular velocity of the optical cavity, $\lambda$ the
laser wavelength and $\theta$ is the relative angle between the direction (or
rotation axis) of $\vec{\Omega}$ and the area vector of the cavity $\vec{A}$.
GINGERINO, see Fig. 1 is rigidly attached to the ground, has square optical
cavity with 3.6m side, lying on the horizontal plane (vertical $\vec{A}$), so
the angular velocity is mainly $\Omega_{\oplus}$, the Earth rotation rate,
$\theta$ is the co-latitude. The sensitivity has been already demonstrated to
be in the range of pico-rad/s, with detection bandwidth of tens of Hz. In the
comparison with GNSS signals the GINGERINO data has been processed in order to
remove environment, tides and laser disturbances, as explained in the
following subsection.

Figure 1: GINGERINO is composed of a cross shaped granite monument, connected
to the floor through a reinforced concrete monument. On top of the granite
cross, four boxes are attached, where the mirrors are located. The four boxes
are connected by pipes where the light travels. In first plane, in between two
mirror boxes the discharge tube is visible, necessary to start and keep
working the laser.

### Details of the procedure for GINGERINO signals reconstruction

The radiation of the two counter-propagating beams transmitted by one mirror
are recombined and their interference is detected on a photodiode. GINGERINO
is connected to the Earth crust and the Earth rotation $\Omega_{\oplus}$ is
the largest signal, the interference is the beat note, $\omega_{m}$.
It is important to note that $\omega_{m}$ cannot be confused with the Sagnac
frequency, $\omega_{s}=2\pi f_{s}$, since the interference signal contains
small perturbations induced by the diffused radiation by the cavity mirrors
and the laser non linear dynamic. It is necessary to reconstruct $f_{s}$ using
the signals exiting the RLG Di Virgilio et al. (2017, 2019): the beat note and
the intensities of the two counter-propagating light beams (mono-beams). In
the present analysis the procedure to evaluate $f_{s}$ follows what described
in previous papers Di Virgilio et al. (2020); Di Virgilio, Angela D. et al.
(2021); Di Virgilio, A. D. V. et al. (2022), where more details can be found
about the theory of analysis procedure and how the terms necessary to
eliminate the typical laser disturbances called back–scatter and null shift
are evaluated and used. It is important to say that those laser systematics
affect $\omega_{m}$ mainly below 1Hz.
The optical cavity of GINGERINO is free running, and because of temperature
variations, the perimeter of the optical cavity changes, and as a consequence,
also the laser emitted wavelength changes. Accordingly, sometimes the light
beams change abruptly its frequency to go back in a position with higher gain.
Two possible effects take place: mode jumps and split mode operations. In the
first one the two counter–propagating beams change optical frequency at the
same time, a very fast transient happens with a few seconds duration. In the
second one only one beam changes frequency and the two beams have difference
in frequency of the free spectral range of the cavity, so that $\omega_{m}$
has very high frequency (MHz), and the GINGERINO apparatus is not able to
acquire. Split mode operation can last for several hours and the data are
lost, but are rare. The data are selected eliminating the portions with split
mode operation or affected by mode jumps, and others short duration large
signals, in the present analysis $20\%$ of the data are eliminated. GINGERINO
is isolated and based on a rather rigid structure, it is assumed that external
variations and signals will affect $\omega_{m}$ at first order, accordingly
the different expected signals and environmental disturbances are estimates
and subtracted with statistical means, basically using standard linear
regression. For the linear regression (LR) it is necessary to elaborate a set
of terms to be used in the linear regression. For the GINGERINO analysis the
following terms are used:

  * •

the laser dynamic has to be subtracted, so several terms are elaborated
following the theory (called LD)

  * •

the polar motions, called $F_{geo}$, the angular rotation around the vertical
axis as expected at the GINGERINO latitude and longitude estimated by the
international system IERS,

  * •

tides affect the apparatus, and the full list of tides elaborated using the
program $GOTIC2\\_{mod}$ , are used as LR terms and called $T$,

  * •

GINGERINO is a mechanical apparatus and accordingly temperature causes
geometry changes, deformation, orientation and the environmental monitors (EM)
of the GINGERINO apparatus, as tilts, temperature and pressure are used as
terms of the LR.

The first step of analysis provides a quantity $\omega_{s0}$, which well
reproduce the true Sagnac frequency above $1mHz$, and evaluates all the
required terms for the LR: the $LD$ (terms to subtract the Laser dynamic),
$F_{geo}$, $T$ and $EM$. The LR procedure is iterative, using each iteration 3
consecutive days, and keeping only the central day, so the coefficients of the
linear regressions changes day by day, but are evaluated using 3 days. It is
important to remind that the linear regression is done proposing about 10
terms, and selecting the ones with p-value below 0.1.

The full analysis releases three different outputs, see Fig. 2, the first
level uses $F_{geo}$ and $LD$, minimising the quantity $F1$ defined as:

$F1=\omega_{s0}-k\cdot F_{geo}-\mu\cdot LD$ (2)

$k$, and $\mu$ are the coefficients found by the linear regression to minimize
$F1$, and

$\omega_{s}=F1/k+F_{geo}$ (3)

Calling $\omega_{s}$ the true angular rotation recorded by GINGERINO, the term
$k$ mitigates for the 3 days used in the LR effects due to scale factor
variations due to temperature variation and change of inclination of the RLG
optical cavity. The first residuals $RES0$ represent the true angular
rotations cleaned by the polar motion effect:

$RES0=\omega_{s}-F_{geo}$ (4)

Last step $RES1$ minimises the residuals subtracting $T$ and $EM$:

$RES1=RES0-b\cdot T-c\cdot EM$ (5)

$b$ and $c$ are the coefficients of the linear regression.

The three quantities $\omega_{s}$, $RES0$ and $RES1$ are interpolated in order
to fill the gaps due to mode jumps and split modes and used in the comparison
with GNSS data.

Figure 2: General scheme of the analysis procedure.

##  III Coherence GINGERINO-GNSS stations signal

To derive a topographic trend over a period of 2 to 256 days, we initially
calculate the coherence between GINGERINO and individual stations. However,
upon examining the plot in Fig.3, which considers the contribution of
individual stations, no apparent common structures or trends are evident.
Consequently, we narrowed our focus to periods of 12-14 days and around 7 days
of interest, regions were typically tides and signals related to anthropic
activity are present.
The results in Figs. 4,5 demonstrate that there is no discernible spatial
pattern either by moving away from GINGERINO’s position or by considering
different areas to the east and west of the Apennines. Additionally, the
coherence value falls short of 40$\%$, a comparatively low value. As we will
delve into later, it’s worth noting that even two random noise signals can
attain a coherence value of 37$\%$ using the same analytical tools at our
disposal.

Figure 3: Coherence between individual stations and GINGERINO

Figure 4: The two plots show the topographic trend of coherence at the choice
of specific frequencies. The GINGERINO signal used is the one in which the
tides were resolved.

Figure 5: At the top the GINGERINO signal in which systematic laser
corrections and terrestrial rotational component, including polar motion and
Chandler wobble (obtained from IERS measurements IER ), were removed. At the
bottom we have the GINGERINO signal, obtained starting from the previous one,
in which we solved and subtracted the tides through the use of the
$GOTIC2\\_{mod}$ program GOT .

##  IV  Evaluate rotations using the GNSS stations

As previously explained, GINGERINO measures angular rotation rate around its
vertical axis, while the GNSS stations measures 3D displacements. To make a
step forward it is necessary to evaluate the rotation using the GNSS array.
For this purpose we employed two different technique.
The crust exhibits both strain and rotational deformation components, implying
that individual stations are spatially correlated. In the subsequent sections,
the two methods employed to determine the rotational component are described.
Firstly, we represent it as the straightforward rotational velocity around a
vertical axis. Secondly, we derive it as anti-symmetric rotational component
of the strain-rate matrix, which can be effortlessly traced back to the rotor
located at the GINGERINO position. One year of data of GINGERINO is compared
with the ones from the GNSS stations. When accompanied by ”$mjd$”, the time is
expressed in the modified Julian date.

### Rotational components approach

The first method to derive the rotational component from a network of GNSS
stations is to directly calculate the rotation vector associated with each
individual station and sum them appropriately; in our specific case, for
reasons mentioned above, the chosen rotation pole is the position of
GINGERINO. Using eq. 6 we can calculate the magnitude of the angular velocity
vector associated with each GNSS station:

$\omega_{1}=\frac{|v_{1}|}{|r_{1}|}\sin{(\alpha_{1}-\theta_{1})}$ (6)

where $|v_{1}|$ is the magnitude of the horizontal velocity of the individual
station seen in the Terrestrial Reference System integral to the Erthe’s crust
in which we decompose it into components of latitude $v_{N}$ and longitude
$v_{E}$ as: $|v_{1}|=\sqrt{v_{E}^{2}+v_{N}^{2}}$, while $\alpha_{1}$ is the
angle that the velocity vector forms with the local parallel, that is:
$\alpha_{1}=\arctan{\left(\frac{v_{N}}{v_{E}}\right)}$. Instead, $|r_{1}|$,
the absolute value of the station’s distance from GINGERINO and $\theta_{1}$,
the angle that the $\hat{r}_{1}$ vector forms with the parallel. They are
derived with the Matlab function ’distance’, specifying the coordinates of
GINGERINO, those of the station and the Earth model, which in our case is the
ellipsoid. In order to evaluate the mean value of the rotational velocities of
the stations array, the weights and thus the errors must be calculated;
details on the error estimation can be found in the Appendix A. We present the
signal obtained through this method in Fig. 6, wherein it is worth noting that
despite having a smaller standard deviation, the signal exhibits amplitudes of
comparable orders to that of GINGERINO and even some similar spikes. However,
to ascertain whether both signals capture the same phenomenon, a coherence
analysis is essential.

Figure 6: Shown in red is the tide-cleaned GINGERINO signal and in blue that
obtained by the rotation vector method, the period under examination is 359
days.

### Curl approach

A second method to derive the rotational component of the area circumscribed
by the GNSS constellations is to calculate the z-component of the rotor at the
GINGERINO position. To do this, we considered the equations linking the
velocities of the stations $v_{i}$ with their positions $x_{i}$ Allmendinger
et al. (2007):

$v_{i}=t_{i}+\frac{\partial v_{i}}{\partial x_{j}}x_{j}=t_{i}+e_{ij}x_{j}$ (7)

where $t_{i}$ is a simple integration constant that takes into account the
displacement at the origin of the coordinate system. The matrix $e_{ij}$ is
the displacement rate gradient tensor, and through the:

$e_{ij}=\epsilon_{ij}+\omega_{ij}=\frac{(e_{ij}+e_{ji})}{2}+\frac{(e_{ij}-e_{ji})}{2}$
(8)

we can divide it into the symmetrical part the infinitesimal strain rate
tensor and the anti-symmetrical part the rotation rate tensor. The latter,
apart from a factor of 2, corresponds to the z-component of the curl:

$\omega_{z}=\left(\frac{\partial v_{x}}{\partial y}-\frac{\partial
v_{y}}{\partial x}\right)$ (9)

Solving these equations involves a well-known inverse problem. To tackle this,
we employ a conventional least squares approach that effectively handles the
overconstrained nature of the problem.
To do this, we minimise the expression:

$\left[\overline{v}^{\prime}-A\cdot\overline{x}^{\prime}\right]^{2}$ (10)

where the matrix $A$, called the design matrix, takes into account how the
velocity components of the individual stations have been arranged, it is the
matrix transposition of the equations 7.
As previously stated, the crust supporting the stations adheres to the
principles of quasi-rigidity. While the stations are free to move, but they
display interconnections and correlations. Unlike the previous approach, we
compute the rotation values by incorporating the signals from the stations
simultaneously. To adequately address the correlation between these signals,
we construct a block matrix consisting of sub-matrices in a diagonal format:

$Corr=\begin{pmatrix}\sigma_{x}^{2}&C_{xy}&C_{xz}\\\
C_{xy}&\sigma_{y}^{2}&C_{yz}\\\ C_{xz}&C_{yz}&\sigma_{z}^{2}\end{pmatrix}$
(11)

where we have $C_{ij}=R_{ij}\sigma_{i}\sigma_{j}$, both the standard
deviations $\sigma_{i}$ and the correlation coefficients $R_{ij}$ are obtained
directly from the data provided by the GNSS stations, for example, in our
case, having used 22 GNSS stations we obtained a block matrix 66x66 and each
block on the diagonal corresponds to a single station with a 3x3 matrix given
by eq. 11. Let us show in Fig. 7 the comparison between the signal of
GINGERINO and the one obtained with the Curl method, we note how it is also
for the previous method that the magnitude of the signal amplitudes is the
same for both GINGERINO and the methods just right ($10^{-12}rad/s$), this is
surprising since, as can be seen from the top part of the Fig. 5, if we had
not resolved the tides for GINGERINO, as they are resolved in GNSS stations,
we would not have obtained the same magnitude but GINGERINO would be an order
of magnitude above. We show the direct comparison between the two signals in
Fig.8.

Figure 7: Shown in red is the tide-cleaned GINGERINO signal and in black that
obtained by the Curl method, the period under examination is the same as the
previous one. Figure 8: Comparison of the two method.

##  V Coherence analysis

Before delving into the analysis of the correlations between the two methods
and the GINGERINO signal, let’s conduct some tests using the MATLAB function
”mscohere.” These tests aim to determine the level of confidence required to
validate a coherence signal. Additionally, we’ll explore how low we can go
below the order of the signals under consideration to find a positive
coherence signal. Lastly, we will perform a fitting procedure on the station
velocities in an attempt to enhance the signal-to-noise ratio for the
implemented methods.

### Baseline for zero coherence: ”mscohere”

To determine the actual degree of coherence between the two signals, we
conducted tests using the ”mscohere” function along with simulated random
noises. Employing a Monte-Carlo simulation approach, we generated 100,000
instances of random noise, using the function $rand$ ******$rand$ generates
uniformly distributed random numbers in the interval [0, 1), but the same
result is achieved by using the $randn$ function, which a instead generates
random numbers distributed according to a standard normal distribution., and
calculated the coherence with respect to: another random noise, the Curl
signal, the rotation vector signal, and the GINGERINO signal. For each of the
aforementioned comparisons, we calculated the upper limit at a 95$\%$
confidence level for every frequency. The resulting data is presented in Fig.
9. Notably, despite considering the comparison with random noise, we observed
non-zero coherence levels. In fact, for the majority of ports, the frequency
band reached 37$\%$, furthermore, by choosing this value as zero we have a
conservative approach, however, it remains uncertain whether structures below
this confidence level are real or not. Hence, we establish this value as the
new baseline for zero coherence, above which we conclude that two signals
exhibit coherence associated with a physical signal; in addition, the
simulation between two random noises has a linear edge, while small structures
appear between the rotation and random noise signals.

Figure 9: We show here the coherence with 95$\%$ confidence intervals for all
simulations made in the comparison with synthetic random noise.

### Detection of a synthetic signal at a known frequency

We validated the analysis method by injecting a synthetic signal into all the
signals we considered. The objective was to assess the effectiveness of the
coherence method in producing positive results when applied to the signals
under analysis. Specifically, we introduced a sinusoidal signal with a
variable amplitude and a fixed frequency corresponding to a 7-day period. We
selected this specific period because the coherence method does not exhibit
any significant signals within that frequency range. By referring to Fig. 10,
we observe that the synthetic signal remains below the noise level established
in the previous test, as long as the amplitude does not exceed $9\cdot
10^{-14}rad/s$. Considering that the signals under examination have an
amplitude of $10^{-12}rad/s$, we can conclude that this analysis method allows
us to identify genuine coherence peaks for signals that are approximately one
order of magnitude lower than the maximum amplitude of the signal being
analyzed.

Figure 10: In the first quadrant, we show the coherence between GINGERINO and
the GNSS data, using the two analysis methods over an interval of 5 to 10
days. In the following quadrants we show by how this coherence varies by
adding to the signal a sinusoidal signal, that exhibited spikes over a
duration of 7 days, with amplitudes: $5\cdot 10^{-14}rad/s$, $7\cdot
10^{-14}rad/s$, $9\cdot 10^{-14}rad/s$, $2\cdot 10^{-13}rad/s$, $4\cdot
10^{-13}rad/s$, $6\cdot 10^{-13}rad/s$.

### Fit of linear velocities

In the previous described analysis we have determined speeds of GNSS stations
by comparing their different positions on consecutive days. The noise can be
reduced by a moving average approach by a linear fitting through multiple data
points.
For this purpose, we utilize the ”polyfit” function in Matlab. This function
takes into account the positions of individual stations along with their
corresponding errors, and outputs the angular coefficient and its associated
error for the best-fit line. We conduct the analysis repeatedly, considering
2, 3, 4, and 5 points at a time. Each time, we shift the analysis by one day
to maintain always the same number of points as the analysis without fitting.
The results of these analyses, incorporating the GINGERINO signal that
includes tides, are presented in Fig. 11. Additionally, we provide the results
in Fig. 12, where the tides have been resolved.

Figure 11: Coherence obtained without subtracting the contribution of the
tides, clear structures are evident with the 2-point fit.  Figure 12: In this
figure, we observe the coherence achieved through the resolution of tides in
GINGERINO. The usual tidal peaks are reduced, revealing previously hidden
structures with periods exceeding 20 days.

We observed that in the coherence analysis of the GINGERINO signal, which
includes tidal effects, a prominent peak emerges at approximately 14 days.
However, this peak vanishes when we perform the coherence analysis for all
signals, taking into account resolved tides. It is likely that the presence of
this peak is attributed to a frequency associated with the rotational dynamics
of tides, which is still present in both the GNSS stations and GINGERINO.
Once the tides are properly resolved, less resolved structures become
apparent, exhibiting a periodicity ranging from 20 to 60 days. These
structures do not originate from tides. Interestingly, they are only visible
when using the Curl method, possibly due to the distinct coordinate systems
employed by the two methods. While the rotation vector (depicted in blue)
solely considers the horizontal components, disregarding the vertical one, the
Curl method utilizes geocentric coordinates, inherently encompassing the
contribution from the radial component (vertical in the former method).
Consequently, these poorly resolved structures, coming out of coherence noise,
might be influenced by vertical phenomena impacting the Earth’s surface. In
the subsequent section, we will propose hypotheses to explicate these signals.
Additionally, it is worth noting that as we increase the number of data points
used for the fitting process, the coherence level of these structures
diminishes until it becomes indistinguishable from the previously calculated
noise. This observation can be attributed to the decrease in signal amplitude
resulting from the speed fit process. Consequently, the magnitudes of both the
Curl and rotation vector signals become comparable to signals that are not
discernible in coherence with the GINGERINO signal, which remains at a
magnitude of approximately $10^{-12}rad/s$. This comparison renders these
structures invisible in the coherence analysis, as demonstrated in the
synthetic signal test.

##  VI REMARKS AND CONCLUSIONS

Our analysis focuses on an initial comparison between two different
instruments: GNSS stations, which measure displacements, and GINGERINO, which
measures rotations. We have developed an analysis approach capable of
determining rotations for a constellation of GNSS satellites using two
different methods. From this comparison, we found that if we do not remove the
influence of tides from the GINGERINO signal, their contribution overwhelms
other structures that only become apparent when we also consider tides in the
RLG signal, similar to how we handled the stations data.
Concerning the coherence structure with a period of 20-60 days, it may be
attributed to a phenomenon acting vertically on the ground. This observation
is evident only in the Curl method, which, unlike the simple rotation vector
around vertical axis, incorporates the vertical component. One potential
factor might involve the alteration of pressure phenomena within atmospheric
loading. The reconstruction of the GNSS stations signal is affected by these
pressure variations, and GINGERINO, being influenced by pressure differences
exerted on the ground, experiences deformations in its surrounding area.
Additionally, these cycles are also impacted by the presence of the aquifer in
the Gran Sasso De Luca et al. (2018), which will be further investigated in
future studies.
The analysis enables us to detect signals with amplitudes one order of
magnitude lower than the average amplitudes of the analyzed signal
($10^{-12}rad/s$). In the GINGERINO RLG apparatus, being a single RLG with
vertical area vector, it is not possible to separate rotations and tilts of
the structure, see variable $\theta$ in eq. 1, which depend on the latitude,
but also on the orientation of the cavity with respect to the rotation axis.
Moreover, being based on a simple mechanical structure, it is limited by
rotations of instrumental origin, due to the fact the the laser cavity frame
is not rigid enough. Based on this experience a new design has been developed
to make tighter the cavity and avoid rotations of instrumental origin Basti et
al. (2021)
At present a more powerful apparatus, called GINGER, is under construction. It
involves the implementation of two larger Ring Laser Gyroscopes (RLGs). One of
these RLGs will operate at maximum Sagnac signal, called RLX, so with area
vector parallel to the Earth rotation axis. For its special alignment, RLX
will allow the identification of the orientation of the others RLG of the
array, allowing the separation of tilts and rotations. The expectation is to
reach 1 part in $10^{11}$ of the Earth rotation rate.
In this article, we have focused on a one-year data, as it represents the
longest period of available data from GINGERINO.
In the future, similarly to GNSS stations, multi-years data will be available
and we hope to further investigate and understand the phenomenon. By
incorporating GNSS stations, we can gain a better understanding of the
surrounding area of the RLG and improve our ability to identify signals that
are not solely related to rotation. Ultimately, this will lead to a
substantial improvement in the signal-to-noise ratio, enhancing the overall
accuracy of our findings.
The positive results of the coherence effectively validate the analysis
procedure of the RLG which is necessary to eliminate the Laser disturbance,
which poses severe limitation for period above 1000s.

## Appendix A: Error Propagation

Following are all the equations to derive angular velocity in the rotation
vector method. Starting from the linear velocities of the individual stations,
the angles defined previously to derive the perpendicular component of the
distance vector and the absolute value of the distance itself.

$\omega_{1}=\frac{|v_{1}|}{|r_{1}|}\sin{(\alpha_{1}-\theta_{1})}$ (12)

In order to do the weighted sum of the individual stations, we must propagate
the error in the appropriate way:

$\displaystyle\sigma_{\omega_{1}}^{2}$
$\displaystyle=\left(\frac{\partial\omega_{1}}{\partial
v_{1}}\right)^{2}\sigma_{v_{1}}^{2}+\left(\frac{\partial\omega_{1}}{\partial
r_{1}}\right)^{2}\sigma_{r_{1}}^{2}+$ (13)
$\displaystyle+\left(\frac{\partial\omega_{1}}{\partial\alpha_{1}}\right)^{2}\sigma_{\alpha_{1}}^{2}+\left(\frac{\partial\omega_{1}}{\partial\theta_{1}}\right)^{2}\sigma_{\theta_{1}}^{2}$

$\displaystyle\left(\frac{\partial\omega_{1}}{\partial v_{1}}\right)^{2}$
$\displaystyle=\left(\frac{\sin{(\alpha_{1}-\theta_{1})}}{|r_{1}|}\right)^{2}$
(14) $\displaystyle\left(\frac{\partial\omega_{1}}{\partial r_{1}}\right)^{2}$
$\displaystyle=\left(-\frac{|v_{1}|\sin{(\alpha_{1}-\theta_{1})}}{|r_{1}|^{2}}\right)^{2}$
(15)
$\displaystyle\left(\frac{\partial\omega_{1}}{\partial\alpha_{1}}\right)^{2}$
$\displaystyle=\left(\frac{|v_{1}|\cos{(\alpha_{1}-\theta_{1})}}{|r_{1}|}\right)^{2}$
(16)
$\displaystyle\left(\frac{\partial\omega_{1}}{\partial\theta_{1}}\right)^{2}$
$\displaystyle=\left(-\frac{|v_{1}|\cos{(\alpha_{1}-\theta_{1})}}{|r_{1}|}\right)^{2}$
(17)

$\sigma_{v_{1}}=\sqrt{\frac{v_{x}^{2}\sigma_{v_{x}}^{2}+v_{y}^{2}\sigma_{v_{y}}^{2}}{v_{x}^{2}+v_{y}^{2}}}$
(18)

$\sigma_{\alpha_{1}}=\frac{\sqrt{v_{y}^{2}\sigma_{v_{x}}^{2}+v_{x}^{2}\sigma_{v_{y}}^{2}}}{v_{x}^{2}+v_{y}^{2}}$
(19)

$\sigma_{v_{x}}=\sqrt{\sigma_{x_{i}}^{2}+\sigma_{x_{i+1}}^{2}}$ (20)

all other parts can be obtained analytically from GNSS measurements, instead,
$\sigma_{r_{1}}$ and $\sigma_{\theta_{1}}$ are evaluated with a parametric
bootstrap method Boos (2003), because they are obtained with the “distance”
function of Matlab.
In our approach, we began by generating normal functions that were centered
around the positions of the stations. The width of these functions was
determined by the sigma value obtained from our measurements. Using the
’distance’ function, we then calculated the distance to GINGERINO as well as
the corresponding angle. We repeated this procedure a million times, randomly
extracting values from the aforementioned normal distribution at each
iteration. This allowed us to obtain distributions of the measurements, and we
associated the standard deviations of these distributions with the errors in
the distance $r_{1}$ and angle $\theta_{1}$ .

## Acknowledgements

The authors express their profound appreciation to Prof. Giancarlo Cella and
Prof. Antonello Ortolan for their stimulating and enlightening discussions,
which significantly improved the quality of this article. Their valuable
insights and thoughtful feedback played an instrumental role in shaping our
overall understanding of the subject matter. Additionally, special thanks are
due to PhD students Alessandro Porcelli and Michele Cardelli, whose thought-
provoking questions and curiosities greatly contributed to the clarity of
fundamental steps in this analysis. Their input has been invaluable to the
development of this work. We also thank Aladino Govoni and Gaetano De Luca for
their contribution to making GINGERINO operational.

## References

  * Angela D.V. Di Virgilio, et al. (2023) Angela D.V. Di Virgilio, et al., to be published in MEMOCS (2023), URL arXiv:2209.09328.
  * Belfi et al. (2017) J. Belfi et al., Rev. Sci. Instrum. 88, 034502 (2017), eprint 1702.02789.
  * Kreemer et al. (2014) C. Kreemer, G. Blewitt, and E. Klein, Geochem. Geophys. Geosyst. 15, 3849–3889 (2014), URL https://doi.org/10.1002/2014GC005407.
  * Altamimi et al. (2023) Z. Altamimi, P. Rebischung, X. Collilieux, L. Métivier, and K. Chanard, J. Geod. 97, 47 (2023), URL https://doi.org/10.1007/s00190-023-01738-w.
  * Okazaki et al. (2021) T. Okazaki, Y. Fukahata, and T. Nishimura, Earth Planets Space 73, 153 (2021), URL https://doi.org/10.1186/s40623-021-01474-5.
  * Serpelloni et al. (2022) E. Serpelloni, A. Cavaliere, L. Martelli, F. Pintori, L. Anderlini, A. Borghi, D. Randazzo, S. Bruni, R. Devoti, P. Perfetti, et al., Front. Earth Sci. 10, 907897 (2022), URL https://doi.org/10.3389/feart.2022.907897.
  * Palano et al. (2023) m. Palano, S. Calcaterra, P. Gambino, B. Porfidia, and F. Sparacino, Results Geophys. Sci. 14, 100056 (2023), URL https://doi.org/10.1016/j.ringps.2023.100056.
  * Leone et al. (2023) G. Leone, N. D’Agostino, L. Esposito, and F. Fiorillo, Environmental Earth Sciences 82, 240 (2023), URL https://doi.org/10.1007/s12665-023-10905-3.
  * Pintori et al. (2021) F. Pintori, E. Serpelloni, L. Longuevergne, A. Garcia, L. Faenza, L. D’Alberto, A. Gualandi, and M. E. Belardinelli, J. Geoph. Res. 126, e2020JB020586 (2021), URL https://doi.org/10.1029/2020JB020586.
  * Michel et al. (2021) A. Michel, A. Santamaría-Gómez, J.-P. Boy, F. Perosanz, and S. Loyer, Remote Sens. 13, 4523 (2021), URL https://doi.org/10.1007/s10291-020-0959-3.
  * Devoti et al. (2015) R. Devoti, D. Zuliani, C. Braitenberg, P. Fabris, and B. Grillo, Earth and Planetary Science Letters 419, 134 (2015), URL https://doi.org/10.1016/j.epsl.2015.03.023.
  * Mémin et al. (2020) A. Mémin, J. Boy, and A. Santamaria-Gomez, GPS Solut. 24, 45 (2020), URL https://doi.org/10.1007/s10291-020-0959-3.
  * Martens et al. (2020) H. Martens, D. Argus, C. Norberg, G. Blewitt, T. Herring, A. Moore, W. Hammond, and C. Kreemer, J. Geod. 94, 115 (2020), URL https://doi.org/10.1007/s00190-020-01445-w.
  * Tercjak and Brzeziński (2017) M. Tercjak and A. Brzeziński, Pure and Applied Geophysics 174, 2719 (2017).
  * Di Virgilio et al. (2017) A. D. V. Di Virgilio, J. Belfi, W.-T. Ni, N. Beverini, G. Carelli, E. Maccioni, and A. Porzio, Eur. Phys. J. Plus 132, 157 (2017).
  * Di Virgilio et al. (2019) A. D. V. Di Virgilio, N. Beverini, G. Carelli, D. Ciampini, F. Fuso, and E. Maccioni, Eur. Phys. J. C 79, 573 (2019), eprint 1904.02533.
  * Di Virgilio et al. (2020) A. D. V. Di Virgilio, A. Basti, N. Beverini, F. Bosi, G. Carelli, D. Ciampini, F. Fuso, U. Giacomelli, E. Maccioni, P. Marsili, et al., Phys. Rev. Res. 2, 032069 (2020), URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.032069.
  * Di Virgilio, Angela D. et al. (2021) Di Virgilio, Angela D., Altucci, Carlo, Bajardi, Francesco, Basti, Andrea, Beverini, Nicolò, Capozziello, Salvatore, Carelli, Giorgio, Ciampini, Donatella, Fuso, Francesco, Giacomelli, Umberto, et al., Eur. Phys. J. C 81, 400 (2021), URL https://doi.org/10.1140/epjc/s10052-021-09199-1.
  * Di Virgilio, A. D. V. et al. (2022) Di Virgilio, A. D. V., Terreni, G., Basti, A., Beverini, N., Carelli, G., Ciampini, D., Fuso, F., Maccioni, E., Marsili, P., Kodet, J., et al., Eur. Phys. J. C 82, 824 (2022), URL https://doi.org/10.1140/epjc/s10052-022-10798-9.
  * (20) _Iers data: Compute time series of the earth orientation matrix / quaternion or derived parameters between two dates_ , URL https://hpiers.obspm.fr/eop-pc/index.php?index=rotation&lang=en.
  * (21) _Gotic2 $\\_{Mod}$: Modified gotic2, software for estimation of ocean tidal loading effects at subsurface observation_, URL https://www.gsj.jp/data/openfile/no0705/openfile0705indexEN.html.
  * Allmendinger et al. (2007) R. Allmendinger, R. Reilinger, and J. Loveless, Tectonics 26 (2007).
  * De Luca et al. (2018) G. De Luca, G. Di Carlo, and M. Tallini, Scientific Reports 8, 15982 (2018), ISSN 2045-2322, URL https://doi.org/10.1038/s41598-018-34444-1.
  * Basti et al. (2021) A. Basti et al., Eur. Phys. J. Plus 136, 537 (2021), eprint 2101.07686.
  * Boos (2003) D. D. Boos, Statistical science 18, 168 (2003).


# Graph signal denoising using $t$-shrinkage priors

Sayantan Banerjee  Operations Management and Quantitative Techniques Area,
Indian Institute of Management Indore, M.P., India  Weining Shen  Department
of Statistics, University of California, Irvine, CA, USA  Corresponding
author<EMAIL_ADDRESS>
###### Abstract

We study the graph signal denoising problem by estimating a piecewise constant
signal over an undirected graph. We propose a new Bayesian approach that first
converts a general graph to a chain graph via the depth first search
algorithm, and then imposes a heavy-tailed $t$-shrinkage prior on the
differences between consecutive signals over the induced chain graph. We show
that the posterior computation can be conveniently conducted by fully
exploring the conjugacy structure in the model. We also derive the posterior
contraction rate for the proposed estimator, and show that this rate is
optimal up to a logarithmic factor, besides automatically adapting to the
unknown edge sparsity level of the graph. We demonstrate the excellent
empirical performance of the proposed method via extensive simulation studies
and applications to stock market data.

> Keywords: Adaptive estimation; Bayesian shrinkage; depth first search; graph
> denoising; heavy tailed prior; posterior contraction rate.

##  1 Introduction

Let $G=(V,E)$ be an undirected graph with node (or, vertex) set
$V=\\{1,\ldots,n\\}$ and edge set $E$. The main focus of this paper is the
estimation of a piecewise constant signal
$\bm{\theta}=(\theta_{1},\ldots,\theta_{n})^{T}\in\mathbb{R}^{n}$ defined over
$G$. In particular, we consider the following model,

$\displaystyle
y_{i}=\theta_{i}+\epsilon_{i},\,i=1,\ldots,n;~{}~{}\epsilon_{i}\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}N(0,\sigma^{2}),$
(1)

where $\theta_{i}$ is the signal associated with the node $i$ for every
$i=1,\ldots,n$, $\epsilon_{1},\ldots,\epsilon_{n}$ are independent and
identically distributed (i.i.d) random Gaussian noises with an unknown
variance parameter $\sigma^{2}$, and $y_{1},\ldots,y_{n}$ are the observed
data points. By assuming that $\bm{\theta}$ has a piecewise constant
structure, we mean that for any pair of indexes $(i,j)\in E$, there is a
chance that the associated signals satisfy $\theta_{i}=\theta_{j}$.

Model (1) is often referred to as a graph denoising model, and it is closely
connected to the vast literature on change point detection (Tartakovsky et
al., 2014; Brodsky and Darkhovsky, 2013) and shape constrained regression
(Guntuboyina et al., 2018), especially when $G$ is a chain graph. In general,
the graph denoising model has been widely studied in various areas such as
image segmentation and denoising (Besag, 1986), signal processing (Gavish et
al., 2010), and network analysis (Crovella and Kolaczyk, 2003; Shuman et al.,
2013). In the literature, regularization has become a popular strategy to
solve the graph denoising problem (Barron et al., 1999; Birgé and Massart,
2007). In particular, Fan and Guan (2018) proposed an $\ell_{0}$-edge
denoising method by considering the following minimization problem,

$\displaystyle\arg\min_{\bm{\theta}\in\mathbb{R}^{n}}\dfrac{1}{2}\|\bm{y}-\bm{\theta}\|_{2}^{2}+\lambda\sum_{(i,j)\in
E}I(\theta_{i}\neq\theta_{j}),$ (2)

where $\bm{y}=(y_{1},\ldots,y_{n})^{T}$ is the observed data, and $\lambda>0$
is the tuning parameter. Also, throughout the paper, we use the notation
$\|\bm{x}\|_{2}$ to denote the $L_{2}$-norm of a vector $\bm{x}$, and
$I(\cdot)$ is the indicator function. Due to the combinatorial nature of the
objective function in (2), it is usually very difficult to find the exact
solution to the minimization for general graphs. Fan and Guan (2018) hence
developed an efficient approximation algorithm for solving (2) within a
polynomial time and showed that the resulting solution achieved the same risk
bound (up to a constant factor) with that of the exact minimizer of (2).
Alternatively, Padilla et al. (2017) proposed a fused Lasso approach by
considering an $L_{1}$ convex relaxation of (2), i.e.,

$\displaystyle\arg\min_{\bm{\theta}\in\mathbb{R}^{n}}\dfrac{1}{2}\|\bm{y}-\bm{\theta}\|_{2}^{2}+\lambda\sum_{(i,j)\in
E}|\theta_{i}-\theta_{j}|.$ (3)

The key idea of that paper is to first perform a depth-first search (DFS)
(Tarjan, 1972) that converts an arbitrary graph to a linear chain, and then
apply the fused Lasso (Rudin et al., 1992; Tibshirani et al., 2005) for
denoising purpose (also known as total variation denoising). The authors also
provided a theoretical justification for using DFS, that is, they showed that
the total variation in the signal defined over the chain graph obtained from
DFS is at most twice as its total variation over the original graph. Based on
this result, the authors obtained an upper bound for the mean squared error
(MSE) of their proposed signal estimators and established minimax optimality
for the obtained bound with respect to tree-structured graphs.

Despite their success in computation and theory, it still remains unclear how
the aforementioned regularization methods can be used to conduct statistical
inference for graph signal denoising problems; and this is the main focus of
our paper. In particular, we propose a simple-yet-useful Bayesian solution to
graph denoising problem. The main idea is to first perform DFS that converts a
general graph to a chain graph, and then conduct Bayesian fusion estimation
over the DFS-induced chain graph by shrinking the difference between
consecutive pairs of signals towards $0$. One advantage of adopting a full
Bayesian framework is that it automatically accounts for the uncertainty from
both DFS and fusion estimation. Moreover, the inference can be conveniently
conducted based on posterior samples. For Bayesian fusion estimation, we
choose to use the recently studied $t$-shrinkage prior (Song and Liang, 2017;
Song and Cheng, 2019), which has been shown to have a better theoretical
(e.g., in terms of posterior contraction rate) and empirical performance
(Bhattacharya et al., 2015; Castillo et al., 2015) than the conventional
Laplace prior in Bayesian fused Lasso (Kyung et al., 2010). We show that the
posterior computation can be conveniently implemented using a Gibbs sampler
thanks to the conjugacy structure in our model. We also derive a posterior
contraction rate for the estimated signals with respect to the class of “edge
sparse” graphs. To be precise, we show that when $G$ satisfies $\sum_{(i,j)\in
E}I(\theta_{i}\neq\theta_{j})\leq s$, then with posterior probability tending
to one, the mean squared error for $\bm{\theta}$ is of order $s\log n/n$,
which only differs with the oracle rate $s/n$ (assuming that the true graph
sparsity structure is known) by a logarithmic factor in $n$. This result
provides a useful theoretical justification for our proposed method, and also
reveals a desirable adaptive estimation property since we do not require the
knowledge of the graph edge sparsity parameter $s$ when constructing the prior
while the posterior can automatically achieve the optimal rate of convergence
with respect to different values of $s$.

Despite the fast growing literature on change point detection and shape-
constrained regression for linear chains, the exact minimax rate for
estimating a piecewise constant or monotonic signal over the linear chain is
only obtained recently by Gao et al. (2020). In the Bayesian domain, several
approaches have been proposed to fit such models, such as the recursion method
(Hutter et al., 2007), Bayesian fused Lasso (Kyung et al., 2010), and Bayesian
trend filtering (Roualdes, 2015). The posterior contraction property for those
methods, however, remains largely unknown until recently. For example, Song
and Cheng (2019) proposed to fit a Bayesian piecewise constant model using a
heavy-tailed shrinkage prior and obtained the posterior contraction rate for
the estimated signals. Liu et al. (2020) considered a Bayesian piecewise
polynomial model using an empirical Bayesian approach and obtained some
posterior contraction rate and structure recovery results. Bayesian piecewise
constant models on general graphs, to the best of our knowledge, still remain
unsolved to date; and our theoretical contribution in this paper is to fill
this gap.

The rest of the paper is organized as follows. We describe the proposed
methodology in Section 2 and discuss an efficient posterior computational
algorithm in Section 3. We then present the theoretical results in Section 4,
and demonstrate the excellent empirical performance of our method via
simulations and real data examples in Sections 5 and 6, respectively. We
discuss several future work directions in Section 7. Additional numerical
results and proof of theorems are given in the Appendix.

##  2 Bayesian Model Specification

We propose a two-step Bayesian approach for solving the piecewise constant
model in (1). We first perform a depth-first search (DFS) over the graph $G$
to convert it into a chain graph, denoted by $G_{C}=(V,E_{C})$, where
$V=\\{1,\ldots,n\\}$ is the same set of nodes in $G$, and $E_{C}$ represents a
new set of $(n-1)$ edges in $G_{C}$. Denote the root node for the DFS search
by $r\in V$ and let $\theta_{r}$ be its associated signal. Since the choice of
the root node $r$ is random, the DFS-induced graph $G_{C}$ is not unique.
However, in our implementation, we can choose a set of different root nodes
and then pool the corresponding posterior samples of $\bm{\theta}$ together
for inference purpose. Doing so will also account for the randomness
associated with DFS in the final results. As discussed in Padilla et al.
(2017), there are two benefits by converting the original graph to a chain
graph using DFS. The first reason is the obvious computational benefit - by
performing DFS, whose computational complexity is of order $O(|E|)$, we only
need to search over a chain graph instead of a potentially complex graph $G$.
The second reason is that the estimated $\bm{\theta}$ obtained from $G_{C}$ is
essentially as accurate as that from the original graph $G$ because the total
variation in $\bm{\theta}$ over the DFS-induced chain graph is bounded by the
twice of that on the original graph. This property provides an intuitive
justification for using DFS and will later help us derive the posterior
contraction property for the proposed Bayesian estimator.

Once we obtain the DFS-induced graph $G_{C}$, the next step is to conduct
Bayesian analysis of $\bm{\theta}$ over $G_{C}$. Note that the fused Lasso
penalty term in (3) has a natural Bayesian analogue as the negative logarithm
of the prior density corresponding to a Laplace distribution. However, it has
been shown in the literature (Bhattacharya et al., 2015; Castillo et al.,
2015) that the induced posterior corresponding to a Laplace prior may have
undesired posterior contraction property (e.g., a sub-optimal contraction
rate). Therefore, motivated by the recent work in Song and Liang (2017) and
Song and Cheng (2019), we consider a heavy-tailed prior for the difference
between signals over edges to facilitate better convergence. In particular, we
use a $t$-shrinkage prior on the successive differences of signals, i.e.,
$\theta_{i}-\theta_{j}$ for every $(i,j)\in E_{C}$ defined on the induced
chain graph. Note that we also need to specify a prior on the signal
corresponding to the root node $\theta_{r}$; and we use a normal prior for
convenience. The prior specification is thus given as:

$\displaystyle\theta_{r}\mid\sigma^{2}$ $\displaystyle\sim$ $\displaystyle
N(0,\lambda_{0}\sigma^{2}),$
$\displaystyle(\theta_{i}-\theta_{j})\mid\sigma^{2}$
$\displaystyle\stackrel{{\scriptstyle\text{ind}}}{{\sim}}$ $\displaystyle
t_{\nu}(m\sigma^{2}),~{}~{}\text{for every~{}}~{}(i,j)\in E_{C},$
$\displaystyle\sigma^{2}$ $\displaystyle\sim$
$\displaystyle\text{IG}(a_{\sigma},b_{\sigma}),$ (4)

where IG denotes an inverse Gamma distribution, $\lambda_{0}$, $a_{\sigma}$
and $b_{\sigma}$ are positive hyperparameters, and $t_{\nu}(\xi)$ denotes a
$t$-distribution with $\nu$ degrees of freedom and scale parameter $\xi$.
Similarly with the Laplace distribution, the $t$-distribution can be expressed
as a Normal-scale mixture with inverse Gamma mixing weights. Therefore we can
further write the prior in (4) as

$\displaystyle\theta_{r}\mid\sigma^{2}$ $\displaystyle\sim$ $\displaystyle
N(0,\lambda_{0}\sigma^{2}),$
$\displaystyle(\theta_{i}-\theta_{j})\mid\sigma^{2},\lambda_{k_{(i,j)}}$
$\displaystyle\stackrel{{\scriptstyle\text{ind}}}{{\sim}}$ $\displaystyle
N(0,\lambda_{k_{(i,j)}}\sigma^{2}),~{}~{}\text{for every~{}}~{}(i,j)\in
E_{C},$ $\displaystyle\lambda_{k_{(i,j)}}$
$\displaystyle\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}$
$\displaystyle\text{IG}(a_{t},b_{t}),~{}~{}\text{for every~{}}~{}(i,j)\in
E_{C},$ $\displaystyle\sigma^{2}$ $\displaystyle\sim$
$\displaystyle\text{IG}(a_{\sigma},b_{\sigma}),$ (5)

where the hyperparameters $a_{t}$ and $b_{t}$ satisfy $\nu=2a_{t}$ and
$m=\sqrt{b_{t}/a_{t}}$, and $k_{(i,j)}$ is the index for the edge that
connects nodes $i$ and $j$.

The choice of hyperparameters in a hiearchical Bayesian model setup plays a
crucial role. In our case, the hyperparamters corresponding to the Inverse
Gamma prior for the variance parameter $\sigma^{2}$ are chosen in such a way
that the prior is almost non-informative, e.g., $a_{\sigma}=b_{\sigma}=0.5.$
For the choice of the hyperparameters $a_{t}$ and $b_{t}$, we follow the
approach in Song and Cheng (2019), i.e., we choose $a_{t}=2$, and the scale
parameter $m$ in the $t$-prior is chosen so as to satisfy
$P(|t_{\nu}(m)|\geq\sqrt{\log(n)/n})\approx n^{-1}$. In the numerical studies,
we also compare our method with Laplace-prior based method, in which we choose
the hyperparameter $\lambda$ as $\lambda=\sqrt{2\log(n)}$ and the value of
$\lambda_{0}$ as $5$ based on the recommendation in Song and Cheng (2019).

##  3 Posterior computation

One of the nice features of the proposed method is that the posterior
calculation can be implemented in a relatively easy way by exploring the
conjugacy structure in the model. Given the Gaussian likelihood and prior
specification in the previous section, the conditional posterior distribution
of the parameters of interest can be conveniently obtained by the Gibbs
sampling as follows,

$\displaystyle\lambda_{k_{(i,j)}}\mid\cdot$ $\displaystyle\sim$
$\displaystyle\text{IG}\left(a_{t}+\dfrac{1}{2},b_{t}+\dfrac{(\theta_{i}-\theta_{j})^{2}}{2\sigma^{2}}\right),\,k_{(i,j)}\in\\{1,\ldots,n-1\\},~{}\text{for
every}~{}(i,j)\in E_{C},$ $\displaystyle\sigma^{2}\mid\cdot$
$\displaystyle\sim$
$\displaystyle\text{IG}\left(a_{\sigma}+n,b_{\sigma}+\dfrac{\theta_{r}^{2}}{2\lambda_{0}}+\dfrac{1}{2}\|\bm{y}-\bm{\theta}\|_{2}^{2}+\sum_{(i,j)\in
E_{C}}\dfrac{(\theta_{i}-\theta_{j})^{2}}{2\lambda_{k_{(i,j)}}}\right),$
$\displaystyle\theta_{i}\mid\cdot$ $\displaystyle\sim$ $\displaystyle
N(\mu_{i},\nu_{i}),\,i=\\{1,\ldots,n\\}\setminus\\{r\\},$ (6)

where $\mid\cdot$ means conditional on the rest of other parameters and the
data, $r$ is the root node for DFS, and

$\displaystyle\nu_{i}^{-1}$ $\displaystyle=$
$\displaystyle\dfrac{1}{\sigma^{2}}+\dfrac{1}{\lambda_{k_{(i,j)}}\sigma^{2}}I((i,j)\in
E_{C})+\dfrac{1}{\lambda_{k_{(j,i)}}\sigma^{2}}I((j,i)\in E_{C}),$
$\displaystyle\mu_{i}$ $\displaystyle=$
$\displaystyle\nu_{i}\left(\dfrac{y_{i}}{\sigma^{2}}+\dfrac{\theta_{j}}{\lambda_{k_{(i,j)}}\sigma^{2}}I((i,j)\in
E_{C})+\dfrac{\theta_{j}}{\lambda_{k_{(j,i)}}\sigma^{2}}I((j,i)\in
E_{C})\right)$

for every $i\in\\{1,\ldots,n\\}\setminus\\{r\\}$, and $k_{(i,j)}$ is the edge
number ranging from $1$ to $(n-1)$ in $G_{C}$ such that $(i,j)\in E_{C}$.

##  4 Theoretical guarantees

In this section, we study the posterior concentration property for the
proposed method. Let
$\bm{\theta}^{*}=(\theta_{1}^{*},\ldots,\theta_{n}^{*})^{T}$ and $\sigma^{*}$
be the respective true values for the signal $\bm{\theta}$ and the noise
standard deviation. We make the following assumptions,

  1. (C1)

Define the number of blocks in $\bm{\theta}^{*}$ by $s=\sum_{(i,j)\in
E}I(\theta_{i}^{*}\neq\theta_{j}^{*})$. We assume $s=o(n/\log n)$.

  2. (C2)

For $a_{t}$ and $b_{t}$ defined in (2), we assume that $a_{t}/b_{t}=n^{c}$ for
some positive constant $c$.

  3. (C3)

For every $i,j$ such that $(i,j)\in E$, we assume
$\log(|\theta_{i}^{*}-\theta_{j}^{*}|/\sigma^{*})=O(\log n)$.

  4. (C4)

For every $i=1,\ldots,n$, we assume
$\theta_{i}^{*2}/(\lambda_{0}\sigma^{*2})+\log(\lambda_{0})=O(\log n)$, where
$\lambda_{0}$ is defined in (4).

Then we have the following theorem that gives the posterior contraction rate
of the proposed estimator. The proof is given in the Appendix.

######  Theorem 1.

Suppose that Conditions (C1)–(C4) hold, then for any constant $M>0$, the
posterior distribution of $\bm{\theta}$, denoted by $\Pi(\cdot\mid\bm{y})$,
satisfies

$\Pi(\|\bm{\theta}-\bm{\theta^{*}}\|_{2}/\sqrt{n}\leq M\sigma^{*}\sqrt{s\log
n/n}\mid\bm{y})\rightarrow 1,~{}~{}\text{as}~{}~{}n\rightarrow\infty,$

where the convergence holds in probability and also with respect to
$L_{1}$-distance.

Theorem 1 shows that the posterior convergence rate for the average
$L_{2}$-estimation error $\|\bm{\theta}-\bm{\theta^{*}}\|_{2}/\sqrt{n}$ is of
order $\sigma^{*}\sqrt{s\log n/n}$. Note that when the graph sparse structure
is known, then the oracle rate of contraction is $O(\sigma^{*}\sqrt{s/n})$.
Therefore the proposed Bayesian shrinkage method achieves the optimal rate of
convergence except the logarithm term in $n$, and this rate is adaptive to the
unknown graph edge sparsity level $s$. This result agrees with the findings in
Song and Cheng (2019) for linear chains, and also connects to the rate
obtained by DFS fused Lasso for general graphs in Padilla et al. (2017) in the
sense that those two rates are equivalent except the $\log n$ factor. It is
worthy discussing and comparing our assumptions (C1)–(C4) with those in
Padilla et al. (2017). In particular, we do not require the minimal signal
strength assumption, but we do need the additional assumptions on the prior
hyperparameters as in (C2) and the maximum signal size as in (C3) and (C4).
Assumptions (C1)-(C4) are in general mild and commonly used in the literature.
For example, (C1) refers to the usual “edge sparse” graphs in Fan and Guan
(2018) and the upper bound on $s$ seems reasonable given that the signal has a
length of $n$ in total. Assumption (C2) ensures the $t$-shrinkage prior to
assign sufficiently amount of prior mass to the tails for optimal convergence
as discussed in Song and Cheng (2019). Assumptions (C3) and (C4) essentially
require the true signals and their paired difference to be bounded by a
diverging sequence of numbers as $n$ increases, which is easily satisfied for
most examples in practice.

##  5 Simulation studies

We carry out extensive simulation studies to evaluate the empirical
performance of our proposed DFS $t$-fusion approach, and compare with two
alternative approaches: (i) Laplace fusion approach, where the Laplace prior
is used instead of the $t$-shrinkage prior in model (4) and the rest steps are
the same with our approach, and (ii) $L_{1}$ fusion, which is the DFS fused
Lasso approach as proposed in Padilla et al. (2017). We consider three types
of the graphs, linear chain graph, 2-D lattice graph, and tree-like graph.
When implementing our method, except for the linear chain graph, we consider
three random root notes for the DFS step and then pool the posterior samples
of $\bm{\theta}$ together for inference purpose. All the computations were
performed in R (R Core Team, 2019) on a computing server with two 512 GB RAM,
AMD Opteron 6378 processors. For evaluating the performance of our proposed
method and comparing with competing approaches, we consider the mean squared
error (MSE) and adjusted MSE, respectively defined as

$\displaystyle\mathrm{MSE}$ $\displaystyle=$
$\displaystyle\|\hat{\bm{\theta}}-\bm{\theta_{0}}\|_{2}^{2}/n,$
$\displaystyle\mathrm{adj.MSE}$ $\displaystyle=$
$\displaystyle\|\hat{\bm{\theta}}-\bm{\theta}_{0}\|_{2}^{2}/\|\bm{\theta}_{0}\|_{2}^{2},$

where $\hat{\bm{\theta}}$ is the estimated signal, and $\bm{\theta}_{0}$ is
the true signal.

###  5.1 Linear Chain graphs

We first consider the estimation of a piecewise constant function on linear
chain graphs. This problem corresponds to the classical change-point detection
problem, in which case we can directly use the true graph structure without
running the DFS algorithm. For a chain with $n=100$ nodes, three different
signals with varyingly spaced pieces are considered – (i) 10 evenly spaced
pieces, each having a length of 10, (ii) 10 unevenly spaced pieces, with the
smaller components having a length of 5, and (iii) 10 very unevenly spaced
pieces, with the smaller components having a length of 2. The response
variable $y$ is simulated from $N(\bm{\theta_{0}},\sigma^{2}I_{n})$, where
$\bm{\theta_{0}}$ is the true signal vector and $\sigma\in\\{0.1,0.3,0.5\\}.$

We first present the estimated signals obtained from three methods in Figure 1
when the true signal is evenly spaced and $\sigma=.1$. It is clear that both
the proposed $t$-fusion prior and fused Lasso are capable of recovering the
underlying piecewise constant block structure very well while Laplace prior
does not have a satisfactory performance in identifying those “blocks”. Also,
the proposed approach can conveniently obtain posterior credible bands that
are much narrower than those from Laplace prior while $L_{1}$ fused Lasso is
only able to provide point estimations. These observations confirm the
excellent performance of the use of $t$ prior in terms of shrinkage effect and
signal recovery. We also present pictures for other cases (unevenly spaced
true signal with higher $\sigma$ values) in Figures 6–10 in the Appendix.
Similar findings are observed.

We then summarize the mean squared error (MSE) and adjusted MSE, and their
associated standard errors (in brackets) over 100 Monte-Carlo replications in
Table 1. We find that our proposed method has achieved the highest accuracy
for most cases when the noise level is reasonably low, i.e., $\sigma=.1$ and
$.3$, but does not perform as well when $\sigma=.5$. One possible explanation
is that as the noise level increases, the observed data deviates more from a
piecewise constant signal and hence it becomes more challenging to recover the
true block structure. As our method uses a heavy-tailed shrinkage prior, it
naturally prefers a more sparse structure with less blocks and hence the
estimation accuracy deteriorates.

|  |  $t$-fusion | Laplace fusion |  $L_{1}$ fusion
---|---|---|---|---
Signal | $\sigma$ | MSE | adj MSE | MSE | adj MSE | MSE | adj MSE
| 0.1 | 0.002(0.001) | 0.001(0.000) | 0.232(0.008) | 0.081(0.016) | 0.005(0.003) | 0.002(0.001)
even | 0.3 | 0.037(0.017) | 0.013(0.006) | 0.265(0.020) | 0.093(0.019) | 0.037(0.020) | 0.013(0.007)
| 0.5 | 0.262(0.308) | 0.092(0.108) | 0.323(0.036) | 0.113(0.024) | 0.098(0.051) | 0.034(0.018)
| 0.1 | 0.002(0.001) | 0.001(0.001) | 0.274(0.008) | 0.113(0.017) | 0.006(0.003) | 0.002(0.001)
uneven | 0.3 | 0.043(0.014) | 0.018(0.007) | 0.309(0.018) | 0.128(0.020) | 0.039(0.019) | 0.016(0.008)
| 0.5 | 0.524(0.457) | 0.216(0.186) | 0.368(0.032) | 0.152(0.025) | 0.099(0.048) | 0.041(0.019)
| 0.1 | 0.007(0.006) | 0.003(0.003) | 0.337(0.006) | 0.156(0.019) | 0.009(0.002) | 0.004(0.001)
v uneven | 0.3 | 0.192(0.240) | 0.088(0.110) | 0.365(0.015) | 0.169(0.022) | 0.063(0.023) | 0.029(0.010)
| 0.5 | 0.989(0.111) | 0.459(0.075) | 0.411(0.025) | 0.191(0.028) | 0.139(0.059) | 0.064(0.027)
Table 1: Linear Chain Graphs: MSE and adjusted MSE (their associated standard errors) for $t$-fusion, Laplace fusion, and $L_{1}$ fusion method, when the true signal is evenly spaced (“even”), unevenly spaced (“uneven”), and very unevenly spaced (“v uneven”). |
---|---
(a) True signal | (b) $t$-fusion estimates
|
(c) Laplace fusion estimates | (d) $L_{1}$ fusion estimates
Figure 1: Signal denoising performance for linear chain graphs with evenly
spaced signals and $\sigma=0.1.$: data in red dots, point estimates in black
dots, and 95% credible bands of the Bayesian procedures in green.

###  5.2 2-D lattice graphs

For the second simulation scenario, we explore 2-D lattice graphs. We consider
a grid graph of dimension 16, with piecewise constant signals given by
$\theta_{0,ij}=\kappa$ for points $(i,j)$ within a distance of 4 units from
the center of the grid, and $\theta_{0,ij}=0$ otherwise. The observations
$y_{ij}$ are simulated independently from $N(\theta_{0,ij},\sigma^{2})$ with
$\sigma=0.3$, and the parameter $\kappa$ that controls the signal strength is
chosen as $\kappa\in\\{1,5,10\\}$.

We first present the plots of the true signal and the recovered signals
obtained from three methods in Figure 2 when $\kappa=10$. Similarly with the
linear chain graphs, we find that both our method and $L_{1}$ fused Lasso
manage to correctly identify the latent structure in the signal while Laplace
prior method has an unsatisfactory performance in identifying the boundary of
the signals. Plots for $\kappa=\\{1,5\\}$ are provided as Figure 11 and 12 in
the Appendix; and we observe a similar pattern there.

|
---|---
(a) True signal | (b) Noisy signal
|
(c) $t$-fusion estimates | (d) Laplace fusion estimates
(e) $L_{1}$ fusion estimates
Figure 2: Signal denoising performance for 2-D lattice graphs with strong
signal strength ($\kappa=10$).

Next we summarize both the MSE and adjusted MSE for all three methods in Table
2 based on 100 Monte-Carlo replications. As the signal becomes stronger, all
three methods have a better accuracy in terms of smaller adjusted MSE values.
Overall, $L_{1}$ fusion has a slightly lower estimation error than our method
although this advantage becomes less noticeable as the signal becomes strong.
The use of $t$ prior has a much better signal estimation accuracy than that of
the Laplace prior, which again confirms the benefit of using a heavy-tailed
shrinkage prior when the true signal itself is indeed piecewise constant.

|  $t$-fusion | Laplace |  $L_{1}$ fusion
---|---|---|---
Signal | MSE | adj MSE | MSE | adj MSE | MSE | adj MSE
weak | 0.031(0.010) | 0.039(0.012) | 0.107(0.007) | 0.131(0.009) | 0.012(0.005) | 0.014(0.006)
medium | 0.047(0.015) | 0.009(0.003) | 0.509(0.015) | 0.100(0.003) | 0.012(0.005) | 0.002(0.001)
strong | 0.028(0.033) | 0.001(0.002) | 1.920(0.031) | 0.094(0.002) | 0.012(0.005) | 0.001(0.000)
Table 2: 2-D lattice Graphs: MSE and adjusted MSE (their associated standard errors) for $t$-fusion, Laplace fusion, and $L_{1}$ fusion method, for weak ($\kappa=1$), medium ($\kappa=5$), and strong ($\kappa=10$) signal strength. |
---|---
(a) True signal | (b) Noisy signal
|
(c) $t$-fusion estimates | (d) Laplace fusion estimates
(e) $L_{1}$ fusion estimates
Figure 3: Signal denoising performance for tree graphs.

###  5.3 Tree graphs

Finally we consider signals defined on tree-structured graphs. We consider
four trees, namely, $T_{1},\ldots,T_{4}$, each having 50 nodes and 3 children
per node. A connected tree-like graph $G_{T}$ having 200 nodes is then formed
by randomly selecting a node from $T_{i}$ and $T_{i+1},\,i=1,2,3$ and placing
an edge between the respective nodes. The true signal values $\theta_{0}$
defined on the trees $T_{1},\ldots,T_{4}$ are given by $1,-1,4,-4$
respectively. We generate the data on the graph $G_{T}$ from
$N(\theta_{0},\sigma^{2})$ with $\sigma=0.3$. The signal denoising performance
is presented in Figure 3 for three methods. Similarly with the previous two
scenarios, both our method and $L_{1}$ fusion have an excellent signal
recovery performance while Laplace fusion tends to produce over-clustered
results. This finding is also confirmed in Table 3, where the estimation
accuracy of $t$-fusion and $L_{1}$-fusion is better than that of the Laplace
prior.

Method | MSE | adj.MSE
---|---|---
$t$-fusion | 0.004 (0.003) | 0.0004 (0.0003)
Laplace | 0.049 (0.006) | 0.0060 (0.0007)
$L_{1}$-fusion | 0.005 (0.002) | 0.0006 (0.0003)
Table 3: MSE and adjusted MSE (standard error) for tree graphs.

##  6 Real data application

Change-point detection in financial data is a widely studied problem. We
consider the problem of change-point detection in the cumulative stock index
returns of two different indices in two different periods, one for the Dow
Jones Industrial Average (DJI) and the other for the Shanghai Stock Exchange
Index (SSE)111The datasets are publicly available at
https://finance.yahoo.com/. We carefully choose the two different time periods
so as to include known stock-market crash phenomena. To be precise, we study
the DJI for the period of January 2007 to January 2010, that encompasses the
2008 Global Recession period, and the SSE for the period January 2013 to April
2016, that encompasses the 2015 China Stock Market crash.

We apply the graph signal denoising method using both the $t$-fusion prior and
the Laplace prior, along with the frequentist fused Lasso method on both data
sets. The results are shown in Figure 4 and 5. We find that the fused Lasso
solution suffers from severe over-fitting problems, and a similar pattern is
observed in the Bayesian Laplace fusion prior as well (see (a) and (b) in
Figure 4 and 5). In contrast, the proposed Bayesian $t$-fusion method produces
almost piecewise constant estimates. Following the suggestions of Song and
Cheng (2019), we “sparsify” the continuous posterior accounting for the prior
concentration. Specifically, we estimate the quantity
$I(\theta_{i}=\theta_{j})$ for $i\neq j$ by
$I(|\hat{\theta}_{i}-\hat{\theta}_{j})|/\hat{\sigma}\leq mt_{1/2n})$, where
$t_{1/2n}$ is the $(1-1/2n)$-quantile of a $t$ distribution with $2a_{t}$
degrees of freedom and $\hat{\theta},\hat{\sigma}$ are Bayes estimates given
by the respective posterior means. It is clear that this sparsification step
is very useful for locating the blocks (or approximate pieces) precisely,
which in turn, provide the change points in the times series.

The change points obtained from the sparsification of the estimates capture
major financial events during the time periods. During the global financial
crisis of 2007-2008, large financial institutions across the United States
failed around mid-September 2008, followed by failure of major banks in
Europe. Price of stocks and other commodities plummeted sharply across the
globe. The-then U.S. President George W. Bush created a Troubled Asset Relief
Program (TARP) for purchasing failing bank assets, via signing the Emergency
Economic Stabilization Act into law. The effects of this 2008 Global Recession
were disastrous, with the world economy taking a serious blow. The resulting
period is referred to as the ‘Great Recession’, that officially spanned the
period of December 2007 to June 2009 in the United States. The DFS $t$-fusion
method captures these major changes in the movement of stocks, with change-
points detected at the beginning of the recession (late 2007), followed by
multiple change-points during late 2008, and at periods when the recession
officially ended.

The 2015-2016 Chinese Stock Market crash started with a stock market bubble
during mid-June 2015, after which the prices crashed sharply throughout July
and August 2015. Following a steep sell-off in the Chinese stock market on
January 2016, global markets felt severe tremors. The proposed method of ours
captures these major change-points in the stock market movements, with several
change-points detected en-route June 2015, indicating the rapid speed at which
the stock market bubble was building up, followed by the series of crashes
post June 2015. The January 2016 global meltdown period has been captured as
well.

In a recent work on dynamic change-point detection in financial networks,
Banerjee and Guhathakurta (2020) analyzed fully connected time-varying
networks based on major stock returns in the Asia-Pacific region along with
stock market indices from the U.S. and the U.K. Their method was able to
detect major financial events in different time periods, including the two
regimes considered above. Our proposed approach works with a single time-
series at a time, but alongside detecting known stock-market crash phenomena,
it additionally captures ripples in the market both apriori and aposteriori
the critical events. Thus, our method may be utilized for serving as potential
warning prior to major disruptions in global financial markets.

---
(a) Fused Lasso solution
(b) Bayesian Laplace-fusion model
(c) Proposed DFS $t$-fusion method
Figure 4: Estimation results for cumulative Stock returns for the Dow Jones
Industrial Average (DJI) in the period of January 2007 to January 2010 by
three methods. The black lines refer to the observed data. The green lines in
(c) indicate the dates where a significant change in the means occur.
---
(a) Fused Lasso solution
(b) Bayesian Laplace-fusion model
(c) Proposed DFS $t$-fusion method
Figure 5: Estimation results for cumulative Stock returns for the Shanghai
Stock Exchange (SSE) in the period of January 2013 to April 2016 by three
methods. The black lines refer to the observed data. The green lines in (c)
indicate the dates where a significant change in the means occur.

##  7 Discussion

In this paper, we propose a new Bayesian approach for analyzing the piecewise
constant model over general graphs. The proposed method makes efficient use of
the depth first search (DFS) algorithm and heavy-tailed shrinkage prior.
Several future working directions remain open. First, our proposed method in
general works for any graph-searching algorithms (instead of restricting to
DFS); and it will be of interest to compare and investigate the effect of
different graph-searching algorithms on the signal estimation performance.
Second, the theory studied in this paper mainly focuses on signal estimation.
It will be of interest to investigate the structure recovery property in the
graph. We expect the results in a recent work (Kim and Gao, 2020) to be
helpful. Lastly, it will be of interest to consider more general residual
distributions other than the normal distribution in model (1).

## Acknowledgements

S.B. is partially supported by DST INSPIRE Faculty Award, Govt. of India,
Grant No. 04/2015/002165, and IIM Indore Young Faculty Research Chair grant.

## References

  * Banerjee and Guhathakurta (2020) Banerjee, S. and Guhathakurta, K. (2020), “Change-point analysis in financial networks,” Stat, 9, e269.
  * Barron et al. (1999) Barron, A., Birgé, L., and Massart, P. (1999), “Risk bounds for model selection via penalization,” Probability theory and related fields, 113, 301–413.
  * Besag (1986) Besag, J. (1986), “On the statistical analysis of dirty pictures,” Journal of the Royal Statistical Society: Series B (Methodological), 48, 259–279.
  * Bhattacharya et al. (2015) Bhattacharya, A., Pati, D., Pillai, N. S., and Dunson, D. B. (2015), “Dirichlet–Laplace priors for optimal shrinkage,” Journal of the American Statistical Association, 110, 1479–1490.
  * Birgé and Massart (2007) Birgé, L. and Massart, P. (2007), “Minimal penalties for Gaussian model selection,” Probability theory and related fields, 138, 33–73.
  * Brodsky and Darkhovsky (2013) Brodsky, E. and Darkhovsky, B. S. (2013), Nonparametric methods in change point problems, vol. 243, Springer Science & Business Media.
  * Castillo et al. (2015) Castillo, I., Schmidt-Hieber, J., Van der Vaart, A., et al. (2015), “Bayesian linear regression with sparse priors,” The Annals of Statistics, 43, 1986–2018.
  * Crovella and Kolaczyk (2003) Crovella, M. and Kolaczyk, E. (2003), “Graph wavelets for spatial traffic analysis,” in IEEE INFOCOM 2003. Twenty-second Annual Joint Conference of the IEEE Computer and Communications Societies (IEEE Cat. No. 03CH37428), IEEE, vol. 3, pp. 1848–1857.
  * Fan and Guan (2018) Fan, Z. and Guan, L. (2018), “Approximate $\ell_{0}$-penalized estimation of piecewise-constant signals on graphs,” The Annals of Statistics, 46, 3217–3245.
  * Gao et al. (2020) Gao, C., Han, F., Zhang, C.-H., et al. (2020), “On estimation of isotonic piecewise constant signals,” Annals of Statistics, 48, 629–654.
  * Gavish et al. (2010) Gavish, M., Nadler, B., and Coifman, R. R. (2010), “Multiscale Wavelets on Trees, Graphs and High Dimensional Data: Theory and Applications to Semi Supervised Learning.” in ICML, pp. 367–374.
  * Guntuboyina et al. (2018) Guntuboyina, A., Sen, B., et al. (2018), “Nonparametric shape-restricted regression,” Statistical Science, 33, 568–594.
  * Hutter et al. (2007) Hutter, M. et al. (2007), “Exact Bayesian regression of piecewise constant functions,” Bayesian Analysis, 2, 635–664.
  * Kim and Gao (2020) Kim, Y. and Gao, C. (2020), “Bayesian model selection with graph structured sparsity,” Journal of Machine Learning Research, 21, 1–61.
  * Kyung et al. (2010) Kyung, M., Gill, J., Ghosh, M., Casella, G., et al. (2010), “Penalized regression, standard errors, and Bayesian lassos,” Bayesian Analysis, 5, 369–411.
  * Liu et al. (2020) Liu, C., Martin, R., and Shen, W. (2020), “Empirical priors and posterior concentration in a piecewise polynomial sequence model,” arXiv e-prints, arXiv–1712.
  * Padilla et al. (2017) Padilla, O., Sharpnack, J., and Scott, J. (2017), “The DFS fused lasso: Linear-time denoising over general graphs,” The Journal of Machine Learning Research, 18, 6410–6445.
  * R Core Team (2019) R Core Team (2019), R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria.
  * Roualdes (2015) Roualdes, E. A. (2015), “Bayesian trend filtering,” arXiv preprint arXiv:1505.07710.
  * Rudin et al. (1992) Rudin, L. I., Osher, S., and Fatemi, E. (1992), “Nonlinear total variation based noise removal algorithms,” Physica D: nonlinear phenomena, 60, 259–268.
  * Shuman et al. (2013) Shuman, D. I., Narang, S. K., Frossard, P., Ortega, A., and Vandergheynst, P. (2013), “The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains,” IEEE signal processing magazine, 30, 83–98.
  * Song and Cheng (2019) Song, Q. and Cheng, G. (2019), “Bayesian Fusion Estimation via t Shrinkage,” Sankhya A, 1–33.
  * Song and Liang (2017) Song, Q. and Liang, F. (2017), “Nearly optimal Bayesian shrinkage for high dimensional regression,” arXiv preprint arXiv:1712.08964.
  * Tarjan (1972) Tarjan, R. (1972), “Depth-first search and linear graph algorithms,” SIAM journal on computing, 1, 146–160.
  * Tartakovsky et al. (2014) Tartakovsky, A., Nikiforov, I., and Basseville, M. (2014), Sequential analysis: Hypothesis testing and changepoint detection, CRC Press.
  * Tibshirani et al. (2005) Tibshirani, R., Saunders, M., Rosset, S., Zhu, J., and Knight, K. (2005), “Sparsity and smoothness via the fused lasso,” Journal of the Royal Statistical Society: Series B (Statistical Methodology), 67, 91–108.

## Appendix

##  1 Proof of Theorem 1

The proof of the results is based upon the proof of Theorem 2.1 and Corollary
2.1 in Song and Cheng (2019) since both papers are using the $t$-shrinkage
prior. After we perform the DFS on the general graph, our proposed model
becomes the piecewise constant model defined on a chain using $t$ shrinkage
prior as presented in Song and Cheng (2019). Under the assumption that
$\sum_{(i,j)\in E}I(\theta_{i}^{*}\neq\theta_{j}^{*})\leq s$, we have
$\sum_{(i,j)\in E_{C}}I(\theta_{i}^{*}\neq\theta_{j}^{*})\leq 2s$ based on
Lemma 1 in Padilla et al. (2017). Therefore we can directly apply Theorem 2.1
and Corollary 2.1 in Song and Cheng (2019) and obtain the rates. Note that
when performing DFS, the root node can be any of the values in
$\theta_{1},\ldots,\theta_{n}$. Therefore the assumption in Song and Cheng
(2019) are generalized to (C3) for arbitrary pairs of connected nodes
$(\theta_{i},\theta_{j})$ and to (C4) for every $\theta_{i}^{*}$ with
$i=1,\ldots,n$.

##  2 Additional numerical results

Here we present additional numerical results for the Simulation Section in the
main paper.

|
---|---
(a) True signal | (b) $t$-fusion estimates
|
(c) Laplace fusion estimates | (d) $L_{1}$ fusion estimates
Figure 6: Figure showing signal denoising performance for linear chain graphs in case of unevenly spaced signals and error sd $\sigma=0.1.$ 95% credible bands (in green) are also provided for the Bayesian procedures. The red dots correspond to observations. |
---|---
(a) True signal | (b) $t$-fusion estimates
|
(c) Laplace fusion estimates | (d) $L_{1}$ fusion estimates
Figure 7: Figure showing signal denoising performance for linear chain graphs in case of very unevenly spaced signals and error sd $\sigma=0.1.$ 95% credible bands (in green) are also provided for the Bayesian procedures. The red dots correspond to observations. |
---|---
(a) True signal | (b) $t$-fusion estimates
|
(c) Laplace fusion estimates | (d) $L_{1}$ fusion estimates
Figure 8: Figure showing signal denoising performance for linear chain graphs in case of evenly spaced signals and error sd $\sigma=0.3.$ 95% credible bands (in green) are also provided for the Bayesian procedures. The red dots correspond to observations. |
---|---
(a) True signal | (b) $t$-fusion estimates
|
(c) Laplace fusion estimates | (d) $L_{1}$ fusion estimates
Figure 9: Figure showing signal denoising performance for linear chain graphs in case of unevenly spaced signals and error sd $\sigma=0.3.$ 95% credible bands (in green) are also provided for the Bayesian procedures. The red dots correspond to observations. |
---|---
(a) True signal | (b) $t$-fusion estimates
|
(c) Laplace fusion estimates | (d) $L_{1}$ fusion estimates
Figure 10: Figure showing signal denoising performance for linear chain graphs in case of very unevenly spaced signals and error sd $\sigma=0.3.$ 95% credible bands (in green) are also provided for the Bayesian procedures. The red dots correspond to observations. |
---|---
(a) True signal | (b) Noisy signal
|
(c) $t$-fusion estimates | (d) Laplace fusion estimates
(e) $L_{1}$ fusion estimates
Figure 11: Figure showing signal denoising performance for 2-d lattice graphs with weak signal strength. |
---|---
(a) True signal | (b) Noisy signal
|
(c) $t$-fusion estimates | (d) Laplace fusion estimates
(e) $L_{1}$ fusion estimates
Figure 12: Figure showing signal denoising performance for 2-d lattice graphs
with medium signal strength.


# Dualities for absolute zeta functions and multiple gamma functions

Nobushige Kurokawa  Nobushige Kurokawara
Tokyo Institute of Technology<EMAIL_ADDRESS>and  Hiroyuki Ochiai
Hiroyuki Ochiai
Kyushu University<EMAIL_ADDRESS>
###### Abstract.

We study absolute zeta functions from the view point of a canonical
normalization. We introduce the absolute Hurwitz zeta function for the
normalization. In particular, we show that the theory of multiple gamma and
sine functions gives good normalizations in cases related to the Kurokawa
tensor product. In these cases, the functional equation of the absolute zeta
function turns out to be equivalent to the simplicity of the associated non-
classical multiple sine function of negative degree.

###### Key words and phrases:

absolute zeta function, multiple gamma function, multiple sine function

###### 2000 Mathematics Subject Classification:

Primary 11M06

##  1\. Introduction

The absolute zeta function of a scheme $X$ over ${\mathbb{F}_{1}}$ was first
studied by Soulé [S] as a “limit of $p\to 1$” of the (congruence) zeta
function over ${\mathbb{F}_{p}}$: see Kurokawa [K2] and Deitmar [D] also.
Then, Connes and Consani [CC1] [CC2] investigated the absolute zeta function
as the following integral

$\zeta_{X}(s)=\exp\left(\mbox{$\displaystyle\int_{1}^{\infty}\frac{N_{X}(u)}{u^{s+1}\log
u}du$}\right),$

where

$N_{X}(u)=\left|X(\mathbb{F}_{1^{u-1}})\right|$

is a suitably interpolated “counting function.” Here we must pay attention to
the needed normalization for the integral near $u=1$: see [CC1] [CC2] for a
discussion. In [CC1, Theorem 4.13] [CC2, Theorem 4.3] Connes and Consani
calculated $\zeta_{X}(s)$ for Noetherian schemes via the Kurokawa tensor
product of [K1].

Our purpose is to introduce the absolute Hurwitz zeta function

$Z_{X}(w;s)=\frac{1}{\Gamma(w)}\int_{1}^{\infty}\frac{N_{X}(u)}{u^{s+1}(\log
u)^{1-w}}du$

to get the canonical normalization:

$\zeta_{X}(s)=\exp\left(\left.\frac{\partial}{\partial
w}Z_{X}(w;s)\right|_{w=0}\right).$

This normalization is essentially due to Riemann (1859) and it is used in the
theory of multiple gamma and sine function as follows.

For each integer $r\geq 1$, the $r$-ple Hurwitz zeta function $\zeta_{r}(w;x)$
is defined in $\operatorname{Re}(w)>r$ as

$\zeta_{r}(w;x)=\sum_{n=0}^{\infty}{}_{r}H_{n}(n+x)^{-w}$ (1)

where ${}_{r}H_{n}={{n+r-1}\choose{n}}$.

The analytic continuation of $\zeta_{r}(w;x)$ to all $w\in\mathbb{C}$ is
obtained via the integral representation of Riemann

$\displaystyle\zeta_{r}(w;x)$
$\displaystyle=\frac{1}{\Gamma(w)}\int_{0}^{\infty}(1-e^{-t})^{-r}e^{-xt}t^{w-1}dt$
$\displaystyle=\frac{1}{\Gamma(w)}\int_{1}^{\infty}(1-u^{-1})^{-r}u^{x-1}(\log
u)^{w-1}du$

by treating the integral around $u=1$ in the usual way.

Thus, by using such analytic continuation we get the $r$-ple gamma function

$\Gamma_{r}(x)=\exp\left(\left.\frac{\partial}{\partial
w}\zeta_{r}(w;x)\right|_{w=0}\right)$

and the $r$-ple sine function

$S_{r}(x)=\Gamma_{r}(x)^{-1}\Gamma_{r}(r-x)^{(-1)^{r}}.$

We refer to Barnes [B] (1904) and Kurokawa-Koyama [KK] (2003) for details,
where more general multiple gamma functions and multiple sine functions were
treated respectively.

We report three results in this introduction. First, for a function
$N:(1,\infty)\rightarrow\mathbb{C}$ we use

$Z_{N}(w;s)=\frac{1}{\Gamma(w)}\int_{1}^{\infty}N(u)u^{-s-1}(\log u)^{w-1}du$

and

$\zeta_{N}(s)=\exp\left(\left.\frac{\partial}{\partial
w}Z_{N}(w;s)\right|_{w=0}\right)$

also.

######  Theorem A.

Let $N(u)=\displaystyle\sum_{\alpha}m(\alpha)u^{\alpha}$ be a finite sum.
Then:

  * (1)

$Z_{N}(w;s)=\displaystyle\sum_{\alpha}m(\alpha)(s-\alpha)^{-w}$.

  * (2)

$\zeta_{N}(s)=\displaystyle\prod_{\alpha}(s-\alpha)^{-m(\alpha)}$.

This result is applicable to calculate many examples (see [K2]) of absolute
zeta functions under our canonical normalization. We note two simple examples.

######  Example 1.

Let $X=\operatorname{Spec}{\mathbb{F}_{1}}$. Then

$\displaystyle N_{X}(u)=1,$ $\displaystyle Z_{X}(w;s)=s^{-w},$
$\displaystyle\zeta_{X}(s)=1/s.$

######  Example 2.

Let $X=\mathbb{SL}(2)$. Then

$\displaystyle N_{X}(u)=u^{3}-u,$ $\displaystyle
Z_{X}(w;s)=(s-3)^{-w}(s-1)^{-w},$ $\displaystyle\zeta_{X}(s)=(s-1)/(s-3).$

Now the following result shows a functoriality.

######  Theorem B.

  * (1)

For $N_{1},N_{2}:(1,\infty)\rightarrow\mathbb{C}$ let

$(N_{1}\oplus N_{2})(u)=N_{1}(u)+N_{2}(u).$

Then

$Z_{N_{1}\oplus N_{2}}(w;s)=Z_{N_{1}}(w;s)+Z_{N_{2}}(w;s)$

and

$\zeta_{N_{1}\oplus N_{2}}(s)=\zeta_{N_{1}}(s)\zeta_{N_{2}}(s).$

  * (2)

Let

$N_{i}(u)=\sum_{\alpha_{i}}m_{i}(\alpha_{i})u^{\alpha_{i}}$

for $i=1,2$. Suppose that both are finite sums. Put

$(N_{1}\otimes N_{2})(u)=N_{1}(u)N_{2}(u).$

Then

$\displaystyle Z_{N_{1}\otimes N_{2}}(w;s)$
$\displaystyle=\sum_{\alpha_{1},,\alpha_{2}}m_{1}(\alpha_{1})m_{2}(\alpha_{2})(s-(\alpha_{1}+\alpha_{2}))^{-w}$

and

$\zeta_{N_{1}\otimes
N_{2}}(s)=\prod_{\alpha_{1},\alpha_{2}}(s-(\alpha_{1}+\alpha_{2}))^{-m_{1}(\alpha_{2})m_{2}(\alpha_{2})}.$

This tensor product is essentially the Kurokawa tensor product originated in
[K1] (see [M], [CC1] and [CC2]) when $\alpha_{j}$’s are real. We remark that
for general $N_{j}$’s (“infinite sums” or “generalized functions”) we must
resolve various difficulties.

For the next result we notice that our construction of $\zeta_{r}(w;x)$,
$\Gamma_{r}(x)$ and $S_{r}(x)$ is valid for negative $r$ also (see the later
explanation).

######  Theorem C.

Let $r$ be a positive integer. Then

  * (1)

$\displaystyle Z_{\mathbb{G}_{\rm m}^{\otimes r}}(w;s)=\zeta_{-r}(w;s-r)$.

  * (2)

$\displaystyle\zeta_{\mathbb{G}_{\rm m}^{\otimes r}}(s)=\Gamma_{-r}(s-r)$
$\displaystyle=\prod_{j=1}^{r}(s-j)^{(-1)^{r-j-1}{r\choose j}}$
$\displaystyle=\left((1-1/s)^{\otimes r}\right)^{-1}$,
where $\otimes r$ is the Kurokawa tensor product.

  * (3)

We have the functional equation

$\zeta_{\mathbb{G}_{\rm m}^{\otimes r}}(s)=\zeta_{\mathbb{G}_{\rm m}^{\otimes
r}}(r-s)^{(-1)^{r}},$

which is equivalent to $S_{-r}(x)=1$.

Our result would suggest that

$\zeta_{\mathbb{G}_{\rm m}^{\otimes r}}(s)=\Gamma_{-r}(s-r)$

holds for $r<0$ also with the functional equation $s\leftrightarrow-r-s$. For
example

$\zeta_{\mathbb{G}_{\rm
m}^{\otimes-1}}(s)=\Gamma_{1}(s+1)=\frac{\Gamma(s+1)}{\sqrt{2\pi}}$

and the functional equation $s\leftrightarrow 1-s$ is the reflection formula
of Euler:

$\Gamma_{1}(s+1)\Gamma(2-s)=S_{1}(s+1)^{-1}=-\frac{1}{2\sin(\pi s)}.$

We remark that Manin [M, §1.7] indicated an idea to consider the gamma
function as the zeta function of the “dual infinite dimensional projective
space over ${\mathbb{F}_{1}}$.”

##  2\. Multiple gamma functions and multiple sine functions

We recall the construction of the multiple Hurwitz zeta function:

$\displaystyle\zeta_{r}(w;x)$ $\displaystyle=\sum_{n=0}^{\infty}{n+r-1\choose
n}(n+x)^{-w}$
$\displaystyle=\frac{1}{\Gamma(w)}\int_{0}^{\infty}(1-e^{-t})^{-r}e^{-xt}t^{w-1}dt$
$\displaystyle=\frac{1}{\Gamma(w)}\int_{1}^{\infty}(1-u^{-1})^{-r}u^{-x-1}(\log
u)^{w-1}du.$

This definition is valid for any $r\in\mathbb{R}$ with sufficiently large
$\operatorname{Re}(x)$ and $\operatorname{Re}(w)$, so we have the analytic
continuation to all $w\in\mathbb{C}$ via the usual method. Thus, we get

$\Gamma_{r}(x)=\exp\left(\left.\frac{\partial}{\partial
w}\zeta_{r}(w;x)\right|_{w=0}\right)$

and

$S_{r}(x)=\Gamma_{r}(x)^{-1}\Gamma_{r}(r-x)^{(-1)^{r}}$

for any $r\in\mathbb{R}$ (or $r\in\mathbb{Z}$ at least without ambiguity of
the meaning of $(-1)^{r}$). For readers interested in the theory of $r<0$, we
refer to [KO].

######  Theorem 1.

Let $r$ be a negative integer. Then

  * (1)

$\displaystyle\Gamma_{r}(x)=\prod_{n=0}^{-r}(x+n)^{(-1)^{n+1}{-r\choose n}}$.

  * (2)

$S_{r}(x)=1$.

###### Proof.

We have

$\displaystyle\zeta_{r}(w;x)$ $\displaystyle=\sum_{n=0}^{\infty}{n+r-1\choose
n}(n+x)^{-w}$ $\displaystyle=\sum_{n=0}^{\infty}(-1)^{n}{-r\choose
n}(n+x)^{-w}.$

Hence

$\displaystyle\Gamma_{r}(x)$
$\displaystyle=\exp\left(\sum_{n=0}^{-r}(-1)^{n+1}{-r\choose
n}\log(n+x)\right)$ $\displaystyle=\prod_{n=0}^{-r}(n+x)^{(-1)^{n+1}{-r\choose
n}}.$

Next,

$\displaystyle S_{r}(x)=\Gamma_{r}(x)^{-1}\Gamma_{r}(r-x)^{(-1)^{r}}$
$\displaystyle=\prod_{n=0}^{-r}(n+x)^{(-1)^{n}{-r\choose
n}}\times\prod_{n=0}^{-r}(n+r-x)^{(-1)^{n-r+1}{-r\choose n}}$
$\displaystyle=\prod_{n=0}^{-r}(n+x)^{(-1)^{n}{-r\choose n}}$
$\displaystyle\qquad\times\prod_{n=0}^{-r}((-r-n)+x)^{(-1)^{(-r-n)+1}{-r\choose
n}},$

where we used

$\sum_{n=0}^{-r}(-1)^{n}{-r\choose n}=0.$

Hence

$\displaystyle S_{r}(x)$
$\displaystyle=\prod_{n=0}^{-r}(n+x)^{(-1)^{n}{-r\choose
n}}\times\prod_{n=0}^{-r}(n+x)^{(-1)^{n+1}{-r\choose n}}$ $\displaystyle=1.$

∎

This result can be generalized to the multi-period case
$\underline{\omega}=(\omega_{1},\ldots,\omega_{r})$ with
$\omega_{1},\ldots,\omega_{r}>0$ as follows, where the above case is contained
as $\underline{\omega}=(1,\ldots,1)$. Put

$\displaystyle\zeta_{-r}(w;x,\underline{\omega})$
$\displaystyle\qquad=\sum_{1\leq i_{1}<\cdots<i_{k}\leq
r}(-1)^{k}(x+\omega_{i_{1}}+\cdots+\omega_{i_{k}})^{-w},$
$\displaystyle\Gamma_{-r}(w,\underline{\omega})=\exp\left(\left.\frac{\partial}{\partial
w}\zeta_{-r}(w;x,\underline{\omega})\right|_{w=0}\right),$ and $\displaystyle
S_{-r}(x,\underline{\omega})=\Gamma_{-r}(x,\underline{\omega})^{-1}$
$\displaystyle\quad\qquad\times\Gamma_{-r}(-(\omega_{1}+\cdots+\omega_{r})-x,\underline{\omega})^{(-1)^{r}}.$

Then we have (see [KO] for more generalizations also)

$\displaystyle\zeta_{-r}(w;x,\underline{\omega})$
$\displaystyle\quad=\frac{1}{\Gamma(w)}\int_{0}^{w}(1-e^{-t\omega_{1}})\cdots(1-e^{-t\omega_{r}})e^{-xt}t^{w-1}dt,$
$\displaystyle\Gamma_{-r}(x,\underline{\omega})=\prod_{1\leq
i_{1}<\cdots<i_{k}\leq r}(x+\omega_{i_{1}}+\cdots+\omega_{i_{k}})^{(-1)^{k}},$
and $\displaystyle S_{-r}(x,\underline{\omega})=1.$

For example, we get

$\zeta_{\mathbb{SL}(2)}(s)=\Gamma_{-1}(s-3,2)=\frac{s-1}{s-3}.$

More generally:

$\displaystyle\zeta_{\mathbb{SL}(r)}(s)=\Gamma_{-(r-1)}(s-(r^{2}-1),(2,3,\cdots,r))$
and
$\displaystyle\zeta_{\mathbb{GL}(r)}(s)=\Gamma_{-r}(s-r^{2},(1,2,3,\cdots,r)),$

where $\left\\{\begin{array}[]{l}r-1=\operatorname{rank}\mathbb{SL}(r)\\\
r^{2}-1=\dim\mathbb{SL}(r)\end{array}\right.$ and
$\left\\{\begin{array}[]{l}r=\operatorname{rank}\mathbb{GL}(r)\\\
r^{2}=\dim\mathbb{GL}(r).\end{array}\right.$ We obtain the functional
equations

$\displaystyle\zeta_{\mathbb{SL}(r)}(s)=\zeta_{\mathbb{SL}(r)}(r(3r-1)/2-1-s)^{(-1)^{r-1}},$
and
$\displaystyle\zeta_{\mathbb{GL}(r)}(s)=\zeta_{\mathbb{GL}(r)}(r(3r-1)/2-s)^{(-1)^{r}}$

from the triviality of the multiple sine function of negative order exactly
similar to Theorem C.

######  Theorem 2.

Let $r$ be a negative real number. Then:

  * (1)

$\zeta_{r}(m;x)=0$ for each integer $m$ satisfying $r<m\leq 0$.

  * (2)

$\displaystyle\Gamma_{r}(x)=\exp\left(\int_{1}^{\infty}(1-u^{-1})^{-r}u^{-x-1}(\log
u)^{-1}du\right)$ for $\operatorname{Re}(x)>0$.

######  Example 3.

$\zeta_{-3}(w;x)=x^{-w}-3(x+1)^{-w}+3(x+2)^{-w}-(x+3)^{-w}$

and

$\zeta_{-3}(0;x)=\zeta_{-3}(-1;x)=\zeta_{-3}(-2;x)=0.$

Notice that $\zeta_{-3}(-3;x)=-6$. (In general $\zeta_{-m}(-m;x)=(-1)^{m}m!$
for integers $m\geq 0$.

######  Example 4.

$\zeta_{-\frac{1}{2}}(w;x)=x^{-w}-\sum_{n=1}^{\infty}\frac{{2n\choose
n}}{(2n-1)4^{n}}(n+x)^{-w}$

and

$\zeta_{-\frac{1}{2}}(0;x)=1-\sum_{n=1}^{\infty}\frac{{2n\choose
n}}{(2n-1)4^{n}}=0,$

that is

$\sum_{n=1}^{\infty}\frac{{2n\choose n}}{(2n-1)4^{n}}=1.$

###### Proof.

The fact (1) follows from the integral representation

$\zeta_{r}(w;x)=\frac{1}{\Gamma(w)}\int_{1}^{\infty}(1-u^{-1})^{-r}u^{-x-1}(\log
u)^{w-1}du,$

since this integral converges for $\operatorname{Re}(w)>-r$ when
$\operatorname{Re}(x)>0$, and $1/\Gamma(w)$ has zeros at $w=0,-1,\ldots,r+1$.
Similarly, (2) is seen by looking at $w=0$. ∎

##  3\. Proof of Theorem A

For a function $N:(1,\infty)\rightarrow\mathbb{C}$ we defined

$Z_{N}(w;s)=\frac{1}{\Gamma(w)}\int_{1}^{\infty}N(u)u^{-s-1}(\log u)^{w-1}du$

and

$\zeta_{N}(s)=\exp\left(\left.\frac{\partial}{\partial
w}Z_{N}(w;s)\right|_{w=0}\right).$

We calculate these functions in the case of a finite sum

$N(u)=\sum_{\alpha}m(\alpha)u^{\alpha}.$

It is sufficient to calculate the following monomial case.

######  Lemma.

Let $N(u)=u^{\alpha}$, then

$Z_{N}(w;s)=(s-\alpha)^{-w}$

and

$\zeta_{N}(s)=\frac{1}{s-\alpha}.$

###### Proof.

$\displaystyle Z_{N}(w;s)$
$\displaystyle=\frac{1}{\Gamma(w)}\int_{1}^{\infty}u^{\alpha-s-1}(\log
u)^{w-1}du$
$\displaystyle=\frac{1}{\Gamma(w)}\int_{0}^{\infty}e^{-(s-\alpha)t}t^{w-1}dt$
$\displaystyle=(s-\alpha)^{-w}.$

Hence

$\left.\frac{\partial}{\partial w}Z_{N}(w,s)\right|_{w=0}=-\log(s-\alpha)$

and

$\zeta_{N}(s)=\frac{1}{s-\alpha}.$

∎

##  4\. Proof of Theorem B

(1) Since

$\displaystyle Z_{N_{1}\oplus N_{2}}(w;s)$
$\displaystyle=\frac{1}{\Gamma(w)}\int_{1}^{\infty}(N_{1}\oplus
N_{2})(u)u^{-s-1}(\log u)^{w-1}du$
$\displaystyle=\frac{1}{\Gamma(w)}\int_{1}^{\infty}(N_{1}(u)+N_{2}(u))u^{-s-1}(\log
u)^{w-1}du$ $\displaystyle=Z_{N_{1}}(w;s)+Z_{N_{2}}(w;s),$

we have

$\zeta_{N_{1}\oplus N_{2}}(s)=\zeta_{N_{1}}(s)\zeta_{N_{2}}(s).$

(2) From

$\displaystyle(N_{1}\oplus N_{2})(u)$ $\displaystyle=N_{1}(u)N_{2}(u)$
$\displaystyle=(\sum_{\alpha_{1}}m_{1}(\alpha_{1})u^{\alpha_{1}})(\sum_{\alpha_{2}}m_{2}(\alpha_{2})u^{\alpha_{2}})$
$\displaystyle=\sum_{\alpha_{1},\alpha_{2}}m_{1}(\alpha_{1})m_{2}(\alpha_{2})u^{\alpha_{1}+\alpha_{2}},$

we have

$\displaystyle Z_{N_{1}\otimes N_{2}}(w;s)$
$\displaystyle=\frac{1}{\Gamma(w)}\int_{1}^{\infty}(N_{1}\otimes
N_{2})(u)u^{-s-1}(\log u)^{w-1}du$
$\displaystyle=\frac{1}{\Gamma(w)}\int_{1}^{\infty}\Big{(}\sum_{\alpha_{1},\alpha_{2}}m_{1}(\alpha_{1})m_{2}(\alpha_{2})u^{\alpha_{1}+\alpha_{2}}\Big{)}$
$\displaystyle\hskip 113.81102pt\times u^{-s-1}(\log u)^{w-1}du$
$\displaystyle=\sum_{\alpha_{1},\alpha_{2}}m_{1}(\alpha_{1})m_{2}(\alpha_{2})(s-(\alpha_{1}+\alpha_{2}))^{-w}.$

Hence

$\zeta_{N_{1}\otimes
N_{2}}(s)=\prod_{\alpha_{1},\alpha_{2}}(s-(\alpha_{1}+\alpha_{2}))^{-m_{1}(\alpha_{1})m_{2}(\alpha_{2})}.\qed$

##  5\. Absolute zeta functions

######  Theorem 3.

Let $r$ be a positive integer. Then

  * (1)

$\displaystyle\zeta_{\mathbb{G}_{\rm m}^{\otimes r}}(s)=\Gamma_{-r}(s-r)$.

  * (2)

$\displaystyle\zeta_{\mathbb{G}_{\rm m}^{\otimes
r}}(s)=\exp\Big{(}\int_{1}^{\infty}N_{\mathbb{G}_{\rm m}^{\otimes
r}}(u)u^{-s-1}(\log u)^{-1}du\Big{)}.$

###### Proof.

(1) Since

$N_{\mathbb{G}_{\rm m}^{\otimes r}}(u)=(u-1)^{r},$

we have

$\displaystyle Z_{\mathbb{G}_{\rm m}^{\otimes r}}(w;s)$
$\displaystyle=\frac{1}{\Gamma(w)}\int_{1}^{\infty}(u-1)^{r}u^{-s-1}(\log
u)^{w-1}du$
$\displaystyle=\frac{1}{\Gamma(w)}\int_{1}^{\infty}(1-u^{-1})^{r}u^{-s+r-1}(\log
u)^{w-1}du$ $\displaystyle=\zeta_{-r}(w;s-r).$

Thus,

$\zeta_{\mathbb{G}_{\rm m}^{\otimes r}}(s)=\Gamma_{-r}(s-r).$

(2) This follows from (1) and Theorem 2(2). ∎

We notice that Theorem 1 and Theorem 3(1) imply Theorem C(1)(2).

##  6\. Functional equations

######  Theorem 4.

Let $r$ be a positive integer. Then

$\frac{\zeta_{\mathbb{G}_{\rm m}^{\otimes r}}(s)}{\zeta_{\mathbb{G}_{\rm
m}^{\otimes r}}(r-s)^{(-1)^{r}}}=S_{-r}(s-r)^{-1}.$

###### Proof.

From Theorem 3(1), we have

$\zeta_{\mathbb{G}_{\rm m}^{\otimes r}}(s)=\Gamma_{-r}(s-r)$

and

$\zeta_{\mathbb{G}_{\rm m}^{\otimes r}}(r-s)=\Gamma_{-r}(-s).$

Hence,

$\displaystyle\zeta_{\mathbb{G}_{\rm m}^{\otimes r}}(s)\zeta_{\mathbb{G}_{\rm
m}^{\otimes r}}(r-s)^{(-1)^{r+1}}$
$\displaystyle=\Gamma_{-r}(s-r)\Gamma_{-r}(-s)^{(-1)^{r+1}}$
$\displaystyle=S_{-r}(s-r)^{-1}.$

∎

We remark that we have the functional equation

$\zeta_{\mathbb{G}_{\rm m}^{\otimes r}}(s)=\zeta_{\mathbb{G}_{\rm m}^{\otimes
r}}(r-s)^{(-1)^{r}}$

from Theorem 1(2) and we know that it is equivalent to $S_{-r}(x)=1$. Thus we
have Theorem C(3).

## References

  * [B] E.W. Barnes, On the theory of the multiple gamma functions. Trans. Cambridge Philos. Soc. 19 (1904) 374–425.
  * [CC1] A. Connes and C. Consani, Schemes over ${\mathbb{F}}_{1}$ and zeta functions. Compositio Mathematica 146 (2010) 1383–1415.
  * [CC2] A. Connes and C. Consani, Characteristic one, entropy and the absolute point. In ”Noncommutative Geometry, Arithmetic, and Related Topics, Proceedings of the JAMI Conference 2009”, Johns Hopkins University Press (2011) 75–140.
  * [D] A. Deitmar, Remarks on zeta functions and $K$-theory over ${\mathbf{F}}_{1}$, Proc. Japan Acad. Ser. A Math. Sci. 82 (2006) 141–146.
  * [K1] N. Kurokawa, Multiple zeta functions: an example. In “Zeta Functions in Geometry” (Tokyo 1990), Adv. Stud. Pure Math. 21, Kinokuniya, Tokyo, 1992, 219–226.
  * [K2] N. Kurokawa, Zeta functions over ${\mathbb{F}}_{1}$, Proc. Japan Acad. Ser. A Math. Sci. 81 (2005) 180–184.
  * [KK] N. Kurokawa and S. Koyama, Multiple sine functions. Forum Math. 15 (2003) 839–876.
  * [KO] N. Kurokawa and H. Ochiai, Multiple gamma functions of negative order, 2013, preprint.
  * [M] Y. Manin, Lectures on zeta functions and motives (according to Deninger and Kurokawa), Astérisque 228 (1995), 121–163.
  * [S] C. Soulé, Les variétés sur le corps à un élément. Mosc. Math. J. 4 (2004) 217–244.