Paper ID #889

Spatial Positioning Token (SPToken) for Smart Parking

Roman Overko${}^{\text{1*}}$, Rodrigo Ordóñez-Hurtado${}^{\text{2}}$, Sergiy
Zhuk${}^{\text{2}}$, Robert Shorten${}^{\text{3}}$

1\. University College Dublin, Dublin, Ireland. Email
<EMAIL_ADDRESS>
2\. IBM Research, Damastown Industrial Park, Dublin 15, Ireland

3\. Imperial College London, South Kensington, London SW7 2AZ, United Kingdom

Abstract
In this paper, we describe an approach to guide drivers searching for a
parking space (PS). The proposed system suggests a sequence of routes that
drivers should traverse in order to maximise the expected likelihood of
finding a PS and minimise the travel distance. This system is built on our
recent architecture SPToken, which combines both Distributed Ledger Technology
(DLT) and Reinforcement Learning (RL) to realise a system for the estimation
of an unknown distribution without disturbing the environment. For this, we
use a number of virtual tokens that are passed from vehicle to vehicle to
enable a massively parallelised RL system that estimates the best route for a
given origin-destination (OD) pair, using crowdsourced information from
participant vehicles. Additionally, a moving window with reward memory
mechanism is included to better cope with non-stationary environments.
Simulation results are given to illustrate the efficacy of our system.

Keywords:
Smart Parking, Reinforcement Learning

##  1 Introduction

A problem of immense importance to drivers is that of finding an available
parking space (PS), which is also a societal problem affecting air pollution,
congestion, and noise emission due to the large number of vehicles currently
on the streets. For illustration, a case study of a small business district in
Los Angeles [1] indicated that 730 tons of CO2 where produced and 47,000 gal
of gasoline were burned in a year by cars searching for parking. Another study
by McKinsey111https://www.mckinsey.com/industries/public-and-social-
sector/our-insights/the-smart-city-solution reported that the average car
owners in Paris spend about four years of their lives searching for PSs. The
emerging paradigm of intelligent transport systems (ITS) offers a great hope
to alleviate these problems by enabling a range of parking guidance systems
that are rapidly becoming an essential part of future sustainable cities. The
most common approach are information boards displaying the available PSs at
various locations around the city. This system provides a valuable guidance
for drivers to avoid areas with potentially limited PS availability (PSA), but
often leads to localised congestion around areas with the largest number of
available PSs which is caused by (i) all the drivers receiving identical data,
and (ii) the trend of most drivers to typically choose the areas with large
PSA. More modern systems [2] provide extensive parking information directly to
the drivers either to their smart devices or to the cars. The wealth of PSA
information is provided to the drivers through the use of PS sensors. Such
systems thus provide drivers with the exact location of currently free PSs
along with e.g. their respective prices. The main disadvantage of such systems
is that a large investment cost is required as sensors are expensive to deploy
and maintain, thus limiting large-scale deployments. On the other hand, if
moving vehicles are used to collect raw environment data (e.g. from dash cams
or LiDAR sensors) to be turned into PSA information in cloud/MEC servers, then
the cost and complexity of maintaining a central database of parking
information built from raw environment data is also prohibitive from the
perspective of the burden on the telecommunications network.

In the context of the EU, the European Commission has created a platform for
the deployment of Cooperative ITS, namely the C-ITS
Platform222https://ec.europa.eu/transport/themes/its/c-its_en, with the aim of
improving the efficiency of road transport at European-wide scale through a
number of identified (essential) applications including probe vehicle data and
parking services. In alignment with this, the smart parking solution proposed
in this paper aims to use moving vehicles as collectors and processors of raw
environment data so that PSA can be shared as floating car data. For instance,
if the ICT architecture proposed in [3] is integrated with the Spatial
Potisioning Token (SPToken) framework [4, 5], then SPToken observers could
include connected automated vehicles (CAVs), road-side units (RSUs), or MEC
servers, and the floating car data can be transmitted via (perhaps extended)
Cooperative Awareness Messages and Collective Perception Messages [6] (or a
new message type to be developed).

This paper is organised as follows. First, we provide a brief summary of
related work on applications in the area of smart parking. Second, we proceed
to describe in detail the integration of SPToken into a smart parking
solution. Third, we provide numerical validation of the proposed approach.
Finally, we close the paper with some concluding remarks.

##  2 Related work

We proceed to present a short overview of relevant smart parking solutions. In
[2], the authors propose a MEC-based architecture using VANETs to improve the
parking experience concerning average parking cost, fuel wastes, and vehicle
exhaust emissions, and the solution gives drivers the option to integrate
their own parking preferences. However, this solution relies on parking
information from interconnected MEC servers deployed at parking lots, thus
substantially restricting its scope to such locations and not allowing for the
integration of unsupervised on-street parking. Concerning privacy-preserving
data management, several architectures for smart parking solutions have
already been proposed. In [7], a blockchain-based descentralised parking
management service is designed, which ensures anonymous authentication and
resistance to data linkability, among others. In [8], the authors introduce a
privacy-preserving smart parking navigation system based on the integration of
cloud services and vehicular communication, that allows for identity/service
authentication and traceability. Even though [7] and [8] provide means for
parking data retrieval (including PSA and navigation guidance) or PS
reservation, their allocation process relies on parking-related databases
maintained with availability data provided by e.g. parking lot owners, and it
does not involve any optimisation problem. The literature on applications of
machine learning to smart parking is vast and we simply refer to a few
relevant publications. In [9], the authors propose a Q-learning-based solution
to find the nearest PS by minimising the total covered distance, time taken
and consumed energy. This approach is evaluated only on stationary
environments, which is restrictive in many real-world applications including
traffic management problems. Situations with non-stationary models may lead to
suboptimal policies provided by RL algorithms. In this context, only a few
existing works have considered to improve RL algorithms for non-stationary
environments. A recent example of such works includes [10], in which the
authors have developed a model-free RL method to effectively detect changes in
the rewards and transition dynamics, and validated the proposed approach on
randomly generated data and traffic signal control problems as well.

##  3 SPToken for Smart Parking

In this work, we make use of the Distributed Ledger Technology (DLT) based
SPToken architecture proposed in [4, 5] to design a solution for the on-street
parking problem. As detailed in [4, 5], the SPToken framework allows to
explore a given environment without perturbing it, this achieved through the
use of (i) virtual entities referred to as tokens in combination with (ii)
crowdsourcing approaches such as multi-agent RL, which makes SPToken suitable
for a wide range of smart mobility applications including smart parking. The
key idea is to use the tokens as virtual containers to be filled with floating
car data, and allow tokens to “jump” from vehicle to vehicle whenever is
required to complete routes determined by the underlying RL algorithm. A
participating vehicle can “collect” a free token via a DLT transaction if its
route coincides with the token route, after which the vehicle with the token
turns into an RL agent that updates the distributed ledger with some
information (e.g. spent travel time, measured roadside parking availability)
whenever it passes an observer. The vehicle must “deposit” the token via
another DLT transaction when it deviates from the token route. Since tokens
are transported by participating vehicles (not forced to follow the entire
token route) already present in an urban scenario, SPToken guarantees that the
environment stays undisturbed during the probing process. Note that the
physical presence of vehicles is imperative to collect/deposit tokens, and so
a Proof-of-Position mechanism is also included (see [5] for details). SPToken
also has all the advantages of DLT such as data privacy preservation, data
ownership retention, and misuse/spamming prevention. As a base requirement for
participating vehicles, we assume they can provide on-street PSA data along
their routes as a result of applying analytics to (raw) real-time environment
data collected by their on-board sensors (e.g. point clouds from LiDAR
sensors). We assume this requirement is reasonably satisfied based on reported
mechanisms for on-board detection of parking availability [11, 12].

Note that the proposed solution is not a straightforward application of the
SPToken architecture to smart parking, and a number of extensions (also to the
underlying RL algorithm [4]) were required: (i) a new mechanism called moving
window with reward memory (MWRM) for MUBEV [5] to improve the performance in
non-stationary environments; (ii) a more advanced action selection strategy to
increase exploration at early stages of the learning process; (iii) an
improved design of the reward function for the RL algorithm to account for a
reduced number of tuning parameters; and (iv) a design parameter $\alpha$ is
introduced as a means for users to provide their preferences in terms of the
their preferred optimisation objectives. Additionally, the application of
SPToken to the on-street parking problem required dealing with a non-trivial
multi-objective optimisation problem: our solution involves the investigation
of optimal routes connecting origin-destination (OD) pairs so that in such
routes (a) travel distance is minimised, and (b) PSA is maximised. While (a)
involves a (widely studied) shortest path (SP) problem, the solution for (b)
requires solving a longest path problem of NP-hard nature. Therefore, we
decided to integrate these two problems into a graph-based, multi-objective
optimisation approach by linearly combining two directed weighted graphs: a
travel Distance Graph (DG), and a Parking Availability Graph (PAG). Both
graphs come from the road network, where each vertex is a road junction, each
edge is a collection of roads, weights $\Lambda$ in DG represent the physical
(time invariant) length of given roads, and weights $K$ in PAG are in function
of the (time-variant) number of unoccupied on-street PSs along given roads at
certain time. As a result, the graph weights of the convex linear combination
$G=(\alpha)\Lambda+(1-\alpha)K$ are computed as

$\mathcal{W}_{G}(\alpha,s,\lambda_{s},\nu_{s}(t))=\left(\alpha\right)*\lambda_{s}+\left(1-\alpha\right)*\kappa_{s}\left(\nu_{s}(t)\right),$
(1)

where: $\alpha$ is a user-defined parameter that satisfies
$\alpha_{min}\leq\alpha\leq 1$ and represents the priority given to the
distance cost; $\lambda_{s}$ is the length of edge $s$ in G;
$\kappa_{s}(\nu_{s}(t))=-\nu_{s}(t)*L_{max}/C_{max}$, with $\nu_{s}(t)$ being
the current roadside PSA along the road links included in edge $s$ at time
$t$, $L_{max}$ the maximum edge length in the road network, and $C_{max}$ the
maximum edge parking capacity in the road network. While SP routing is widely
used as the de-facto default solution (DS) in most navigation systems, SP does
not always guarantees maximum PSA. However, the total PS capacity along road
links is indeed time-invariant, and thus in our approach we compute DSs using
Equation 1 with the value of $\alpha$ provided by the users. In addition,
$\alpha_{min}$ is determined empirically such that the graph $G$ does not
contain negative cycles for weights computed as
$\mathcal{W}_{G}(\alpha,s,\lambda_{s},c_{s})$, where $c_{s}$ is the total
(time-invariant) roadside parking capacity along the road segments included in
edge $s$. Consequently, if $\alpha\in[\alpha_{min},1]$, then we can use an SP
algorithm to compute DSs to be used as the initial policy of the RL algorithm.
Clearly, DSs degenerate into SPs when $\alpha=1$ (as per Equation 1).

###  3.1 Optimal policy search algorithm

As previously mentioned, we use an RL strategy to solve our target problem. A
full RL problem is usually modeled as a Markov Decision Process (MDP). Our
decision problem is a stationary Finite Horizon MDP (FHMDP), which is a
discrete-time stochastic control process defined by a tuple
$\langle\mathcal{S},\mathcal{A},P,\mathcal{R},H\rangle$, where: $\mathcal{S}$
is the state space, $\mathcal{A}$ is the action space, $P$ is the tensor of
transition probabilities, $\mathcal{R}$ is the reward matrix, and $H$ is the
length of the time horizon. Let $S=|\mathcal{S}|$, $A=|\mathcal{A}|$ and
$A_{s}=|\mathcal{A}_{s}|$, that is, $S$ is the number of states, $A$ is the
total number of actions, and $A_{s}$ is the number of actions allowable in
state $s\in\mathcal{S}$. We assume that all the states are fully observable.
We also adopt this general approach: each state $s$ has a different set of
allowable actions $\mathcal{A}_{s}$, hence
$\mathcal{A}=\bigcup_{s\in\mathcal{S}}\mathcal{A}_{s}$. Finally,
$P(s^{\prime}|s,a)$ is the probability of transition to state $s^{\prime}$ if
action $a$ is chosen in state $s$, and $\mathcal{R}(s,a)$ is the reward of
playing action $a$ in state $s$.

An RL agent, i.e. the decision maker, interacts with the environment in a
sequence of episodes, each episode with length $H$. Within an episode $k$ at
each time step $t\in[H]$ (where $[\bullet]=\\{1,2,...,\bullet\\}$ with
$\bullet$ an integer number), the agent plays action $a$ based on the
observation of state $s$ and receives a reward $r(s,a)$. The next state
$s^{\prime}$ is drawn from a distribution $P$ which defines the trajectory of
the MDP, that is $s^{\prime}\sim P(\cdot|s,a)$. Over each episode $k$, the
agent selects actions according to a policy $\pi_{k}$, which maps states and
time steps to actions. $\pi_{k}$ is updated after $H$ interactions with the
environment, i.e. at the end of each episode. Note that in stationary MDPs,
the transition probabilities and reward distribution do not vary with time
step $t$. The policy $\pi_{k}$, however, is generally time-step-dependent for
FHMDPs. The expected return until the end of an episode is represented by a
value function (for state $s$, time step $t$, and policy $\pi_{k}$) defined as

$V_{t}^{(\pi_{k})}(s_{t})\coloneqq\mathbb{E}\left[\sum\nolimits_{i=t}^{H}\mathcal{R}\big{(}s_{i},\pi_{k}(s_{i},i)\big{)}\right],$
(2)

where the expectation is taken with respect to states $s^{\prime}$ encountered
in the MDP. The quality of a policy at episode $k$ is characterised by the
total expected reward defined as
$u^{(\pi_{k})}=p_{0}^{\top}V^{(\pi_{k})}_{1}$, where $p_{0}$ is the
distribution of the initial states. The value function represented in Equation
2 can be rewritten as follows:

$V_{t}^{(\pi_{k})}(s_{t})=\mathcal{R}\big{(}s_{t},\pi_{k}(s_{t},t)\big{)}+P\big{(}\cdot|s_{t},\pi_{k}(s_{t},t)\big{)}^{\top}V_{t+1}^{(\pi_{k})}.$
(3)

Equation 3 is called the Bellman equation or the optimality equation, which is
often solved using the backward induction process with boundary condition
$V_{H+1}^{(\pi_{k})}\coloneqq 0$. The goal of an RL agent is then to find an
optimal trajectory which maximises the total expected reward. Thus, the
optimal policy $\pi^{*}$ is calculated through the backward induction
procedure as follows:

$\pi^{*}(s,t)=\operatorname*{argmax}_{a\in\mathcal{A}_{s}}\left\\{\mathcal{R}(s,a)+P\big{(}\cdot|s,a\big{)}^{\top}V_{t+1}^{\pi^{*}}\right\\},\,\,\,\pi^{*}(s,H)=\operatorname*{argmax}_{a\in\mathcal{A}_{s}}\mathcal{R}(s,a).$
(4)

We are now in a position to present the Modified UBEV for Stationary FHMDPs
with Moving window (MUBEV-SM) algorithm, which is the tool we use in our
approach. Algorithm 1 represents a new variant of the MUBEV algorithm
introduced in [13] (see [4] for the most recent version of MUBEV), which also
includes the key modifications of the UBEV-S algorithm presented in [14] in
order to achieve a better performance on stationary MDPs. Namely, UBEV-S
incorporates a slightly modified exploration bonus and produces tighter regret
bounds in the settings of time-independent rewards and transition dynamics
[14].

It is worth noting that MUBEV-SM includes all the modifications introduced in
the MUBEV algorithm [5, 13]. That is, (i) MUBEV-SM does not require the
estimation of transition probabilities due to the specific design of the
state-action space; (ii) in MUBEV-SM, observation of the environment is
achieved using multi-agent approach; and (iii) MUBEV-SM incorporates a better
action selection strategy. However, we introduce a few new modifications to
the MUBEV-SM algorithm as explained below.

First, the action selection strategy is slightly improved now. Specifically,
MUBEV is forced to use the default policy (e.g. SP policy) whenever all the
entries of $Q_{s}$ are equal to each other, where $Q_{s}$ denotes a vector of
$Q$-function values for state $s$. In MUBEV-SM, we explicitly compute and use
the maximum value of $Q_{s}$, i.e. $Q_{max}$ (Algorithm 1, line 14). When
$Q_{max}$ appears multiple times in $Q_{s}$, the selection strategy is as
follows: if one of the repetitions corresponds to the default action, then
action is chosen according to the default policy; otherwise, an entry is
randomly selected from those repetitions using uniform distribution, and the
associated action is chosen (Algorithm 1, line 16). Such a randomness boosts
the exploration at the beginning of the learning process.

Second, MUBEV-SM involves a new mechanism not used in [4, 5, 13], referred to
as MWRM, which is introduced as a mitigation method against stagnation that
may occur when a change in the environment happens after a long period of
stationarity, such as a long sequence of episodes with similar traffic
conditions. For the calculation of the estimated reward $\hat{r}$ in UBEV and
UBEV-S [14, 15], all reward values from the first episode are accumulated,
while in MUBEV-SM this only happens in the context of a moving window. That
is, in MUBEV-SM, a maximum of past $J$ rewards are accumulated (with $J$ the
length of the moving window), and the oldest reward is dropped every time a
new reward is collected. Additionally, the accumulation is done using weighted
rewards, and all the collected statistics are reset once the confidence bound
reaches a certain minimum value. As a result, MUBEV-SM can better cope with
environments whose reward distribution dynamically changes or includes long
periods of stationarity333Note that reward is independent of time step $t$ in
stationary MDPs, although it can vary from one episode to another.. However,
the MWRM mechanism can be disabled if required ($MWRM=\text{False}$ in
Algorithm 1).

(a)  (b)

Figure 1:  Left: The MWRM mechanism used in MUBEV-SM. Gray cells represent
initial values of $\rho$. A green cell represents the current observed reward.
Rewards observed at previous interactions are stored sequentially (blue
cells). Values of $\rho$ in orange and green cells (if moving window is
reached) are subject to discounts. Right: Schematic representation of the
state weight $w$ and the route weight $\omega$ used for computation of
rewards. $w_{s^{\prime}}$ and $w_{\bar{s}^{\prime}}$ are the state weights for
nodes $s^{\prime}$ and $\bar{s}^{\prime}$, respectively. $\omega_{s^{\prime}}$
and $\omega_{\bar{s}^{\prime}}$ are the route weights to the destination node
$s_{f}$ from nodes $s^{\prime}$ and $\bar{s}^{\prime}$, respectively.

For a better understanding of the MWRM mechanism, see Figure 1-Left in which
we illustrate this technique for a given state-action pair $(s,a)$ observed at
the $i$-th interaction with the environment (without loss of generality).
Recall that $J$ is the length of the moving window. For each $(s,a)$ pair, we
introduce an additional transition count $n^{\prime}$, and a reward memory
vector $\rho$ of length $J$ to store the collected rewards. Initial values of
$n$, $n^{\prime}$, and $\rho_{j}$ (gray cells) are set to zero for all $(s,a)$
pairs and $\forall j\in[J]$, where $n$ and $n^{\prime}$ represent the number
of times action $a$ is played in state $s$. Note that both $n$ and
$n^{\prime}$ are equal to each other and incremented by $1$ only if $i<J$. At
interaction $i$, reward $r_{i}$ is observed and $\rho$ is updated such that
$\rho_{i}=r_{i}$ (green cells) if $i\leq J$, i.e. within the moving window,
and the estimated reward $\hat{r}$ is the mean of the collected rewards. When
the moving window length is reached, i.e. $n=n^{\prime}=J$, we stop
incrementing $n$. Thus, $n(s,a)\leq J$ for all state-action pairs $(s,a)$. The
estimated reward $\hat{r}$ at interaction $i=J$ is the mean of $\rho$. For
further interactions $i>J$, the count $n^{\prime}$ is determined as a function
of the confidence bound $\phi$ as follows: $n^{\prime}(s,a)$ is reset to $1$
for all $a\in\mathcal{A}_{s}$ whenever the confidence bound for state-action
pair $(s,a)$ decreases below some predefined threshold $\phi_{min}$ (which
ensures a faster return to the default policy). Additionally, at interactions
$i>J$, the oldest value of the collected reward is dropped, all other entries
of $\rho$ are shifted one position to the left, and the new collected reward
(green cell) is inserted at the $J$-th position. Then, $\hat{r}$ is computed
as the weighted mean of $\rho$ using the corresponding weighting coefficients
$\gamma_{i}$. We refer to $\gamma$ as the vector of “forgetting” factors (or
discount vector), and it is a tuning parameter along with $\phi_{min}$.

Algorithm 1 Modified Upper Bounding the Expected Next State Value (UBEV)
Algorithm for Stationary FHMDPs with Moving Window—MUBEV-SM

1:
$\mathcal{S},\,\mathcal{A},\,P,\,\Pi_{DP};\,\,\varepsilon,\delta\in\left(0,1\right];\,\,H,\,J,\,M\in\mathbb{N};\,\,r_{max},\phi_{min}\in\mathbb{R};\,\,\gamma\in\mathbb{R}^{J};\,\,MWRM\in\\{\text{True,
False}\\}$.

2:$n(s,a)=n^{\prime}(s,a)=\hat{r}(s,a)=R(s,a)=Q(s,a)=\phi(s,a)=\hat{V}(s,t)=\rho(s,a,j)=0\,\,\forall
s\in\mathcal{S},\,\,a\in\mathcal{A}_{s},\,\,t\in[H+1],\,\,j\in[J]$.

3:$\delta^{\prime}=\delta/9,\,V_{max}=H*r_{max},\,\phi_{+}=0,\,\pi=\Pi_{DP}$.

4:for $k=1,2,3...$ do

5: for $t=H$ to $1$ do$\triangleright$ Optimistic planning loop

6:
$\hat{V}_{t+1}=\hat{V}(\cdot,t+1);\,\,\tilde{V}_{max}=\min\big{(}\max(\hat{V}_{t+1}),V_{max}\big{)}$

7: for $s\in\mathcal{S}$ do

8: $\eta=\ln\big{(}27*S*A_{s}*H/\delta^{\prime}\big{)}$

9: for $a\in\mathcal{A}_{s}$ do

10: $r=r_{max};\,\,EV=\tilde{V}_{max}$

11: if $n^{\prime}(s,a)>0$ then

12:
$\phi(s,a)=\varepsilon*\sqrt{\frac{2\ln\ln\big{(}\max\big{(}e,\,n^{\prime}(s,\,a)\big{)}\big{)}+\eta}{n^{\prime}(s,\,a)}};\,\,\hat{V}_{next}=P(\cdot,s,a)\times\hat{V}_{t+1};\,\,\bar{v}=\max(\hat{V}_{t+1})-\min(\hat{V}_{t+1})$

13:
$\bar{v}^{\prime}=\min\big{(}(H-t)*r_{max},\bar{v}+\phi_{+}\big{)};\,\,r=\min\big{(}r_{max},\hat{r}(s,a)+\phi(s,a)\big{)};\,\,EV=\min\left(\tilde{V}_{max},\hat{V}_{next}+\bar{v}^{\prime}*\phi(s,a)\right)$

14: $Q(s,a)=r+EV$

15: $Q_{max}=\max_{a\in\mathcal{A}_{s}}\big{(}Q(s,a)\big{)}$

16: if $\,\,Q_{max}$ duplicates in $Q(s,\cdot)\,\,$ then

17: if $\,\,Q_{max}=Q\big{(}s,\Pi_{DP}(s,t)\big{)}\,\,$ then
$\,\,\tilde{a}=\Pi_{DP}(s,t)\,\,$ else
$\,\,\tilde{a}\sim\mathcal{U}\big{(}\\{a\,|\,Q(s,a)=Q_{max}\\}\big{)}$

18: else$\,\,\tilde{a}=\operatorname*{argmax}_{a\in\mathcal{A}_{s}}Q(s,a)$

19: $\pi_{k}(s,t)=\tilde{a};\,\,\hat{V}(s,t)=Q(s,\tilde{a})$

20: if $\,\,n^{\prime}(s,\tilde{a})>0\,\,$ then
$\,\,\tilde{\phi}=\varepsilon*\sqrt{\frac{2\ln\ln\big{(}\max\big{(}e,\,n^{\prime}(s,\,\tilde{a})\big{)}\big{)}+\eta}{n^{\prime}(s,\,\tilde{a})}};\,\,\phi_{+}=\max\big{(}4*\sqrt{S}*H^{2}*\tilde{\phi},\phi_{+}\big{)}$

21:
$\tilde{s}=\big{[}s_{1}^{(1)},...,s_{1}^{(M)}]\sim\mathcal{U}\big{(}1,S\big{)},\,s_{1}^{(i)}\neq
s_{1}^{(j)}\,\,\,\forall i,j\in[1,M]$

22: for $m=1$ to $M$, $t=1$ to $H$ do $\triangleright$ Execute policy for one
episode

23:
$a_{t}^{(m)}=\pi_{k}(s_{t}^{(m)},t);\,\,\bar{a}_{t}^{(m)}=\Pi_{DP}(s_{t}^{(m)},t);\,\,s_{t+1}^{(m)}\sim
P(\cdot|s_{t}^{(m)},a_{t}^{(m)});\,\,\bar{s}_{t+1}^{(m)}\sim
P(\cdot|s_{t}^{(m)},\bar{a}_{t}^{(m)})$

24:
$r^{(m)}_{t}=\mathcal{R}\big{(}s_{t}^{(m)},s_{t+1}^{(m)},\bar{s}_{t+1}^{(m)});\,\,R\big{(}s_{t}^{(m)},a_{t}^{(m)}\big{)}+=r^{(m)}_{t}$
$\triangleright$ Call to Function 1, accumulate rewards

25: if $MWRM=\text{True}$ then$\triangleright$ Employ moving window

26: $j_{m}\coloneqq
n\big{(}s_{t}^{(m)},a_{t}^{(m)}\big{)};\,\,j_{m}^{\prime}\coloneqq
n^{\prime}\big{(}s_{t}^{(m)},a_{t}^{(m)}\big{)}$

27: if $j_{m}\leq J$ then

28:
$\rho\big{(}s_{t}^{(m)},a_{t}^{(m)},j_{m}\big{)}=r^{(m)}_{t};\,n\big{(}s_{t}^{(m)},a_{t}^{(m)}\big{)}\tiny{++};\,n^{\prime}\big{(}s_{t}^{(m)},a_{t}^{(m)}\big{)}\tiny{++};\,\hat{r}\big{(}s_{t}^{(m)},a_{t}^{(m)}\big{)}=R\big{(}s_{t}^{(m)},a_{t}^{(m)}\big{)}\big{/}n\big{(}s_{t}^{(m)},a_{t}^{(m)}\big{)}$

29: else

30:
$R\big{(}s_{t}^{(m)},a_{t}^{(m)}\big{)}=R\big{(}s_{t}^{(m)},a_{t}^{(m)}\big{)}-\rho\big{(}s_{t}^{(m)},a_{t}^{(m)},1\big{)}$
$\triangleright$ “Forget” the oldest reward

31:
$\rho\big{(}s_{t}^{(m)},a_{t}^{(m)},\cdot\big{)}=\rho\big{(}s_{t}^{(m)},a_{t}^{(m)},\cdot\big{)}_{2:J}+\\!\\!\\!+\big{\\{}r^{(m)}_{t}\big{\\}}$
$\triangleright$ Incorporate the newest reward

32:
$\hat{r}\big{(}s_{t}^{(m)},a_{t}^{(m)}\big{)}=\Big{(}\gamma\times\rho\big{(}s_{t}^{(m)},a_{t}^{(m)}\big{)}\Big{)}\big{/}J$

33: if $\,\,\phi\big{(}s_{t}^{(m)},a_{t}^{(m)}\big{)}\leq\phi_{min}\,\,$ then
$\,\,n^{\prime}\big{(}s_{t}^{(m)},a\big{)}=1\,\,\forall
a\in\mathcal{A}_{s}\,\,$ else
$\,\,n^{\prime}\big{(}s_{t}^{(m)},a_{t}^{(m)}\big{)}++$

34: else$\triangleright$ Do not apply moving window

35:
$n\big{(}s_{t}^{(m)},a_{t}^{(m)}\big{)}++;\,\,n^{\prime}\big{(}s_{t}^{(m)},a_{t}^{(m)}\big{)}++;\,\,\hat{r}\big{(}s_{t}^{(m)},a_{t}^{(m)}\big{)}=R\big{(}s_{t}^{(m)},a_{t}^{(m)}\big{)}\big{/}n\big{(}s_{t}^{(m)},a_{t}^{(m)}\big{)}$

The notation for MUBEV-SM (Algorithm 1) is as follows. $\mathcal{S}$ is the
set of states; $\mathcal{A}_{s}$ is the set of allowable actions in state $s$,
and thus, $\mathcal{A}=\bigcup_{s\in\mathcal{S}}\mathcal{A}_{s}$ is the total
action set; $S$, $A_{s}$ and $A$ denote cardinality of finite sets
$\mathcal{S}$, $A_{s}$ and $\mathcal{A}$, respectively; $P$ is a 3-dimensional
tensor of predefined transition probabilities; $\Pi_{DP}$ is the default path
(DP) policy, computed using an SP algorithm as the initial estimate for $\pi$;
$\varepsilon$ is the exploration bonus, and $\delta$ is the failure
probability (see [15] for details); $H$ is the length of the MDP’s time
horizon; $J$ is the size of moving window; $M$ is the number of MUBEV-SM
tokens; $r_{max}$ is the maximum reward an agent can receive per transition;
$\phi_{min}$ is the minimum value of the confidence bound which determines
when to reset the collected statistics (i.e. counts $n^{\prime}$); $\gamma$ is
a vector of discounts for the rewards stored in the reward memory $\rho$;
$MWRM$ is a boolean flag to enable/disable MWRM; $n(s,a)$ and
$n^{\prime}(s,a)$ count the times action $a$ is played in state $s$ (see
previous paragraph for details); $\hat{r}(s,a)$ and $R(s,a)$ are the
normalised and accumulated rewards in state $s$ under action $a$,
respectively; $Q(s,a)$ is the Q-function [15] for state $s$ and action $a$;
$\phi(s,a)$ is the confidence bound for state $s$ and action $a$;
$\hat{V}(s,t)$ is the value function from time step $t$ for state $s$;
$\rho(s,a,j)$ represents a reward stored in a 3-dimensional tensor of reward
memory $\rho$ for state $s$ and action $a$ at memory epoch $j$;
$\delta^{\prime}$ represents a scaled failure tolerance (see [15] for
details); $V_{max}$ is the maximum value of the value function for next
states; $\phi_{+}$ is the correction term which keeps track of the largest
confidence bound for the least visited state-action pair under the policy
$\pi$ [14]; $\hat{V}(\cdot,t+1)$ and $P(\cdot,s,a)$ denote vectors of length
$S$, and $Q(s,\cdot)$ is interpreted as a vector of length $A_{s}$; $e$ is the
Euler’s number; $\eta$, $r$, $EV$, $\bar{v}$ and $\bar{v}^{\prime}$ are
auxiliary variables; $\bar{v}$ represents the range of vector $\hat{V}_{t+1}$
[14]; vector $\tilde{s}$ is a vector of initial states of MUBEV-SM tokens,
which is uniformly sampled in range from $1$ to $S$ with no repeated entries;
$m$ is the index of an agent (a vehicle with a token) that interacts with the
environment each time step $t$, and receives reward $r_{t}^{(m)}$ determined
by the reward function defined in Function 1.

###  3.2 Reward function

The reward function (Function 1) returns the total reward $r_{t}$ at time $t$,
which implicitly incorporates two objectives, namely travel distance and
parking availability.

1:$\alpha\in\left[\alpha_{min},1\right];s_{t},s_{t+1},\bar{s}_{t+1}\in\mathcal{S};\beta_{1},\beta_{2},r_{max}\in\mathbb{R}$.

2:$r_{t}$.

3:$s_{t}\coloneqq s,s^{\prime}\coloneqq
s_{t+1},\bar{s}^{\prime}\coloneqq\bar{s}_{t+1}$

4:function $\mathcal{R}$($s,s^{\prime},\bar{s}^{\prime}$)

5: Get the state lengths $\lambda_{s^{\prime}}$ and
$\lambda_{\bar{s}^{\prime}}$, and the total roadside parking capacity
$c_{\bar{s}^{\prime}}$.

6: Observe the current roadside parking availability $\nu_{s^{\prime}}$.

7:
$w_{\bar{s}^{\prime}}=\mathcal{W}_{G}(\alpha,\bar{s}^{\prime},\lambda_{\bar{s}^{\prime}},c_{\bar{s}^{\prime}}),\,\,$
$w_{s^{\prime}}=\mathcal{W}_{G}(\alpha,s^{\prime},\lambda_{s^{\prime}},\nu_{s^{\prime}})$
$\triangleright$ Compute weights for $\bar{s}^{\prime},s^{\prime}$,
respectively, using Equation 1

8: if $w_{\bar{s}^{\prime}}\neq w_{s^{\prime}}$ then$\triangleright$
Uncertainty has been encountered

9: Get route total weights $\omega_{s^{\prime}}$ and
$\omega_{\bar{s}^{\prime}}$.
$\Omega_{s,\bar{s}^{\prime}}=\omega_{\bar{s}^{\prime}}+w_{\bar{s}^{\prime}};\,\,\Omega_{s,s^{\prime}}=\omega_{s^{\prime}}+w_{s^{\prime}};\,\,r_{t}=1-\frac{\Omega_{s,s^{\prime}}}{\Omega_{s,\bar{s}^{\prime}}}$

10: if $\,s^{\prime}=\bar{s}^{\prime}\,$ then $\,r_{t}=\beta_{1}*r_{t}$
$\triangleright$ Agent takes the DP action

11: else$\,\,\,r_{t}=\beta_{2}*r_{t}$ $\triangleright$ Agent takes an
alternative action

12: else$\,\,\,r_{t}=r_{max}$ $\triangleright$ There is no uncertainty return
$r_{t}$

Function 1 The Reward Function

Function 1 receives three arguments: (i) $s_{t}$, the previous state of an
agent; (ii) $s_{t+1}$, the current state; and (iii) $\bar{s}_{t+1}$, the state
to which the agent would “jump” if the default action in state $s_{t}$ were
taken. To begin with, $\lambda_{s^{\prime}}$ and $\lambda_{\bar{s}^{\prime}}$
are the lengths of road links included in states $s^{\prime}$ and
$\bar{s}^{\prime}$, respectively, $c_{\bar{s}^{\prime}}$ is the total on-
street parking capacity for state $\bar{s}^{\prime}$, and $\nu_{s^{\prime}}$
is the current parking availability in state $s^{\prime}$. Note that the state
length $\lambda$ can be easily computed for a given road network. Also, as we
assume that total roadside PS capacity for the road network is available (e.g.
via digital map provider), and thus the total state capacity $c$ is known for
each state. However, PSA is subject to real-time measurements (e.g. obtained
via [11, 12]). $w_{\bar{s}^{\prime}}$ is the state weight if the agent would
transit from state $s$ to state $\bar{s}^{\prime}$ in which all the on-street
parking stops are treated as available. In contrast, $w_{s^{\prime}}$
represents the state weight of the actual transition from state $s$ to state
$s^{\prime}$ at time $t$. If uncertainties occur in the environment, then the
weights computed in line 7 are not equal to each other, and thus the DP policy
has to be amended. In this case, the reward is computed using the total route
weight of two (in general different) routes. Further,
$\Omega_{s,\bar{s}^{\prime}}$ and $\Omega_{s,s^{\prime}}$ are the total
weights of routes from $s$ to the destination state $s_{f}$ via
$\bar{s}^{\prime}$ and $s^{\prime}$, respectively, and
$\omega_{\bar{s}^{\prime}}$ and $\omega_{s^{\prime}}$ are the total weights of
routes from $\bar{s}^{\prime}$ and $s^{\prime}$ to the destination $s_{f}$,
respectively (see Figure 1-Right). As a result,
$\Omega_{s,\bar{s}^{\prime}}=\omega_{s}$, and since
$\Omega_{s,\bar{s}^{\prime}}<\Omega_{s,s^{\prime}}$, reward computed in line 7
always satisfies $r_{t}<0$. Finally, $\beta_{1}$ and $\beta_{2}$ are used to
speed up the learning process of searching for a detour and returning to the
DP route respectively, and their values may depend on many factors including
values of $\alpha$ and $r_{max}$ or the road network itself. To achieve a
better and faster performance of the algorithm, $\beta_{1}$ and $\beta_{2}$
need to be increasing functions of $\alpha$ as illustrated in the following
section. Finally, if a given state is not affected by uncertainties, the
agents receive the maximum reward $r_{max}$ available per transition.

##  4 Numerical evaluation

In this section we are interested in evaluating a route recommender system for
smart parking, for which a number of MUBEV-SM tokens is distributed to emulate
real vehicles probing the uncertain environment. These tokens are passed from
vehicle to vehicle using the same DLT framework introduced in [4, 5], where
vehicles in possession of tokens are permitted to write data to the DLT. The
token passing mechanism is dictated by both the MUBEV-SM algorithm and the DLT
architecture and can be implemented using a MEC/cloud-based service (also
described in [4, 5]). Once the unknown environment has been ascertained, route
recommendations can be accessed via a smart parking app by a variety of users
interested in optimal routing.

Figure 2: The road network used in the experiments comes from an urban area
close to the city centre in Würzburg, Germany. In the OpenStreetMap network of
this area555https://www.openstreetmap.org/#map=14/49.7830/9.9553&layers=T, PSs
have been tagged, and thus on-street PSs can be easily imported to a SUMO road
network file. Some road segments of interest are highlighted, representing
either origins (O), destinations (D), or sources of full parking areas (F).

To evaluate MUBEV-SM, a number of numerical experiments was designed using the
traffic simulator SUMO [16]. Interactions with running simulations (including
those to measure roadside parking occupancy) are accomplished using the Python
programming language and SUMO packages Sumolib and TraCI. The generic settings
for our simulations are as follows:

  * •

The urban area in Würzburg, Germany shown in Figure 5 is used for traffic
simulations in all our experiments. Once all the U-turns are removed666U-turns
may lead to undesirable recurrent attempts to use DP policy. from the network
file, the road network results in a graph with 861 nodes (states).

  * •

A number of road links are selected as origins, destinations, and sources of
full on-street parking areas (see Figure 5 for details). Specifically, we use:
origin $O1$ and full parking area (Full-PA) $F1$ in Experiment 1 and 2; origin
$O2$ and Full-PA $F2$ in Experiment 3; destination $D1$ in all our
experiments.

  * •

In all our simulations, we create and use a vehicle type based on the default
SUMO vehicle type with maximum speed 118.8 km/h and
impatience777https://sumo.dlr.de/wiki/Definition_of_Vehicles_Vehicle_Types_and_Routes.
equal to 0.5. When these vehicles are in possession of a token, they become
virtual MUBEV-SM vehicles. In Experiment 1 and 2, these vehicles have no
parking stops assigned. In Experiment 3, every 100th vehicle added to the
simulation is assigned to a random parking stop along its route. For the
generation of Full-PAs, we release a number of cars of the aforementioned
vehicle type and populate the selected on-street parking lanes with them.

  * •

For all our experiments, we add an average of 250 background cars with random
routes, which can carry tokens if required.

  * •

DPs are calculated using Dijkstra’s SP algorithm
(SciPy888https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csgraph.dijkstra.html
implementation).

  * •

A state of an RL agent is considered to be a set of road links joined into one
state using the merging technique introduced in [4]. We refer to a token trip
as an RL episode. Before policy execution, the tokens are uniformly
distributed on the road network. The length of the MDP’s time horizon $H$ is
set to 50. $H$ determines a maximum number of allowed road links that each
token can traverse over a given episode. If tokens do not reach a specified
destination within this restriction, their trips will be declared incomplete
(i.e. unsuccessful).

Regarding the input parameters of Algorithm 1 and Function 1, in all our
experiments we set $\varepsilon=\delta=r_{max}=1$, $J=10$, $\phi_{min}=1.25$.
The components of vector $\gamma$ are computed as follows:
$\gamma_{j}=1-(J-j)*(1-\gamma_{min})/(J-1)$, $j\in[J]$, with $\gamma_{min}$
empirically set to $0.1$. In Experiment 1, we set $\beta_{1}=\beta_{2}=20$ for
simulations with $\alpha=0.66$, and $\beta_{1}=\beta_{2}=60$ for
$\alpha=0.80$. In Experiment 2 and 3, we set $\beta_{1}=\beta_{2}=20$ and
$\alpha=\alpha_{min}=0.66$. Any specific additional setting for each
individual experiment will be described in the corresponding subsection below.

###  4.1 Experiment 1: Optimal route estimation under uncertainty using
different values of $\alpha$

The objectives of this experiment is to illustrate that our DLT-enabled RL
approach is able to determine a simple uncertain environment using a single
token for different values of $\alpha$. For a chosen OD pair, we artificially
simulate a Full-PA on some road links along DP at various time intervals. In
this scenario, our goal is to show that the token-enabled MUBEV-SM algorithm
can distinguish between Free Parking Intervals (Free-PIs) and Full Parking
Intervals (Full-PIs), and, in the latter case, find the next optimal route for
the selected OD pair.

Particularly, this illustrative experiment is designed as follows. A single
MUBEV-SM token is used over each episode of the learning process to collect
data used to update the MDP’s policy. This token has a fixed OD pair, namely
$\\{O1,D1\\}$ as marked in Figure 5. Additionally, we choose $F1$ (a set of
road segments which belong to DP for the specified OD pair) and generate a
Full-PA in it at different time intervals. Over each episode, we start the
token from $O1$ and ask it to travel to $D1$, keeping a record of its
performance in terms of route characteristics, such as the total weight
(computed using Equation 1), the number of available on-street PSs, and the
travel distance (route length), regardless of its success in attempting to
reach $D1$. During Full-PIs, the total weight of the token route enlarges
while parking availability decreases (this is annotated with “DP full” for the
first Full-PI in Figure 3). MUBEV-SM is expected to provide an alternative
path (AP) which is an optimal solution leading to a minimum possible value of
total weight during Full-PIs (as also annotated in Figure 3), and that DP
routing would eventually be advised in all other situations.

(a)  (b)

Figure 3:  Experiment 1: Performance indices of MUBEV-SM with MWRM during the
learning process on a changing environment using one token and fixed OD pair
$\\{O1,D1\\}$ for $\alpha=0.66$ (Left) and $\alpha=0.80$ (Right). Annotations
are only provided for the first part of the learning process (DP: default
path, AP: alternative path).

The results of a single realisation of the previously described experiment are
depicted in Figure 3 for different values of $\alpha$. As it can be noted, the
token succeeds in avoiding Full-PAs once they are created (which can be better
observed in the travel distance plots), and learns an AP route with a smaller
total weight than that of DP during Full-PIs (annotated as “DP full”). The
token also successfully returns to DP once Full-PAs are cleared. Both actions,
avoiding and returning, happen within a reasonably small number of episodes.
Additionally, thanks to MWRM, this rapid adaptation is present even at the
beginning of the learning process, and the performance of MUBEV-SM is nearly
uniform as time passes. These observations validate our expectations about
MUBEV-SM concerning its ability to adapt rapidly to uncertain environments.

When comparing the right and left panels of Figure 3, we can notice that the
length of DP, and the total parking capacity along it, are larger for a
smaller value of $\alpha$ (left panel) than the corresponding values for a
larger $\alpha$ (right panel). This indicates, as expected, that smaller
values of $\alpha$ generally lead to longer DP routes. Consequently, it would
be impractical to consider an extremely small value of $\alpha$ (even if it
satisfies $\alpha\in[\alpha_{min},1]$), as it could lead to a dramatic
increase in travel distances/times. As it can also be seen in Figure 3, our
MUBEV-SM performs as expected for both values of $\alpha$.

###  4.2 Experiment 2: Comparative analysis of MUBEV-SM with and without MWRM

In this second experiment, we want to evaluate the performance of MUBEV-SM
with and without MWRM using a single token and non-uniform Full-PIs. For this,
we use a similar setup as in Experiment 1 (i.e. same OD pair, Full-PA on $F1$,
and one token probing the environment). To begin with, the token starts at
$O1$ and is asked to travel to $D1$ each episode, keeping a record of the
performance indices as mentioned in the previous subsection. We intentionally
do not introduce a Full-PA for a long period of time, allowing the token to
collect lots of “positive” rewards while traveling along DP, after which we
generate a Full-PA over a reasonable short time interval. Afterwards, we allow
a short Free-PI followed by a Full-PA on $F1$ again, but over a significantly
longer period of time. The results of this comparative analysis for a single
realisation are shown in Figure 4.

(a)  (b)

Figure 4: Experiment 2: Comparative performance of MUBEV-SM without (Left) and
with (Right) MWRM on a changing environment with non-uniform Full-PIs, using a
single token and fixed OD pair. Annotations are only provided for the first
part of the learning process (DP: default path, AP: alternative path).

From Figure 4-Left, we can notice that the learning process using MUBEV-SM
without MWRM is notably slow during the first Full-PI (better observed in the
travel distance plot, left panel), meaning that a large number of episodes are
required to learn an AP. On the contrary, the algorithm is able to provide an
optimal detour AP much faster (in terms of fewer episodes) if MWRM is applied
as seen in Figure 4-Right. The latter configuration also allows for a better
performance during the second Full-PI: MUBEV-SM is able to learn an AP route
in a small number of episodes, and such an AP solution is preferred over the
Full-PI. It must be noted that without MWRM, the system mostly selects DP
routing over the whole learning process, which can be seen in the lowest
subplot in Figure 4-Left. The reason of such performance is that the system
without MWRM has collected lots of “positive” rewards and statistics during
the first (long) Free-PI.

Finally, it can be observed that returning to DP (once the Full-PAs are
cleared) takes slightly longer with MWRM than without. However, with MWRM, the
system does not require much tuning of its parameters for each particular
uncertainty, and performs significantly better in general. Note that
experiments with a single token are useful to analyse the performance of
MUBEV-SM in uncertain environments, but more elaborated experiments must be
performed for a comprehensive validation. Experiment 3 addresses this.

###  4.3 Experiment 3: Route recommendations from a MUBEV-SM-based system, and
speedup in learning

Along with an illustration of how MUBEV-SM route recommendations can be
delivered to a number of drivers, this experiment is also designed to evaluate
the performance of the system as a function of the number of tokens used to
collect the data to update the MDP’s policy. The setup for this experiment is
as follows. Over each episode, we release a set of tokens starting at
different origins according to a uniform spatial distribution. All these
tokens have a common destination $D1$. In parallel, we analyse the performance
of a test (non-MUBEV-SM) vehicle starting at origin $O2$ every single episode.
The test vehicle is asked to travel to destination $D1$ using recommendations
provided by a MUBEV-SM-based routing system, whose recommendations are defined
as follows: (i) a route connecting $\\{O2,D1\\}$ is calculated based on the
estimated MDP’s policy; and (ii) whenever a complete route for $\\{O2,D1\\}$
cannot be calculated using the MUBEV-SM policy (e.g. the destination is not
reached within the length horizon $H$), the latest valid recommendation is
reused. A Full-PA is generated on the set of road links $F2$ (which belong to
the DP for $\\{O2,D1\\}$) during episodes 25 to 45. It is worth noting that a
new test car is released at the end of each episode as soon as the policy is
updated. Figure 5 depicts the results for this experiment obtained from 50
different realisations.

(a)  (b)

Figure 5: Experiment 3: Left: Median values of route attributes for a test
vehicle using recommendations from a MUBEV-SM-based routing system involving
multiple tokens. Each datapoint represents the median value of 50 different
realisations of the experiment. Right: Average learning speed using multiple
tokens. The error bars depict the 95% confidence interval for the mean. Dashed
lines were obtained using exponential curve fitting.

As it can be concluded from Figure 5-Left, the recommender system has a
remarkable performance for different number of participating tokens. However,
the number of tokens directly affects the convergence rate of the MUBEV-SM
algorithm: the more tokens involved, the faster the learning process. A more
conclusive relationship between the number of participating tokens and the
average number of episodes required to learn a new traffic condition (either
free or full parking areas) is depicted in Figure 5-Right. This relationship
reflects an exponential-like decay for learning of Full-PAs, and a nearly
linear decrease for returning to DP.

##  5 Conclusion

We propose a routing system design based on DLT and multi-agent RL technique
to solve the on-street parking problem. The proposed approach is applied to
estimate the best routes in terms of short travel distance and large number of
available roadside PSs, and the design allows drivers to choose their
preferred objectives. A moving window mechanism is also incorporated into the
underlying RL algorithm to better cope with non-stationary environments.
Experimental results have proven the efficiency of our system in improving the
parking problems. For future work, we plan to include more objectives into our
optimisation problem (e.g. travel time, walking distance), and extend the
range of allowable values of $\alpha$.

##  6 Acknowledgements

This work has been supported by Science Foundation Ireland (SFI) under grant
No. 16/IA/4610. Part of this work has received funding from the EU’s Horizon
2020 research & innovation programme via the ICT4CART project under grant
agreement No. 768953. Content reflects only the authors’ view and the European
Commission is not responsible for any use that may be made of the information
it contains.

## References

  * [1] D. C. Shoup, “Cruising for parking,” Transport Policy, vol. 13, no. 6, pp. 479–486, 2006.
  * [2] C. Tang et al., “Towards Smart Parking Based on Fog Computing,” IEEE Access, vol. 6, pp. 70172–70185, 2018.
  * [3] M. Buchholz et al., “Enabling automated driving by ICT infrastructure: A reference architecture,” in Proceedings of TRA’20, Helsinki (conference cancelled), 2020.  doi: 10.18725/OPARU-26023.
  * [4] R. Overko et al., “Spatial Positioning Token (SPToken) for Smart Mobility,” IEEE Transactions on Intelligent Transportation Systems, pp. 1–14, 2020.
  * [5] R. Overko et al., “Spatial Positioning Token (SPToken) for Smart Mobility,” in Proceedings of IEEE ICCVE’19, pp. 1–6, 2019.
  * [6] K. Garlichs et al., “Leveraging the Collective Perception Service for CAM Information Aggregation at Intersections,” 2020.  Accepted at IEEE VTC2020-Fall.
  * [7] W. Al Amiri et al., “Privacy-Preserving Smart Parking System Using Blockchain and Private Information Retrieval,” in Proceedings of SmartNets’19, pp. 1–6, 2019.
  * [8] J. Ni et al., “Privacy-Preserving Smart Parking Navigation Supporting Efficient Driving Guidance Retrieval,” IEEE Transactions on Vehicular Technology, vol. 67, no. 7, pp. 6504–6517, 2018.
  * [9] M. Khalid et al., “A Reinforcement Learning based Path Guidance Scheme for Long-range Autonomous Valet Parking in Smart Cities,” in Proceedings of ComNet’20, pp. 1–7, 2020.
  * [10] S. Padakandla et al., “Reinforcement Learning in Non-Stationary Environments,” Applied Intelligence, vol. 50, pp. 3590–3606, 2020.
  * [11] M.-C. Wu et al., “Early detection of vacant parking spaces using dashcam videos,” in Proceedings of AAAI’19, vol. 33, pp. 9613–9618, 2019.
  * [12] F. Bock et al., “On-street parking statistics using lidar mobile mapping,” in Proceedings of ITSC’15, pp. 2812–2818, IEEE, 2015.
  * [13] R. Overko et al., “Reinforcement Learning Augmented Optimization for Smart Mobility,” in Proceedings of IEEE CDC’19, pp. 1286–1292, 2019.
  * [14] A. Zanette et al., “Problem Dependent Reinforcement Learning Bounds Which Can Identify Bandit Structure in MDPs,” in Proceedings of MLR’18, vol. 80, pp. 5747–5755, PMLR, 10–15 Jul 2018.
  * [15] C. Dann et al., “Unifying PAC and regret: Uniform PAC bounds for episodic reinforcement learning,” in Proceedings of NeurIPS’17, pp. 5713–5723, 2017.
  * [16] P. A. Lopez et al., “Microscopic Traffic Simulation using SUMO,” in Proceedings of ITSC’18, pp. 2575–2582, IEEE, 2018.


# Projection predictive inference beyond Exponential Family models

Alejandro Catalina

###### Abstract

Given a reference model that includes all the available variables, projection
predictive inference replaces its posterior with a constrained projection
including only a subset of all variables. We extend projection predictive
inference to enable computationally efficient variable and structure selection
in models outside the exponential family. By adopting a latent space
projection predictive perspective we are able to: 1) propose a unified and
general framework to do variable selection in complex models while fully
honouring the original model structure, 2) properly identify relevant
structure and retain posterior uncertainties from the original model, and 3)
provide an improved approach also for non-Gaussian models in the exponential
family. We demonstrate the superior performance of our approach by thoroughly
testing and comparing it against popular variable selection approaches in a
wide range of settings, including realistic data sets. Our results show that
our approach successfully recovers relevant terms and model structure in
complex models, selecting less variables than competing approaches for
realistic datasets.

Latent Space Projection Predictive Inference

Alejandro Catalina1 Paul Bürkner2 Aki Vehtari1 1 Helsinki Institue for
Information Technology, HIIT, Aalto University, Finland 2 Cluster of
Excellence SimTech, University of Stuttgart, Germany

##  1 Introduction

Variable and structure selection plays an important role in a robust Bayesian
workflow (Gelman et al., 2020). While variable selection has been extensively
studied and successfully applied for models in the exponential family
(Koopman, 1936), it has not received the same attention regarding models
outside of the exponential family (e.g., advanced ordinal, count or time-to-
event (also known as survival) data distributions), despite their important
applications in different fields (Kelter, 2020; Nagler, 1994; Bürkner and
Vuorre, 2019; Barron, 1992).

We propose an efficient, stable, and information theoretically justified
method to make variable selection for non-normal observation models in or
beyond the exponential family. The main benefits of the proposed latent space
projection predictive inference are:

  1. 1.

we enable the projection predictive variable and structure selection for
models outside the exponential family while honouring the original model
structure and its predictive uncertainty,

  2. 2.

we obtain more stable projections for non-Gaussian exponential family models,

  3. 3.

we demonstrate the superior performance of our method as compared to state-of-
the-art competitors in both simulated and real-world scenarios,

  4. 4.

we provide a ready-to-use open source implementation of the new methods.

##  2 Projection predictive inference

Given posterior draws $\\{\lambda_{*}^{(s)}\\}_{s=1}^{S}\sim
p(\lambda\mid\mathcal{D})$ from a reference model with data
$\mathcal{D}=\\{X,y\\}$, outcome observations $y$ and predictor variables $X$,
projection predictive inference (Piironen et al., 2020b; Catalina et al.,
2020) learns a projection $q_{\bot}(\lambda)$ containing only a subset of
variables that matches the reference predictive performance as close as
possible. The solution is given by the minimiser of the Kullback-Leibler (KL)
divergence from the reference model to the projection predictive distributions
(Dupuis and Robert, 2003). Let $p(\tilde{y}\mid{\cal D})$ be the reference
predictive distribution and $\\{\lambda_{\bot}^{(s)}\\}_{s=1}^{S}\sim
q_{\bot}(\lambda)$ be draws from the projection, then:

KL $\displaystyle\left(p\left(\tilde{y}\mid{\cal D}\right)\parallel
q_{\bot}\left(\tilde{y}\right)\right)$
$\displaystyle=-\mathbb{E}_{\lambda_{*}}\left(\mathbb{E}_{\tilde{y}\mid\lambda_{*}}\left(\log\mathbb{E}_{\lambda_{\bot}}\left(p\left(\tilde{y}\mid\lambda_{\bot}\right)\right)\right)\right)+\text{C}.$
(1)

As the integrals involved are in most cases computationaly infeasible, Goutis
and Robert (1998) suggest to change the order of integration and minimisation
and solve the optimisation for each posterior draw separately

$\text{arg}\max_{\lambda}\mathbb{E}_{\tilde{y}\mid\lambda_{*}^{(s)}}\left(\log\mathbb{E}_{\lambda_{\bot}}\left(p\left(\tilde{y}\mid\lambda\right)\right)\right).$
(2)

Piironen et al. (2020b) proposed a further speed-up by first clustering the
posterior draws $\lambda_{*}^{(s)}$ and then solving the optimisation
individually for each resulting cluster centre
$\\{\lambda_{*}^{(c)}\\}_{c=1}^{C}$.

For models in the exponential family distribution, Eq. 2 coincides with
computing maximum likelihood estimates under the projection model as

$\lambda_{\bot}=\arg\max_{\lambda}\sum_{i=1}^{N}\mu_{i}^{*}\xi_{i}(\lambda)-B(\xi_{i}(\lambda)),$
(3)

where $\mu^{*}=\mathbb{E}_{\tilde{y}\mid\lambda_{*}}(\tilde{y})$ are mean
predictions of the reference model, $\xi_{i}$ are the natural parameters for
the $i$th observation and $B(\cdot)$ is a function of the natural parameters.
These maximum likelihood estimates can be efficiently computed by solving
penalised iteratively reweighted least squares (PIRLS; Marx, 1996). If non-
constant, the projected scale parameter of the exponential family model is
then obtained as

$\phi_{\bot}=\text{arg}\max_{\phi}\sum_{i=1}^{N}\left(\dfrac{r_{i}(\lambda_{\bot})}{A(\phi)}+\mathrm{E}_{\tilde{y}_{i}\mid\lambda_{*}}(H(\tilde{y}_{i},\phi))\right),$
(4)

where $A,H$ are family-specific functions and
$r_{i}(\lambda_{\bot})=\mu_{i}^{*}\xi_{i}(\lambda_{\bot})-B(\xi_{i}(\lambda_{\bot}))$
does not depend on $\phi$.

##  3 Latent space projective inference

The equivalence between maximum likelihood estimates and KL minimiser does not
hold for models outside the exponential family. For data ${\cal D}=\\{X,y\\}$,
parameters $\lambda$ and inverse link function $g$, we assume the general
model formulation

$y\sim p\left(\mu,\phi\right),\quad\mu=g(\eta),\quad\eta\sim
p\left(\lambda\mid X\right),$ (5)

where $\eta$ is commonly called the latent predictor (McCullagh and Nelder,
1989). In Fig. 1 we show the difference between the linear and the transformed
predictor spaces for two common models, a Poisson model with a log link
function and an ordinal cumulative model (Bürkner and Vuorre, 2019) with logit
link function.

(a) Latent linear predictor, transformed predictor and response space for a
Poisson model with log link function.

(b) Latent linear predictor, transformed predictor and response space for an
Ordinal cumulative model with probit link function.

Figure 1: Latent, transformed and response space representation for different
models.

We propose to reformulate the projection problem by solving the KL
minimisation in the latent predictive space $p(\tilde{\eta}\mid\lambda_{*})$
as:

KL $\displaystyle\left(p\left(\tilde{\eta}\mid{\cal D}\right)\parallel
q_{\bot}\left(\tilde{\eta}\right)\right)=$
$\displaystyle-\mathbb{E}_{\lambda_{*}}\left(\mathbb{E}_{\tilde{\eta}\mid\lambda_{*}}\left(\log\mathbb{E}_{\lambda_{\bot}}\left(p\left(\tilde{\eta}\mid\lambda_{\bot}\right)\right)\right)\right)+\text{C}.$

###  3.1 Non exponential family case

The distribution $p(\eta\mid\lambda)$ is model-dependent and in non-
exponential family case doesn’t, in general, have a nice closed form. We
propose to approximate $p(\eta\mid\lambda)$ with a Gaussian distribution, a
choice motivated by several reasons: 1) as potential boundaries on $\mu$ are
enforced only via the inverse link $g$, $\eta$ itself is unbounded so that its
support matches the support of a Gaussian, 2) the Gaussian distribution
belongs to the exponential family and is thus computationally easy to handle,
and 3) the model on $\eta$ is typically additive, such that commonly chosen
Gaussian priors on the latent parameters $\lambda$ make the implied
distribution of $\eta$ closer to Gaussian as well. Now that the transformed
projection is again in the exponential family, we can apply the approach
outlined in Section 2 to solve the projection.

Although the dispersion of the latent space is typically not a model
parameter, we still need to approximate it to compute its projection in Eq. 4.
The approximate dispersion is model dependent, and may sometimes even be known
analytically.

###  3.2 Exponential family case

The latent projection formulation also works for models in the exponential
family. For Gaussian observation models, it coincides exactly with the
original framework, and therefore no improvement is gained with the latent
approach.

For non-Gaussian exponential family models, the original framework computes an
approximate solution to Eq. 1 via PIRLS (as in Eq. 3), which often results in
unstable solutions for models with complex structure or link functions
(Catalina et al., 2020). In contrast, the latent approach computes an
approximation on the latent space, removing the complexity of the link
function and the response model. Experimental results show that the latent
approach results in significant improvements also for non-Gaussian exponential
family models.

##  4 Related work

For models in the exponential family, variable selection has a strong presence
in the literature. Some methods perform variable selection by optimising a
penalised likelihood formulation (Tibshirani, 1996; Zou and Hastie, 2005;
Friedman et al., 2010a; Candes and Tao, 2007; Breiman, 1995; Fan and Li,
2001), while at the same time trying to select a subset of relevant variables.
These methods impose a penalisation on large coefficients, effectively driving
some of them towards zero, depending on the choice of regularization. For more
details, we refer the reader to Hastie (2015). These methods suffer from
several drawbacks. First, by dealing with the estimation of the model and the
selection of variables at the same time, they often result in suboptimal
solutions (Piironen and Vehtari, 2017a). Second, they are derived for specific
likelihoods, and generalising them to other models is difficult, if at all
possible. Most notably, these methods cannot perform variable selection on
group-specific parameters in hierarchical models (Catalina et al., 2020).

These approaches have been generalized to specific likelihoods outside of the
exponential family, such as ordinal or Cox models (Wurm et al., 2017; Archer
and Williams, 2012; Archer et al., 2014; Fan et al., 2005), but no general
framework exists. Penalised likelihood approaches can be applied to models
outside of the exponential family by approximating the likelihood with an
exponential family distribution. This enables variable selection but the
resulting model is likely to underperform in terms of predictive performance.
A model in the original response space can be obtained by fitting a Bayesian
model including only the selected variables. We compare our method against
baselines that follow this approach in Section 5.2.

Variable selection is also addressed from the Bayesian perspective (O’Hara and
Sillanpää, 2009). This is typically done by imposing so called sparsity
priors, such as the horseshoe (Carvalho et al., 2010; Piironen and Vehtari,
2017b) or the spike-and-slab (Ishwaran and Rao, 2005). As the posterior itself
is not fully sparse, a sparse solution is then obtained by thresholding based
on posterior expectations.

If a Markov chain Monte Carlo (MCMC) algorithm (Robert and Casella, 2013) is
used for inference, no strong assumptions on the likelihood or overall
structure are needed. This makes Bayesian inference applicable to models
outside of the exponential family with multilevel (or other complex)
structure. This helps with the generality of the inference, but still falls
short on the selection, since the user still needs to manually decide which
variables should be selected.

Bayesian reference models have been used for variable and structure selection
tasks in the context of exponential family models, including (additive)
multilevel models (Piironen et al., 2020b; Catalina et al., 2020; Pavone et
al., 2020; Piironen and Vehtari, 2016). However, its application to models
outside of the exponential family has remained unexplored.

Our approach solves variable selection on non exponential family models by
performing the selection on the latent space of the original model, therefore
keeping the original structure. We approximate the unknown latent distribution
with a Gaussian. The problem of learning an implicit distribution has been
approached in the Bayesian inference literature from different angles. Some
authors tackle this problem with variational inference (Blei et al., 2017), by
learning an implicit mapping between samples of the implicit distribution and
a powerful and expressive learner, typically normalising flows (Rezende and
Mohamed, 2016; Huszár, 2017; Pequignot et al., 2020; Titsias and Ruiz, 2019).
Optimising these flexible models is typically expensive, requiring many
iterations and diagnostics to have any guarantee of convergence (Dhaka et al.,
2020). On top of that, the main bottleneck in projection predictive inference
is not solving this projection _once_ , but possibly many times, as complex
models require many posterior draws to be projected for the projections to
fully capture the posterior uncertainty in the reference model.

(a) Residuals histogram plots for three selected projections for a reference
model of size $D=50$. In all cases, the residuals distribution is close to
normal with a decreasing standard deviation as more terms are introduced in
the projection. Notice the shrinking $x$-axis as the projections include more
terms.

subfigure]fig:residuals_combined

(b) KL divergence between reference and projection predictive distributions.
Different linetypes indicate size of the reference model. As projections
include a higher proportion of the reference model, the KL-divergence
approaches 0 faster.

subfigure]fig:kl_combined

Figure 2: Residuals histogram and KL-divergence plots for projections on a
simulated time-to-event survival analysis model and an ordinal cumulative
model, both with $N=100$ observations, and correlation factor $\rho=0.3$.

##  5 Experiments

We evaluate the performance of the proposed method with simulated and real-
world data experiments. In particular, we must account for complex high
dimensional posterior geometries, the role of correlated terms, and the size
of the selected subset of variables in different models outside the
exponential family.

First, we demonstrate the use of different diagnostics for the latent
approximation. Second, we demonstrate the variable selection itself in high
dimensional non-exponential family examples. Then, we extend the experiments
to real datasets. Finally, we show the benefits that our method brings to
models in the exponential family, too.

For the implementation of our new approach, we use modified projpred (Piironen
et al., 2020a). The implementation of the latent projection predictive
inference approach is publicly available at https://github.com/stan-
dev/projpred/tree/latent_projection.

###  5.1 Diagnosing the quality of the latent approximation

We use two diagnostics to asses the quality of the approximate projections:

  * •

We Perform projection predictive checks on the residuals
$\tilde{\eta}_{*}-\tilde{\eta}_{\bot}$, where $\tilde{\eta}_{*}$ corresponds
to the latent predictions of the reference model and $\tilde{\eta}_{\bot}$ to
the predictions of the projection.

  * •

We check the Kullback-Leibler divergence between the reference’s and
projection’s predictions after convergence, which should approach $0$ as more
terms are included in the projection.

We use simulated data from 1) an ordinal cumulative model with probit link
function and 2) a time-to-event survival model with a log link function. The
outcome in both cases is generated as a function of $D\in[50,100,250,500]$
sampled predictors with a uniform correlation factor of $0.3$, where only
$60$% of the predictors have a non-zero effect on the response. The full
generation process for these data is detailed and further analysed in Section
5.2.

##### Projection predictive checks assess normality assumptions.

Histograms of projection residuals (LABEL:fig:residuals_combined) for various
model sizes are a practical diagnostic for the normality assumption. As more
terms enter the projection, the residuals get smaller and more concentrated,
as indicated by the shrinking $x$-axis in the figure.

##### Kullback-Leibler divergence shows that the latent projections eventually
match the reference model predictions.

Even for the most challenging scenarios, the KL-divergence of the latent
predictions shows that the projection predictive distribution gets closer (the
KL-divergence approaches $0$) to the reference predictive distribution as more
terms enter the projection (LABEL:fig:kl_combined).

###  5.2 Non-exponential family models with simulated data

We compare the predictive performance of the optimal submodels in terms of
held-out expected log predictive density (ELPD) (Vehtari and Ojanen, 2012) for
two types of models: an ordinal cumulative model and a time-to-event survival
analysis model with a Weibull hazard process. Additionally, we examine the
performance regarding the selection of truly relevant variables.

For the simulated high-dimensional data, we compare the performance of our
approach to other popular sparsifying solutions in the literature:

  * •

Elastic net regularization as implemented by glmnet (Friedman et al., 2010b),
abbreviated as glmnet in the figures,

  * •

Spike-and-slab sparsifying priors as implemented by spikeSlabGAM (Scheipl,
2011), abbreviated as ss in the figures,

  * •

Spike-and-slab LASSO priors as implemented by SSLASSO (Rockova and George,
2018), abbreviated as sslasso,

  * •

Projection predictive inference on the approximate response space as
implemented by projpred (Piironen et al., 2020a), abbreviated as projpred in
the figures.

Since these approaches do not exist for the specific likelihoods we use, or
for models outside the exponential family in general, we run variable
selection on an approximate model, where we assume a normal likelihood of the
response rather than the appropriate likelihood (cumulative or Weibull in our
examples). Note that this normal approximation is different from our latent
approach, since the latter approximates the latent predictor, not the
response.

Figure 3: Median and 90% confidence interval ELPD difference to the reference
model for the selected subset of predictors for an ordinal cumulative model
(left) and a time-to-event survival model (right) with $N=100,\rho=0$. Higher
values are better. Competing methods are described in Section 5.2.  Figure 4:
Median and 90% confidence interval of the area under the curve (AUC) for the
ROC for the selection of truly relevant terms in an ordinal cumulative model
(left) and a time-to-event survival model (right) with $N=100,\rho=0$ for all
competing methods. Higher values are better. Performance of pure random
selection is shown as dashed dark green line. Competing methods are described
in Section 5.2.

We first fit ordinal and time-to-event survival regression models on extensive
simulation conditions. The generative process is

$\displaystyle x_{n}$ $\displaystyle\sim\mathrm{Normal}(0,\Sigma_{\rho}),\quad
z_{d}\sim\mathrm{Bernoulli}(0.6),$ $\displaystyle\beta_{d}$ $\displaystyle\sim
z_{d}\cdot\mathrm{Normal}(0,1.5),$ $\displaystyle\eta_{n}$ $\displaystyle\mid
x_{n},\beta=\beta^{T}x_{n},\quad y_{n}\mid\eta_{n}\sim p(g(\eta_{n}),\phi),$

where $y\in\mathbb{R}$, $x_{n}\in\mathrm{R}^{D}$ are the outcome and
predictors respectively, $\beta$ denotes the unknown coefficients, number of
variables $D$ is varied in $[50,100,250,500]$ and the number of observations
$N$ is also varied in $[100,200,300]$. $\rho$ indicates the uniform
correlation between predictors, $g$ is the model-specific inverse link
function and $p$ is either an ordinal cumulative likelihood (with $5$ ordered
categories) or a time-to-event likelihood with Weibull base hazard process
(maximum time-to-event is capped to $5$), and $\phi$ are model specific
parameters. Here, we show results for $N=100$ observations and $\rho=0$. We
show more simulations in the Appendix.

We fit the reference model making use of brms (Bürkner, 2018) with a
regularized horseshoe prior (Piironen and Vehtari, 2017b) accounting for the
sparsity in the predictors. The regularized horseshoe prior helps in avoiding
overfitting, particularly for the higher dimensionality models.

Both the KL and residuals diagnostics indicate a good approximation for these
models as shown in Section 5.1. Here we analyse the predictive and selection
performance.

##### Latent projections achieve the best held-out ELPD performance.

The latent approach always achieves the best performance (Fig. 3), even
surpassing the full reference model with models that include substantially
fewer terms. This can be explained by a slight overfitting in the reference
model. In the most complex scenarios, we see that the performance of all
methods suffers from slightly larger variance, including the full reference
model.

##### Approximate likelihood methods result in underperforming selections
(Fig. 3).

While the performance of glmnet is the closest to our proposed method on
average, its performance distribution is wider and less reliable, while other
methods suffer greatly from the approximate likelihood representation. Spike-
and-slab priors (ss in the figures) impose a very strong penalisation, which,
together with the arbitrary choice of threshold ($0.5$ in our experiments, as
is common, e.g. in Scheipl (2011)), can result in suboptimal selections. The
ill-informed approximate likelihood results in very few terms crossing the
selection threshold, particularly in ss and naive projpred. Even though
sslasso does not need an arbitrary threshold for the selection, it suffers
from similar problems.

##### Other approaches require refitting a Bayesian model.

Other approaches operate on an approximate likelihood model that does not
result in a sensible predictive model. To have a functional reduced model in
the same domain as the reference model, all approaches except ours require
refitting a full Bayesian model with the selected predictors and correct (non-
normal) likelihood, increasing their cost. Our latent approach fully honours
the reference model structure, and predictions in the original model space
only require evaluating the inverse link function for the latent predictions.

##### Latent approach retains reference model uncertainty.

The uncertainty present in the reference model is projected as well,
contributing to better informed solutions than the other approaches, which
cannot represent the original model likelihood. As seen in Fig. 4, both
projpred approaches achieve the highest AUC, even though naive projpred fails
to assess the predictive performance of the projections. For the simpler
ordinal experiments, glmnet achieves similar AUC. Interestingly, ss has
trouble identifying any relevant features, and therefore its selection is not
accurate, whereas sslasso manages to identify at least some.

###  5.3 Non-exponential family models with real data

In this section, we assess the selection and predictive capabilities of our
latent approach on real datasets and compare it against glmnet as the best
competitor in our simulated examples. We focus on the quality of the selection
path and the properties of the latent projections. A brief description of each
dataset can be found in the Appendix.

(a) Eyedisease

(b) HCCFrame

(c) Cancer

(d) Rotterdam

Figure 5: Full data zoomed-in ELPD performance comparison in the response
model space. For latent_projpred the projections are computed in the latent
space and the predictions are then transformed. Red shows glmnet and blue
latent_projpred. Error bars indicate ELPD 95% quantiles. Larger error bars
correspond to datasets with fewer observations.

(a) Eyedisease

(b) HCCFrame

(c) Cancer

(d) Rotterdam

Figure 6: Full data ELPD performance in the latent predictor space for
latent_projpred. Reference performance in dark red dashed horizontal line.

##### Latent projections identify superior solution path.

We show the comparison of the latent approach against glmnet in Fig. 5. From a
predictive performance point of view, the latent approach identifies a better
solution path for all datasets we tested, except for eyedisease, where both
approaches are equal. Note that, along the path, there are multiple model
sizes for which the latent approach’s projections’ performance is superior to
their glmnet counterpart. Datasets with few observations (cancer and rotterdam
particularly) suffer from higher variance in performance.

##### The response space is noisier.

Our latent approach allows us to perform the variable selection on either the
response scale or the latent predictor scale. When compared to the selection
in response scale (as seen in Fig. 5), the latent predictor space is
significantly less noisy (Fig. 6), mostly due to the complex observation model
and sometimes relatively small datasets. The reference model filters noise
from the data and thus the latent variables have less noise, resulting in
reduced variance in the selection criterion (Fig. 6 _Latent_ label), even in
data with few observations (Cancer and Rotterdam in particular).

Figure 7: Median and $90$% confidence interval relative wall time of
latent_projpred over projpred for $50$ Bernoulli data realizations and varying
number of observations (different shapes). Lower than $1$ means
latent_projpred is faster. Figure 8: Bootstrap inclusion frequencies for a
simulated Bernoulli dataset with $N=100$ observations and $D=50$ variables
computed from $50$ boostrap samples.

###  5.4 Non-Gaussian exponential family models with simulated data

We simulate $50$ realizations of a Bernoulli data with $N=300$ observations
and $D=50,100,250,500$ uncorrelated variables.

##### The latent approach is faster to compute.

The latent normal approximation is solved efficiently in closed form, which is
faster to compute than the iterative solution of the original projection
approach, resulting in significantly faster solutions (Fig. 7).

##### Latent approach has less variability and selects fewer variables.

We simulate $50$ bootstrap samples from the same Bernoulli data generation
process with $N=100$ observations and $D=50$ variables and compare the
selection of the original framework against the latent approach. As seen in
Fig. 8, the latent approach results in smaller subsets of relevant variables
with a smaller variability across bootstrap samples, while the original
framework fails at identifying a suitable subset and overselects variables in
every case.

Figure 9: Comparison of full data expected log predictive density (ELPD)
performance for a hierarchical Bernoulli model of fertility among kākāpō
population. Reference performance in dark red dashed horizontal line. The
original projpred solution does not reach the reference model predictive
performance due to convergence issues in the underlying projections when
including the varying intercept per individual.

###  5.5 Non-Gaussian exponential family models with real data

For complex non-Gaussian models, the original framework provides a solution
that might suffer from unstable projections and often fail to converge to
sensible results. We show that the latent approach obtains a superior solution
for a real model in the context of endangered species conservation taken from
an ongoing collaboration with domain experts (Digby et al., 2021).

These data consist of a set of 211 kākāpō individuals (it is not possible to
obtain more observations due to the very nature of the data set) and the aim
of the study is to find a model for the fertility of upcoming eggs. With the
domain experts, we constructed a hierarchical Bernoulli model (the response is
either 0 or 1 for fertile and infertile eggs) with a varying intercept per
individual. Such group effects can soak most of the outcome variance in the
model, which paired with the nonlinear link function, can often result in
PIRLS solver convergence issues.

##### Latent approach provides superior performing projections (Fig. 9).

Since the latent predictor is naturally continuous and unbounded, the latent
approach results in more stable and efficient projections for non-Gaussian
exponential family models, as seen in Fig. 9. The latent approach accurately
captures the reference model predictive performance, and also a better
ordering in the solution path.

##  6 Discussion

We proposed a novel latent space projection predictive approach to perform
variable and structure selection in models with non-exponential family
observation model. As shown in the experiments, the approach is also
beneficial for non-Gaussian exponential family models.

We have shown that our approach offers superior performance in extensive
experiments on high dimensional generalized linear models and several real
datasets. The proposed method is not limited to latent linear models, as it
can readily be used for variable and structure selection in more complicated
models, including linear and non-linear hierarchical models (as demonstrated
in Fig. 9).

Our approach not only improves variable and structure selection in a whole new
set of models while respecting their original structure, but also offers
superior performance in currently supported models in projpred.

The method may have worse performance for models whose latent space is far
from normal. We provide diagnostics that alert the user in these cases, but
further research is needed to improve the performance in such cases.

#### Acknowledgments

The authors thank Aalto Science-IT for computational resources and FCAI for
funding and support.

## References

  * Archer and Williams [2012] K. J. Archer and A. A. Williams.  L1 penalized continuation ratio models for ordinal response prediction using high-dimensional datasets.  _Statistics in Medicine_ , 31:1464–1474, 2012.  URL https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3718008/.
  * Archer et al. [2014] K. J. Archer, J. Hou, Q. Zhou, K. Ferber, J. G. Layne, and A. E. Gentry.  ordinalgmifs: An R package for ordinal regression in high-dimensional data settings.  _Cancer informatics_ , 13:187—195, 2014.  ISSN 1176-9351.  doi: 10.4137/cin.s20806.  URL https://europepmc.org/articles/PMC4266195.
  * Barron [1992] D. N. Barron.  The analysis of count data: Overdispersion and autocorrelation.  _Sociological methodology_ , pages 179–220, 1992.
  * Blei et al. [2017] D. M. Blei, A. Kucukelbir, and J. D. McAuliffe.  Variational inference: A review for statisticians.  _Journal of the American statistical Association_ , 112(518):859–877, 2017.
  * Breiman [1995] L. Breiman.  Better subset regression using the nonnegative garrote.  _Technometrics_ , 37(4):373–384, 1995.  ISSN 00401706.  URL http://www.jstor.org/stable/1269730.
  * Bürkner [2018] P.-C. Bürkner.  Advanced Bayesian multilevel modeling with the r package brms.  _The R Journal_ , 10(1):395, 2018.  ISSN 2073-4859.  doi: 10.32614/rj-2018-017.  URL http://dx.doi.org/10.32614/rj-2018-017.
  * Bürkner and Vuorre [2019] P.-C. Bürkner and M. Vuorre.  Ordinal regression models in psychology: A tutorial.  _Advances in Methods and Practices in Psychological Science_ , 2(1):77–101, 2019.
  * Candes and Tao [2007] E. Candes and T. Tao.  The dantzig selector: Statistical estimation when p is much larger than n.  _The Annals of Statistics_ , 35(6):2313–2351, Dec 2007.  ISSN 0090-5364.  doi: 10.1214/009053606000001523.  URL http://dx.doi.org/10.1214/009053606000001523.
  * Carvalho et al. [2010] C. M. Carvalho, N. G. Polson, and J. G. Scott.  The horseshoe estimator for sparse signals.  _Biometrika_ , 97(2):465–480, Apr 2010.  ISSN 1464-3510.  doi: 10.1093/biomet/asq017.  URL http://dx.doi.org/10.1093/biomet/asq017.
  * Catalina et al. [2020] A. Catalina, P.-C. Bürkner, and A. Vehtari.  Projection predictive inference for generalized linear and additive multilevel models.  _arxiv preprint:2010.06994_ , 2020.
  * Dhaka et al. [2020] A. K. Dhaka, A. Catalina, M. R. Andersen, M. Magnusson, J. H. Huggins, and A. Vehtari.  Robust, accurate stochastic optimization for variational inference.  In _NeurIPS_ , 2020.  URL https://proceedings.neurips.cc/paper/2020/hash/7cac11e2f46ed46c339ec3d569853759-Abstract.html.
  * Digby et al. [2021] A. Digby, D. Eason, A. Catalina, K. R. Programme, M. Lierz, S. Galla, G. A. authors TBC, T. Steeves, and D. Vercoe.  Captive rearing reduces fertility in the critically endangered kakapo, 2021.
  * Dupuis and Robert [2003] J. A. Dupuis and C. P. Robert.  Variable selection in qualitative models via an entropic explanatory power.  _Journal of Statistical Planning and Inference_ , 111, 2003.  ISSN 03783758.  doi: 10.1016/S0378-3758(02)00286-0.
  * Fan and Li [2001] J. Fan and R. Li.  Variable selection via nonconcave penalized likelihood and its oracle properties.  _Journal of the American Statistical Association_ , 96(456):1348–1360, 2001.  doi: 10.1198/016214501753382273.  URL https://doi.org/10.1198/016214501753382273.
  * Fan et al. [2005] J. Fan, G. Li, and R. Li.  An overview on variable selection for survival analysis.  In _Contemporary Multivariate Analysis And Design Of Experiments: In Celebration of Professor Kai-Tai Fang’s 65th Birthday_ , pages 315–336. World Scientific, 2005.
  * Friedman et al. [2010a] J. Friedman, T. Hastie, and R. Tibshirani.  Regularization paths for generalized linear models via coordinate descent.  _Journal of Statistical Software_ , 33(1), 2010a.  ISSN 1548-7660.  doi: 10.18637/jss.v033.i01.  URL http://dx.doi.org/10.18637/jss.v033.i01.
  * Friedman et al. [2010b] J. Friedman, T. Hastie, and R. Tibshirani.  Regularization paths for generalized linear models via coordinate descent.  _Journal of Statistical Software_ , 33(1):1–22, 2010b.  URL http://www.jstatsoft.org/v33/i01/.
  * Gelman et al. [2020] A. Gelman, A. Vehtari, D. Simpson, C. C. Margossian, B. Carpenter, Y. Yao, L. Kennedy, J. Gabry, P.-C. Bürkner, and M. Modrák.  Bayesian workflow, 2020.
  * Goutis and Robert [1998] C. Goutis and C. P. Robert.  Model choice in generalised linear models: A Bayesian approach via Kullback-Leibler projections.  _Biometrika_ , 85(1):29–37, Mar 1998.  ISSN 1464-3510.  doi: 10.1093/biomet/85.1.29.  URL http://dx.doi.org/10.1093/biomet/85.1.29.
  * Hastie [2015] T. Hastie.  _Statistical Learning with Sparsity_.  Chapman and Hall/CRC, May 2015.  ISBN 9781498712170.  doi: 10.1201/b18401.  URL http://dx.doi.org/10.1201/b18401.
  * Huszár [2017] F. Huszár.  Variational inference using implicit distributions.  _CoRR_ , abs/1702.08235, 2017.  URL http://arxiv.org/abs/1702.08235.
  * Ishwaran and Rao [2005] H. Ishwaran and J. S. Rao.  Spike and slab variable selection: Frequentist and Bayesian strategies.  _The Annals of Statistics_ , 33(2):730–773, Apr 2005.  ISSN 0090-5364.  doi: 10.1214/009053604000001147.  URL http://dx.doi.org/10.1214/009053604000001147.
  * Kelter [2020] R. Kelter.  Bayesian survival analysis in stan for improved measuring of uncertainty in parameter estimates.  _Measurement: Interdisciplinary Research and Perspectives_ , 18(2):101–109, 2020.  doi: 10.1080/15366367.2019.1689761.  URL https://doi.org/10.1080/15366367.2019.1689761.
  * Koopman [1936] B. O. Koopman.  On distributions admitting a sufficient statistic.  _Transactions of the American Mathematical Society_ , 39(3):399–399, Mar 1936.  ISSN 0002-9947.  doi: 10.1090/s0002-9947-1936-1501854-3.  URL http://dx.doi.org/10.1090/s0002-9947-1936-1501854-3.
  * Marx [1996] B. D. Marx.  Iteratively reweighted partial least squares estimation for generalized linear regression.  _Technometrics_ , 38(4):374–381, 1996.
  * McCullagh and Nelder [1989] P. McCullagh and J. A. Nelder.  _Generalized Linear Models_.  Springer US, 1989.  ISBN 9781489932426.  doi: 10.1007/978-1-4899-3242-6.  URL http://dx.doi.org/10.1007/978-1-4899-3242-6.
  * Nagler [1994] J. Nagler.  Scobit: An alternative estimator to logit and probit.  _American Journal of Political Science_ , 38(1):230–255, 1994.  ISSN 00925853, 15405907.  URL http://www.jstor.org/stable/2111343.
  * O’Hara and Sillanpää [2009] R. B. O’Hara and M. J. Sillanpää.  A review of Bayesian variable selection methods: what, how and which.  _Bayesian Analysis_ , 4(1):85 – 117, 2009.  doi: 10.1214/09-BA403.  URL https://doi.org/10.1214/09-BA403.
  * Pavone et al. [2020] F. Pavone, J. Piironen, P.-C. Bürkner, and A. Vehtari.  Using reference models in variable selection.  _arXiv preprint arXiv:2004.13118_ , 2020.
  * Pequignot et al. [2020] Y. Pequignot, M. Alain, P. Dallaire, A. Yeganehparast, P. Germain, J. Desharnais, and F. Laviolette.  Implicit variational inference: the parameter and the predictor space, 2020.
  * Piironen and Vehtari [2016] J. Piironen and A. Vehtari.  Projection predictive model selection for Gaussian processes.  In _2016 IEEE 26th International Workshop on Machine Learning for Signal Processing (MLSP)_ , pages 1–6, 2016.  doi: 10.1109/MLSP.2016.7738829.
  * Piironen and Vehtari [2017a] J. Piironen and A. Vehtari.  Comparison of Bayesian predictive methods for model selection.  _Statistics and Computing_ , 27(3):711–735, 2017a.  ISSN 0960-3174, 1573-1375.  doi: 10.1007/s11222-016-9649-y.  URL http://link.springer.com/10.1007/s11222-016-9649-y.
  * Piironen and Vehtari [2017b] J. Piironen and A. Vehtari.  Sparsity information and regularization in the horseshoe and other shrinkage priors.  _Electronic Journal of Statistics_ , 11(2):5018–5051, 2017b.  ISSN 1935-7524.  doi: 10.1214/17-EJS1337SI.  URL https://projecteuclid.org/euclid.ejs/1513306866.
  * Piironen et al. [2020a] J. Piironen, M. Paasiniemi, A. Catalina, and A. Vehtari.  _projpred: Projection Predictive Feature Selection_ , 2020a.  https://mc-stan.org/projpred, https://discourse.mc-stan.org/.
  * Piironen et al. [2020b] J. Piironen, M. Paasiniemi, and A. Vehtari.  Projective inference in high-dimensional problems: Prediction and feature selection.  _Electronic Journal of Statistics_ , 14(1):2155–2197, 2020b.  doi: 10.1214/20-EJS1711.  URL https://doi.org/10.1214/20-EJS1711.
  * Rezende and Mohamed [2016] D. J. Rezende and S. Mohamed.  Variational inference with normalizing flows, 2016.
  * Robert and Casella [2013] C. Robert and G. Casella.  _Monte Carlo statistical methods_.  Springer Science & Business Media, 2013.
  * Rockova and George [2018] V. Rockova and E. George.  The spike-and-slab lasso.  _Journal of the American Statistical Association_ , 113(521):431–444, 2018.
  * Scheipl [2011] F. Scheipl.  spikeSlabGAM: Bayesian variable selection, model choice and regularization for generalized additive mixed models in R.  _Journal of Statistical Software_ , 43(14):1–24, 2011.  URL http://www.jstatsoft.org/v43/i14/.
  * Tibshirani [1996] R. Tibshirani.  Regression shrinkage and selection via the Lasso.  _Journal of the Royal Statistical Society. Series B (Methodological)_ , 58(1):267–288, 1996.  ISSN 00359246.  URL http://www.jstor.org/stable/2346178.
  * Titsias and Ruiz [2019] M. K. Titsias and F. Ruiz.  Unbiased implicit variational inference.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , pages 167–176. PMLR, 2019.
  * Vehtari and Ojanen [2012] A. Vehtari and J. Ojanen.  A survey of Bayesian predictive methods for model assessment, selection and comparison.  _Statistics Surveys_ , 6(none):142 – 228, 2012\.  doi: 10.1214/12-SS102.  URL https://doi.org/10.1214/12-SS102.
  * Wurm et al. [2017] M. J. Wurm, P. J. Rathouz, and B. M. Hanlon.  Regularized Ordinal Regression and the ordinalNet R Package, 2017.
  * Zou and Hastie [2005] H. Zou and T. Hastie.  Regularization and variable selection via the elastic net.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 67(2):301–320, Apr 2005.  ISSN 1467-9868.  doi: 10.1111/j.1467-9868.2005.00503.x.  URL http://dx.doi.org/10.1111/j.1467-9868.2005.00503.x.