# Commonsense Knowledge Aware Concept Selection For
Diverse and Informative Visual Storytelling

Hong Chen1,3, Yifei Huang1, Hiroya Takamura2,3, Hideki Nakayama 1
With help from the AAAI Publications Committee.

###### Abstract

Visual storytelling is a task of generating relevant and interesting stories
for given image sequences. In this work we aim at increasing the diversity of
the generated stories while preserving the informative content from the
images. We propose to foster the diversity and informativeness of a generated
story by using a concept selection module that suggests a set of concept
candidates. Then, we utilize a large scale pre-trained model to convert
concepts and images into full stories. To enrich the candidate concepts, a
commonsense knowledge graph is created for each image sequence from which the
concept candidates are proposed. To obtain appropriate concepts from the
graph, we propose two novel modules that consider the correlation among
candidate concepts and the image-concept correlation. Extensive automatic and
human evaluation results demonstrate that our model can produce reasonable
concepts. This enables our model to outperform the previous models by a large
margin on the diversity and informativeness of the story, while retaining the
relevance of the story to the image sequence.

## Introduction

Telling a story based on a sequence of images is a natural task for humans and
a fundamental problem for machine intelligence for various scenarios such as
assisting the visually impaired people. Also known as visual storytelling
(VST), the task has raised extensive research attention, since VST requires
the model to not only understand the complex content within one image but also
reason about the event across images as they occur and change. Since image
sequences contain rich and diverse information, it is especially difficult for
a model to tell a relevant story that is both informative of the image content
and diverse in story style.

Most previous works on VST construct end-to-end frameworks (Yang et al. 2019;
Wang et al. 2018; Jung et al. 2020; Yu, Bansal, and Berg 2017). However,
although these methods can produce legitimate stories with high score in
automatic metrics like BLEU (Papineni et al. 2002), it is shown that the
stories tend to be monotonous which contains limited lexical diversity and
knowledge (Hsu et al. 2019a) (see the example in Figure 1). Recently, two-
stage generation methods, also known as plan-write strategy, aroused much
research attention in story generation tasks (Yao et al. 2019; Martin et al.
2017; Ammanabrolu et al. 2020). When adopted to the task of VST, Hsu et al.
(2019a) shows that this strategy is capable of generating more diverse stories
compared with end-to-end methods. However, their method directly generates
concepts from the images using sequence-to-sequence models. Since the concept
is selected from the full vocabulary, this kind of direct generation often
produces concepts of low quality which affects the informativeness of the
story.

Figure 1: Stories generated by the state-of-the-art model (Jung et al. 2020)
and our proposed model using concept selection (red). The state-of-the-art
model tends to generate similar stories (blue) for different inputs. Compared
with it, our model can generate more informative and diverse stories. Figure
2: An overview of our visual storytelling model. The image features are
obtained by a pretrained CNN combined with a bi-LSTM layer. The concepts are
obtained from a concept detection model and enriched by ConceptNet (Liu and
Singh 2004). These concepts from the nodes in a graph and are connected
according to the relationship in the knowledge base. Initialized by the word
embedding vector, the concept features are then updated by a Graph Attention
Network. Our proposed concept selection module is then applied to select exact
concept words using the image features and the concept features. Finally, both
image features and concept words are used to generate a full story.

In this work we aim to generate stories that are both diverse and informative
for a given input image sequence. Taking the advantage of the previous two-
stage models, we detect image concepts and construct concept graphs for
proposing a set of concept candidates, and propose two novel methods for
better selecting the appropriate concept for the second generation stage.
After detecting the concept in each input image, we first extend the concepts
into a larger commonsense graph using ConceptNet (Liu and Singh 2004). This
extension step increases the informativeness of generated stories. Since
selecting appropriate candidates from the concept graph is critical for
generating stories of good quality, a natural way is to use a graph attention
network (Veličković et al. 2017) to refine the node features. Using Graph
Attention Network can allow message passing along the graph structure, so that
information of related concepts can be updated and integrated. This would
allow us to get a better pool of concept candidates.

For selecting the most adequate concept from the candidates as the input to
the second stage of our model, two novel modules are proposed in this work.
The first one, named Sequential Selection Module (SSM), operates in a
straightforward manner that uses an encoder-decoder for selecting concepts for
each image. Differently from SSM, the second module called Maximal Clique
Selection Module (MCSM) processes the concept graph as a whole. It learns a
probability for each concept in the training phase, and during inference it
finds a maximal clique using the Bron Kerbosch algorithm (Bron and Kerbosch
1973). The concepts within the clique are used for the next story generation
step. Our experiments show that improved quality of concept selection can
greatly help to increase the diversity of the generated stories while keeping
the relevance with the input images.

The second stage of our model generates a story with the image features and
the selected concepts. Other than using the same module for fair comparison
with existing works, we also propose to modify the large scale pre-trained
model BART (Lewis et al. 2019) to input the images and concepts and output the
full stories.

We conduct extensive experiments on the public VIST dataset (Huang et al.
2016). Our experiments demonstrate that using our proposed concept selection
modules, our generated stories can achieve better performance on both
automatic metric and multiple human evaluation metrics using the same
generation module. When equipped with BART, the quality of the stories can be
remarkably improved, with the generated story diversity similar to human
writing.

In summary, our main contributions are listed as follows:

  * •

We propose two novel modules SSM and MCSM to select concepts from the given
candidates concepts under a plan-write two-stage visual storytelling system.
The experiments show that our proposed methods can output more appropriate
concepts than the previous work.

  * •

We exploit modified BART as our story generation module to mitigate the
problem caused by limited vocabulary and knowledge in the dataset. To the best
of our knowledge, this is the first work to use a large scale pre-trained
language model in a visual storytelling task.

  * •

Large scale experiments using automatic metrics and human evaluation show that
our model can outperform previous models by a large margin in both diversity
and informativeness, while retaining the relevance and logicality as the
previous work.

Figure 3:  (a) Sequential Selection Module: this module sequentially use
current hidden state as a query to select the concept from the commonsense
graph. (b) Maximal Clique Selection Module: this module calculates concept-to-
concept and image-to-concept correlation maps. The correlation maps can be
viewed as a fully connected graph of concepts. We set a threshold to prune the
edges. Then, a maximal clique algorithm is applied to find the maximal cliques
in the remaining graph. Finally, those cliques are scored and the one with
highest score will be selected and the concepts to be used in the next
generation stage.

## Related Work

Visual storytelling aims at generating stories for image sequences. Many
existing methods focused on generating relevant and coherent stories (Hu et
al. 2020; Wang et al. 2020; Chandu et al. 2019; Li et al. 2019; Hsu et al.
2019b). These works can be separated into two lines: one line is to construct
end-to-end models to generate the stories directly from the images. The other
line is to build a two-stage generation model that first outputs a mid-level
abstraction and then generates the full story.

### End-to-End Methods

Wang et al. (2018) proposed a visual storytelling framework which is widely
used as a base model in the coming-up studies. This framework uses an end-to-
end structure that first convert the image into features and then transfer its
information to the adjacent images by a BiLSTM layer. Finally, a decoder
decodes the features separately and merge the sentences into a story. While
many succeeding works (Huang et al. 2019; Jung et al. 2020) can achieve high
automatic scores, the story may not be interesting and informative (Hu et al.
2020) for human as they often contain repetitive texts and limited
information. One of the main reasons for the low diversity and informativeness
is that these model are trained end-to-end under the maximum likelihood
estimation (Yang et al. 2019).

### Two-stage Methods

To alleviate the low diversity problem, Hsu et al. (2019a) proposed to
generate several concepts before outputting the full stories. The discrete
concept words can guide the decoder to produce more diverse stories. This
plan-and-write strategy (Yao et al. 2019) can substantially increase the
diversity of the stories. During the planning, to enhance the concepts that
models can obtain, some researchers (Hsu et al. 2019a; Yang et al. 2019)
introduce external commonsense knowledge database such as OpenIE (Angeli,
Premkumar, and Manning 2015), Visual Genome (Krishna et al. 2017) or
ConceptNet in the VST task. Their results show that using external knowledge
base helps to generate more informative sentences. These works also show that
the concepts are critical for the quality of generated stories because they
can control the story flow.

In this work we aim to improve the concept selection for increasing the
diversity and informativeness of VST. We propose two concept selection modules
that carefully selects concept from a well-designed pool of concept
candidates. The stories generated using our selected concepts thus become more
diverse and informative. We further introduce to modify the pretrained model
BART and use it to generate even better stories.

## Method

Figure 2 depicts an overview of our proposed model. Given a sequence of $N$
images as input, our model 1) encode image features, 2) construct a large
commonsense graph for the images, 3) update concept feature in the graph, 4)
select the concepts from the graph and 5), send concepts and image features
into the decoder to output the full story. The details of each step are as
follows:

### Image Feature Encoding

We send the images into ResNet152 (He et al. 2016) to obtain image features
$I=\\{I_{1},...,I_{N}\\}$. Following Wang et al. (2018), a bidirectional GRU
further encodes the high-level visual features to capture sequential
information. The encoded feature for each image contains both the information
of itself and the information from the adjacent images. Note that position
embedding is applied before sending the image features into GRU to identify
the order of the images.

### Commonsense Graph Construction

To build our commonsense knowledge graph for image sequences, we need some
seed concepts. Following Yang et al. (2019), we use clarifai (Sood 2015) to
obtain the top 10 seed concepts from each image. Each concept is used as a
query to select relative commonsense concepts in the ConceptNet (Liu and Singh
2004). Since the number of the commonsense concepts is usually very large
$(>500)$, we make several rules to filter some concepts which are less useful:

  * •

Remove the commonsense concepts that appear less than 5 times in the whole
training corpus.

  * •

Remove the commonsense concepts that do not co-occur with the seed concepts
either in the sentence or in the training corpus.

  * •

If the concept number is still larger than $K$ for one image. We simply
randomly sample $K$ words from it.

After the filtering process, each image contains $K$ concepts. Note that while
different images can obtain the same concept in one image sequence, they can
represent different semantic meanings in different positions of the story.
Each concept forms a node in the commonsense graph. An undirected edge is
established between concepts if they are related in ConceptNet. Also, a
concept in one image will connect to the related concepts in the adjacent
images to allow information flow between images. Like (Yang et al. 2019), we
do not use the specific relation (e.g., isA, has) between concepts. Till now,
we build a graph which is a graph structure that is both connected within an
image and between images.

### Concept Features Update

The concept features are initialized with word embedding vectors. To
incorporate the visual information into the concepts, we also connect the
image feature to its corresponding concept features in the graph. These
features are updated by a two-layer Graph Attention Network, which passes
information between connected concepts and image using attention mechanism.

### Concept Selection Module

We propose two methods to select concepts given the concept features and the
image features.

To better formalize the procedure in the methods, we denote $c^{i,j}$ as the
j-th concept of the i-th ($1\leq i\leq N$) image. we let
$\mathcal{C}_{S}=\\{c_{S}^{1,1},...,c_{S}^{N,K}\\}$ and
$\mathcal{C}_{G}=\\{c_{G}^{1,1},...\\}$ denote the concepts set in the source
candidate concepts and the full word set in the gold story, respectively. The
target concepts are their intersection:
$\mathcal{C}_{T}=\mathcal{C}_{S}\cap\mathcal{C}_{G}$.

#### Sequential Selection Module (SSM)

One straightforward way of selecting concepts is to adopt an encoder-decoder
model where we can forward the updated concept features into the encoder, and
the decoder will output the selected concepts. Inspired by the Copy Mechanism
(Gu et al. 2016), instead of generating a probability distribution with
vocabulary size in each step, the SSM outputs are directly chosen from the
inputs $\mathcal{C}_{S}$. As shown in Figure 3(a), we use a GRU (Cho et al.
2014) to first encode the concept embedding feature $v_{S}^{t-1}$ and the
hidden state into a new hidden state $h^{t}$. We then use $h^{t}$ to query all
the concepts in $\mathcal{C}_{S}$ to get a probability $p_{S}$ for each
concept in the source set. Finally the concept with the highest probability is
selected as the output concept, while its feature is directly copied for the
generation of the next step:

$\displaystyle\begin{split}h^{t}&=GRU\left(h^{t-1},v^{t-1}_{S}\right)\\\
p_{S}&=softmax\left((W_{h}h^{t})^{\mathrm{T}}W_{c}V_{S}\right)\\\
c_{S}^{t}&=argmax(p_{S})\end{split}$ (1)

Here $W_{h}$ and $W_{c}$ are trainable projection matrices. The objective
function is to maximize the probability score of the concepts which locate in
$C_{T}$.

$\displaystyle\mathcal{L}_{ssm}$ $\displaystyle=-\Sigma y_{S,T}\log(p_{S}),$
(2)

where $y_{S,T}$ is an indicator of whether a concept in $\mathcal{C}_{S}$ is
in $\mathcal{C}_{T}$. The sequence selection step stops when the module
generates $<$end$>$ token. This $<$end$>$ token is added to the set of
candidate concepts with a uniform random initialized feature without any
update during the training phase. The same procedure is done to the
$<$start$>$ token except that it is not involved in the candidates.

#### Maximal Clique Selection Module (MCSM)

Different from SSM, this method aims to calculate the co-occurrence
probability of all candidate concepts $c_{s}$ in the graph. An illustration of
MCSM is shown in Figure 3(b). In the beginning, we calculate self-attention to
compute a correlation matrix $M_{C}\in(NK\times NK)$ which contains the
correlation between each pair of nodes. We also calculate another correlation
matrix for each image $M_{I}\in(N\times K)$ indicating the correlation between
the concept embedding feature ($v_{S}$) and image features ($I$).

$\displaystyle\begin{split}M_{C}&=\sigma(v_{S}^{\mathrm{T}}W_{a}^{\mathrm{T}}W_{b}v_{S})\\\
M_{I}&=\sigma(I^{\mathrm{T}}W_{c}^{\mathrm{T}}W_{d}v_{S})\end{split}$ (3)

Here, $W_{a},W_{b},W_{c},W_{d}$ is trainable weights, $\sigma$ denotes sigmoid
activation function. Intuitively, the concepts that appear in a gold story
should own high correlations with each other, and the image should be highly
related to the gold concepts to describe it. Thus, our target correlation maps
can be written as follow:

$\displaystyle\begin{split}\hat{M}_{C}[i,j]&=\begin{cases}1,&c_{i}\in\mathcal{C}_{T}\wedge
c_{j}\in\mathcal{C}_{T}\\\ 0,&\text{otherwise}\\\ \end{cases}\\\
\hat{M}_{I}[i,j]&=\begin{cases}1,&c_{j}\in\mathcal{C}^{i}_{T}\\\
0,&\text{otherwise}\\\ \end{cases}\end{split}$ (4)

Then, the objective is to minimize the difference between predicted and target
correlation maps:

$\displaystyle\mathcal{L}_{mcsm}=\lambda_{1}\left\|M_{C}-\hat{M}_{C}\right\|_{2}^{2}+\lambda_{2}\left\|M_{I}-\hat{M}_{I}\right\|_{2}^{2}$
(5)

In testing phase, $M_{C}$ can be viewed as a fully connected graph in which
the edge weights correspond to the values in the matrix. Therefore, a low edge
weight means less co-occurrence probability between two concepts. Based on
this assumption, we set a threshold $\tau$ to remove the edges whose weight is
less than $\tau$. Then we apply Bron Kerbosch algorithm (Bron and Kerbosch
1973) to find all maximal cliques from the remaining sub-graph. Finally, we
score each of them with equation 6 and select a clique with maximum score $s$.
The output concepts are the nodes within the selected cliques.

$\displaystyle\begin{split}s&=s_{C}+s_{I}\\\
s_{C}&=\frac{1}{(\left\|\mathcal{C}_{P}\right\|-1)\left\|\mathcal{C}_{P}\right\|}\sum_{i}\sum_{j\neq
i}log(M_{C}[i,j])\\\
s_{I}&=\frac{1}{\left\|\mathcal{C}_{P}\right\|}\sum_{i=1}^{N}\sum_{c_{j}\in\mathcal{C}_{P}^{i}}log(M_{I}[i,j]).\end{split}$
(6)

where $\mathcal{C}_{P}$ denotes the concepts in a clique, and
$\mathcal{C}_{P}^{i}$ denotes the concept of the i-th image in the clique.

### Concept to Story Module

The selected concepts are assigned to its corresponding image to generate the
sentences. We tried two kinds of encoder-decoder to decode the story: 1) a
simple encoder-decoder module that uses multi-head pooling to encode the
concept embeddings, and decode the sentences with a RNN decoder. 2) a large
scale encoder-decoder which both can encode the input and output the
sentences.

RNN Decoder decodes sentences separately for each image and then concatenates
into a story, while BART accepts all the images and concepts at once to output
a full story.

#### RNN

To get the concept feature from the concept words, we apply Multi-head Pooling
(Liu and Lapata 2019) to calculate multi head self-attention score on the
input concept embedding and do weighted summation on the concept embedding.
Each image and corresponding concepts are decoded separately with the same
decoder as Jung et al. (2020). The decoder accepts the image features $I_{i}$
and the pooled concept feature $v_{i}$ as input. Formally, the generation
process can be written as:

$\displaystyle h^{t}_{i}$
$\displaystyle=RNN\left(h^{t-1}_{i},\left[w_{t-1}^{i};I_{i};v_{i}\right]\right)$
(7) $\displaystyle\pi^{t}$ $\displaystyle=softmax\left(W_{s}h^{t}_{i}\right)$
(8)

Here $W_{s}$ is a projection matrix for outputting the probability
distribution $\pi^{t}$ on the full vocabulary. “;” denotes the channel-wise
concatenation operation. Finally, all of the corresponding words are merged
into a complete story.

#### BART

Although there exists many large scale pretrained language models, none of
them are used on VST in previous works. In this paper, we propose to use a
modified version of BART as our base decoder model. BART is pretrained on a
large dataset and their arbitrary noise function helps to mitigate the
exposure bias. For inputting the image features into BART, We add a new
randomly initialized projection matrix to transform the size of image features
to fit BART’s input. For the concepts we simply use the pretrained BART word
embedding. Since BART is powerful enough to generate a long story. Different
with RNN decoder, we input N images and all concepts into the BART together
and use a $<$SEP$>$ token to separate the images and concepts, and to separate
the concepts from different images. The transformed image feature and concept
embedding are sent to BART and generate the final stories.

## Experiment

### Dataset

We conduct experiments on the widely used VIST dataset (Huang et al. 2016),
which consists of 40101 samples for training 5050 test samples. Each sample
contains 5 images and one corresponding gold story. For fair comparison, we
follow the same experiment setting as (Jung et al. 2020) except that we set
the vocabulary size to 28000. All models use the same fixed random seed.

### Concept Detector

In this paper, we use the Concept Detection API produced by clarifai
111www.clarifai.com. This Detector can predict 11,000 concepts which is larger
than any other detection model. This powerful pretrained detector helps us to
precisely find out the concepts inside the images, so that the word in the
gold sentences can be easily involved in our knowledge enhanced candidate
concepts.

### Implementation Details

When training SSM, since we assume that there is no order relationship between
the concepts in one image, so during the training phase, we randomly shuffle
the target concept in one image. When training BART, we conduct a two-stage
training: 1) freeze all BART parameters except for the image projection
matrix. 2) finetune all parameters. We use Adam (Kingma and Ba 2014) optimizer
with an initial learning rate of 4e-4 in the first stage, then the learning
rate is decreased to 4e-5 in the fine tuning stage. Each stage is trained for
5 epochs. All the other parts of our model are trained with an Adam Optimizer
with learning rate 4e-4. During training, we follow Jung et al. (2020) to
blind one of the images starting from the 50-th epoch and increase the
blinding into two images from epoch 80. The training stops at epoch 100. Our
model uses gold concepts extracted from gold stories to train the concept to
story model. This step is similar to the common auto-regressive models that
use the target token as the input to generate the next token. As has been
discussed in the previous sections, this kind of generation often meets with
the problem caused by the train-test discrepancy that we cannot see the gold
concepts in the testing phase. To alleviate, a simple and effective way is to
add noise to the inputs. In this work, we add two kinds of noise into the
inputs in the story generation module: masking and random replacement. We mask
30% concepts and replace 20% of them into other similar words in training.

###  $\tau$ search in Maximal Clique Selection Module

$\tau$ is set as the threshold in pruning edges for the maximal clique
selection. Larger $\tau$ leads to fewer concepts that can be selected and will
further result in the lack of imagination (diversity) in the generated story.
However, smaller $\tau$ would lead to too many concepts that may mislead the
model to generate irrelevant stories. To make a trade-off, we initialize
$\tau$ as 0.3 and continual decreasing the number until Bron Kerbosch
algorithm (Bron and Kerbosch 1973) can produce at least 5 candidate cliques
that contains 7 to 15 candidate concepts in each clique.

Method | Precision | Recall | F measure
---|---|---|---
Rand | 2.68 | 2.45 | 2.56
C_Attn | 30.38 | 43.37 | 35.86
I2C | 31.32 | 20.75 | 24.96
SSM | 40.43 | 40.30 | 40.36
MCSM | 45.30 | 40.90 | 42.99
Table 1: Concept selection performance of different methods. The results show
that our MCSM achieved the best f-score among all methods.

### Decoding strategy

During model inference, usually beam search is used in decoding the sentences
from the decoder (Wang et al. 2018; Jung et al. 2020; Yu, Bansal, and Berg
2017). However, it has been proved that beam search can result in bland,
incoherent stories, or gets stuck in repetitive outputs (Holtzman et al.
2019). In this paper, to further improve diversity in generated stories, we
apply nucleus sampling (Holtzman et al. 2019) which outperforms top-k sampling
method. However, nucleus sampling does not perform well without a well-trained
language model. We show this by an experiment where we train our model with
RNN as the story generation module on the VIST dataset and decode using
nucleus sampling. Note the VIST dataset is relatively small compared with
other natural language datasets like Book Corpus (Zhu et al. 2015). We conduct
human evaluation and let the annotators to rank the overall story qualities
(lower the better) generated with different softmax temperatures
[0.5,0.6,0.7,0.8,0.9] for the same image sequence. We randomly pick 100
samples for this experiment. Figure 4 shows that with the temperature
increases, the story quality would drop. Especially, when the temperature is
0.9 which is recommended in Holtzman et al. (2019), the generated sentences
are incomprehensible and full of grammatical errors, as shown in the example
of Figure 4. This result motivates us to use a large scale pre-trained
language model in our experiment.

Figure 4: The y axis shows the ranking score (lower the better) and Distinct-4
score (higher the better) while we change the temperature in nucleus sampling
with the model trained from scratch on VIST dataset. We ask the workers to
rank the overall score for five stories generated by 5 different temperatures.
As we can see, with the increasing of the temperature, the stories become more
diverse, however, the quality of them become lower.

Thus, for fair comparison with previous works, the same RNN with beam search
is used as the story generation module to validate the effectiveness of our
concept planning model. While, we use nucleus sampling (temperature=0.9,
p=0.9) when BART (Lewis et al. 2019) is applied to demonstrate a higher upper
bound of our model’s ability.

### Experiment on Concept Selection

We first test the ability to select appropriate concepts for different models.
Since each image sequence in the test sample corresponds to 3 to 5 gold
stories, we report the highest performance of the output selection result with
respect to all gold stories in each sample. Similar as Keyphrase generation
tasks, we apply precision, recall and f measure to evaluate the efficiency of
concept selection methods. The precision is defined as the number of correct
concepts selected over the target concepts, and recall is defined as the
number of correct selections over all candidates.

We compare among several methods:

  * •

Rand: A simple baseline where we randomly pick 3 concepts from the candidates
for each image. On average, each image contains 2.65 gold concepts.

  * •

C_Attn: We extract the attended concepts where the attention score is larger
than a threshold from the model of Yang et al. (2019). This is an end-to-end
model with sigmoid attention on concept words. We choose 0.8 as the threshold
since this contributes the best f-score.

  * •

Image to concept(I2C): This is a straightforward version of concept selection
where the concepts are directly generated from the images. We simply add a
projection layer on each hidden state to predict the concept words from the
vocabulary size of the concepts, which is very similar to the model of Hsu et
al. (2019a).

  * •

SSM: Our proposal which uses a copy mechanism in each step of selection.

  * •

MCSM: Our proposal which calculates the correlation score for concept-concept
and image-concept and uses maximal clique selection.

Method | Dist-2 | Dist-3 | Dist-4
---|---|---|---
INet$\smwhitestar$ | 8.36 | 18.92 | 31.02
KS$\smwhitestar$ | 10.84 | 22.90 | 36.51
KG-Story† | 18.73 | 38.65 | 57.22
Our(MCSM) | 13.98 | 34.01 | 54.11
Image+BART† | 21.63 | 46.23 | 67.57
Our(MCSM)+BART† | 34.95 | 69.88 | 88.74
Gold | 47.76 | 82.27 | 95.05
Table 2: Diversity of generated stories by different methods. Two-stage
generation methods can produce more diverse stories. Using BART, we can
achieve the diversity close the human writing, while achieving same level
story quality in other aspect. †denotes the story generation module is pre-
trained with other dataset. $\smwhitestar$ denotes end-to-end methods.

Qualitative results are shown in Table 1. We can see that our proposed SSM and
MCSM can achieve significantly higher f-score than other methods. This helps
our model to keep the story relevance to the input images while generating
diverse stories.

### Experiment on Visual Storytelling

Here we show the results of the visual storytelling. We use the following
baselines for comparison:

INet (Jung et al. 2020) This is a recent work which uses a “hide-and-tell”
strategy to train an end-to-end model. In this method no concept is used.

KS (Yang et al. 2019) This method uses sigmoid attention to incorporate
concept features into the model. We change the structure of the visual encoder
and decoder the same as INet for fair comparison.

KG-Story† (Hsu et al. 2019a) is a strong baseline that use two stage plan-
write strategy and pretrain the decoder on RocStories Corpora (Mostafazadeh et
al. 2017). †indicates the model introduces external large-scale pretraining
data.

Image+BART†is an end-to-end baseline that uses BART on top of image features
to directly generate the story. This baseline is one-stage that does not
generate concepts.

We also change the concept selection module and story generation module in our
model to validate the effectiveness of each component. Specifically, we
compare: Rand+RNN, C_Attn+RNN, SSM+RNN, MCSM+RNN, and MCSM+BART† .

#### Comparison on diversity

We first compare the ability of generating diverse stories of different
models. Quantitative comparison is shown in Table 2. We report Distinct-n
(Dist-n) scores (Li et al. 2015) that calculate the percentage of unique
n-gram in all generated stories in the test dataset. Higher score means less
inter-story repetition. From the table, two stage generation methods (KG-Story
and ours) can achieve significantly higher diversity scores. Our MCSM can
generate the most diverse stories among all the methods without using external
pretrained models. When equipped with BART, we can even achieve diversity
close to human writing. We show in the following that the increased diversity
also improves the overall quality of the story.

Method | B-3 | B-4 | R | M | C
---|---|---|---|---|---
INet$\smwhitestar$ | 23.5 | 14.4 | 29.7 | 35.3 | 9.5
KS$\smwhitestar$ | 24.7 | 15.0 | 31.0 | 35.0 | 9.8
Rand+RNN | 13.3 | 6.1 | 27.2 | 31.1 | 2.2
C_Attn+RNN | 20.7 | 11.2 | 29.7 | 34.5 | 7.8
SSM+RNN | 22.1 | 12.0 | 30.0 | 35.4 | 10.5
MCSM+RNN | 23.1 | 13.0 | 30.7 | 36.1 | 11.0
Table 3: Automatic metric in story quality. We report BLEU (B), METEOR (M), ROUGH-L (R), and CIDEr (C) scores. The two-stage generation can achieve higher METEOR and CIDER scores. Choices(%) | MCSM vs INet | MCSM vs KS | MCSM vs SSM | MCSM+BART† vs KS | MCSM+BART† vs MCSM
---|---|---|---|---|---
| MCSM | INet | MCSM | KS | MCSM | SSM | MCSM+BART | KS | MCSM+BART | MCSM
Revelence | 47.4 | 35.6 | 26.3 | 31.6 | 50.5 | 40.0 | 28.8 | 33.6 | 35.2 | 35.2
Informativeness | 51.0* | 31.6 | 46.3* | 28.9 | 44.7 | 41.2 | 62.5** | 18.8 | 58.8** | 23.5
Logicality | 35.5 | 34.3 | 34.2 | 29.0 | 32.9 | 42.3 | 35.3 | 33.3 | 40.2 | 37.5
Overall | 55.0** | 30.0 | 44.7 | 34.2 | 48.3 | 37.1 | 43.5** | 23.0 | 47.0* | 31.6
Table 4: Human evaluation. Numbers indicate the percentage of annotators
believe that a model outperforms its opponent. Methods without (+BART) means
using RNN as the story generation module. Cohen’s Kappa coefficients
($\kappa$) for all evaluations are in Moderate (0.4-0.6) or Fair (0.2-0.4)
agreement, which ensures inter-annotator agreement. We also conduct a sign
test to check the significance of the differences. The scores marked with *
denotes $p<0.05$ and ** indicates $p<0.01$ in sign test.

#### Automatic evaluation

In Table 3, for comparing the quality of generated stories, we use automatic
metrics BLEU (B) (Papineni et al. 2002), METEOR (M) (Banerjee and Lavie 2005),
ROUGH-L (R) (Lin 2004), and CIDEr (C) (Vedantam, Lawrence Zitnick, and Parikh
2015). Note that it remains tricky for automatic scores to appropriately
evaluate story qualities. Using concept, KS can achieve better performance
than INet that does not use concept. From the comparison of the variants of
our model, we can see that better concept selection can lead to better
automatic scores. With reasonable concept selection, our SSM and MCSM can
achieve highest METEOR and CIDER scores.

#### Human Evaluation

To better evaluate the quality of generated stories, we conduct human
evaluation to compare pair-wise outputs with several models via the Amazon
Mechanical Turk (AMT). We randomly sample 200 image sequences from the test
set and generate stories using each model. For each sample pair, two
annotators participate in the judgement and decide their preference on either
story (or tie) in terms of Relevance, Informativeness, Logicality and Overall.
Relevance evaluates how relevant the stories and the images are.
Informativeness assesses how much information can be achieved in the generated
stories, and this score from one side reflects the diversity of stories.
Logicality evaluates the logical coherence in the stories. Overall is a
subjective criterion that shows the preference of workers.

Table 4 shows the human evaluation result. Since there exists randomness in
human evaluation, we compute the Cohen’s Kappa coefficient and found that all
evaluations are in Moderate Agreement and Fair agreement, which indicates the
evaluation result is reasonable as good inner agreement between evaluators is
reached. We also conduct a Sign test to illustrate the significance of the
evaluation difference: if p is below 0.05 it would indicate that the two
compared models have a significant performance difference. From the comparison
between MCSM and INet and the comparison between MCSM and KS, we can see that
our two-stage planning method greatly outperforms the end-to-end models,
especially in the informativeness score. The MCSM module also outperforms the
SSM module, which indicates positive correlation between the quality of
concept selection and the overall quality of generated stories. Finally, using
BART with MCSM can help to achieve further informativeness and generate even
better stories.

Figure 5: We calculate the Pearson Correlation Coefficient for each criteria
in human evaluation. R, I, L, O denotes Relevance, Informativeness, Logicality
and Overall, respectively. We can see that Informativeness is almost
independent to Relevance and Logicality, while is highly correlated to Overall
score.

### Importance of informativeness in story quality

We calculate the Pearson Correlation Coefficient on four criteria in human
evaluation. In Figure 5, R, I, L, O denotes Relevance, Informativeness,
Logicality and Overall, respectively. Ranged from -1 to 1, the Informativeness
score has low correlation score with Logicality (0.011) and Relevance (0.08),
while a high correlation score with Overall (0.47). This indicates that
Informativeness is almost independent on Relevance and Logicality, but highly
dependent on the Overall score. This suggests that humans tend to choose
stories with more interesting information. This phenomenon proves the
significance of informativeness and diversity in visual storytelling.

Figure 6: The examples of generated stories by different methods. Our MCSM and
SSM can generate better stories compared with other baselines that do not use
BART. When using the pretrained BART, the concept selection with MCSM can
produce vivid and informative story.

### Case Study

We show a qualitative result of a random test sample in Figure 6. This is a
hard example because the last three images are very similar and the objects in
all images are hard to recognize. We can see that INet generates monotonous
and even irrelevant sentences. KS can generate better sentences but still low
in lexical diversity. For the stories generated by two-stage strategy with RNN
(SSM+RNN, MCSM+RNN), we can see that the story follows the selected concepts
and the stories seem more reasonable than that of end-to-end training methods.
When using BART, we compare three methods that represent no concept selection
(Image+BART), bad concept selection (Rand+BART) and ours concept selection
(MCSM+BART). We can see that without using concepts or using randomly selected
concepts, the generated stories are of low quality and to a certain extent
irrelevant to the images. However, when guided by the selected concept, the
story becomes vivid, relevant and logical.

## Conclusion

In this work we exploit concept selection for improving the diversity and
informativeness of stories generated from image sequences. By constructing a
commonsense graph and two novel modules for concept selection, our proposed
model outperforms all previous works in diversity by a large margin while
still preserving the relevance and logical consistency on the VIST dataset.
Our future direction aims to increase the relevance of the generated story by
better leveraging the visual information.

## Acknowledgements

We thank the anonymous reviewers for the useful comments. This paper was based
on results obtained from a project, JPNP20006, commissioned by the New Energy
and Industrial Technology Development Organization (NEDO) and was also
supported by JSPS KAKENHI Grant Number JP19H04166.

## References

  * Ammanabrolu et al. (2020) Ammanabrolu, P.; Tien, E.; Cheung, W.; Luo, Z.; Ma, W.; Martin, L. J.; and Riedl, M. O. 2020.  Story Realization: Expanding Plot Events into Sentences.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , 7375–7382.
  * Angeli, Premkumar, and Manning (2015) Angeli, G.; Premkumar, M. J. J.; and Manning, C. D. 2015.  Leveraging linguistic structure for open domain information extraction.  In _Proceedings of the 53rd Annual Meeting of the Association for Computational Linguistics and the 7th International Joint Conference on Natural Language Processing (Volume 1: Long Papers)_ , 344–354.
  * Banerjee and Lavie (2005) Banerjee, S.; and Lavie, A. 2005.  METEOR: An automatic metric for MT evaluation with improved correlation with human judgments.  In _Proceedings of the acl workshop on intrinsic and extrinsic evaluation measures for machine translation and/or summarization_ , 65–72.
  * Bron and Kerbosch (1973) Bron, C.; and Kerbosch, J. 1973.  Algorithm 457: Finding All Cliques of an Undirected Graph.  _Commun. ACM_ 16(9): 575–577.  ISSN 0001-0782.  doi:10.1145/362342.362367.  URL https://doi.org/10.1145/362342.362367.
  * Chandu et al. (2019) Chandu, K.; Prabhumoye, S.; Salakhutdinov, R.; and Black, A. W. 2019.  “My Way of Telling a Story”: Persona based Grounded Story Generation.  In _Proceedings of the Second Workshop on Storytelling_ , 11–21.
  * 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_ .
  * Gu et al. (2016) Gu, J.; Lu, Z.; Li, H.; and Li, V. O. 2016.  Incorporating copying mechanism in sequence-to-sequence learning.  _arXiv preprint arXiv:1603.06393_ .
  * He et al. (2016) He, K.; Zhang, X.; Ren, S.; and Sun, J. 2016.  Deep Residual Learning for Image Recognition.  In _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 770–778.
  * Holtzman et al. (2019) Holtzman, A.; Buys, J.; Du, L.; Forbes, M.; and Choi, Y. 2019.  The curious case of neural text degeneration.  _arXiv preprint arXiv:1904.09751_ .
  * Hsu et al. (2019a) Hsu, C.-C.; Chen, Z.-Y.; Hsu, C.-Y.; Li, C.-C.; Lin, T.-Y.; Huang, T.-H.; and Ku, L.-W. 2019a.  Knowledge-Enriched Visual Storytelling.  _arXiv preprint arXiv:1912.01496_ .
  * Hsu et al. (2019b) Hsu, T.-Y.; Huang, C.-Y.; Hsu, Y.-C.; and Huang, T.-H. 2019b.  Visual story post-editing.  _arXiv preprint arXiv:1906.01764_ .
  * Hu et al. (2020) Hu, J.; Cheng, Y.; Gan, Z.; Liu, J.; Gao, J.; and Neubig, G. 2020.  What Makes A Good Story? Designing Composite Rewards for Visual Storytelling.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , 7969–7976.
  * Huang et al. (2019) Huang, Q.; Gan, Z.; Celikyilmaz, A.; Wu, D.; Wang, J.; and He, X. 2019.  Hierarchically structured reinforcement learning for topically coherent visual story generation.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 33, 8465–8472.
  * Huang et al. (2016) Huang, T.-H.; Ferraro, F.; Mostafazadeh, N.; Misra, I.; Agrawal, A.; Devlin, J.; Girshick, R.; He, X.; Kohli, P.; Batra, D.; et al. 2016.  Visual storytelling.  In _Proceedings of the 2016 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_ , 1233–1239.
  * Jung et al. (2020) Jung, Y.; Kim, D.; Woo, S.; Kim, K.; Kim, S.; and Kweon, I. S. 2020.  Hide-and-Tell: Learning to Bridge Photo Streams for Visual Storytelling.  _arXiv preprint arXiv:2002.00774_ .
  * Kingma and Ba (2014) Kingma, D. P.; and Ba, J. 2014.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ .
  * Krishna et al. (2017) Krishna, R.; Zhu, Y.; Groth, O.; Johnson, J.; Hata, K.; Kravitz, J.; Chen, S.; Kalantidis, Y.; Li, L.-J.; Shamma, D. A.; et al. 2017.  Visual genome: Connecting language and vision using crowdsourced dense image annotations.  _International journal of computer vision_ 123(1): 32–73.
  * Lewis et al. (2019) Lewis, M.; Liu, Y.; Goyal, N.; Ghazvininejad, M.; Mohamed, A.; Levy, O.; Stoyanov, V.; and Zettlemoyer, L. 2019.  Bart: Denoising sequence-to-sequence pre-training for natural language generation, translation, and comprehension.  _arXiv preprint arXiv:1910.13461_ .
  * Li et al. (2015) Li, J.; Galley, M.; Brockett, C.; Gao, J.; and Dolan, B. 2015.  A diversity-promoting objective function for neural conversation models.  _arXiv preprint arXiv:1510.03055_ .
  * Li et al. (2019) Li, N.; Liu, B.; Han, Z.; Liu, Y.-S.; and Fu, J. 2019.  Emotion reinforced visual storytelling.  In _Proceedings of the 2019 on International Conference on Multimedia Retrieval_ , 297–305.
  * Lin (2004) Lin, C.-Y. 2004.  Rouge: A package for automatic evaluation of summaries.  In _Text summarization branches out_ , 74–81.
  * Liu and Singh (2004) Liu, H.; and Singh, P. 2004.  ConceptNet—a practical commonsense reasoning tool-kit.  _BT technology journal_ 22(4): 211–226.
  * Liu and Lapata (2019) Liu, Y.; and Lapata, M. 2019.  Hierarchical transformers for multi-document summarization.  _arXiv preprint arXiv:1905.13164_ .
  * Martin et al. (2017) Martin, L. J.; Ammanabrolu, P.; Wang, X.; Hancock, W.; Singh, S.; Harrison, B.; and Riedl, M. O. 2017.  Event representations for automated story generation with deep neural nets.  _arXiv preprint arXiv:1706.01331_ .
  * Mostafazadeh et al. (2017) Mostafazadeh, N.; Roth, M.; Louis, A.; Chambers, N.; and Allen, J. 2017.  Lsdsem 2017 shared task: The story cloze test.  In _Proceedings of the 2nd Workshop on Linking Models of Lexical, Sentential and Discourse-level Semantics_ , 46–51.
  * Papineni et al. (2002) Papineni, K.; Roukos, S.; Ward, T.; and Zhu, W.-J. 2002.  BLEU: a method for automatic evaluation of machine translation.  In _Proceedings of the 40th annual meeting of the Association for Computational Linguistics_ , 311–318.
  * Sood (2015) Sood, G. 2015.  _clarifai: R Client for the Clarifai API_.  R package version 0.2.
  * Vedantam, Lawrence Zitnick, and Parikh (2015) Vedantam, R.; Lawrence Zitnick, C.; and Parikh, D. 2015.  Cider: Consensus-based image description evaluation.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 4566–4575.
  * Veličković et al. (2017) Veličković, P.; Cucurull, G.; Casanova, A.; Romero, A.; Lio, P.; and Bengio, Y. 2017.  Graph attention networks.  _arXiv preprint arXiv:1710.10903_ .
  * Wang et al. (2020) Wang, R.; Wei, Z.; Li, P.; Zhang, Q.; and Huang, X. 2020.  Storytelling from an Image Stream Using Scene Graphs.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , 9185–9192.
  * Wang et al. (2018) Wang, X.; Chen, W.; Wang, Y.-F.; and Wang, W. Y. 2018.  No metrics are perfect: Adversarial reward learning for visual storytelling.  _arXiv preprint arXiv:1804.09160_ .
  * Yang et al. (2019) Yang, P.; Luo, F.; Chen, P.; Li, L.; Yin, Z.; He, X.; and Sun, X. 2019.  Knowledgeable Storyteller: A Commonsense-Driven Generative Model for Visual Storytelling.  In _IJCAI_ , 5356–5362.
  * Yao et al. (2019) Yao, L.; Peng, N.; Weischedel, R.; Knight, K.; Zhao, D.; and Yan, R. 2019.  Plan-and-write: Towards better automatic storytelling.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 33, 7378–7385.
  * Yu, Bansal, and Berg (2017) Yu, L.; Bansal, M.; and Berg, T. L. 2017.  Hierarchically-attentive rnn for album summarization and storytelling.  _arXiv preprint arXiv:1708.02977_ .
  * Zhu et al. (2015) Zhu, Y.; Kiros, R.; Zemel, R.; Salakhutdinov, R.; Urtasun, R.; Torralba, A.; and Fidler, S. 2015.  Aligning Books and Movies: Towards Story-Like Visual Explanations by Watching Movies and Reading Books.  In _The IEEE International Conference on Computer Vision (ICCV)_.


# Semi-convergence of the APSS method for a class of nonsymmetric three-by-
three singular saddle point problems

Hamed Aslani†, Davod Khojasteh Salkuyeh‡
†Faculty of Mathematical Sciences, University of Guilan, Rasht, Iran
‡Faculty of Mathematical Sciences, and Center of Excellence for Mathematical
Modelling,
Optimization and Combinational Computing (MMOCC), University of Guilan, Rasht,
Iran
Corresponding author.
Emails<EMAIL_ADDRESS>(H. Aslani<EMAIL_ADDRESS>(D.K.
Salkuyeh)

Abstract. For nonsymmetric block three-by-three singular saddle point problems
arising from the Picard iteration method for a class of mixed finite element
scheme, recently Salkuyeh et al. in (D.K. Salkuyeh, H. Aslani, Z.Z. Liang, An
alternating positive semi-definite splitting preconditioner for the three-by-
three block saddle point problems, Math. Commun. 26 (2021) 177-195)
established an alternating positive semi-definite splitting (APSS) method. In
this work, we analyse the semi-convergence of the APSS method for solving a
class of nonsymmetric block three-by-three singular saddle point problems. The
APSS induced preconditioner is applied to improve the semi-convergence rate of
the flexible GMRES (FGMRES) method. Experimental results are designated to
support the theoretical results. These results show that the served
preconditioner is efficient compared with FGMRES without a preconditioner.

Keywords: iterative methods, sparse matrices, saddle point, semi-convergence ,
preconditioning, Krylov methods.
AMS Subject Classification: 65F10, 65F50, 65F08.

##  1 Introduction

We are interested in solving the following large and sparse block three-by-
three saddle point system by iteration methods,

$\left(\begin{array}[]{ccc}{A}&{B^{T}}&{0}\\\ {B}&{0}&{C^{T}}\\\
{0}&{C}&{0}\end{array}\right)\left(\begin{array}[]{l}{x}\\\ {y}\\\
{z}\end{array}\right)=\left(\begin{array}[]{l}{f}\\\ {g}\\\
{h}\end{array}\right)$ (1)

where $A\in\mathbb{R}^{n\times n},B\in\mathbb{R}^{m\times n}$ and
$C\in\mathbb{R}^{l\times m}.$ Here $f\in\mathbb{R}^{n},g\in\mathbb{R}^{m}$ and
$h\in\mathbb{R}^{l}$ are known and ${\bf x}=\left(x;y;z\right)$ is an unknown
vector that has to be determined. The coefficient matrix of the system (2) is
of order $\textbf{n}=m+n+l$. Many practical applications produce linear
systems of the form (2); for example in application of the Picard iteration
method for a class of mixed finite element scheme for stationary
magnetohydrodynamics models [17] and finite element method to solve the time-
dependent Maxwell equations having discontinuous coefficients in polyhedral
domains with Lipschitz boundary. References [26, 14] provide additional
instances and references therein.

In this paper, we make an attempt to solve the following system,

$\mathcal{A}{\bf x}\equiv\left(\begin{array}[]{ccc}{A}&{B^{T}}&{0}\\\
{-B}&{0}&{-C^{T}}\\\
{0}&{C}&{0}\end{array}\right)\left(\begin{array}[]{l}{x}\\\ {y}\\\
{z}\end{array}\right)=\left(\begin{array}[]{l}{f}\\\ {g}\\\
{h}\end{array}\right)\equiv\mathbf{b},$ (2)

which is an equivalent form of (1). Note that the coefficient matrix in (1) is
symmetric, however, $\mathcal{A}$ is nonsymmetric. Nevertheless, $\mathcal{A}$
has some good properties. For instance, $\mathcal{A}$ is positive semi-
definite. It means that $\mathcal{A}+\mathcal{A}^{T}$ is Symmetric Positive
Semi-Definite (SPSD). This property can greatly improve the performance of the
GMRES method for solving the system. Some other features are available in
[20].

Consider the two-by-two block saddle point problem

$\mathcal{\tilde{A}}{\bf\tilde{x}}\equiv\left(\begin{array}[]{ccc}\tilde{A}&\tilde{B}^{T}\\\
-\tilde{B}&0\end{array}\right)\left(\begin{array}[]{l}\tilde{x}\\\
\tilde{y}\end{array}\right)=\left(\begin{array}[]{l}\tilde{f}\\\
\tilde{g}\end{array}\right)\equiv\mathbf{\tilde{b}}.$ (3)

The matrix $\tilde{A}\in\mathbb{R}^{p\times p}$ is assumed to be a Symmetric
Positive Definite (SPD), $\tilde{B}\in\mathbb{R}^{p\times q}$ with $p>q$ is a
matrix that $rank(B)=r<q$ (i.e., $\tilde{B}$ is a rank-deficient matrix),
$\tilde{f}\in\mathbb{R}^{p}$ and $\tilde{g}\in\mathbb{R}^{q}$. Thereby, the
linear system (3) leads to a two-by-two singular saddle point problem. In the
context of the two-by-two singular saddle point problem, Bai [4] established
the HSS method for singular saddle point problem and also derived some
conditions for guaranteeing the semi-convergence of the HSS method. Li et al.
[22] generalized the HSS method for solving non-Hermitian, positive semi-
definite, and singular linear systems. They studied semi-convergence analysis
of the Generalized HSS (GHSS) method. In addition, an upper bound for the
semi-convergence factor was derived. In [12] a Generalized Preconditioned
Hermitian and skew-Hermitian splitting method (GPHSS) was considered to solve
the singular saddle point problems. The semi-convergence of the GPHSS scheme
was proved under some conditions. In addition, the authors obtained the
induced preconditioner and discussed the eigenvalues of the preconditioned
matrix. Then, the Local Hermitian and skew-Hermitian (LHSS) method and the
Modified LHSS (MLHSS) method were established in [20]. They also gave the
semi-convergence conditions. Motivated by the Uzawa method, Chao et al. [11]
designed the Uzawa-SOR method for singular saddle point problems.
Generalization of the Uzawa-SOR method was investigated in [31]. The authors
established the Uzawa-AOR scheme to solve (3). The distribution of the
eigenvalues of the iteration matrix and the semi-convergence properties were
given. Numerical results indicate that the Uzawa-AOR method outperforms the
Uzawa method [2], the parameterized Uzawa method [35], and the Uzawa-SOR
method [11]. Some other efficient methods for the singular saddle point
problem of the form (3) were studied in [24, 34, 33].

Note that some partitioning techniques can be employed to represent the three-
by-three coefficient matrix in (2) as:

$\left(\begin{array}[]{cccc}A&B^{T}&\vdots&0\\\ B&0&\vdots&C^{T}\\\
\cdots&\cdots&\vdots&\cdots\\\
0&C&\vdots&0\end{array}\right)\qquad\text{or}\qquad\left(\begin{array}[]{cccc}A&\vdots&B^{T}&0\\\
\cdots&\vdots&\cdots&\cdots\\\ B&\vdots&0&C^{T}\\\
0&\vdots&C&0\end{array}\right),$

which are the standard two-by-two saddle point problem.

Many available preconditioning schemes for (3), can not be implemented for
solving (1). This is because of different properties in two problems. Note
that the matrix $\mathcal{A}$ described in (2) is nonsingular if $A$ is SPD
and the matrices $B$ and $C$ are of full row rank (see [18, 19, 28]). In
recent years, many researchers have considered the three-by-three saddle point
problem. Huang et al. [20] arranged a block diagonal preconditioner for
solving the nonsingular system of the form (2). The exact and inexact versions
of the preconditioner were also studied. Then, in [10] the shift splitting
(SS) and the relaxed shift splitting (RSS) method were designated. Xie et al.
[30] considered three efficient preconditioners. The authors analyzed the
eigenvalues of the corresponding preconditioned matrices. Aslani et al. [3]
presented a new method for solving (2) when $A$ is SPD and $B,C$ are full row
rank matrices. Convergence properties of the method was derived. Moreover, the
spectral properties of the preconditioned matrix were discussed. Abdolmaleki
et al. [1] from another perspective organized a block three-by-three diagonal
preconditioner for (2). A suitable estimation strategy for lower and upper
bounds of eigenvalues of the preconditioned matrix was considered. In [29] an
exact parameterized block SPD preconditioner and its inexact version for a
class of block three-by-three saddle point problems was proposed. The authors
also estimated the eigenvalue bounds for the preconditioned matrix.

In [21], Liang and Zhang proposed the Alternating Positive Semi-definite
Splitting (APSS) iteration method for double saddle point problems. Using the
idea of [21], Salkuyeh et al. in [27] applied the APSS method for solving the
problem (2) and proved its convergence. In the case of $C$ being rank-
deficient, the coefficient matrix (2) is singular. Accordingly, the linear
system (2) is labeled as a singular three-by-three saddle point problem. In
this work, the three-by-three large, sparse, and singular saddle point problem
is considered and the semi-convergence analysis of the APSS method is
discussed.

The structure of this paper is as follows. The paper is started with a review
of the APSS method and its corresponding induced preconditioner. In Section 2,
we focus on the semi-convergence properties of the APSS method for solving
(2). Unconditionally semi-convergent for the APSS iteration method are derived
in Section 3. A strategy is given to estimate the parameter of the proposed
method in Section 4. Section 5 is devoted to give some numerical tests to
support the theoretical results. The end of the paper will be companied by
some brief conclusions.

##  2 Review of the APSS method

Let us first give a brief study of the APSS method. Consider the following
decomposition for the coefficient matrix $\mathcal{A}$ in (2):

$\mathcal{A}=\mathcal{A}_{1}+\mathcal{A}_{2},$ (4)

where

$\mathcal{A}_{1}=\left(\begin{array}[]{ccc}{A}&{B^{T}}&{0}\\\ {-B}&{0}&{0}\\\
{0}&{0}&{0}\end{array}\right),\quad\mathcal{A}_{2}=\left(\begin{array}[]{ccc}{0}&{0}&{0}\\\
{0}&{0}&{-C^{T}}\\\ {0}&{C}&{0}\end{array}\right).$ (5)

Let $\alpha>0$ be a given parameter. Based on the decomposition (4), the
following splittings for the matrix $\mathcal{A}$ can be stated

$\mathcal{A}=(\alpha\mathcal{I}+\mathcal{A}_{1})-(\alpha\mathcal{I}-\mathcal{A}_{2})=(\alpha\mathcal{I}+\mathcal{A}_{2})-(\alpha\mathcal{I}-\mathcal{A}_{1}),$

where $\mathcal{I}$ is the identity matrix of order $\mathbf{n}$. Now, by
using these splittings the APSS method can be written as

$\displaystyle(\alpha\mathcal{I}+\mathcal{A}_{1})x^{\left(k+\frac{1}{2}\right)}=(\alpha\mathcal{I}-\mathcal{A}_{2})x^{(k)}+b,$
$\displaystyle(\alpha\mathcal{I}+\mathcal{A}_{2})x^{(k+1)}=(\alpha\mathcal{I}-\mathcal{A}_{1})x^{\left(k+\frac{1}{2}\right)}+b,$

where $x^{(0)}\in\mathbb{R}^{\mathbf{n}}$ is an initial guess. By eliminating
$x^{(k+\frac{1}{2})},$ the iteration scheme can be rewritten as the stationary
form

$x^{k+1}=\mathcal{T_{\alpha}}x^{k}+f,$ (6)

with

$\mathcal{T_{\alpha}}=(\alpha\mathcal{I}+\mathcal{A}_{2})^{-1}(\alpha\mathcal{I}-\mathcal{A}_{1})(\alpha\mathcal{I}+\mathcal{A}_{1})^{-1}(\alpha\mathcal{I}-\mathcal{A}_{2}),$
(7)

and

$f=2\alpha(\alpha\mathcal{I}+\mathcal{A}_{2})^{-1}(\alpha\mathcal{I}+\mathcal{A}_{1})^{-1}b.$

Similar to the Hermitian and Skew-Hermitian splitting (HSS) iteration method
[5], if we set

$\tilde{\mathcal{M}_{\alpha}}=\frac{1}{2\alpha}(\alpha\mathcal{I}+\mathcal{A}_{1})(\alpha\mathcal{I}+\mathcal{A}_{2}),\quad\tilde{\mathcal{N}_{\alpha}}=\frac{1}{2\alpha}(\alpha\mathcal{I}-\mathcal{A}_{1})(\alpha\mathcal{I}-\mathcal{A}_{2}),$

then $\mathcal{A}=\tilde{\mathcal{M}_{\alpha}}-\tilde{\mathcal{N}_{\alpha}}$
and

$\mathcal{T}_{\alpha}=\tilde{\mathcal{M}_{\alpha}}^{-1}\tilde{\mathcal{N}_{\alpha}}=\mathcal{I}-\tilde{\mathcal{M}_{\alpha}}^{-1}{\mathcal{A}}.$

From now on, we use
$\mathcal{M}_{\alpha}=(\alpha\mathcal{I}+\mathcal{A}_{1})(\alpha\mathcal{I}+\mathcal{A}_{2})$
as the APSS preconditioner, since the pe-factor $\frac{1}{2\alpha}$ has no
effect on the preconditioned matrix. So, the saddle point system (2) can be
preconditioned from the left as ${\mathcal{M}_{\alpha}}^{-1}\mathcal{A}{\bf
x}={\mathcal{M}_{\alpha}}^{-1}{\bf b}$. In this case, we have

${\mathcal{M}_{\alpha}}^{-1}\mathcal{A}{\bf
x}=(\mathcal{I}-{\mathcal{M}_{\alpha}}^{-1}{\mathcal{N}_{\alpha}}){\bf
x}={\mathcal{M}_{\alpha}}^{-1}{\bf b},$ (8)

where

${\mathcal{N}_{\alpha}}=(\alpha\mathcal{I}-\mathcal{A}_{1})(\alpha\mathcal{I}-\mathcal{A}_{2})=\alpha\mathcal{I}.\alpha\mathcal{I}-\alpha(\mathcal{A}_{1}+\mathcal{A}_{2})+\mathcal{A}_{1}\mathcal{A}_{2}.$
(9)

##  3 The semi-convergence of APSS iteration method

In this section, we will analyze the semi-convergence properties of the APSS
iteration method for solving the double saddle point problem (2). First, let
us give some related definitions, lemmas and theorems.

######  Definition 1.

The iteration method (6) is said to be semiconvergent if for any initial guess
${\bf{x}}^{(0)}$, the iteration sequence $\\{\bf{x}^{(k)}\\}$ produced by (6)
converges to a solution of linear systems $\mathcal{A}\bf{x}=\bf{b}.$

######  Theorem 1.

[7] The iteration method (6) is semi-convergent if and only if

$index(\mathcal{I}-\mathcal{T}_{\alpha})=1\quad\text{and}\quad\vartheta(\mathcal{T}_{\alpha})<1,$

where

$\vartheta(\mathcal{T}_{\alpha})\equiv\max\\{|\lambda|:\lambda\in\sigma(\mathcal{T}_{\alpha}),\lambda\neq
1\\}<1,$

in which $\vartheta(\mathcal{T}_{\alpha})$ is called the pseudo-spectral
radius of $\mathcal{T}_{\alpha}.$ Moreover, the limit of the sequence can be
written as

${\bf{x}}^{(*)}=(\mathcal{I}-\mathcal{T}_{\alpha})^{D}f+[\mathcal{I}-(\mathcal{I}-\mathcal{T}_{\alpha})^{D}(\mathcal{I}-\mathcal{T}_{\alpha})]{\bf{x}}^{(0)},$

where $(\mathcal{I}-{\mathcal{T}_{\alpha}})^{D}$ denotes the Drazin inverse of
$(\mathcal{I}-{\mathcal{T}_{\alpha}})$ [8].

######  Lemma 1.

[25] (Kellogg’s lemma) If $A\in\mathbb{C}^{n\times n}$ is positive semi-
definite, then

$\left\|(\alpha I+A)^{-1}(\alpha I-A)\right\|_{2}\leq 1,$

for all $\alpha>0$. Moreover, if $A\in\mathbb{C}^{n\times n}$ is positive
definite, then

$\left\|(\alpha I+A)^{-1}(\alpha I-A)\right\|_{2}<1,$

for all $\alpha>0.$

######  Remark 1.

For the saddle point matrix of the form

$\mathcal{A}=\left(\begin{array}[]{ccc}{A}&{B^{T}}&{0}\\\ {-B}&{0}&{-C^{T}}\\\
{0}&{C}&{0}\end{array}\right),$

if $\mathcal{A}$ is singular, we can easily see that at least one of the sets
$null(B^{T})\cap null(C)$ and $null(C^{T})$ is nontrivial, i.e., the dimension
of one of the two sets is at least one.

Noticing that

$\displaystyle(\alpha\mathcal{I}-\mathcal{A}_{1})(\alpha\mathcal{I}+\mathcal{A}_{1})^{-1}$
$\displaystyle=$
$\displaystyle(\alpha\mathcal{I}+\mathcal{A}_{1})^{-1}(\alpha\mathcal{I}-\mathcal{A}_{1})$
(10)
$\displaystyle(\alpha\mathcal{I}-\mathcal{A}_{2})(\alpha\mathcal{I}+\mathcal{A}_{2})^{-1}$
$\displaystyle=$
$\displaystyle(\alpha\mathcal{I}+\mathcal{A}_{2})^{-1}(\alpha\mathcal{I}-\mathcal{A}_{2}),$
(11)

it is easy to see that the matrix $\mathcal{T}_{\alpha}$ is similar to

$\mathcal{L_{\alpha}}=(\alpha\mathcal{I}+\mathcal{A}_{1})^{-1}(\alpha\mathcal{I}-\mathcal{A}_{1})(\alpha\mathcal{I}+\mathcal{A}_{2})^{-1}(\alpha\mathcal{I}-\mathcal{A}_{2}).$
(12)

Now, since the matrices $\mathcal{A}_{1}$ and $\mathcal{A}_{2}$ are both
positive semi-definite, then using the Kellogg’s lemma we have

$\|\mathcal{L_{\alpha}}\|_{2}\leq\|(\alpha\mathcal{I}+\mathcal{A}_{1})^{-1}(\alpha\mathcal{I}-\mathcal{A}_{1})\|_{2}\|(\alpha\mathcal{I}+\mathcal{A}_{2})^{-1}(\alpha\mathcal{I}-\mathcal{A}_{2})\|_{2}\leq
1.$

Thus, it holds that $\mathcal{L}_{\alpha}\mathbf{x}=\mathbf{x}$ if and only if
${}\mathcal{L}_{\alpha}^{*}\mathbf{x}=\mathbf{x},$ for any
$\mathbf{x}\in\mathbb{C}^{n\times n}$ [4]. As a result the index of
$\text{index}(\mathcal{I}-\mathcal{L}_{\alpha})=1$. Eventually, since two
similar matrices have the same index, we see that

$\text{index}(\mathcal{I}-\mathcal{T}_{\alpha})=1.$ (13)

Next, we give the conditions for $\vartheta(\mathcal{T}_{\alpha})<1.$

######  Theorem 2.

Suppose that $A\in\mathbb{R}^{n\times n}$ is symmetric positive definite,
$B\in\mathbb{R}^{m\times n}$ is an arbitrary matrix and
$C\in\mathbb{R}^{l\times m}$ is rank-deficient. Then, for the decomposition
(4) of the matrix $\mathcal{A}$, the pseudo-spectral radius of
$\mathcal{T}_{\alpha}$ is less than one, i.e.,
$\vartheta(\mathcal{T}_{\alpha})<1,\forall\alpha>0.$

###### Proof.

First of all note that the matrix $\alpha\mathcal{I}-\mathcal{A}_{2}$ is non-
singular, and the equalities 10 and 11 hold. Hence, the matrix
$\mathcal{T}_{\alpha}$ is similar to

$\mathcal{P}_{\alpha}=(\alpha\mathcal{I}+\mathcal{A}_{2})^{-1}(\alpha\mathcal{I}-\mathcal{A}_{2})(\alpha\mathcal{I}+\mathcal{A}_{1})^{-1}(\alpha\mathcal{I}-\mathcal{A}_{1}),$

So, $\vartheta(\mathcal{T}_{\alpha})=\vartheta(\mathcal{P}_{\alpha}).$ Let
$x\in\mathbb{C}^{n}$ be any eigenvector of the matrix $\mathcal{P}_{\alpha}$
and $\lambda$ be the eigenvalue of matrix $\mathcal{P}_{\alpha}$ corresponding
to eigenvector $x$, i.e., $\mathcal{P}_{\alpha}x=\lambda x.$ Without loss of
generality, we assume that $\|x\|_{2}=1.$ Now, we prove the theorem according
to the following four cases:

Case 1. $x\in null(\mathcal{A}_{1})\cap null(\mathcal{A}_{2}).$ So, we have
$\mathcal{A}_{1}x=\mathcal{A}_{2}x=0$ and it follows that

$(\alpha\mathcal{I}+\mathcal{A}_{1})x=(\alpha\mathcal{I}-\mathcal{A}_{1})x\Longrightarrow
x=(\alpha\mathcal{I}+\mathcal{A}_{1})^{-1}(\alpha\mathcal{I}-\mathcal{A}_{1})x,$
(14)

and

$(\alpha\mathcal{I}+\mathcal{A}_{2})x=(\alpha\mathcal{I}-\mathcal{A}_{2})x\Longrightarrow
x=(\alpha\mathcal{I}+\mathcal{A}_{2})^{-1}(\alpha\mathcal{I}-\mathcal{A}_{2})x.$
(15)

By combining (14) and (15), we have $\mathcal{P_{\alpha}}x=x.$ Thus,
$\lambda=1.$

Case 2. $x\in null(\mathcal{A}_{2}),$ but $x\notin null(\mathcal{A}_{1}).$ It
means that $\mathcal{A}_{1}x\neq 0$ and $\mathcal{A}_{2}x=0.$ From
$\mathcal{P}_{\alpha}x=\lambda x,$ by easy computations we can obtain

$\lambda
x=(\alpha\mathcal{I}+\mathcal{A}_{1})^{-1}(\alpha\mathcal{I}-\mathcal{A}_{1})x.$
(16)

Since $\mathcal{A}_{1}$ is positive semi-definite, it follows from the
Kellogg’s lemma that

$\|(\alpha\mathcal{I}+\mathcal{A}_{1})^{-1}(\alpha\mathcal{I}-\mathcal{A}_{1}\|_{2}\leq
1.$

Therefore $|\lambda|\leq 1.$ In the following, we further prove that
$|\lambda|<1$ for any $\alpha>0.$ We will argue it by contradiction. If
$|\lambda|=1,$ then there exists $\theta\in(-\pi,\pi]$ such that

$(\alpha\mathcal{I}-\mathcal{A}_{1})(\alpha\mathcal{I}+\mathcal{A}_{1})^{-1}x=e^{\imath\theta}x.$
(17)

If we set
$\mathcal{V}_{\alpha}=(\alpha\mathcal{I}-\mathcal{A}_{1})(\alpha\mathcal{I}+\mathcal{A}_{1})^{-1},$
then we can rewritten (17) as follows

$\mathcal{V}_{\alpha}x=e^{\imath\theta}x.$

It then follows that

$\|\mathcal{V}_{\alpha}x\|_{2}=\|x\|_{2}.$ (18)

Now, letting

$w:=(u;v;p)=(\alpha\mathcal{I}+\mathcal{A}_{1})^{-1}x,\quad\text{ with }\quad
u\in\mathbb{C}^{n},\quad v\in\mathbb{C}^{m},\text{ and }\quad
p\in\mathbb{C}^{l},$

we get

$\|(\alpha\mathcal{I}-\mathcal{A}_{1})w\|_{2}=\|(\alpha\mathcal{I}+\mathcal{A}_{1})w\|_{2}.$
(19)

From (19), it holds that

$w^{*}(\mathcal{A}_{1}+\mathcal{A}_{1}^{*})w=0,$

or equivalently,

$\left(\begin{array}[]{ccc}u^{*}&v^{*}&p^{*}\end{array}\right)\left(\begin{array}[]{ccc}{A}&{0}&{0}\\\
{0}&{0}&{0}\\\ {0}&{0}&{0}\end{array}\right)\left(\begin{array}[]{l}u\\\ v\\\
p\end{array}\right)=0,$

which leads $u=0$, due to the symmetric positive definiteness of $A.$ On one
hand, from (17) we have

$(\alpha\mathcal{I}-\mathcal{A}_{1})w=e^{\imath\theta}(\alpha\mathcal{I}+\mathcal{A}_{1})w,$

which gives

$\left(\begin{array}[]{ccc}\alpha I-A&-B^{T}&{0}\\\ B&\alpha I&{0}\\\
{0}&{0}&\alpha I\\\ \end{array}\right)\left(\begin{array}[]{c}0\\\ v\\\
p\end{array}\right)=e^{\imath\theta}\left(\begin{array}[]{ccc}\alpha
I+A&B^{T}&{0}\\\ -B&\alpha I&{0}\\\ {0}&{0}&\alpha I\\\
\end{array}\right)\left(\begin{array}[]{l}0\\\ v\\\ p\end{array}\right),$

and so

$\left(\begin{array}[]{c}-B^{T}v\\\ \alpha v\\\ \alpha
p\end{array}\right)=e^{\imath\theta}\left(\begin{array}[]{c}B^{T}v\\\ \alpha
v\\\ \alpha p\end{array}\right).$

It leads to

$\displaystyle B^{T}v=-e^{\imath\theta}B^{T}v,$ (20) $\displaystyle\alpha
v=e^{\imath\theta}\alpha v,$ (21) $\displaystyle\alpha
p=e^{\imath\theta}\alpha p.$ (22)

Here, if $B^{T}v=0,$ then we see that $v=0$, because $B$ is a full row rank
matrix. It is immediate to conclude that $\mathcal{A}_{1}x=0,$ which is a
contradiction with the assumption $x\notin null(\mathcal{A}_{1}).$ Thereby,
from (20) it holds $e^{\imath\theta}=-1.$ Substituting the equality
$e^{\imath\theta}=-1$ into (21) and (22), results in $v=0$ and $p=0,$
respectively. So $w=0.$ Consequently, we have
$x=(\alpha\mathcal{I}+\mathcal{A}_{1})w=0,$ which contradicts with the fact
that $x$ is an eigenvector. Therefore, $|\lambda|<1.$

Case 3. $x\in null(\mathcal{A}_{1}),$ but $x\notin null(\mathcal{A}_{2}).$ So,
we have

$\mathcal{A}_{1}x=0\Longrightarrow(\alpha\mathcal{I}+\mathcal{A}_{1})x=\alpha
x\Longrightarrow(\alpha\mathcal{I}+\mathcal{A}_{1})^{-1}x=\frac{1}{\alpha}x.$

Hence, using $\mathcal{A}_{1}x=0$ and the above equation, we get

$\displaystyle\mathcal{P}_{\alpha}x$ $\displaystyle=$
$\displaystyle\alpha(\alpha\mathcal{I}+\mathcal{A}_{2})^{-1}(\alpha\mathcal{I}-\mathcal{A}_{2})(\alpha\mathcal{I}+\mathcal{A}_{1})^{-1}x$
$\displaystyle=$
$\displaystyle(\alpha\mathcal{I}+\mathcal{A}_{2})^{-1}(\alpha\mathcal{I}-\mathcal{A}_{2})x.$

Therefore, using the fact that $\mathcal{P}_{\alpha}x=\lambda x$, positive
semi-definiteness of $\mathcal{A}_{2}$ and the Kellogg’s lemma, we deduce that

$|\lambda|\leq\|\mathcal{P}_{\alpha}x\|\leq\|(\alpha\mathcal{I}+\mathcal{A}_{2})^{-1}(\alpha\mathcal{I}-\mathcal{A}_{2})\|_{2}\leq
1.$ (23)

Moreover, we claim that $|\lambda|=1$ never happens. By contradiction we
assume that $|\lambda|=1.$ So, there exists $\theta\in(-\pi,\pi]$ so that

$(\alpha\mathcal{I}-\mathcal{A}_{2})x=e^{\imath\theta}(\alpha\mathcal{I}+\mathcal{A}_{2})x,$

which is equivalent to

$\displaystyle\alpha x_{1}=e^{\imath\theta}\alpha x_{1},$ (24)
$\displaystyle\alpha x_{2}+C^{T}x_{3}=e^{\imath\theta}(\alpha
x_{2}-C^{T}x_{3}),$ (25) $\displaystyle-Cx_{2}+\alpha
x_{3}=e^{\imath\theta}(Cx_{2}+\alpha x_{3}).$ (26)

From Eq. (24), either $e^{\imath\theta}=1$ or $x_{1}=0.$ First, suppose that
$e^{\imath\theta}=1.$ From $\mathcal{A}_{1}x=0,$ we have

$\displaystyle Ax_{1}+B^{T}x_{2}=0,$ (27) $\displaystyle Bx_{2}=0.$ (28)

Substituting $x_{1}=-A^{-1}B^{T}x_{2}$ which is deduced from (27) into (28),
gives $BA^{-1}B^{T}x_{2}=0,$ consequently, $x_{2}=0,$ and then by (27),
$x_{1}=0.$ Now, from Eq. (25), we have $C^{T}x_{3}=0.$ Thus,
$\mathcal{A}_{2}x=0$ which is in contradiction with the assumption. Therefore,
$e^{\imath\theta}=1$ never happens. Now, we discuss $x_{1}=0.$ Clearly, the
assumption $\mathcal{A}_{1}x=0,$ gives

$\left(\begin{array}[]{ccc}{A}&{B^{T}}&{0}\\\ {-B}&{0}&{0}\\\
{0}&{0}&{0}\end{array}\right)\left(\begin{array}[]{l}0\\\ x_{2}\\\
x_{3}\end{array}\right)=\left(\begin{array}[]{l}B^{T}x_{2}\\\ 0\\\
0\end{array}\right)=0.$

Since $B$ has full row rank, we deduce that $x_{2}=0.$ Substituting $x_{2}=0$
into (26) gives $x_{3}=0.$ So, $x=(x_{1};x_{2};x_{3})=0,$ which is impossible.
thereupon, $|\lambda|=1$ is unacceptable.

Case 4. $x\notin null(\mathcal{A}_{1}),x\notin null(\mathcal{A}_{2}).$ For
this case, we have $\mathcal{A}_{1}x\neq 0$ and $\mathcal{A}_{2}x\neq 0.$ If
we define
$\mathcal{Q}_{\alpha}=(\alpha\mathcal{I}+\mathcal{A}_{2})^{-1}(\alpha\mathcal{I}-\mathcal{A}_{2})$
and
$\mathcal{V}_{\alpha}=(\alpha\mathcal{I}+\mathcal{A}_{1})^{-1}(\alpha\mathcal{I}-\mathcal{A}_{1}),$
we can see that
$\mathcal{P}_{\alpha}=\mathcal{Q}_{\alpha}\mathcal{V}_{\alpha}.$ Since
$\mathcal{Q}_{\alpha}$ is a unitary matrix, we see that
$\left\|\mathcal{Q}_{\alpha}\right\|_{2}=1$. On the other hand, since
$\mathcal{A}_{1}$ is positive semi-definite matrix from the Kellogg’s lemma,
we have $\left\|\mathcal{V}_{\alpha}\right\|_{2}\leq 1$ for any $\alpha>0.$ It
leads to

$\|\mathcal{P}_{\alpha}\|_{2}\leq\|\mathcal{Q}_{\alpha}\|_{2}\|\mathcal{V}_{\alpha}\|_{2}=\|\mathcal{V}_{\alpha}\|_{2}\leq
1,\quad\forall\alpha>0.$

Thus, $|\lambda|\leq 1.$ In the following, we further prove that
$|\lambda|<1,$ for any $\alpha>0.$ We will argue it by contradiction. If
$|\lambda|=1,$ then there exists $\theta\in(-\pi,\pi]$ such that

$\mathcal{P}_{\alpha}x=e^{\imath\theta}x,$

which is equivalent to,

$\mathcal{V}_{\alpha}x=e^{\imath\theta}\mathcal{Q}_{\alpha}^{*}x.$ (29)

Cosquently,

$\|\mathcal{V}_{\alpha}x\|_{2}=\|x\|_{2}.$ (30)

Substituting $\mathcal{V}_{\alpha}$ into (30) and using the change of variable

$w:=(u;v;p)=(\alpha\mathcal{I}+\mathcal{A}_{1})^{-1}x,\quad\text{ with }\quad
u\in\mathbb{C}^{n}\quad,v\in\mathbb{C}^{m}\text{ and }\quad
p\in\mathbb{C}^{l},$

gives

$\|(\alpha\mathcal{I}-\mathcal{A}_{1})w\|_{2}=\|(\alpha\mathcal{I}+\mathcal{A}_{1})w\|_{2}.$
(31)

From (31), it holds

$w^{*}(\mathcal{A}_{1}+\mathcal{A}_{1}^{*})w=0,$

or equivalently,

$\left(\begin{array}[]{ccc}u^{*}&v^{*}&p^{*}\end{array}\right)\left(\begin{array}[]{ccc}{A}&{0}&{0}\\\
{0}&{0}&{0}\\\ {0}&{0}&{0}\end{array}\right)\left(\begin{array}[]{l}u\\\ v\\\
p\end{array}\right)=0,$

which leads to $u=0$, due to the symmetric positive definiteness of $A.$ On
one hand, from (29) we have

$(\alpha\mathcal{I}-\mathcal{A}_{1})w=e^{\imath\theta}\mathcal{Q}_{\alpha}^{*}(\alpha\mathcal{I}+\mathcal{A}_{1})w,$

which results in

$\left(\begin{array}[]{ccc}\alpha I-A&-B^{T}&{0}\\\ B&\alpha I&{0}\\\
{0}&{0}&\alpha I\\\ \end{array}\right)\left(\begin{array}[]{c}0\\\ v\\\
p\end{array}\right)=e^{\imath\theta}\mathcal{Q}_{\alpha}^{*}\left(\begin{array}[]{ccc}\alpha
I+A&B^{T}&{0}\\\ -B&\alpha I&{0}\\\ {0}&{0}&\alpha I\\\
\end{array}\right)\left(\begin{array}[]{l}0\\\ v\\\ p\end{array}\right).$

So

$\left(\begin{array}[]{c}-B^{T}v\\\ \alpha v\\\ \alpha
p\end{array}\right)=e^{\imath\theta}\mathcal{Q}_{\alpha}^{*}\left(\begin{array}[]{c}B^{T}v\\\
\alpha v\\\ \alpha p\end{array}\right).$

The above equality can be rewritten as

$(\alpha\mathcal{I}-\mathcal{A}_{2})\left(\begin{array}[]{c}-B^{T}v\\\ \alpha
v\\\ \alpha
p\end{array}\right)=e^{\imath\theta}(\alpha\mathcal{I}+\mathcal{A}_{2})\left(\begin{array}[]{c}B^{T}v\\\
\alpha v\\\ \alpha p\end{array}\right).$

It leads to

$\left(\begin{array}[]{ccc}\alpha I&{0}&{0}\\\ {0}&\alpha I&C^{T}\\\
{0}&-C&\alpha I\\\ \end{array}\right)\left(\begin{array}[]{c}-B^{T}v\\\ \alpha
v\\\ \alpha
p\end{array}\right)=e^{\imath\theta}\left(\begin{array}[]{ccc}\alpha
I&{0}&{0}\\\ {0}&\alpha I&-C^{T}\\\ {0}&C&\alpha I\\\
\end{array}\right)\left(\begin{array}[]{l}B^{T}v\\\ \alpha v\\\ \alpha
p\end{array}\right),$

that is implies that

$\displaystyle B^{T}v=-e^{\imath\theta}B^{T}v,$ (32) $\displaystyle\alpha
v+C^{T}p=e^{\imath\theta}(\alpha v-C^{T}p),$ (33) $\displaystyle-Cv+\alpha
p=e^{\imath\theta}(Cv+p).$ (34)

If $B^{T}v=0,$ then we can easily see that $\mathcal{A}_{1}v=0.$ It is
immediate to conclude that $\mathcal{A}_{1}x=0,$ which is a contradiction with
the consideration $x\notin null(\mathcal{A}_{1}).$ Thereby, from (32), it
holds $e^{\imath\theta}=-1.$ Substituting the identity $e^{\imath\theta}=-1$
into (33) and (34), derives $v=0$ and $p=0,$ respectively. So $w=0.$
Consequently, we have $x=(\alpha\mathcal{I}+\mathcal{A}_{1})w=0,$ which
contradicts with the fact that $x$ is a non-zero vector. Therefore,
$|\lambda|<1.$

In summary, from the Cases 1-4, we see that
$\vartheta(\mathcal{T}_{\alpha})<1.$ ∎

######  Theorem 3.

Suppose that the assumptions of Theorem 2 hold. Then, the APSS iteration
method (6) is semi-convergent for any $\alpha>0$.

###### Proof.

The proof immediately follows from Eq. (13) and Theorem 2. ∎

##  4 Estimation strategy for the parameter $\alpha$

Finding the optimal parameter $\alpha$ of the APSS method is generally hard.
In this section, we propose an appropriate strategy for estimating $\alpha$ in
the APSS method, which has been studied by Cao [9].

Notably, from (8) it is anticipated that $\mathcal{M}_{\alpha}$ is as close as
conceivable to $\mathcal{A}$ when $\mathcal{N}_{\alpha}\approx 0.$ In this
way, having Eq. (9) in mind, the function

$\displaystyle\Psi(\alpha)$
$\displaystyle=\alpha\|\mathcal{I}\|_{F}\cdot\alpha\|\mathcal{I}\|_{F}-\alpha\left(\left\|\mathcal{A}_{1}\right\|_{F}+\left\|\mathcal{A}_{2}\right\|_{F}\right)+\left\|\mathcal{A}_{1}\right\|_{F}\left\|\mathcal{A}_{2}\right\|_{F}$
$\displaystyle=(n+m+l)\alpha^{2}-\alpha\left(\left\|\mathcal{A}_{1}\right\|_{F}+\left\|\mathcal{A}_{2}\right\|_{F}\right)+\left\|\mathcal{A}_{1}\right\|_{F}\left\|\mathcal{A}_{2}\right\|_{F}$
$\displaystyle={\bf{n}}\alpha^{2}-\alpha\left(\left\|\mathcal{A}_{1}\right\|_{F}+\left\|\mathcal{A}_{2}\right\|_{F}\right)+\left\|\mathcal{A}_{1}\right\|_{F}\left\|\mathcal{A}_{2}\right\|_{F},$

can be characterized. Minimizing $\Psi(\alpha)$ with respect to $\alpha$ leads
to the estimation parameter $\alpha_{est}$ being set to

$\alpha_{est}=\frac{\left\|\mathcal{A}_{1}\right\|_{F}+\left\|\mathcal{A}_{2}\right\|_{F}}{2{\bf{n}}}.$

Since the matrix $A$ is SPD (in this case, $\mathcal{A}_{1}\neq 0$), one can
conclude that $\alpha_{est}>0.$ In the following section, the efficiency of
this choice will be verified.

##  5 Numerical experiments

To indicate the efficiency of the APSS preconditioner, we have conducted some
numerical tests. We provide two examples and each one is scaled by a symmetric
diagonal matrix. To this end, the coefficient matrix $\mathcal{A}$ is replaced
by the matrix
$\mathcal{D}^{-\frac{1}{2}}\mathcal{A}\mathcal{D}^{-\frac{1}{2}},$ in which
${\cal D}=\textrm{diag}(\|{\mathcal{A}}_{1}\|_{2},\ldots,\|{\mathcal{A}}_{\bf
n}\|_{2})$. In addition, the $j$th column of the matrix $\mathcal{A}$ is
represented by $\mathcal{A}_{j}.$

For all the examples, a zero vector was used as the initial guess, and the
iteration was terminated once

$Res=\frac{\left\|{\mathbf{b}}-\mathcal{A}{\bf{x}}^{(k)}\right\|_{2}}{\left\|{\mathbf{b}}\right\|_{2}}\leq
10^{-7},$

or the specified number of iteration steps, $maxit=2000$, is exceeded. Note
that ${\bf{x}}^{(k)}$ stands for the computed solution at the $k$th iteration.
The right hand side vector $\bf{b}$ was chosen such that the exact solution of
(2) was a vector of all ones.

In the following, we apply the complete version of FGMRES method with right
preconditioning technique in conjunction with APSS preconditioner
$\mathcal{M}_{\alpha}.$ For the APSS preconditioner, we need to solve two SPD
linear systems including $\alpha I+A+\frac{1}{\alpha}B^{T}B$ and
$\alpha^{2}I+CC^{T}$ as sub-tasks. To solve these linear systems, we employ
the conjugate gradient (CG) method without preconditioning and the iteration
was stopped when the residual 2-norm is reduced by a factor of $10^{3}$ or the
maximum number of inner iteration is reached to 200.

In the tables below, we use “CPU”, “IT” and “RES” to represent the elapsed CPU
time to converge in second, iteration counts and relative residual,
respectively. In addition, Degree of freedom (DOF) is defined as
$DOF=n+m+l={\bf n}$. Finally, a dagger (†), means that more than the maximum
number of iterations are needed to converge. All the examples were performed
in Matlab-R2019a. All the numerical results were obtained by a Laptop by the
following features:

intel (R), Core(TM) i5-8265U, CPU @ 1.60 GHz, 8 GB.

######  Example 1.

Consider the two dimensional leaky lid-driven cavity problem

$\begin{dcases}-\Delta\mathbf{u}+\nabla p=0,&\text{ in }\Omega,\\\
\nabla\cdot\mathbf{u}=0,&\text{ in }\Omega,\end{dcases}$ (35)

in which a suitable boundary conditions are applied on the side and bottom
points. $\mathbf{u}$ and $p$ symbolized the velocity vector field and the
pressure scalar field, respectively. In addition, $\Delta$ and $\nabla$ refer
to the vector Laplacian in $\mathbb{R}^{2}$ and the gradient, respectively.
This problem is called the Stokes problem. For the saddle point problem (2),
the matrices $A$ and $B$ come from the Stokes problem. To obtain these
matrices, we use the IFISS package by Elman et al. [15]. Note that the
stablized $Q_{1}-P_{0}$ and $Q_{2}-P_{1}$ finite element method (FEM) is
employed to discretize the Stokes equation (35). The grid parameters are
chosen as $h=\frac{1}{8},\frac{1}{16},\ldots,\frac{1}{256},$ for all uniform
or stretched grid points. Now, the matrix $C$ is taken to be the form

$C=\left(\begin{array}[]{l}C_{1}\\\
C_{2}\end{array}\right)=\left(\begin{array}[]{l}C_{1}\\\ c_{2}\\\
c_{2}\end{array}\right)\in\mathbb{R}^{(l+2)\times m},$

where
$C_{1}=\left(\begin{array}[]{ll}\operatorname{diag}\\{1,3,5,\ldots,2l-1\\}&\operatorname{randn}(l,m-l)\end{array}\right)\in\mathbb{R}^{l\times
m},$ and $c_{1}=\left(\begin{array}[]{ll}e^{T};&0^{T}\end{array}\right)C_{1},$
$c_{2}=\left(\begin{array}[]{ll}0^{T};&e^{T}\end{array}\right)C_{1},$
$e^{T}=\left(\begin{array}[]{llll}1;&1;&1;\cdots,&1\end{array}\right)\in\mathbb{R}^{\frac{l}{2}}.$
Note that $l=m-2.$ Here, $\operatorname{randn}(l,m-l)$ is a normally
distributed random matrix of order $l\times m-l.$

This example is a technical variant of Example 1 in [29]. The results are
given in Tables 4-4. From these tables, we can observe that the FGMRES method
in conjunction with the APSS preconditioner $\mathcal{M}_{\alpha}$ is strongly
more efficient than the FGMRES without a preconditioner. As we can see from
Figure 1, the preconditioner $\mathcal{M}_{\alpha}$ is efficient to cluster
the eigenvalues of the original coefficient matrix. Another observation which
can be posed here is that the number of iterations of the FGMRES method with
the APSS preconditioner remains almost constant when the size of the problem
is increased, whereas this is not the case for FGMRES without preconditioning.

Figure 1:  Eigenvalues distribution of $\mathcal{A}$ and
$\mathcal{M}_{\alpha}^{-1}\mathcal{A}$ (from the left to right) for
$h=\frac{1}{16}$ for Example 1 .  Table 1: Numerical results for FGMRES to
solve Example 1 with $Q_{1}-P_{0}$ FEM on uniform grids.

Prec. |  $h(DOF)$ |  $\frac{1}{8}$ (286) |  $\frac{1}{16}(1086)$ |  $\frac{1}{32}(4222)$ |  $\frac{1}{64}(16638)$ |  $\frac{1}{128}(66046)$ |  $\frac{1}{256}(263166)$
---|---|---|---|---|---|---|---
$\mathcal{I}$ |  IT |  49 |  103 |  552 |  188 |  474 |  1492
CPU |  0.0416 |  0.0600 |  0.4866 |  0.5810 |  4.1456 |  135.3632
RES |  2.2e-08 |  9.9e-08 |  1.0e-07 |  9.3e-08 |  9.8e-08 |  1.1e-07
$\mathcal{M}_{\alpha}$ |  $\alpha_{est}$ |  0.0396 |  0.0201 |  0.0101 |  0.0051 |  0.0025 |  0.0013
IT |  11 |  13 |  11 |  9 |  10 |  11
CPU |  0.0795 |  0.1009 |  0.0915 |  0.1807 |  0.6871 |  6.6438
RES |  9.9e-08 |  4.6e-08 |  8.9e-08 |  3.1e-08 |  2.3e-08 |  6.6e-08

Table 2: Numerical results for FGMRES to solve Example 1 with $Q_{1}-P_{0}$
FEM on stretched grids.

Prec. |  $h(DOF)$ |  $\frac{1}{8}$ (286) |  $\frac{1}{16}(1086)$ |  $\frac{1}{32}(4222)$ |  $\frac{1}{64}(16638)$ |  $\frac{1}{128}(66046)$ |  $\frac{1}{256}(263166)$
---|---|---|---|---|---|---|---
$\mathcal{I}$ |  IT |  48 |  133 |  93 |  187 |  406 |  974
CPU |  0.0391 |  0.0625 |  0.1089 |  0.5716 |  3.5570 |  118.1586
RES |  3.2e-08 |  7.5e-08 |  9.9e-08 |  9.9e-08 |  9.9e-08 |  1.0e-07
$\mathcal{M}_{\alpha}$ |  $\alpha_{est}$ |  0.3962 |  0.2001 |  0.1011 |  0.0051 |  0.0025 |  0.0013
IT |  11 |  13 |  7 |  7 |  7 |  9
CPU |  0.0739 |  0.0846 |  0.0807 |  0.1851 |  0.5726 |  8.6066
RES |  6.4e-08 |  3.5e-08 |  3.1e-08 |  6.5e-08 |  9.0e-08 |  3.1e-08

Table 3: Numerical results for FGMRES to solve Example 1 with $Q_{2}-P_{1}$
FEM on uniform grids.

Prec. |  $h(DOF)$ |  $\frac{1}{8}(254)$ |  $\frac{1}{16}(958)$ |  $\frac{1}{32}(3710)$ |  $\frac{1}{64}(14590)$ |  $\frac{1}{128}(57854)$ |  $\frac{1}{256}(230398)$
---|---|---|---|---|---|---|---
$\mathcal{I}$ |  IT |  51 |  130 |  186 |  253 |  947 |  $\dagger$
CPU |  0.0391 |  0.0653 |  0.1494 |  0.6904 |  7.8228 |
RES |  8.8e-08 |  9.0e-08 |  9.8e-08 |  9.9e-08 |  1.0e-07 |
$\mathcal{M}_{\alpha}$ |  $\alpha_{est}$ |  0.0419 |  0.0214 |  0.0108 |  0.0054 |  0.0027 |  0.0013
IT |  11 |  14 |  13 |  10 |  12 |  12
CPU |  0.0737 |  0.0771 |  0.0978 |  0.2462 |  1.0375 |  9.0643
RES |  3.5e-08 |  2.6e-08 |  7.7e-08 |  7.3e-08 |  5.8e-08 |  9.0e-08

Table 4: Numerical results for FGMRES to solve Example 1 with $Q_{2}-P_{1}$
FEM on stretched grids.

Prec. |  $h(DOF)$ |  $\frac{1}{8}(254)$ |  $\frac{1}{16}(958)$ |  $\frac{1}{32}(3710)$ |  $\frac{1}{64}(14590)$ |  $\frac{1}{128}(57854)$ |  $\frac{1}{256}(230398)$
---|---|---|---|---|---|---|---
$\mathcal{I}$ |  IT |  55 |  131 |  121 |  198 |  512 |  1137
CPU |  0.0380 |  0.0782 |  0.1123 |  0.5466 |  4.1606 |  129.2754
RES |  6.5-08 |  1.0e-07 |  1.0e-07 |  9.5e-08 |  9.9e-08 |  1.0e-07
$\mathcal{M}_{\alpha}$ |  $\alpha_{est}$ |  0.4195 |  0.0214 |  0.0108 |  0.0054 |  0.027 |  0.014
IT |  11 |  13 |  8 |  9 |  9 |  10
CPU |  0.0732 |  0.0722 |  0.0833 |  0.2363 |  0.9544 |  10.7704
RES |  9.6e-08 |  9.5e-08 |  9.8e-08 |  2.1e-08 |  4.2e-08 |  8.0e-08

######  Example 2.

Consider a technical modification of Example 1 in [20, 30] as follows:

$\displaystyle A=\left(\begin{array}[]{cc}{I\otimes T+T\otimes I}&{0}\\\
{0}&{I\otimes T+T\otimes I}\end{array}\right)\in\mathbb{R}^{2p^{2}\times
2p^{2}},$ $B=(I\otimes F\quad F\otimes I)\in\mathbb{R}^{p^{2}\times 2p^{2}}$

and

$C=\left(\begin{array}[]{l}C_{1}\\\
C_{2}\end{array}\right)=\left(\begin{array}[]{l}C_{1}\\\ c_{2}\\\
c_{2}\end{array}\right)\in\mathbb{R}^{(p^{2}+2)\times p^{2}},$

where

$C_{1}=E\otimes F\in\mathbb{R}^{p^{2}\times p^{2}},$

and

$c_{1}=\left(\begin{array}[]{ll}e^{T};&0^{T}\end{array}\right)C_{1},\qquad
c_{2}=\left(\begin{array}[]{ll}0^{T};&e^{T}\end{array}\right)C_{1},\qquad
e^{T}=\left(\begin{array}[]{llll}1;&1;&1;\cdots,&1\end{array}\right)\in\mathbb{R}^{\frac{p^{2}}{2}}.$

where

$\displaystyle
T=\frac{1}{h^{2}}\cdot\operatorname{tridiag}(-1,2,-1)\in\mathbb{R}^{p\times
p},\quad
F=\frac{1}{h}\cdot\operatorname{tridiag}(0,1,-1)\in\mathbb{R}^{p\times p},$

and $E=\operatorname{diag}\left(1,p+1,2p+1,\ldots,p^{2}-p+1\right)$ in which
the Kronecker product is denoted by $\otimes$, while the discretization mesh
size is represented by $h={1}/{(p+1)}$.

Table 5 reports the result of the FGMRES method and the FGMRES in conjunction
with the proposed preconditioner $\mathcal{M}_{\alpha},$ with respect to IT,
CPU and RES. Table 5 shows that the suggested preconditioner requires
significantly less iteration numbers and CPU time than the FGMRES without a
preconditioner. We also see that when the size of the problem is increased,
the number of iterations of the FGMRES without preconditioning increase
drastically, whereas for the FGMRES with the preconditioner
$\mathcal{M}_{\alpha}$ increases moderately. Figure 2 illustrates that the
preconditioned matrix works well in clustering the eigenvalues.

Figure 2: Eigenvalue distributions of $\mathcal{A}$ and $\mathcal{M}_{\alpha}^{-1}\mathcal{A}$ (from the left to right) for $h=\frac{1}{16}$ for Example 2 . Table 5: Numerical results for FGMRES to solve Example 2 . Prec. |  $p(DOF)$ |  8 (258) |  16 (1026) |  32 (4098) |  64 (16386) |  128 (65538)
---|---|---|---|---|---|---
$\mathcal{I}$ |  IT |  659 |  1999 |  $\dagger$ |  $\dagger$ |  $\dagger$
CPU |  0.1371 |  0.5270 |  |  |
RES |  9.9e-08 |  1.0e-07 |  |  |
$\mathcal{M}_{\alpha}$ |  $\alpha_{est}$ |  0.0434 |  0.0219 |  0.0110 |  0.0055 |  0.0027
IT |  13 |  14 |  15 |  17 |  27
CPU |  0.0700 |  0.0857 |  0.1448 |  0.5156 |  2.5788
RES |  2.4e-08 |  3.6e-08 |  4.5e-08 |  8.3e-08 |  7.3e-08

##  6 Conclusions

In this paper, the APSS method was employed to solve a class of nonsymmetric
three-by-three singular saddle point problems. We have applied the induced
preconditioner, $\mathcal{M}_{\alpha},$ to improve the semi-convergence rate
,when it is conjucted with FGMRES. We have proved that if $C$ is rank-
deficient the APSS method is unconditionally semi-convergent. Numerical tests
prove our theoretical claims.

##  7 Declaration

There is no conflict of interest.

##  8 Acknowledgements

The authors would like to that Prof. Fatemeh Panjeh Ali Beik (Vali-e-Asr
University of Rafsanjan) and Dr. Mohsen Masoudi (University of Guilan) for
their useful comments on an earlier version of this paper.

## References

  * [1] M. Abdolmaleki, S. Karimi, D.K. Salkuyeh, A new block-diagonal preconditioner for a class of 3$\times$ 3 block saddle point problems, Mediterr. J. Math. 19 (2022) 1-15.
  * [2] K. Arrow, L. Hurwicz, H. Uzawa, Studies in nonlinear programming, Stanford University Press, Stanford, 1958.
  * [3] H. Aslani, D.K. Salkuyeh, F.P.A. Beik, On the preconditioning of three-by-three block saddle point problems, Filomat 15 (2021) 5181-5194.
  * [4] Z.Z. Bai, On semi-convergence of Hermitian and skew-Hermitian splitting methods for singular linear systems, Computing, 89 (2010) 171-197.
  * [5] Z.Z. Bai, G.H. Golub, M.K. Ng, Hermitian and skew-Hermitian splitting methods for non-Hermitian positive definite linear systems, SIAM J. Matrix Anal. Appl. 24 (2003) 603-626.
  * [6] Z.-Z. Bai, B.N. Parlett, Z.-Q. Wang, On generalized successive overrelaxation methods for augmented linear systems, Numer. Math. 102 (2005) 1-38.
  * [7] A. Berman, R. Plemmons, Nonnegative Matrices in Mathematical Science. Academic Press, NewYork, 1979.
  * [8] S.L. Campbell, C.D. Meyer, Generalized Inverses of Linear Transformations, Pitman, London, 1979.
  * [9] Y. Cao, A block positive-semidefinite splitting preconditioner for generalized saddle point linear systems, J. Comput. Appl. Math. 374 (2020) 112787.
  * [10] Y. Cao, Shift-splitting preconditioners for a class of block three-by-three saddle point problems, Appl. Math. Lett. 96 (2019) 40-46.
  * [11] Z. Chao, G. Chen, Semi-convergence analysis of the Uzawa-SOR methods for singular saddle point problems, Appl. Math. Lett. 35 (2014) 52–57.
  * [12] Z. Chao, N.-M. Zhang, A generalized preconditioned HSS method for singular saddle point problems, Numer. Algor. 66 (2014) 203-221.
  * [13] Z. Chen, Q. Du, J. Zou, Finite element methods with matching and nonmatching meshes for Maxwell equations with discontinuous coefficients, SIAM J. Numer. Anal. 37 (2000) 1542-1570.
  * [14] P. Degond, P.A. Raviart, An analysis of the Darwin model of approximation to Maxwell’s equations, Forum Math. 4 (1992) 13-44.
  * [15] H.C. Elman, A. Ramage, D.J. Silvester, Algorithm 866: IFISS, a Matlab toolbox for modelling incompressible flow, ACM Trans. Math. Softw. 33 (2007) Article 14.
  * [16] D. Han, X. Yuan, Local linear convergence of the alternating direction method of multipliers for quadratic programs, SIAM J. Numer. Anal. 51 (2013) 3446-3457.
  * [17] K.B. Hu, J.C. Xu, Structure-preserving finite element methods for stationary MHD models, arXiv:1503.06160.
  * [18] N. Huang, Variable parameter Uzawa method for solving a class of block three-bythree saddle point problems, Numer. Algor. 85 (2020), 1233-1254.
  * [19] N. Huang, Y.-H. Dai, Q. Hu, Uzawa methods for a class of block three-by-three saddle-point problems, Numer Linear Algebra Appl. 26 (2019), e2265.
  * [20] N. Huang, C.-F. Ma, Spectral analysis of the preconditioned system for the $3\times 3$ block saddle point problem, Numer. Algor. 81 (2019) 421-444.
  * [21] Z.-Z. Liang, G.-F. Zhang, Alternating positive semidefinite splitting preconditioners for double saddle point problems, Calcolo 56 (2019) 26.
  * [22] W. Li, Y.-P. Liu, X.-F. Peng, The generalized HSS method for solving singular linear systems. J. Comput. Appl. Math. 236 (2012) 2338-2353.
  * [23] J.-L. Li, Q.-N. Zhang, S.-L. Wu, Semi-convergence of the local Hermitian and SkewHermitian splitting iteration methods for singular generalized saddle point problems, Appl. Math. E-Notes 11 (2011) 82-90.
  * [24] H.-F. Ma, N.-M. Zhang, A note on block-diagonally preconditioned PIU methods for singular saddle point problems, Int. J. Comput. Math. 88 (2011) 3448-3457.
  * [25] G.I. Marchuk, Methods of Numerical Mathematics, Springer-Verlag, New York, 1984.
  * [26] P. Monk, Analysis of a finite element method for Maxwell’s equations, SIAM J. Numer. Anal. 29 (1992) 714-729.
  * [27] D.K. Salkuyeh, H. Aslani, Z.Z. Liang, An alternating positive semidefinite splitting preconditioner for the three-by-three block saddle point problems, Math. Commun. 26 (2021) 177-195.
  * [28] L. Wang and K. Zhang, Generalized shift-splitting preconditioner for saddle point problems with block three-by-three structure, Open Access Library Journal 6 (2019) e5968.
  * [29] N.N. Wang, J.C. Li, On parameterized block symmetric positive definite preconditioners for a class of block three-by-three saddle point problems, J. Comput. Appl. Math. 405 (2022) 113-959.
  * [30] X. Xie, H.-B. Li, A note on preconditioning for the 3 $\times$ 3 block saddle point problem, Comput. Math. Appl. 79 (2020) 3289-3296.
  * [31] J.S. Xiong, X.B. Gao, Semi-convergence analysis of Uzawa-AOR method for singular saddle point problems, Comput. Appl. Math. 36 (2017) 383-395.
  * [32] J.-Y. Yuan, Numerical methods for generalized least squares problems, J. Comput. Appl. Math. 66 (1996) 571-584.
  * [33] N.-M. Zhang, T.-T. Lu, Y.-M. Wei, Semi-convergence analysis of Uzawa methods for singular saddle point problems, J. Comput. Appl. Math. 255 (2014) 334-345.
  * [34] G.-F. Zhang, S.-S. Wang, A generalization of parameterized inexact Uzawa method for singular saddle point problems, Appl. Math. Comput. 219 (2013) 4225-4231.
  * [35] B. Zheng, Z.-Z. Bai, X. Yang, On semi-convergence of parameterized Uzawa methods for singular saddle point problems, Linear Algebra Appl. 431 (2009) 808-817.