# Entropy and its conservation in expanding Universe

Sinya Aoki<EMAIL_ADDRESS>Center for Gravitational Physics and
Quantum Information, Yukawa Institute for Theoretical Physics, Kyoto
University, Kitashirakawa Oiwakecho, Sakyo-ku, Kyoto 606-8502, Japan  Kiyoharu
Kawana<EMAIL_ADDRESS>School of Physics, KIAS, Seoul 02455, Korea

###### Abstract

We investigate properties of the conserved charge in general relativity,
recently proposed by one of the present authors with his collaborators, in the
inflation era, the matter dominated era and the radiation dominated era of the
expanding Universe. We show that the conserved charge in the inflation era
becomes the Bekenstein-Hawking entropy for de Sitter space, and it becomes the
matter entropy and the radiation entropy in the matter and radiation dominated
eras, respectively, while the charge itself is always conserved. These
properties are qualitatively confirmed by a numerical analysis of a model with
a scalar field and radiations. Results in this paper provide more evidences on
the interpretation that the conserved charge in general relativity corresponds
to entropy.

††preprint: YITP-22-115

##  I Introduction

Entropy is one of the most fundamental quantities in nature and plays an
important role to connect microscopic and macroscopic physics. For example, in
the micro-canonical statistical mechanics, it is given by
$S(E,V,N)=\ln\Omega(E,V,N)$ where $\Omega(E,V,N)$ is a number of micro states
in a system having definite $E,V,N$, suggesting that entropy is defined from
some microscopic degree of freedom (dof). In thermodynamic limit
$V\rightarrow\infty$, this entropy correctly reproduces various properties of
macroscopic thermodynamic entropy.

Another famous example is of course a black hole Bekenstein (1973); Hawking
(1975); Bekenstein (1974), whose entropy is given by

$S_{H}={A_{H}\over 4G_{N}},$ (1)

where $A_{H}$ is an area of the black hole horizon and $G_{N}$ is the Newton
constant. This expression suggests that the microscopic dof of a black hole
may be localized at its surface just like (mem)brane objects Strominger and
Vafa (1996); Maldacena et al. (1999); Thorne et al. (1986). Although its
microscopic origin is still not fully understood yet, this formula seems a key
to quantum gravity.

In recent works Aoki et al. (2021a, b); Aoki and Onogi (2022); Aoki (2022),
one of the present authors and his collaborators proposed a new definition of
a conserved charge $Q[\xi]$ for a generic metric field $g_{\mu\nu}$, even when
a system has no Killing vectors. (See Section II for details.) In particular,
it was found for several cases that the conserved charge constructed by a
particular vector field $\xi^{\mu}=-\beta\delta^{\mu}_{0}$ corresponds to
entropy of a system with $\beta$ proportional to an inverse temperature. In
fact, they explicitly showed that the new conserved charge correctly
reproduces the Bekenstein-Hawking entropy Eq. (1) for several black hole
systems such as the Schwarzschild black hole and the BTZ blackhole Aoki et al.
(2021b). Interestingly, their definition of “entropy” contains the volume
integral of the energy-momentum tensor and the entropy density is localized at
the singularity not at the surface. This may also open a new insight for the
black hole information paradox Hawking (1976, 1975); Raju (2022).

One of trivial but important aspects of $Q[\zeta]$ is that it is always
conserved no matter what happens during the dynamical evolution of a system.
In particular, it is quite common in cosmology that dominated matter content
of the Universe is changing from inflation era to radiation (matter) dominated
era. If the interpretation of $Q[\zeta]$ as entropy is actually true, it
should represent some “entropy” for each different eras. A purpose of this
paper is to clarify the physical meaning of $Q[\zeta]$ in such dynamical
transition processes in the expanding Universe. After a brief review on the
new proposal for the conserved charge in general relativity in the next
section, we first analytically study a few transition processes (inflation
$\to$ matter $\to$ radiation) and find that $Q[\zeta]$ can be actually
interpreted as entropy for each epoch in the following way in Sec. III: During
the inflation era, $Q[\zeta]$ reproduces the Bekenstein-Hawking entropy for de
Sitter space Gibbons and Hawking (1977) by taking $\beta_{0}$ (an initial
value of $\beta$) equal to the inverse de Sitter temperature. On the other
hand, $Q[\zeta]$ simply corresponds to the total matter entropy as
$Q[\zeta]=E/T_{M}$ with a constant matter temperature (energy) $T_{M}$ in the
matter dominated era, while $Q[\zeta]$ agrees with the radiation entropy in
the radiation dominated era taking ${d\over
d-1}\beta^{-1}(x^{0}):=T_{R}(x^{0})$ as radiation temperature. In particular,
it is notable that the radiation energy density satisfies the Stefan-Boltzmann
law $\rho_{R}\propto\beta^{-d}$ without assuming thermal equilibrium of
radiations. While $\beta$ is continuous, the temperature is discontinuous at
transition points (from inflation to matter and from matter to radiation).

After the analytical studies in Sec. IV, we numerically solve a dynamics of a
scalar field theory, which emulate the above dynamics, to confirm the
conservation of $Q[\zeta]$ and above interpretations. In Sec. V, we conclude
our paper.

##  II A matter conserved charge

In this section, we briefly review a new definition on a matter conserved
charge in general relativity(GR) proposed in Refs. Aoki et al. (2021a, b);
Aoki and Onogi (2022).

We first define a matter energy in a $d$-dimensional curve spacetime as

$\displaystyle E(x^{0})$ $\displaystyle:=$
$\displaystyle\int_{\Sigma(x^{0})}[d^{d-1}x]_{\mu}\,T^{\mu}{}_{\nu}(x)\xi_{E}^{\nu}(x)~{},$
(2)

where $T^{\mu}{}_{\nu}$ is a (matter) energy momentum tensor (EMT),
$\Sigma(x^{0})$ is a constant $x^{0}$ (space-like) hyper-surface with a $d-1$
dimensional hyper-surafce element $[d^{d-1}x]_{\mu}$, and
$\xi_{E}^{\mu}(x)=-\delta_{0}^{\mu}$ is a generator of the time translation.
Here we take a coordinate system which satisfies $g_{0j}=0$ for
$j=1,2,\cdots,d-1$ without loss of generality. This energy is not guaranteed
to be conserved in general, except special cases such that $\xi_{E}$ becomes a
Killing vector Aoki et al. (2021a, b); Aoki and Onogi (2022).

It has been pointed out Aoki et al. (2021b), however, that there always exists
a conserved charge in GR, given by

$\displaystyle Q(x^{0})$ $\displaystyle:=$
$\displaystyle\int_{\Sigma(x^{0})}[d^{d-1}x]_{\mu}\,T^{\mu}{}_{\nu}(x)\beta(x)n^{\nu}(x)~{},$
(3)

where the scalar function $\beta$ is chosen to satisfy
$\nabla_{\mu}(T^{\mu}{}_{\nu}\beta n^{\nu})=0$, and $n^{\mu}$ is a time-like
vector. While an appropriate choice of $n^{\mu}$ is needed to be clarified for
a general case in future studies, it is reasonable to take
$n^{\mu}=\xi^{\mu}_{E}$ in this paper as we will see. The conservation
condition $\nabla_{\mu}(T^{\mu}{}_{\nu}\beta n^{\nu})=0$ implies

$\displaystyle
Q(x^{0}_{a})-Q(x_{b}^{0})=\int_{x^{0}_{a}}^{x^{0}_{b}}dx^{0}\,\int[d^{d-2}x]_{i}T^{i}{}_{\nu}(x)\beta(x)n^{\nu}(x),~{}~{}~{}~{}$
(4)

so that $Q(x^{0})$ is $x^{0}$-independent if contributions at spatial
boundaries in the right-hand side vanish. An absence of such boundary
contributions is obviously true if $T^{\mu}{}_{\nu}=0$ at the boundary. In the
presence of non-zero EMT at spatial boundaries,
$[d^{d-2}x]_{i}T^{i}{}_{\nu}=0$ is always satisfied: If
$[d^{d-2}x]_{i}T^{i}{}_{\nu}\not=0$, the spacetime must develop in the normal
direction to the boundary according to the Einstein equation. It is impossible
since there is no spacetime beyond the boundary.

Using $\nabla_{\mu}T^{\mu}{}_{\nu}=0$, the condition for $\beta$ becomes

$\displaystyle
T^{0}{}_{0}\partial_{0}\beta+T^{j}{}_{0}\partial_{j}\beta+T^{\mu}{}_{\nu}\Gamma^{\nu}_{\mu
0}\beta=0~{},$ (5)

which is a linear partial differential equation, whose solution can be easily
obtained once an initial condition of $\beta$ is given at some $x^{0}$ Aoki et
al. (2021b). This charge can be regarded as a conserved charge from the
Noether’s 1st theorem for a global symmetry of the matter action in GR, called
a hidden (global) matter symmetry Aoki (2022).

In Ref. Aoki et al. (2021b), it has been pointed out that $Q$ can be
interpreted as entropy with $\beta(x)$ identified as a (time-dependent local)
inverse temperature, since they satisfy the first law of thermodynamics for
some cases including the expanding Universe. Furthermore, it was shown that
$Q$ also reproduces the Bekenstein-Hawking entropy Eq. (1) for Schwarzschild
blackhole and BTZ blackhole Aoki et al. (2021b).

In the next section, we will investigate the conserved charge $Q$ for the
expanding Universe in more detail, in order to add more evidences for our
interpretation that $Q$ is entropy and $\beta(x)$ is an inverse temperature.

##  III Conservation in expanding Universe

###  III.1 Setup

We consider the expanding Universe described by

$\displaystyle ds^{2}=-(dx^{0})^{2}+a^{2}(x^{0})\delta_{ij}dx^{i}dx^{j}~{},$
(6)

and the EMT for a perfect fluid as

$\displaystyle T^{0}{}_{0}=-\rho(x^{0})~{},\quad
T^{i}{}_{j}=\delta^{i}_{j}P(x^{0}):=\delta^{i}_{j}\omega(x^{0})\rho(x^{0})~{},~{}~{}~{}$
(7)

whose covariant conservation implies

$\displaystyle\dot{\rho}$ $\displaystyle=$
$\displaystyle-(d-1)(\rho+P)H=-(d-1)(1+\omega)\rho H~{},$ (8)

where $\omega(x^{0})$ is an equation of state, and $H=\dot{a}/a$. Throughout
this letter, we take $c=\hbar=k_{B}=1$. In this setup, the condition (5)
becomes

$\displaystyle{\dot{\beta}\over\beta}$ $\displaystyle=$
$\displaystyle(d-1)H{P\over\rho}~{}.$ (9)

Energy and charge are evaluated as

$\displaystyle E(x^{0})$ $\displaystyle=$
$\displaystyle-\int_{\Sigma_{x^{0}}}d^{d-1}x\,T^{0}{}_{0}=V_{d-1}[\rho
a^{d-1}](x^{0}),$ (10) $\displaystyle Q(x^{0})$ $\displaystyle=$
$\displaystyle-\int_{\Sigma_{x^{0}}}d^{d-1}x\,T^{0}{}_{0}\beta=V_{d-1}[\rho
a^{d-1}\beta](x^{0}),~{}~{}~{}$ (11)

where $V_{d-1}=\int d^{d-1}x$ is a co-moving volume. While the energy is not
conserved ($\dot{E}\not=0$), it is easy to see from eq. (8) and eq. (9) that
the charge is indeed conserved ($\dot{Q}=0$ ). Furthermore, using $Q=E\beta$
and eq. (9), we have

$\displaystyle{dQ\over dx^{0}}$ $\displaystyle=$ $\displaystyle{dE\over
dx^{0}}\beta+E{d\beta\over dx^{0}}=\left({dE\over dx^{0}}+P{dV\over
dx^{0}}\right)\beta,$ (12)

where $V:=V_{d-1}a^{d-1}$ is a space volume. The 1st law of thermodynamics,
$TdS=dE+PdV$, suggests that $Q$ and $\beta$ may be regarded as entropy and a
(time-dependent) inverse temperature, respectively Aoki et al. (2021b).

###  III.2 Solution for a constant $\omega$

When $\omega=$ constant, Einstein equation reads

$\displaystyle\dot{H}$ $\displaystyle=$ $\displaystyle-{(d-1)(1+\omega)\over
2}H^{2}~{},$ (13)

whose solution is given by

$\displaystyle H(x^{0})$ $\displaystyle=$ $\displaystyle{H_{0}\over
1+C_{0}H_{0}x^{0}}~{},\quad C_{0}:={(d-1)(1+\omega)\over 2}~{}.$ (14)

We then obtain

$\displaystyle a(x^{0})$ $\displaystyle=$ $\displaystyle
a_{0}(1+C_{0}H_{0}x^{0})^{1/C_{0}}~{},$ (15) $\displaystyle\rho(x^{0})$
$\displaystyle=$
$\displaystyle{\rho_{0}\over(1+C_{0}H_{0}x^{0})^{2}}=\rho_{0}\left(\frac{a_{0}}{a(x^{0})}\right)^{(d-1)(1+\omega)}~{},$
(16) $\displaystyle\beta(x^{0})$ $\displaystyle=$
$\displaystyle\beta_{0}(1+C_{0}H_{0}x^{0})^{2\omega\over
1+\omega}=\beta_{0}\left(\frac{a(x^{0})}{a_{0}}\right)^{(d-1)\omega}~{},$ (17)

where $H_{0}$, $a_{0}$, $\rho_{0}$ and $\beta_{0}$ are initial values. Since
$Q$ is a conserved charge, it is determined by its initial value as

$\displaystyle Q$ $\displaystyle=$ $\displaystyle
V_{d-1}\rho_{0}a_{0}^{d-1}\beta_{0}~{}.$ (18)

In the $\omega+1\to 0~{}(C_{0}\rightarrow 0)$ limit, we have $H(x^{0})=H_{0}$,
$\rho(x^{0})=\rho_{0}$,

$\displaystyle
a(x^{0})=a_{0}e^{H_{0}x^{0}},\quad\beta(x^{0})=\beta_{0}e^{-(d-1)H_{0}x^{0}}.$
(19)

###  III.3 Model of expanding Universe

We consider a simplified model of the expanding Universe, whose time
evaluation is given as follows: (1) Inflation era at $0\leq x^{0}<t_{M}$ where
$P=-\rho$ ($\omega=-1$), (2) Matter dominated era at $t_{M}\leq x^{0}<t_{R}$
where $P=0$ ($\omega=0$), and (3) Radiation dominated era at $t_{R}\leq t$
where $P=\rho/(d-1)$ ($\omega=1/(d-1)$). This history captures the essence of
inflaton dynamics after the inflation. We determine the time dependences of
$a,H,\rho,\beta,E$ and $Q$ at each epochs, and fix initial constants $H_{0}$,
$a_{0}$, $\rho_{0}$ and $\beta_{0}$ by requiring that $Q$ represents the de-
Sitter entropy during the inflation.

####  III.3.1 Inflation era

Solutions with $\omega=-1$ read

$\displaystyle H(x^{0})$ $\displaystyle=$ $\displaystyle H_{I},\
\rho(x^{0})=\rho_{I}={(d-1)(d-2)\over 16\pi G_{N}}H_{I}^{2},$ (20)
$\displaystyle a(x^{0})$ $\displaystyle=$ $\displaystyle a_{0}e^{H_{I}x^{0}},\
\beta(x^{0})=\beta_{0}e^{-(d-1)H_{I}x^{0}},$ (21) $\displaystyle E(x^{0})$
$\displaystyle=$ $\displaystyle
V_{d-1}\rho_{I}a_{0}^{d-1}e^{(d-1)H_{I}x^{0}},Q=a_{0}^{d-1}V_{d-1}\rho_{I}\beta_{0}.~{}~{}~{}~{}~{}~{}$
(22)

While the energy is exponentially increasing, $Q$ is constant.

As for the volume factor $a_{0}^{d-1}V_{d-1}$, we can relate it to the total
number of the Hubble patches $N_{H}$ in the following way. First, the maximal
length light would travel from $x^{0}=0$ to $x^{0}=\infty$ 111Here, the choice
of initial time $x^{0}=t_{0}$ is arbitrary, and $Q$ does not depend on it.  is
determined by $ds^{2}=0\leftrightarrow
dr/dx^{0}=1/a(x^{0})=e^{-H_{I}x^{0}}/a_{0}$, which leads to

$\displaystyle r_{\rm max}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}{dr\over dx^{0}}dx^{0}={1\over a_{0}H_{I}}.$
(23)

This means that the Hubble volume is given by

$\displaystyle V_{H}$ $\displaystyle:=$ $\displaystyle
a_{0}^{d-1}\Omega_{d-2}\int_{0}^{r_{\rm max}}r^{d-2}dr={H_{I}^{-(d-1)}\over
d-1}\Omega_{d-2}~{},$ (24)

where $\Omega_{d-2}$ is a solid angle of $d-2$ dimensional polar coordinates.
Then, the total number of the Hubble patches $N_{H}$ is

$\displaystyle N_{H}=\frac{a_{0}^{d-1}V_{d-1}}{V_{H}}~{}.$ (25)

By using this relation, we can rewrite $Q$ as

$\displaystyle Q$ $\displaystyle=$ $\displaystyle N_{H}\times{A_{H}\over
4G_{N}}\times{(d-2)\over 2}T_{H}\beta_{0}~{},$ (26) $\displaystyle A_{H}:$
$\displaystyle=$ $\displaystyle R_{H}^{d-2}\Omega_{d-2}~{},$ (27)

where $A_{H}$ is an area of the Hubble horizon, which is equivalent to the
cosmological horizon for de Sitter spacetime with a radius $R_{H}=1/H_{I}$,
since the metric with $a(x^{0})=a_{0}e^{H_{I}x^{0}}$ describes (a part of) de
Sitter spacetime, and thus $T_{H}=H_{I}/(2\pi)$ is a temperature of Hubble or
de Sitter horizon. Under this equivalence, the constant $\rho_{I}$ for $w=-1$
is identified with the cosmological constant as
$\Lambda_{dS}=(d-1)(d-2)/(2R_{H}^{2})=8\pi G_{N}\rho_{I}$, which implies
$R_{H}=1/H_{I}$. Thus, taking an initial condition of $\beta_{0}$ as

$\displaystyle\beta_{0}={1\over\tilde{T}_{H}},\quad\tilde{T}_{H}:={(d-2)\over
2}T_{H}~{},$ (28)

and focusing on one Hubble patch i.e. $N_{H}=1$, $Q$ reproduces the
Bekenstein-Hawking entropy formula Bekenstein (1973, 1974); Hawking (1975) for
de Sitter spacetime Gibbons and Hawking (1977) as $Q={A_{H}\over 4G_{N}}$.
This fact justifies our interpretation that the conserved charge $Q$ is
entropy and $\beta(x^{0})$ is proportional to an inverse temperature, so that
the temperature during the inflation era is given by
$T_{H}(x^{0})={2\over(d-2)\beta(x^{0})}$. In other words, we can derive the
Bekenstein-Hawking entropy formula for a (static) de Sitter spacetime Gibbons
and Hawking (1977) from the entropy $Q(x^{0})$ in the inflation era with an
appropriate choice of the initial inverse temperature $\beta_{0}$. Note
however that the entropy $Q$ is uniformly distributed inside the Hubble
horizon but is not concentrated on the horizon. An overall factor $(d-2)/2$
for $T_{H}$ in $\beta_{0}$ becomes unity at $d=4$ and this factor also appears
in the Schwarzschild blackhole Aoki et al. (2021b).

####  III.3.2 Matter dominated era

At $x^{0}=t_{M}$, the inflation era ends, and the matter dominated era starts
typically due to the oscillations of inflaton around the origin.

Solutions with $\omega=0$ lead to time-dependent $H$ and $\rho$ during the
matter dominated era as

$\displaystyle H(x^{0})$ $\displaystyle=$ $\displaystyle
H_{I}\left(\frac{a(t_{M})}{a(x^{0})}\right)^{\frac{d-1}{2}},~{}\rho(x^{0})=\rho_{I}\left(\frac{a(t_{M})}{a(x^{0})}\right)^{d-1},~{}~{}~{}~{}~{}$
(29)

while $\beta(x^{0})=\beta(t_{M}):=T_{M}^{-1}$, $Q$ and

$\displaystyle
E(x^{0})=V_{d-1}\rho_{I}a_{0}^{d-1}e^{(d-1)H_{I}t_{M}}=V_{d-1}\rho_{I}a(t_{M})^{d-1}~{}~{}~{}$
(30)

are all constant during the matter dominated era. Using the law of
equipartition, the energy can be also expressed as

$\displaystyle E(t_{M})$ $\displaystyle=$ $\displaystyle{g_{M}N_{M}\over
2}\varepsilon_{M}~{},$ (31)

where $N_{M}$ is a number of matter particles, $g_{M}$ is an effective degrees
of freedom for a particle, and $\varepsilon_{M}$ is the effective temperature
(energy) of matter particle. Thus a number of particle $N_{M}$ does not change
during the matter dominated era. While $\varepsilon_{M}$ is different from
$T_{M}$ in general, $t_{M}$ and $T_{M}$ are related as

$\displaystyle t_{M}$ $\displaystyle=$
$\displaystyle{1\over(d-1)H_{I}}\log{T_{M}\over\tilde{T}_{H}}~{},$ (32)

where a positivity of $t_{M}>0$ implies $T_{M}>\tilde{T}_{H}$, which is easily
satisfied, since the temperature increases during the inflation. Note also
that while $\beta(x^{0})$ is continuous at $x^{0}=t_{M}$, the temperature
$T_{H}(t_{M})$ at the end of the inflation era and the matter temperature
$T_{M}$ are different as $T_{M}={(d-2)\over 2}T_{H}(t_{M})$, due to the new
relation $T_{M}=1/\beta(t_{M})$, which is required for $Q(t_{M})$ to be the
matter entropy during the matter dominated era.

####  III.3.3 Radiation dominated era

After $x^{0}=t_{R}$, the radiation era starts. Solutions with $\omega=1/(d-1)$
are

$\displaystyle H(x^{0})$ $\displaystyle=$ $\displaystyle
H(t_{R})\left(\frac{a(t_{R})}{a(t)}\right)^{\frac{d}{2}}~{},$ (33)
$\displaystyle\rho(x^{0})$ $\displaystyle=$
$\displaystyle\rho(t_{R})\left(\frac{a(t_{R})}{a(x^{0})}\right)^{d}~{},\
\beta(x^{0})=T_{M}^{-1}\frac{a(x^{0})}{a(t_{R})}~{}.$ (34)

If the conserved charge density $q(x^{0}):=\rho(x^{0})\beta(x^{0})$ is
identified with the entropy density for radiations as

$\displaystyle\rho(x^{0})\beta(x^{0})$ $\displaystyle=$
$\displaystyle{\rho(x^{0})+P(x^{0})\over T_{R}(x^{0})}={d\over
d-1}{\rho(x^{0})\over T_{R}(x^{0})}~{},$ (35)

we can read the relation between radiation temperature $T_{R}(x^{0})$ and
$\beta(x^{0})$ as222 Due to the direct transition from the matter era to the
radiation era imposed by hand in our simplified model, this definition of
$T_{R}(x^{0})$ is different from the usual definition of radiation temperature
$\rho(x^{0})=\sigma_{d}T^{d}(x^{0})$ in cosmology. In particular, when the
decay rate of inflaton $\Gamma_{\phi}$ satisfies $\Gamma_{\phi}\ll H_{I}$, the
ordinary reheating temperature is given by
$T_{R}:=T(t_{R})\sim\sqrt{\Gamma_{\phi}M_{\rm Pl}}$, which is not related to
$T_{M}$. We expect that this mismatch will be fixed in more realistic models.

$\displaystyle\beta(x^{0})={d\over d-1}{1\over T_{R}(x^{0})}\Rightarrow
T_{R}(x^{0})={d\over d-1}T_{M}\frac{a(t_{R})}{a(x^{0})}~{},~{}~{}~{}~{}$ (36)

which shows an extra factor $d/(d-1)$ between $\beta(x^{0})$ and the inverse
temperature. Thus, while $\beta(x^{0})$ is continuous, the matter temperature
$T_{M}$ and the radiation temperature $T_{R}(t_{R})$ at $x^{0}=t_{R}$ are
different as $(d-1)T_{R}(t_{R})=dT_{M}$, as in the case at $x^{0}=t_{M}$.

It is non-trivial and interesting that this time dependent temperature
satisfies the Stefan-Boltzmann law as

$\displaystyle\rho(x^{0})$ $\displaystyle=$ $\displaystyle
C_{R}T_{R}^{d}(x^{0}),\quad C_{R}:=\rho(t_{R})\left({d-1\over
dT_{M}}\right)^{d},~{}~{}$ (37)

which again justifies our interpretation on $Q$ and $\beta(x^{0})$. Note that
we have not assumed thermal equilibrium for radiation to obtain the above
result.

In addition, a starting time of radiation dominated era $t_{R}$ can be
estimated by assuming that $C_{R}$ reproduces the Stefan-Boltzmann constant
for black-body radiations in $d$ dimensions Cardoso and de Castro (2005) as

$\displaystyle
C_{R}=\sigma_{d}:=g_{R}{(d-1)\Gamma(d/2)\zeta(d)\over\pi^{d\over 2}}~{},$ (38)

where $g_{R}$ represents a number of degrees of freedom for radiations
($g_{R}=d-2$ for photons), and $\zeta(x)$ is the Riemann zeta function. In
such a case, we obtain

$\displaystyle t_{R}-t_{M}$ $\displaystyle=$ $\displaystyle{2\over
d-1}\left[\left({d-1\over dT_{M}}\right)^{d\over 2}\sqrt{(d-1)(d-2)\over 16\pi
G_{N}\sigma_{d}}-{1\over H_{I}}\right],~{}~{}~{}~{}~{}$ (39)

and a positivity of $t_{R}-t_{M}$ implies

$\displaystyle T_{M}^{d\over 2}<\sqrt{\pi}M_{\rm
pl}\tilde{T}_{H}\left({d-1\over d}\right)^{d\over 2}\sqrt{(d-1)(d-2)\over
4\sigma_{d}},$ (40)

which gives

$\displaystyle T_{M}^{2}<{27\over 32}\sqrt{10\over\pi}M_{\rm pl}T_{H}\simeq
1.5M_{\rm pl}T_{H}$ (41)

at $d=4$, where $M_{\rm pl}=1/\sqrt{G_{N}}$ is the Planck mass.

##  IV A comparison with a scalar model

Here we numerically solve the inflaton dynamics at $d=4$ and compare it with
the analytical results in the previous section. Equations are

$\displaystyle\ddot{\phi}+(3H+\Gamma)\dot{\phi}+\frac{\partial
V}{\partial\phi}=0~{},$ (42)
$\displaystyle\dot{\rho}_{R}+4H\rho_{R}=\Gamma(\rho_{\phi}+p_{\phi})~{},$ (43)
$\displaystyle\dot{\beta}-3H\frac{P}{\rho}\beta=0~{},\quad H^{2}=\frac{8\pi
G_{N}}{3}\rho~{},$ (44)
$\displaystyle\rho=\rho_{\phi}+\rho_{R}~{},~{}\rho_{\phi}=\frac{1}{2}\dot{\phi}^{2}+V(\phi)~{},$
(45) $\displaystyle
P=p_{\phi}+\frac{1}{3}\rho_{R}~{},~{}p_{\phi}=\frac{1}{2}\dot{\phi}^{2}-V(\phi)~{},$
(46)

where $\Gamma$ represents the decay rate of inflaton. As for the inflaton
potential, we choose

$\displaystyle
V(\phi)=\frac{m_{\phi}^{2}}{2}\phi^{2}+\frac{\lambda}{4!}\phi^{4}$ (47)

as a toy model though its CMB predictions are already excluded by Planck2018
Akrami et al. (2020). We also introduce the radiation entropy

$\displaystyle S=V_{3}a^{3}(x^{0})\rho_{R}(x^{0})\beta(x^{0})~{},$ (48)

which must coincide with $Q$ after the reheating.

Numerical results are obtained with following parameters and initial
conditions:

$\displaystyle\lambda$ $\displaystyle=10^{-2}~{},\quad m_{\phi}=0.1M_{\rm
Pl}~{},\quad\Gamma=10^{-2}m_{\phi}~{},$ (49) $\displaystyle\phi(0)$
$\displaystyle=5M_{\rm
Pl}~{},\quad\dot{\phi}(0)=0~{},\quad\beta(0)=T_{H}^{-1}~{}.$ (50)

Fig. 1 (1st) shows the EoS $\omega(x^{0})$ as a function of $m_{\phi}t$
(logarithmic scale in the horizontal axis), which counts the number of
oscillations of $\phi(x^{0})$ around the origin. The plot indicates that a
slow rolling during inflation ends and an oscillation in the matter dominated
era starts around $m_{\phi}t\simeq O(1)$, while it ends at $m_{\phi}t\simeq
m_{\phi}\Gamma^{-1}=100$.

Figure 1: (1st) EoS $\omega(x^{0})$. (2nd) Energy densities $\rho_{\phi}$
(blue) and $\rho_{R}$(orange). (3rd) The inverse temperature $\beta$ by
analytic (orange) and numerical (blue) calculations. (4th) The conserved
charge $Q$ (blue), the radiation entropy $S$ (orange) and energy $E$ (green)
for this model, together with an analytic prediction of $E$ (red).

Fig. 1 (2nd) shows a log plot of energy densities of matters $\rho_{\phi}$ and
radiations $\rho_{R}$, which also shows that the radiation dominated era
starts at $m_{\phi}t\sim 100$, where $\rho_{R}$ becomes comparable to
$\rho_{\phi}$.

In Fig. 1 (3rd), we compare the inverse temperature $\beta$ between the
numerical calculation (blue) and the analytic result (orange). A qualitative
behavior of the inverse temperature in the numerical calculation is well
captured by the analytic prediction, whose constant behavior tells us that the
matter dominated era starts at around $m_{\phi}t\simeq 3$ and ends at
$m_{\phi}t\simeq 100$. The 4th panel of Fig. 1 shows $Q$ (blue), $S$ (orange),
and total energy $E$ (green), together with the analytic prediction of $E$
(red) for a comparison. One can see that $Q$ is conserved during the
transition dynamics and that it is converted to the radiation entropy $S$ in
the end. In addition, the analytic prediction of $E$ well describes the
numerical result qualitatively.

##  V Summary

We have investigated properties of a conserved charge obtained by a vector
field $\xi^{\mu}=-\beta\delta^{\mu}_{0}$ during dynamical evolutions of a
homogeneous and isotropic expanding Universe with several different constant
equations of state. While the conserve charge in the inflation era with an
appropriate choice of the initial temperature agrees with the Bekenstein-
Hawking entropy for de Sitter spacetime, it is identified as the matter
entropy with a constant temperature (energy) in the matter dominated era.
Finally the charge becomes the radiation entropy and the time dependent
temperature is shown to satisfy the Stefan-Boltzmann law during the radiation
dominated era. As a concrete example of such a dynamical transitions, we have
numerically studied a scalar (inflation) model with radiations and found that
time evolutions of the conserved charge and temperature are qualitatively
reproduced. Our results give strong evidences on the interpretation that the
conserved charge in general relativity Aoki et al. (2021b) is indeed entropy
of the system together with the spacetime dependent temperature.

###### Acknowledgements.

This work is supported in part by the Grant-in-Aid of the Japanese Ministry of
Education, Sciences and Technology, Sports and Culture (MEXT) for Scientific
Research (Nos. JP22H00129) and by “Program for Promoting Researches on the
Supercomputer Fugaku” (Simulation for basic science: from fundamental laws of
particles to creation of nuclei). K.K. would like to thank Yukawa Institute
for Theoretical Physics, Kyoto University for the support and the hospitality
during his stay by the long term visiting program.

## References

  * Bekenstein (1973) J. D. Bekenstein, Phys. Rev. D 7, 2333 (1973).
  * Hawking (1975) S. W. Hawking, Commun. Math. Phys. 43, 199 (1975), [Erratum: Commun.Math.Phys. 46, 206 (1976)].
  * Bekenstein (1974) J. D. Bekenstein, Phys. Rev. D 9, 3292 (1974).
  * Strominger and Vafa (1996) A. Strominger and C. Vafa, Phys. Lett. B 379, 99 (1996), eprint hep-th/9601029.
  * Maldacena et al. (1999) J. M. Maldacena, G. W. Moore, and A. Strominger (1999), eprint hep-th/9903163.
  * Thorne et al. (1986) K. S. Thorne, R. H. Price, and D. A. Macdonald, eds., _BLACK HOLES: THE MEMBRANE PARADIGM_ (1986), ISBN 978-0-300-03770-8.
  * Aoki et al. (2021a) S. Aoki, T. Onogi, and S. Yokoyama, Int. J. Mod. Phys. A 36, 2150098 (2021a), eprint 2005.13233.
  * Aoki et al. (2021b) S. Aoki, T. Onogi, and S. Yokoyama, Int. J. Mod. Phys. A 36, 2150201 (2021b), eprint 2010.07660.
  * Aoki and Onogi (2022) S. Aoki and T. Onogi, Int. J. Mod. Phys. A 37, 2250129 (2022), eprint 2201.09557.
  * Aoki (2022) S. Aoki, PTEP 2022, 123A02 (2022), eprint 2206.00283.
  * Hawking (1976) S. W. Hawking, Phys. Rev. D 14, 2460 (1976).
  * Raju (2022) S. Raju, Phys. Rept. 943, 1 (2022), eprint 2012.05770.
  * Gibbons and Hawking (1977) G. W. Gibbons and S. W. Hawking, Phys. Rev. D 15, 2738 (1977).
  * Cardoso and de Castro (2005) T. R. Cardoso and A. S. de Castro, Rev. Bras. Ens. Fis. 27, 559 (2005), eprint quant-ph/0510002.
  * Akrami et al. (2020) Y. Akrami et al. (Planck), Astron. Astrophys. 641, A10 (2020), eprint 1807.06211.


marginparsep has been altered.
topmargin has been altered.
marginparwidth has been altered.
marginparpush has been altered.
The page layout violates the ICML style. Please do not change the page layout,
or include packages like geometry, savetrees, or fullpage, which change it for
you. We’re not able to reliably undo arbitrary changes to the style. Please
remove the offending package(s), or layout-changing commands and try again.

Bayesian Low-Rank Interpolative Decomposition for Complex Datasets

Anonymous Authors1

††footnotetext:  1Anonymous Institution, Anonymous City, Anonymous Region,
Anonymous Country.

###### Abstract

In this paper, we introduce a probabilistic model for learning interpolative
decomposition (ID), which is commonly used for feature selection, low-rank
approximation, and identifying hidden patterns in data, where the matrix
factors are latent variables associated with each data dimension. Prior
densities with support on the specified subspace are used to address the
constraint for the magnitude of the factored component of the observed matrix.
Bayesian inference procedure based on Gibbs sampling is employed. We evaluate
the model on a variety of real-world datasets including CCLE $EC50$, CCLE
$IC50$, CTRP $EC50$, and MovieLens 100K datasets with different sizes, and
dimensions, and show that the proposed Bayesian ID GBT and GBTN models lead to
smaller reconstructive errors compared to existing randomized approaches.

#### Keywords:

Interpolative decomposition, Low-rank approximation, Bayesian inference,
Hierarchical model.

##  1 Introduction

Matrix factorization methods such as singular value decomposition (SVD),
factor analysis, principal component analysis (PCA), and independent component
analysis (ICA) have been used extensively over the years to reveal hidden
structures of matrices in many fields of science and engineering such as deep
learning (Liu et al., 2015), recommendation systems (Comon et al., 2009; Lu,
2021a; 2022), computer vision (Goel et al., 2020), clustering and
classification (Li et al., 2009; Wang et al., 2013; Lu, 2021b), collaborative
filtering (Marlin, 2003; Lim & Teh, 2007; Mnih & Salakhutdinov, 2007; Raiko et
al., 2007; Chen et al., 2009; Brouwer & Lio, 2017; Lu, 2022) and machine
learning in general (Lee & Seung, 1999).

Moreover, low-rank approximations are essential in modern data science. Low-
rank matrix approximation with respect to the Frobenius norm - minimizing the
sum squared differences to the target matrix - can be easily solved with
singular value decomposition. For many applications, however, it is sometimes
advantageous to work with a basis that consists of a subset of the columns of
the observed matrix itself (Halko et al., 2011; Martinsson et al., 2011). The
interpolative decomposition (ID) provides one such approximation. Its
distinguishing feature is that it reuses columns from the original matrix.
This enables it to preserve matrix properties such as sparsity and
nonnegativity that also help save space in memory.

Interpolative decomposition is widely used as a feature selection tool that
extracts the essence and allows dealing with big data which is originally too
large to fit into the RAM. In addition, it removes the non-relevant parts of
the data which consist of error and redundant information (Liberty et al.,
2007; Halko et al., 2011; Martinsson et al., 2011; Arı et al., 2012; Lu,
2021a). In the meantime, finding the indices associated with the spanning
columns is frequently valuable for the purpose of data interpretation and
analysis, it can be very useful to identify a subset of the columns that
distills the information in the matrix. If the columns of the observed matrix
have some specific interpretations, e.g., they are transactions in a
transaction dataset, then the columns of the factored matrix in ID will have
the same meaning as well. The factored matrix obtained from interpolative
decomposition is also numerically stable since its maximal magnitude is
limited to a certain range.

On the other hand, matrix factorizations can also be thought of as statistical
models in which we seek the factorization that offers the maximum marginal
likelihood (MML) for the underlying data (Arı et al., 2012; Brouwer & Lio,
2017). Probabilistic interpretations are investigated for many popular matrix
decompositions in the literature such as real-valued matrix factorization,
nonnegative matrix factorization (NMF) (Brouwer & Lio, 2017; Lu & Ye, 2022),
principal component analysis (Tipping & Bishop, 1999), and generalized to
tensor factorizations (Schmidt & Mohamed, 2009). In the meantime,
probabilistic models can easily accommodate constraints on the specific range
of the factored matrix.

In this light, we focus on the Bayesian ID (BID) of underlying matrices. The
ID problem of observed matrix $\bm{A}$ can be stated as
$\bm{A}=\bm{C}\bm{W}+\bm{E}$, where
$\bm{A}=[\bm{a}_{1},\bm{a}_{2},\ldots,\bm{a}_{N}]\in{}^{M\times N}$ is
approximately factorized into an $M\times K$ matrix $\bm{C}\in{}^{M\times K}$
containing $K$ basis columns of $\bm{A}$ and a $K\times N$ matrix
$\bm{W}\in{}^{K\times N}$ with entries no larger than 1 in magnitude 111A
weaker construction is to assume no entry of $\bm{Y}$ has an absolute value
greater than 2. See proof of the existence of the decomposition in Lu
(2021a).; the noise is captured by matrix $\bm{E}\in{}^{M\times N}$. Training
such models amounts to finding the best rank-$K$ approximation to the observed
$M\times N$ target matrix $\bm{A}$ under the given loss function. Let
$\bm{r}\in\\{0,1\\}^{N}$ be the state vector with each element indicating the
type of the corresponding column, basis column or interpolated (remaining)
column: if $r_{n}=1$, then the $n$-th column $\bm{a}_{n}$ is a basis column;
if $r_{n}=0$, then $\bm{a}_{n}$ is interpolated using the basis columns plus
some error term. Suppose further $J$ is the set of the indices of the selected
basis columns, $I$ is the set of the indices of the interpolated columns such
that $J\cup I=\\{1,2,\ldots,N\\}$, $J=J(\bm{r})=\\{n|r_{n}=1\\}_{n=1}^{N}$,
and $I=I(\bm{r})=\\{n|r_{n}=0\\}_{n=1}^{N}$. Then $\bm{C}$ can be described as
$\bm{C}=\bm{A}[:,J]$ where the colon operator implies all indices. The
approximation $\bm{A}\approx\bm{C}\bm{W}$ can be equivalently stated that
$\bm{A}\approx\bm{C}\bm{W}=\bm{X}\bm{Y}$ where $\bm{X}\in{}^{M\times N}$ and
$\bm{Y}\in{}^{N\times N}$ such that $\bm{X}[:,J]=\bm{C}$,
$\bm{X}[:,I]=\mathbf{0}\in{}^{M\times(N-K)}$; $\bm{W}=\bm{Y}[J,:]$. We also
notice that there exists an identity matrix $\bm{I}\in{}^{K\times K}$ in
$\bm{W}$ and $\bm{Y}$:

$\bm{I}=\bm{W}[:,J]=\bm{Y}[J,J].$ (1)

To find the low-rank ID of $\bm{A}\approx\bm{C}\bm{W}$ then can be transformed
into the problem of finding the $\bm{A}\approx\bm{X}\bm{Y}$ with state vector
$\bm{r}$ recovering the submatrix $\bm{C}$ (see Figure 1). To evaluate the
approximation, reconstruction error measured by mean squared error (MSE or
Frobenius norm) is minimized (assume $K$ is known):

$\mathop{\min}_{\bm{W},\bm{Z}}\,\,\frac{1}{MN}\sum_{n=1}^{N}\sum_{m=1}^{M}\left(a_{mn}-\bm{x}_{m}^{\top}\bm{y}_{n}\right)^{2},$
(2)

where $\bm{x}_{m}$, $\bm{y}_{n}$ are the $m$-th row and $n$-th column of
$\bm{X}$, $\bm{Y}$ respectively. In this paper, we approach the magnitude
constraint in $\bm{W}$ and $\bm{Y}$ by considering the Bayesian ID model as a
latent factor model and we describe a fully specified graphical model for the
problem and employ Bayesian learning methods to infer the latent factors. In
this sense, explicit magnitude constraints are not required on the latent
factors, since this is naturally taken care of by the appropriate choice of
prior distribution; here we use general-truncated-normal prior.

The main contribution of this paper is to propose a novel Bayesian ID method
which is flexible called the GBT model. To further favor flexibility and
insensitivity on the hyperparameter choices, we also propose the hierarchical
model known as the GBTN algorithm, which has simple conditional density forms
with little extra computation, Meanwhile, the method is easy to implement. We
show that our method can be successfully applied to the large, sparse, and
very imbalanced Movie-User dataset, containing 100,000 user/movie ratings.

##  2 Related Work

The most popular algorithm to compute the low-rank ID approximation is the
randomized ID (RID) algorithm (Liberty et al., 2007). At a high level, the
algorithm randomly samples $M>K$ columns from $\bm{A}$, uses column-pivoted QR
(Lu, 2021a) to select $K$ of those $M$ columns for basis matrix $\bm{C}$, and
then computes $\bm{W}$ via least squares. $M$ is usually set to $M=1.2K$ to
oversample the columns so as to capture a large portion of the range of matrix
$\bm{A}$.

We propose a new Bayesian approach for interpolative decomposition, which is
designed to have a strict constraint on the magnitude of the factored matrix
$\bm{W}$. The proposed GBT and GBTN models introduce a general-truncated-
normal density over the factored matrix $\bm{W}$ that finds interpolated
columns to reconstruct the observed matrix.

Figure 1: Demonstration of the interpolative decomposition of a matrix where
the yellow vector denotes the basis columns of matrix $\bm{A}$, white entries
denote zero, purple entries denote one, blue and black entries denote elements
that are not necessarily zero. The Bayesian ID models find the approximation
$\bm{A}\approx\bm{X}\bm{Y}$ and the post processing procedure finds the
approximation $\bm{A}\approx\bm{C}\bm{W}$.

###  2.1 Probability Distributions

We introduce all notations and probability distribution in this section.

$\mathcal{N}(x|\mu,\tau^{-1})=\sqrt{\frac{\tau}{2\pi}}\exp\\{-\frac{\tau}{2}(x-\mu)^{2}\\}$
is a Gaussian distribution with mean $\mu$ and precision $\tau$ (variance
$\sigma^{2}=\tau^{-1}$).

$\mathcal{G}(x|\alpha,\beta)=\frac{\beta^{\alpha}}{\Gamma(\alpha)}x^{\alpha-1}\exp\\{-\beta
x\\}u(x)$ is a Gamma distribution where $\Gamma(\cdot)$ is the gamma function
and $u(x)$ is the unit step function that has a value of $1$ when $x\geq 0$
and 0 otherwise.

$\mathcal{G}^{-1}(x|\alpha,\beta)=\frac{\beta^{\alpha}}{\Gamma(\alpha)}x^{-\alpha-1}\exp\\{-\frac{\beta}{x}\\}u(x)$
is an inverse-Gamma distribution.

$\mathcal{TN}(x|\mu,\tau^{-1})=\frac{\sqrt{\frac{\tau}{2\pi}}\exp\\{-\frac{\tau}{2}(x-\mu)^{2}\\}}{1-\Phi(-\mu\sqrt{\tau})}u(x)$
is a truncated-normal (TN) with zero density below $x=0$ and renormalized to
integrate to one. $\mu$ and $\tau$ are known as the “parent mean” and “parent
precision”. $\Phi(\cdot)$ is the cumulative distribution function of standard
normal density $\mathcal{N}(0,1)$.

$\mathcal{GTN}(x|\mu,\frac{1}{\tau},a,b)=\frac{\sqrt{\frac{\tau}{2\pi}}\exp\\{-\frac{\tau}{2}(x-\mu)^{2}\\}}{\Phi((b-\mu)\cdot\sqrt{\tau})-\Phi((a-\mu)\cdot\sqrt{\tau})}$$u(x|a,b)$
is a general-truncated-normal (GTN) with zero density below $x=a$ or above
$x=b$ and renormalized to integrate to one where $u(x|a,b)$ is a step function
that has a value of 1 when $a\leq x\leq b$ and 0 otherwise. Similarly, $\mu$
and $\tau$ are known as the “parent mean” and “parent precision” of the normal
distribution. When $a=0$ and $b=\infty$, the GTN distribution reduces to a TN
density.

##  3 Bayesian Interpolative Decomposition

###  3.1 Bayesian GBT and GBTN Models for ID

(a) GBT.

(b) GBTN.

Figure 2: Graphical representation of GBT and GBTN models. Orange circles
represent observed and latent variables, green circles denote prior variables,
and plates represent repeated variables. “/” in the variable represents “or”,
and comma “,” in the variable represents “and”. Parameters $a,b$ are fixed
with $a=-1,b=1$ in our case; while a weaker construction can set them to
$a=-2,b=2$.

We view the data $\bm{A}$ as being produced according to the probabilistic
generative process shown in Figure 2. The observed $(m,n)$-th entry $a_{mn}$
of matrix $\bm{A}$ is modeled using a Gaussian likelihood function with
variance $\sigma^{2}$ and mean given by the latent decomposition
$\bm{x}_{m}^{\top}\bm{y}_{n}$ (Eq. (2)),

$p(a_{mn}|\bm{x}_{m}^{\top}\bm{y}_{n},\sigma^{2})=\mathcal{N}(a_{mn}|\bm{x}_{m}^{\top}\bm{y}_{n},\sigma^{2}).$
(3)

We choose a conjugate prior over the data variance, an inverse-Gamma
distribution with shape $\alpha_{\sigma}$ and scale $\beta_{\sigma}$,

$p(\sigma^{2}|\alpha_{\sigma},\beta_{\sigma})=\mathcal{G}^{-1}(\sigma^{2}|\alpha_{\sigma},\beta_{\sigma}).$
(4)

While it can also be equivalently given a conjugate Gamma prior over the
precision and we shall not repeat the details.

We treat the latent variables $y_{kl}$’s as random variables. And we need
prior densities over these latent variables to express beliefs for their
values, e.g., constraint with magnitude smaller than 1 in this context though
there are many other constraints (nonnegativity in Lu & Ye (2022), semi-
nonnegativity in Ding et al. (2008), discreteness in Gopalan et al. (2014;
2015)). Here we assume further that the latent variable $y_{kl}$’s are
independently drawn from a general-truncated-normal prior

$p(y_{kl}|\cdot)=\mathcal{GTN}(y_{kl}|\mu_{kl},(\tau_{kl})^{-1},a=-1,b=1).$
(5)

This prior serves to enforce the constraint on the components $\bm{Y}$ with no
entry of $\bm{Y}$ having an absolute value greater than 1, and is conjugate to
the Gaussian likelihood. The posterior density is also a general-truncated-
normal distribution. While in a weaker construction of interpolative
decomposition, the constraint on the magnitude can be loosened to 2; the prior
is flexible in that the parameters can be then set to $a=-2,b=2$ accordingly.

#### Hierarchical prior

To further favor flexibility, we choose a convenient joint hyperprior density
over the parameters $\\{\mu_{kl},\tau_{kl}\\}$ of GTN prior in Eq. (5),
namely, the GTN-scaled-normal-Gamma (GTNSNG) prior,

$\displaystyle\,\,\,\,\,\,\,\,p(\mu_{kl},\tau_{kl}|\cdot)=\mathcal{GTNSNG}(\mu_{kl},\tau_{kl}|\mu_{\mu},\frac{1}{\tau_{\mu}},\alpha_{t},\beta_{t})$
(6)
$\displaystyle=\big{\\{}\Phi((b-\mu_{\mu})\cdot\sqrt{\tau_{\mu}})-\Phi((a-\mu_{\mu})\cdot\sqrt{\tau_{\mu}})\big{\\}}$
$\displaystyle\,\,\,\,\,\,\,\,\cdot\mathcal{N}(\mu_{kl}|\mu_{\mu},(\tau_{\mu})^{-1})\cdot\mathcal{G}(\tau_{kl}|a_{t},b_{t}).$

This prior can decouple parameters $\mu_{kl},\tau_{kl}$, and the posterior
conditional densities of them are normal, and Gamma respectively due to this
convenient scale.

#### Terminology

There are three types of choices we make that determine the specific type of
matrix decomposition model we use, namely, the likelihood function, the priors
we place over the factor matrices $\bm{W}$ and $\bm{Z}$, and whether we use
any further hierarchical priors. We will call the model by the density
function in the order of the types of likelihood and priors. For example, if
the likelihood function for the model is chosen to be a Gaussian density, and
the two prior density functions are selected to be exponential density and
Gaussian density functions respectively, then the model will be denoted as
Gaussian Exponential-Gaussian (GEG) model. Sometimes, we will put a hyperprior
over the parameters of the prior density functions, e.g., we put a Gamma prior
over the Gaussian density, then it will further be termed as a Gaussian
Exponential-Gaussian Gamma (GEGA) model. In this sense, the proposed
hierarchical model is named the GBT and GBTN model where $B$ stands for Beta-
Bernoulli density intrinsically.

###  3.2 Gibbs Sampler

Algorithm 1 Gibbs sampler for GBT and GBTN ID models. The procedure presented
here may not be efficient but is explanatory. A more efficient one can be
implemented in a vectorized manner. By default, uninformative priors are
$a=-1,b=1,\alpha_{\sigma}=0.1,\beta_{\sigma}=1$,
($\\{\mu_{kl}\\}=0,\\{\tau_{kl}\\}=1$) for GBT, ($\mu_{\mu}=0$,
$\tau_{\mu}=0.1,\alpha_{t}=\beta_{t}=1$) for GBTN.

1: for $t=1$ to $T$ do

2: Sample state vector $\bm{r}$ from Eq. (9);

3: Update matrix $\bm{X}$ by $\bm{A}[:,J]$ where index vector $J$ is the index
of $\bm{r}$ with value 1 and set $\bm{X}[:,I]=\mathbf{0}$ where index vector
$I$ is the index of $\bm{r}$ with value 0;

4: Sample $\sigma^{2}$ from $p(\sigma^{2}|\bm{X},\bm{Y},\bm{A})$ in Eq. (10);

5: for $k=1$ to $N$ do

6: for $l=1$ to $N$ do

7: Sample $y_{kl}$ from Eq. (7);

8: (GBTN only) Sample $\mu_{kl}$ from Eq. (11);

9: (GBTN only) Sample $\tau_{kl}$ from Eq. (12);

10: end for

11: end for

12: Report loss in Eq. (2), stop if it converges.

13: end for

14: Report mean loss in Eq. (2) after burn-in iterations.

In this paper, we use Gibbs sampling since it tends to be very accurate at
finding the true posterior. Other than this method, variational Bayesian
inference can be an alternative way but we shall not go into the details. We
shortly describe the posterior conditional density in this section. A detailed
derivation can be found in Appendix A. The conditional density of $y_{kl}$ is
a general-truncated-normal density. Denote all elements of $\bm{Y}$ except
$y_{kl}$ as $\bm{Y}_{-kl}$, we then have

$\displaystyle\,\,\,\,\,\,\,\,p(y_{kl}|\bm{A},\bm{X},\bm{Y}_{-kl},\mu_{kl},\tau_{kl},\sigma^{2})$
(7) $\displaystyle\propto p(\bm{A}|\bm{X},\bm{Y},\sigma^{2})\cdot
p(y_{kl}|\mu_{kl},\tau_{kl})$
$\displaystyle=\prod_{i,j=1}^{M,N}\mathcal{N}\left(a_{ij}|\bm{x}_{i}^{\top}\bm{y}_{j},\sigma^{2}\right)\cdot\mathcal{GTN}(y_{kl}|\mu_{kl},\frac{1}{\tau_{kl}},-1,1)$
$\displaystyle\propto\mathcal{GTN}(y_{kl}|\widetilde{\mu},(\widetilde{\tau})^{-1},a=-1,b=1),$

where $\widetilde{\tau}=\frac{\sum_{i}^{M}x_{ik}^{2}}{\sigma^{2}}+\tau_{kl}$
is the posterior “parent precision” of the GTN distribution, and

$\widetilde{\mu}=\bigg{(}\frac{1}{\sigma^{2}}\sum_{i}^{M}x_{ik}\big{(}a_{il}-\sum_{j\neq
k}^{N}x_{ij}y_{jl}\big{)}+{\color[rgb]{0,0,0}\tau_{kl}\mu_{kl}}\bigg{)}\big{/}\widetilde{\tau}$

is the posterior “parent mean” of the GTN distribution.

Given the state vector $\bm{r}=[r_{1},r_{2},\ldots,r_{N}]^{\top}\in{}^{N}$,
the relation between $\bm{r}$ and the index set $J$ is simple;
$J=J(\bm{r})=\\{n|r_{n}=1\\}_{n=1}^{N}$ and
$I=I(\bm{r})=\\{n|r_{n}=0\\}_{n=1}^{N}$. A new value of state vector $\bm{r}$
is to select one index $j$ from index set $J$ and another index $i$ from index
sets $I$ (we note that $r_{j}=1$ and $r_{i}=0$ for the old state vector
$\bm{r}$) such that

$\displaystyle o_{j}$
$\displaystyle=\frac{p(r_{j}=0,r_{i}=1|\bm{A},\sigma^{2},\bm{Y},\bm{r}_{-ji})}{p(r_{j}=1,r_{i}=0|\bm{A},\sigma^{2},\bm{Y},\bm{r}_{-ji})}$
(8) $\displaystyle=\frac{p(r_{j}=0,r_{i}=1)}{p(r_{j}=1,r_{i}=0)}$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\times\frac{p(\bm{A}|\sigma^{2},\bm{Y},\bm{r}_{-ji},r_{j}=0,r_{i}=1)}{p(\bm{A}|\sigma^{2},\bm{Y},\bm{r}_{-ji},r_{j}=1,r_{i}=0)},$

where $\bm{r}_{-ji}$ denotes all elements of $\bm{r}$ except $j$-th and $i$-th
entries. Trivially, we can set $p(r_{j}=0,r_{i}=1)=p(r_{j}=1,r_{i}=0)$. Then
the full conditionally probability of
$p(r_{j}=0,r_{i}=1|\bm{A},\sigma^{2},\bm{Y},\bm{r}_{-ji})$ can be calculated
by

$p(r_{j}=0,r_{i}=1|\bm{A},\sigma^{2},\bm{Y},\bm{r}_{-ji})=\frac{o_{j}}{1+o_{j}}.$
(9)

Finally, the conditional density of $\sigma^{2}$ is an inverse-Gamma
distribution by conjugacy,

$\displaystyle\,\,\,\,\,\,\,\,p(\sigma^{2}|\bm{X},\bm{Y},\bm{A})=\mathcal{G}^{-1}(\sigma^{2}|\widetilde{\alpha_{\sigma}},\widetilde{\beta_{\sigma}}),$
(10)

where $\widetilde{\alpha_{\sigma}}=\frac{MN}{2}+\alpha_{\sigma}$,
$\widetilde{\beta_{\sigma}}=\frac{1}{2}\sum_{i,j=1}^{M,N}(a_{ij}-\bm{x}_{i}^{\top}\bm{y}_{j})^{2}+\beta_{\sigma}$.

#### Extra updates for the GBTN model

Following the graphical representation of the GBTN model in Figure 2, the
conditional density of $\mu_{kl}$ can be derived,

$\displaystyle\,\,\,\,\,\,\,\,p(\mu_{kl}|\tau_{kl},\mu_{\mu},\tau_{\mu},a_{t},b_{t},y_{kl})$
(11)
$\displaystyle\propto\mathcal{GTN}(y_{kl}|\mu_{kl},(\tau_{kl})^{-1},a=-1,b=1)$
$\displaystyle\,\,\,\,\,\,\,\,\times\mathcal{GTNSNG}(\mu_{kl},\tau_{kl}|\mu_{\mu},(\tau_{\mu})^{-1},\alpha_{t},\beta_{t})$
$\displaystyle\propto\mathcal{N}(\mu_{kl}|\widetilde{m},(\,\widetilde{t}\,)^{-1}),$

where $\widetilde{t}=\tau_{kl}+\tau_{\mu}$, and
$\widetilde{m}=(\tau_{kl}y_{kl}+\tau_{\mu}\mu_{\mu})/\widetilde{t}$ are the
posterior precision and mean of the normal density. Similarly, the conditional
density of $\tau_{kl}$ is,

$\displaystyle\,\,\,\,\,\,\,\,p(\tau_{kl}|\mu_{kl},\mu_{\mu},\tau_{\mu},a_{t},b_{t},y_{kl})$
(12)
$\displaystyle\propto\mathcal{GTN}(y_{kl}|\mu_{kl},(\tau_{kl})^{-1},a=-1,b=1)$
$\displaystyle\,\,\,\,\,\,\,\,\times\mathcal{GTNSNG}(\mu_{kl},\tau_{kl}|\mu_{\mu},(\tau_{\mu})^{-1},\alpha_{t},\beta_{t})$
$\displaystyle\propto\mathcal{G}(\tau_{kl}|\widetilde{a},\widetilde{b}),$

where $\widetilde{a}=a_{t}+1/2$ and
$\widetilde{b}=b_{t}+\frac{(y_{kl}-\mu_{kl})^{2}}{2}$ are the posterior
parameters of the Gamma density. The full procedure is then formulated in
Algorithm 1.

###  3.3 Aggressive Update

In Algorithm 1, we notice that we set $\bm{X}[:,I]=\mathbf{0}$ when the new
state vector $\bm{r}$ is sampled. However, in the next iteration step, the
index set $I$ may be updated in which case one entry $i$ of $I$ may be altered
to have value of 1:

$r_{i}=0\rightarrow r_{i}=1.$

This may cause problems in the update of $y_{kl}$ in Eq. (7), a zero $i$-th
column in $\bm{X}$ cannot update $\bm{Y}$ accordingly. One solution is to
record a proposal state vector $\bm{r}_{2}$ and a proposal factor matrix
$\bm{X}_{2}$ from the $\bm{r}_{2}$ vector. When the update in the next
iteration selects the old state vector $\bm{r}$, the factor matrix $\bm{X}$ is
adopted to finish the updates; while the algorithm chooses the proposal state
vector $\bm{r}_{2}$, the proposal factor matrix $\bm{X}_{2}$ is applied to do
the updates. We call this procedure the aggressive update. The aggressive
sampler for the GBT model is formulated in Algorithm 2. For the sake of
simplicity, we don’t include the sampler for the GBTN model because it may be
done in a similar way.

Algorithm 2 Aggressive Gibbs sampler for GBT ID model. The procedure presented
here may not be efficient but is explanatory. A more efficient one can be
implemented in a vectorized manner. By default, uninformative priors are
$a=-1,b=1,\alpha_{\sigma}=0.1,\beta_{\sigma}=1$,
($\\{\mu_{kl}\\}=0,\\{\tau_{kl}\\}=1$) for GBT.

1: for $t=1$ to $T$ do

2: Sample state vector $\bm{r}$ from $\\{\bm{r}_{1},\bm{r}_{2}\\}$ by Eq. (9);

3: Decide $\bm{Y}$: $\bm{Y}=\bm{Y}_{1}$ if $\bm{r}$ is $\bm{r}_{1}$;
$\bm{Y}=\bm{Y}_{2}$ if $\bm{r}$ is $\bm{r}_{2}$;

4: Update state vector $\bm{r}_{1}=\bm{r}$;

5: Sample proposal state vector $\bm{r}_{2}$ based on $\bm{r}$;

6: Update matrix $\bm{X}$ by $\bm{r}=\bm{r}_{1}$;

7: Update proposal $\bm{X}_{2}$ by $\bm{r}_{2}$;

8: Sample $\sigma^{2}$ from $p(\sigma^{2}|\bm{X},\bm{Y},\bm{A})$ in Eq. (10);

9: Sample $\bm{Y}_{1}=\\{y_{kl}\\}$ using $\bm{X}$;

10: Sample $\bm{Y}_{2}=\\{y_{kl}\\}$ using $\bm{X}_{2}$;

11: Report loss in Eq. (2), stop if it converges.

12: end for

13: Report mean loss in Eq. (2) after burn-in iterations.

###  3.4 Post Processing

The Gibbs sampling algorithm finds the approximation
$\bm{A}\approx\bm{X}\bm{Y}$ where $\bm{X}\in{}^{M\times N}$ and
$\bm{Y}\in{}^{N\times N}$. As stated above, the redundant columns in $\bm{X}$
and redundant rows in $\bm{Y}$ can be removed by the index vector $J$:

$\displaystyle\bm{C}$ $\displaystyle=\bm{X}[:,J]=\bm{A}[:,J],$
$\displaystyle\bm{W}$ $\displaystyle=\bm{Y}[J,:].$

Since the submatrix $\bm{Y}[J,J]=\bm{W}[:,J]$ (Eq. (1)) from the Gibbs
sampling procedure is not enforced to be an identity matrix (as required in
the interpolative decomposition). We need to set it to be an identity matrix
manually. This will basically reduce the reconstructive error further. The
post processing procedure is shown in Figure 1.

Dataset | Rows | Columns | Fraction obs.
---|---|---|---
CCLE $EC50$ | 502 | 48 | 0.632
CCLE $IC50$ | 504 | 48 | 0.965
CTRP $EC50$ | 887 | 1090 | 0.801
MovieLens 100K | 943 | 2946 | 0.072
Table 1: Overview of the CCLE $EC50$, CCLE $IC50$, CTRP $EC50$, and MovieLens
100K datasets, giving the number of rows, columns (after copying twice for
redundancy), and the fraction of entries that are observed. The MovieLens 100K
dataset has a larger fraction of unobserved entries.

Figure 3: Data distribution of CCLE $EC50$, CCLE $IC50$, CTRP $EC50$, and
MovieLens 100K datasets. The former three datasets have complex distributions
and the last one has a Gaussian-like distribution.

##  4 Experiments

To evaluate the strategy and demonstrate the main advantages of the proposed
Bayesian ID method, we conduct experiments with different analysis tasks; and
different data sets including Cancer Cell Line Encyclopedia (CCLE $EC50$ and
CCLE $IC50$ datasets, (Barretina et al., 2012)) and Cancer Therapeutics
Response Portal (CRTP $EC50$, (Seashore-Ludlow et al., 2015)) from
bioinformatics, and MovieLens 100K from movie ratings for different users
(Harper & Konstan, 2015). Following Brouwer & Lio (2017), we preprocess these
datasets by capping high values to 100 and undoing the natural log transform
for the former three datasets. Then we standardize to have zero mean and unit
variance for all datasets and filling missing entries by 0. Moreover, the user
vectors or movie vectors with less than 3 observed entries are cleaned in the
MovieLens 100K dataset. Finally, we copy every column twice in order to
increase redundancy in the matrix. A summary of the four datasets can be seen
in Table 1 and their distributions are shown in Figure 3. The MovieLens 100K
dataset is unbalanced in that it has a larger fraction of unobserved entries
with only $7.2\%$ being observed so that the algorithm may tend to simply
predict all entries to be zero.

In all scenarios, the same parameter initialization is adopted when conducting
different tasks. Experimental evidence reveals that post processing procedure
can increase performance to a minor level, and that the outcomes of the GBT
and GBTN models (aggressive and non-aggressive versions) are relatively
similar. For clarification, we only provide the findings of the GBT model
after an aggressive update and post processing. We compare the results with
randomized ID (RID) algorithm (Liberty et al., 2007). In a wide range of
scenarios across various datasets, the GBT model improves reconstructive error
and leads to performances that are as good or better than the randomized ID
method in low-rank ID approximation. We also notice that all the four datasets
contain missing entries that are filled by zero values. When the fraction of
unobserved data is large, the model can simply predict the entries with zero.
So we also report the mean squared error for observed entries, termed as GBT
(Observed) for the proposed GBT model, and RID (Observed) for the randomized
ID model.

In order to measure overall decomposition performance, we use mean square
error (MSE, Eq. (2)), which measures the similarity between the true and
reconstructive; the smaller the better performance.

###  4.1 Hyperparameters

In this experiments, we use $a=-1,b=1,\alpha_{\sigma}=0.1,\beta_{\sigma}=1$,
($\\{\mu_{kl}\\}=0,\\{\tau_{kl}\\}=1$) for GBT, ($\mu_{\mu}=0$,
$\tau_{\mu}=0.1,\alpha_{t}=\beta_{t}=1$) for GBTN. These are very weak prior
choices and the models are not sensitive to them. As long as the
hyperparameters are set, the observed or unobserved variables are initialized
from random draws as this initialization procedure provides a better initial
guess of the right patterns in the matrices. In all experiments, we run the
Gibbs sampler 500 iterations with a burn-in of 100 iterations and a thinning
of 5 iterations as the convergence analysis shows the algorithm can converge
in less than 50 iterations.

###  4.2 Convergence Analysis

Firstly we show the convergence in terms of iterations on the CCLE $EC50$,
CCLE $IC50$, CTRP $EC50$ and MovieLens 100K datasets. We run each model with
$K=5,10,15,20$ for the CCLE $EC50$ and CCLE $IC50$ datasets, and
$K=20,100,150,200$ for the CTRP $EC50$ and MovieLens 100K datasets; and the
loss is measured by mean squared error (MSE). Figure 4(a) shows the average
convergence results. Figure 4(b) shows autocorrelation coefficients of samples
computed using Gibbs sampling. The coefficients are less than 0.1 when the
lags are more than 10 showing the mixing of the Gibbs sampler is good. In all
experiments, the algorithm converges in less than 50 iterations. On the CCLE
$EC50$ and MovieLens 100K datasets, the sampling are less noisy; while on the
CCLE $IC50$ and CTRP $EC50$ datasets, the sampling seems to be noisy.

Data | $K$ |  GBT |  RID
---|---|---|---
CCLE $EC50$ | 5 | 0.35 | 0.50
10 | 0.22 | 0.36
15 | 0.13 | 0.22
| 20 | 0.07 | 0.08
CCLE $IC50$ | 5 | 0.30 | 0.45
10 | 0.23 | 0.30
15 | 0.16 | 0.17
| 20 | 0.13 | 0.09
CTRP $EC50$ | 20 | 0.58 | 1.11
100 | 0.48 | 1.91
150 | 0.47 | 1.79
MovieLens 100K | 20 | 0.06 | 0.12
100 | 0.05 | 0.24
150 | 0.05 | 0.26
Table 2: Mean squared error measure with various latent dimension $K$ parameters for CCLE $EC50$, CCLE $IC50$, CTRP $EC50$, and MovieLens 100K datasets. Dataset | $K=5$ | $K=10$ | $K=15$
---|---|---|---
CCLE $EC50$ | 0.0 | 4.44e-16 | 4.44e-16
CCLE $IC50$ | 1.99e-2 | 4.44e-16 | 6.66e-16
Dataset | $K=20$ | $K=100$ | $K=150$
CTRP $EC50$ | 4.44e-16 | 2.69e-2 | 1.99e-15
MovieLens 100K | 2.22e-16 | 6.66e-16 | 8.88e-16
Table 3: Magnitude that exceeds 1 in randomized ID algorithm with various
latent dimensions $K$ for different datasets.

(a) Convergence of the models on the CCLE $EC50$, CCLE $IC50$, CTRP $EC50$,
and MovieLens 100K datasets, measuring the training data fit (mean squared
error). The algorithm almost converges in less than 50 iterations.

(b) Averated autocorrelation coefficients of samples of $y_{kl}$ computed
using Gibbs sampling on the CCLE $EC50$, CCLE $IC50$, CTRP $EC50$, and
MovieLens 100K datasets.

(c) Reconstructive analysis on the CCLE $EC50$, CCLE $IC50$, CTRP $EC50$, and
MovieLens 100K datasets with increasing latent dimension $K$.

Figure 4: Convergence results (upper), sampling mixing analysis (middle), and
reconstructive results (lower) on the CCLE $EC50$, CCLE $IC50$, CTRP $EC50$,
and MovieLens 100K datasets for various latent dimensions.

###  4.3 Comparison with Randomized ID

Comparative results for the proposed GBT model and the RID model on the four
datasets are given in Figure 4(c) and Table 2. We also observe that, in all
experiments, the MSEs of the observed entries are worse than those of all the
entries. In the results of CCLE $EC50$ and CCLE $IC50$ datasets, the proposed
GBT model recovers smaller MSE when the latent dimension is lower than 15 (for
both the observed ones and all entries); However, when $K$ tends to be closer
to the rank (24 in this case), the RID becomes better. Similar results can
also be observed on the CTRP $EC50$ dataset whose rank is 545, and on the
MovieLens 100K dataset whose rank is 943. This demonstrates that our proposed
models are more competitive when we need a low-rank ID approximation that can
distill data to a great extent.

For the proposed Bayesian GBT and GBTN models, in all circumstances, the
maximal magnitude of the factored matrix $\bm{W}$ is no greater than 1 because
the prior densities in Bayesian ID models guarantee the magnitude constraint.
The randomized ID, on the other hand, finds a weaker interpolative
decomposition with a maximal magnitude of no more than 2. Table 3 shows the
magnitudes of the randomized ID model for different datasets and different
latent dimensions. We may notice that when $K=100$ on the CTRP $EC50$ dataset,
the maximal magnitude is about $1.0269$. Though this is not a big issue in
many applications, it has potential problems for other applications that have
a high requirement for numerical stability. Advani & O’Hagan (2021) also
reports that the randomized ID may have a magnitude larger than 167 making it
less numerical stable.

##  5 Conclusion

The aim of this paper is to solve the numerical stability issue of the
randomized algorithms in computing the low-rank ID approximation. We propose a
simple and computationally efficient algorithm that requires little extra
computation and is easy to implement for interpolative decomposition. Overall,
we show that the proposed GBT and GBTN models are versatile algorithms that
have good convergence results and better reconstructive performances on both
sparse and dense datasets, especially for low-rank approximation ID. GBT and
GBTN are able to force the magnitude of the factored matrix to be no greater
than 1 such that numerical stability is guaranteed.

## References

  * Advani & O’Hagan (2021) Advani, Rishi and O’Hagan, Sean.  Efficient algorithms for constructing an interpolative decomposition.  _arXiv preprint arXiv:2105.07076_ , 2021.
  * Arı et al. (2012) Arı, Ismail, Cemgil, A Taylan, and Akarun, Lale.  Probabilistic interpolative decomposition.  In _2012 IEEE International Workshop on Machine Learning for Signal Processing_ , pp. 1–6. IEEE, 2012.
  * Barretina et al. (2012) Barretina, Jordi, Caponigro, Giordano, Stransky, Nicolas, Venkatesan, Kavitha, Margolin, Adam A, Kim, Sungjoon, Wilson, Christopher J, Lehár, Joseph, Kryukov, Gregory V, Sonkin, Dmitriy, et al.  The cancer cell line encyclopedia enables predictive modelling of anticancer drug sensitivity.  _Nature_ , 483(7391):603–607, 2012.
  * Brouwer & Lio (2017) Brouwer, Thomas and Lio, Pietro.  Prior and likelihood choices for Bayesian matrix factorisation on small datasets.  _arXiv preprint arXiv:1712.00288_ , 2017.
  * Chen et al. (2009) Chen, Gang, Wang, Fei, and Zhang, Changshui.  Collaborative filtering using orthogonal nonnegative matrix tri-factorization.  _Information Processing & Management_, 45(3):368–379, 2009.
  * Comon et al. (2009) Comon, Pierre, Luciani, Xavier, and De Almeida, André LF.  Tensor decompositions, alternating least squares and other tales.  _Journal of Chemometrics: A Journal of the Chemometrics Society_ , 23(7-8):393–405, 2009.
  * Ding et al. (2008) Ding, Chris HQ, Li, Tao, and Jordan, Michael I.  Convex and semi-nonnegative matrix factorizations.  _IEEE transactions on pattern analysis and machine intelligence_ , 32(1):45–55, 2008.
  * Goel et al. (2020) Goel, Abhinav, Tung, Caleb, Lu, Yung-Hsiang, and Thiruvathukal, George K.  A survey of methods for low-power deep learning and computer vision.  In _2020 IEEE 6th World Forum on Internet of Things (WF-IoT)_ , pp. 1–6. IEEE, 2020.
  * Gopalan et al. (2014) Gopalan, Prem, Ruiz, Francisco J, Ranganath, Rajesh, and Blei, David.  Bayesian nonparametric poisson factorization for recommendation systems.  In _Artificial Intelligence and Statistics_ , pp. 275–283. PMLR, 2014.
  * Gopalan et al. (2015) Gopalan, Prem, Hofman, Jake M, and Blei, David M.  Scalable recommendation with hierarchical poisson factorization.  In _UAI_ , pp. 326–335, 2015.
  * Halko et al. (2011) Halko, Nathan, Martinsson, Per-Gunnar, and Tropp, Joel A.  Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions.  _SIAM review_ , 53(2):217–288, 2011.
  * Harper & Konstan (2015) Harper, F Maxwell and Konstan, Joseph A.  The movielens datasets: History and context.  _Acm transactions on interactive intelligent systems (tiis)_ , 5(4):1–19, 2015.
  * Lee & Seung (1999) Lee, Daniel D and Seung, H Sebastian.  Learning the parts of objects by non-negative matrix factorization.  _Nature_ , 401(6755):788–791, 1999.
  * Li et al. (2009) Li, Tao, Zhang, Yi, and Sindhwani, Vikas.  A non-negative matrix tri-factorization approach to sentiment classification with lexical prior knowledge.  In _Proceedings of the Joint Conference of the 47th Annual Meeting of the ACL and the 4th International Joint Conference on Natural Language Processing of the AFNLP_ , pp. 244–252, 2009.
  * Liberty et al. (2007) Liberty, Edo, Woolfe, Franco, Martinsson, Per-Gunnar, Rokhlin, Vladimir, and Tygert, Mark.  Randomized algorithms for the low-rank approximation of matrices.  _Proceedings of the National Academy of Sciences_ , 104(51):20167–20172, 2007.
  * Lim & Teh (2007) Lim, Yew Jin and Teh, Yee Whye.  Variational Bayesian approach to movie rating prediction.  In _Proceedings of KDD cup and workshop_ , volume 7, pp. 15–21. Citeseer, 2007.
  * Liu et al. (2015) Liu, Baoyuan, Wang, Min, Foroosh, Hassan, Tappen, Marshall, and Pensky, Marianna.  Sparse convolutional neural networks.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 806–814, 2015.
  * Lu (2021a) Lu, Jun.  Numerical matrix decomposition and its modern applications: A rigorous first course.  _arXiv preprint arXiv:2107.02579_ , 2021a.
  * Lu (2021b) Lu, Jun.  A survey on Bayesian inference for Gaussian mixture model.  _arXiv preprint arXiv:2108.11753_ , 2021b.
  * Lu (2022) Lu, Jun.  Matrix decomposition and applications.  _arXiv preprint arXiv:2201.00145_ , 2022.
  * Lu & Ye (2022) Lu, Jun and Ye, Xuanyu.  Flexible and hierarchical prior for Bayesian nonnegative matrix factorization.  _arXiv preprint arXiv:2205.11025_ , 2022.
  * Marlin (2003) Marlin, Benjamin M.  Modeling user rating profiles for collaborative filtering.  _Advances in neural information processing systems_ , 16, 2003.
  * Martinsson et al. (2011) Martinsson, Per-Gunnar, Rokhlin, Vladimir, and Tygert, Mark.  A randomized algorithm for the decomposition of matrices.  _Applied and Computational Harmonic Analysis_ , 30(1):47–68, 2011.
  * Mnih & Salakhutdinov (2007) Mnih, Andriy and Salakhutdinov, Russ R.  Probabilistic matrix factorization.  _Advances in neural information processing systems_ , 20, 2007.
  * Raiko et al. (2007) Raiko, Tapani, Ilin, Alexander, and Karhunen, Juha.  Principal component analysis for large scale problems with lots of missing values.  In _European Conference on Machine Learning_ , pp. 691–698. Springer, 2007.
  * Schmidt & Mohamed (2009) Schmidt, Mikkel N and Mohamed, Shakir.  Probabilistic non-negative tensor factorization using Markov chain Monte Carlo.  In _2009 17th European Signal Processing Conference_ , pp. 1918–1922. IEEE, 2009.
  * Seashore-Ludlow et al. (2015) Seashore-Ludlow, Brinton, Rees, Matthew G, Cheah, Jaime H, Cokol, Murat, Price, Edmund V, Coletti, Matthew E, Jones, Victor, Bodycombe, Nicole E, Soule, Christian K, Gould, Joshua, et al.  Harnessing connectivity in a large-scale small-molecule sensitivity dataset.  _Cancer discovery_ , 5(11):1210–1223, 2015.
  * Tipping & Bishop (1999) Tipping, Michael E and Bishop, Christopher M.  Probabilistic principal component analysis.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 61(3):611–622, 1999.
  * Wang et al. (2013) Wang, Jim Jing-Yan, Wang, Xiaolei, and Gao, Xin.  Non-negative matrix factorization by maximizing correntropy for cancer clustering.  _BMC bioinformatics_ , 14(1):1–11, 2013.

##  Appendix A Derivation of Gibbs Sampler for Bayesian Interpolative
Decomposition

In this section, we give the derivation of the Gibbs sampler for the proposed
Bayesian ID models. As shown in the main paper, the observed data $(m,n)$-th
entry $a_{mn}$ of matrix $\bm{A}$ is modeled using a Gaussian likelihood with
variance $\sigma^{2}$ and mean given by the latent decomposition
$\bm{x}_{m}^{\top}\bm{y}_{n}$ (Eq. (2)),

$\displaystyle p(a_{mn}|\bm{x}_{m}^{\top}\bm{y}_{n},\sigma^{2})$
$\displaystyle=\mathcal{N}(a_{mn}|\bm{x}_{m}^{\top}\bm{y}_{n},\sigma^{2});$
$\displaystyle
p(\bm{A}|\bm{\theta})=\prod_{m,n=1}^{M,N}\mathcal{N}\left(a_{mn}|(\bm{X}\bm{Y})_{mn},\sigma^{2}\right)$
$\displaystyle=\prod_{m,n=1}^{M,N}\mathcal{N}\left(a_{mn}|(\bm{X}\bm{Y})_{mn},\tau^{-1}\right),$

where $\bm{\theta}=\\{\bm{X},\bm{Y},\sigma^{2}\\}$ denotes all parameters in
the model, $\sigma^{2}$ is the variance, $\tau^{-1}=\sigma^{2}$ is the
precision,

#### Prior

We choose a conjugate prior on the data variance, an inverse-Gamma
distribution with shape $\alpha_{\sigma}$ and scale $\beta_{\sigma}$,

$p(\sigma^{2}|\alpha_{\sigma},\beta_{\sigma})=\mathcal{G}^{-1}(\sigma^{2}|\alpha_{\sigma},\beta_{\sigma}).$

We assume further that the latent variable $y_{kl}$’s are drawn from a
general-truncated-normal prior,

$p(y_{kl}|\cdot)=\mathcal{GTN}(y_{kl}|\mu_{kl},\tau_{kl}^{-1},a=-1,b=1)=\frac{\sqrt{\frac{\tau_{kl}}{2\pi}}\exp\\{-\frac{\tau_{kl}}{2}(y_{kl}-\mu_{kl})^{2}\\}}{\Phi((b-\mu_{kl})\cdot\sqrt{\tau_{kl}})-\Phi((a-\mu_{kl})\cdot\sqrt{\tau_{kl}})}u(y_{kl}|a,b),$
(13)

where $u(x|a,b)$ is a step function with a value of 1 when $a\leq x\leq b$ and
0 otherwise.

Let $\bm{r}\in\\{0,1\\}^{N}$ be the state vector with each element indicating
the type of the corresponding column. If $r_{n}=1$, then $\bm{a}_{n}$ is a
basis column; otherwise, $\bm{a}_{n}$ is interpolated using the basis columns
plus some error term.

#### Hyperprior

To further favor flexibility, we choose a convenient joint hyperprior density
over the parameters $\\{\mu_{kl},\tau_{kl}\\}$ of GTN prior in Eq. (13),
namely, the GTN-scaled-normal-Gamma (GTNSNG) prior,

$\displaystyle\,\,\,\,\,\,\,\,p(\mu_{kl},\tau_{kl}|\cdot)=\mathcal{GTNSNG}(\mu_{kl},\tau_{kl}|\mu_{\mu},(\tau_{\mu})^{-1},\alpha_{t},\beta_{t})$
(14)
$\displaystyle=\big{\\{}\Phi((b-\mu_{\mu})\cdot\sqrt{\tau_{\mu}})-\Phi((a-\mu_{\mu})\cdot\sqrt{\tau_{\mu}})\big{\\}}\cdot\mathcal{N}(\mu_{kl}|\mu_{\mu},(\tau_{\mu})^{-1})\cdot\mathcal{G}(\tau_{kl}|a_{t},b_{t}).$

This prior can decouple parameters $\mu_{kl},\tau_{kl}$, and the posterior
conditional densities of them are normal and Gamma respectively due to this
convenient scale.

#### Posterior

Following the graphical representation of the GBT (or the GBTN) model in
Figure 2, the conditional density of $y_{kl}^{W}$ can be derived,

$\displaystyle\,\,\,\,\,\,\,\,p(y_{kl}|\bm{A},\bm{X},\bm{Y}_{-kl},\mu_{kl},\tau_{kl},\sigma^{2})\propto
p(\bm{A}|\bm{X},\bm{Y},\sigma^{2})\cdot p(y_{kl}|\mu_{kl},\tau_{kl})$
$\displaystyle=\prod_{i,j=1}^{M,N}\mathcal{N}\left(a_{ij}|\bm{x}_{i}^{\top}\bm{y}_{j},\sigma^{2}\right)\times\mathcal{GTN}(y_{kl}|\mu_{kl},(\tau_{kl})^{-1},a=-1,b=1)$
$\displaystyle\propto\exp\left\\{-\frac{1}{2\sigma^{2}}\sum_{i,j=1}^{M,N}(a_{ij}-\bm{x}_{i}^{\top}\bm{y}_{j})^{2}\right\\}\exp\\{-\frac{\tau_{kl}}{2}(y_{kl}-\mu_{kl})^{2}\\}u(y_{kl}|a,b)$
$\displaystyle\propto\exp\left\\{-\frac{1}{2\sigma^{2}}\sum_{i}^{M}(a_{il}-\bm{x}_{i}^{\top}\bm{y}_{l})^{2}\right\\}\exp\\{-\frac{\tau}{2}(y_{kl}-\mu_{kl})^{2}\\}u(y_{kl}|a,b)$
$\displaystyle\propto\exp\left\\{-\frac{1}{2\sigma^{2}}\sum_{i}^{M}\bigg{(}x_{ik}^{2}y_{kl}^{2}+2x_{ik}y_{kl}(\sum_{j\neq
k}^{N}x_{ij}y_{jl}-a_{il})\bigg{)}\right\\}\exp\\{-\frac{\tau_{kl}}{2}(y_{kl}-\mu_{kl})^{2}\\}u(y_{kl}|a,b)$
$\displaystyle\propto\exp\left\\{-y_{kl}^{2}\big{(}\frac{\sum_{i}^{M}x_{ik}^{2}}{2\sigma^{2}}+{\color[rgb]{0,0,0}\frac{\tau_{kl}}{2}}\big{)}+y_{kl}\underbrace{\bigg{(}\frac{1}{\sigma^{2}}\sum_{i}^{M}x_{ik}\big{(}a_{il}-\sum_{j\neq
k}^{N}x_{ij}y_{jl}\big{)}+{\color[rgb]{0,0,0}\tau_{kl}\mu_{kl}}\bigg{)}}_{{\color[rgb]{0,0,0}\widetilde{\tau}\widetilde{\mu}}}\right\\}u(y_{kl}|a,b)$
$\displaystyle\propto\mathcal{N}(y_{kl}|\widetilde{\mu},(\widetilde{\tau})^{-1})u(y_{kl}|a,b)\propto\mathcal{GTN}(y_{kl}|\widetilde{\mu},(\widetilde{\tau})^{-1},a=-1,b=1),$

where again, for simplicity, we assume the rows of $\bm{X}$ are denoted by
$\bm{x}_{i}$’s and columns of $\bm{Y}$ are denoted by $\bm{y}_{j}$’s,
$\widetilde{\tau}=\frac{\sum_{i}^{M}x_{ik}^{2}}{\sigma^{2}}+\tau_{kl}$ is the
posterior “parent precision” of the general-truncated-normal distribution, and

$\widetilde{\mu}=\bigg{(}\frac{1}{\sigma^{2}}\sum_{i}^{M}x_{ik}\big{(}a_{il}-\sum_{j\neq
k}^{N}x_{ij}y_{jl}\big{)}+{\color[rgb]{0,0,0}\tau_{kl}\mu_{kl}}\bigg{)}\big{/}\widetilde{\tau}$

is the posterior “parent mean” of the general-truncated-normal distribution.

Given the state vector $\bm{r}=[r_{1},r_{2},\ldots,r_{N}]^{\top}\in{}^{N}$,
the relation between $\bm{r}$ and the index sets $J$ is simple;
$J=J(\bm{r})=\\{n|r_{n}=1\\}_{n=1}^{N}$ and
$I=I(\bm{r})=\\{n|r_{n}=0\\}_{n=1}^{N}$. A new value of state vector $\bm{r}$
is to select one index $j$ from index sets $J$ and another index $i$ from
index sets $I$ (we note that $r_{j}=1$ and $r_{i}=0$ for the old state vector
$\bm{r}$) such that

$o_{j}=\frac{p(r_{j}=0,r_{i}=1|\bm{A},\sigma^{2},\bm{Y},\bm{r}_{-ji})}{p(r_{j}=1,r_{i}=0|\bm{A},\sigma^{2},\bm{Y},\bm{r}_{-ji})}=\frac{p(r_{j}=0,r_{i}=1)}{p(r_{j}=1,r_{i}=0)}\times\frac{p(\bm{A}|\sigma^{2},\bm{Y},\bm{r}_{-ji},r_{j}=0,r_{i}=1)}{p(\bm{A}|\sigma^{2},\bm{Y},\bm{r}_{-ji},r_{j}=1,r_{i}=0)},$

where $\bm{r}_{-ji}$ denotes all elements of $\bm{r}$ except $j$-th and $i$-th
entries. Trivially, we can set $p(r_{j}=0,r_{i}=1)=p(r_{j}=1,r_{i}=0)$. Then
the full conditionally probability of
$p(r_{j}=0,r_{i}=1|\bm{A},\sigma^{2},\bm{Y},\bm{r}_{-ji})$ can be calculated
by

$p(r_{j}=0,r_{i}=1|\bm{A},\sigma^{2},\bm{Y},\bm{r}_{-ji})=\frac{o_{j}}{1+o_{j}}.$

Finally, the conditional density of $\sigma^{2}$ is an inverse-Gamma
distribution by conjugacy,

$\displaystyle\,\,\,\,\,\,\,\,p(\sigma^{2}|\bm{X},\bm{Y},\bm{A})=\mathcal{G}^{-1}(\sigma^{2}|\widetilde{\alpha_{\sigma}},\widetilde{\beta_{\sigma}}),$

where $\widetilde{\alpha_{\sigma}}=\frac{MN}{2}+\alpha_{\sigma}$,
$\widetilde{\beta_{\sigma}}=\frac{1}{2}\sum_{i,j=1}^{M,N}(a_{ij}-\bm{x}_{i}^{\top}\bm{y}_{j})^{2}+\beta_{\sigma}$.

#### Extra update for GBTN model

Following the graphical representation of the GBTN model in Figure 2, the
conditional density of $\mu_{kl}$ can be derived,

$\displaystyle\,\,\,\,\,\,\,\,p(\mu_{kl}|\tau_{kl},\mu_{\mu},\tau_{\mu},a_{t},b_{t},y_{kl})\propto\mathcal{GTN}(y_{kl}|\mu_{kl},(\tau_{kl})^{-1},a=-1,b=1)\cdot\mathcal{GTNSNG}(\mu_{kl},\tau_{kl}|\mu_{\mu},(\tau_{\mu})^{-1},\alpha_{t},\beta_{t})$
$\displaystyle\propto\mathcal{GTN}(y_{kl}|\mu_{kl},(\tau_{kl})^{-1},a=-1,b=1)\cdot\big{\\{}\Phi((b-\mu_{\mu})\cdot\sqrt{\tau_{\mu}})-\Phi((a-\mu_{\mu})\cdot\sqrt{\tau_{\mu}})\big{\\}}\cdot{\mathcal{N}(\mu_{kl}|\mu_{\mu},(\tau_{\mu})^{-1})}\cdot\cancel{\mathcal{G}(\tau_{kl}|a_{t},b_{t})}$
$\displaystyle\propto\sqrt{\tau_{kl}}\cdot\exp\left\\{-\frac{\tau_{kl}}{2}(y_{kl}-\mu_{kl})^{2}\right\\}\cdot\exp\left\\{-\frac{\tau_{\mu}}{2}(\mu_{\mu}-\mu_{kl})^{2}\right\\}$
$\displaystyle\propto\exp\left\\{-\frac{\tau_{kl}+\tau_{\mu}}{2}\mu_{kl}^{2}+\mu_{kl}\underbrace{(\tau_{kl}y_{kl}+\tau_{\mu}\mu_{\mu})}_{\widetilde{m}\cdot\widetilde{t}}\right\\}\propto\mathcal{N}(\mu_{kl}|\widetilde{m},(\,\widetilde{t}\,)^{-1}),$

where $\widetilde{t}=\tau_{kl}+\tau_{\mu}$, and
$\widetilde{m}=(\tau_{kl}y_{kl}+\tau_{\mu}\mu_{\mu})/\widetilde{t}$ are the
posterior precision and mean of the normal density. Similarly, the conditional
density of $\tau_{kl}$ is,

$\displaystyle\,\,\,\,\,\,\,\,p(\tau_{kl}|\mu_{kl},\mu_{\mu},\tau_{\mu},a_{t},b_{t},y_{kl})\propto\mathcal{GTN}(y_{kl}|\mu_{kl},(\tau_{kl})^{-1},a=-1,b=1)\cdot\mathcal{GTNSNG}(\mu_{kl},\tau_{kl}|\mu_{\mu},(\tau_{\mu})^{-1},\alpha_{t},\beta_{t})$
$\displaystyle\propto\mathcal{GTN}(y_{kl}|\mu_{kl},(\tau_{kl})^{-1},a=-1,b=1)\cdot\big{\\{}\Phi((b-\mu_{\mu})\cdot\sqrt{\tau_{\mu}})-\Phi((a-\mu_{\mu})\cdot\sqrt{\tau_{\mu}})\big{\\}}\cdot\cancel{\mathcal{N}(\mu_{kl}|\mu_{\mu},(\tau_{\mu})^{-1})}\cdot{\mathcal{G}(\tau_{kl}|a_{t},b_{t})}$
$\displaystyle\propto\exp\left\\{-\tau_{kl}\frac{(y_{kl}-\mu_{kl})^{2}}{2}\right\\}\tau_{kl}^{1/2}\tau_{kl}^{a_{t}-1}\exp\left\\{-b_{t}\tau_{kl}\right\\}\propto\exp\left\\{-\tau_{kl}\left[b_{t}+\frac{(y_{kl}-\mu_{kl})^{2}}{2}\right]\right\\}\cdot\tau_{kl}^{(a_{t}+1/2)-1}$
$\displaystyle\propto\mathcal{G}(\tau_{kl}|\widetilde{a},\widetilde{b}),$

where $\widetilde{a}=a_{t}+1/2$ and
$\widetilde{b}=b_{t}+\frac{(y_{kl}-\mu_{kl})^{2}}{2}$ are the posterior
parameters of the Gamma density.