# Wideband Effects on Near-Field Pose Estimation of Target-Lodged RIS
††thanks: This work was partially supported by the European Union under the
Italian National Recovery and Resilience Plan (NRRP) of NextGenerationEU,
partnership on “Telecommunications of the Future” (PE00000001 - program
“RESTART”).

Dario Tagliaferri

###### Abstract

Reconfigurable intelligent surfaces (RISs) have been recently considered in
sensing context to either improve localization performance or to extend the
coverage in non-line-of-sight scenarios. While most of the literature
considers RISs as statically placed in the environment, the usage of target-
lodged RISs is relatively new and could be of interest where the target’s
reflectivity can/must be increased to improve its detection or parameters’
estimation. This letter derives the Cramér-Rao bound (CRB) on the estimation
of position and orientation (pose) of a target-mounted RIS, in generic
conditions: near-field, bistatic and wideband operation (i.e., when the
wavefront across the RIS is curved and the employed sensing bandwidth is large
enough to obtain a frequency-dependent RIS behavior). In particular, we focus
on the wideband effect, that implies a pose-dependent filtering on the
impinging signal decreasing or increasing the CRB depending on the RIS size
and the employed signal bandwidth.

###### Index Terms:

Cramér-Rao bound, near-field, wideband, position, orientation, RIS

##  I Introduction

Electromagnetic (EM) metasurfaces are rapidly surging in both academia and
industry as one of the promising technologies for the smart propagation in the
next generation of communication systems (6G) [1]. In particular,
reconfigurable intelligent surfaces (RISs) are discrete metasurfaces made of
passive regularly spaced, sub-wavelength-sized elements, whose reflection
coefficient can be dynamically tuned to manipulate the direction of reflection
and refraction of an impinging wave [2].

The usage of RIS for sensing is relatively recent [3]. In [4], the authors
proposed to use a RIS to assist radar in non-line-of-sight conditions,
focusing on the maximization of the probability of correct detection of a
given target, via signal-to-noise (SNR) maximization with constrained phase
design. The results analyze the Cramér-Rao bound (CRB) on position estimation
and suggest the placement of the RIS nearby either the radar or the target.
The work [5] considers indoor localization by means of wall-placed RIS. The
paper [6] compares discrete and continuous metasurfaces in terms of CRB for
position estimation in the integrated sensing and communication context. With
the advent of large RISs, that enable ultra-high precision localization, the
near-field (NF) propagation condition applies, i.e., the wavefront across the
RIS is no more flat and allows proper energy focusing through reflection. NF
operation opened for further research on fundamental localization bounds and
practical RIS phase configuration. The work [7] derives the CRB for position
and orientation estimation in some selected practical cases, assuming NF
operating conditions and discussing the phase configuration. The authors of
[8] propose a practical NF RIS phase configuration and analyze the CRB on
position estimation.

In all the the aforementioned literature, RISs are statically placed in the
environment and used to aid localization, by magnifying the Tx/Rx aperture to
decrease the CRB. In some scenarios, however, the RIS can be directly lodged
on targets, to enhance their "visibility" in the sensing data and improve
their localization performance. This setup is only considered in [9] and [10].
In the former, the authors investigate the CRB on position and orientation
estimation for a target-mounted RIS under far-field (FF) assumption, while the
latter addresses the CRB evaluation on position estimation with a perfect RIS
phase configuration and bistatic settings. However, none of the previous works
evaluate the wideband effects, namely the change of the RIS response due to a
change in the frequency within the employed signal bandwidth. The only work
addressing this issue is [11], for a single-input-single output system and RIS
again placed in the environment. In FF conditions, the wideband effect gives
rise to reflection beam squinting, i.e, a frequency-dependent spatial
spreading of the reflected signal around the desired direction diminishing the
amount of received energy. Similarly, a RIS under NF and wideband operation is
subject to defocusing.

This letter tackles the problem of pose (3D position plus 3D orientation)
estimation of a target-lodged RIS, in a generic bistatic, NF and wideband
(frequency-selective) setting. We derive the CRB for perfect RIS phase
configuration (although imperfect phase configuration can be considered as
well). Further, we outline the importance of considering the wideband RIS
behavior in practical cases, where the RIS acts as a pose-dependent filter for
the impinging sensing signal. Remarkably, the latter effect allows increasing
or reducing the available information w.r.t. the narrowband case (where the
RIS is assumed to behave like a frequency-independent target) depending on the
RIS size and the specific employed bandwidth value.

The letter is organized as follows: Sect. II presents the system model, in
Sect. III we evaluate the CRB, Sect. IV presents the numerical results while
Sect. V concludes the paper. We adopt the following notation: bold upper- and
lower-case denote matrices and column vectors. The L2-norm of a vector is
denoted with $\|\cdot\|$. Matrix transposition and conjugate transposition are
indicated respectively as $\mathbf{A}^{T}$ and $\mathbf{A}^{H}$.
$\mathbf{I}_{n}$ is the identity matrix of size $n$.
$\mathbf{a}\sim\mathcal{CN}(\boldsymbol{\mu},\mathbf{C})$ denotes a multi-
variate circularly complex Gaussian random variable with mean
$\boldsymbol{\mu}$ and covariance $\mathbf{C}$. $\mathrm{trace}(\mathbf{A})$
extracts the trace of $\mathbf{A}$. $\mathbb{R}$ and $\mathbb{C}$ stand for
the set of real and complex numbers, respectively. $\delta_{n}$ is the
Kronecker delta.

##  II System Model

Figure 1: System model.

Let us consider a generic sensing terminal in a bistatic setting (e.g., a
multiple-input multiple-output (MIMO) radar) that aims at estimating the
position and orientation of a target equipped with a RIS (Fig. 1). In a global
reference system, the phase centers of transmitting (Tx) and receiving (Rx)
antenna arrays are located in $\mathbf{x}_{\mathrm{T}}\in\mathbb{R}^{3\times
1}$ and $\mathbf{x}_{\mathrm{R}}\in\mathbb{R}^{3\times 1}$, respectively. The
sensing terminal implements $L$ measurement channels, for all Tx-Rx antenna
pairs (say, each single antennas is located in $\mathbf{x}_{\mathrm{T},\ell}$
and $\mathbf{x}_{\mathrm{R},\ell}$, respectively). At each Tx antenna, the
sensing terminal emits the signal

$s(t)=g(t)e^{j2\pi f_{0}t}$ (1)

where the base-band waveform $g(t)$ has bandwidth $B$ and the carrier
frequency is $f_{0}$. For simplicity, we assume $g(t)$ to be the same on each
Tx antenna, as for a MIMO radar operating in time-division multiplexing. The
RIS is in the generic position $\mathbf{x}\in\mathbb{R}^{3\times 1}$ and is
made by $N\times M$ elements, displaced along the $x$ and $y$ axes of the of
the local reference system. The local axes are rotated w.r.t. the global ones
by the Euler angles
$\boldsymbol{\gamma}=[\psi_{x},\psi_{y},\psi_{z}]^{T}\in\mathbb{R}^{3\times
1}$ (roll $\psi_{x}$, pitch $\psi_{y}$ and yaw $\psi_{z}$). These latter
represent the RIS orientation, that is the objective of the estimation
together with position $\mathbf{x}$. The position of the $(n,m)$-th element of
the RIS in global coordinates is

$\displaystyle\mathbf{x}_{nm}$
$\displaystyle=\mathbf{x}+\mathbf{Q}(\boldsymbol{\gamma})\mathbf{p}_{nm}$ (2)

for $n=-N/2,...,N/2-1$, $m=-M/2,...,M/2-1$, where
$\mathbf{p}_{nm}=\left[nd,\,md,\,0\right]^{\mathrm{T}}$ is the local position
of the $(n,m)$-th element in local coordinates (for inter-element spacing
$d$), matrix
$\mathbf{Q}(\boldsymbol{\gamma})\in\mathrm{SO}(3)\overset{\triangle}{=}\left\\{\mathbf{Q}|\mathrm{det}(\mathbf{Q})=1,\mathbf{Q}\mathbf{Q}^{T}=\mathbf{I}_{3}\right\\}$
defines the relative counterclockwise rotation of the local reference system
around $z$, $y$ and $x$ axis by angles $\psi_{z}$, $\psi_{y}$, $\psi_{x}$,
respectively. The expression of $\mathbf{Q}(\boldsymbol{\gamma})$ is:

$\mathbf{Q}(\boldsymbol{\gamma})=\mathbf{Q}_{z}(\psi_{z})\mathbf{Q}_{y}(\psi_{y})\mathbf{Q}_{x}(\psi_{x})$
(3)

where single cascaded rotation matrices are defined according to the
convention indicated in Fig. 1.

The Rx signal for the $\ell$-th measurement channel of the sensing terminal
(after demodulation) is reported in (4),

$\begin{split}y_{\ell}(t)&=\rho\sum_{n,m}e^{j\Phi_{nm}}e^{-j2\pi
f_{0}(\tau^{i}_{\ell,nm}+\tau^{o}_{\ell,nm})}g\left(t-\tau^{i}_{\ell,nm}-\tau^{o}_{\ell,nm}\right)+z_{\ell}(t)\\\
&\overset{(a)}{\approx}\rho\,e^{-j2\pi
f_{0}(\tau^{i}_{0}+\tau^{o}_{0}+\Delta\tau^{i}_{\ell}+\Delta\tau^{o}_{\ell})}\sum_{n,m}e^{j\Phi_{nm}}e^{-j2\pi
f_{0}(\Delta\tau^{i}_{nm}+\Delta\tau^{o}_{nm})}g\left(t-\tau^{i}_{0}-\tau^{o}_{0}-\Delta\tau^{i}_{\ell}-\Delta\tau^{o}_{\ell}-\Delta\tau^{i}_{nm}-\Delta\tau^{o}_{nm}\right)+z_{\ell}(t)\\\
&\overset{(b)}{\approx}\rho\,g(t-\tau^{i}_{0}-\tau^{o}_{0})\,e^{-j2\pi
f_{0}(\tau^{i}_{0}+\tau^{o}_{0}+\Delta\tau^{i}_{\ell}+\Delta\tau^{o}_{\ell})}\sum_{n,m}e^{j\Phi_{nm}}e^{-j2\pi
f_{0}(\Delta\tau^{i}_{nm}+\Delta\tau^{o}_{nm})}+z_{\ell}(t)\end{split}$ (4)

$\begin{split}Y_{\ell}(f)&=\rho\,G(f)\sum_{n,m}e^{j\Phi_{nm}}e^{-j2\pi(f_{0}+f)(\tau^{i}_{\ell,nm}+\tau^{o}_{\ell,nm})}+Z_{\ell}(f)\\\
&\overset{(a)}{\approx}G(f)e^{-j2\pi(f_{0}+f)(\tau^{i}_{0}+\tau^{o}_{0}+\Delta\tau^{i}_{\ell}+\Delta\tau^{o}_{\ell})}\beta(f|\mathbf{\Phi})+Z_{\ell}(f)\overset{(b)}{\approx}G(f)e^{-j2\pi(f_{0}+f)(\tau^{i}_{0}+\tau^{o}_{0})}e^{-j2\pi
f_{0}(\Delta\tau^{i}_{\ell}+\Delta\tau^{o}_{\ell})}\beta(0|\mathbf{\Phi})+Z_{\ell}(f)\end{split}$
(5)

where $\rho$ denotes geometrical energy losses, $\Phi_{nm}$ is the phase
applied at the $(n,m)$-th element of the RIS,

$\tau^{i}_{\ell,nm}=\frac{\|\mathbf{x}_{nm}-\mathbf{x}_{\mathrm{T},\ell}\|}{c},\,\,\,\,\,\tau^{o}_{\ell,nm}=\frac{\|\mathbf{x}_{\mathrm{R},\ell}-\mathbf{x}_{nm}\|}{c},$
(6)

are the absolute delays between the $\ell$-th Tx and Rx antennas and the
$(n,m)$-th element of the RIS, respectively. Term
$z_{\ell}(t)\in\mathcal{CN}(0,\sigma_{z}^{2}\delta_{\ell-k}\delta(t))$ is the
additive noise, assumed as white and uncorrelated over different channels. The
expression of $\rho$ follows from the radar equation [12]:

$\rho=\sqrt{\frac{c^{2}}{(4\pi)^{3}f^{2}_{0}\|\mathbf{x}-\mathbf{x}_{\mathrm{T}}\|^{2}\|\mathbf{x}_{\mathrm{R}}-\mathbf{x}\|^{2}}\Gamma_{\mathrm{elem}}}\;e^{j\delta}$
(7)

where $\Gamma_{\mathrm{elem}}$ is the radar cross section of the single RIS
element of size $d^{2}$, and $\delta$ models residual phase uncertainties
about the target (the RIS) arising, for instance, from Tx-Rx circuitry and
Doppler effects from motion. We assume that Tx and Rx antennas, as well as RIS
elements, are isotropic for simplicity, although more rigorous single-element
models can be used [13]. The overall path-loss depends only on the macroscopic
Tx-RIS and RIS-Rx distances.

Approximation $(a)$ in (4) assumes FF at both RIS and sensing terminal, i.e,
uniform planar wavefronts, thus the delays can be linearized as

$\displaystyle\tau^{i}_{\ell,nm}$
$\displaystyle\overset{(a)}{\simeq}\underbrace{\frac{\|\mathbf{x}-\mathbf{x}_{\mathrm{T}}\|}{c}}_{\tau^{i}_{0}}+\underbrace{\frac{\mathbf{x}_{\mathrm{T},\ell}^{T}\mathbf{u}(\boldsymbol{\zeta}_{\mathrm{T}})}{c}}_{\Delta\tau^{i}_{\ell}}+\underbrace{\frac{\mathbf{p}^{T}_{nm}\mathbf{u}(\boldsymbol{\xi}_{i})}{c}}_{\Delta\tau^{i}_{nm}}$
(8) $\displaystyle\tau^{o}_{\ell,nm}$
$\displaystyle\overset{(a)}{\simeq}\underbrace{\frac{\|\mathbf{x}_{\mathrm{R}}-\mathbf{x}\|}{c}}_{\tau^{o}_{0}}+\underbrace{\frac{\mathbf{x}_{\mathrm{R},\ell}^{T}\mathbf{u}(\boldsymbol{\zeta}_{\mathrm{R}})}{c}}_{\Delta\tau^{o}_{\ell}}+\underbrace{\frac{\mathbf{p}^{T}_{nm}\mathbf{u}(\boldsymbol{\xi}_{o})}{c}}_{\Delta\tau^{o}_{nm}},$
(9)

where $\tau^{i}_{0}$ and $\tau^{o}_{0}$ are the macroscopic delays between the
Tx and Rx phase centers and the RIS phase center, respectively, while
$\Delta\tau^{i}_{\ell}$ and $\Delta\tau^{o}_{\ell}$ are the excess delays at
the sensing terminal, function of angles of departure and arrivals
$\boldsymbol{\zeta}_{\mathrm{T}}$ and $\boldsymbol{\zeta}_{\mathrm{R}}$
respectively, whereas $\Delta\tau^{i}_{nm}$ and $\Delta\tau^{o}_{nm}$ are the
excess delays at the RIS, function of incidence and reflection angles
$\boldsymbol{\xi}_{i}=[\phi_{i},\theta_{i}]^{T}$ and
$\boldsymbol{\xi}_{o}=[\phi_{o},\theta_{o}]^{T}$ (elevation and azimuth). The
latter angles are computed as:

$\boldsymbol{\xi}_{i}=J\left(\mathbf{Q}(\boldsymbol{\gamma})(\mathbf{x}_{\mathrm{T}}-\mathbf{x})\right),\,\,\,\boldsymbol{\xi}_{o}=J\left(\mathbf{Q}(\boldsymbol{\gamma})(\mathbf{x}_{\mathrm{R}}-\mathbf{x})\right)$
(10)

where $J:\mathbb{R}^{3\times 1}\rightarrow\mathbb{R}_{2\pi}^{2\times 1}$ is
the transformation between global Cartesian coordinates to RIS-local spherical
coordinates, according to the convention in Fig. 1. Unit vector
$\mathbf{u}(\cdot)\in\mathbb{R}^{3\times 1}$ defines the direction in local or
global coordinates. For the planar RIS considered in this work, the excess
delays are:

$\displaystyle\Delta\tau^{i}_{nm}$
$\displaystyle=\frac{d}{c}\left(n\sin\phi_{i}\cos\theta_{i}+m\sin\phi_{i}\sin\theta_{i}\right)$
(11) $\displaystyle\Delta\tau^{o}_{nm}$
$\displaystyle=\frac{d}{c}\left(n\sin\phi_{o}\cos\theta_{o}+m\sin\phi_{o}\sin\theta_{o}\right).$
(12)

The further approximation $(b)$ in (4) assumes that

$g\left(t-\tau^{i}_{0}-\tau^{o}_{0}-\Delta\tau^{i}_{\ell}-\Delta\tau^{o}_{\ell}-\Delta\tau^{i}_{nm}-\Delta\tau^{o}_{nm}\right)\overset{(b)}{\approx}g\left(t-\tau^{i}_{0}-\tau^{o}_{0}\right)$
(13)

i.e., the base-band Tx waveform is not affected by the residual (excess)
delays at the RIS and at the sensing terminal. The narrowband operation
described by (13) is currently assumed in all the literature [9, 10] but does
not practically match the reality in most of sensing acquisitions, where the
employed bandwidth $B$ is large enough to induce reflection beam squinting
from the RIS (i.e., a frequency-dependent RIS reflection coefficient). To gain
insight on the latter phenomenon, we can analyze the Rx signal in the
frequency domain, as reported in (5). Here, $G(f)$ is the Fourier transform of
$g(t)$,

$\beta(f|\mathbf{\Phi})=\rho\,\sum_{n,m}e^{j\Phi_{nm}}e^{-j2\pi(f_{0}+f)(\Delta\tau^{i}_{nm}+\Delta\tau^{o}_{nm})}$
(14)

is the overall frequency-dependent RIS reflection coefficient (including path-
loss) depending on the phase configuration
$\mathbf{\Phi}\in\mathbb{C}^{NM\times 1}$, and
$Z_{\ell}(f)\in\mathcal{CN}(0,N_{0}\delta_{\ell-k}\delta(f))$ is the noise in
the frequency domain, with power spectral density $N_{0}$. The base-band
frequency $f\in[-B/2,B/2]$ models the wideband response of both Tx/Rx arrays
and the RIS (e.g., see [14] for further details) through
$\beta(f|\boldsymbol{\Phi})$, that filters the impinging signal $G(f)$. Notice
that the narrowband approximation in the frequency domain implies that the RIS
has the same reflection coefficient for all the frequencies $f$ (corresponding
to $f_{0}$). In the following, we derive and compare the CRB for the model
(5), highlighting the wideband effects.

##  III CRB calculation

The model (5) for the $\ell$-th measurement channel can be stacked into a
$L\times 1$ vector:

$\mathbf{y}(f)=\mathbf{a}(f,\boldsymbol{\theta}|\overline{\mathbf{\Phi}})+\mathbf{z}(f)$
(15)

where $\mathbf{a}(f,\boldsymbol{\theta}|\overline{\mathbf{\Phi}})$ is the non-
linear model relating the parameters to be estimated
$\boldsymbol{\theta}=[\mathbf{x}^{T},\boldsymbol{\gamma}^{T}]^{T}\in\mathbb{R}^{6\times
1}$ with the observation, conditioned to the specific phase configuration of
the RIS $\overline{\mathbf{\Phi}}$, while
$\mathbf{z}(f)\sim\mathcal{CN}(\mathbf{0},N_{0}\mathbf{I}_{L}\delta(f))$.
Notice that we assumed the knowledge of the reflection coefficient (including
path-loss) $\rho$, that allows obtaining an optimistic lower bound on pose
estimation. If $\rho$ is unknown, its real and imaginary parts can be included
in the parameters to be estimated. The RIS can be configured according to:

$[\overline{\mathbf{\Phi}}]_{nm}=\begin{dcases}2\pi
f_{0}\left[\frac{\|\mathbf{x}_{nm}-\mathbf{x}_{\mathrm{T}}\|}{c}+\frac{\|\mathbf{x}_{\mathrm{R}}-\mathbf{x}_{nm}\|}{c}\right]&\text{NF}\\\
2\pi
f_{0}\left[\frac{\mathbf{p}^{T}_{nm}\mathbf{u}(\boldsymbol{\xi}_{i})}{c}+\frac{\mathbf{p}^{T}_{nm}\mathbf{u}(\boldsymbol{\xi}_{o})}{c}\right]&\text{FF}\end{dcases}$
(16)

where in the first case (NF) the RIS is able to perfectly focus the impinging
radiation towards the phase center of the sensing terminal, while in the
second case (FF) the RIS is configured to reflect the signal from
$\boldsymbol{\xi}_{i}$ to $\boldsymbol{\xi}_{o}$. The perfect RIS
configuration represents an upper performance bound modeling the situation in
which the RIS follows some focusing alignment procedure such as the one in
[15], not necessarily implying the perfect knowledge of $\mathbf{x}$ and
$\boldsymbol{\gamma}$ (and thus of $\mathbf{x}_{nm}$), which is the objective
of the estimation. A similar consideration can be made for FF configuration.
In any case, the NF/FF phase configuration is optimal only at $f_{0}$.

Based on model (15), the Fisher information matrix (FIM) $\mathbf{F}$ is
block-partitioned as follows:

$\mathbf{F}=\begin{bmatrix}\mathbf{F}_{\mathbf{x}\mathbf{x}}&\mathbf{F}_{\mathbf{x}\boldsymbol{\gamma}}\\\
\mathbf{F}_{\boldsymbol{\gamma}\mathbf{x}}&\mathbf{F}_{\boldsymbol{\gamma}\boldsymbol{\gamma}}\end{bmatrix}$
(17)

with straightforward dimensions, whose entries are

$\displaystyle\mathbf{F}_{\boldsymbol{\mu}\boldsymbol{\nu}}=\frac{2}{N_{0}}\Re\left\\{\int\limits_{-B/2}^{B/2}\left(\frac{\partial\mathbf{a}(f,\boldsymbol{\theta}|\overline{\mathbf{\Phi}})}{\partial\boldsymbol{\mu}}\right)^{H}\frac{\partial\mathbf{a}(f,\boldsymbol{\theta}|\overline{\mathbf{\Phi}})}{\partial\boldsymbol{\nu}}df\right\\}$
(18)

where $\boldsymbol{\mu}$ and $\boldsymbol{\nu}$ can be any set of the
parameters to be estimated, i.e., $\mathbf{x}$, $\boldsymbol{\gamma}$. Single
FIM entries can be evaluated as in (19)-(20), where
$\boldsymbol{\tau}^{i}_{nm}=[\tau^{i}_{1,nm},...,\tau^{i}_{L,nm}]^{T}$,
$\boldsymbol{\tau}^{o}_{nm}=[\tau^{o}_{1,nm},...,\tau^{o}_{L,nm}]^{T}$ are the
delay vectors, whose gradient w.r.t. $\mathbf{x}$ and $\boldsymbol{\gamma}$ is
computed as (exemplary for $\boldsymbol{\tau}^{i}_{nm}$):

$\displaystyle\frac{\partial\mathbf{a}(f,\boldsymbol{\theta})}{\partial\mathbf{x}}$
$\displaystyle=\rho\,G(f)\sum_{n,m}e^{j\overline{\Phi}_{nm}}[-j2\pi(f_{0}+f)]e^{-j2\pi(f_{0}+f)(\boldsymbol{\tau}^{i}_{nm}+\boldsymbol{\tau}^{o}_{nm})}\frac{\partial}{\partial\mathbf{x}}\left(\boldsymbol{\tau}^{i}_{nm}+\boldsymbol{\tau}^{o}_{nm}\right)+$
$\displaystyle+G(f)\frac{\partial\rho}{\partial\mathbf{x}}\sum_{n,m}e^{j\overline{\Phi}_{nm}}e^{-j2\pi(f_{0}+f)(\boldsymbol{\tau}^{i}_{nm}+\boldsymbol{\tau}^{o}_{nm})}\in\mathbb{C}^{L\times
3}$ (19)
$\displaystyle\frac{\partial\mathbf{a}(f,\boldsymbol{\theta})}{\partial\boldsymbol{\gamma}}$
$\displaystyle=\rho\,G(f)\sum_{n,m}e^{j\overline{\Phi}_{nm}}[-j2\pi(f_{0}+f)]e^{-j2\pi(f_{0}+f)(\boldsymbol{\tau}^{i}_{nm}+\boldsymbol{\tau}^{o}_{nm})}\frac{\partial}{\partial\boldsymbol{\gamma}}\left(\boldsymbol{\tau}^{i}_{nm}+\boldsymbol{\tau}^{o}_{nm}\right)\in\mathbb{C}^{L\times
3}$ (20)

$\displaystyle\left[\frac{\partial\boldsymbol{\tau}^{i}_{nm}}{\partial\mathbf{x}}\right]_{\ell}$
$\displaystyle=\frac{1}{c}\frac{(\mathbf{x}_{nm}-\mathbf{x}_{\mathrm{T},\ell})^{T}}{\|\mathbf{x}_{nm}-\mathbf{x}_{\mathrm{T},\ell}\|}$
(21)
$\displaystyle\left[\frac{\partial\boldsymbol{\tau}^{i}_{nm}}{\partial\boldsymbol{\gamma}}\right]_{\ell}$
$\displaystyle=\frac{1}{c}\frac{(\mathbf{x}_{nm}-\mathbf{x}_{\mathrm{T},\ell})^{T}}{\|\mathbf{x}_{nm}-\mathbf{x}_{\mathrm{T},\ell}\|}\frac{\partial\left(\mathbf{Q}(\boldsymbol{\gamma})\mathbf{p}_{nm}\right)}{\partial\boldsymbol{\gamma}}.$
(22)

The components of the FIM under the FF approximation (and possibly narrowband
settings) follow from (19)-(20) under the approximation of the delays in
(8)-(9) and condition (13).

##  IV Numerical Results

(a)

(b)

Figure 2: PEB vs. bandwidth $B$ for (a) $10\times 10$ cm2 RIS and (b)
$20\times 20$ cm2 RIS.  Figure 3: OEB of $\psi_{y}$ vs. RIS size (side size).

This section shows numerical results quantifying the wideband effects for the
pose estimation of a target-lodged RIS. We consider a monostatic sensing
terminal operating at $f_{0}=78.5$ GHz, on a variable bandwidth $B\in[1,10]$
GHz, equipped with $1$ Tx antenna located in the origin of the global
coordinate system and $20\times 20$ Rx antennas along $y$ and $z$ axes. The
RIS size is variable, from $A_{\mathrm{ris}}=5\times 5$ cm2 to
$A_{\mathrm{ris}}=35\times 35$ cm2. In all the evaluations, the RIS is in
$\mathbf{x}=[5,0,-5.5]^{T}$ m, located on the $xy$ plane, that may represent a
situation where the RIS is on a vehicle rooftop illuminated by a radar on a
road side unit above. The Tx power is $23$ dBm at the single Tx antenna, and
the Rx signal is corrupted by thermal noise, $N_{0}=-173$ dBm/Hz. Notice that,
as the power is kept fixed, increasing the bandwidth $B$ means decreasing the
energy spectral density of the Tx signal $G(f)$.

The first set of results is summarized in Fig. 2. We show the position error
bound (PEB), evaluated as

$\displaystyle\mathrm{PEB}=\sqrt{\frac{\mathrm{trace([\mathbf{F}^{-1}]_{1:3,1:3})}}{3}}$
(23)

where the FIM can be either the NF or the FF one, under either wideband or
narrowband operation. The latter, being the current state of the art [9, 10],
assumes that the RIS behavior is independent on the base-band frequency $f$,
as for approximation $(b)$ in (4) and (5). We also report as a benchmark the
lower bound in which the orientation $\boldsymbol{\gamma}$ is known (or its
estimation is assumed to be decoupled from position), thus
$\mathbf{F}=\mathbf{F}_{\mathbf{xx}}$. Let us first consider the case of a
medium sized RIS, i.e., $10\times 10$ cm2 (Fig. 2a). For increasing bandwidth
$B$, the PEB increases as well, as long as the wideband modeling (both NF and
FF) is concerned. PEB for narrowband FF is invariant, and it lower bounds the
wideband one, underestimating the true PEB. The former effect is counter
intuitive, as increasing $B$ is expected to decrease the PEB thanks to an
increase in the sensing resolution. This not true, as the RIS in wideband
conditions operates as a filter on the Tx signal spectrum $G(f)$ that actually
decreases the reflected energy w.r.t. the narrowband (frequency-flat) case.
Notice that even when the RIS is assumed to be not frequency-selective (FF,
narrowband curve), the PEB does not improve with $B$, as the effective
bandwidth of $g(t)$ (1) is [16]:

$B_{\mathrm{eff}}^{2}=f_{0}^{2}+\frac{B^{2}}{12}\approx f_{0}^{2},$ (24)

largely dominated by the carrier component for the considered settings, and
almost insensitive to $B$ (for $B=10$ GHz and $f_{0}=78.5$ GHz, the first term
of (24) is 3 orders of magnitude higher than the second). Fig. 2a also
highlights the difference between the NF and FF modeling of the RIS behavior.
When the orientation $\boldsymbol{\gamma}$ is known, e.g., a RIS is deployed
on the rooftop of a vehicle, the FF shows higher PEB, as expected, as NF
enables energy focusing at the Rx side.

Interestingly, for known $\boldsymbol{\gamma}$, the wideband PEB can be lower
than the narrowband one. Thus, the information brought in the FIM by the pose-
dependent filtering operated by the RIS overcomes the energy loss due to beam
squinting and yields lower PEB. As $B$ increases, the energy loss becomes
dominant and wideband PEB starts to be larger than the narrowband one (i.e,,
the NF/FF wideband PEB curves cross the FF narrowband one). This latter
threshold effect depends on the size of the RIS and on its distance from the
sensing terminal. Although not reported for brevity, it can be shown that as
the latter increases, the information loss due to spatial energy spreading
dominates over the pose-dependent filtering effect, shifting the threshold
towards lower values of $B$. This phenomenon is also observed by increasing
the size of the RIS (see Fig. 2b). In this latter case, the PEB gap between NF
and FF—when $\boldsymbol{\gamma}$ is known—is higher, whereas again no
difference between NF and FF is observed for the pose estimation problem.

The second result, in Fig. 3, reports the orientation error bound (OEB) on the
pitch angle $\psi_{y}$ for the NF model, versus the RIS size, varying the
bandwidth $B=1,4,8$ GHz. The OEB is defined as

$\mathrm{OEB}_{y}=\sqrt{[\mathbf{F}^{-1}]_{5,5}}$ (25)

where the lower bound corresponding to known position $\mathbf{x}$ is obtained
by plugging $\mathbf{F}=\mathbf{F}_{\boldsymbol{\gamma}\boldsymbol{\gamma}}$
in (25). We select the pitch angle $\psi_{y}$ to be analyzed as it is the one
providing the lower OEB among $\psi_{x},\psi_{y},\psi_{z}$, thus the higher
model sensitivity. We notice that increasing the size of the RIS means
increasing its frequency-selectivity for a given bandwidth $B$ (narrowing its
reflection beamwidth and exacerbating the wideband effect), thus penalizing
the Rx signal energy. The result is a OEB that worsens with $B$ except for
very small RIS. For known $\mathbf{x}$, instead, this latter effect is only
observed above a given RIS size, approximately $20\times 20$ cm2 in the
considered settings. Below, the lower OEB value is attained by $B=4$ GHz, as a
trade-off between frequency-selectivity (dominating for large $B$) and a
beneficial pose-dependent filtering (dominating for small $B$). This is a
further confirmation of what observed on the PEB in Fig. 2, and underlines the
importance of the wideband modeling in the estimation of pose of a target-
lodged RIS. As a final remark, the frequency-dependent RIS behavior opens for
further research on the optimal spectral shaping of $G(f)$ to optimize (lower)
the PEB/OEB (not covered here).

##  V Conclusion

This letter analyzes the wideband effects on the pose (position+orientation)
estimation of a target-lodged RIS, namely the frequency-dependent (filter) RIS
behavior when the employed sensing signal bandwidth is large. We derive the
CRB in NF/FF and wideband/narrowband modeling assumptions, discussing the
trend of PEB and OEB varying the employed bandwidth and RIS size. The results
highlight the importance of a wideband modeling of the RIS in the pose
estimation problem, as both PEB and OEB worsen as the bandwidth increases.
Under some circumstances and typical NF conditions, the wideband PEB is
tighter than the narrowband one, as the RIS behaves like a pose-dependent
filter.

## References

  * [1] M. Di Renzo, A. Zappone, M. Debbah, M.-S. Alouini, C. Yuen, J. De Rosny, and S. Tretyakov, “Smart radio environments empowered by reconfigurable intelligent surfaces: How it works, state of research, and the road ahead,” _IEEE Journal on Selected Areas in Communications_ , vol. 38, no. 11, pp. 2450–2525, 2020.
  * [2] Y. Zhang, J. Zhang, M. Di Renzo, H. Xiao, and B. Ai, “Performance analysis of ris-aided systems with practical phase shift and amplitude response,” _IEEE Transactions on Vehicular Technology_ , 2021.
  * [3] H. Zhang, B. Di, K. Bian, Z. Han, H. V. Poor, and L. Song, “Toward ubiquitous sensing and localization with reconfigurable intelligent surfaces,” _Proceedings of the IEEE_ , vol. 110, no. 9, pp. 1401–1422, 2022.
  * [4] S. Buzzi, E. Grossi, M. Lops, and L. Venturino, “Foundations of mimo radar detection aided by reconfigurable intelligent surfaces,” _IEEE Transactions on Signal Processing_ , vol. 70, pp. 1749–1763, 2022.
  * [5] H. Zhang, H. Zhang, B. Di, K. Bian, Z. Han, and L. Song, “Metalocalization: Reconfigurable intelligent surface aided multi-user wireless indoor localization,” _IEEE Transactions on Wireless Communications_ , vol. 20, no. 12, pp. 7743–7757, 2021.
  * [6] Z. Wang, Z. Liu, Y. Shen, A. Conti, and M. Z. Win, “Location awareness in beyond 5g networks via reconfigurable intelligent surfaces,” _IEEE Journal on Selected Areas in Communications_ , vol. 40, no. 7, pp. 2011–2025, 2022\.
  * [7] D. Dardari, N. Decarli, A. Guerra, and F. Guidi, “Los/nlos near-field localization with a large reconfigurable intelligent surface,” _IEEE Transactions on Wireless Communications_ , vol. 21, no. 6, pp. 4282–4294, 2022\.
  * [8] M. Luan, B. Wang, Y. Zhao, Z. Feng, and F. Hu, “Phase design and near-field target localization for ris-assisted regional localization system,” _IEEE Transactions on Vehicular Technology_ , vol. 71, no. 2, pp. 1766–1777, 2022.
  * [9] P. Wang, W. Mei, J. Fang, and R. Zhang, “Target-mounted intelligent reflecting surface for joint location and orientation estimation,” 2023. [Online]. Available: https://arxiv.org/abs/2301.09248
  * [10] R. Ghazalian, K. Keykhosravi, H. Chen, H. Wymeersch, and R. Jäntti, “Bi-static sensing for near-field ris localization,” in _GLOBECOM 2022 \- 2022 IEEE Global Communications Conference_ , 2022, pp. 6457–6462.
  * [11] K. Keykhosravi, M. F. Keskin, G. Seco-Granados, P. Popovski, and H. Wymeersch, “Ris-enabled siso localization under user mobility and spatial-wideband effects,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 16, no. 5, pp. 1125–1140, 2022.
  * [12] M. Skolnik, _Introduction to Radar Systems_. London: McGraw-Hill Education, 2002.
  * [13] S. Ellingson, “Path loss in reconfigurable intelligent surface-enabled channels,” in _2021 IEEE 32nd Annual International Symposium on Personal, Indoor and Mobile Radio Communications (PIMRC)_ , 2021, pp. 829–835.
  * [14] B. Wang, F. Gao, S. Jin, H. Lin, and G. Y. Li, “Spatial- and frequency-wideband effects in millimeter-wave massive mimo systems,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 13, pp. 3393–3406, 2018.
  * [15] A. Elzanaty, A. Guerra, F. Guidi, and M.-S. Alouini, “Reconfigurable intelligent surfaces for localization: Position and orientation error bounds,” _IEEE Transactions on Signal Processing_ , vol. 69, pp. 5386–5402, 2021.
  * [16] U. Spagnolini, _Statistical Signal Processing in Engineering_. John Wiley & Sons Ltd, 2018.


# Of Spiky SVDs and Music Recommendation

Darius Afchar<EMAIL_ADDRESS>,  Romain Hennequin  Deezer
ResearchParisFrance and  Vincent Guigue  AgroParisTech, MLIA – Sorbonne
UniversitésFrance

(2023)

###### Abstract.

The truncated singular value decomposition is a widely used methodology in
music recommendation for direct similar-item retrieval and downstream tasks
embedding musical items. This paper investigates a curious effect that we show
naturally occurring on many recommendation datasets: spiking formations in the
embedding space. We first propose a metric to quantify this spiking
organization’s strength, then mathematically prove its origin tied to
underlying communities of items of varying internal popularity. With this new-
found theoretical understanding, we finally open the topic with an industrial
use case of estimating how music embeddings’ top-k similar items will change
over time under the addition of data.

††copyright: acmcopyright††journalyear: 2023††doi: TBD††conference: 17th ACM
Conference on Recommender Systems; Sept 18–22, 2023; Singapore

##  1\. Introduction

There is no unique definition of music recommendation, but rather a range of
tasks that fall under this name: track sequence recommendation, context-aware
recommendation, playlist continuation or generation, similar track, artist, or
playlist retrieval (Schedl et al., 2018). These settings represent the many
use cases of recommendation found in the wild (e.g., in a streaming service).
Despite this proteiformity, many recommenders leverage a model of similar item
retrieval as a basis for their computation. A standard methodology is to rely
on some high-dimensional representations of the musical items to recommend.
That is why we see more and more work aiming to find general-purpose
representations of music tracks (e.g., self-supervision on music spectrograms
(Wang et al., 2022)) to be used in various downstream recommendation
scenarios. In the realm of collaborative filtering, the literature is
dominated by matrix factorization techniques to create such representations.
Compared to other research fields of machine learning, it has been reported in
recent years that modern approaches based on neural networks did not bring
significant benefit over simpler baselines of factorization and also suffered
many reproducibility issues (Rendle et al., 2019; Ferrari Dacrema et al.,
2019; Rendle et al., 2020; Ferrari Dacrema et al., 2021). In this context, the
truncated Singular Value Decomposition (SVD) is still a strong contender for
factorization. It has demonstrated consistent performances across a variety of
datasets (Cremonesi et al., 2010). SVD also benefits from rich theoretical
literature and has been implemented in efficient ways suitable to an
industrial scale (e.g., distributed RSVD (Halko et al., 2011; Constantine and
Gleich, 2011; Briand et al., 2021)). That is why SVD is still widely used in
music recommendation.

Nevertheless, despite not requiring much hyperparameter tuning, some blind
spots still exist in using SVD as embeddings (Levy et al., 2015). For
instance, after approximating a dataset $M$ with the low-rank decomposition
$U\Sigma V^{*}$, should we use $U$, $U\Sigma$, or $U\Sigma^{p}$ as a basis of
similarity? Should the dot product or cosine similarity be used? Many
practices for item retrieval rely on experimental results but have yet to find
theoretical justifications. In this paper, we submit the reader with a curious
geometry that tends to emerge when applying SVD to multiple recommendation
datasets: spikes. We observe that embedding vectors tend to self-organize
along lines that pass through the origin (see Figure 1). This effect is
understudied and goes unnoticed when using the cosine distance for similarity
since the normalization squashes all points on a hypersphere (Levy et al.,
2015). However, we show that leveraging the norm of vectors may be more
insightful than it first appears and that they should not be discarded. In
detail, we prove that spikes represent communities, and embeddings’ norms stem
from their varying importance within that community (i.e., intra-popularity).
These two latter features are notably prevalent in music data (Cross, 2001;
Jacobson and Sandler, 2009) and substantiate the need for theoretical support
in this context. As an opening to our results, we show that the norm is
strongly indicative of the stability of embeddings, which has industrial
implications on how music track representations evolve through time. To our
knowledge, this spiking behavior has yet to be leveraged for recommendation
before.

The paper proceeds as follows: we first propose a unified metric to quantify
this spiking effect, then we mathematically prove the equivalence of spiking
geometry with the presence of communities and intra-popularity within each
community. We discuss and open our formalization with a practical case of
evaluating the stability of music embeddings.

All experiments may be reproduced via our code repository at
https://github.com/deezer/spiky_svd.

Figure 1. Naturally emerging ”spikes” formations across several independently
computed SVDs. Since the embedding vectors are high-dimensional, we can only
display some projections of the points on slices of dimensions (indicated
above).

Four subplots

##  2\. Geometry of the spikes

We present and compute SVDs on several datasets and show the emergence of
spike formations. We then propose establishing a measure of the ”spikiness” of
a distribution of embeddings, to quantify the strength of this effect.

###  2.1. Background on the SVD

We first provide the reader with a brief theoretical background. The Truncated
Singular Value Decomposition (SVD111A full-rank decomposition is, of course,
intractable to use since it would create embeddings in millions of dimensions.
So we abuse notation and drop the T in the acronym. As a reference, the
truncated method can also be found under the designation pureSVD in the
recommenders community.) is a factorization technique that aims to approximate
a matrix $M\in\mathbb{R}^{n\times m}$ as the product of three matrices
$U\in\mathbb{R}^{n\times f}$, $\Sigma\in\textrm{diag}(\mathbb{R}^{f}_{+})$,
and $V\in\mathbb{R}^{m\times f}$, such that $\|M-\hat{M}\|_{F}^{2}$ is
minimal, with $\hat{M}=U\Sigma V^{*}$ and $U$ and $V$ are orthonormal by
columns (i.e., $U^{*}U=V^{*}V=I_{f}$). $\hat{M}$ is a matrix of rank $f$, with
$f$ chosen such that $f\ll\min(n,m)$. $\hat{M}$ is thus said to be a low-rank
approximation of $M$. The ”truncated” term refers to the optimal solution for
$\hat{M}$ where only the $f$ largest singular values and eigenvectors of $M$
are retained for the decomposition (Eckart-Young-Mirsky theorem (Eckart and
Young, 1936)).

In the context of music collaborative filtering, several remarks can be made.
$M$ is often a large sparse matrix ($n$ and $m$ are frequently counted in
millions) and often either denotes track-user, artist-user, or track-playlist
interactions. $M$ is often positive, and sometimes additionally binary (e.g.,
with the formalism of implicit feedbacks (Rendle et al., 2009)). Since the
goal is often to leverage $\hat{M}$ for similar item retrieval rather than
rebuilding $M$, the term $V$ is often ignored, and only the similarities
between each row of $U$ or $U\Sigma$ are considered. The close-form solution
to the SVD is such that $U$ is the $f$ dominant eigenvectors of $MM^{*}$.
Therefore, it does not matter whether we compute the SVD on $M$ or $MM^{*}$ to
retrieve $U$. With an abuse of notations, we identify $M$ to its symmetrized
version $MM^{*}$ since it does not impact the resulting $U$.

###  2.2. Datasets

In order to demonstrate that spiking SVDs are not tied to a particular
setting, we use six recommendation datasets: Spotify Million Playlist
(Spotify, 2020) (Spotify), The Art of the Mix (McFee and Lanckriet, 2012)
(AoTM-2011), Cornell’s playlists (Chen et al., 2012), LastFM-360K (Celma,
2010). We also include the classic MovieLens-25M dataset (Harper and Konstan,
2015) to hint at an extension of our result to other settings, and a private
dataset from Deezer containing twelve SVDs computed over one year in an
industrial context. This latter dataset will be further exploited in section 4
to evaluate the stability of embedding representations.

We compute SVDs using the same setting for all six datasets. We use the
symmetrized version $MM^{*}$ of the matrix. To curb popularity biases,
denoise, and sparsify the matrix, we filter the top-$k$ highest non-null
interactions per item. As commonly done in information retrieval, we apply a
positive pointwise mutual information (PPMI) normalization to the resulting
matrix, which further helps with popularity biases and retrieval performances
(Levy et al., 2015) and has been shown to be equivalent to modern skip-gram
formulations with negative sampling (Levy and Goldberg, 2014), which is common
for music embeddings (Cheng et al., 2017; Caselles-Dupré et al., 2018; Briand
et al., 2021). We choose $f=128$ for the decomposition and take $E=U\Sigma$ as
the resulting embeddings222Any other choice of normalization as $U\Sigma^{p}$
with $p\in[0,1]$ simply scales the embedding space differently but does not
impact the presence of spikes.. The rows of $E=(e_{1},...e_{n})$ provide a
latent representation for each item.

The obtained SVDs are displayed in Figure 1. Since we cannot lay out all 128
dimensions in a single figure, we report curious readers to our repository for
a detailed visualization. As a general comment, it seems that the spikiness of
embeddings is stronger on dominant dimensions (i.e., eigenvectors associated
with the higher singular values and hence the first indexes) and tends to
degrade to a blurry cloud of points in the higher indexes. This seems to align
with the literature on robust SVD (Stewart, 1998) that dominant components
first reconstruct underlying structures and then reconstruct noisy
perturbations. The Deezer and Spotify datasets exhibit the cleanest spikes
overall.

###  2.3. Spikiness metric

We could not find an off-the-shelf definition of spikes of points that worked
in high-dimension333The closest applicable notion is probably the measure of
kurtosis. We have tried several experiments with this statistic. However, it
turned out to be very sensitive to hyper-parameter choice and does not have a
straightforward extension to multi-dimensional variables.. In our context, it
seems suited to study the embeddings with the highest norms as many candidate
peaks for spikes. We can then measure whether many other points with a lower
norm are collinear to them444Note that the spikes all seem to pass through the
origin.. If a small number of peaks are collinear to most of the rest of the
distribution, it means that the distribution is spiky overall.

Since the distributions are noisy, we must instead rely on approximate
collinearity. We set some thresholds: $\theta$ and $\rho$. We then denote by
$e^{*}$ a given embedding with a high norm and say that any other vector $e$
belongs to the spike if $\cos(e^{*},e)>\cos(\theta)$. We iterate for candidate
peaks $e^{*}$ in descending order of norm until we have captured a ratio
$\rho$ of points in the distribution. This can be seen as a greedy heuristic
to the partial set cover problem555… which is notoriously NP-complete in high
dimension and thus could not be solved exactly with our number of considered
points.. The obtained number of spikes is divided by $n$ – the total number of
embeddings – and denoted Spk. We have $0<\texttt{Spk}\leq\rho$, where the
lower bound occurs when all vectors form a single spike, and the upper bound
means that no two points are collinear.

The results are given in Table 1 for $\cos\theta=0.9$ and $\rho=50\%$. Though
this measure is computed with a heuristic and constitutes an upper bound to
possibly more optimal choices of peaks to cover the distribution, the obtained
results seem consistent with our intuitions on the spikiness of each SVD. The
Deezer and Spotify come up as the most spiky. To fully spell out their
obtained result, half of the millions of embedding vectors’ directions – that
could be arbitrarily anything in 128 dimensions – fall under one of 361 and
1395 spikes. We also confirm that the embeddings are spikier in their first
dimensions and then exhibit a sudden increase of Spk when the noise dominates
in later dimensions.

Table 1. Estimated Spk on several datasets: needed ratio of points that are peaks of spikes to capture $\rho=50\%$ of the total distribution of $n$ embeddings, with $cos\theta=0.9$. The parameter @$f$ denotes several choices of truncation for the SVD. We add a Gaussian sampling as a reference on the improbability of finding collinear vectors in $\mathbb{R}^{f}$ – even with a $\theta$ angle tolerance. Dataset | Deezer | Spotify | AoTM-2011 | Cornell | LastFM-360K | Movielens-25M | $\mathcal{N}(0,1)$
---|---|---|---|---|---|---|---
Spk@128 | 0.04% | 0.14% | 18.1% | 4.55% | 2.42% | 7.09% | 50% $=\rho$
Spk@64 | 150/$n$ ¡ 0.01% | 228/$n$ ¡ 0.03% | 11.4% | 1.59% | 0.49% | 3.05% | 50%
Spk@32 | 50/$n$ ¡ 0.01% | 55/$n$ ¡ 0.01% | 2.31% | 0.49% | 0.14% | 1.12% | 50%
Spk@16 | 23/$n$ ¡ 0.01% | 23/$n$ ¡ 0.01% | 0.93% | 0.09% | 0.06% | 0.02% | 49.4%

We have shown a systematic spiking effect happening over multiple computed
SVDs. The result that a significantly low number of spikes may cluster half of
the music tracks is highly reminiscent of spectral clustering (but with
spikes) and hints at the presence of structured groups in the data.

##  3\. Modeling

We next demonstrate a formal equivalence between the spiking structure in the
SVD-embedding space and the presence of communities of varying node degrees in
the graph associated with $\hat{M}$. This enables us to get network properties
from the embeddings without explicitly building a graph, for which we show a
practical impact for recommendations.

###  3.1. Preliminaries

We define some notations and outline some basic properties to be used in the
next section.

####  3.1.1. Spikes in $M$

As argued in section 2.1, we may assume $M$ to be symmetric since we only seek
to retrieve similarities from the ”left side” of $M$. We denote by
$\textbb{S}^{n}(\mathbb{R}_{+})$ the space of symmetric matrices of
$\mathbb{R}^{n\times n}_{+}$. Applying the spectral theorem, $M$ is
diagonalizable as $PDP^{*}$, with $P$ an orthonormal basis of eigenvectors. We
rank the eigenvalues $(d_{i})$ and corresponding vectors in decreasing order
of absolute value. We then have two useful lemmas:

###### Lemma 3.1.

Matrix $M$ has a truncated SVD given by $\hat{M}=U\Sigma V^{*}$ such that:

$U=P_{1..f}\quad\quad\Sigma=|D_{1..f,1..f}|=\textrm{diag}(\sigma_{1},...\sigma_{f})\quad\quad
V=P_{1..f}\textrm{sign}(D_{1..f,1..f})$

where $|.|$ denote the term-by-term absolute value and $\textrm{sign}()$ the
sign function, since the SVD requires $\Sigma$ to be positive.

###### Proof.

We identify a desired decomposition. Note that the unitary eigenvectors of $P$
are not uniquely defined, and thus, so does the decomposition. ∎

###### Definition 3.2 (Spike formation).

A set of vectors $(e_{1},...e_{k})$ is said to form a spike – as studied in
2.2 – if we may rewrite them as $(\alpha_{1}\vec{s},...\alpha_{k}\vec{s})$
with $\vec{s}$ a unitary vector and
$(\alpha_{1},...\alpha_{k})\in\mathbb{R}_{+}^{k}$ their distance from the
origin.

###### Lemma 3.3.

Let $M$ have a decomposition $\hat{M}=U\Sigma V^{*}$. If the rows
$(u_{i_{1}},...u_{i_{k}})$ of $U$ belong to a spike, then, we have that the
rows $(i_{1},...i_{k})$ of $U\Sigma^{p}$ and of $V$ also form a spike:

$u_{i}\Sigma^{p}=\alpha_{i}\left(\vec{s}\cdot\Sigma^{p}\right)=\alpha_{i}\vec{s}^{\prime}\quad\quad\quad
v_{i}=u_{i}\cdot\textrm{sign}(D_{..f,..f})=\alpha_{i}\left(\vec{s}\cdot\textrm{sign}(D_{..f,..f})\right)=\alpha_{i}\vec{s}^{\prime\prime}$

The spiking effect is insensitive to any choice of embedding as $U\Sigma^{p}$
or $V$: only the directions of the spikes change. In particular, the vector
$\vec{s}^{\prime\prime}$ of $V$ remains unitary. In the rest of this section,
we arbitrarily fix the embeddings as $U\Sigma$.

####  3.1.2. Community graph models

The Stochastic Block Model (SBM) is an important model introduced in (Holland
et al., 1983) as a means to generate random graphs according to an underlying
structure in communities. Given a set of $n$ nodes indexed as $\llbracket
1,n\rrbracket=\\{1,2,...n\\}$, a partition $C:\llbracket
1,n\rrbracket\rightarrow\llbracket 1,K\rrbracket$ of the nodes into $K$
communities, and a matrix $B\in\mathcal{S}^{K}([0,1])$ of edge probabilities,
the SBM model samples an adjacency matrix $A\in\mathcal{S}^{n}(\\{0,1\\})$
such that $A_{i,j}=A_{j,i}\sim\mathcal{B}(B_{C(i),C(j)})$, with $\mathcal{B}$
the Bernouilli law. On average, we thus obtain a matrix with a constant edge
degree within and between any pair of communities of nodes. If there is not a
unique way of defining communities (Fortunato, 2010; Lehmann, 2014), it is
usually agreed that the edge density within each community must be superior to
that with other ones (Radicchi et al., 2004), which translates as
$B_{i,i}>B_{i,j}$ for all $i\neq j$. This last property is known as
assortivity. Music data has been shown to be particularly assortive (Jacobson
and Sandler, 2009).

The Degree-Corrected SBM (DCBM) (Karrer and Newman, 2011) introduces an
additional set of parameters $(\alpha_{i})_{1}^{n}\in[0,1]^{n}$ allowing the
edge probability to vary inside communities. Concisely, the DCBM proposes to
sample edges between $i$ and $j$ with a probability
$\alpha_{i}\alpha_{j}B_{C(i),C(j)}$. An example is depicted in Figure 2. In
order to make the parameters identifiable, we may renormalize $B$ to such that
$B_{i,i}=1$. For further properties on the DCBM, we refer readers to (Karrer
and Newman, 2011; Lei and Rinaldo, 2015).

###  3.2. Equivalence of spikes and degree-varying communities

We may finally spell out our main theorem linking recommendation embeddings
and modeling of communities:

###### Theorem 3.4.

Let $M\in\mathcal{S}^{n}(\mathbb{R}_{+})$ have a SVD $\hat{M}=U\Sigma V^{*}$
as defined in section 2.1. Let $E=U\Sigma=(e_{1},...e_{n})$ be embeddings that
exhibit $K$ distinct spikes within its rows whose directions are given by
$\vec{s}_{1},...\vec{s}_{K}$. If row $e_{i}$ belongs to a spike, let us define
$C(i)$ its corresponding spike index. We additionally denote $n^{\prime}\leq
n$ the total number of points belonging to spikes. For ease of notation, we
rearrange $E$ and $V$ such that the $n^{\prime}$ first rows belong to spikes.
Then, the submatrix $\tilde{M}=|\hat{M}_{1..n^{\prime},1..n^{\prime}}|$ is
proportional to the mean value of a DCBM.

###### Proof.

For all $i\in\llbracket
1,n^{\prime}\rrbracket,e_{i}=\alpha_{i}\vec{s}_{C(i)}$. Similarly,
$v_{i}=\alpha_{i}\vec{s}^{\prime}_{C(i)}$, thanks to lemma (3.3). Then, for
all $i,j\in\llbracket 1,n^{\prime}\rrbracket^{2}$,

$(\hat{M})_{i,j}=(EV^{*})_{i,j}=\alpha_{i}\alpha_{j}\langle\vec{s}_{C(i)},\vec{s}^{\prime}_{C(j)}\rangle\quad\quad\quad|\hat{M}_{1..n^{\prime},1..n^{\prime}}|=K^{2}\cdot\textrm{DCBM}\left((C(i))_{i=1}^{n^{\prime}}\,;\,B\,;\,\left(\frac{\alpha_{i}}{K}\right)_{i=1}^{n^{\prime}}\right)$

with $K=\max_{i}\alpha_{i}$, $B\in[0,1]^{K\times K}$ the matrix
$B_{k,l}=|\langle\vec{s}_{k},\vec{s}^{\prime}_{l}\rangle|$, and
$\frac{\alpha_{i}}{K}\in[0,1]$ as required. ∎

In practice, there are good reasons to believe that $\hat{M}$ will be
(roughly) non-negative and that the absolute value can be dropped in
$\tilde{M}$. When $f=1$, the Perron-Frobenius theorem (Perron, 1907) ensures
the non-negativity, which is also trivially true when $f\geq\textrm{rank}(M)$.
In between, we may expect negative values to principally arise from punctual
approximation artifacts of null values of $M$, rather than large negative
blocks having $\langle\vec{s}_{k},\vec{s}^{\prime}_{l}\rangle<0$.

Figure 2. Illustration of theorem (3.4) with a synthetic matrix with $n=1000$
and $n^{\prime}=3\times 200$ points forming three spikes. We highlight the
three direction vectors $\vec{s}_{\\{1,2,3\\}}$ in $E$ and their corresponding
communities $C_{\\{1,2,3\\}}$ in $\hat{M}\,$ (reindexed for clarity).

Two subplots

#### Interpretation

We have shown that the sub-matrix $\hat{M}_{1..n^{\prime},1..n^{\prime}}$
follows a structure in communities given by a DCBM, as illustrated in Figure 2
with a synthetic example. If we think geometrically, the DCBM is merely a
matrix with constant blocks that is multiplied element-wise by the rank-1
matrix
$[\alpha_{1},...\alpha_{n^{\prime}}]^{*}[\alpha_{1},...\alpha_{n^{\prime}}]$.
Our interpretation of this theorem is the following. Since $M$ is often sparse
(c.f., section 2.2), rebuilding as many communities as possible with rank-1
matrices is intuitively more optimal to minimize $\|M-\hat{M}\|_{F}^{2}$ than
”spending” many eigenvectors on a single community to fine-tune it with a
matrix of higher rank. The widespread emergence of spike formations we
uncovered across multiple datasets would thus stem from this general
underlying prioritization of rank-1 approximation of communities. This
coincidentally results in a DCBM with convenient theoretical properties that
we can now use. We leave the detection and analysis of eventual rank ¿ 1
approximation out of the scope of this paper.

#### Related work

We found theorem (3.4) to be related to a result from community detection that
the diagonalization of DCBM matrices leads to the appearance of spikes (Jin,
2012; Lei and Rinaldo, 2015). Let us stress the difference. We prove the
reciprocate (i.e., spikes $\Rightarrow$ DCBM), which is novel and in a more
general case since a portion of $M$ is noisy ($n^{\prime}<n$). In spirit also,
the former literature aim to correct this effect – seen as noise – while we
instead try to exploit it to gain knowledge on $\hat{M}$.

###  3.3. Practical advice

The duality of spike formations with DCBM communities has many convenient
properties for recommendation. We highlight some corollaries of our theorem
that provide practical advice when using embeddings computed with a SVD.

####  3.3.1. Fast degree estimation

As a first remark, the $(\alpha)_{i}$ parameters of our theorem’s DCBM may be
easily retrieved as the norm $\|e_{i}\|_{2}$ of embeddings that belong to a
spike. Access to the degree of a node of the denoised matrix $\hat{M}$ without
explicitly building a graph may open doors for many benefits. For instance,
the topology of a graph may be estimated from its degree distribution (e.g.,
whether the graph is scale-free), which has been suggested to be indicative of
a rather social or information nature of the music graphs (Cano et al., 2006).
Degree centrality has also been related to popularity biases in music data
(Celma and Cano, 2008; South et al., 2020) or used to create hierarchies
(Afchar et al., 2022a). More generally, the degree of a node is a strong
predictor of its robustness to change in a graph (Gasteiger et al., 2019;
Geisler et al., 2020). In section 4, we will further exploit such properties
of stability.

####  3.3.2. Cosine similarity or dot product?

The literature on community detection often seeks to curb the effect of the
degree heterogeneity of the DCBM (e.g., (Karrer and Newman, 2011; Qin and
Rohe, 2013; Gao et al., 2018)). One popular method is to normalize the
eigenvectors with a spherical clustering (Lei and Rinaldo, 2015). Normalizing
on a hypersphere turns out to correspond to the widespread use of cosine
similarity in recommendation settings. It effectively removes the
$(\alpha)_{i}$ terms, turning the DCBM into a simple SBM and thus attributing
a unique value to all the items of a given community (the cosine similarity
between their community direction $\vec{s}_{i}$ and that of the reference
embedding to retrieve items from). From the recommendation lens, both
unpopular and popular items will be reweighted with an equal relevance value.
This may be useful to improve discovery objectives – e.g., improving coverage
and serendipity metrics (Ge et al., 2010), but conversely promotes filter
bubbles by design (Pariser, 2011; Zhang et al., 2012).

On the flip side of the coin, measuring the internal popularity of a node
within a community – i.e., its ”importance” (Gao et al., 2018) – may be
helpful to navigate between communities while displaying their most
representative items. This setting is, in particular, more suited for artist
similarity and cold-start applications (Salha et al., 2019; Salha-Galvan et
al., 2021). There, a dot product will always return the item with the highest
norm of the community as top-similar, but may also merge popular items from
nearby communities in the retrieved items. The dot product thus seems more
suited to propose a panopticon view of the many available items, but with the
adverse risk of being biased by popularity and thus making long-tail items
harder to access.

As always, this choice of hyper-parameter depends on the task at hand and
desired criteria to optimize.

####  3.3.3. Tuning of $f$

The Spk metric we proposed in 2.3 may constitute a promising offline metric on
optimally choosing the truncation parameter $f$ in the SVD, i.e., finding when
recommendation performances are saturated despite adding more embedding
dimensions. In the light of theorem (3.4), we can now reinterpret our metric
as an approximation of $K$, where we have implicitly assumed that
$\hat{n}^{\prime}=\rho n$ points should be considered as being modeled by a
DCBM. We have additionally used the empirical heuristic that noise was
prevalent in the lower norm values of $E$. When the chosen $\hat{n}^{\prime}$
comes close to the ground-truth $n^{\prime}$, and the hypothesis starts to
fail, noise starts to permeate, and the number of spikes grows linearly with
the number of points (similar to the baseline Normal distribution where
collinear points are rare).

In Table 1, this effect was striking between Deezer and Spotify, where the Spk
metric was proportional between the two datasets for low values of $f$, then
exploded at $f=128$ for the latter. We conclude that $f=64$ for Spotify and
$f=128$ for Deezer constitutes an appropriate choice to capture clean
communities in the embeddings. This metric makes the assumption to treat non-
spikes as uninformative noise, which dismisses potential higher-ranking
approximations that may hold further information on the data. However, we
leave this last aspect for future work.

##  4\. Opening: industrial use-case

We open our results with an industrial use case to illustrate the few
previously highlighted properties. Specifically, designing an industrial
recommender system involves many challenges (e.g., scalability constraints)
that can render adopting a complex machine learning model cumbersome (Afchar
et al., 2022b). Insights like the one we propose can thus be valuable to gain
knowledge to better tune and evaluate recommenders leveraging simple item
similarity models, such as the SVD, without changing the computation pipeline.
We had access to a history of size twelve of music track embeddings from
Deezer, spanning between May 2022 and May 2023, computed on user data with a
truncated SVD as described in 2.2. These SVDs contain 1.5 million common music
tracks and have $f=128$.

Figure 3. Stability through time of the top-$500$ of Deezer embeddings
starting in May 2022 and compared with eleven successive embeddings within
five-week steps until May 2023. We condition the recommendation on a partition
of the embedding norms (indicated on the right). Both the cosine similarity
and dot product are considered. 95% confidence intervals are displayed.

Two plots showing decreasing similarity curves.

Music data is particularly prone to exhibiting community effects – linked to
its many historical and cultural groundings (Cross, 2001; Jacobson and
Sandler, 2009), but it is also ever-changing. New trends emerge every day
(e.g., lo-fi music, phonk), others vanish (e.g., emo, jumpstyle), or are
sometimes revivified (e.g., fans of Kate Bush thanks to Stranger Things’ OST),
or shape one of many types of evolution (Bródka et al., 2013). It is thus
natural to wonder how relevant some embeddings computed at a given date will
remain in the following months under the change of the underlying data.

In section 3.3, we stressed that the estimated degree is considered a good
predictor of the robustness of a node in a graph. We thus propose to link both
notions and show that embeddings of higher norms have better stability in time
than their lower-norm counterpart. Using the oldest SVD of May 2022 as a
reference, we partition the embeddings equally into five according to their
norm. We sample $1000$ points in each partition and compute their top-$500$
similar items in each SVDs – using either the cosine similarity or the dot
product – among the almost 300.000 track candidates in each partition. We
finally compute the mean Jaccard index (IoU) (Greene et al., 2010) to quantify
the relative similarity between the top items of the reference and subsequent
compared date.

The results are given in Figure 3. We can see that conditioning our results on
several norm ranges reveals strikingly different dynamics that could not have
been inferred solely from a mean analysis. The dot product notably displays
the most disparate curves where low-norm embeddings quickly become irrelevant
while higher-norm ones remain over 80% similar after almost a year. Note that
our point is not to state that change is necessarily undesirable in the
embeddings, though we could expect well-established music genres to remain
stable (e.g., jazz music tracks). Instead, we envision this type of
conditioned analysis as a means to better evaluate music recommenders by
discriminating between various dynamics of embeddings – and thus some expected
associated properties.

##  5\. Conclusion

Embeddings are commonly used in music recommendation. Yet, there are mostly so
in a black-box manner, and many practices rely on empirical knowledge. In this
paper, we have proven that the widespread emergence of spikes in the
embeddings of music items was tied to an underlying graph structure of degree-
varying communities. Our insights are particularly relevant for music data,
where such structure is found naturally due to cultural and historical
groundings. Reframing a community detection model through the recommendation
lens enabled us to formulate novel practical insights for music embeddings
(e.g., cosine or dot product, stability) that are better grounded in theory.
We hope to see more work bridging the gap between graph theory and
recommendation. Our future work includes the analysis of rank ¿ 1
approximation effects on embeddings.

## References

  * (1)
  * Afchar et al. (2022a) Darius Afchar, Romain Hennequin, and Vincent Guigue. 2022a.  Learning Unsupervised Hierarchies of Audio Concepts. In _23rd International Society for Music Information Retrieval Conference (ISMIR)_.
  * Afchar et al. (2022b) Darius Afchar, Alessandro Melchiorre, Markus Schedl, Romain Hennequin, Elena Epure, and Manuel Moussallam. 2022b.  Explainability in music recommender systems.  _AI Magazine_ 43, 2 (2022), 190–208.
  * Briand et al. (2021) Léa Briand, Guillaume Salha-Galvan, Walid Bendada, Mathieu Morlon, and Viet-Anh Tran. 2021\.  A semi-personalized system for user cold start recommendation on music streaming apps. In _27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining_. 2601–2609.
  * Bródka et al. (2013) Piotr Bródka, Stanisław Saganowski, and Przemysław Kazienko. 2013.  GED: the method for group evolution discovery in social networks.  _Social Network Analysis and Mining_ 3 (2013), 1–14.
  * Cano et al. (2006) Pedro Cano, Òscar Celma, Markus Koppenberger, and Javier M. Buldú. 2006.  Topology of music recommendation networks.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 16, 1 (01 2006).  https://doi.org/10.1063/1.2137622
  * Caselles-Dupré et al. (2018) Hugo Caselles-Dupré, Florian Lesaint, and Jimena Royo-Letelier. 2018.  Word2vec applied to recommendation: Hyperparameters matter. In _12th ACM Conference on Recommender Systems_ _(RecSys)_. 352–356.
  * Celma (2010) Òscar Celma. 2010\.  _Music Recommendation and Discovery in the Long Tail_.  Springer.  http://ocelma.net/MusicRecommendationDataset/lastfm-360K.html
  * Celma and Cano (2008) Òscar Celma and Pedro Cano. 2008.  From hits to niches? or how popular artists can bias music recommendation and discovery. In _Proceedings of the 2nd KDD Workshop on Large-Scale Recommender Systems and the Netflix Prize Competition_. 1–8.
  * Chen et al. (2012) Shuo Chen, Josh L Moore, Douglas Turnbull, and Thorsten Joachims. 2012.  Playlist prediction via metric embedding. In _18th ACM SIGKDD international conference on Knowledge discovery and data mining (KDD)_. 714–722.  https://www.cs.cornell.edu/~shuochen/lme/data_page.html
  * Cheng et al. (2017) Zhiyong Cheng, Jialie Shen, Lei Zhu, Mohan S Kankanhalli, and Liqiang Nie. 2017.  Exploiting Music Play Sequence for Music Recommendation. In _IJCAI_ , Vol. 17. 3654–3660.
  * Constantine and Gleich (2011) Paul G Constantine and David F Gleich. 2011.  Tall and skinny QR factorizations in MapReduce architectures. In _2nd international workshop on MapReduce and its applications_. 43–50.
  * Cremonesi et al. (2010) Paolo Cremonesi, Yehuda Koren, and Roberto Turrin. 2010\.  Performance of Recommender Algorithms on Top-n Recommendation Tasks. In _4th ACM Conference on Recommender Systems (RecSys)_. 39–46.  https://doi.org/10.1145/1864708.1864721
  * Cross (2001) Ian Cross. 2001\.  Music, cognition, culture, and evolution.  _Annals of the New York Academy of sciences_ 930, 1 (2001), 28–42.
  * Eckart and Young (1936) Carl Eckart and Gale Young. 1936.  The approximation of one matrix by another of lower rank.  _Psychometrika_ 1, 3 (01 Sep 1936), 211–218.  https://doi.org/10.1007/BF02288367
  * Ferrari Dacrema et al. (2021) Maurizio Ferrari Dacrema, Simone Boglio, Paolo Cremonesi, and Dietmar Jannach. 2021.  A troubling analysis of reproducibility and progress in recommender systems research.  _ACM Transactions on Information Systems (TOIS)_ 39, 2 (2021), 1–49.
  * Ferrari Dacrema et al. (2019) Maurizio Ferrari Dacrema, Paolo Cremonesi, and Dietmar Jannach. 2019\.  Are we really making much progress? A worrying analysis of recent neural recommendation approaches. In _13th ACM conference on recommender systems (RecSys)_. 101–109.
  * Fortunato (2010) Santo Fortunato. 2010\.  Community detection in graphs.  _Physics reports_ 486, 3-5 (2010), 75–174.
  * Gao et al. (2018) Chao Gao, Zongming Ma, Anderson Y. Zhang, and Harrison H. Zhou. 2018.  Community Detection in Degree-Corrected Block Models.  _The Annals of Statistics_ 46, 5 (2018), 2153–2185.  https://www.jstor.org/stable/26542860
  * Gasteiger et al. (2019) Johannes Gasteiger, Stefan Weißenberger, and Stephan Günnemann. 2019\.  Diffusion improves graph learning.  _Advances in neural information processing systems (NeurIPS)_ 32 (2019).
  * Ge et al. (2010) Mouzhi Ge, Carla Delgado-Battenfeld, and Dietmar Jannach. 2010\.  Beyond accuracy: evaluating recommender systems by coverage and serendipity. In _4th ACM conference on Recommender systems (RecSys)_. 257–260.
  * Geisler et al. (2020) Simon Geisler, Daniel Zügner, and Stephan Günnemann. 2020\.  Reliable graph neural networks via robust aggregation.  _Advances in Neural Information Processing Systems (NeurIPS)_ 33 (2020), 13272–13284.
  * Greene et al. (2010) Derek Greene, Donal Doyle, and Padraig Cunningham. 2010\.  Tracking the evolution of communities in dynamic social networks. In _2010 international conference on advances in social networks analysis and mining_. IEEE, 176–183.
  * Halko et al. (2011) Nathan Halko, Per-Gunnar Martinsson, and Joel A Tropp. 2011\.  Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions.  _SIAM review_ 53, 2 (2011), 217–288.
  * Harper and Konstan (2015) F Maxwell Harper and Joseph A Konstan. 2015.  The movielens datasets: History and context.  _Acm transactions on interactive intelligent systems (tiis)_ 5, 4 (2015), 1–19.  https://grouplens.org/datasets/movielens/25m/
  * Holland et al. (1983) Paul W Holland, Kathryn Blackmond Laskey, and Samuel Leinhardt. 1983.  Stochastic blockmodels: First steps.  _Social networks_ 5, 2 (1983), 109–137.
  * Jacobson and Sandler (2009) Kurt Jacobson and Mark Sandler. 2009.  Musically meaningful or just noise? an analysis of on-line artist networks. In _Computer Music Modeling and Retrieval. Genesis of Meaning in Sound and Music: 5th International Symposium_. Springer, 107–118.
  * Jin (2012) Jiashun Jin. 2012\.  Fast network community detection by SCORE.  _The Annals of Statistics_ 43 (11 2012).  https://doi.org/10.1214/14-AOS1265
  * Karrer and Newman (2011) Brian Karrer and Mark EJ Newman. 2011.  Stochastic blockmodels and community structure in networks.  _Physical review E_ 83, 1 (2011), 016107.
  * Lehmann (2014) Sune Lehmann. 2014\.  _Community Detection, Current and Future Research Trends_.  Springer New York, New York, NY, 214–220.  https://doi.org/10.1007/978-1-4614-6170-8_27
  * Lei and Rinaldo (2015) Jing Lei and Alessandro Rinaldo. 2015.  Consistency of spectral clustering in stochastic block models.  _The Annals of Statistics_ 43, 1 (2015), 215–237.
  * Levy and Goldberg (2014) Omer Levy and Yoav Goldberg. 2014.  Neural word embedding as implicit matrix factorization. In _Advances in neural information processing systems (NIPS)_ , Vol. 27.
  * Levy et al. (2015) Omer Levy, Yoav Goldberg, and Ido Dagan. 2015.  Improving distributional similarity with lessons learned from word embeddings.  _Transactions of the association for computational linguistics (ACL)_ 3 (2015), 211–225.
  * McFee and Lanckriet (2012) B. McFee and G. R. G. Lanckriet. 2012.  Hypergraph models of playlist dialects. In _13th International Symposium for Music Information Retrieval (ISMIR)_.  https://brianmcfee.net/data/aotm2011.html
  * Pariser (2011) Eli Pariser. 2011\.  _The filter bubble: What the Internet is hiding from you_.  penguin UK.
  * Perron (1907) Oskar Perron. 1907\.  Zur theorie der matrices.  _Math. Ann._ 64, 2 (1907), 248–263.
  * Qin and Rohe (2013) Tai Qin and Karl Rohe. 2013\.  Regularized spectral clustering under the degree-corrected stochastic blockmodel.  _Advances in neural information processing systems (NIPS)_ 26 (2013).
  * Radicchi et al. (2004) Filippo Radicchi, Claudio Castellano, Federico Cecconi, Vittorio Loreto, and Domenico Parisi. 2004\.  Defining and identifying communities in networks.  _Proceedings of the national academy of sciences_ 101, 9 (2004), 2658–2663.
  * Rendle et al. (2009) Steffen Rendle, Christoph Freudenthaler, Zeno Gantner, and Lars Schmidt-Thieme. 2009.  BPR: Bayesian Personalized Ranking from Implicit Feedback. In _25th Conference on Uncertainty in Artificial Intelligence_ _(UAI ’09)_. 452–461.
  * Rendle et al. (2020) Steffen Rendle, Walid Krichene, Li Zhang, and John Anderson. 2020.  Neural collaborative filtering vs. matrix factorization revisited. In _14th ACM Conference on Recommender Systems (RecSys)_. 240–248.
  * Rendle et al. (2019) Steffen Rendle, Li Zhang, and Yehuda Koren. 2019.  On the Difficulty of Evaluating Baselines: A Study on Recommender Systems.  _CoRR_ abs/1905.01395 (2019).  arXiv:1905.01395
  * Salha et al. (2019) Guillaume Salha, Stratis Limnios, Romain Hennequin, Viet-Anh Tran, and Michalis Vazirgiannis. 2019\.  Gravity-inspired graph autoencoders for directed link prediction. In _28th ACM international conference on information and knowledge management (CIKM)_. 589–598.
  * Salha-Galvan et al. (2021) Guillaume Salha-Galvan, Romain Hennequin, Benjamin Chapus, Viet-Anh Tran, and Michalis Vazirgiannis. 2021.  Cold start similar artists ranking with gravity-inspired graph autoencoders. In _15th ACM Conference on Recommender Systems (RecSys)_. 443–452.
  * Schedl et al. (2018) Markus Schedl, Hamed Zamani, Ching-Wei Chen, Yashar Deldjoo, and Mehdi Elahi. 2018.  Current challenges and visions in music recommender systems research.  _International Journal of Multimedia Information Retrieval_ 7 (2018), 95–116.
  * South et al. (2020) Tobin South, Matthew Roughan, and Lewis Mitchell. 2020\.  Popularity and centrality in Spotify networks: critical transitions in eigenvector centrality.  _Journal of Complex Networks_ 8, 6 (2020), cnaa050.
  * Spotify (2020) Spotify. 2020.  _Spotify Million Playlist_.  https://research.atspotify.com/2020/09/the-million-playlist-dataset-remastered/
  * Stewart (1998) Gilbert W Stewart. 1998\.  _Perturbation theory for the singular value decomposition_.  Technical Report.
  * Wang et al. (2022) Luyu Wang, Pauline Luc, Yan Wu, Adria Recasens, Lucas Smaira, Andrew Brock, Andrew Jaegle, Jean-Baptiste Alayrac, Sander Dieleman, Joao Carreira, et al. 2022\.  Towards learning universal audio representations. In _IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. 4593–4597.
  * Zhang et al. (2012) Yuan Cao Zhang, Diarmuid Ó Séaghdha, Daniele Quercia, and Tamas Jambor. 2012\.  Auralist: introducing serendipity into music recommendation. In _5th ACM international conference on Web search and data mining (WSDM)_. 13–22.


MIDLMedical Imaging with Deep Learning 2019 MIDL 2019 – Extended Abstract
Track Éric Fournié<EMAIL_ADDRESS>
Matthias Baer-Beck<EMAIL_ADDRESS>
Karl Stierstorfer<EMAIL_ADDRESS>
Siemens Healthcare GmbH, Computed Tomography, Siemensstr. 3, 91301 Forchheim,
Germany

# CT Field of View Extension Using Combined Channels Extension and Deep
Learning Methods

###### Abstract

This paper proposes a method to extend the field of view of computed
tomography images. In a first step, the field of view is increased by
extrapolating linearly the outer channels in the sinogram space. The modified
sinogram is then used to reconstruct extended field of view (EFoV) images
containing artifacts due to the channels extension. In a second step, those
artifacts are reduced by a deep learning network in image space.
The proposed method has been evaluated on a collection of clinical scans. The
resulting volumes have been checked for consistency and plausibility and
compared to an existing state of the art EFoV method.

###### keywords:

Extended field of view, CT, sinogram, U-Net, artifact reduction.

##  1 Introduction

X-ray computed tomography (CT) reconstructs slices of the scanned structure
from the sinograms acquired by rotating a detector array and an X-ray emitting
tube around the scanner gantry. One of its fundamental limitation is that only
data that is located inside the scan field of view (SFoV) of the CT scanner
can be reconstructed without artifacts. The size of the SFoV is determined by
the geometrical properties of the CT scanner, i.e. by the size of the detector
and by the distance of the X-ray tube and the detector from the isocenter of
the CT scanner. Sinograms nonetheless contains information from regions
outside of the SFoV. Artificially extending their angular width, for example
by padding or extrapolating the channels as in fig. 1, will allow for a
partially reconstructed EFoV region polluted by artifacts.

Figure 1: Field of view extension by channel padding or extrapolation.

The proposed method first extends the sinograms in channel direction, then
uses a U-Net [Ronneberger et al.(2015)Ronneberger, P.Fischer, and Brox] to
clean up the artifacts. The result is a realistic and computationally cheap
extension of the field of view. Compared to existing methods [Hsieh et
al.(2004)Hsieh, Chao, Thibault, Grekowicz, and Horst, Kunze et al.(2007)Kunze,
Härer, and Stierstorfer, Bruder et al.(2008)Bruder, Suess, and Stierstorfer],
the proposed approach delivers improved results at a low computational cost.

##  2 Method

The proposed method consists of a two-step approach. Let us consider a scanned
object $y$ and its sinogram obtained by the Radon transform $\Re(y)$. In a
first step, the angular width of the sinogram is increased by linear
extrapolation $g(\Re(y))$ of the outermost channels linearly towards zero,
thus increasing the number of channels by a factor corresponding to
$EFoV/SFoV$. We can now use any method $\Re^{+}$ such as filtered
backprojection (FBP) to reconstruct an EFoV image $X=\Re^{+}(g(\Re(y)))$
containing artifacts. The SFoV region is only slightly influenced by the
extrapolation via the convolution step of the FBP and therefore can be used
for diagnostic purposes without limitation.

The artifacts produced outside of the SFoV region in $X$ are then cleaned up
by a U-Net $F$ trained to match $(X,y)$ pairs. The output of the network
$\hat{y}=F(X)$ can then be used as a realistic estimation of $y$ for the EFoV
region. With this approach, the problem faced is the reduction of artifacts
instead of inpainting large portions of the image, thus largely reducing the
space of acceptable solutions.

Figure 2: Data generation pipeline and network architecture

To remediate possible modifications of the SFoV region by the network, it is
replaced by the corresponding SFoV of $X$. In order to improve the final
SFoV/EFoV transition, the network output is adapted to the input distribution:

$\hat{y}=\frac{\hat{y}-\mu_{\hat{y}}}{\sigma_{\hat{y}}}\cdot\sigma_{X}+\mu_{X}$

To train the network, a collection of clinical CT scan slices $y$ are
projected to their Radon transform to obtain sinograms. The sinograms are
extended by extrapolation of the outer channels as described previously (cf.
fig. 1 and 2) and backprojected to images $X$ used as inputs of the network.
The loss function is only computed in the EFoV region and is defined as the
weighted sum of the structural dissimilarity $DSSIM=(1-SSIM)/2$ [Wang et
al.(2004)Wang, Bovik, Sheikh, and Simoncelli] and the mean squared error.

##  3 Results and Conclusions

To evaluate performance on clinical data, a plugin implementing the various
steps of the method including execution of the trained network has been
developed for ReconCT, a proprietary offline reconstruction software developed
by Siemens Healthineers. ReconCT can also reconstruct EFoV images using the
existing state of the art HDFoV method.

Scanner raw data corresponding to 11 patients and 4 phantoms has been
collected and EFoV images reconstructed using the proposed method and HDFoV.
While no ground truth exists for the patient data sets, reconstructions have
been checked for artifacts errors and plausibility in axial, coronal and
sagittal slices. Some qualitative results are presented in fig. 3. When
comparing the results depicted in figure 3 it shows up that the results with
the proposed method are superior to the HDFoV results for cases 2, 4 and 5.
For cases 1 and 3 the proposed method cannot outperform the HDFoV method as it
either introduces artificial structures (case 1) or fails to reconstruct some
parts of the patient anatomy (arm in case 3).

Figure 3: Comparison between our method (top) and HDFoV (bottom)

Extrapolation by linear extrapolation followed by the removal of artifacts in
the corresponding reconstruction using a deep learning network proved to be an
efficient method to extend the field of view of CT scans and shows potential
to improve quality of EFoV reconstructions. Like all methods aiming for the
reconstruction of volume areas using incomplete projection data, the EFoV
regions cannot be used for diagnostic but they are sufficient for an array of
clinical applications where the normal SFoV is too narrow, such as radio
therapy planning or the imaging of obese patients.

## References

  * [Bruder et al.(2008)Bruder, Suess, and Stierstorfer] H. Bruder, C. Suess, and K. Stierstorfer.  Efficient extended field of view (efov) reconstruction techniques for multi-slice helical ct.  6913:81–, 04 2008.  10.1117/12.768477.
  * [Hsieh et al.(2004)Hsieh, Chao, Thibault, Grekowicz, and Horst] J. Hsieh, E. Chao, J. Thibault, B. Grekowicz, and A. Horst.  A novel reconstruction algorithm to extend the CT scan field-of-view.  _Med. Phys._ , 31:2385 – 2391, 2004.
  * [Kunze et al.(2007)Kunze, Härer, and Stierstorfer] H. Kunze, W. Härer, and K. Stierstorfer.  Iterative extended field of view reconstruction.  _Proceedings of SPIE - The International Society for Optical Engineering_ , 03 2007.  10.1117/12.707336.
  * [Ronneberger et al.(2015)Ronneberger, P.Fischer, and Brox] O. Ronneberger, P.Fischer, and T. Brox.  U-net: Convolutional networks for biomedical image segmentation.  In _Medical Image Computing and Computer-Assisted Intervention (MICCAI)_ , volume 9351 of _LNCS_ , pages 234–241. Springer, 2015.  URL http://lmb.informatik.uni-freiburg.de/Publications/2015/RFB15a.  (available on arXiv:1505.04597 [cs.CV]).
  * [Wang et al.(2004)Wang, Bovik, Sheikh, and Simoncelli] Z. Wang, A. Bovik, H. Sheikh, and E. Simoncelli.  Image quality assessment: From error visibility to structural similarity.  _IEEE TRANSACTIONS ON IMAGE PROCESSING_ , 13(4):600–612, 2004.