# Enhanced Winter Road Surface Condition Monitoring with Computer Vision

Risto Ojala, Alvari Seppänen Related code is available at
https://github.com/ojalar/SIWNet

###### Abstract

Winter conditions pose several challenges for automated driving applications.
A key challenge during winter is accurate assessment of road surface
condition, as its impact on friction is a critical parameter for safely and
reliably controlling a vehicle. This paper proposes a deep learning regression
model, SIWNet, capable of estimating road surface friction properties from
camera images. SIWNet extends state of the art by including an uncertainty
estimation mechanism in the architecture. This is achieved by including an
additional head in the network, which estimates a prediction interval. The
prediction interval head is trained with a maximum likelihood loss function.
The model was trained and tested with the SeeingThroughFog dataset, which
features corresponding road friction sensor readings and images from an
instrumented vehicle. Acquired results highlight the functionality of the
prediction interval estimation of SIWNet, while the network also achieved
similar point estimate accuracy as the previous state of the art. Furthermore,
the SIWNet architecture is several times more lightweight than the previously
applied state-of-the-art model, resulting in more practical and efficient
deployment.

###### Index Terms:

Computer vision, convolutional neural networks, intelligent vehicles, vehicle
safety

##  I Introduction

Friction between the road and vehicle tyres plays a key role in defining how a
vehicle should be controlled and manoeuvred in winter conditions. The impacts
of friction on driving safety are substantial, directly affecting factors such
as braking distance and slip angle. The friction between a vehicle tyre and
the road depends on both the tyre and road surface properties. Whereas the
tyre properties remain mostly static, the road surface properties can greatly
vary, especially during winter in countries with below-freezing temperatures.
It has been explicitly noted in the literature that accident rates are
strongly affected by road surface condition [1]. In winter conditions, the
road surface friction properties are mainly dependent on the amount of snow,
ice, and water on the road. Quantifying these road friction properties is
essential for modelling, estimating, and predicting the friction between the
tyre and the road. Consequently, vehicle control needs to adapt to the
prevailing road surface conditions to ensure safe operation. As the
availability and development of automated driving features is on the rise, the
estimation methods for road surface friction properties have an increasingly
important role. Automated driving solutions must be capable of adapting to
different friction conditions, modifying their control based on the
environment.

Road surface condition can be analysed with several on-board methods [2, 3]. A
common approach has been to utilise wheel dynamics information for the
estimation task. However, these methods generally have difficulties
quantifying the road surface friction properties before severe slippage has
occurred. Special optical sensors have also been developed for analysing road
surface condition, yet these are generally too expensive for consumer vehicle
applications. Recently, computer vision solutions with typical visible light
cameras have been a popular approach for the task. Computer vision approaches
offer the convenience of utilising the existing windshield camera equipment.
Additionally, they have the potential of assessing the road surface condition
in front of the vehicle, which enables predictive actions.

This paper expands computer vision-based road surface condition monitoring in
winter conditions by presenting a deep learning model, SIWNet, for the task.
The regression model is developed to predict a scalar estimate for the road
surface friction properties, summarising the effects of visible snow, ice, and
water on tyre-road friction. Hence the name of the model, SIWNet (Snow, Ice,
Water Network). Additionally, the model estimates the uncertainty of the
prediction by providing a prediction interval with a multi-head architecture.
Such feature has not been previously proposed for regression-based road
surface condition monitoring models. Furthermore, SIWNet has been designed
with practical on-board deployment in mind, the model being several times
smaller than the previous state of the art [4], yet achieving similar
accuracy. Similarly to [4], SIWNet is trained in an automated manner by
matching images and corresponding friction information. However, the work
presented here utilises the SeeingThroughFog-dataset [5], where the friction
values have been acquired with an optical road friction sensor.

The novel contributions of this work can be summarised as:

$\bullet$

SIWNet is the first road surface friction regression model to feature
prediction intervals in the estimates.

$\bullet$

SIWNet is more computationally lightweight than previous state of the art.

$\bullet$

This paper is the first work training a road surface friction estimation
computer vision model for winter conditions based on optical road friction
sensor data.

##  II Related Work

###  II-A Road surface condition monitoring

Road surface condition monitoring is essential for different vehicle
applications, and the topic has been studied in the past with several
approaches. Current interest in automated driving solutions has increased the
importance of the field, as vehicle control algorithms require information of
road surface condition to properly assess the situation. Road surface
condition monitoring is especially important in winter conditions, where the
road surface friction properties greatly vary. The existing methods can be
divided to contact and non-contact based methods.

Contact methods are not analysed here in-detail, as this work is focused on
computer vision methodology. For a more thorough overview of contact-based
methods, the review by Acosta et al. [2] is recommended. Fundamentals of
contact-based estimation methods are still briefly presented here. Contact-
based methods generally measure the actual friction between the tyre and the
road. Road surface condition can be consequently estimated from this
information.

Of the contact-based methods, slip-based methods are the most common.
Measuring wheel rotation information and inertial measurement unit readings,
friction between the tyre and the road can be estimated [6]. Utilising this
information, the road surface condition can be defined based on the friction
characteristics. In a commercial vehicle, this information is readily
available in the anti-lock braking unit and electronic stability control.
However, slip-based methods typically cannot accurately estimate the road
surface conditions until severe slippage has already occured. The review by
Acosta et al. [2] notes that slip-based approaches are typically considered
inadequate for reliably improving the functionality of advanced driver
assistance systems. Other contact-based methods for road surface monitoring
are based on vibration. Low frequency methods can utilise signals such as the
vehicle rotation speed [7]. However, the approaches rely largely on slip-slope
assumptions, and consequently lack robustness. High frequency methods have
yielded impressive results, yet these methods rely on additional sensor
equipment. Works on the topic have utilised microphones installed on the
vehicle to monitor the tyres [8]. This might not be a commercially feasible
solution due to the burden of additional sensor installations.

To enable alternative means for estimating road surface condition, a number of
non-contact methods have been developed. Non-contact methods typically exhibit
different operation characteristics and attributes compared to the contact-
based methods. They offer alternative ways for the measurement, or they could
be fused with contact-based methods. Review for non-contact friction
estimation presented by Ma et al. [3]. The non-contact approaches can be
roughly divided into methods that utilise special optical sensors, and methods
which utilise computer vision algorithms to process images captured with
traditional visible light cameras. Additionally, there have also been some
studies exploring road surface condition monitoring with automotive radars
[9].

With dedicated optical sensors, non-contact approaches commonly utilise
infrared spectroscopy [10]. The approach is based on the different reflectance
characteristics of water, ice, and snow. This measurement approach is commonly
utilised, and multiple commercial products applying the method are available
on the market [11, 12, 13]. Another optical technique for road surface
condition monitoring is based on analysing polarisation of light [14].

Computer vision applied on regular visible light cameras has been a popular
topic for non-contact road surface condition estimation. The approach is
lucrative from a practical point-of-view, since modern vehicles are equipped
with forward-facing cameras. Review of computer vision-based estimation has
been prepared by Wu et al. [15]. Generally, computer vision approaches utilise
machine learning models such as convolutional neural networks (CNNs) to
analyse images of the road.

Several works utilising computer vision have applied classification techniques
to assess the road surface condition. Nolte et al. [16] proposed applying CNNs
to perform this classification task. They trained ResNet50 [17] and
InceptionV3 [18] models to recognise six distinct categories of road surface:
asphalt, dirt, grass, wet asphalt, cobblestone, and snow. Similar studies have
been conducted by Šabanovič et al. [19], who developed a CNN capable of
classifying the road pavement type (asphalt, cobblestone, gravel) as well as
the surface condition (dry, wet). The effectiveness of the developed road
surface monitoring system was demonstrated with vehicle braking tests. The
tests highlighted that the stopping distance was shorter when utilising an
adaptive anti-lock braking system control strategy, which was tuned based on
the classification results. To extend the capabilities of road surface
classification, Wang et al. [20] proposed applying a segmentation CNN to
perform the classification task. Their classification task featured a total of
nine different pavement and road surface condition combinations, including
categories from winter conditions.

To further enhance the development of computer vision-based road surface
monitoring in summer conditions, Cordes et al. published an open
classification dataset called RoadSaw [21]. They collected the dataset in an
automated manner, mounting an optical road surface monitoring sensor as well
as a camera on a vehicle. As a result, a large dataset of images in realistic
road conditions was captured, along with the corresponding road surface
information. The dataset features three different pavement types (asphalt,
cobblestone, concrete), which have four different surface conditions (dry,
damp, wet, very wet). The authors evaluated the performance of a MobileNetV2
[22] classifier on the data. They also proposed adding deterministic
uncertainty quantification [23] functionality to the network, allowing
assessment of the uncertainty of the classification predictions.

Classification-based road surface monitoring has also been extended to model
the road surface in finer spatial detail. In the camera view, several road
surface conditions may be visible simultaneously, such as patches of snow on
an asphalt road. Providing a single classification for the entire visible road
may result in inadequate assessment of the prevailing situation. Roychowdhury
et al. [24] proposed a two-stage approach for detailed analysis of the road
surface in winter conditions. After classifying the overall road surface
condition with a CNN, they split the image of the road into a grid, in which
each cell was classified separately.

The resolution of road surface condition estimation has been also extended by
applying regression models to the task, instead of utilising classification
models. In terms of friction-related information, regression models allow for
more accurate representation of the road surface properties. This is due the
model predicting a continuous value, instead of utilising a discrete number of
categories labelled with certain values. Vosahlik et al. [4] proposed
automated training of a CNN regression model based on corresponding friction
information derived from a slip-based contact method. Based on data acquired
with a sub-scale car model, they created a dataset of matching friction values
and images, which also included samples from winter conditions. They trained a
ResNet50 network on the data to perform the regression task. Regression-based
road surface condition estimation was also developed by Du et al. [25], who
studied the problem in summer conditions. They collected road surface
condition information with a commercial contact-based monitoring system,
utilising a trailer pulled behind the vehicle. With matching image
information, they trained ResNet50 and InceptionV3 [18] regression models to
predict the road surface friction properties. The vision-based estimators were
fused with another neural network, which processed the vehicle dynamics
information. The fusion model seemed to offer clear benefits, as several key
performance values were improved.

###  II-B Prediction intervals for deep learning regression models

Uncertainty assessment has been a key topic in deep learning due to the
obscurity of the data-driven the models. In order to evaluate the reliability
of the estimates produced by deep learning models, several approaches have
been proposed to quantify related uncertainties. An extensive review on the
topic has been written by Abdar et al. [26]. Most extensive uncertainty
quantification methodologies include Bayesian neural networks and ensemble
techniques. Bayesian neural networks model the network weights as
distributions, allowing the predictions to be accurately modelled as posterior
distributions. Ensemble models rely on utilising multiple neural network
models to process the input, and determining the uncertainty based on the
outputs of the networks. Bayesian and ensemble methods have been proven to
generally provide reliable uncertainty information. However, these methods
typically require immense computational resources to operate, and applying
them is often not practically feasible.

A common approach for assessing uncertainty in regression problems is the
estimation of prediction intervals [27]. Several approaches have been proposed
for generating prediction intervals with deep learning regression models,
including the previously mentioned Bayesian and ensemble methods. However,
more lightweight approaches for estimating prediction intervals have also been
developed. Typically this has been achieved by modifying the neural network
architecture to feature an additional output node. In order to quantify
uncertainty, Nix and Weigend [28] proposed adding a separate output node with
its own hidden layers to a fully connected network architecture. The
additional output node was responsible for estimating the variance of each
prediction on each forward pass, whereas the other output node remained
responsible for producing the point estimate. Modelling each network
prediction as a probability distribution, the network could be trained with a
maximum likelihood loss function. In such architecture, uncertainty is
quantified by the value of the variance output, which can also be used to
generate prediction intervals. Somewhat similar prediction interval estimation
has been proposed by Khosravi et al. [29]. Their approach was based on two
output nodes, which were responsible for predicting the lower and upper bound
of the prediction interval. A special loss function was used for training,
which determined the target coverage of the prediction interval.

###  II-C Research gap

This paper aims to enhance the existing state of the art of winter road
surface condition evaluation. A novel CNN architecture, SIWNet, is proposed
for the task of computer vision-based estimation of road friction properties.
Similarly to the work of Vosahlik et al. [4], SIWNet is implemented as a
regression model.

SIWNet expands state of the art by including an additional prediction head in
the network architecture for uncertainty quantification. This enables
representation of the model output as a prediction interval. The approach is
based on the work of Nix and Weigend [28]. Uncertainty quantification is a
vital feature for computer vision-based road surface monitoring systems, as
visual estimation of road surface friction properties is bound to feature
varying levels of uncertainty. Cordes et al. [21] have previously proposed an
uncertainty estimation approach in their classification-based work, yet this
is the first paper to propose such feature for a regression approach.

In addition to including an uncertainty quantification mechanism, SIWNet is
also computationally extremely lightweight. The findings of this paper
highlight that winter road surface condition estimation does not require an
extensively large neural network architecture. This is a notable benefit,
considering the generally limited vehicle on-board computational resources.

##  III Methods

###  III-A Problem formulation and dataset

The goal of the research was to develop a model capable of assessing the road
surface friction properties based on images captured from a typical vehicle
windshield camera. To acquire relevant data for training such a model, the
publicly available SeeingThroughFog dataset [5] was utilised. The dataset has
been gathered by driving an instrumented vehicle in central European and
Nordic countries during winter. Part of the dataset has also been recorded in
the springtime, to include summer-like conditions. The instrumented vehicle
contained a forward-facing camera, as well as an optical road friction sensor.
The sensor setup is illustrated in Figure 1. In this paper, the camera images
and corresponding road friction sensor readings were utilised for developing a
neural network model, SIWNet, for assessing the road surface condition. SIWNet
is trained to process an image of the road, and predict a value corresponding
to the road friction sensor reading.

Figure 1: Illustration of the sensor installations, with the camera at the
windshield and the road friction sensor at the bumper.

The optical friction measurement unit used in the dataset is a prototype from
a widely recognised manufacturer, seemingly similar to those of the
manufacturer’s other models [11, 12, 13]. Similarly to the commercial
equivalents of the sensor, the unit measures the amount of snow, ice, and
water on the road surface. By manufacturer software, these road surface
friction properties are summarised into a single factor, grip factor, within
range 0.09…0.82. This factor effectively indicates the slipperiness of the
road. Here, this grip factor was normalised to a range 0.00…1.00, and called
friction factor, $\hat{f}$. It should be noted that the friction factor is not
the friction coefficient between the tyre and the road, as determining this
value depends on the tyre properties. The friction factor represents a value
that can be used to estimate slipperiness of the road, and consequently the
actual friction coefficient in case relevant tyre parameters are known.

The images in the dataset have a resolution of 1920x1024 pixels. To focus the
analysis on the most relevant portion of the image, the road directly ahead, a
predefined static section of the images was cropped and transformed to
bird’s-eye-view. The predefined section was chosen manually based on a single
image, and the bird’s-eye-view transformation was achieved by stretching the
images to a square. This process is visualised in Figure 2. These pre-
processing steps are similar to those presented in [4].

Figure 2: Transformation of the area in front to bird’s-eye-view.

The used dataset contained a total of 4330 samples, as this is the number of
measurements in the SeeingThroughFog-dataset which contain readings from the
road friction sensor. These samples were collected in different locations on
12 different days. Figure 3 displays the number of samples collected on each
date. The distribution of the friction factor values in the dataset is
visualised in Figure 4. Image samples from the dataset with corresponding
friction factor readings are provided in Figure 5.

Figure 3: Number of samples per date in the dataset.
---
Figure 4: Friction factor values in the dataset, with a zoomed view providing a clearer depiction of the under-represented values. |  |  |
---|---|---|---
$f=0.19$ | $f=0.38$ | $f=1.0$ | $f=0.67$
|  |  |
$f=0.69$ | $f=1.0$ | $f=0.94$ | $f=0.50$
|  |  |
$f=0.38$ | $f=0.98$ | $f=0.81$ | $f=0.39$
Figure 5: Samples from the utilised data, featuring images with the
corresponding ground truth friction factor values.

###  III-B Model architecture

The architecture of SIWNet consists of a feature backbone, as well as a point
estimate head and a prediction interval head. SIWNet architecture is presented
in detail in Figure 6. The feature backbone is responsible for processing
relevant features from the images. Based on this information, the point
estimate head outputs the predicted friction factor $\hat{f}$. The prediction
interval head is responsible for assessing the uncertainty related to the
point estimate. Based on the features and the point estimate, the prediction
interval head outputs a predicted standard deviation $\hat{\sigma}$, which is
utilised to establish a prediction interval around the point estimate.

Figure 6: SIWNet architecture, with tensor sizes reported for processing a
single image of 324x324 pixels. Each block shows the number of features in the
output.

The feature backbone is based on the ResNet [17] architecture, applying the
same basic residual building blocks. Contrary to typical ResNet
implementations, each block is used only once. This design results in an
extremely lightweight architecture, ensuring practical applicability in
embedded on-board systems. In the point estimate head, a single fully
connected layer is applied, which is typical for regression tasks. A sigmoid
activation function was added to the final output, as the friction factor
values were bound between 0 and 1. The prediction interval head is also
constructed of fully connected layers. These layers analyse the features
provided by the feature backbone, as well as the prediction provided by the
point estimate head. At the end of the prediction interval head, a sigmoid
activation was added to enhance the stability of the output. This conveniently
limited the predicted standard deviation $\hat{\sigma}$ to positive numbers.
The upper limit of 1 for $\hat{\sigma}$ was deemed reasonable, since with the
probabilistic model applied to train the network, this corresponded to a
nearly uniform distribution within the range 0…1.

###  III-C Loss function for prediction interval head

The prediction interval head of SIWNet is trained to quantify the uncertainty
related to the prediction of the point estimate head. Based on the quantified
uncertainty, a prediction interval can be generated. In order to perform this
task, the output of SIWNet is modelled as a truncated normal distribution. The
goal of the training process is to maximise the likelihood of the training
labels with regard to the predicted distributions. Previous work has applied
similar methodology with a regular normal distribution [28]. A truncated
normal distribution of a random variable $x$ has a probability density
function (PDF) of the form [30]

$\displaystyle
p(\mu,\sigma,a,b;x)=\frac{\phi(\frac{x-\mu}{\sigma})}{\Phi(\frac{b-\mu}{\sigma})-\Phi(\frac{a-\mu}{\sigma})}$
(1)

where $\mu$ and $\sigma$ denote the mean and standard deviation of the
underlying normal distribution, respectively. The lower and upper truncation
bounds are denoted by $a$ and $b$, respectively. The truncated normal
distribution PDF representation is based on the underlying normal distribution
PDF, defined as

$\displaystyle\phi(\frac{x-\mu}{\sigma})=\frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(x-\mu)^{2}}{2\sigma^{2}}}$
(2)

as well as the underlying normal distribution cumulative distribution
function, defined as

$\displaystyle\Phi(\frac{b-\mu}{\sigma})=\frac{1}{2}(1+erf(\frac{b-\mu}{\sigma\sqrt{2\pi}})),$
(3)

where $erf(\cdot)$ is the error function. Consequently, the negative log-
likelihood of the truncated normal distribution has the form

$-\ln{p(\mu,\sigma,a,b;x)}=\ln\sigma+\frac{(\mu-x)^{2}}{2\sigma^{2}}\\\
+\ln{(erf(\frac{\mu-b}{\sigma\sqrt{2}})-erf(\frac{\mu-a}{\sigma\sqrt{2}}))}.$
(4)

The training process of the prediction interval head of SIWNet is based on
maximising the likelihood of the corresponding ground truths and predicted
truncated normal distributions. Consequently, this means minimising the
negative log-likelihood with the predictions. The utilised loss function for a
training batch has the form

$\displaystyle\mathcal{L}=\sum_{i=i}^{n}-\ln{p(\hat{f}_{i},\hat{\sigma}_{i},a,b;f_{i})}$
(5)

where $n$ denotes the number of samples in the batch and $f$ denotes the
ground truth friction factor. The underlying normal distribution mean is the
predicted friction factor $\hat{f}$ from the point estimate head, whereas the
predicted standard deviation $\hat{\sigma}$ is the output of the prediction
interval head. To improve stability, $\hat{\sigma}$ is thresholded to a
minimum value of $1\cdot 10^{-4}$ when computing the loss.

###  III-D Training and testing

For training SIWNet and conducting the presented experiments, the following
steps were taken. The utilised dataset was randomly divided into train-
validation-test sets with a respective 70%-15%-15% split. When splitting the
data, timestamps were utilised to ensure that data samples gathered from the
same location were not included in different sets. Hyperparameters were
optimised by training the network on the training set and finding the best
result on the validation set. For the evaluation on the test set, the model
was re-trained on a combination set containing both the training and
validation sets.

During the training, the feature backbone and point estimate head were first
trained for 60 epochs utilising regular mean squared error as the loss
function. After these parts of the network were trained, their weights were
frozen. Afterwards, the prediction interval head was trained for 60 epochs
with the loss function presented in Equation 5. Dropout with a probability of
0.5 was applied when feeding the feature backbone output to the prediction
interval head during training.

Training was carried out with a batch size of 32, and stochastic gradient
descent with a momentum of 0.9 was used as the optimisation method. During
training, the learning rate was decayed with a step-based scheduler. Every
twenty epochs the learning rate was reduced to one tenth of the previous
value. As for data augmentation during training, the images were randomly
flipped horizontally as well as rotated with a value from [-4,4] degrees.
Furthermore, the pixel values were slightly scaled by random color jitter in
the range [0.9, 1.1]. For both training and testing, images were reshaped to
324x324 pixels before being fed to the network. Additionally, the pixel values
were normalised with the mean and standard deviation of the pixel values in
the training set.

The presented experiments demonstrate comparison of SIWNet to the ResNet50
model previously applied in the literature [4]. The same training and testing
procedures were applied on the ResNet50 model, except for steps related to the
prediction interval head, which the ResNet50 model does not include. A sigmoid
activation was also added to the output of the ResNet50 model, since this was
noted to boost performance. Both models were implemented on the PyTorch deep
learning framework [31]. Number of trainable parameters in the models as well
as floating point operations executed during inference were analysed with the
ptflops-tool [32].

In order to benchmark the prediction intervals and the probabilistic
properties of the predictions in the experiments, interval score and
continuous ranked probability score (CRPS) were utilised [33]. The 90%
prediction interval was used for defining the interval score. For constructing
the prediction interval based on the truncated normal distribution output of
SIWNet, the 90% interval of the distribution with the highest likelihood was
used. Since the compared ResNet50 model produces only point estimates, there
was no straightforward approach to define interval scores for the model. Thus,
the ResNet50 model was given a static prediction interval surrounding the
point estimate. The boundaries were set at a distance from the point estimate,
which corresponded to the 90% error threshold $e_{90\%}$ on the validation
set. The 90% error was defined as the value below which 90% of absolute errors
on the validation set were located. The prediction interval boundaries were
clamped by the plausible friction factor values, 0 to 1.

##  IV Results

SIWNet was evaluated with several experiments, to assess the efficacy of the
model in the task of predicting road surface friction properties. The
computational demands of SIWNet were also analysed to determine the
applicability of the model in embedded on-board usecases. The accuracy of
SIWNet was evaluated on the test set in terms of point estimates, as well as
prediction intervals. Samples of SIWNet predictions plotted next to
corresponding test set images are presented in Fig. 7. Each prediction plot
features the point estimate of the friction factor $\hat{f}$, as well as the
truncated normal distribution and prediction interval estimated by the
prediction interval head. Results of SIWNet were compared to those of
ResNet50, which has in previous literature been applied for a similar
regression task.

|
---|---
|
|
|
Figure 7: Samples of SIWNet predictions on the test set.

###  IV-A Model size and computational load

The neural network models were evaluated in terms of computational load by
analysing the number of trainable parameters in the architectures, as well as
the floating point operations (FLOPs) required to perform inference on a
single input image. The SIWNet architecture was designed to require minimal
computational resources, as evident in the results presented in Table I.

TABLE I: Model sizes and computational loads. |  |
---|---|---
Method |  Number of parameters ($\cdot 10^{6}$) |  Inference load (GFLOPs)
SIWNet |  5.04 |  4.23
ResNet50 |  23.5 |  18.3

Compared to the previously applied ResNet50 model, SIWNet features several
times fewer parameters. Additionally, the SIWNet model requires notably fewer
FLOPs to process an image. The results highlight that compared to ResNet50,
SIWNet requires less computational resources to operate. This indicates that
SIWNet is well-fitted for on-board utilisation with limited embedded hardware.

###  IV-B Point estimate head performance

Training the feature backbone and point estimate heads were the initial steps
for training SIWNet. Hyperparameters were tuned with grid search, resulting in
the initial learning rate and weight decay values of 0.1 and $1\cdot 10^{-3}$,
respectively. Identical hyperparameter tuning was performed for the ResNet50,
resulting in the same optimal parameter values. Table II reports the point
estimate accuracies achieved by the networks when predicting friction factors
on the test data. Accuracies are reported in terms of mean absolute error
(MAE) and root-mean-square error (RMSE).

TABLE II: Point estimate errors on the test set. |  |
---|---|---
Method |  MAE |  RMSE
SIWNet |  0.078 |  0.132
ResNet50 |  0.079 |  0.132

The presented results highlight that SIWNet achieved a nearly identical point
estimate accuracy as ResNet50. Considering the minimal size of the SIWNet,
this indicates that a larger architecture does not offer added benefit for the
prediction task.

###  IV-C Prediction interval head performance

The key feature of SIWNet is its capability of assessing the uncertainty of
the friction factor prediction, enabled by the prediction interval head. The
prediction interval head was trained while keeping the feature backbone and
point estimate head frozen. Hyperparameters were optimised with grid search,
tuning the initial learning rate, weight decay, as well as number of neurons
in the fully connected layers of the prediction head. Based on the
optimisation, the initial learning rate and weight decay were set at values of
$5\cdot 10^{-4}$ and $1\cdot 10^{-3}$, respectively. Acquired average interval
score results are presented in Table III for SIWNet and ResNet50. As described
previously, the prediction interval for ResNet50 was formulated based on the
90% error on the validation set.

TABLE III: Average interval scores of the networks on the test set. |
---|---
Method |  Average interval score*
SIWNet |  0.482
ResNet50 with $e_{90\%}$ |  0.656
*Lower is better

Based on the comparison, SIWNet clearly achieved a more favourable average
interval score. This indicates that the prediction interval head was capable
of assessing the uncertainty related to the friction factor predictions of the
point estimate head.

As another test to evaluate the performance of the prediction interval head,
the CRPS metric was used to assess the quality of the distributions predicted
by SIWNet. The point estimates of ResNet50 were scored on the metric for
comparison. For point estimates, the CRPS is equal to the MAE score.

TABLE IV: Average CRPS values of the networks on the test set. |
---|---
Method |  Average CRPS*
SIWNet |  0.060
ResNet50 |  0.079
*Lower is better

The acquired results highlight that the probabilistic forecasts of SIWNet were
capable of more accurate prediction than the point estimates of ResNet50. This
indicates that the prediction interval head of SIWNet did learn a functional
strategy for assessing the uncertainty of the friction factor estimates.

###  IV-D Ablation study

In order to ensure that the design decisions behind the SIWNet architecture
were sensible and beneficial for the performance, an ablation study was
carried out. Similarly to the previously presented results, the ablation study
analysed the point estimate head and prediction interval head separately.

The point estimate head was modified by removing the sigmoid activation, and
the related training and testing procedures were repeated. Identical actions
were taken with the ResNet50, removing the sigmoid activation from its output.
During testing, the outputs of the networks were clamped between 0 and 1.
Resulting accuracies on the test set are presented in Table V. Based on the
results, it is evident that the sigmoid activation was a beneficial addition
for both networks.

TABLE V: Ablation study of the point estimate head. |  |
---|---|---
Method |  MAE |  RMSE
SIWNet w/o sigmoid |  0.099 |  0.151
ResNet50 w/o sigmoid |  0.144 |  0.190

Ablation study of the prediction interval head also investigated the effect of
the sigmoid activation. The sigmoid activation was removed from the prediction
interval output, and training and testing were repeated. In another test, the
prediction interval head was trained and tested without using dropout
regularisation on the feature backbone output. Additionally, the SIWNet model
was studied by training the model without the proposed technique of
pretraining the feature and backbone and point estimate head, and freezing
their weights. Instead the whole network was optimised simultaneously
utilising the loss function presented in Equation 5. Finally, the overall
efficacy of the prediction interval head was evaluated by formulating a static
prediction interval around the point estimates based on the 90% error on the
validation set, similarly to the prediction interval used for the ResNet50
predictions. The ablations were evaluated with the interval score, and results
from the tests are presented in Table VI. All tests concluded that the
originally proposed SIWNet architecture was capable of scoring higher interval
scores.

TABLE VI: Ablation study of the prediction interval head. |
---|---
Method |  Average interval score*
SIWNet w/o dropout |  0.621
SIWNet with $e_{90\%}$ |  0.710
SIWNet w/o pretrain & freeze |  0.748
SIWNet w/o PI Sigmoid |  0.904
*Lower is better

##  V Discussion

SIWNet was capable of producing point estimates with equivalent accuracy as
the ResNet50. This is an impressive result, considering that the model has
approximately 79% fewer parameters, and required roughly 77% fewer FLOPs to
process an input image. The presented experiments also highlight that SIWNet
was capable of successfully assessing the uncertainty related to its friction
factor predictions. This was evident in the CPRS as well as interval score
results in Tables IV and III. However, exact quantification of the prediction
interval capabilities of SIWNet is difficult, since only a naive method was
used for comparison. The static prediction interval applied with ResNet50
point estimates was by no means an optimal method for uncertainty
quantification. However, since no previous works applying uncertainty
quantification to road surface friction regression exist, there were no clear
alternatives for comparison. The implementation of SIWNet is published as open
source, and future works on the topic are encouraged to compare their
approaches to SIWNet.

Concerning the training procedure of SIWNet, the chosen approach was likely
not optimal. During training of the prediction interval head with the loss
function presented in Equation 5, the point estimate head as well as the
feature backbone were frozen. This was done to stabilise the learning process,
as the loss function can produce quite extreme values. The result of a direct
training approach was demonstrated in the ablation study. However, the
training stability comes at the cost of not finding optimal weights in the
feature backbone for the uncertainty estimation. Finding a strategy to
optimise the feature backbone accordingly might result in better uncertainty
quantification and prediction interval estimates.

Regarding the validity of the presented results overall, some unreliability
was likely caused by the used dataset and its characteristics. Overall, a
larger dataset would provide more definitive results, as the models could
learn richer representations from data, and the testing would provide a more
thorough investigation of the prediction capabilities. In addition, the used
dataset contained a clear over-representation of certain values, as visible in
Figure 4. A more balanced dataset would likely result in better learning
results of the detectors, as an unbalanced dataset can result in bias in the
models. These factors may have contributed to the SIWNet and ResNet50 models
achieving such similar point estimate accuracies. Increased variation in the
used dataset might have resulted in more notable differences in the prediction
results.

Another minor disadvantage of the used dataset was the fact that the road
friction sensor readings did not exactly correspond to the visible road in the
camera view. This can be seen in Figure 1, which illustrates the sensor
placement. The road friction sensor was measuring the road area directly below
the bumper of the vehicle. This specific road area was not actually visible in
the camera view simultaneously, as the camera was monitoring the road slightly
ahead of the bumper. Consequently, in order to use the data, an assumption had
to made that the road friction properties were the same below the bumper and
in the visible road area. If full time series of the recordings were
available, this issue could be mitigated by matching the sensor reading to a
previously captured image. However, assumption of same road condition below
the bumper and slightly ahead in the camera view should hold most of the time.
Effects of this drawback were likely minimal in the presented results.

Perhaps a more impactful feature of the used dataset was the utilisation of a
single road friction sensor reading for determining the ground truth value for
an entire image. The road friction sensor measured only a single point on the
road, whereas the camera image captured a large portion of the road ahead. In
some scenarios, the visible road conditions may have varied greatly within the
image. For example, the road may have been partially covered in snow, while
also featuring tyre tracks with clear asphalt or water. This naturally led to
a certain degree of ambiguity and noise in the data, as the ground truth value
has not always fully represented the entire information in the image. This was
an inherent problem in the dataset, which likely resulted in the models
learning approximate solutions that on average provided the highest scoring
results despite the noise in the data. However, despite this limitation,
similar datasets with limited ground truth coverage have been previously
successfully applied in road condition monitoring development [4, 21]. An
optimal dataset for road surface friction estimation purposes should include
several ground truth labels for the visible road area, allowing more fine-
grained estimation. This could be achieved by for example mounting a vehicle
with several measurement units. To the best of the authors’ knowledge, no such
public dataset unfortunately exists.

It should also be noted that the utilised dataset, and consequently the
trained models, only consider road friction properties related to snow, ice,
and water. This was due to the focus on winter conditions, in which these
factors largely determine the overall friction properties of the road.
Consequently, analysis of road surface pavement type was ignored, which can be
a drawback in certain road conditions. For better generalisability, future
development of SIWNet could aim to include the road pavement type as a factor
in the analysis.

##  VI Conclusion

This work enhanced camera-based winter road condition monitoring by presenting
the SIWNet model. Based on an image of the road, the deep learning regression
model predicts a friction factor, which summarises the friction properties of
the visible road. SIWNet also includes built-in uncertainty quantification
capabilities, allowing the model to generate prediction intervals instead of
traditional point estimates. Due to this feature, SIWNet enables robust road
condition monitoring via computer vision. This is a key finding, as reliable
road condition monitoring is vital for proper tuning of controllers in
automated vehicle applications. Additionally, the SIWNet model is
computationally lightweight, allowing convenient implementation of the network
in embedded on-board systems.

Future research efforts on the topic should focus on more fine-grained spatial
analysis of winter road surfaces. Computer vision solutions for road condition
monitoring could be extended to pixel-level analysis of the road surface.
Furthermore, spatio-temporal data analysis might also provide clear benefits
to the accuracy at which road condition can be monitored. Integrating such
features to the SIWNet architecture could allow for even more reliable
modelling and estimation of the road surface conditions in the future.
Considering the wide variety of winter road conditions and their effect on
optimal automated driving strategies, additional research efforts should be
concentrated on developing increasingly robust solutions for these demanding
circumstances.

In order to improve the reliability of automated driving solutions, the
technology must be capable of operating in all weather and road conditions.
Proper situational awareness of vehicles in adverse conditions is a key factor
in enabling safe and robust future automated vehicle solutions. SIWNet
advances the road condition monitoring capabilities of vehicles in winter
conditions, potentially enhancing the operation of automated driving
functionalities and bringing reliable automated driving functions closer to
reality.

## Acknowledgments

The authors wish to acknowledge the funding provided by Henry Ford Foundation
Finland.

## References

  * [1] C.-G. Wallman and H. Åström, “Friction measurement methods and the correlation between road friction and traffic safety: A literature review,” 2001\.
  * [2] M. Acosta, S. Kanarachos, and M. Blundell, “Road friction virtual sensing: A review of estimation techniques with emphasis on low excitation approaches,” _Applied Sciences_ , vol. 7, no. 12, p. 1230, 2017.
  * [3] Y. Ma, M. Wang, Q. Feng, Z. He, and M. Tian, “Current non-contact road surface condition detection schemes and technical challenges,” _Sensors_ , vol. 22, no. 24, p. 9583, 2022.
  * [4] D. Vosahlik, J. Cech, T. Hanis, A. Konopisky, T. Rurtle, J. Svancar, and T. Twardzik, “Self-supervised learning of camera-based drivable surface friction,” in _2021 IEEE International Intelligent Transportation Systems Conference (ITSC)_. IEEE, 2021, pp. 2773–2780.
  * [5] M. Bijelic, T. Gruber, F. Mannan, F. Kraus, W. Ritter, K. Dietmayer, and F. Heide, “Seeing through fog without seeing fog: Deep multimodal sensor fusion in unseen adverse weather,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 11 682–11 692.
  * [6] F. Gustafsson, “Slip-based tire-road friction estimation,” _Automatica_ , vol. 33, no. 6, pp. 1087–1099, 1997.
  * [7] D. Pavković, J. Deur, J. Asgari, and D. Hrovat, “Experimental analysis of potentials for tire friction estimation in low-slip operating mode,” _SAE Transactions_ , pp. 369–380, 2006.
  * [8] J. Alonso, J. López, I. Pavón, M. Recuero, C. Asensio, G. Arcas, and A. Bravo, “On-board wet road surface identification using tyre/road noise and support vector machines,” _Applied acoustics_ , vol. 76, pp. 407–415, 2014.
  * [9] V. V. Viikari, T. Varpula, and M. Kantanen, “Road-condition recognition using 24-ghz automotive radar,” _IEEE transactions on intelligent transportation systems_ , vol. 10, no. 4, pp. 639–648, 2009.
  * [10] P. Jonsson and M. Riehm, “Infrared thermometry in winter road maintenance,” _Journal of Atmospheric and Oceanic Technology_ , vol. 29, no. 6, pp. 846–856, 2012.
  * [11] _Remote Road Surface State Sensor DSC111_ , VAISALA, 2022 (Accessed 20.03.2023), Available: https://docs.vaisala.com/v/u/B210470EN-E/en-US.
  * [12] _Remote Road State Sensor DSC211_ , VAISALA, 2022 (Accessed 20.03.2023), Available: https://docs.vaisala.com/v/u/B210470EN-E/en-US.
  * [13] _Mobile Detector MD30_ , VAISALA, 2021 (Accessed 06.09.2023), Available: https://docs.vaisala.com/v/u/B210470EN-E/en-US.
  * [14] J. Casselgren and M. Sjödahl, “Polarization resolved classification of winter road condition in the near-infrared region,” _Applied optics_ , vol. 51, no. 15, pp. 3036–3045, 2012.
  * [15] Y. Wu, F. Liu, L. Guan, and X. Yang, “A survey of vision-based road parameter estimating methods,” in _Intelligent Computing Methodologies: 16th International Conference, ICIC 2020, Bari, Italy, October 2–5, 2020, Proceedings, Part III 16_. Springer, 2020, pp. 314–325.
  * [16] M. Nolte, N. Kister, and M. Maurer, “Assessment of deep convolutional neural networks for road surface classification,” in _2018 21st International Conference on Intelligent Transportation Systems (ITSC)_. IEEE, 2018, pp. 381–386.
  * [17] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 770–778.
  * [18] C. Szegedy, V. Vanhoucke, S. Ioffe, J. Shlens, and Z. Wojna, “Rethinking the inception architecture for computer vision,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 2818–2826.
  * [19] E. Šabanovič, V. Žuraulis, O. Prentkovskis, and V. Skrickij, “Identification of road-surface type using deep neural networks for friction coefficient estimation,” _Sensors_ , vol. 20, no. 3, p. 612, 2020.
  * [20] Z.-h. Wang, S.-b. Wang, L.-r. Yan, and Y. Yuan, “Road surface state recognition based on semantic segmentation,” _Journal of Highway and Transportation Research and Development (English Edition)_ , vol. 15, no. 2, pp. 88–94, 2021.
  * [21] K. Cordes, C. Reinders, P. Hindricks, J. Lammers, B. Rosenhahn, and H. Broszio, “Roadsaw: A large-scale dataset for camera-based road surface and wetness estimation,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2022, pp. 4440–4449.
  * [22] M. Sandler, A. Howard, M. Zhu, A. Zhmoginov, and L.-C. Chen, “Mobilenetv2: Inverted residuals and linear bottlenecks,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2018, pp. 4510–4520.
  * [23] J. Van Amersfoort, L. Smith, Y. W. Teh, and Y. Gal, “Uncertainty estimation using a single deep deterministic neural network,” in _International conference on machine learning_. PMLR, 2020, pp. 9690–9700.
  * [24] S. Roychowdhury, M. Zhao, A. Wallin, N. Ohlsson, and M. Jonasson, “Machine learning models for road surface and friction estimation using front-camera images,” in _2018 International Joint Conference on Neural Networks (IJCNN)_. IEEE, 2018, pp. 1–8.
  * [25] Z. Du, A. Skar, M. Pettinari, and X. Zhu, “Pavement friction evaluation based on vehicle dynamics and vision data using a multi-feature fusion network,” _Transportation Research Record_ , p. 03611981231165029, 2023.
  * [26] M. Abdar, F. Pourpanah, S. Hussain, D. Rezazadegan, L. Liu, M. Ghavamzadeh, P. Fieguth, X. Cao, A. Khosravi, U. R. Acharya _et al._ , “A review of uncertainty quantification in deep learning: Techniques, applications and challenges,” _Information fusion_ , vol. 76, pp. 243–297, 2021.
  * [27] N. Dewolf, B. D. Baets, and W. Waegeman, “Valid prediction intervals for regression problems,” _Artificial Intelligence Review_ , vol. 56, no. 1, pp. 577–613, 2023.
  * [28] D. Nix and A. Weigend, “Estimating the mean and variance of the target probability distribution,” in _Proceedings of 1994 IEEE International Conference on Neural Networks (ICNN’94)_ , vol. 1, 1994, pp. 55–60 vol.1.
  * [29] A. Khosravi, S. Nahavandi, D. Creighton, and A. F. Atiya, “Lower upper bound estimation method for construction of neural network-based prediction intervals,” _IEEE transactions on neural networks_ , vol. 22, no. 3, pp. 337–346, 2010.
  * [30] C. P. Robert, “Simulation of truncated normal variables,” _Statistics and computing_ , vol. 5, pp. 121–125, 1995.
  * [31] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga _et al._ , “Pytorch: An imperative style, high-performance deep learning library,” _Advances in neural information processing systems_ , vol. 32, 2019.
  * [32] V. Sovrasov, _ptflops: a flops counting tool for neural networks in pytorch framework_ , 2023 (Accessed 25.05.2023), Available: https://github.com/sovrasov/flops-counter.pytorch.
  * [33] T. Gneiting and A. E. Raftery, “Strictly proper scoring rules, prediction, and estimation,” _Journal of the American statistical Association_ , vol. 102, no. 477, pp. 359–378, 2007.


title paragraph section subsection subsubsection descriptionlabel caption
captionlabel numbers=noenddot [ type = fltalgorithm , float , ]fltalgorithm
author date

# Performance of $\mathcal{H}$-Matrix-Vector Multiplication with Floating
Point Compression

Ronald Kriemann
MPI for Mathematics i.t.S.
Leipzig Germany
<EMAIL_ADDRESS>
###### Abstract

Matrix-vector multiplication forms the basis of many iterative solution
algorithms and as such is an important algorithm also for hierarchical
matrices. However, due to its low computational intensity, its performance is
typically limited by the available memory bandwidth. By optimizing the storage
representation of the data within such matrices, this limitation can be lifted
and the performance increased. This applies not only to hierarchical matrices
but for also for other low-rank approximation schemes, e.g. block low-rank
matrices.

AMS Subject Classification: 65Y05, 65Y20, 68W10, 68W25, 68P30
Keywords: hierarchical matrices, low-rank arithmetic, data compression,
matrix-vector multiplication

##  1 Introduction

Introduced in [15] was a version of (hierarchical) low-rank arithmetic where
the matrix data, i.e., dense and low-rank blocks, was compressed using
floating point compression methods. As the standard $\mathcal{H}$-arithmetic
is typically based on dense arithmetic functions defined by the BLAS and
LAPACK function set [5], the modified $\mathcal{H}$-arithmetic in [15] was
based on the idea of decompressing all input data of arithmetic kernel
functions, executing the arithmetic kernel in standard double precision and
then compressing the output data. This way, the actual arithmetic functions
remain unchanged. Already this approach showed superior performance for the
$\mathcal{H}$-matrix-vector multiplication ($\mathcal{H}$-MVM), which is often
memory bandwidth limited and as such, any reduction of the memory size will
increase performance.

Another reason for this _semi-on-the-fly_ approach was the general compression
approach in [15], i.e., for floating point data _any_ compressor could be
used. This prevents direct arithmetic within the (unknown) compression format.

However, some of the compression schemes in [15] allow random access of
entries in the compressed storage and hence, special arithmetic functions can
be implemented. The aim of this work is to investigate the benefit of such an
approach for $\mathcal{H}$-matrix vector multiplication.

An analog strategy was used in [6] with the idea of a _memory accessor_ ,
i.e., transparent conversion between a storage and a computation format within
a sparse matrix computation. This work is therefore the application of this
concept for $\mathcal{H}$-matrix arithmetic.

A different strategy is used by mixed precision schemes ([20, 1]) where
combinations of hardware provided floating point formats are used to reduce
the memory footprint and to increase performance, partly due to faster
execution of such smaller data formats. However, these approaches are limited
in the reduction of the $\mathcal{H}$-matrix memory which is crucial on
computer systems with memory bandwidth limitations. Also, not all floating
point formats are (yet) hardware supported on all platforms, e.g., half
precision formats like BF16 or FP16.

The rest of this work is structured as follows: in Section 2 basic definitions
and algorithms for $\mathcal{H}$-matrices are introduced together with the
introduction of compression schemes for dense and low-rank data. Section 3
will discuss different strategies for $\mathcal{H}$-matrix-vector
multiplication with and without compression. Numerical experiments will be
presented in Section 4, followed by a conclusion in Section 5.

##  2 $\mathcal{H}$-Matrices

For an indexset $I$ we define the _cluster tree_ (or $\mathcal{H}$-tree) as
the hierarchical partitioning of $I$ into disjoint sub-sets of $I$:

######  Definition 2.1 (Cluster Tree)

Let $T_{I}=(V,E)$ be a tree with $V\subset\mathcal{P}(I)$. $T_{I}$ is called a
_cluster tree_ over $I$ if

  1. 1.

$I=\operatorname{root}(T_{I})$ and

  2. 2.

for all $v\in V$ with
$\operatorname{sons}(v)\neq\emptyset:v=\dot{\cup}_{v^{\prime}\in\operatorname{sons}(v)}v^{\prime}$.

A node in $T_{I}$ is also called a _cluster_ and we write $\tau\in T_{I}$ if
$\tau\in V$. The set of leaves of $T_{I}$ is denoted by $\mathcal{L}(T_{I})$.

Similar to a cluster tree we can extend the hierarchical partitioning to the
product $I\times J$ of two index sets $I,J$, while restricting the possible
set of nodes by given cluster trees $T_{I}$ and $T_{J}$ over $I$ and $J$,
respectively. Furthermore, the set of leaves will be defined by an
_admissibility condition_. In the literature, various examples of
admissibility can found, e.g. standard [11], weak [12] or off-diagonal
admissibility [9, 2].

######  Definition 2.2 (Block Tree)

Let $T_{I},T_{J}$ be two cluster trees and let $\operatorname{adm}:T_{I}\times
T_{J}\to\mathbbm{B}$. The _block tree_ $T=T_{I\times J}$ is recursively
defined starting with $\operatorname{root}(T)=(I,J)$:

$\displaystyle\operatorname{sons}(\tau,\sigma)=$
$\displaystyle\begin{cases}\emptyset,\textnormal{ if
}\operatorname{adm}(\tau,\sigma)=\textnormal{true}\vee\operatorname{sons}(\tau)=\emptyset\vee\operatorname{sons}(\sigma)=\emptyset,\\\
\left\\{(\tau^{\prime},\sigma^{\prime})\,:\,\tau^{\prime}\in\operatorname{sons}(\tau),\sigma^{\prime}\in\operatorname{sons}(\sigma)\right\\}\textnormal{
else}.\end{cases}$

A node in $T$ is also called a _block_. Again, the set of leaves of $T$ is
denoted by $\mathcal{L}(T):=\left\\{b\in
T\;:\;\operatorname{sons}(b)=\emptyset\right\\}$.

The admissibility condition is used to detect blocks in $T$ which can be
efficiently approximated by low-rank matrices with a predefined rank $k$,
i.e., blocks $b$ with $\operatorname{adm}(b)=\textnormal{true}$. The set of
all such matrices forms the set of $\mathcal{H}$-matrices:

######  Definition 2.3 ($\mathcal{H}$-Matrix)

For a block tree $T$ over cluster trees $T_{I},T_{J}$ and $k\in\mathbbm{N}$,
the set of $\mathcal{H}$-matrices $\mathcal{H}(T,k)$ is defined as

$\displaystyle\mathcal{H}(T,k):=\\{$ $\displaystyle M\in\mathbbm{R}^{I\times
J}\;:\;\forall(\tau,\sigma)\in\mathcal{L}(T):$
$\displaystyle\operatorname{rank}(M_{\tau,\sigma})\leq
k\vee\tau\in\mathcal{L}(T_{I})\vee\sigma\in\mathcal{L}(T_{J})\\}$

Here, $M_{\tau,\sigma}$ refers to the sub-block $M|_{\tau\times\sigma}$.

In practice the constant rank $k$ is typically replaced by a fixed low-rank
approximation accuracy $\varepsilon>0$ as the resulting $\mathcal{H}$-matrices
are often more memory efficient. For this we assume for an admissible block
$M_{\tau,\sigma}$:

$||M_{\tau,\sigma}-U_{\tau,\sigma}V_{\tau,\sigma}^{H}||\leq\varepsilon||M_{\tau,\sigma}||\
.$ (1)

In an analog way to $\mathcal{H}(T,k)$, the set $\mathcal{H}(T,\varepsilon)$
can be defined as the set of $\mathcal{H}$-matrices with _local_ low-rank
approximation error of $\varepsilon$. We will also use $\mathcal{H}(T)$ if
either a fixed rank or a fixed accuracy is used.

###### Remark 2.4

The set $\mathcal{H}(T)$ also includes various other formats like _block low-
rank_ (BLR) [3] or _hierarchical off-diagonal low-rank_ (HODLR) [2], as only
the clustering or the admissibility has to be chosen appropriately.

###### Remark 2.5

For $\mathcal{H}$-matrices with a full hierarchy, the set
$M_{\tau}:=\left\\{M|_{\tau,\sigma}\,:\,(\tau,\sigma)\in\mathcal{L}(T)\wedge\operatorname{rank}(M_{\tau,\sigma})\leq
k\right\\}$ of low-rank blocks for a cluster $\tau\in T_{I}$ is bounded by the
constant $c_{\textnormal{sp}}$ [10] for a particular application.

###  2.1 Compressed $\mathcal{H}$-Matrices

Floating point data in $\mathcal{H}$-matrices appears in inadmissible blocks
as dense matrices holding the coefficients and in low-rank blocks in the form
of the low-rank factors. Often these are stored in FP64 (or FP32) format.
However, due to low-rank approximation with accuracy $\varepsilon$, already an
error is introduced which is typically much larger than the unit roundoff of
FP64 (or even FP32).

In [16, 15] the FP64 storage was replaced by error adaptive floating point
compression, i.e., an optimized storage format was chose with a representation
error depending on $\varepsilon$. Different compressors are available to
implement such a direct compression of floating point data, e.g., ZFP [19] or
BLOSC [8]. Furthermore, different storage schemes based on the IEEE-754
floating point standard were examined, were the number of mantissa bits
$m_{\varepsilon}$ is chosen based on the low-rank approximation error
$\varepsilon$ as $m:=\lceil-\log_{2}\varepsilon\rceil$. Different choices for
the number of exponent bits $e$ were also examined, e.g., with 8 bits as in
the FP32 or BF16 formats (called BFL), 11 bits as in FP64 (called DFL) and an
adaptive choice based on the _dynamic range_ of the data, i.e., the base 10
logarithm of the ratio between the largest and smallest (absolute) value
(called AFLP). In all cases $m_{\varepsilon}$ was increased such that the
number of bits per value $1+e+m_{\varepsilon}$111With the additional sign bit.
is a multiple of 8 for fast byte aligned storage.

Independent on the particular choice of the compression scheme, this _direct_
compression mode is then applied to the dense data of inadmissible blocks
$M_{\tau,\sigma}\in\mathcal{L}_{\textnormal{inadm}}$ and the low-rank factors
$U_{\tau^{\prime},\sigma^{\prime}},V_{\tau^{\prime},\sigma^{\prime}}$ of
admissible blocks
$M_{\tau^{\prime},\sigma^{\prime}}\in\mathcal{L}_{\textnormal{adm}}$.

Furthermore, as described in [15], low-rank matrices permit an advanced
compression scheme with an adaptive accuracy choice for each column in the
low-rank factors. This _adaptive precision_ compression (_APLR_) is based on
the mixed precision approach described in [4]. For a block $M_{\tau,\sigma}$
we assume a rank-$k$ approximation $U\cdot V^{H}$ with
$\|M_{\tau,\sigma}-UV^{H}\|\leq\delta$. Using the singular value decomposition
we can find orthogonal matrices $W$ and $X$ and a diagonal matrix
$\Sigma=\operatorname{diag}(\sigma_{0},\ldots,\sigma_{k-1})$ with the singular
values $\sigma_{0}>\sigma_{1}>\ldots\sigma_{k-1}$ of $UV^{H}$.

If the $i$’th column $w_{i}$ of $W$ and $x_{i}$ of $X$ is stored with
precision $\delta/\sigma_{i}$, then the total approximation error is (see [15,
Section 4])

$\displaystyle\|M_{\tau,\sigma}-\widetilde{W}\Sigma\widetilde{X}^{H}\|\leq\delta+\left(2\delta
k+\delta^{2}\sum_{i=1}^{k}\frac{1}{\sigma_{i}}\right).$

With this, any direct floating point compression method $\mathcal{Z}$ can be
used to yield an improved storage method for low-rank matrices, denoted
APLR-$\mathcal{Z}$, e.g., APLR-AFLP, APLR-BFL or APLR-DFL. The main advantage
of this scheme compared to direct compression is, that in the latter case the
chosen precision is applied to the full data whereas with APLR even for a high
accuracy, a low precision may be used for some part of the data.

##  3 $\mathcal{H}$-Matrix-Vector Multiplication

We consider the update operation

$y:=\alpha Ax+y$

with an $\mathcal{H}$-matrix $A\in\mathcal{H}(T,k)$ and vectors $x$ and $y$.
The product is computed by looping over the leaf blocks of $A$ and performing
local matrix-vector multiplications, either with a dense matrix for
inadmissible blocks or in low-rank format, i.e.,
$t:=V_{\tau,\sigma}^{H}x|_{\sigma}$ followed by $y|_{\tau}:=y|_{\tau}+\alpha
U_{\tau,\sigma}t$. The full procedure is shown in Algorithm 3.

[htb]  Algorithm 1: $\mathcal{H}$-Matrix-Vector Multiplication procedure
hmvm($\alpha,A,x,y$)  for $(\tau,\sigma)\in\mathcal{L}$ do if $(\tau,\sigma$
is admissible then $y|_{\tau}:=y|_{\tau}+\alpha
U_{\tau,\sigma}V_{\tau,\sigma}^{H}x|_{\sigma}$;  else
$y|_{\tau}:=y|_{\tau}+\alpha D_{\tau,\sigma}x|_{\sigma}$;

Versions of the $\mathcal{H}$-MVM for parallel systems need to consider load
balancing due to different, not a priori known ranks in different low-rank
blocks of the $\mathcal{H}$-matrix if a fixed accuracy $\varepsilon$ is used.
Also the block structure is typically not equal throughout the matrix. This
poses a serious scalability issue for the distributed memory case (see [7, 17]
or systems with a NUMA architecture.

On shared memory systems a task-based approach can avoid these problems if the
scheduling algorithm is able to assign ready tasks to idle processors.
However, this may lead to other problems as the memory layout of the blocks
handled by a single processor may not be optimal for efficient execution. This
is of special importance because of the low computational intensity of matrix-
vector multiplication, which normally leads to a memory bandwidth limited
performance. Different optimization strategies are discussed in [13], where
especially the memory layout of the $\mathcal{H}$-matrix data is adjusted such
that memory loads are faster.

Another issue with shared memory programming is handling potential collisions
when writing to the same memory positions, e.g., with matrix blocks
$A_{\tau,\sigma}$ and $A_{\tau,\sigma^{\prime}}$ handled by different
processors writing simultaneously to $y|_{\tau}$. Solutions to this problem
involve _atomic_ updates [14] or reduction of thread local results [13]. A
reduction approach of local results is also the default choice for the
distributed memory case [7, 17].

An alternative approach is a collision free design in which the memory blocks
are scheduled to the processors in a way to prevent simultaneous writing to
the same memory positions. Such a method is used in the following.

Let
$\mathcal{A}_{\tau}:=\left\\{A_{\tau,\sigma}\;:\;(\tau,\sigma)\in\mathcal{L}(T)\right\\}$
be the set of all matrix blocks in $A$ with identical row cluster $\tau$ and
let $\mathcal{A}:=\left\\{\mathcal{A}_{\tau}\;:\>\tau\in T_{I}\right\\}$ be
the set of all such block lists. Since $\mathcal{A}$ is defined based on
$T_{I}$, it can be considered to be structurally identical to the cluster
tree. Due to its definition, the number of matrix blocks in any
$\mathcal{A}_{\tau}$ is bounded by $c_{\textnormal{sp}}$ (see 2.5) and
therefore independent on the dimension of the matrix..

Now let $\tau_{0},\ldots,\tau_{\ell}$ be clusters of $T_{I}$ with identical
level, i.e., $0pt{\tau_{i}}=0pt{\tau_{j}},0\leq i,j\leq\ell$. Then, for any
$0\leq i,j\leq\ell$ the matrix-vector products in the corresponding sets
$\mathcal{A}_{\tau_{i}}$ and $\mathcal{A}_{\tau_{j}}$ can be computed in
parallel since $\tau_{i}\cap\tau_{j}=\emptyset$.

For any $\tau,\sigma\in T_{I}$ with $0pt{\tau}\neq\sigma$ the sets
$\mathcal{A}_{\tau}$ and $\mathcal{A}_{cls}$ can only be executed in parallel
if $\tau\cap\sigma=\emptyset$. However, due to the definition of $T_{I}$ if
$\tau\cap\sigma\neq\emptyset$ then either $\tau\subseteq\sigma$ or
$\sigma\subseteq\tau$ holds. Therefore, if $T_{I}$ is traversed from root to
bottom with execution of matrix blocks in a given $\mathcal{A}_{\tau}$ before
proceeding to the sons in $\mathcal{S}(\tau)$, any race condition when
accessing $y$ is prevented. This procedure is implemented in Algorithm 3.

[htb]  Algorithm 2: Parallel $\mathcal{H}$-Matrix-Vector Multiplication
procedure phmvm($\alpha,\tau,\mathcal{A},x,y$)  for all
$A_{\tau,\sigma}\in\mathcal{A}_{\tau}$ do if $(\tau,\sigma$ is admissible then
$y|_{\tau}:=y|_{\tau}+\alpha U_{\tau,\sigma}V_{\tau,\sigma}^{H}x|_{\sigma}$;
else $y|_{\tau}:=y|_{\tau}+\alpha D_{\tau,\sigma}x|_{\sigma}$;  parallel for
$(\tau^{\prime}\in\mathcal{S}(\tau)$ do
phmvm($\alpha,\tau^{\prime},\mathcal{A},x,y$)

Figure 1: Roofline plot for $\mathcal{H}$-MVM.

When using this algorithm with the architecture used in Section 4, almost
optimal performance is achieved as shown in the (empirical) roofline plot in
Figure 1. Please note, that the plot shows values obtained for different
problem sizes and therefore data sizes, demonstrating its performance
consistency.

In principle, the computation of all products for matrix blocks in
$\mathcal{A}_{\tau}$ in Algorithm 3 can be further parallelized using a
reduction scheme. One could also combine memory layout optimizations from [13]
easily with this approach. However, as already the performance limit is
reached by the above procedure, only minor improvements by such modifications
are possible.

In Algorithm 3 no parallelism of the per matrix block products was considered,
which is a major drawback of this procedure if a block structure like HODLR
[2] is used. There, the most time consuming computations are performed on the
upper levels of the block cluster tree, where only a few processors may be
used in parallel. However, as the numerical results in Section 4 demonstrate,
typical $\mathcal{H}$-matrix block structures do not show these problems, at
least not for the number of processors cores considered in this work.

###  3.1 Compressed $\mathcal{H}$-Matrix-Vector Multiplication

The main interest of this work is in the performance of $\mathcal{H}$-MVM when
using floating point compression. This was already the topic of [15]. There,
for a dense or low-rank block, the compressed data was first fully converted
into the computation format and only then the local matrix-vector
multiplication was performed in double precision. This way, the standard
arithmetic kernels, typically optimized by hardware vendors, could be reused.

In [6] the concept of a _memory-accessor_ is described which implements on-
the-fly conversion between the storage format and the computation format
during the arithmetic. Since $\mathcal{H}$-MVM only uses decompression this
approach is easier to apply in this case compared to the full
$\mathcal{H}$-arithmetic. Also, since $\mathcal{H}$-MVM is often memory
bandwidth limited it may be more forgiving for a (potentially) less heavily
optimized implementation.

In any case, such an approach requires fast access to the compressed values,
which holds for the above described compression schemes AFLP, BFL and DFL but
much less so for ZFP or BLOSC. Furthermore, especially ZFP showed a low
performance compared to other formats (see [15]). However, in principle by
tightly coupling the compression scheme with the matrix-vector multiplication,
any compression format could be used. As such, the restriction to AFLP, BFL
and DFL in this work should be considered a proof-of-concept.

Since Algorithm 3 only uses dense matrix-vector multiplication, one only needs
to focus on this function. The actual implementation of the compressed version
is straight-forward, without particular optimizations, aside from standard
code reorganization due to the used column-major storage scheme and shown in
Algorithm 3.1 for the application of a non-transposed $n\times m$ matrix $D$.
Only the access to the coefficients in $D$ is replaced by the corresponding
coefficient decompression.

[htb]  Algorithm 3: Matrix-vector multiplication with compressed dense
$n\times m$ matrix procedure zmvm(in: $D,x$, inout: $y$)  for $0\leq j<m$ do
for $0\leq i<n$ do $y_{i}:=y_{i}+\textnormal{{{decompress}($D_{ij}$)}}x_{j}$;

For the mixed precision approach in [4] in the predecessor in [20] (using FP64
and FP32) which are based on hardware supported floating point formats, one
has the advantage of performing the computations for the corresponding
floating point hardware natively, i.e., without converting to FP64, thereby
potentially increasing the performance. As this had no negative side effects
on the error for the model problems used in Section 4 this was also used in
the following.

##  4 Numerical Experiments

Figure 2: Compression rates compared to uncompressed $\mathcal{H}$-matrices
for fixed accuracy of $\varepsilon=10^{-6}$ (left) and a fixed problem size of
$n=2.097.152$ (right).

The model problem is based on a boundary element discretization for the
Laplace single layer potential (Laplace SLP) while the domain is defined by
the unit sphere:

$\int_{\Omega}\frac{1}{\|x-y\|}u(x)dy=f(x),\quad x\in\Omega$ (2)

with $\Omega=\left\\{x\in\mathbbm{R}^{3}:\|x\|_{2}=1\right\\}$. Piecewise
constant ansatz functions are used for the discretization. Furthermore,
standard admissibility

$\min\left\\{\operatorname{diam}(t),\operatorname{diam}(s)\right\\}\leq\eta\operatorname{dist}(\tau,\sigma)$

is applied for setting up the block tree.

All experiments are performed on an AMD Epyc 9554 CPU with 64 cores in total
and 12 32GB DDR5-4800 memory DIMMs. For parallelization Intel TBB v2021.11 was
used while Intel oneMKL v2024.0 provided the BLAS and LAPACK functions for the
uncompressed case. Please note, that the sequential version was chosen as all
parallelization is performed within the $\mathcal{H}$-arithmetic itself.
Furthermore, the AVX512 code path in MKL was activated. All code was compiled
using GCC v12.3.

The algorithms described in this work are implemented in the open source
software HLR222http://libhlr.org, program: _mpmvm_. For the numerical
experiments version _9cdb804_ was used.

For runtime results the median of ten runs is presented, while only little
variations between each run was observed.

Aside from the compressors AFLP, BFL, DFL the mixed precision formats _MP-3_
using FP64, FP32 and BF16333BF16 was preferred over FP16 due to much faster
conversion from FP64/FP32. (as in [4]) and _MP-2_ with FP64 and FP32 (as in
[20]) are used.

First, the compression ratio for the different schemes is shown in Figure 2
for a comparison on the memory savings. As can be seen the compression ration
is slightly improving with a growing problem size if a fixed accuracy is
chosen. This corresponds to the increasing portion of low-rank memory compared
to memory associated to inadmissible blocks in $\mathcal{H}$-matrices with
larger problems and is more pronounced for the mixed precision formats as
there inadmissible blocks are not compressed. Due to a larger number of
exponent bits in BFL and DFL, these formats have a higher memory consumption
compared to AFLP.

For a fixed problem size and varying accuracy, APLR compression yields a high
compression for low and high accuracies. In contrast to this MP-3 and MP-2
show a much worse compression if low accuracy is used, again due to the fact
that it is applied to low-rank blocks only. However, for higher accuracies
especially MP-3 is able to come closer to the other formats.

Figure 3: Roofline plot for $\mathcal{H}$-MVM without (red) and with
compression using AFLP (blue) for a fixed accuracy of $\varepsilon=10^{-6}$.

Figure 4: Speedup of compressed $\mathcal{H}$-MVM vs. uncompressed
$\mathcal{H}$-MVM (base line) for fixed accuracy of $\varepsilon=10^{-6}$
(left) and a fixed problem size of $n=2.097.152$ (right).

Figure 3 shows the the performance of uncompressed and compressed
$\mathcal{H}$-MVM in an empirical roofline plot with a maximal memory
bandwidth of 390 GB/s and maximal floating point performance of 3.8 TFlop/s,
measured by the Likwid tool [18]. The arithmetic intensity of
$\mathcal{H}$-MVM is well within the bandwidth limited regime. While the
uncompressed multiplication is close to the limit, the compressed
$\mathcal{H}$-MVM is slightly less optimal. However, the latter does not make
use of the optimized matrix-vector product from the Intel MKL library and also
has the additional overhead of the coefficient decompression.

Nevertheless, the performance is significantly increased by lowering the
memory bandwidth requirements. The runtime speedup compared to the
uncompressed multiplication is shown in Figure 4 for all compression schemes.
There a better compression ratio directly translates into a better performance
with AFLP yielding best results.

To demonstrate, that this performance improvement is not restricted to
$\mathcal{H}$-matrices, the same problem was computed using the Block Low-Rank
approach from [3]. Here, the weak admissibility [12] was used. The results for
BLR matrix-vector multiplication (BLR-MVM) are shown in Figure 5. The maximal
problem sizes are smaller due to the higher memory demands of the BLR format.

Figure 5: Speedup of compressed BLR-MVM vs. uncompressed BLR-MVM (base line)
for fixed accuracy of $\varepsilon=10^{-6}$ (left) and a fixed problem size of
$n=1.048.576$ (right).

While the general picture is similar to the $\mathcal{H}$-matrix case, the
runtime improvements compared to the uncompressed case are even slightly
bigger. Only the mixed precision formats do show a worse behavior for small
problem sizes.

##  5 Conclusion

Matrix-vector multiplication for $\mathcal{H}$-matrices can benefit
significantly from an optimized memory representation of the dense and low-
rank data on platforms with slow memory access as is the case for many CPU
based computations. With APLR an efficient compressor is available which
permits such optimizations if combined with floating point compression schemes
with fast access to individual values. The such presented algorithm not only
performs well for $\mathcal{H}$-matrices but may by used also for other forms
of low-rank storage like BLR. It will be interesting to see if these concepts
can be applied to more complex forms of $\mathcal{H}$-arithmetic, e.g. matrix-
multiplication or LU factorization.

## References

  * [1] S. Abdulah et al.  “Accelerating Geostatistical Modeling and Prediction With Mixed-Precision Computations: A High-Productivity Approach With PaRSEC”  In _IEEE Transactions on Parallel and Distributed Systems_ 33.4, 2022, pp. 964–976  DOI: 10.1109/TPDS.2021.3084071
  * [2] S. Ambikasaran and E. Darve  “An $\mathcal{O}(N\log N)$ Fast Direct Solver for Partial Hierarchically Semi-Separable Matrices”  In _Journal of Scientific Computing_ 57.3, 2013, pp. 477–501
  * [3] P. Amestoy et al.  “Improving Multifrontal Methods by Means of Block Low-Rank Representations”  In _SIAM Journal on Scientific Computing_ 37.3, 2015, pp. A1451–A1474  DOI: 10.1137/120903476
  * [4] P. Amestoy et al.  “Mixed precision low-rank approximations and their application to block low-rank LU factorization”  In _IMA Journal of Numerical Analysis_ , 2022  DOI: 10.1093/imanum/drac037
  * [5] E. Anderson et al.  “LAPACK Users’ Guide”  Philadelphia, PA: Society for IndustrialApplied Mathematics, 1999
  * [6] Hartwig Anzt, Thomas Grützmacher and Enrique S. Quintana-Ortí  “Toward a modular precision ecosystem for high-performance computing”  In _International Journal of High Performance Computing Applications_ 33.6, 2019, pp. 1069–1078  DOI: 10.1177/1094342019846547
  * [7] M. Bebendorf and R. Kriemann  “Fast parallel solution of boundary integral equations and related problems”  In _Computing and Visualization in Science_ 8.3, 2005, pp. 121–134  DOI: 10.1007/s00791-005-0001-x
  * [8] Blosc Development Team  “A fast, compressed and persistent data store library”, https://blosc.org, 2009-2023
  * [9] S. Chandrasekaran et al.  “Some Fast Algorithms for Sequentially Semiseparable Representations”  In _SIAM Journal on Matrix Analysis and Applications_ 27, 2005, pp. 341–364
  * [10] L. Grasedyck and W. Hackbusch  “Construction and arithmetics of ${\mathcal{H}}$-matrices”  In _Computing_ 70, 2003, pp. 295–334
  * [11] W. Hackbusch and B.N. Khoromskij  “A sparse $\mathcal{H}$-matrix arithmetic. Part II. Application to multi-dimensional problems”  In _Computing_ 64.1, 2000, pp. 21–47
  * [12] W. Hackbusch, B.N. Khoromskij and R. Kriemann  “Hierarchical Matrices Based on a Weak Admissibility Criterion”  In _Computing_ 73.3, 2004, pp. 207–243  DOI: 10.1007/s00607-004-0080-4
  * [13] T. Hoshino, A. Ida and T. Hanawa  “Optimizations of H-matrix-vector Multiplication for Modern Multi-core Processors”  In _2022 IEEE International Conference on Cluster Computing (CLUSTER)_ , 2022, pp. 462–472  DOI: 10.1109/CLUSTER51413.2022.00056
  * [14] A. Ida, T. Iwashita, T. Mifune and Y. Takahashi  “Parallel Hierarchical Matrices with Adaptive Cross Approximation on Symmetric Multiprocessing Clusters”  In _Journal of Information Processing_ 22.4, 2014, pp. 642–650  DOI: 10.2197/ipsjjip.22.642
  * [15] R. Kriemann  “Hierarchical Lowrank Arithmetic with Binary Compression”, 2023  arXiv:2308.10960 [cs.MS]
  * [16] R. Kriemann et al.  “High-Performance Spatial Data Compression for Scientific Applications”  In _Euro-Par 2022: Parallel Processing_ Cham: Springer International Publishing, 2022, pp. 403–418
  * [17] Y. Li, J. Poulson and L. Ying  “Distributed-Memory $\mathcal{H}$-Matrix Algebra I: Data Distribution and Matrix-Vector Multiplication”  In _CSIAM Transactions on Applied Mathematics_ 2.3, 2021, pp. 431–459  DOI: 10.4208/csiam-am.2020-0206
  * [18] “LIKWID”, 2024  URL: https://github.com/RRZE-HPC/likwid/
  * [19] P. Lindstrom  “Fixed-Rate Compressed Floating-Point Arrays”  In _IEEE Transactions on Visualization and Computer Graphics_ 20.12, 2014, pp. 2674–2683
  * [20] R. Ooi et al.  “Effect of Mixed Precision Computing on H-Matrix Vector Multiplication in BEM Analysis”  In _Proceedings of the International Conference on High Performance Computing in Asia-Pacific Region_ , HPCAsia ’20  , Fukuoka, Japan,: Association for Computing Machinery, 2020, pp. 92–101  DOI: 10.1145/3368474.3368479