# Joint Multimodal Transformer for Emotion Recognition in the Wild

Paul Waligora1, Haseeb Aslam1, Osama Zeeshan1, Soufiane Belharbi1, Alessandro
Lameiras Koerich1,
Marco Pedersoli1, Simon Bacon2, and Eric Granger1
1 LIVIA, Dept. of Systems Engineering, ETS Montreal, Canada
2 Dept. of Health, Kinesiology & Applied Physiology, Concordia University,
Montreal, Canada
<EMAIL_ADDRESS>
###### Abstract

Systems for multimodal emotion recognition (MMER) can typically outperform
unimodal systems by leveraging the inter- and intra-modal relationships
between, e.g., visual, textual, physiological, and auditory modalities. In
this paper, an MMER method is proposed that relies on a joint multimodal
transformer for fusion with key-based cross-attention. This framework aims to
exploit the diverse and complementary nature of different modalities to
improve predictive accuracy. Separate backbones capture intra-modal
spatiotemporal dependencies within each modality over video sequences.
Subsequently, a joint multimodal transformer fusion architecture integrates
the individual modality embeddings, allowing the model to capture inter-modal
and intra-modal relationships effectively. Extensive experiments111Code:
https://github.com/PoloWlg/Joint-Multimodal-Transformer-6th-ABAW. on two
challenging expression recognition tasks: (1) dimensional emotion recognition
on the Affwild2 dataset (with face and voice), and (2) pain estimation on the
Biovid dataset (with face and biosensors), indicate that the proposed method
can work effectively with different modalities. Empirical results show that
MMER systems with our proposed fusion method allow us to outperform relevant
baseline and state-of-the-art methods.

Keywords: Affective Computing, Multi-Modal Fusion, Audio-Visual Fusion,
Transformers, Cross-Attention, Valence-Arousal Estimation.

Figure 1: (Top) Illustration of the vanilla multimodal transformer fusion
architecture. (Bottom) Multimodal transformer fusion with joint representation
(in red) for two input sources, A and B.

##  1 Introduction

Human-computer interaction is applied in a wide range of real-world scenarios,
e.g., health care, IoT, and autonomous driving. Researchers have classified
human emotions in different ways, most notably according to discrete
categories, ordinal intensity levels, and the valence/arousal circumplex [1].
With the recent advancements in deep learning and sensor technologies,
research in affective computing has evolved from lab-controlled to real-world
(in the wild) scenarios. In the latter, human emotions are usually expressed
over a broader spectrum beyond the six basic categorical expressions - anger,
disgust, fear, happy, sad, and surprise [9]. Therefore, there is much interest
in analyzing and modeling complex and subtle expressions of emotions in real-
world scenarios. For instance, the spectrum of emotions can be formulated as
dimensional emotional recognition (ER), where complex human emotions are
represented along the axes of arousal (intensity) and valence (positiveness).

Multimodal fusion has been widely explored for the problem of video-based ER
in the literature [50, 28]. For instance, audio and visual modalities may
provide complementary and redundant information over a video sequence. These
relationships must be captured to model the intricacies of human emotions
effectively. Furthermore, effectively capturing both the intra-modal temporal
dependencies within the audio and visual modalities and the inter-modal
association across the audio and visual modalities is crucial in developing an
effective AER system[47].

Several methods have been proposed for video-based ER and recurrent networks
have been employed to capture the intra-modal temporal dependencies from video
sequences [50, 28, 5]. Recently, attention-based methods have been introduced
to extract features that are the most relevant to downstream tasks. Cross-
attention-based methods have also been [48] employed to capture the inter-
modal association between the audio, visual, and other modalities. Lu et al.
[35] proposed ViLBERT, which was the seminal work in multimodal co-attention.
Since then, many transformer-based cross-attention methods have been proposed
[43, 57]. These methods, however, cannot effectively capture the intra-modal
temporal dynamics. Further, they specialize in capturing the complementary
information among the modalities but do not have the mechanism to capture the
redundant information explicitly.

The proposed method introduces the third branch with the joint representation
of the multiple modalities, as shown in Figure 1. By incorporating a joint
representation branch, the model can access additional contextual information
that might not be fully captured by cross-attention alone. This can help
improve the model’s understanding of complex relationships between the input
sequences. Further, the proposed method becomes more robust to noise or
irrelevant information present in individual sequences. This helps in
mitigating the sensitivity of cross-attention to noisy inputs and improves the
overall performance of the system.

Our main contributions are summarized as follows.
(1) This paper proposes a joint multimodal transformer fusion framework that
leverages joint modality representations. It captures inter- and intra-modal
information in videos using key-based cross-attention. (2) A joint multimodal
transformer multimodal fusion that exploits the redundant and complementary
associations among modalities. (3) An extensive set of experiments on the
challenging emotion recognition (pain estimation on BioVid and dimension
valence-arousal assessment on Affwild2) datasets indicate that our proposed
multimodal fusion architecture can outperform relevant baseline and state-of-
the-art methods.

##  2 Related Work in Emotion Recognition

###  2.1 Multimodal Methods

MMER refers to the integration of multiple sources of information (modalities)
to improve the accuracy and robustness of automated emotion recognition
systems, at the expense of complexity. These modalities typically include
visual, audio, textual, and physiological. The seminal work in multimodal deep
learning was proposed by Ngiam et al. [39], where the features from the audio
and visual modalities were extracted separately, and then autoencoders and
Restricted Boltzmann Machines were used to feature fusion. Tzirakis et al.
[55] proposed one of the early approaches for A-V fusion for dimensional
emotion recognition, in which the visual features were extracted using a
ResNet50 and the audio features were obtained using a 1D convolutional neural
network (CNN). The modality-specific features were concatenated and fed to a
recurrent net for simultaneous temporal modeling and modality fusion. An
empirical study was presented by Juan et al. [42], where the authors studied
the impact of fine-tuning multiple layers in a pretrained CNN for the visual
modality. A two-stream autoencoder with a long short-term memory (LSTM)
network was proposed by Nguyen et al. [40] to jointly learn and compact
representative features from the visual and audio modalities. A knowledge
distillation-based approach was investigated by Schonevald et al. [50] for
visual modality. For the audio modality, spectrograms were obtained and fed to
a CNN model, and the two modalities were fused using a recurrent net. A novel
self-distillation scheme was put forward by Deng et al. [6] to overcome the
problem of noisy labels in a multitasking setting. A two-stream aural visual
(TSAV) network was proposed by Kuhnke et al. [28], in which the audio features
were extracted using a ResNet18, and the visual features were extracted using
a 3D-CNN. The obtained embeddings were fed to a specially designed TSAV
network for information fusion.

Pain classification is one of the primary problems in affective computing.
Researchers have proposed many multimodal datasets for the pain estimation
task. The facial activity descriptors method for pain estimation was
introduced by Werner et al. [59]. Dragomir et al. [7] propose a subject-
independent method from facial images with a residual learning technique. A
Sparse LSTM-based method was proposed by Zhi et al. [65] to solve the problem
of vanishing gradients in temporal learning. Morabit et al. [38] proposed a
data-efficient image transformer. To process multiscale electrodermal activity
signals, a SE-Net-based network was proposed by Lu et al. [36]. Multimodal
solutions to fuse the physiological and visual modalities were proposed by
Werner et al. [59], Kachele et al. [13], and Zhu et al. [65]. Physiological
signals are more discriminative for pain classification than the visual
modality.

###  2.2 Attention-Based and Transformer Methods

Since its inception, attention models have shown extraordinary performance in
many applications. These models have been extensively investigated for
capturing the inter and intra-modal associations between the audio and visual
modalities for tasks like action localization [32], A-V event localization
[8], and multimodal emotion recognition [43]. An attention-based fusion
mechanism was proposed by Zhang et al. [64], 3D-CNNs and 2D-CNNs were used to
extract multi-features in the visual modality, and for the audio modality, a
2D-CNN was used to learn representation from spectrograms. Specialized scoring
functions were used to re-weight the audio and visual features. Recently,
cross-modal attention has shown promising results because of its ability to
model inter-modal relationships. Srinivas et al. [43] explored a transformer
network with encoder layers, where cross-modal attention is used to fuse audio
and visual features for continuous arousal/valence prediction in the wild.
Tzirakis et al. [54] explored the idea of cross-attention in conjunction with
self-attention. The authors proposed a transformer-based fusion architecture.
Although the methods mentioned above have used cross-modal attention with
transformers, they do not have any explicit mechanism to capture semantic
relevance between the A-V features, particularly the intra-modal correlations.
Zhang et al. [63] proposed a method for A-V fusion using leader-follower
attentive fusion for continuous arousal/valence prediction. Attention weights
are combined with the encoded visual and audio features. Cross attention
presented in Praveen et al. [48] has shown a substantial increase in
performance by using cross-correlation across the individual features. In
contrast, our proposed method uses key-based cross-attention in multimodal
transformers and explores the idea of feeding the joint A-V feature vector. By
feeding the joint A-V feature representation, the proposed method effectively
captures the inter- and intra-modal relationships simultaneously by
interacting across itself and the other modalities.

Huang et al. [11] investigated the idea of multi-head attention in
transformer-based fusion architecture, which was further combined with LSTM to
capture the high-level representations. Tran et al. [53] proposed a cross-
modal transformer architecture that consisted of a multimodal cross-modal
attention block, where the Queries were generated from one modality and the
key-values were generated from the other modality. Le et al. [31] put forward
an end-to-end transformer-based fusion mechanism for multilabel multimodal
emotion classification; the model consisted of three parts: i) three backbone
networks for visual, audio, and textual feature extractor, ii) a transformer
network for information fusion, and iii) classification network. Zhou et al.
[66] proposed a transformer-based fusion scheme along with the temporal
convolutional network (TCN); the audio and visual features were extracted
using pretrained backbones followed by a TCN, the output of TCN was
concatenated and fed to a transformer encoder block. A multilayer perceptron
(MLP) was then used for the final prediction.

All the aforementioned transformer-based fusion architectures primarily focus
on intermodality correlation. In contrast, the proposed method, in addition to
modeling the intermodality relationships to capture the complementarity
between modalities, we explicitly feed the joint (combined) features to the
multimodal transformer to introduce redundancy. By incorporating this third
joint representation branch, the proposed model has access to enhanced
contextual information that might not be fully captured by cross-attention
alone. Doing this improves the model’s understanding of complex relationships
between the input sequences. Further, the proposed method becomes more robust
to noise or irrelevant information present in individual sequences. Having
this third joint representation allows the model to dynamically focus on this
newly introduced information in sequences where both modalities are noisy
simultaneously. This helps in mitigating the sensitivity of cross-attention to
noisy inputs and improves the overall performance of the system.

Figure 2: An overview of the proposed joint multimodal transformer model for
A-V fusion. The audio and visual modalities are cross-attended using
transformer blocks. The JMT block also takes in the joint representation
(shown with red arrows). The output of the cross-attended features is
concatenated, and an FC layer is used for valence/arousal prediction.

##  3 Proposed Approach

The proposed method is a hierarchical fusion mechanism, where the intra-
modality features are combined using transformer-based self-attention, and
cross-modality features are fused using transformer-based cross-attention.
Further, we feed a third joint representation to the joint transformer module
(JTM) to enhance robustness. The $K$ (key matrix), $V$ (value matrix), and $Q$
(query matrix) vectors are shared among the six transformer blocks. In the
end, the output of these six blocks is again fed to a transformer self-
attention block to dynamically weigh the most relevant representations. The
final prediction is made using fully connected (FC) layers.

###  3.1 Modality Specific Feature Extraction

In the first step, modality-specific features are extracted using backbones.
The proposed method allows to combine multiple backbones for each modality to
improve system robustness. The extracted feature vectors from each backbone
are fed to a transformer self-attention block. The combined feature vector is
the representation of the particular modality. For example, to capture
information about a person’s emotional state, we can use an R(2+1)D CNN
pretrained on the Kinetics-400 dataset [29] to extract visual features. For
the audio modality, we could extract features using a ResNet18 [10] CNN with a
GRU [4].

###  3.2 Multi-transformer Attention Fusion

We define $\bm{F_{A}}$ as the deep features extracted from backbone A, and
$\bm{F_{B}}$ as the deep features extracted from the backbone B in Figure 2.
Given $\bm{F_{B}}$ and $\bm{F_{A}}$, the joint feature representation is
obtained by concatenating $\bm{F_{B}}$ and $\bm{F_{A}}$ feature vectors:

${\bm{J}}=[{\bm{F}}_{\mathbf{B}};{\bm{F}}_{\mathbf{A}}]\;,$ (1)

where ${[\cdot;\cdot]}$ denotes a concatenation operation. The concatenated
${\bm{J}}$ features are then fed to an FC layer for dimensionality reduction
of the joint feature representation to yield $\bm{J_{AB}}$. We now have three
key sources of information: ${\bm{F}}_{\mathbf{B}}$, ${\bm{F}}_{\mathbf{A}}$
and $\bm{J_{AB}}$, which have the same dimensionality.

Each representation is then fed to a specific encoder. Our model is composed
of three different encoders, one for each type of feature
${\bm{F}}_{\mathbf{B}}$, ${\bm{F}}_{\mathbf{A}}$, and $\bm{J_{AB}}$. Each
encoder consists of a multi-head-self-attention (Eq. 2) followed by a fully
connected feed-forward network. Residual connection and layer normalization
are performed around both of these layers. The key is used to associate a
sequence to a key value, the value matrix holds information that is ultimately
used to compute the output of the attention mechanism, and the query matrix
represents a set of vectors used to query the key-value pairs. The ${K}$,
${V}$, ${Q}$ matrix are calculated this way: ${K={X}{W_{K}}}$,
${Q={X}{W_{Q}}}$, ${V={X}{W_{V}}}$. ${X}$ corresponds to one of the sources of
information and ${W_{K}}$, ${W_{Q}}$, ${W_{V}}$ are the weights of the key,
query, value matrices respectively. The output values of the self-attention
layers are given by:

${Attention(Q,K,V)=softmax({K}{Q^{T}}/{d_{k}})V}$ (2)

With self-attention layers, each encoder focuses independently on important
cues related to its respective source of information. Afterward, each encoder
embedding is combined by utilizing six cross-modal attention layers where the
query matrix Q is shared with the key K and value V matrix of the other source
of information. Sharing this matrix between each source of information helps
the model add redundancy and complement the visual and audio modalities, thus
improving its performance. At the output of each of the six cross-attention
modules, the feature vector of dimension 512 is output. These six feature
vectors are then stacked to form a sequence, which is then fed to a
transformer self-attention block. This block dynamically selects and weighs
these feature vectors. The final attended features are fed to an FC layer for
final prediction.

The model aims to maximize the Concordance Correlation Coefficient (CCC) [30],
which is common in dimensional emotion recognition. To achieve this, we
minimize the following loss:

$\\\ {L}_{c}=1-\rho_{c}=\
1-\frac{2\rho_{xy}^{2}}{\rho_{x}^{2}+\rho_{y}^{2}+(\mu_{x}-\mu_{y})^{2}}\ $
(3)

where ${\rho_{xy}^{2}}$ is the covariance between the predictions and the
ground truth, ${\rho_{x}^{2}}$ and ${\rho_{y}^{2}}$ the variance of the
prediction and the ground truth, ${\mu_{x}}$ and ${\mu_{y}}$ the mean of
prediction and the ground truth.

Table 1: Comparison with the state-of-the-art pain estimation methods on the Biovid Heat Pain Database. The highest score is indicated in bold.  Method | Modality | CV Scheme | Accuracy (%)
---|---|---|---
Werner et al. [60] ICPR 2014 | EDA, ECG, EMG + Video | 5-FOLD | 80.6
Werner et al. [58] IEEE TAC 2016 | Video | LOSO | 72.4
Kachele et al. [12] IEEE IJSTSP 2016 | EDA, ECG, EMG | LOSO | 82.73
Lopez et al. [33] ACIIW 2017 | EDA, ECG | 10-FOLD | 82.75
Lopez et al. [34] EMBC 2018 | EDA | LOSO | 74.21
Thiam et al. [51] Sensors 2019 | EDA | LOSO | 84.57
Wang et al. [56] EMBC 2020 | EDA, ECG, EMG | LOSO | 83.3
Pouromran et al. [45] PLoS ONE 2021 | EDA | LOSO | 83.3
Thiam et al. [52] Frontiers 2021 | EDA, ECG, EMG | LOSO | 84.25
Phan et al. [44] IEEE Access 2023 | EDA, ECG, EMG | LOSO | 84.8
1D CNN | EDA | 5-FOLD | 77.2
R3D model | Video | 5-FOLD | 72.9
Feature concatenation | EDA + Video | 5-FOLD | 83.5
Vanilla transformer fusion | EDA + Video | 5-FOLD | 87.8
Joint Multimodal Transformer (ours) | EDA + Video | 5-FOLD | 89.1

##  4 Experimental Methodology

###  4.1 Datasets

#### Affwild2

was put forward by [23]. It is one of the largest and most comprehensive
datasets for affective computing focused in-the-wild scenarios. The dataset is
also a part of the Affective Behavior Analysis in the Wild (ABAW) challenge
[2, 22, 21, 17, 16, 3, 15, 27, 26, 20, 18, 25, 24, 19, 61]. The dataset
comprises $564$ videos that were collected from YouTube. These videos are
labeled for three main affective computing tasks: i) categorical expression
recognition, ii) continuous valence arousal prediction, and iii) action unit
detection. The dataset comes with a train, validation, and test split with
341, 71, and 152 videos, respectively. The continuous valence/arousal
annotation set is used to validate the proposed method.

#### Biovid Heat Pain Database:

The BioVid Heat Pain Database [60] comprises 87 subjects where heat pain was
induced experimentally at the right arm with four different intensities. The
data comprises video and depth map video from a Kinect camera, galvanic skin
response (EDA), electromyograph (EMG) on trapezius muscle, and
electrocardiogram (ECG). The Biovid dataset has various partitions from Part A
through E. These partitions differ in the modalities, annotations, and tasks.
We use Part A of the dataset and utilize the video and raw EDA data to
validate our proposed method.

###  4.2 Implementation Details

#### Affwild2:

In the visual modality, cropped and aligned facial images provided with the
dataset are used [27]. Black frames (zero pixels) replace missing frames in
the visual modality. These images are then fed to a 3D network that takes the
input size of 224$\times$224\. A clip length of 8 is used, which makes up a
sub-sequence of 64 frames. Each sub-sequence contains eight clips. A dropout
with a value of 0.8 was used in the linear layers for network regularization.
The network was optimized using the stochastic gradient descent (SGD)
iterative method with an initial learning rate (LR) of $10^{-3}$. The batch
size used for the visual modality was 8. For further generalization, random
cropping and random horizontal flips were added as data augmentation. The
maximum number of epochs is 50, with early stopping for model selection. In
the audio modality, the audio from each video is separated and resampled for
44 kHz. Following the segmentation in the visual modality, small vocal
segments are extracted corresponding to the length of the clip. A ResNet18 is
used to extract the features from the spectrograms. A Discrete Fourier
transform with a length of 1024, a hop length of 10 msec, and a window length
of 20 msec is used to obtain spectrograms of each vocal segment. The
spectrograms have a resolution of 64$\times$107 pixels, corresponding to a
single clip in the visual modality. Other preprocessing of spectrograms
include mean and variance normalization, as well as conversion to log-power
spectrum. The first convolutional layer in the pretrained ResNet model is
adapted to take in the single-channel spectrograms. The learning rate of
$1\times 10^{-2}$ is used and Adam optimizer is used to optimize the network.
64 batch size is set for the audio modality.

The audio and visual backbones are frozen to train the A-V fusion network and
only train the whole transformer-based fusion model. Each of these backbones
outputs deep feature vectors of dimension 512. These features are concatenated
(to obtain the joint feature representation, a feature vector of dimension
1024) and then fed to an FC layer to reduce dimensionality to 512. Each of
these features is fed to the fusion model, as we can see in Figure 2. At the
output of each of the six cross-attention modules, the feature vector of
dimension 512 is outputted. These six feature vectors are then stacked to form
a sequence, which is then fed to a transformer self-attention block. This
block dynamically selects and weighs these feature vectors. The final attended
features are fed to a FC layer for final prediction. We perform a grid search
to find the optimal learning rate and batch size for our fusion network. Thus,
we use a list of learning rates: [$8\times 10^{-4}$, $6\times 10^{-4}$,
$3\times 10^{-4}$] and an SGD optimizer to optimize the fusion network. We use
a batch size of 32, and the maximum number of epochs is 5, with early stopping
for fusion model selection. We select the best model among the learning rates
listed before.

#### Biovid:

In the visual modality, the faces are cropped and aligned using an MTCNN. We
apply the frame retention strategy for missing frames where the face is not
visible, and the MTCNN cannot capture any frame. Further, to ensure noise-free
input to the visual model, we clip the first 2 seconds of the video and the
last 0.5 seconds at the end because the subjects show no signs of pain during
this duration. The total number of frames is 75. The extracted faces are fed
to an R3D model for a visual feature extractor. The batch size is set to 64.
The network is separately optimized using the SGD optimizer and the learning
rate of $10^{-3}$.

Table 2: Description of the architecture of the custom 1D-CNN for the physiological backbone. Layer | No. of | Size of |  |
---|---|---|---|---
type | Filters | Kernel | Stride | Output
Input | - | - | - |  $2816$ $\times$ $1$
1st Conv | $32$ | $5$ | $2$ |  $1406$ $\times$ $32$
ReLU | - | - | - |  $1406$ $\times$ $32$
Max Pooling | - | $2$ | - |  $703$ $\times$ $32$
2nd Conv | $64$ | $5$ | $1$ |  $699$ $\times$ $64$
ReLU | - | - | - |  $699$ $\times$ $64$
Max Pooling | - | $2$ | - |  $349$ $\times$ $64$
1st FC | - | - | - | $512$
ReLU | - | - | - | $512$
2nd FC | - | - | - | $2$

For the physiological modality, the EDA is used. The signal is clipped to
correspond to the visual modality and fed to a 1D CNN. The architecture of the
custom 1D-CNN is shown in Table 2. The CNN outputs a 512-dimensional feature
vector. The model is optimized using the SGD optimizer with a learning rate of
$10^{-4}$. The batch size for the physiological backbone is set to 1024.

Figure 3: Illustration of the proposed joint multimodal transformer
architecture used for the Biovid pain estimation task. The blue branch shows
the visual backbone, and the yellow branch is the physiological backbone. The
joint representation is shown with a red block. The three feature vectors are
fed into the joint transformer block.

For the modality fusion, the two feature vectors and the joint representation
are fed to the joint transformer block, as show in Figure 3. The FC layers are
removed from both backbones that were added in the backbone training phase.
512-dimensional feature vectors from visual and physiological backbones are
obtained and fed to the joint multimodal transformer module. The backbones are
frozen, and the joint transformer block is optimized using the ADAM optimizer
with a learning rate of $5\times 10^{-6}$, and the batch size is set to 128.

##  5 Results and Discussion

###  5.1 Comparison with the State-of-the-Art

Table 1 presents the performance of our fusion model on the Biovid Heat Pain
Database. We performed 5-fold cross-validation to pick up the best average
fusion model. It can be seen from the table that the proposed model can
achieve state-of-the-art performance with multimodal input while using 5-fold
cross-validation. The physiological modality is stronger in the Biovid
database. Many studies have validated the models on physiological modality.
Our empirical results show that the EDA-only accuracy is 77.2%, whereas the
visual-only accuracy is 72.9%. The proposed model can improve over unimodal
performance and achieves state-of-the-art performance on the Biovid dataset.
We also compare it with standard fusion techniques like feature concatenation.
For a fair comparison, we keep all the parameters the same. The proposed model
improves 6% over simple feature concatenation and 1.3% over a vanilla
multimodal transformer i.e. without a joint representation. Figure 4 shows the
visualization of the attention weights generated by the joint transformer
model. It can be seen that the model gives more weightage to the physiological
modality.

Figure 4:  Visualization of attention weights for visual and physiological
modalities on the Biovid heat pain database. The facial frames are taken 1400
msec each.

Table 3 shows the valence and arousal CCC values on the Affwild2 official
validation set and custom-defined folds to increase generalizability. On the
official validation set, the proposed method achieves a 0.666 average with
0.717 valence and 0.614 arousal.

Table 3: CCC for valence and arousal of the fusion model trained on different folds of the Affwild2 validation set. Highest scores are indicated in bold. Validation Set | Valence | Arousal | Average
---|---|---|---
Official | 0.717 | 0.614 | 0.666
fold-1 | 0.705 | 0.683 | 0.694
fold-2 | 0.741 | 0.623 | 0.682
fold-3 | 0.657 | 0.637 | 0.647
fold-4 | 0.760 | 0.666 | 0.713
fold-5 | 0.684 | 0.629 | 0.657

Table 4 compares the proposed method with the state-of-the-art on Affwild2
development set. Our method achieves a 0.722 average with 0.76 valence and
0.68 arousal.

Table 4: CCC performance of the proposed and state-of-the-art methods for A-V fusion on the Affwild2 development set. Method | Valence | Arousal | Mean
---|---|---|---
Joint Multimodal Transformer (ours) | 0.760 | 0.684 | 0.722
Situ-RUCAIM3 [37] | 0.655 | 0.709 | 0.682
Joint Cross- Attention [46] | 0.670 | 0.590 | 0.630
FlyingPigs [62] | 0.559 | 0.671 | 0.615
PRL [41] | 0.390 | 0.540 | 0.465
AU-NO [14] | 0.361 | 0.551 | 0.456
HSE-NN [49] | 0.408 | 0.477 | 0.443

Table 5 compares the proposed method with the state-of-the-art on the Affwild2
test set. The proposed method can achieve 0.42 valence and 0.46 arousal. The
average of valence and arousal is 0.443. The proposed method significantly
improves over the baseline, which is 0.17 average. It is essential to mention
here that the other methods that achieve higher performance are due to
extensive pertaining, the use of additional modalities like text, and the use
of more robust backbones. For a fair comparison, we use a similar setting to
Praveen et al. [46], which includes similar pretraining and identical
backbones for audio and visual modalities. They used joint cross-attention to
fuse the two modalities and achieve a 0.369 average. On the other hand, the
proposed model uses a joint multimodal transformer-based fusion method and can
achieve a 0.443 average. The proposed method improves 7% over the joint cross-
attention-based method.

Table 5: CCC performance of the proposed method and state-of-the-art methods for A-V fusion on the Affwild2 test set. Method | Valence | Arousal | Mean
---|---|---|---
Situ-RUCAIM3 [37] | 0.606 | 0.596 | 0.601
FlyingPigs [62] | 0.520 | 0.602 | 0.561
PRL [41] | 0.450 | 0.445 | 0.448
Joint Multimodal Transformer (ours) | 0.420 | 0.467 | 0.443
HSE-NN [49] | 0.417 | 0.454 | 0.436
AU-NO [14] | 0.418 | 0.407 | 0.413
Joint Cross- Attention [46] | 0.374 | 0.363 | 0.369
Baseline | 0.180 | 0.170 | 0.175

###  5.2 Ablations

To see the effectiveness of the joint representation of the proposed method,
we performed an ablation study. Table 6 shows the results of the proposed
method with and without the joint representation. On the Biovid dataset, the
joint multimodal transformer improves by 1.3% over the vanilla multimodal
transformer.

Table 6: Performance of the proposed approach compared to the vanilla multimodal transformer on test set. Valence and arousal for the Affwild2 dataset and accuracy for the Biovid dataset. I3D + ResNet18 backbones are used for the Affwild2 dataset and R3D + 1D CNN are used for the BioVid dataset. Default training/validation split is used in the Affwild2 dataset and 5-fold cross validation is performed on the BioVid dataset. Database | Method | Accuracy
---|---|---
Affwild2 |  Vanilla Multimodal Transformer |  V: 0.432 A: 0.410 Avg: 0.421
Joint Multimodal Transformer |  V: 0.425 A: 0.450 Avg: 0.438
Biovid |  Vanilla Multimodal Transformer | 87.8
Joint Multimodal Transformer | 89.1

Table 7 shows that the proposed model improves accuracy over the vanilla
multimodal transformer by 1.8% with a R2D1 and ResNet18 backbone on the
Affwild2 dataset. With an I3D and a ResNet18 backbone, results are improved by
1.7%, as shown in Table 8. Notice that the results on the test set are
different from Table 5 because we used a different train/validation split.

Table 7: CCC performance of the backbones alone, and combined using baselines and the proposed A-V fusion method on the Affwild2 test set. Experiments are performed on the default training-validation split. R2D1 (visual) and ResNet18 (audio) backbones are used for all cases. Fusion Model | Valence | Arousal | Mean
---|---|---|---
R2D1 model only (visual) | 0.194 | 0.310 | 0.252
ResNet18 model only (audio) | 0.273 | 0.246 | 0.260
Concat + FC layers | 0.320 | 0.327 | 0.323
Vanilla Transformer | 0.376 | 0.334 | 0.355
Joint Multimodal Transformer (ours) | 0.366 | 0.379 | 0.373
Table 8: CCC performance of the backbones alone, and combined using baselines and the proposed A-V fusion method on the Affwild2 test set. Experiments are performed on the default training-validation split. I3D (visual) and ResNet18 (audio) backbones are used for all cases. Fusion Model | Valence | Arousal | Mean
---|---|---|---
I3D model only (visual) | 0.336 | 0.422 | 0.379
ResNet18 model only (audio) | 0.273 | 0.246 | 0.260
Concat + FC layers | 0.387 | 0.453 | 0.420
Vanilla Transformer | 0.432 | 0.410 | 0.421
Joint Multimodal Transformer (ours) | 0.425 | 0.450 | 0.438
Joint Cross Attention [46] | 0.374 | 0.363 | 0.369

##  6 Conclusion

Multimodal emotion recognition systems outperform their unimodal counterparts,
especially in the wild environment. The missing and noisy modality is a
prevalent issue with in-the-wild emotion recognition systems. Many attention-
based methods have been proposed in the literature to overcome this problem.
These methods aim to weigh the modalities dynamically. This paper introduces a
joint multimodal transformer for dimensional emotional recognition in the
wild. This transformer-based architecture introduces a joint feature
representation in order to add more redundancy and complementary between audio
and visual data. The proposed model outperforms the state of the art on the
Biovid dataset and improves 6% on the Affwild2 dataset. In the future, we aim
to introduce more modalities and sophisticated backbones for effective feature
extraction.

## References

  * Anagnostopoulos et al., [2015] Anagnostopoulos, C., Iliou, T., and Giannoukos, I. (2015).  Features and classifiers for emotion recognition from speech: a survey from 2000 to 2011.  Artif Intell Rev, 43(6):155–177.
  * Aslam et al., [2024] Aslam, M. H., Zeeshan, M. O., Belharbi, S., Pedersoli, M., Koerich, A. L., Bacon, S., and Granger, E. (2024).  Distilling privileged multimodal information for expression recognition using optimal transport.  In 18th IEEE International Conference on Automatic Face and Gesture Recognition. arXiv preprint arXiv:2401.15489.
  * Aslam et al., [2023] Aslam, M. H., Zeeshan, O., Pedersoli, M., Koerich, A. L., Bacon, S., and Granger, E. (2023).  Privileged knowledge distillation for dimensional emotion recognition in the wild.  In CVPRw 2023: IEEE/CVF Conference on Computer Vision and Pattern Recognition, Vancouver, Canada.
  * Cho et al., [2014] Cho, K., Van Merriënboer, B., Gulcehre, C., Bahdanau, D., Bougares, F., Schwenk, H., and Bengio, Y. (2014).  Learning phrase representations using rnn encoder-decoder for statistical machine translation.  arXiv preprint arXiv:1406.1078.
  * de Melo et al., [2021] de Melo, W. C., Granger, E., and Lopez, M. B. (2021).  Mdn: A deep maximization-differentiation network for spatio-temporal depression detection.  IEEE Tran. on Affective Computing.
  * Deng et al., [2021] Deng, D., Wu, L., and Shi, B. E. (2021).  Iterative distillation for better uncertainty estimates in multitask emotion recognition.  In ICCVW, pages 3550–3559.
  * Dragomir et al., [2020] Dragomir, M.-C., Florea, C., and Pupezescu, V. (2020).  Automatic subject independent pain intensity estimation using a deep learning approach.  In 2020 International Conference on e-Health and Bioengineering (EHB), pages 1–4.
  * Duan et al., [2021] Duan, B., Tang, H., Wang, W., Zong, Z., Yang, G., and Yan, Y. (2021).  Audio-visual event localization via recursive fusion by joint co-attention.  In WACV, pages 4012–4021.
  * Ekman, [1992] Ekman, P. (1992).  An argument for basic emotions.  Cognition and Emotion, 6(3-4):169–200.
  * He et al., [2016] He, K., Zhang, X., Ren, S., and Sun, J. (2016).  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778.
  * Huang et al., [2020] Huang, J., Tao, J., Liu, B., Lian, Z., and Niu, M. (2020).  Multimodal transformer fusion for continuous emotion recognition.  In ICASSP 2020 - 2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 3507–3511.
  * Kächele et al., [2016] Kächele, M., Thiam, P., Amirian, M., Schwenker, F., and Palm, G. (2016).  Methods for person-centered continuous pain intensity assessment from bio-physiological channels.  IEEE Journal of Selected Topics in Signal Processing, 10(5):854–864.
  * Kächele et al., [2015] Kächele, M., Werner, P., Al-Hamadi, A., Palm, G., Walter, S., and Schwenker, F. (2015).  Bio-visual fusion for person-independent recognition of pain intensity.  In International Workshop on Multiple Classifier Systems.
  * Karas et al., [2022] Karas, V., Tellamekala, M. K., Mallol-Ragolta, A., Valstar, M., and Schuller, B. W. (2022).  Continuous-time audiovisual fusion with recurrence vs. attention for in-the-wild affect recognition.  arXiv preprint arXiv:2203.13285.
  * Kollias, [2022] Kollias, D. (2022).  Abaw: Valence-arousal estimation, expression recognition, action unit detection & multi-task learning challenges.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 2328–2336.
  * [16] Kollias, D. (2023a).  Abaw: learning from synthetic data & multi-task learning challenges.  In European Conference on Computer Vision, pages 157–172. Springer.
  * [17] Kollias, D. (2023b).  Multi-label compound expression recognition: C-expr database & network.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 5589–5598.
  * Kollias et al., [2020] Kollias, D., Schulc, A., Hajiyev, E., and Zafeiriou, S. (2020).  Analysing affective behavior in the first abaw competition.  In FG 2020.
  * [19] Kollias, D., Sharmanska, V., and Zafeiriou, S. (2019a).  Face behavior a la carte: Expressions, affect and action units in a single network.  arXiv preprint arXiv:1910.11111.
  * Kollias et al., [2021] Kollias, D., Sharmanska, V., and Zafeiriou, S. (2021).  Distribution matching for heterogeneous multi-task learning: a large-scale face study.  arXiv preprint arXiv:2105.03790.
  * Kollias et al., [2023] Kollias, D., Tzirakis, P., Baird, A., Cowen, A., and Zafeiriou, S. (2023).  Abaw: Valence-arousal estimation, expression recognition, action unit detection & emotional reaction intensity estimation challenges.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 5888–5897.
  * Kollias et al., [2024] Kollias, D., Tzirakis, P., Cowen, A., Zafeiriou, S., Shao, C., and Hu, G. (2024).  The 6th affective behavior analysis in-the-wild (abaw) competition.  arXiv preprint arXiv:2402.19344.
  * [23] Kollias, D., Tzirakis, P., Nicolaou, M. A., Papaioannou, A., Zhao, G., Schuller, B., Kotsia, I., and Zafeiriou, S. (2019b).  Deep affect prediction in-the-wild: Aff-wild database and challenge, deep architectures, and beyond.  IJCV, 127:907–929.
  * [24] Kollias, D., Tzirakis, P., Nicolaou, M. A., Papaioannou, A., Zhao, G., Schuller, B., Kotsia, I., and Zafeiriou, S. (2019c).  Deep affect prediction in-the-wild: Aff-wild database and challenge, deep architectures, and beyond.  International Journal of Computer Vision, pages 1–23.
  * Kollias and Zafeiriou, [2019] Kollias, D. and Zafeiriou, S. (2019).  Expression, affect, action unit recognition: Aff-wild2, multi-task learning and arcface.  arXiv preprint arXiv:1910.04855.
  * [26] Kollias, D. and Zafeiriou, S. (2021a).  Affect analysis in-the-wild: Valence-arousal, expressions, action units and a unified framework.  arXiv preprint arXiv:2103.15792.
  * [27] Kollias, D. and Zafeiriou, S. (2021b).  Analysing affective behavior in the second abaw2 competition.  In ICCVW, pages 3652–3660.
  * [28] Kuhnke, F., Rumberg, L., and Ostermann, J. (2020a).  Two-stream aural-visual affect analysis in the wild.  In FG Workshop, pages 600–605.
  * [29] Kuhnke, F., Rumberg, L., and Ostermann, J. (2020b).  Two-stream aural-visual affect analysis in the wild.  In 2020 15th IEEE International Conference on Automatic Face and Gesture Recognition (FG 2020), pages 600–605. IEEE.
  * Lawrence and Lin, [1989] Lawrence, I. and Lin, K. (1989).  A concordance correlation coefficient to evaluate reproducibility.  Biometrics, pages 255–268.
  * Le et al., [2023] Le, H.-D., Lee, G.-S., Kim, S.-H., Kim, S., and Yang, H.-J. (2023).  Multi-label multimodal emotion recognition with transformer-based fusion and emotion-level representation learning.  IEEE Access, 11:14742–14751.
  * Lee et al., [2021] Lee, J.-T., Jain, M., Park, H., and Yun, S. (2021).  Cross-attentional audio-visual fusion for weakly-supervised action localization.  In ICLR.
  * Lopez-Martinez and Picard, [2017] Lopez-Martinez, D. and Picard, R. (2017).  Multi-task neural networks for personalized pain recognition from physiological signals.  In 2017 Seventh International Conference on Affective Computing and Intelligent Interaction Workshops and Demos (ACIIW), pages 181–184. IEEE.
  * Lopez-Martinez and Picard, [2018] Lopez-Martinez, D. and Picard, R. (2018).  Continuous pain intensity estimation from autonomic signals with recurrent neural networks.  In 2018 40th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), pages 5624–5627. IEEE.
  * Lu et al., [2019] Lu, J., Batra, D., Parikh, D., and Lee, S. (2019).  Vilbert: Pretraining task-agnostic visiolinguistic representations for vision-and-language tasks.  Advances in neural information processing systems, 32.
  * Lu et al., [2023] Lu, Z., Ozek, B., and Kamarthi, S. (2023).  Transformer encoder with multiscale deep learning for pain classification using physiological signals.  arXiv preprint arXiv:2303.06845.
  * Meng et al., [2022] Meng, L., Liu, Y., Liu, X., Huang, Z., Jiang, W., Zhang, T., Deng, Y., Li, R., Wu, Y., Zhao, J., et al. (2022).  Multi-modal emotion estimation for in-the-wild videos.  arXiv preprint arXiv:2203.13032.
  * Morabit and Rivenq, [2022] Morabit, S. E. and Rivenq, A. (2022).  Pain detection from facial expressions based on transformers and distillation.  In 2022 11th International Symposium on Signal, Image, Video and Communications (ISIVC), pages 1–5.
  * Ngiam et al., [2011] Ngiam, J., Khosla, A., Kim, M., Nam, J., Lee, H., and Ng, A. Y. (2011).  Multimodal deep learning.  In Proceedings of the 28th International Conference on International Conference on Machine Learning, ICML’11, page 689–696, Madison, WI, USA. Omnipress.
  * Nguyen et al., [2021] Nguyen, D., Nguyen, D. T., Zeng, R., Nguyen, T. T., Tran, S., Nguyen, T. K., Sridharan, S., and Fookes, C. (2021).  Deep auto-encoders with sequential learning for multimodal dimensional emotion recognition.  IEEE Trans. on Multimedia, pages 1–1.
  * Nguyen et al., [2022] Nguyen, H.-H., Huynh, V.-T., and Kim, S.-H. (2022).  An ensemble approach for facial expression analysis in video.  arXiv preprint arXiv:2203.12891.
  * Ortega et al., [2019] Ortega, J. D. S., Cardinal, P., and Koerich, A. L. (2019).  Emotion recognition using fusion of audio and video features.  In SMC, pages 3847–3852.
  * Parthasarathy and Sundaram, [2020] Parthasarathy, S. and Sundaram, S. (2020).  Training strategies to handle missing modalities for audio-visual expression recognition.  In ICMI, page 400–404.
  * Phan et al., [2023] Phan, K. N., Iyortsuun, N. K., Pant, S., Yang, H.-J., and Kim, S.-H. (2023).  Pain recognition with physiological signals using multi-level context information.  IEEE Access, 11:20114–20127.
  * Pouromran et al., [2021] Pouromran, F., Radhakrishnan, S., and Kamarthi, S. (2021).  Exploration of physiological sensors, features, and machine learning models for pain intensity estimation.  Pone Journal, 16(7):e0254108.
  * [46] Praveen, R. G., de Melo, W. C., Ullah, N., Aslam, H., Zeeshan, O., Denorme, T., Pedersoli, M., Koerich, A. L., Bacon, S., Cardinal, P., et al. (2022a).  A joint cross-attention model for audio-visual fusion in dimensional emotion recognition.  In IEEE/CVF CVPR, pages 2486–2495.
  * [47] Praveen, R. G., de Melo, W. C., Ullah, N., Aslam, H., Zeeshan, O., Denorme, T., Pedersoli, M., Koerich, A. L., Bacon, S., Cardinal, P., and Granger, E. (2022b).  A joint cross-attention model for audio-visual fusion in dimensional emotion recognition.  In CVPRW, pages 2485–2494.
  * Rajasekhar et al., [2021] Rajasekhar, G. P., Granger, E., and Cardinal, P. (2021).  Cross attentional audio-visual fusion for dimensional emotion recognition.  In FG, pages 1–8.
  * Savchenko, [2022] Savchenko, A. V. (2022).  Frame-level prediction of facial expressions, valence, arousal and action units for mobile devices.  arXiv preprint arXiv:2203.13436.
  * Schoneveld et al., [2021] Schoneveld, L., Othmani, A., and Abdelkawy, H. (2021).  Leveraging recent advances in deep learning for audio-visual emotion recognition.  Pattern Recognition Letters, 146:1–7.
  * Thiam et al., [2019] Thiam, P., Bellmann, P., Kestler, H. A., and Schwenker, F. (2019).  Exploring deep physiological models for nociceptive pain recognition.  Sensors, 19(20):4503.
  * Thiam et al., [2021] Thiam, P., Hihn, H., Braun, D. A., Kestler, H. A., and Schwenker, F. (2021).  Multi-modal pain intensity assessment based on physiological signals: A deep learning perspective.  Frontiers in Physiology, 12:720464.
  * Tran and Soleymani, [2022] Tran, M. and Soleymani, M. (2022).  A pre-trained audio-visual transformer for emotion recognition.  In ICASSP 2022 - 2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 4698–4702.
  * Tzirakis et al., [2021] Tzirakis, P., Chen, J., Zafeiriou, S., and Schuller, B. (2021).  End-to-end multimodal affect recognition in real-world environments.  Information Fusion, 68:46–53.
  * Tzirakis et al., [2017] Tzirakis, P., Trigeorgis, G., Nicolaou, M. A., Schuller, B. W., and Zafeiriou, S. (2017).  End-to-end multimodal emotion recognition using deep neural networks.  IEEE J. of Selected Topics in Signal Processing, 11(8):1301–1309.
  * Wang et al., [2020] Wang, R., Xu, K., Feng, H., and Chen, W. (2020).  Hybrid rnn-ann based deep physiological network for pain recognition.  In 2020 42nd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC), pages 5584–5587. IEEE.
  * Wei et al., [2020] Wei, X., Zhang, T., Li, Y., Zhang, Y., and Wu, F. (2020).  Multi-modality cross attention network for image and sentence matching.  In 2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 10938–10947.
  * Werner et al., [2016] Werner, P., Al-Hamadi, A., Limbrecht-Ecklundt, K., Walter, S., Gruss, S., and Traue, H. C. (2016).  Automatic pain assessment with facial activity descriptors.  IEEE Transactions on Affective Computing, 8(3):286–299.
  * [59] Werner, P., Al-Hamadi, A., Niese, R., Walter, S., Gruss, S., and Traue, H. C. (2014a).  Automatic pain recognition from video and biomedical signals.  In 2014 22nd International Conference on Pattern Recognition, pages 4582–4587.
  * [60] Werner, P., Al-Hamadi, A., Niese, R., Walter, S., Gruss, S., and Traue, H. C. (2014b).  Automatic pain recognition from video and biomedical signals.  In 2014 22nd international conference on pattern recognition, pages 4582–4587. IEEE.
  * Zafeiriou et al., [2017] Zafeiriou, S., Kollias, D., Nicolaou, M. A., Papaioannou, A., Zhao, G., and Kotsia, I. (2017).  Aff-wild: Valence and arousal ‘in-the-wild’challenge.  In Computer Vision and Pattern Recognition Workshops (CVPRW), 2017 IEEE Conference on, pages 1980–1987. IEEE.
  * Zhang et al., [2022] Zhang, S., An, R., Ding, Y., and Guan, C. (2022).  Continuous emotion recognition using visual-audio-linguistic information: A technical report for abaw3.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 2376–2381.
  * Zhang et al., [2021] Zhang, S., Ding, Y., Wei, Z., and Guan, C. (2021).  Continuous emotion recognition with audio-visual leader-follower attentive fusion.  In ICCVW, pages 3560–3567.
  * Zhang et al., [2020] Zhang, Y.-H., Huang, R., Zeng, J., and Shan, S. (2020).  Multi-modal continuous valence-arousal estimation in the wild.  In IEEE FG, pages 632–636.
  * Zhi et al., [2021] Zhi, R., Zhou, C., Yu, J., Li, T., and Zamzmi, G. (2021).  Multimodal-based stream integrated neural networks for pain assessment.  IEICE Trans. Inf. Syst., 104-D:2184–2194.
  * Zhou et al., [2023] Zhou, W., Lu, J., Xiong, Z., and Wang, W. (2023).  Leveraging tcn and transformer for effective visual-audio fusion in continuous emotion recognition.  In 2023 IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW), pages 5756–5763.


# Half-Precision Kronecker Product SVD Preconditioner
for Structured Inverse Problems

Yizhou Chen Department of Mathematics, Emory University, Atlanta, GA 30322,
USA. E-mail<EMAIL_ADDRESS>Xiang Ji Department of Mathematics, Emory
University, Atlanta, GA 30322, USA. Corresponding author.
E-mail<EMAIL_ADDRESS>James Nagy Department of Mathematics, Emory
University, Atlanta, GA 30322, USA. E-mail<EMAIL_ADDRESS>Research
supported in part by the U.S. National Science Foundation, grant codes
DMS-2038118 and DMS-2208294.

###### Abstract

In this paper we investigate the use of half-precision Kronecker product
singular value decomposition (SVD) approximations as preconditioners for
large-scale Tikhonov regularized least squares problems. Half precision
reduces storage requirements and has the potential to greatly speedup
computations on certain GPU architectures. We consider both standard PCG and
flexible PCG algorithms, and investigate, through numerical experiments on
image deblurring problems, the trade-offs between potentially faster
convergence with the additional cost per iteration when using this
preconditioning approach. Moreover, we also investigate the use of several
regularization parameter choice methods, including generalized cross
validation and the discrepancy principle.

##  1 Introduction

Inverse problems can occur in various fields, when observations are acquired
to infer the original data. They can be described as the following linear
system:

${\bm{b}}=A{\bm{x}}_{{}_{\mathrm{true}}}+{\bm{\varepsilon}}$ (1)

where ${\bm{b}}$ is the outside measurements containing some noise
${\bm{\varepsilon}}$ and $A$, usually ill-conditioned and large-scale,
represents some mapping of the true solution ${\bm{x}}_{{}_{\mathrm{true}}}$
[1]. Since the noise is unknown, we can only calculate an approximate solution
${\bm{x}}$, instead of ${\bm{x}}_{{}_{\mathrm{true}}}$. A naïve approach
computes an approximation by ignoring the noise, and either solving
$A{\bm{x}}={\bm{b}}$, if $A$ is $n\times n$ and non-singular, or solving the
least squares problem,

$\min_{\bm{x}}\|{\bm{b}}-A{\bm{x}}\|^{2}_{2}\,.$ (2)

In either case, when $A$ is large, iterative methods are needed. An additional
difficulty is that $A$ is typically very ill-conditioned, with singular values
clustering at $0$. In this case, noise in the measured data is greatly
amplified when computing the naïve solution, resulting in very poor
approximations of ${\bm{x}}_{{}_{\mathrm{true}}}$. To get around this
difficulty, regularization methods are usually incorporated into the solution
process.

One of the most common approaches is Tikhonov regularization, which is
implemented to balance the smoothness of the solution and the residual term.
It can be written as the least squares problem,

$\min_{\bm{x}}\|{\bm{b}}-A{\bm{x}}\|^{2}_{2}+\lambda^{2}\|{\bm{x}}\|^{2}_{2}=\min_{\bm{x}}\left\|\left[\begin{array}[]{c}{\bm{b}}\\\
{\bm{0}}\\\ \end{array}\right]-\left[\begin{array}[]{c}A\\\ \lambda I\\\
\end{array}\right]{\bm{x}}\right\|^{2}_{2}\,.$ (3)

The regularization parameter, $\lambda$, balances the residual norm with the
solution norm [8]; large $\lambda$ favors small norm (smoother) solutions,
while small $\lambda$ favors large residual norms (and hence, less smooth
solutions). There are various ways to choose the regularization parameter
$\lambda$, such as generalized cross validation (GCV), weighted GCV, and the
discrepancy principle [5, 9].

The conjugate gradient (CG) algorithm is a well-known and efficient iterative
method for solving large-scale sparse (or structured) linear systems. For
least squares problems, CG can be applied to the normal equations, which for
the Tikhonov regularization problem (3)takes the form

$\left[\begin{array}[]{cc}A^{T}&\lambda I\\\
\end{array}\right]\left[\begin{array}[]{c}A\\\ \lambda I\\\
\end{array}\right]{\bm{x}}=\left[\begin{array}[]{cc}A^{T}&\lambda I\\\
\end{array}\right]\left[\begin{array}[]{c}{\bm{b}}\\\ {\bm{0}}\\\
\end{array}\right]\,.$ $(A^{T}A+\lambda^{2}I){\bm{x}}=A^{T}{\bm{b}}\,.$ (4)

Normally, it is recommended to use either CGLS or LSQR implementations [3].
The overall cost of this algorithm depends on the number of iterations and the
cost of matrix-vector multiplications involving $A$ and $A^{T}$, which can be
done efficiently if $A$ is sparse or if it has an exploitable structure (e.g.,
Toeplitz or Kronecker product).

Preconditioning can be used in the CG algorithm (PCG) to accelerate
convergence. The basic idea, for equation (4), is to construct a matrix $M$
that has three properties. The first is that $M\approx A^{T}A+\lambda^{2}I$,
such that the eigenvalues of $M^{-1}(A^{T}A+\lambda^{2}I)$ are clustered
around a point bounded away from zero. The second is that solving the linear
system $M{\bm{z}}={\bm{r}}$ is not too expensive, which is required at each
iteration of the PCG algorithm. The third property is that it is not too
expensive to construct the approximating matrix $M$.

There exists various ways to construct $M$. In this paper we use Kronecker
product singular value decomposition (SVD) approximations. To further
accelerate the process, we use mixed-precision computations in the PCG
algorithm. Specifically, we use 16 bit half-precision arithmetic to store $M$
and to compute solutions of $M{\bm{z}}={\bm{r}}$, and 64 bit double precision
for the rest of the computations. Half precision reduces storage requirements
and has the potential to greatly speedup computations on certain GPU
architectures [14]. We consider both standard PCG and flexible PCG algorithms,
and investigate the trade-offs between potentially faster convergence with the
additional cost per iteration when using this preconditioning approach.
Numerical experiments to illustrate performance are reported for deconovlution
and image deblurring problems. To have full control over numerical precision,
all half precision computations are simulated in MATLAB using Higham’s chop
function [11].

We begin by providing some background and establish notation in Section 2. The
basic idea of our preconditioning approach is described in Section 3, and we
discuss some issues related to low-precision arithmetic simulations in Section
4. Numerical results are presented in Section 5, and we summarize with
concluding remarks in Section 6.

##  2 Background and Notation

This section is intended to establish notation used throughout the paper, as
well as to provide background material on inverse problems, regularization and
preconditioning.

###  2.1 Tikhonov Regularization

Because $A$ is severely ill-conditioned, ignoring the noise (even if it is
small) in (1) and naïvely solving $A{\bm{x}}={\bm{b}}$ (or
$\min\||{\bm{b}}-A{\bm{x}}\|_{2}$) typically produces very poor approximations
of ${\bm{x}}_{{}_{\mathrm{true}}}$. This can be seen from a singular value
decomposition (SVD).

Specifically, suppose111To simplify notation, we assume $A$ is an $n\times n$
matrix, but everything we discuss in the paper is applicable to $m\times n$
matrices. $A=U\Sigma V^{T}$ is the SVD of $A\in\mathds{R}^{n\times n}$, where
$U$ and $V$ are orthogonal matrices, and
$\Sigma=\,$diag$(\sigma_{1},\sigma_{2},\ldots,\sigma_{n})$ is a diagonal
matrix, with the singular values satisfying
$\sigma_{1}\geq\sigma_{2}\geq\cdots\geq\sigma_{n}\geq 0$. For inverse
problems, the singular values of $A$ decay rapidly toward zero, so even if we
assume $A$ is nonsingular and compute the inverse solution
${\bm{x}}_{{}_{\mathrm{inv}}}=A^{-1}{\bm{b}}$, we obtain

${\bm{x}}_{{}_{\mathrm{inv}}}=A^{-1}(A{\bm{x}}_{{}_{\mathrm{true}}}+{\bm{\varepsilon}})={\bm{x}}_{{}_{\mathrm{true}}}+A^{-1}{\bm{\varepsilon}}={\bm{x}}_{{}_{\mathrm{true}}}+\sum_{i=1}^{n}\frac{{\bm{u}}_{i}^{T}{\bm{\varepsilon}}}{\sigma_{i}}{\bm{v}}_{i}\,,$

where ${\bm{u}}_{i}$ is the $i$-th right singular vector of $A$ (i.e., $i$-th
column of $U$), and ${\bm{v}}_{i}$ is the $i$-th left singular vector of $A$
(i.e., $i$-th column of $V$). That is,
$\displaystyle{\bm{x}}_{{}_{\mathrm{inv}}}={\bm{x}}_{{}_{\mathrm{true}}}+$
error, and we observe that noise components
${\bm{u}}_{i}^{T}{\bm{\varepsilon}}$ can be greatly magnified when divided by
tiny $\sigma_{i}$ in the error term,
$\displaystyle\sum_{i=1}^{n}\frac{{\bm{u}}_{i}^{T}{\bm{\varepsilon}}}{\sigma_{i}}{\bm{v}}_{i}$.
Another important property of inverse problems is that the singular vectors
${\bm{v}}_{i}$ corresponding to large singular values tend to be smooth (low
frequency), and ${\bm{v}}_{i}$ corresponding to small singular values tend to
highly oscillatory (high frequency). Thus, division by tiny singular values,
$\sigma_{i}$, magnifies the high frequency components in the error,
${\bm{v}}_{i}$, and these high frequency components dominate any information
about ${\bm{x}}_{{}_{\mathrm{true}}}$ in the inverse solution
${\bm{x}}_{{}_{\mathrm{inv}}}$.

There are many regularization approaches that can be used to mitigate noise
magnification [4, 6, 9, 15, 21]. In this paper we use Tikhonov regularization,
which can be written as

$\min_{\bm{x}}\|{\bm{b}}-A{\bm{x}}\|^{2}_{2}+\lambda^{2}\|{\bm{x}}\|^{2}_{2}=\min_{\bm{x}}\left\|\left[\begin{array}[]{c}{\bm{b}}\\\
{\bm{0}}\\\ \end{array}\right]-\left[\begin{array}[]{c}A\\\ \lambda I\\\
\end{array}\right]{\bm{x}}\right\|^{2}_{2}\,.$ (5)

The regularization parameter $\lambda$ needs to be chosen to balance fit-to-
data ($\|{\bm{b}}-A{\bm{x}}\|_{2}^{2}$) with smoothness of the solution
($\|{\bm{x}}\|_{2}^{2}$); if $\lambda$ is tiny, then we get good fit-to-data
with a solution that is highly magnified by noise, and if we choose $\lambda$
too large, then we get a smooth solution that does not provide a very good
fit-to-data. It is also interesting to see that if we replace $A$ with its
SVD, and compute the least squares solution of (5), we obtain

${\bm{x}}=\sum_{i=1}^{n}\phi_{i}\frac{{\bm{u}}_{i}^{T}{\bm{b}}}{\sigma_{i}}{\bm{v}}_{i}=\sum_{i=1}^{n}\phi_{i}\frac{{\bm{u}}_{i}^{T}{\bm{b}}_{{}_{\mathrm{true}}}}{\sigma_{i}}{\bm{v}}_{i}+\sum_{i=1}^{n}\phi_{i}\frac{{\bm{u}}_{i}^{T}{\bm{\varepsilon}}}{\sigma_{i}}{\bm{v}}_{i}\,,$
(6)

where
$\displaystyle\phi_{i}=\frac{\sigma_{i}^{2}}{\sigma_{i}^{2}+\lambda^{2}}$ acts
as a filter; when $\sigma_{i}$ is large compared to $\lambda$,
$\phi_{i}\approx 1$ allowing low frequency components of the inverse solution
(those corresponding to large singular values) to be reconstructed, but when
$\sigma_{i}$ is small compared to $\lambda$, $\phi_{i}\approx 0$ and the high
frequency components in the error term (those corresponding to small singular
values) are suppressed.

###  2.2 Finding Regularization Parameters

Choosing the regularization parameter $\lambda$ is crucial to obtaining a good
approximate solution. In this subsection we review a few approaches that can
be used, provided we are able to compute the SVD, $A=U\Sigma V^{T}$.

For experimental purposes, it can be helpful to use test problems where the
true solution, ${\bm{x}}_{{}_{\mathrm{true}}}$ is known, and find a parameter
that minimizes the relative error; that is, find

$\lambda_{{}_{\mathrm{opt}}}=\arg\min_{\lambda}\frac{\|{\bm{x}}_{\lambda}-{\bm{x}}_{{}_{\mathrm{true}}}\|_{2}}{\|{\bm{x}}_{{}_{\mathrm{true}}}\|_{2}}\,,\quad\mbox{where}\quad{\bm{x}}_{\lambda}=\arg\min_{\bm{x}}\left\|\left[\begin{array}[]{c}{\bm{b}}\\\
{\bm{0}}\end{array}\right]-\left[\begin{array}[]{c}A\\\ \lambda
I\end{array}\right]{\bm{x}}\right\|^{2}_{2}\,.$ (7)

Substituting $A=U\Sigma V^{T}$ into (7), it is not difficult to show that
$\lambda_{{}_{\mathrm{opt}}}$ can be found by minimizing a one-dimensional
function,

$\lambda_{{}_{\mathrm{opt}}}=\arg\min_{\lambda}{\cal
O}(\lambda)=\arg\min_{\lambda}\sum_{i=1}^{n}\left(\frac{\sigma_{i}\hat{b}_{i}}{\sigma_{i}^{2}+\lambda^{2}}-\hat{x}_{i}\right)^{2}\,,$
(8)

where $\hat{\bm{b}}=U^{T}{\bm{b}}$ and
$\hat{\bm{x}}=V^{T}{\bm{x}}_{{}_{\mathrm{true}}}$. This minimization problem
can be found in MATLAB using, for example, the built-in function fminbnd,
assuming $\sigma_{n}\leq\lambda\leq\sigma_{1}$. We remark that this approach
is not practical for real problems, but it can be a good reference point (or
benchmark) in numerical tests, when comparing with other approaches that do
not require knowledge of ${\bm{x}}_{{}_{\mathrm{true}}}$.

The generalized cross validation (GCV) method is a popular approach that does
not require ${\bm{x}}_{{}_{\mathrm{true}}}$. The goal of GCV is to predict any
missing data from the data set. This means that the difference between an
arbitrary $b_{i}$ left-out of the original data set and the corresponding
solution filtered by $\lambda$ should be minimized [5]. This minimization is
carried out for all data points in ${\bm{b}}$. Although it is not obvious from
this description, it can be shown (see, e.g., [9]) that an estimate of the
regularization parameter using the GCV approach can be found by minimizing the
one-dimensional function

${\cal
G}(\lambda)=\frac{n\|(I-AA_{\lambda}^{\dagger}){\bm{b}}\|_{2}^{2}}{\left(\text{trace}(I-AA_{\lambda}^{\dagger})\right)^{2}}\,.$
(9)

Here $A_{\lambda}^{\dagger}=(A^{T}A+\lambda^{2}I)^{-1}A^{T}$ and trace means
summing up the diagonal elements of a matrix. Using $A=U\Sigma V^{T}$, it can
be shown [9] that the GCV chosen regularization parameter can be found by the
minimization

$\lambda_{{}_{\mathrm{gcv}}}=\arg\min_{\lambda}{\cal
G}(\lambda)=\arg\min_{\lambda}\frac{n\displaystyle\sum_{i=1}^{n}\left(\frac{\hat{b}_{i}}{\sigma_{i}^{2}+\lambda^{2}}\right)^{2}}{\left(\displaystyle\sum_{i=1}^{n}\frac{1}{\sigma_{i}^{2}+\lambda^{2}}\right)^{2}}\,.$
(10)

Note that although equation (10) appears to be complicated, it is a simple
one-dimensional minimization problem for $\lambda$, which can be solved, for
example, using MATLAB’s fminbnd function.

According to [5], GCV can sometimes choose a $\lambda$ that is too small, and
that a weighted version of GCV (wGCV) can overcome this problem:

${\cal
G}_{\omega}(\lambda)=\frac{n\|(I-AA_{\lambda}^{\dagger}){\bm{b}}\|_{2}^{2}}{\left(\text{trace}(I-\omega
AA_{\lambda}^{\dagger})\right)^{2}}\,.$ (11)

Using $A=U\Sigma V^{T}$, the wGCV chosen regularization parameter (for
Tikhonov regularization) can be found by the minimization

$\lambda_{{}_{\mathrm{wgcv}}}=\arg\min_{\lambda}{\cal
G}_{\omega}(\lambda)=\arg\min_{\lambda}\frac{n\displaystyle\sum_{i=1}^{n}\left(\frac{\lambda^{2}\hat{b}_{i}}{\sigma_{i}^{2}+\lambda^{2}}\right)^{2}}{\left(\displaystyle\sum_{i=1}^{n}\frac{(1-\omega)\sigma_{i}^{2}+\lambda^{2}}{\sigma_{i}^{2}+\lambda^{2}}\right)^{2}}\,.$
(12)

When the weight parameter $\omega=1$ we get the standard GCV function,
$\omega>1$ gives us smoother solutions, and solutions for $\omega<1$ are less
smooth.

Another well-known method for estimating regularization parameters is the
discrepancy principle. The main idea is to set an upper bound for the
residuals $\|A{\bm{x}}_{\lambda}-{\bm{b}}\|_{2}$, where ${\bm{x}}_{\lambda}$
is the solution computed with regularization parameter $\lambda$. That is, we
look for a regularization parameter $\lambda$ such that

${\cal D}(\lambda)=\|A{\bm{x}}_{\lambda}-{\bm{b}}\|_{2}^{2}-\epsilon^{2}=0\,,$
(13)

where $\epsilon=\eta\|{\bm{\varepsilon}}\|$, and $\eta\approx 1$. As with the
previous strategies for estimating regularization parameters, if we use
$A=U\Sigma V^{T}$ then the function ${\cal D}(\lambda)$ can be simplified.
Specifically, the discrepancy principle reduces to a root finding problem,
where we find $\lambda=\lambda_{{}_{\mathrm{disc}}}$ such that

${\cal
D}(\lambda)=\sum_{i=1}^{n}\left(\frac{\lambda^{2}\hat{b}_{i}}{\sigma_{i}^{2}+\lambda^{2}}\right)^{2}-\epsilon^{2}=0\,.$
(14)

The built-in MATLAB function fzero can be used to find
$\lambda_{{}_{\mathrm{disc}}}$.

We remark that to use the discrepancy principle, it is necessary to have a
good estimate of the noise, $\epsilon=\|{\bm{\varepsilon}}\|_{2}$.

###  2.3 Iterative approaches for Tikhonov regularization

For large-scale problems, it may not be computationally feasible to compute an
SVD unless $A$ has an exploitable structure. For example, if $A$ is a
circulant matrix, then fast Fourier transforms (FFTs) can be used to compute a
spectral decomposition, and if $A$ is a block circulant matrix with circulant
blocks, then 2-dimensional FFTs can be used. Another example occurs if $A$ can
be written as a Kronecker product of two smaller matrices; we discuss this in
more detail in Section 3.

If it is not computationally feasible to compute an SVD of $A$, then iterative
methods, such as the conjugate gradient (CG) algorithm, are often recommended.
For least squares problems, CG can be applied to the normal equations, which
for the Tikhonov regularization problem (3)takes the form

$\left[\begin{array}[]{cc}A^{T}&\lambda I\\\
\end{array}\right]\left[\begin{array}[]{c}A\\\ \lambda I\\\
\end{array}\right]\mathbf{x}=\left[\begin{array}[]{cc}A^{T}&\lambda I\\\
\end{array}\right]\left[\begin{array}[]{c}\mathbf{b}\\\ \mathbf{0}\\\
\end{array}\right]\,.$ $(A^{T}A+\lambda^{2}I)\mathbf{x}=A^{T}\mathbf{b}\,.$
(15)

Normally, it is recommended to use either CGLS or LSQR implementations [3].
The overall cost of this algorithm depends on the number of iterations and the
cost of matrix-vector multiplications involving $A$ and $A^{T}$, which can be
done efficiently if $A$ is sparse or if it has an exploitable structure (e.g.,
Toeplitz or a sum of Kronecker products).

Preconditioning can be used in the CG algorithm (PCG) to accelerate
convergence. The basic idea, for equation (15), is to construct a matrix $M$
that has three properties. The first is that $M\approx A^{T}A+\lambda^{2}I$,
such that the eigenvalues of $M^{-1}(A^{T}A+\lambda^{2}I)$ are clustered
around a point bounded away from zero (the more eigenvalues that are
clustered, and tighter clusters, leads to faster convergence). The second is
that solving the linear system $M\mathbf{z}=\mathbf{r}$ is not too expensive,
because it is required at each iteration of the PCG algorithm. The third
property is that it is not too expensive to construct the approximating matrix
$M$.

In this paper we consider preconditioners constructed from Kronecker product
SVD approximations of $A$. We discuss Kroneker products later, but for now it
is helpful to understand the basic idea. Suppose were are able to efficiently
compute the decomposition
$\widehat{A}=\widehat{U}\widehat{\Sigma}\widehat{V}^{T}\approx A$, and define

$\displaystyle M$ $\displaystyle=\widehat{A}^{T}\widehat{A}+\lambda^{2}I$
$\displaystyle=\widehat{V}\widehat{\Sigma}^{T}\widehat{U}^{T}\widehat{U}\widehat{\Sigma}\widehat{V}^{T}+\lambda^{2}I$
$\displaystyle=\widehat{V}\widehat{\Sigma}^{T}\widehat{\Sigma}\widehat{V}^{T}+\lambda^{2}I$
$\displaystyle=\widehat{V}(\widehat{\Sigma}^{T}\widehat{\Sigma}+\lambda^{2}I)\widehat{V}^{T}\,.$

If we can compute such an approximation, then the solution of the linear
system $M{\bm{z}}={\bm{r}}$ can be computed as

${\bm{z}}=M^{-1}{\bm{r}}=\widehat{V}\big{(}D^{-1}\big{(}\widehat{V}^{T}{\bm{r}}\big{)}\big{)}\,,\quad\mbox{where}\quad
D=\widehat{\Sigma}^{T}\widehat{\Sigma}+\lambda^{2}I\,.$ (16)

That is, solving $M{\bm{z}}={\bm{r}}$ requires matrix-vector multiplication
with $\widehat{V}$ and $\widehat{V}^{T}$, and multiplication by the diagonal
matrix $D^{-1}$.

##  3 Kronecker Product SVD Approximations for Preconditioning

In this section we discuss the details of the low precision, Kronecker product
SVD approximation preconditioner we propose to use in this work.

We first begin by supposing $A\in\mathds{R}^{N\times N}$ where $N=n^{2}$ can
be approximated by a Kronecker product

$A\approx\widehat{A}=A_{r}\otimes
A_{c}=\left[\begin{array}[]{ccc}a^{{}^{(r)}}_{11}A_{c}&\cdots&a^{{}^{(r)}}_{1n}A_{c}\\\\[3.0pt]
\vdots&&\vdots\\\\[6.0pt]
a^{{}^{(r)}}_{n1}A_{c}&\cdots&a^{{}^{(r)}}_{nn}A_{c}\end{array}\right]$

where $A_{r}$ and $A_{c}$ are $n\times n$ matrices, and $a^{{}^{(r)}}_{ij}$ is
the $(i,j)$ entry of $A_{r}$. There are many convenient properties of
Kronecker products that can be exploited for efficient computations [19].
Specifically, for general matrices $C,D,E,F\in\mathds{R}^{n\times n}$ and
${\bm{y}}\in\mathds{R}^{n^{2}}$,

  * •

$(C\otimes D)(E\otimes F)=CE\otimes DF$.

  * •

$(C\otimes D)^{T}=C^{T}\otimes D^{T}$.

  * •

If $C$ and $D$ are nonsingular, then so is $C\otimes D$, and $(C\otimes
D)^{-1}=C^{-1}\otimes D^{-1}$.

  * •

${\bm{x}}=(C\otimes D){\bm{y}}={\tt vec}(DYC^{T})$, where

${\bm{y}}=\mbox{vec}(Y)=\mbox{vec}\left(\left[\begin{array}[]{cccc}{\bm{y}}_{1}&{\bm{y}}_{2}&\cdots&{\bm{y}}_{n}\end{array}\right]\right)=\left[\begin{array}[]{c}{\bm{y}}_{1}\\\
{\bm{y}}_{2}\\\ \vdots\\\
{\bm{y}}_{n}\end{array}\right]\,,\quad{\bm{y}}_{i}\in\mathds{R}^{n}\,.$

Using these properties for $A\approx A_{r}\otimes A_{c}$, with
$A_{r}=U_{r}\Sigma_{r}V^{T}_{r}$ and $A_{c}=U_{c}\Sigma_{c}V^{T}_{c}$, then

$A\approx\widehat{A}=A_{r}\otimes
A_{c}=(U_{r}\Sigma_{r}V^{T}_{r})\otimes(U_{c}\Sigma_{c}V^{T}_{c})=(U_{r}\otimes
U_{c})(\Sigma_{r}\otimes\Sigma_{c})(V_{r}\otimes V_{c})^{T}\,.$ (17)

The matrices $U_{r}\otimes U_{c}$ and $V_{r}\otimes V_{c}$ do not need to be
formed explicitly, thus required storage for the decomposition is only $O(N)$.
(We do not mention anything about $\Sigma_{r}\otimes\Sigma_{c}$ because it
involves diagonal matrices, and so contributes only trivially to the
computational costs.) In addition, the decomposition and computations with
$U_{r}\otimes U_{c}$ and $V_{r}\otimes V_{c}$ require only $O(N^{3/2})$
arithmetic operations. This is very efficient compared to the standard
approaches applied directly to $A$, which require $O(N^{2})$ storage and
$O(N^{3})$ operations.

A general approach for constructing a Kronecker product approximation was
proposed by Van Loan and Pitsianis [20], and an alternative approach for
matrices that have a block Toeplitz with Toeplitz block (BTTB) structure,
which arise in image deblurring problems, was described in [12, 13, 16, 17].
An implementation to efficiently construct this approximation for image
deblurring problems, which we use for this work, is in the IR Tools MATLAB
software package [7].

To solve linear systems (16) with the Kronecker product SVD approximation
(17), we exploit the properties of Kronecker products. Specifically, the
vector ${\bm{z}}$ in equation (16) can be computed as follows,

$Z=V_{c}\Big{(}S\,\mathchoice{\mathbin{{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}*\;\big{(}V_{c}^{T}R\,V_{r}\big{)}\Big{)}V_{r}^{T}\,,$
(18)

where we use the MATLAB notation
$\,\mathchoice{\mathbin{{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}*\;$
for element-wise multiplication, $\mbox{vec}(S)=\mbox{diag}(D^{-1})$,
${\bm{r}}=\mbox{vec}(R)$, and ${\bm{z}}=\mbox{vec}(Z)$.

Another important aspect of our preconditioner is that we store the components
of $M$ (i.e., $V_{r}$, $V_{c}$ and $S$), and perform the computations in (18),
using 16-bit (half) precision arithmetic. Although we do not compute $S$,
$V_{r}$ and $V_{c}$ in half precision, this should be considered in future
work. We note that using half precision for preconditioner solves means that,
due to roundoff errors, the preconditioner does not remain constant at each
iteration, and could lead to a loss of orthogonality of the generated Krylov
search directions in PCG [2]. A workaround to this problem is to use a
flexible PCG algorithm [18], which essentially replaces how an update
parameter of the search direction (standard PCG uses Fletcher-Reeves) with an
alternative formula (Polak-Ribieŕe). This workaround requires an additional
inner product and additional storage [2].

##  4 Low Precision Simulation in MATLAB

Before proceeding to the numerical experiments, we discuss a few important
issues related to simulating low precision arithmetic in MATLAB. Higham and
Pranesh presented a MATLAB function named chop that can be used to simulate
lower precision arithmetic by rounding entries of vectors and matrices to the
target precision level [11]. It supports various precision formats, for
example IEEE half precision, or even user-customized formats. The input of the
function is in either double or single precision, and the output is a value at
the specified precision; the number is stored in higher precision type with
extra bits set to zero.

When simulating low-precision arithmetic, chop needs to be applied for each
floating point operation. For example, the scalar operation $x+y\times z$
requires two calls to chop:

${\rm chop}(x+{\rm chop}(y\times z))\,.$

When simulating low-precision computations for large-scale problems, there are
two important issues to consider:

  * •

Each call to chop is costly.

  * •

Proper implementations are needed to reduce buildup of roundoff error when
replacing built-in computatioanl kernels (e.g., norm) with user defined
functions.

Consider, for example, the simple computation of an inner product,
$s={\bm{x}}^{T}{\bm{y}}$. Since every floating point operation needs to be
properly rounded, one might consider the implementation:

  * s = 0;
for i = 1:n
s = chop(s + chop(x(i)*y(i)))
end

This requires $2n$ calls to chop. We can reduce this to $n+1$ calls by
exploiting vectorization:

  * z = chop(x$\,\mathchoice{\mathbin{{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}*\;$y);
s = 0;
for i = 1:n
s = chop(s + z)
end

However, it is well-known that the above approach to compute the sum of
floating point numbers can lead to a buildup of rounding errors. Alternative
summation approaches include various ordering schemes, blocking, Kahan’s
compensated summation, and pairwise summation [10]. In our implementations, we
use pairwise summation because, in addition to having good rounding error
properties, it also significantly reduces the number of calls to chop.

To visualize the pairwise summation approach to compute
$s={\bm{x}}^{T}{\bm{y}}$, first compute the element-wise product
${\bm{z}}={\tt
chop}({\bm{x}}\,\mathchoice{\mathbin{{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}*\;{\bm{y}})$,
so that only one vectorized call to chop is needed to compute all entries in
${\bm{z}}$. Pairwise summation then proceeds with a divide-and-conquer
approach to sum all of the entries in ${\bm{z}}$, two at a time. This can be
illustrated with a small example, ${\bm{z}}\in\mathds{R}^{8}$, noticing that
we can do pairwise summation as a vectorized operation, followed by a single
vectorized call to chop at each stage:

$\displaystyle{\bm{z}}^{T}={\bm{z}}^{(0)T}$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{cccccccc}{\bm{z}}^{(0)}_{1}&{\bm{z}}^{(0)}_{2}&{\bm{z}}^{(0)}_{3}&{\bm{z}}^{(0)}_{4}&{\bm{z}}^{(0)}_{5}&{\bm{z}}^{(0)}_{6}&{\bm{z}}^{(0)}_{7}&{\bm{z}}^{(0)}_{8}\end{array}\right]$
$\displaystyle\begin{array}[]{cccc}\;\underbrace{\hskip
42.0pt}_{\text{add}}&\;\underbrace{\hskip
42.0pt}_{\text{add}}&\;\underbrace{\hskip
42.0pt}_{\text{add}}&\;\underbrace{\hskip 42.0pt}_{\text{add}}\end{array}$
$\displaystyle\hskip 20.0pt\underbrace{\hskip 182.0pt}_{\text{vectorized
chop}}$ $\displaystyle=$ $\displaystyle\left[\begin{array}[]{cccc}\hskip
14.0pt{\bm{z}}^{(1)}_{1}&\hskip 14.0pt{\bm{z}}^{(1)}_{2}&\hskip
14.0pt{\bm{z}}^{(1)}_{3}&\hskip 14.0pt{\bm{z}}^{(1)}_{4}\end{array}\right]$
$\displaystyle\begin{array}[]{cc}\hskip 16.0pt\underbrace{\hskip
78.0pt}_{\text{add}}&\hskip 24.0pt\underbrace{\hskip
78.0pt}_{\text{add}}\end{array}$ $\displaystyle\hskip 50.0pt\underbrace{\hskip
130.0pt}_{\text{vectorized chop}}$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{cc}\hskip 46.0pt{\bm{z}}^{(2)}_{1}&\hskip
41.0pt{\bm{z}}^{(2)}_{2}\end{array}\right]$ $\displaystyle\hskip
54.0pt\underbrace{\hskip 130.0pt}_{\text{add \& final chop}}$

With this approach, a MATLAB implementation of a low-precision inner product
${\bm{x}}^{T}{\bm{y}}$, with ${\bm{x}},{\bm{y}}\in\mathds{R}^{n}$ and
$n=2^{p}$, could look like:

  * z = chop(x$\,\mathchoice{\mathbin{{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}*\;$y)’;
for k = 1:p
K = kron(speye(n/2),[1;1]);
z = chop(z*K);
n = n/2;
end

Some comments about the above implementation:

  * •

The statement K = kron(speye(n/2),[1;1]); constructs a sparse representation
of the matrix

$K=I\otimes\left[\begin{array}[]{c}1\\\
1\end{array}\right]=\left[\begin{array}[]{ccc}1&&\\\ 1&&\\\ &1&\\\ &1&\\\
&&\ddots\end{array}\right]$

  * •

The statement z*K computes pair-wise sums of entries in (the row vector) z,
and resizes to half the length. For example, using notation in the previous
diagram, ${\bm{z}}^{(1)T}={\bm{z}}^{(0)T}K$.

  * •

When completed, z has one entry, the inner product ${\bm{x}}^{T}{\bm{y}}$.

  * •

This divide-and-conquer approach requires $\log_{2}(n)+1$ calls to chop, which
is significantly less than the $2n$ or even $n+1$ calls required by the
previously discribed direct summation approaches.

  * •

The code can easily be modified for the case when $n$ is not a power of 2.

  * •

An advantage of this approach is that it uses a natural blocking that reduces
the buildup of roundoff errors when $n$ is large [10].

This pairwise summation procedure can be used for matrix-vector
multiplication. Specifically, if computing $A{\bm{v}}$, we can first scale the
columns of $A$ with entries of ${\bm{v}}$ and call chop once using
vectorization (e.g., C =
chop(A$\,\mathchoice{\mathbin{{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}*\;$v’)),
and then use pairwise summation to add the columns of C, which again requires
only $\log(n)$ additional calls to chop. Similarly, matrix-matrix
multiplication can be performed using pairwise summation of outer products,
and is the approach we use for our low precision preconditioner solve
computation (18).

##  5 Numerical Results

In this section we provide some numerical experiments to illustrate the
effectiveness of using a half-precision Kronecker product SVD approximation as
a preconditioner for Tikhonov regularized least squares problems. We use
Higham’s chop package [11] to simulate half-precision for the preconditioner
solves, and the IR Tools package [7] to generate image deblurring test
problems. Specifically, we consider Gaussian (often used for image deblurring
test problems), speckle (used to simulate blurred observations in telescope
imaging of space objects), defocus (i.e., out-of-focus lens imaging), and
shake (simulating random camera movements) blurs.

The IR Tools package allows to easily generate $A$ (stored in a structured
format that allows for fast matrix-vector multiplications),
${\bm{x}}_{{}_{\mathrm{true}}}$ and
${\bm{b}}_{{}_{\mathrm{true}}}=A{\bm{x}}_{{}_{\mathrm{true}}}$. For each test
problem we add noise ${\bm{\varepsilon}}$ to ${\bm{b}}_{{}_{\mathrm{true}}}$
to obtain the simulated measured data,

${\bm{b}}=A{\bm{x}}_{{}_{\mathrm{true}}}+{\bm{\varepsilon}}\,,$

and for the various tests we report the amount of noise as

$\mbox{noise
level}=\frac{\|{\bm{\varepsilon}}\|_{2}}{\|{\bm{b}}_{{}_{\mathrm{true}}}\|_{2}}\,.$

For each test problem, we use $128\times 128$ images, which means
${\bm{x}}_{{}_{\mathrm{true}}},{\bm{b}}\in\mathds{R}^{16384}$ and
$A\in\mathds{R}^{16384\times 16384}$.

For each blur, it is possible to efficiently decompose $A$ into a sum of
Kronecker products [12, 13],

$A=\sum_{i=1}^{k}A_{r}^{(i)}\otimes A_{c}^{(i)}\,$ (24)

where

$A_{r}^{(1)}\otimes A_{c}^{(1)}=\arg\min_{A_{r},A_{c}}\|A-A_{r}\otimes
A_{c}\|_{F}\,.$

We can therefore use the Frobenius relative error,
$\|A-\widehat{A}\|_{F}/\|A\|_{F}$ as a rough measure of the quality of the
preconditioner’s approximation of $A$, where $\widehat{A}$ is given in
equation (17). These norms, for each of the test problems, is given in Table
1. For the Gaussian blur, $A$ is exactly decomposed into a single Kronecker
product, so the relative error between $A$ and $\widehat{A}$ is due only to
the half-precision approximation.

Table 1: Preconditioner approximations for various test problems. | # of terms in Kron sum |
---|---|---
Blur type | (i.e., $k$ in equation (24)) | $\|A-\widehat{A}\|_{F}/\|A\|_{F}$
Gauss |  1 | $7.9195\times 10^{-8}$
Defocus |  6 | $3.1006\times 10^{-1}$
Speckle | 18 | $4.3102\times 10^{-1}$
Shake | 10 | $6.2834\times 10^{-1}$
Motion |  9 | $7.9212\times 10^{-1}$

###  5.1 Comparison between Different Parameter Choice Methods

We first consider the issue of determining regularization parameters. The
techniques discussed in Section 2.2 require the SVD of $A$, which is very
expensive to compute for large-scale problems. However, because we are using a
half-precision Kronecker product SVD approximation for the preconditioner, we
will use that to compute estimates of the regularization parameter. While the
parameter given by GCV works well for the Gaussian blur type, $\lambda$ tended
to be too small for other blur types, such as defocus and speckle, and the
relative errors blew up at the second iteration. Figure 1(a) was a typical
example of the instant degeneration.

|
---|---
(a) | (b)
Figure 1: A defocus problem of size 128 and 1% noise (a) with GCV; (b) with
weight = $3$ for wGCV

Since GCV chose $\lambda$ too small and under-smoothed the solutions, we
instead used wGCV and set the weight $\omega>1$. For the same problem as in
Figure 1(a), after applying wGCV with $\omega=3$, we obtained a solution that
was comparable to the one using $\lambda_{opt}$. The vertical axis was
rescaled in Figure 1(b) so that the difference in convergence rate is
presented more clearly. The algorithm was not sensitive to the weight: good
results could still be obtained if we set $\omega=5$ or $8$, and the
differences in relative errors were negligible, as demonstrated in Figure 2.

|  |
---|---|---
(a) | (b) | (c)
Figure 2: Result for a defocus problem of size 128 and 1% noise with wGCV (a)
with weight = $3$; (b) with weight = $5$; (c) with weight = $8$

As for the $\|{\bm{\varepsilon}}\|_{2}$-based discrepancy principle, this did
not work well with our approximate SVD preconditioner unless we used a large
value for $\eta$. Figure 3(b) demonstrated that the discrepancy principle with
$\eta=2$ resulted in a $\lambda$ that was comparable to the performance of
$\lambda_{opt}$ for the same problem in Figures 1 and 3(a). Here the vertical
axis was also rescaled for better presentation. However, the discrepancy
principle was more sensitive to $\eta$ than the weight parameter in wGCV;
there were obvious visible differences in the relative errors and the
resulting reconstructed images if $\eta$ was not chosen optimally. Figure 4
showed the results for the same defocus problem with different values of
$\eta$, and the quality of the solution deteriorated quickly as $\omega$
deviates from its optimal value.

|
---|---
(a) | (b)
Figure 3: A defocus problem of size 128 and 1% noise (a) with the discrepancy principle; (b) with $\omega=2$ for the modified discrepancy principle |  |
---|---|---
(a) | (b) | (c)
Figure 4: Result for a defocus problem of size 128 and 1% noise with the
discrepancy principle (a) with $\eta=2$; (b) with $\eta=3$; (c) with $\eta=5$

Moreover, the optimal $\eta$ varied a lot across different problems for the
discrepancy principle. Whereas a suitable weight $\omega$ usually ranged from
$2$ to $8$ for wGCV, the optimal $\eta$ could be as large as $15$ for certain
problem sizes, blur types, and noise levels in the case of the discrepancy
principle. A more systematic study of choosing appropriate weight parameters
is desired, both for wGCV and the discrepancy principle.

###  5.2 Accuracy and Work Unit Evaluation

It is certain that reducing the preconditioning to low precision can achieve
faster computational speed, but the potential loss in accuracy is also to be
considered. Based on the numerical experiments, the difference in relative
errors is negligible. For the same test problem as in Section 5.1, we compared
the relative errors from algorithms that used double-precision preconditioning
and half-precision preconditioning, respectively, as shown in Figure 5. While
the double precision preconditioning achieved smaller relative errors before
convergence, both precision levels converged to the same relative errors in
approimately the same number of iterations; it is difficult to tell the two
lines apart in the plot. Therefore, although a small increase in the relative
errors is incurred when switching to lower precision in the convergence
process, the precision level of the preconditioning did not affect the
relative errors that the algorithm could reach.

|
---|---
(a) | (b)
Figure 5: Comparison between double and half precision preconditioning for a
defocus problem of size 128 and 1% noise (a) with the optimal method; (b) with
wGCV

As previously mentioned, using half precision for preconditioner solves means
that, due to roundoff errors, the preconditioner does not remain constant at
each iteration, and could lead to a loss of orthogonality of the generated
Krylov search directions in PCG [2]. We tested a flexible PCG method [18],
which requires an additional inner product and additional storage [2], but did
not see any difference in convergence behavior. This is probably due to the
fact that we do not need many iterations, so loss of orthogonality is not
significant.

To get a rough comparison of the additional cost of preconditioning, and to
determine if the extra cost is worth preconditioning, we note that for each
PCG iteration, the most expensive computations are matrix-vector
multiplications with $A$ and $A^{T}$, and the preconditioner solve with $M$.
Since we are simulating half-precision in software, a direct timing comparison
between the preconditioned and unpreconditioned algorithms would not be
meaningful. However, we can get a rough estimate if we assume that each
preconditioner solve (in half-precision) is at most $\frac{1}{4}$ the cost of
a double-precision matrix-vector multiplication with $A$ or $A^{T}$. This will
be true if we exploit Kronecker product structure, or if we can use fast
Fourier transforms (FFTs), which is possible for certain spatially invariant
image deblurring problems.

Let $t$ be the time required for a double-precision matrix-vector
multiplication, and let $m_{N}$ and $m_{P}$ denote the total number of
iterations needed, respectively, for the unpreconditioned and preconditioned
algorithms. Then in order for the extra cost of preconditioning to result in a
more efficient algorithm, we need

$\left(2t+\frac{1}{4}t\right)m_{P}<2tm_{N}\quad\Rightarrow\quad
m_{P}<\frac{8}{9}m_{N}\,.$

In the next section, we see that in most cases we achieve this requirement.

###  5.3 Additional Experiments

So far, we have only tested the algorithm on the defocus problem with 1%
noise. In order to fully understand the advantages and pitfalls of our
preconditioning approach, we further applied it to different blur types and
noise levels. In Section 5.1, we investigated several parameter choice methods
and concluded that both wGCV and the discrepancy principle could be effective
with suitably chosen $\omega$ (for wGCV) and $\eta$ (for the discrepancy
principle). Therefore, we directly used $\lambda_{opt}$ in this section to
assess the maximum potential performance of our algorithm. However, we
emphasize that this $\lambda_{opt}$ is computed using the Kronecker product
approximate SVD, and not the SVD of the matrix $A$.

First consider the Gaussian blur, where $A$ can be written as a Kronecker
product, so the only approximation of the preconditioner comes from using
half-precision. As demonstrated by Figure 6, in this ideal situation, our
algorithm converges at the second iteration, a significant improvement on the
no-preconditioning case. In Figure 6(a), there is a jump of the relative error
in the second iteration. Because we are using the same regularization
parameter for both the preconditioner solve and the inverse problem, the
parameter is too small for the preconditioner solve in the first few
iterations. We could use a larger parameter for the preconditioner solve and
eliminate the jump, but since the algorithm quickly recovers from it, we
decided to use the same parameter for both scenarios.

|
---|---
(a) | (b)
Figure 6: A Gaussian problem of size 128 (a) with 1% noise; (b) with 10% noise

As its name implies, the shake blur applies random, shaking motion to each
pixel [7]. As summarized in Table 1, in this case $A$ is a sum of $10$
Kronecker product terms. For this type of blur, our preconditioning algorithm
converges much faster than the no-preconditioning scenario when the noise
level is low ($1\%$ in the test case of Figure 7(a)), but for a high noise
level ($10\%$ in the test case of Figure 7(b)), there is no obvious advantage
in preconditioning. One reason behind the lack of improvement from the
preconditioning is that we can get little information about the true solution
because of the high noise, and therefore the preconditioner is not so useful.

|
---|---
(a) | (b)
Figure 7: A shake problem of size 128 (a) with 1% noise; (b) with 10% noise

The motion blur creates a linear motion with a $45$-degree angle [7]. As
summarized in Table 1, in this case $A$ is a sum of $9$ Kronecker product
terms. However, the performance of the proconditioner on this blur type is
very different from the shake blur: for both noise levels, the preconditioning
provides us with no advantage. We are uncertain of the reason behind such
disparity, it is possible that different regularization parameters could have
an affect on the convergence behavior.

|
---|---
(a) | (b)
Figure 8: A motion problem of size 128 (a) with 1% noise; (b) with 10% noise

The last blur we investigated is the speckle blur, which is used to model blur
caused atmospheric turbulence [7]. As summarized in Table 1, in this case $A$
is a sum of $18$ Kronecker products. Similar to the shake blur, the
preconditioning greatly improves the rate of convergence for low noise level,
but the method is not very helpful when the noise level is high, as shown in
Figure 9. Despite the large number of Kronecker product terms in $A$, our
half-precision, 1-term approximation works effectively.

|
---|---
(a) | (b)
Figure 9: A speckle problem of size 128 (a) with 1% noise; (b) with 10% noise

##  6 Concluding Remarks

In this project, we investigated the use of a low-precision one-term Kronecker
product approximation as a preconditioner for the PCG algorithm to solve ill-
conditioned image deblurring problems. For the Gaussian blur, our
preconditioning method converges significantly faster than the no
preconditioning case since the Kronecker product approximation is exact, and
thus the only approximation comes from using half-precision arithmetic. For
other types of blur, our algorithm is more efficient when the noise level is
low; with high noise, PCG has the same convergence rate as the no
preconditioning method. Thus, our algorithm is more suitable for image
blurring problems with low noise level. Additionally, in our experiments, we
found that flexible PCG has similar performance as PCG. In any case, our low-
precision algorithm is much cheaper than its double-precision counterpart
without loss of accuracy, as discussed in Section 5.2, demonstrating the
potentials of low-precision algorithms.

We also compared different regularization parameter choice methods, using our
half-precision Kronecker product SVD approximation. While the normal GCV and
discrepancy principle tend to choose $\lambda$ too small and blow up the
relative error, we were able to get a reasonable parameter given a suitable
weights ($\omega$ for wGCV and $\eta$ for the discrepancy principle). With
suitably chosen weights, both methods can produce good regularization
parameters, but we found that wGCV was less sensitive to the weight parameter.
We manually chose an appropriate weight through trial and error in this
project, and further research can be done on systematic ways of choosing a
suitable regularization parameter.

## References

  * [1] B. M. Afkham, J. Chung, and M. Chung.  Learning regularization parameters of inverse problems via deep neural networks.  Inverse Problems, 37(10):105017, 2021.
  * [2] H. Anzt, J. Dongarra, G. Flegar, N. J. Higham, and E. S. Quintana-Ortí.  Adaptive precision in block-jacobi preconditioning for iterative sparse linear system solvers.  Concurrency and Computation: Practice and Experience, 31(6):e4460, 2019.
  * [3] Å. Björck.  Least squares methods.  Handbook of numerical analysis, 1:465–652, 1990.
  * [4] D. Calvetti and E. Somersalo.  Introduction to Bayesian Scientific Computing.  Springer-Verlag, New York, 2007.
  * [5] J. Chung, J. G. Nagy, and D. P. O’Leary.  A weighted-gcv method for lanczos-hybrid regularization.  Electronic Transactions on Numerical Analysis, 28:149–167, 2008\.
  * [6] H. W. Engl, M. Hanke, and A. Neubauer.  Regularization of Inverse Problems.  Kluwer Academic Publishers, Dordrecht, 2000.
  * [7] S. Gazzola, P. C. Hansen, and J. G. Nagy.  IR tools: a matlab package of iterative regularization methods and large-scale test problems.  Numerical Algorithms, 81:773 – 811, 2019.
  * [8] G. H. Golub, P. C. Hansen, and D. P. O’Leary.  Tikhonov regularization and total least squares.  SIAM journal on matrix analysis and applications, 21(1):185–194, 1999.
  * [9] P. C. Hansen.  Discrete Inverse Problems: Insight and Algorithms.  SIAM, 2010.
  * [10] N. J. Higham.  The accuracy of floating point summation.  SIAM Journal on Scientific Computing, 14(4):783–799, 1993.
  * [11] N. J. Higham and S. Pranesh.  Simulating low precision floating-point arithmetic.  SIAM Journal on Scientific Computing, 41(5):C585–C602, 2019.
  * [12] J. Kamm and J. G. Nagy.  Kronecker product and SVD approximations in image restoration.  Linear Algebra Appl., 284:177–192, 1998.
  * [13] J. Kamm and J. G. Nagy.  Optimal Kronecker product approximation of block Toeplitz matrices.  SIAM J. Matrix Anal. Appl., 22:155–172, 2000.
  * [14] P. Luszczek, J. Kurzak, I. Yamazaki, and J. Dongarra.  Towards numerical benchmark for half-precision floating point arithmetic.  In 2017 IEEE High Performance Extreme Computing Conference (HPEC), pages 1–5. IEEE, 2017.
  * [15] J. L. Mueller and S. Siltanen.  Linear and Nonlinear Inverse Problems with Practical Applications.  SIAM, Philadelphia, PA, 2012.
  * [16] J. G. Nagy and M. E. Kilmer.  Kronecker product approximation in three-dimensional imaging applications.  IEEE Trans. Image Proc., pages 604–613, 2006.
  * [17] J. G. Nagy, M. K. Ng, and L. Perrone.  Kronecker product approximation for image restoration with reflexive boundary conditions.  SIAM J. Matrix Anal. Appl., 25:829–841, 2004.
  * [18] Y. Notay.  Flexible conjugate gradients.  SIAM J. Sci. Comput., 22(4):1444–1460, 2000.
  * [19] C. F. Van Loan.  The ubiquitous Kronecker product.  J. Comp. Appl. Math., 123:85–100, 2000.
  * [20] C. F. Van Loan and N. P. Pitsianis.  Approximation with Kronecker products.  In M. S. Moonen and G. H. Golub, editors, Linear Algebra for Large Scale and Real Time Applications, pages 293–314. Kluwer Publications, 1993\.
  * [21] C. R. Vogel.  Computational Methods for Inverse Problems.  SIAM, Philadelphia, PA, 2002.