# Deep Non-Monotonic Reasoning for Visual Abstract Reasoning Tasks

Yuan Yang<EMAIL_ADDRESS>
Deepayan Sanyal<EMAIL_ADDRESS>
Joel Michelson<EMAIL_ADDRESS>
James Ainooson<EMAIL_ADDRESS>
Maithilee Kunda<EMAIL_ADDRESS>
Vanderbilt University Department of Computer Science, 400 24th Ave S,
Nashville, TN 37212 USA

###### Abstract

While achieving unmatched performance on many well-defined tasks, deep
learning models have also been used to solve visual abstract reasoning tasks,
which are relatively less well-defined, and have been widely used to measure
human intelligence. However, current deep models struggle to match human
abilities to solve such tasks with minimum data but maximum generalization.
One limitation is that current deep learning models work in a monotonic way,
i.e., treating different parts of the input in essentially fixed orderings,
whereas people repeatedly observe and reason about the different parts of the
visual stimuli until the reasoning process converges to a consistent
conclusion, i.e., non-monotonic reasoning. This paper proposes a non-monotonic
computational approach to solve visual abstract reasoning tasks. In
particular, we implemented a deep learning model using this approach and
tested it on the RAVEN dataset—a dataset inspired by the Raven’s Progressive
Matrices test. Results show that the proposed approach is more effective than
existing monotonic deep learning models, under strict experimental settings
that represent a difficult variant of the RAVEN dataset problem.

Keywords: visual reasoning; abstract reasoning; Raven’s Progressive Matrices;
deep learning.

## Introduction

Visual abstract reasoning (VAR) tasks have long played a prominent role in
understanding and measuring human intellignece. For instance, VAR task
variants have been incorporated into many standard human intelligence tests,
for example: Wechsler, Stanford-Binet, Woodcock–Johnson, Kaufmann, etc.
Psychometrics research has showed that VAR tasks often significantly load on
core factors of intelligence, such as general intelligence and fluid
intelligence within different theoretical models of human intelligence, such
as Spearman’s $g$, the three-stratum model of Cattell, Horn, & Carroll, and
the $g$-VPR model of Johnson & Bouchard.

A canonical example of VAR tasks is Raven’s Progressive Matrices (RPM), which
has been widely used in intelligence testing since the 1930s, and was shown to
be the best single-format measure of general intelligence across a wide range
of psychometric tests (?, ?). We thus use RPM-style matrix reasoning problems
as being a representative and important variant of VAR tasks in the following
discussion.

A matrix reasoning item consists of a matrix of visual figures with the last
entry missing, along with several answer choices, with one answer choice that
correctly completes the patterns in the matrix. While the published RPM tests
contain on the order of several dozen items manually designed by test author
John Raven, recent efforts in AI has produced several large datasets of RPM-
like problems to support machine learning research for solving VAR tasks.
Figure 1 shows an item from one such dataset called RAVEN that was inspired by
the RPM (?, ?).

Figure 1: A sample item from the RAVEN dataset (?, ?). The correct answer is
the first image in the first row. It satisfies underlying four patterns in the
matrix: (1) three spatial arrangements repeat in each row, but in different
orders; (2) shapes in each column are the same; (3) the (relative) sizes of
shapes in each row are constant; (4) two dark entries and one white entry are
in each row.

Despite this extremely simple format, there is a lot we can tell about RPM.
First, how could one determine that an RPM item is correctly solved? In other
words, what are the criteria for an answer choice to be correct?
Unfortunately, this is usually unclear until one sees the specific item. The
designer of RPM items may use any patterns to create the items as long as the
patterns can be visually expressed. The only rule is that, as intelligence
tests, the underlying patterns should involve as less non-intelligence factors
as possible, for example, languages and education experiences. This does not
leave too many choices for item designers. In most cases, they can only use
common geometric shapes, natural things (e.g., sun, moon, trees, and rivers),
and body parts (hands and foot). Nonetheless, they created many novel and
diversified designs for RPM tests and other similar intelligence tests. And
more importantly, the correct answer is unambiguously favored by most human
subjects when they find by themselves or are told the underlying patterns.
Interesting! Isn’t it? No universal criteria for being correct or wrong, but
most people’s understanding converges to the same answer. This requires much
experience, especially when designing difficult items. For example, the item
in Figure 1 looks very much in disarray at first glance. But, if one observes
the underlying patterns of shape, color, and spatial arrangement (see the
caption of Figure 1), the first answer choice is obviously better than others.
On the contrary, if a human subject has difficulty in any cognitive process
related to extracting these underlying patterns, this item poses a great
intellectual challenge.

Second, how are these RPM items solved differently by different people? There
has been a long line of studies discussing various aspects of the human
problem solving process, such as high-level strategies, e.g., constructive
matching and response elimination (?, ?), information processing components
such as working memory capacity and selective attention (?, ?), and perceptual
complexity a(?, ?). These works provide rich inspiration for building AI
systems for solving RPM. Taking deep learning models as an example, these AI
systems can be divided into single-choice and multiple-choice models, roughly
corresponding to the constructive matching and response elimination strategies
used by human subjects; some of these systems mimic the structure mapping
process of analogy making to organize the information processing; and there
also exists models with built-in memory and attention mechanisms, albeit not
exactly the mechanisms of human cognition. However, these works, both in human
and AI, are mostly focusing and stressing the importance of certain aspects or
levels of the solving process. This has resulted in many ad hoc AI systems for
solving RPM with less focus on flexible and integrative strategies. In this
paper, we re-examine this question of more integrative solving processes for
RPM-like tasks from scratch (i.e., with system inputs in the form of images at
the pixel level).

Considering inputs at the pixel level poses an interesting challenge for AI
systems in terms of inferential complexity: If each entry in Figure 1 is a
$160\times 160$ gray scale image, there exist 409,600 pixels in the image of
the entire item. What pixels should an AI system process at each step of its
computation? In what order, are these pixels processed? How are they
processed? There are many possible answers to these questions. Eye tracking
studies(?, ?) show that human subjects switch between matrix entries, rows,
and columns during the entire solving process, but without a fixed order.
Generally speaking, it is hard to tell when they are inspecting a matrix
entry, when they are comparing multiple ones, and when they are drawing
conclusions of the whole matrix. In contrast, every existing AI system for
solving RPM assumes a certain order to process different parts of the visual
input. These orders vary from simple ones, like sequentially processing every
entry, to complex ones, like processing entries, pairs of entries, rows,
columns, double rows, double columns, and combining the processings in a even
more complex hierarchy. In particular, these AI systems are designed with an
assumption that the solving process—from pixel values to final predictions—is
a monotonic one, where we use the term monotonic to refer to a systematic
ordering in how the reasoner processes elements of the problem input (?, ?, ?,
?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?). While this assumption
is effective for many deep learning tasks, it is limiting and not
representative of human-like reasoning flexibility for solving VAR tasks. For
example, the process could be non-monotonic by moving back and forth between
perceptual and conceptual processings on different parts of raw visual input
until it settles on a balanced point. Such a non-monotonic approach seems more
reflective of human-like VAR reasoning, but such approaches have not received
sustained attention in AI research on such tasks.

Again, take the RPM item in Figure 1 as an example. If we insert each answer
choice into the matrix, we have 8 3$\times$3 matrices with each matrix entry a
160$\times$160 images. Since there are limited number of items in RPM-like
datasets, we can assume that we have $N$ underlying patterns and each
underlying pattern has $M$ instances (i.e., items) in this dataset. Then, this
RPM-like dataset can be represented by a 7D array of shape
$(N,M,a=8,r=3,c=3,h=160,w=160)$111$a$ for answer choices, $r$ for rows, $c$
for columns, $h$ for height of matrix entry images, and $w$ for width of
them.. If we want our AI systems to learn the underlying patterns in the
dataset, how does it process this huge array? If the $(h,w)$ and $(N)$
dimensions correspond to the lowest and highest cognitive levels,
respectively, the the monotonic AI system would move the processing from the
right to the left of $(N,M,c=8,r=3,c=3,h=160,w=160)$, without backtracking.
The monotonic assumption here means that correctly deciphering a dimension is
conditioned on correctly deciphering the ones on its right. Now, imagine a
situation where all the entry images, and thus all the items, in the RPM-like
dataset are the same. In this case, processing the first five dimensions
first, rather than the last two, would be more efficient and easier for the AI
system to detect the underlying pattern—everything is constant. This might be
a very extreme example, but the general idea is clear. There should not be any
fixed order to process input dimensions; at a certain point of processing for
solving certain items, some dimensions are more informative than others.
Therefore, the non-monotonic approach has its advantage over the monotonic
approach as it sees all dimensions and allows to select the informative ones
at each step.

Although the non-monotonic approach seems appealing, such freedom of
processing poses another problem—we do not know what dimensions should be
processed at each step and so many choices would make the model intractable.
In this paper, we provide one solution to this problem:

  1. 1.

We describe a general non-monotonic computational approach for solving VAR
tasks.

  2. 2.

We implement a deep learning model to verify the feasibility of our non-
monotonic approach.

  3. 3.

We evaluate our model on the RAVEN dataset, using the strictest experimental
settings that involve the least training data but the most underlying
patterns, without any meta information of item generation or prior knowledge
about the test, and using single-choice prediction, and we show that our model
outperforms other deep learning models on this most challenging training and
evaluation regime.

## Related Work

There is a long line of research of computational models solving VAR tasks,
especially RPM-like tasks. A technical taxonomy classifies these works into
four categories. We briefly describe them here. More extensive reviews can be
found in (?, ?) and (?, ?).

### Imagery-Based Approach

Visual mental imagery refers to mental images (?, ?). It plays a crucial role
in human’s ability to solve VAR tasks. The most important characteristic of
mental imagery is that human can experience mental imagery in the absence of
the concurrent sensory input. The imagery-based approach simulates human
mental imagery by directly operating on the raw visual input and, through
mental operations, such as rotation, addition, and subtraction, it can solve a
substantial portion of original RPM tests (?, ?).

### Logical Reasoning

The computational models using logical reasoning works on symbolic
representations of RPM-like items and reason in formal logic. For example, a
entry image $A$ in a matrix could be described by a set of propositions:
“triangle($A$)=True, triangle-large($A$)=False, triangle-on-the-
left($A$)=True, square($A$)=True, square-small($A$)=True, etc”. The symbolic
representations in these models are either manually constructed or obtained
through a preprocessing module. Representatives of this approach are ANALOGY
(?, ?), FAIRMAN ,and BETTERMAN (?, ?).

### Neuro-Symbolic Approach

The neuro-symbolic model usually consist of two modules—a neural perception
frontend and a symbolic reasoning backend. The neural perception frontend
(usually implemented as neural networks) extracts/approximates the
distributions over the values of objects in each entry in a predefined formal
representation system. The symbolic reasoning backend performs probability
calculation according to a predefined set of rules. Neuro-symbolic approach
and the next approach are data-driven approach, whereas the first two
approaches are knowledge-based approaches. Examples of neuro-symbolic models
for solving RPM-like tasks include ALANS2, PrAE, VAE-GPP, TRIVR, LoGe, and
NVSA (?, ?, ?, ?, ?, ?, ?).

### Learning Approach

Unlike the previous approach, learning approach does not rely on predefined
representation systems of geometric objects and abstract patterns. Instead,
the representations are learned from raw perceptual data. This approach has
become more popular because of large datasets of VAR tasks created in the last
five years (?, ?, ?, ?, ?). When this paper is written, almost of all the
popular deep learning architectures have been applied to solve VAR tasks, such
as CNN, ResNet family, recurrent neural networks, and attention models. The
proposed model in the following sections also falls into this category.

## Method

Figure 2: Non-monotonic processing of a multi-dimensional array $D$. This
figure visualizes the computation of Equation (1), (2), (3), and (4).

Let’s continue the example of the RPM item in the introduction section, which
is a 7D array of shape $(N,M,a=8,r=3,c=3,h=160,w=160)$. Each dimension
represent a distinct type of information about this huge input array (although
the information in different dimensions might be entangled). For example, the
$(r,c)$ dimensions are used to realize the underlying patterns; the difference
between answer choices, i.e., how well they complete the underlying patterns,
is realized through the $a$ dimension; and different types of underlying
patterns are organized in the $N$ dimension. All the first five dimensions are
related to some information that is more abstract than the last two dimensions
$(h,w)$. Thus, we refer to them as conceptual and perceptual dimensions,
respectively. To further simplify the discussion, we consider the 4D sub-array
of shape $(r=3,c=3,h=160,w=160)$, which represents a single matrix completed
by an (correct or wrong) answer choice, with $(h=160,w=160)$ and $(r=3,c=3)$
as perceptual and conceptual dimensions, respectively. The task is then
simplified to taking as input such a 4D array and providing as output
indicating whether or not the answer choice is the correct one. So, it becomes
a binary classification problem.

As mentioned above, monotonic systems process perceptual and conceptual
dimensions sequentially, whereas non-monotonic systems are allowed to process
in any order. Generally, if we want to process conceptual dimensions for $a$
times and perceptual dimensions for $b$ times, there are ${a+b\choose a}$ ways
to perform these processings sequentially. Since none of them is guaranteed to
be more effective than others, we can try to perform all of them in parallel.

But how is this different from using multiple monotonic systems simultaneously
in parallel? Well, let’s take a close look at our example of the 4D array. We
first insert a trivial dimension as an abstract feature dimension and it is
thus of shape $(r=3,c=3,h=160,w=160,f=1)$. Consider two possible processing
paths.

  * •

Path 1:
$(r=3,c=3,h=160,w=160,f=1)\xRightarrow{f_{h,w}}(r=3,c=3,h=1,w=1,f=16)\xRightarrow{f_{r,c}}(r=1,c=1,h=1,w=1,f=32)$

  * •

Path 2:
$(r=3,c=3,h=160,w=160,f=1)\xRightarrow{f^{\prime}_{r,c}}(r=1,c=1,h=160,w=160,f=16)\xRightarrow{f^{\prime}_{h,w}}(r=1,c=1,h=1,w=1,f=32)$

$f_{h,w}$ and $f^{\prime}_{h,w}$ are perceptual processings on the $(h,w)$
dimensions, with $(r,c)$ dimensions being transparent to them. Similarly,
$G_{r,c}$ and $G^{\prime}_{r,c}$ are conceptual processings on the $(r,c)$
dimensions, with $(h,w)$ dimensions being transparent to them. In this case,
Path 1 represents the model’s bottom-up “understanding” of underlying patterns
of the matrix, whereas Path 2 represents the model’s top-down “understanding”
of underlying pattern of the matrix. Which one is correct? There is no
universal answer to this question. The only condition one can tell is that
only when the outputs of the two paths are consistent (or even identical) with
each other can they be correct222This is very common in optimization and
machine learning models. For example, many optimization problem is optimized
when it meets its dual problem, i.e., they are optimized simultaneously., for
that two different understandings of the same thing are both correct is
impossible 333This is true when the RPM item is well-designed for intelligence
testing. The ambiguous case of VAR tasks is beyond the scope of this paper..
Therefore, to make non-monotonic approach work and work better than an
ensemble of multiple monotonic processing paths in parallel, we can try to
force the consistency between the multiple paths, through a contrasting
module.

Now, we formalize the non-monotonic approach. Suppose that we have an input
array $D$ of shape $(d_{1},d_{2},...,d_{l})$. Let
$\Omega_{1},\Omega_{2},...,\Omega_{T}\subset\\{1,2,...,l\\}$ be subsets of
dimensions that we are going to process sequentially. Thus, there exist $K=T!$
possible processing paths and we create $K$ copies
$\\{D^{(0)}_{1},D^{(0)}_{2},...,D^{(0)}_{K}\\}$ of $D$, with each copy for a
path. A non-monotonic step, say Step $j$, transforms
$\\{D^{(j)}_{1},D^{(j)}_{2},...,D^{(j)}_{K}\\}$ to
$\\{D^{(j+1)}_{1},D^{(j+1)}_{2},...,D^{(j+1)}_{K}\\}$ through Equation (1),
(2), (3), and (4). For all path $i\in\\{1,2,...,K\\}$:

$\displaystyle f^{(j)}_{i}=f^{(j)}_{\Omega_{i_{T}}}\circ
f^{(j)}_{\Omega_{i_{T-1}}}\circ\cdots\circ f^{(j)}_{\Omega_{i_{1}}}$ (1)
$\displaystyle P^{(j)}_{i}=f^{(j)}_{i}\left(D^{(j)}_{i}\right)$ (2)
$\displaystyle
C^{(j)}=g^{(j)}\left(P^{(j)}_{1},P^{(j)}_{2},...,P^{(j)}_{K}\right)$ (3)
$\displaystyle D^{(j+1)}_{i}=P^{(j)}_{i}-C^{(j)}$ (4)

where $f^{(j)}=(f^{(j)}_{1},f^{(j)}_{2},...,f^{(j)}_{K})$ and $g^{(j)}$ are
processing paths and contrasting module of Step $j$. The processing paths
$f^{(j)}$ and the contrasting module $g^{(j)}$ can be implemented using
standard deep learning modules, such as convolutions, linear, and pooling, or
other custom operations, as long as the input and output dimensions are
compatible with neighboring modules.

Figure 3: 7 spatial configurations of the RAVEN dataset.

Figure 2 depicts the non-monotonic approach in a more extended way. In
practice, we can stack multiple non-monotonic steps together, as we stack
multiple convolutional layers in monotonic models. At Step $0$, $f^{(0)}$
extracts information from certain dimensions of the input data array $D$ and
optionally reduces the sizes of the dimensions; $g^{(0)}$ then contrasts the
outputs from different paths. Since the parallel paths of $f^{(0)}$ are
processing the same dimensions in different orders, the outputs of different
paths are expected to be consistent (or even identical) with each other, as we
have discussed in the introduction section. Therefore, the contrasting module
$g^{(0)}$ computes how the outputs from different paths are not consistent.
Intuitively, the contrasting process will lead to two outcomes. First, if the
outputs are consistent, the contrasting module will produce a certain type of
output to give the following modules a signal that “I have had a good
understanding of the input”, and the following modules then respond to this
signal by producing the desired final output, e.g., $Y=0$ in the Figure 2.
Second, if the outputs from different paths are not consistent, the
contrasting module will filter out the components that causes the
inconsistency, for which the following modules are more and more likely to be
able to “figure out” why they are inconsistent, and, if enough non-monotonic
steps are stacked, achieve the status of being consistent. To sum up, in
either case, the model aims to achieve the consistency between different ways
to decipher the raw input.

## Experiment

### Dataset

To verify feasibility of the non-monotonic approach, we implemented a deep
learning model to solve the RAVEN dataset (?, ?)—an RPM-like dataset that was
procedurally generated according to 4 predefined underlying patterns
(constant, progression, arithmetic, distribute of three) and 5 predefined
geometric elements (number/position, type, size, color)444The values of
geometric elements are limited in finite discrete sets. Thus, the RAVEN
dataset is not as diversified as real-world intelligence tests.. All the items
in RAVEN share the same format in Figure 1, i.e., a $3\times 3$ matrix with 8
answer choices. The goemetric elements in RAVEN are organized in 7 predefined
spatial configurations, as in Figure 3. Each configuration was used to
generated 10,000 items. Thus, there are 70,000 items in total. We use the
common splits of 60%, 20%, and 20% for training, validation, and testing,
respectively.

### Model Implementation

Our model takes as input arrays of shape $(r=3,c=3,h=160,w=160)$, which
represents a matrix completed by each of its answer choices, and provides as
output the probability that the answer choice is correct. For each item, the
answer choice with the highest probability is selected as the model’s response
to the item, which is used to compute the accuracy.

Table 1: The implementation of our model. Note that the parameters“k”, “s”, and “p” are kernel size, stride, and padding, respectively. The “$\times$2” imples that there are two data arrays of the shape, each being processed by a path, as shown in Figure 2. Layer Name | Output Shape | Implementation
---|---|---
Input | (3$\times$3,1,160,160) |
Conv1 | (3$\times$3,64,40,40) | Conv2d(1,64,k=7,s=2,p=3), BatchNorm2d(64), Relu, Maxpool2d(k=3,s=2,p=1)
Conv2 | (3$\times$3,64,20,20) | Conv2d(64,64,k=3,s=1,p=1), BatchNorm2d(64), Relu, Maxpool2d(k=3,s=2,p=1)
$f^{(0)}$ | (3,3,128,10,10)$\times$2 |  | Path 1: ResidualBlock2d(64,64,k=3,s=1), ResidualBlock2d(64,128,k=3,s=2)
---
Path 2: ResidualBlock2d(64,128,k=3,s=2), ResidualBlock2d(128,128,k=3,s=1)
$g^{(0)}$ | (3,3,128,10,10)$\times$2 | Avg of its two inputs, Conv3d(128,k=3,s=1,p=1),BatchNorm3d(128), Subtraction
$f^{(1)}$ | (3,3,256,5,5)$\times$2 |  | Path 1: ResidualBlock2d(128,128,k=3,s=1), ResidualBlock2d(128,256,k=3,s=2)
---
Path 2: ResidualBlock2d(128,256,k=3,s=2), ResidualBlock2d(256,256,k=3,s=1)
$g^{(1)}$ | (3,3,256,5,5)$\times$2 | Avg of its two inputs, Conv3d(256,k=3,s=1,p=1),BatchNorm3d(256), Subtraction
Pool | (256)$\times$2 | AvgPool4d(1)
MLP | (1) | Multi-Layer Perceptron of sizes [512, 256, 1]

Table 1 describes the specific implementation of our model in the experiment.
We use two convolutional layers before non-monotonic steps to reduce the sizes
of $(h,w)$ dimensions, because they are too large and cost too much
computation if we use only non-monotonic steps. In this model, we use two non-
monotonic steps to process the $(r,c)$ and $(h,w)$ dimensions
alternatively555The contrasting modules should have been 4D convolutions. But,
it is too inconvenient to implement it in Pytorch and we used 3D convolutions
instead.. The model is implemented in Pytorch and trained with ADAM optimizer.

Table 2: Model Model Performance on RAVEN (?, ?, ?). CoPINet is evaluated on a variant of the RAVEN dataset since its result on the original RAVEN is not reliable(?, ?). Model | Avg. Acc. | Center | 2x2Grid | 3x3Grid | L-R | U-D | O-IC | O-IG
---|---|---|---|---|---|---|---|---
CNN | 36.97% | 33.58% | 30.30% | 33.53% | 39.43% | 41.26% | 43.20% | 37.54%
ResNet | 53.43% | 52.82% | 41.86% | 44.29% | 58.77% | 60.16% | 63.19% | 53.12%
LSTM | 13.07% | 13.19% | 14.13% | 13.69% | 12.84% | 12.35% | 12.15% | 12.99%
WReN | 14.69% | 13.09% | 28.62% | 28.27% | 7.49% | 6.34% | 8.38% | 10.56%
CoPINet | 46.1% | 54.4% | 36.8% | 31.9% | 52.2% | 42.8% | 51.9% | 52.5%
Ours | 79.32% | 95.6% | 55.25% | 56.95% | 97.55% | 96.15% | 95.95% | 58.5%

### RAVEN Evaluation Settings

Our experiment is by no mean comprehensive since there are several different
variants of RAVEN(?, ?) and several different evaluation protocols (e.g.,
access to meta information of item generation, prior knowledge about the test,
single-choice or multi-choice, and generation or decision). The combination of
all these exponentially increases the cost of experiments, which we cannot
afford in this single study. And the task difficulty also varies greatly in
different combinations. Since our goal is to solve general VAR tasks, we chose
the most ordinary one without using the distinctive features of RAVEN, i.e.,
the experimental setting of the least data but the most underlying patterns,
without any meta information of item generation or prior knowledge about the
test, and using single-choice prediction, in which the model needs to assess
each answer choice separately, i.e., not allowed to compare them. Therefore,
our experiment setting is the hardest one among all the settings for solving
RAVEN, which has never been used before to the best of our knowledge.

### Results

Table 2 gives the accuracies of our model and other important baseline models
on RAVEN. Our model achieves the best accuracies on both the entire dataset
and each subset of spatial configuration. The results of WReN and CoPINet were
obtained in similar settings to ours. The other three models—CNN, ResNet, and
LSTM—are common deep learning models to provide a clue of how RAVEN compares
to other deep learning tasks. Admittedly, there exist works, such as (?, ?, ?,
?, ?), that can achieve higher performance than ours by using easier
experimental settings, such as using meta information of item generation, such
as using the underlying rules and geometric elements as additional supervising
signal, and prior knowledge about RAVEN, such as decomposing a matrix into
rows and processing each row individually (because the underlying patterns in
RAVEN are represented only by rows). Because these works use significantly
more input information and prior knowledge, we do not include them for direct
comparison here.

Table 2 shows that our model does not work as well on the 2x2Grid, 3x3Grid,
and Out-InGrid configurations, compared to other configurations. All three of
these configurations contain grid structures, which allow more variation in
number and position of geometric elements. These three configurations and
other four configurations, respectively, represent two common types of
perceptual organization in VAR tasks—classical view and normal view (?, ?)—in
which the correspondence between geometric elements are established in very
different ways. Our model is good at one of them but not the other. However,
it still achieves more than 50% on its worst configurations, which is strong
given the random baseline of 12.5%. This raises an interesting question for
further investigation—because the model is learning abstract patterns and
already achieved 50% accuracy, it must have learned some representations of
the abstract concepts and be able to detect them in some cases; but what stops
it from detecting them in other cases?

## Concluding Remarks

Through RPM, we can see some common characteristics of visual abstract
reasoning. First, there must be multiple parts in an VAR item because the
abstract concepts are in nature relations between these multiple parts. In
most cases, it requires additional dimensions to organize the multiple parts.
The more complex the abstract concepts, the more additional dimensions needed.
This nature of VAR task calls for effective multi-dimensional processing
methods. The proposed non-monotonic approach is such a method.

Second, the name “visual abstract reasoning” has been misleading for building
AI systems as VAR tasks in intelligence tests are in no way pure reasoning
problems. RPM tests and other VAR tests are for testing fluid intelligence,
which is largely a superset of reasoning ability, or a precondition for the
reasoning ability. This point can be seen through the administration
procedures of RPM tests, in which very sparse instructions are given. The
human subject only knows that she needs to point at one of the answer choices
to complete the matrix, but have no idea of what “complete” means and what is
the criterion of being correct or wrong. Therefore, VAR as intelligence tests
is primarily about discovering and secondarily about reasoning. In contrast,
the AI systems for solving VAR is mainly working on the reasoning part, with
the discovering part ideally formulated by human designers. The proposed
approach intends to capture the discovering part by pursuing the consistency
between different processing paths. This approach is by no means the only way
to do so, and it might not be a very efficient one. But the crux of building
capable AI is definitely on the “discovering” part.

## References

  * An ChoAn Cho An, J., Cho, S.  (2020).  Hierarchical transformer encoder with structured representation for abstract reasoning.  _IEEE Access_ , _8_ , 200229–200236.
  * Arendasy SommerArendasy Sommer Arendasy, M., Sommer, M.  (2005).  The effect of different types of perceptual manipulations on the dimensionality of automatically generated figural matrices.  _Intelligence_ , _33_(3), 307–324.
  * Barrett et al.Barrett et al. Barrett, D., Hill, F., Santoro, A., Morcos, A., Lillicrap, T.  (2018).  Measuring abstract reasoning in neural networks.  In _International conference on machine learning_ (pp. 511–520).
  * Benny et al.Benny et al. Benny, Y., Pekar, N., Wolf, L.  (2021).  Scale-localized abstract reasoning.  In _Proceedings of the ieee/cvf conference on computer vision and pattern recognition_ (pp. 12557–12565).
  * Bethell-Fox et al.Bethell-Fox et al. Bethell-Fox, C. E., Lohman, D. F., Snow, R. E.  (1984).  Adaptive reasoning: Componential and eye movement analysis of geometric analogy performance.  _Intelligence_ , _8_(3), 205–238.
  * Carpenter et al.Carpenter et al. Carpenter, P. A., Just, M. A., Shell, P.  (1990).  What one intelligence test measures: a theoretical account of the processing in the raven progressive matrices test.  _Psychological review_ , _97_(3), 404.
  * EvansEvans Evans, T. G.  (1964).  _A program for the solution of a class of geometric-analogy intelligence-test questions._ (Tech. Rep.).  AIR FORCE CAMBRIDGE RESEARCH LABS LG HANSCOM FIELD MASS.
  * Hahne et al.Hahne et al. Hahne, L., Lüddecke, T., Wörgötter, F., Kappel, D.  (2019).  Attention on abstract visual reasoning.  _arXiv preprint arXiv:1911.05990_.
  * He et al.He et al. He, W., Ren, J., Bai, R.  (2021).  Two-stage rule-induction visual reasoning on rpms with an application to video prediction.  _arXiv preprint arXiv:2111.12301_.
  * Hersche et al.Hersche et al. Hersche, M., Zeqiri, M., Benini, L., Sebastian, A., Rahimi, A.  (2022).  A neuro-vector-symbolic architecture for solving raven’s progressive matrices.  _arXiv preprint arXiv:2203.04571_.
  * Hill et al.Hill et al. Hill, F., Santoro, A., Barrett, D. G., Morcos, A. S., Lillicrap, T.  (2019).  Learning to make analogies by contrasting abstract relational structure.  _arXiv preprint arXiv:1902.00120_.
  * Hoshen WermanHoshen Werman Hoshen, D., Werman, M.  (2017).  Iq of neural networks.  _arXiv preprint arXiv:1710.01692_.
  * Hu et al.Hu et al. Hu, S., Ma, Y., Liu, X., Wei, Y., Bai, S.  (2020).  Stratified rule-aware network for abstract visual reasoning.  _arXiv preprint arXiv:2002.06838_.
  * Jahrens MartinetzJahrens Martinetz Jahrens, M., Martinetz, T.  (2020).  Solving raven’s progressive matrices with multi-layer relation networks.  In _International joint conference on neural networks._
  * Kiat et al.Kiat et al. Kiat, N. Q. W., Wang, D., Jamnik, M.  (2020).  Pairwise relations discriminator for unsupervised raven’s progressive matrices.  _arXiv preprint arXiv:2011.01306_.
  * Kosslyn et al.Kosslyn et al. Kosslyn, S. M., Thompson, W. L., Ganis, G.  (2006).  _The case for mental imagery_ (1st. ed.).  New York, NY: Oxford University Press.
  * Małkiński MańdziukMałkiński Mańdziuk Małkiński, M., Mańdziuk, J.  (2022).  Deep learning methods for abstract visual reasoning: A survey on raven’s progressive matrices.  _arXiv preprint arXiv:2201.12382_.
  * Pekar et al.Pekar et al. Pekar, N., Benny, Y., Wolf, L.  (2020).  Generating correct answers for progressive matrices intelligence tests.  _Advances in Neural Information Processing Systems_ , _33_ , 7390–7400.
  * PrimiPrimi Primi, R.  (2001).  Complexity of geometric inductive reasoning tasks: Contribution to the understanding of fluid intelligence.  _Intelligence_ , _30_(1), 41–70.
  * Rahaman et al.Rahaman et al. Rahaman, N., Gondal, M. W., Joshi, S., Gehler, P., Bengio, Y., Locatello, F., Schölkopf, B.  (2021).  Dynamic inference with neural interpreters.  _Advances in Neural Information Processing Systems_ , _34_.
  * Shekhar TaylorShekhar Taylor Shekhar, S., Taylor, G. W.  (2021).  Neural structure mapping for learning abstract visual analogies.  In _Svrhm 2021 workshop@ neurips._
  * Shi et al.Shi et al. Shi, F., Li, B., Xue, X.  (2021).  Raven’s progressive matrices completion with latent gaussian process priors.  In _Proceedings of the aaai conference on artificial intelligence_ (Vol. 35, pp. 9612–9620).
  * Sinha et al.Sinha et al. Sinha, I., Webb, T. W., Cohen, J. D.  (2020).  A memory-augmented neural network model of abstract rule learning.  _arXiv preprint arXiv:2012.07172_.
  * Snow et al.Snow et al. Snow, R. E., Kyllonen, P. C., Marshalek, B., et al.  (1984).  The topography of ability and learning correlations.  _Advances in the psychology of human intelligence_ , _2_(S 47), 103.
  * Spratley et al.Spratley et al. Spratley, S., Ehinger, K., Miller, T.  (2020).  A closer look at generalisation in raven.  In _Computer vision–eccv 2020: 16th european conference, glasgow, uk, august 23–28, 2020, proceedings, part xxvii 16_ (pp. 601–616).
  * Steenbrugge et al.Steenbrugge et al. Steenbrugge, X., Leroux, S., Verbelen, T., Dhoedt, B.  (2018).  Improving generalization for abstract reasoning tasks using disentangled feature representations.  _arXiv preprint arXiv:1811.04784_.
  * van Steenkiste et al.van Steenkiste et al. van Steenkiste, S., Locatello, F., Schmidhuber, J., Bachem, O.  (2019).  Are disentangled representations helpful for abstract visual reasoning?  _arXiv preprint arXiv:1905.12506_.
  * Wang et al.Wang et al. Wang, D., Jamnik, M., Lio, P.  (2020).  Abstract diagrammatic reasoning with multiplex graph networks.  _arXiv preprint arXiv:2006.11197_.
  * Wu et al.Wu et al. Wu, Y., Dong, H., Grosse, R., Ba, J.  (2021).  The scattering compositional learner: Discovering objects, attributes, relationships in analogical reasoning.  _arXiv preprint arXiv:2007.04212_.
  * Yang et al.Yang et al. Yang, Y., McGreggor, K., Kunda, M.  (2020).  Not quite any way you slice it: How different analogical constructions affect raven’s matrices performance.  In _Proceedings of the eighth annual conference on advances in cognitive systems (acs)._
  * Yang et al.Yang et al. Yang, Y., Sanyal, D., Michelson, J., Ainooson, J., Kunda, M.  (2022a).  Automatic item generation of figural analogy problems: A review and outlook.  In _The ninth annual conference on advances in cognitive systems._
  * Yang et al.Yang et al. Yang, Y., Sanyal, D., Michelson, J., Ainooson, J., Kunda, M.  (2022b).  A conceptual chronicle of solving raven’s progressive matrices computationally.
  * Yu et al.Yu et al. Yu, S., Mo, S., Ahn, S., Shin, J.  (2021).  Abstract reasoning via logic-guided generation.  _arXiv preprint arXiv:2107.10493_.
  * Zhang, Gao, et al.Zhang, Gao, et al. Zhang, C., Gao, F., Jia, B., Zhu, Y., Zhu, S.-C.  (2019).  Raven: A dataset for relational and analogical visual reasoning.  In _Proceedings of the ieee/cvf conference on computer vision and pattern recognition_ (pp. 5317–5327).
  * Zhang, Jia, et al.Zhang, Jia, et al. Zhang, C., Jia, B., Gao, F., Zhu, Y., Lu, H., Zhu, S.-C.  (2019).  Learning perceptual inference by contrasting.  In (Vol. 32).
  * Zhang et al.Zhang et al. Zhang, C., Jia, B., Zhu, S.-C., Zhu, Y.  (2021).  Abstract spatial-temporal reasoning via probabilistic abduction and execution.  In _Proceedings of the ieee/cvf conference on computer vision and pattern recognition_ (pp. 9736–9746).
  * Zhang et al.Zhang et al. Zhang, C., Xie, S., Jia, B., Zhu, Y., Wu, Y. N., Zhu, S.-C.  (2020).  Learning algebraic representation for abstract spatial-temporal reasoning.
  * Zheng et al.Zheng et al. Zheng, K., Zha, Z.-J., Wei, W.  (2019).  Abstract reasoning with distracting features.  _Advances in Neural Information Processing Systems_ , _32_.


# Inexact Adaptive Cubic Regularization Algorithms on Riemannian Manifolds and
Application

Z. Y. Lia and X. M. Wangb aEmail<EMAIL_ADDRESS>author,
Email<EMAIL_ADDRESS>a,bCollege of Mathematics and Statistics, Guizhou
University, Guiyang 550025, P. R. China

###### Abstract

The adaptive cubic regularization algorithm employing the inexact gradient and
Hessian is proposed on general Riemannian manifolds, together with the
iteration complexity to get an approximate second-order optimality under
certain assumptions on accuracies about the inexact gradient and Hessian. The
algorithm extends the inexact adaptive cubic regularization algorithm under
true gradient in [Math. Program., 184(1-2): 35-70, 2020] to more general cases
even in Euclidean settings. As an application, the algorithm is applied to
solve the joint diagonalization problem on the Stiefel manifold. Numerical
experiments illustrate that the algorithm performs better than the inexact
trust-region algorithm in [Advances of the neural information processing
systems, 31, 2018].

###### keywords:

Riemannian manifolds; retraction; inexact adaptive cubic regularization
algorithm; iteration complexity; joint diagonalization

##  1 Introduction

We consider the large-scale separable unconstrained optimization problem on
general Riemannian manifolds:

$\min_{x\in\mathcal{M}}f(x)=\frac{1}{n}\sum_{i=1}^{n}f_{i}(x),$ (1)

where $\mathcal{M}$ is a Riemannian manifold, $n\gg 1$, and each
$f_{i}:\mathcal{M}\rightarrow\mathbb{R}$ is continuously differentiable
($i:=1,2,\dots n$). Such problems frequently appear in machine learning and
scientific computing, where each $f_{i}$ is a loss (or misfit) function
corresponding to $i-$th observation (or measurement); see, e.g., [1, 2, 3]. In
such “large-scale” settings, since the evaluations of the gradient or the
Hessian of $f$ can be computationally expensive, some inexact techniques are
used for approximating the first and the second derivatives, which
particularly includes the random sampling technique. As a result, many popular
first-order and second-order stochastic algorithms are proposed to solve
problem (1). The stochastic gradient descent (SGD) algorithm is one of the
first-order stochastic/inexact algorithms, which uses only the gradient
information. The idea of the SGD comes from [4] and one can find some variants
in [5, 6, 27, 10, 11] and the references therein (the Riemannian versions of
the SGD can be found in [10, 11]). The second-order stochastic/inexact
algorithms always use the inexact Hessian information, which includes the
popular sub-sampling Newton method, sub-sampling/inexact trust region
algorithm, and the sub-sampling/inexact cubic regularization algorithm; see
e.g, [28, 29, 30, 31, 32, 17] (the Riemannian version of the inexact trust
region algorithm is in [17]). Particularly, we note the inexact trust-region
algorithm and the inexact adaptive cubic regularization algorithms introduced
in [32] for solving problem (1). Both algorithms employ the true gradient and
the inexact Hessian (by sub-sampling) of $f$ and solve the corresponding sub-
problems approximately every iteration. These two algorithms are respectively
shown to own iteration complexity
$\mathcal{O}(\max\\{\varepsilon_{g}^{-2}\varepsilon_{H}^{-1},\varepsilon_{H}^{-3}\\})$
and $\mathcal{O}(\max\\{\varepsilon_{g}^{-2},\varepsilon_{H}^{-3}\\})$ of
achieving the $(\varepsilon_{g},\varepsilon_{H})$-optimality (see definition
of the $(\varepsilon_{g},\varepsilon_{H})$-optimality in Definition 2.2).

It is known that the Riemannian geometry framework has some advantages in many
applications, such as translating some nonconvex (constrained) optimization
problems in Euclidean space into convex (unconstrained) ones over a Riemannian
manifold; see, e.g., [12, 34, 35, 36]. As a result, some classical numerical
methods for solving optimization problems on the Euclidean space, such as
Newton’s method, BFGS algorithm, trust region method, gradient algorithm,
subgradient algorithm, etc., have been extended to the Riemannian manifold
setting; see, e.g., [16, 12, 13, 14, 15, 34, 33, 35, 36]. As we have noted
above, for solving large-scale problem (1), the stochastic gradient descent
algorithm and its variants have been extended from Euclidean spaces to
Riemannian manifolds in the literature; see, e.g., [10, 11]. Recently, the
Riemannian inexact trust region algorithm is studied in [17], which uses the
inexact gradient and the inexact Hessian every iteration, and in addition the
inexact solution of the trust region sub-problem. Noting that the algorithm
employs the inexact gradient instead of the true gradient, it extends the
corresponding one in [32] to more general cases even in Euclidean settings.
Similarly, for achieving the $(\varepsilon_{g},\varepsilon_{H})$-optimality,
the iteration complexity of the algorithm is also proven to be
$\mathcal{O}(\max\\{\varepsilon_{g}^{-2}\varepsilon_{H}^{-1},\varepsilon_{H}^{-3}\\})$.
At the same time, some numerical results are provided which illustrate the
algorithm performs significantly better than state-of-the-art deterministic
and stochastic algorithms in some applications.

Inspired by the prior works (in particular, the works in [32] and [17]), we
propose the Riemannian inexact adaptive cubic regularization algorithm for
solving problem (1) on general Riemannian manifolds, which employs the inexact
gradient and the inexact Hessian every iteration. The algorithm is proven to
own the similar iteration complexity
$\mathcal{O}(\max\\{\varepsilon_{g}^{-2},\varepsilon_{H}^{-3}\\})$ for
obtaining the $(\varepsilon_{g},\varepsilon_{H})$-optimality under certain
conditions on accuracies about the inexact gradient and the inexact Hessian.
Particularly, iteration complexities of the Riemannian (deterministic)
adaptive cubic regularization algorithm and the Riemannian inexact adaptive
cubic regularization algorithm under the true gradient are established. As an
application, the proposed algorithms are applied to solve the joint
diagonalization problem on the Stiefel manifold. Numerical results indicate
that inexact algorithms are more efficient than the deterministic algorithm at
the same accuracy. Meanwhile, the inexact Riemannian adaptive cubic
regularization algorithm outperforms the inexact Riemannian trust-region
algorithm (in [17, Algorithm 1]).

The paper is organized as follows. As usual, some basic notions and notation
on Riemannian manifolds, together with some related properties about the sub-
sampling, are introduced in the next section. The Riemannian inexact adaptive
cubic regularization algorithm and its iteration complexity are presented in
section 3, and the application to joint diagonalization on the Stiefel
manifold is shown in the last section.

##  2 Notation and Preliminaries

Notation and terminologies used in the present paper are standard; the readers
are referred to some textbooks for more details; see, e.g., [12, 18, 19, 22].

Let $\mathcal{M}$ be a connected and complete $n$-dimensional Riemannian
manifold. Let $x\in\mathcal{M}$, and let ${T}_{x}\mathcal{M}$ stand for the
tangent space at $x$ to $\mathcal{M}$.
$T\mathcal{M}:=\cup_{x\in\mathcal{M}}T_{x}\mathcal{M}$ is called the tangent
bundle on $\mathcal{M}$. We denote by $\langle,\rangle_{x}$ the scalar product
on $T_{x}\mathcal{M}$ with the associated norm $\|\cdot\|_{x}$, where the
subscript $x$ is sometimes omitted. For $y\in\mathcal{M}$, let
$\gamma:[0,1]\rightarrow\mathcal{M}$ be a piecewise smooth curve joining $x$
to $y$. Then, the arc-length of $\gamma$ is defined by
$l(\gamma):=\int_{0}^{1}\|{\gamma}^{\prime}(t)\|dt$, while the Riemannian
distance from $x$ to $y$ is defined by ${\rm d}(x,y):=\inf_{\gamma}l(\gamma)$,
where the infimum is taken over all piecewise smooth curves
$\gamma:[0,1]\rightarrow\mathcal{M}$ joining $x$ to $y$.

We use $\nabla$ to denote the Levi-Civita connection on $\mathcal{M}$. A
vector field $V$ is said to be parallel along $\gamma$ if
$\nabla_{{\gamma}^{\prime}}V=0$. In particular, for a smooth curve $\gamma$,
if ${\gamma}^{\prime}$ is parallel along itself, then $\gamma$ is called a
geodesic; thus, a smooth curve $\gamma$ is a geodesic if and only if
$\nabla_{{\gamma}^{\prime}}{{\gamma}^{\prime}}=0$. A geodesic
$\gamma:[0,1]\rightarrow M$ joining $x$ to $y$ is minimal if its arc-length
equals its Riemannian distance between $x$ and $y$. By the Hopf-Rinow theorem
[18], $(\mathcal{M},{\rm d})$ is a complete metric space, and there is at
least one minimal geodesic joining $x$ to $y$ for any points $x$ and $y$.

Let $f:\mathcal{M}\rightarrow{{\mathbb{R}}}$ be a twice continuously
differentiable real-valued function defined on $\mathcal{M}$. The Remannian
gradient and Hessian of $f$ at $x\in\mathcal{M}$ are denoted by ${\rm
grad}f(x)$ and ${\rm Hess}f(x)$, which are respectively defined as

$\langle{\rm grad}f(x),\xi\rangle=\xi(f)\quad\mbox{and}\quad{\rm
Hess}f(x)(\xi,\eta)=\langle\nabla_{\xi}{\rm
grad}f(x),\eta\rangle\quad\forall\xi\in T_{x}\mathcal{M}.$

On manifolds $\mathcal{M}$, we use the retraction defined as below to move in
the direction of a tangent vector, which can be found in [12, Definition 4.1.1
and Proposition 5.5.5].

######  Definition 2.1 (retraction and seconder-order retraction).

A retraction on a manifold $\mathcal{M}$ is a smooth mapping
$R:T\mathcal{M}\rightarrow\mathcal{M}$ with the following properties, where
$R_{x}$ is the restriction of $R$ on $T_{x}\mathcal{M}$:

(1) $R_{x}(0_{x})=x$, where $0_{x}$ denotes the zero element of
$T_{x}\mathcal{M}$.

(2) $DR_{x}(0_{x})=\mathrm{Id}_{T_{x}\mathcal{M}}$, where
$\mathrm{Id}_{T_{x}\mathcal{M}}$ is the identity mapping on $T_{x}\mathcal{M}$
with the canonical identification $T_{0_{x}}T_{x}\mathcal{M}\simeq
T_{x}\mathcal{M}$.

Moreover, $R$ is said to be a seconder-order retraction if it further
satisfies

$\frac{D^{2}}{dt^{2}}R_{x}(t\eta)\mid_{t=0}=0\quad\forall
x\in\mathcal{M}\mbox{ and }\forall\eta\in T_{x}\mathcal{M},$

where $\frac{D^{2}}{dt^{2}}\gamma$ denotes acceleration of the curve $\gamma$.

######  Remark 1.

(i) $\mathcal{M}$ admits a seconder-order retraction defined by the
exponential mapping, and one can find more general retractions on matrix
manifolds in [12].

(ii) In general, the Euclidean Hessian $\nabla^{2}f\circ R_{x}(0_{x})$ differs
from the Riemannian Hessian $\mathrm{Hess}f(x)$, while they are identical
under second-order retractions (see[21, Lemma 3.9]).

Let $J:=\\{1,2,\cdots,n\\}$. We consider the large-scale minimization (1):

$\min_{x\in\mathcal{M}}f(x)=\frac{1}{n}\sum_{i=1}^{n}f_{i}(x),$

where $n\gg 1$, and for each $i\in J$,
$f_{i}:\mathcal{M}\rightarrow\mathbb{R}$ is a twice continuously
differentiable real-valued function. Note that $f$ may be non-convex, some
authors refer to find a proximate
$(\varepsilon_{g},\varepsilon_{H})$-optimality in practice, which is defined
by [32, Definition 1] (see also [17, Definition 2.1]). As usual,
$\lambda_{\min}(A)$ stands the minimum eigenvalue of matrix $A$.

######  Definition 2.2 ($(\varepsilon_{g},\varepsilon_{H})$-optimality).

Let $\varepsilon_{g},\varepsilon_{H}\in(0,1)$. A point $x^{*}\in\mathcal{M}$
is said to be an $(\varepsilon_{g},\varepsilon_{H})$-optimality of (1) if

$\|{\rm{grad}}f(x^{*})\|\leq\varepsilon_{g}\quad\mbox{and}\quad\lambda_{\min}({\rm{Hess}}f(x^{*}))\geq-\varepsilon_{H}.$

In practice, we adopt the sub-sampled inexact gradient and Hessian. To
proceed, let $S_{g},S_{H}\subset J$ be the sample collections with or without
replacement from $J$, and their cardinalities are denoted as $|S_{g}|$ and
$|S_{H}|$, respectively. The sub-sampled inexact gradient and Hessian are
defined as

$G:=\frac{1}{|S_{g}|}\sum_{i\in S_{g}}\mathrm{grad}f_{i}\quad\mbox{and}\quad
H:=\frac{1}{|S_{H}|}\sum_{i\in S_{H}}\mathrm{Hess}f_{i},$ (2)

respectively. The following lemma provides sufficient sample sizes of $S_{g}$
and $S_{H}$ to guarantee that the inexact $G$ and $H$ approximate ${\rm
grad}f$ and ${\rm Hess}f$ in a probabilistic way, which is taken form [17,
Theorem 4.1]. Here we set

$K^{\max}_{g}:=\max_{i\in
J}\sup_{x\in\mathcal{M}}\|\mathrm{grad}f_{i}(x)\|\quad\mbox{and}\quad
K^{\max}_{H}:=\max_{i\in J}\sup_{x\in\mathcal{M}}\|\mathrm{Hess}f_{i}(x)\|.$

######  Lemma 2.3.

Let $\delta,\delta_{g},\delta_{H}\in(0,1)$ and let $R$ be a seconder-order
retraction. Assume that the sampling is done uniformly at random to generate
$S_{g}$ and $S_{H}$, and let $G$ and $H$ be defined by (2). If the sample
sizes $|S_{g}|$ and $|S_{H}|$ satisfy

$|S_{g}|\geq\frac{32(K_{g}^{\max})^{2}(\log\frac{1}{\delta})+\frac{1}{4}}{\delta_{g}^{2}}\quad\mbox{and}\quad|S_{H}|\geq\frac{32(K_{H}^{\max})^{2}(\log\frac{1}{\delta})+\frac{1}{4}}{\delta_{H}^{2}},$

then the following estimates hold for any $x\in\mathcal{M}$ and any $\eta\in
T_{x}\mathcal{M}$:

$\mathrm{Pr}(\|G(x)-\mathrm{grad}f(x)\|\leq\delta_{g})\geq 1-\delta,$
$\mathrm{Pr}(\|(H(x)-\nabla^{2}f\circ
R_{x}(0_{x}))[\eta]\|\leq\delta_{H}\|\eta\|)\geq 1-\delta.$

##  3 Inexact Riemannian adaptive cubic regularization algorithm

In this section, we propose the following inexact Riemannian adaptive cubic
regularization algorithm for solving problem (1), which employs the inexact
Hessian and gradient at each iteration.

######  Algorithm 3.1.

(Inexact Riemannian adaptive cubic regularization)

Step 0. Choose $\varepsilon_{g},\varepsilon_{H},\rho_{TH}\in(0,1)$, and
$\gamma\in(1,+\infty)$.

Step 1. Choose $x_{0}\in\mathcal{M}$ and $\sigma_{0}\in(0,+\infty)$, and set
$k:=0$.

Step 2. Construct inexact gradient $G_{k}$ and approximate Hessian $H_{k}$ of
$f$ at $x_{k}$, respectively.

Step 3. If $\|G_{k}\|\leq\varepsilon_{g}$ and
$\lambda_{\min}(H_{k})\geq-\varepsilon_{H}$, then return $x_{k}$.

Step 4. Solve the following sub-problem approximately:

$\eta_{k}\approx{\underset{\eta\in
T_{x_{k}}\mathcal{M}}{{\arg\min}\,m_{k}(\eta)}:=\langle
G_{k},\eta\rangle+\frac{1}{2}\langle
H_{k}[\eta],\eta\rangle+\frac{1}{3}\sigma_{k}\|\eta\|^{3}.}$ (3)

Step 5. Set $\rho_{k}:=\frac{f(x_{k})-f\circ
R_{x_{k}}(\eta_{k})}{-m_{k}(\eta_{k})}$.

Step 6. If $\rho_{k}\geq\rho_{TH}$, then set $x_{k+1}=R_{x_{k}}(\eta_{k})$ and
$\sigma_{k+1}=\frac{\sigma_{k}}{\gamma}$; otherwise set $x_{k+1}=x_{k}$ and
$\sigma_{k+1}=\gamma\sigma_{k}$.

Step 7. Set $k:=k+1$ and go to Step 2.

######  Remark 2.

Algorithm 3.1 is an extension version of the adaptive cubic regularization
algorithm in Euclidean spaces. As we known, in Euclidean settings, the
algorithm was proposed and extensively studied by Cartis et al. in [23, 24].
Recently, for solving the “large-scale” separable problem (1), Wu et al.
proposed the adaptive cubic regularization algorithm with inexact Hessian (and
true gradient) in [32, Algorithm 2] and the iteration complexity of the
algorithm is proven to be
$\mathcal{O}(\max\\{\varepsilon_{g}^{-2},\varepsilon_{H}^{-3}\\})$.

Let $\\{x_{k}\\}$ be a sequence generated by Algorithm 3.1. We make the
following blanket assumptions throughout the paper:

(A1) There exists constant $L_{H}>0$ such that for all
$k\in\mathbb{N}$$,\emph{}f\circ R$ satisfies

$|f\circ
R_{x_{k}}(\eta_{k})-f(x_{k})-\langle\mathrm{grad}f(x_{k}),\eta_{k}\rangle-\frac{1}{2}\langle\nabla^{2}f\circ
R_{x_{k}}(0_{x_{k}})[\eta_{k}],\eta_{k}\rangle|\leq\frac{1}{2}L_{H}\|\eta_{k}\|^{3}.$

(A2) There exists constant $K_{H}>0$ such that

$\|H_{k}\|:=\sup_{\eta\in T_{x_{k}}\mathcal{M},\|\eta\|\leq
1}\langle\eta,H_{k}[\eta]\rangle\leq K_{H}\quad\forall k\in\mathbb{N}.$

(A3) There exist constants $\delta_{g},\delta_{H}\in(0,1)$ such that

$\begin{array}[]{lll}&\mbox{\bf(A3)-a}\quad\|G_{k}-\mathrm{grad}f(x_{k})\|\leq\delta_{g}\quad\forall
k\in\mathbb{N},\\\ &\mbox{\bf(A3)-b}\quad\|(H_{k}-\nabla^{2}f\circ
R_{x_{k}}(0_{x_{k}}))[\eta_{k}]\|\leq\delta_{H}\|\eta_{k}\|\quad\forall
k\in\mathbb{N}.\end{array}$

######  Remark 3.

(i) As shown in [21, Appendix B], Assumption (A1) is satisfied in the case
when $\mathcal{M}$ is compact and $R$ is a second-order retraction.

(ii) In practice, according to Lemma 2.3, we shall construct $\\{G_{k}\\}$ and
$\\{H_{k}\\}$ by sub-sampling to ensure Assumption (A3) in a probabilistic
way.

As stated in Step 4, we solve the sub-problem (3) approximately every
iteration. The most popular conditions for the approximate solution in the
literature are the Cauchy and Eigenpoint conditions; see, e.g., [23, 24, 32].
To ensure the convergence of the algorithm, we make the following similar
assumptions on the approximate solutions $\\{\eta_{k}\\}$:

(A4) $-m_{k}(\eta_{k})\geq-m_{k}(\eta_{k}^{C})$, and $-m_{k}(\eta_{k})\geq-
m_{k}(\eta_{k}^{E})\quad if\enspace\lambda_{\min}(H_{k})<0$ $\forall
k\in\mathbb{N}$,

where $\eta_{k}^{C}$ and $\eta_{k}^{E}$ are the approximate optimal solutions
of (3) along the negative gradient and the negative curvature directions,
respectively, and $\eta_{k}^{E}$ satisfies

$\langle\eta_{k}^{E},H_{k}[\eta_{k}^{E}]\rangle\leq\nu\lambda_{\min}(H_{k})\|\eta_{k}^{E}\|^{2}<0,\quad\nu\in(0,1);$

see more details in ([23, 24]).

Now, let $k\in\mathbb{N}$. Noting that the sub-problem (3) of Algorithm 3.1 is
actually posed on Euclidean spaces, the follow properties of $\eta_{k}^{C}$
and $\eta_{k}^{E}$ are valid (see [32, Lemmas 6,7]).

######  Lemma 3.1.

The following estimates for $\eta_{k}^{C}$ and $\eta_{k}^{E}$ hold:

$-m_{k}(\eta_{k}^{C})\geq\max\\{\frac{1}{12}||\eta_{k}^{C}||^{2}(\sqrt{K_{H}^{2}+4\sigma_{k}||G_{k}||}-K_{H}),\frac{||G_{k}||}{2\sqrt{3}}\min\\{\frac{||G_{k}||}{K_{H}},\sqrt{\frac{||G_{k}||}{\sigma_{k}}}\\}\\},$

(4)
$-m_{k}(\eta_{k}^{E})\geq\frac{\nu|\lambda_{\min}(H_{k})|}{6}\max\\{\|\eta_{k}^{E}\|^{2},\frac{\nu^{2}|\lambda_{\min}(H_{k})|^{2}}{\sigma_{k}^{2}}\\},$
(5)
$\|\eta_{k}^{C}\|_{x_{k}}\geq\frac{1}{2\sigma_{k}}(\sqrt{K_{H}^{2}+4\sigma_{k}\|G_{k}\|}-K_{H}),$
(6) $\|\eta_{k}^{E}\|\geq\frac{\nu|\lambda_{\min}(H_{k})|}{\sigma_{k}}.$ (7)

Lemma 3.2 below estimates the sufficient decrease in the objective function.

######  Lemma 3.2.

Assume (A1) and (A3). Then we have

$f\circ
R_{x_{k}}(\eta_{k})-f(x_{k})-m_{k}(\eta_{k})\leq\delta_{g}\|\eta_{k}\|+\frac{1}{2}\delta_{H}\|\eta_{k}\|^{2}+(\frac{L_{H}}{2}-\frac{\sigma_{k}}{3})\|\eta_{k}\|^{3}.$
(8)

###### Proof.

In view of (3), we have that

$\begin{split}&|f\circ
R_{x_{k}}(\eta_{k})-f(x_{k})-m_{k}(\eta_{k})+\frac{\sigma_{k}}{3}\|\eta_{k}\|^{3}|\\\
=&|f\circ R_{x_{k}}(\eta_{k})-f(x_{k})-\langle
G_{k},\eta_{k}\rangle-\frac{1}{2}\langle H_{k}[\eta_{k}],\eta_{k}\rangle|\\\
\leq&|f\circ R_{x_{k}}(\eta_{k})-f(x_{k})-\langle{\rm
gradf}(x_{k}),\eta_{k}\rangle-\frac{1}{2}\langle\nabla^{2}f\circ
R_{x_{k}}(0_{x_{k}})[\eta_{k}],\eta_{k}\rangle|\\\ &+|\langle G_{k}-{\rm
gradf}f(x_{k}),\eta_{k}\rangle|+\frac{1}{2}|\langle(\nabla^{2}f\circ
R_{x_{k}}(0_{x_{k}})-H_{k})[\eta_{k}],\eta_{k}\rangle|\\\
\leq&\delta_{g}\|\eta_{k}\|+\frac{\delta_{H}}{2}\|\eta_{k}\|^{2}+\frac{L_{H}}{2}\|\eta_{k}\|^{3},\end{split}$

where the second inequality is from Assumptions (A1) and (A3). Then, (8)
follows immediately. ∎

######  Lemma 3.3.

Assume (A1)-(A4). Assume further that $\sigma_{k}\geq 2L_{H}$ and
$\delta_{g}\leq\frac{9\delta_{H}^{2}}{4\sigma_{k}}$. If one of the following
conditions is satisfied, then parameter $\sigma_{k}$ decreases, that is,
$\sigma_{k+1}=\frac{\sigma_{k}}{\gamma}\;(\gamma>1)$:

  1. (a)

$\|G_{k}\|>\varepsilon_{g}$,
$\delta_{H}\leq\min\\{\frac{1}{18},\frac{1-\rho_{TH}}{9}\\}(\sqrt{K_{H}^{2}+4\sigma_{k}\varepsilon_{g}}-K_{H})$.

  2. (b)

$\lambda_{\min}(H_{k})<-\varepsilon_{H}$,
$\delta_{H}\leq\min\\{\frac{1}{9},\frac{2(1-\rho_{TH})}{9}\\}\nu\varepsilon_{H}$.

###### Proof.

In light of (8) and $\sigma_{k}\geq 2L_{H}$, we get that

$\displaystyle f\circ R_{x_{k}}(\eta_{k})-f(x_{k})-m_{k}(\eta_{k})$
$\displaystyle\leq(\frac{L_{H}}{2}-\frac{\sigma_{k}}{3})\|\eta_{k}\|^{3}+\frac{1}{2}\delta_{H}\|\eta_{k}\|^{2}+\delta_{g}\|\eta_{k}\|$
(9)
$\displaystyle\leq-\frac{\sigma_{k}}{12}\|\eta_{k}\|^{3}+\frac{1}{2}\delta_{H}\|\eta_{k}\|^{2}+\delta_{g}\|\eta_{k}\|.$

We claim that the following implication holds:

$\mbox{[$\|\eta_{k}\|\geq\frac{9\delta_{H}}{\sigma_{k}}]\Longrightarrow[f\circ
R_{x_{k}}(\eta_{k})-f(x_{k})-m_{k}(\eta_{k})\leq 0]$}.$ (10)

In fact, note that

$-\sigma_{k}t^{2}+6\delta_{H}t+12\delta_{g}\leq 0$ for all
$t\geq\frac{6\delta_{H}+\sqrt{36\delta_{H}^{2}+48\sigma_{k}\delta_{g}}}{2\sigma_{k}}=\frac{3\delta_{H}}{\sigma_{k}}(1+\sqrt{1+\frac{4\sigma_{k}\delta_{g}}{3\delta_{H}^{2}}})$.
(11)

In view of $\delta_{g}\leq\frac{9\delta_{H}^{2}}{4\sigma_{k}}$ (by
assumption), we see that $\frac{4\sigma_{k}\delta_{g}}{3\delta_{H}^{2}}\leq
3$, and then have

$\frac{3\delta_{H}}{\sigma_{k}}(1+\sqrt{1+\frac{4\sigma_{k}\delta_{g}}{3\delta_{H}^{2}}})\leq\frac{9\delta_{H}}{\sigma_{k}}.$

If $\|\eta_{k}\|\geq\frac{9\delta_{H}}{\sigma_{k}}$, we get from (11) that
$-\sigma_{k}\|\eta_{k}\|^{2}+6\delta_{H}\|\eta_{k}\|+12\delta_{g}\leq 0$.
Therefore, implication (10) is valid by (9).

Below, we show that the following inequality holds under either condition (a)
or (b):

$\rho_{k}\geq\rho_{TH}.$ (12)

Granting this and Step 6 of Algorithm 3.1, we conclude that
$\sigma_{k+1}=\frac{\sigma_{k}}{\gamma}$, completing the proof

Assume that condition (a) is satisfied. Then, we have

$\|\eta_{k}^{C}\|\geq\frac{1}{2\sigma_{k}}(\sqrt{K_{H}^{2}+4\sigma_{k}\varepsilon_{g}}-K_{H})\geq\frac{9\delta_{H}}{\sigma_{k}},$
(13)

where the first inequality is by (6) and the first item of condition (a), and
the last inequality thanks to the second item of condition (a). We claim that

$f\circ
R_{x_{k}}(\eta_{k})-f(x_{k})-m_{k}(\eta_{k})\leq\frac{1}{2}\delta_{H}\|\eta_{k}^{C}\|^{2}+\delta_{g}\|\eta_{k}^{C}\|.$
(14)

Indeed, in the case when $\|\eta_{k}\|\leq\|\eta_{k}^{C}\|$, (14) is evident
from (9). Otherwise, $\|\eta_{k}\|>\|\eta_{k}^{C}\|$. This, together with (13)
and implication (10), implies that $f\circ
R_{x_{k}}(\eta_{k})-f(x_{k})-m_{k}(\eta_{k})\leq 0$, and so claim (14) is
trivial. Using (13), (14) and
$\delta_{g}\leq\frac{9\delta_{H}^{2}}{4\sigma_{k}}$, we estimate that

$f\circ
R_{x_{k}}(\eta_{k})-f(x_{k})-m_{k}(\eta_{k})\leq\frac{1}{2}\delta_{H}\|\eta_{k}^{C}\|^{2}+\frac{\delta_{H}}{4}\frac{9\delta_{H}}{\sigma_{k}}\|\eta_{k}^{C}\|\leq\frac{3}{4}\delta_{H}\|\eta_{k}^{C}\|^{2}.$
(15)

Due to assumption (A4) and (4), there holds that

$-m_{k}(\eta_{k})\geq-
m_{k}(\eta_{k}^{C})\geq\frac{1}{12}\|\eta_{k}^{C}\|^{2}(\sqrt{K_{H}^{2}+4\sigma_{k}\varepsilon_{g}}-K_{H}).$

This together with (15), implies that

$1-\rho_{k}=\frac{f\circ
R_{x_{k}}(\eta_{k})-f(x_{k})-m_{k}(\eta_{k})}{-m_{k}(\eta_{k})}\leq\frac{\frac{3}{4}\delta_{H}\|\eta_{k}^{C}\|^{2}}{\frac{1}{12}\|\eta_{k}^{C}\|^{2}(\sqrt{K_{H}^{2}+4\sigma_{k}\varepsilon_{g}}-K_{H})}\leq
1-\rho_{TH},$

where the second inequality thanks to the second item of condition (a), which
shows (12).

Assume that condition (b) is satisfied. Then, we get from (7) and condition
(b) that

$\|\eta_{k}^{E}\|\geq\frac{\nu|\lambda_{\min}(H_{k})|}{\sigma_{k}}\geq\frac{\nu\varepsilon_{H}}{\sigma_{k}}\geq\frac{9\delta_{H}}{\sigma_{k}},$
(16)

Thus, we claim that

$f\circ
R_{x_{k}}(\eta_{k})-f(x_{k})-m_{k}(\eta_{k})\leq\frac{1}{2}\delta_{H}\|\eta_{k}^{E}\|^{2}+\delta_{g}\|\eta_{k}^{E}\|.$
(17)

Indeed, if $\|\eta_{k}\|\leq\|\eta_{k}^{E}\|$, claim (17) is clear by (9).
Otherwise $\|\eta_{k}\|>\|\eta_{k}^{E}\|$. This, together with (16) and
implication (10), implies that $f\circ
R_{x_{k}}(\eta_{k})-f(x_{k})-m_{k}(\eta_{k})\leq 0$. Therefore, claim (17)
holds trivially. In view of (16), (17) and
$\delta_{g}\leq\frac{9\delta_{H}^{2}}{4\sigma_{k}}$, we get that

$f\circ
R_{x_{k}}(\eta_{k})-f(x_{k})-m_{k}(\eta_{k})\leq\frac{1}{2}\delta_{H}\|\eta_{k}^{E}\|^{2}+\frac{\delta_{H}}{4}\frac{9\delta_{H}}{\sigma_{k}}\|\eta_{k}^{E}\|\leq\frac{3}{4}\delta_{H}\|\eta_{k}^{E}\|^{2}.$
(18)

Moreover, in view of $\lambda_{\min}(H_{k})<0$ by condition (b), we have from
assumption (A4) and (5) that

$-m_{k}(\eta_{k})\geq-
m_{k}(\eta_{k}^{E})\geq\frac{\nu|\lambda_{\min}(H_{k})|}{6}\|\eta_{k}^{E}\|^{2}(\sqrt{K_{H}^{2}+4\sigma_{k}\varepsilon_{g}}-K_{H}).$

Combining this and (18), there holds that

$1-\rho_{k}=\frac{f\circ
R_{x_{k}}(\eta_{k})-f(x_{k})-m_{k}(\eta_{k})}{-m_{k}(\eta_{k})}\leq\frac{\frac{3}{4}\delta_{H}\|\eta_{k}^{E}\|^{2}}{\frac{\nu|\lambda_{\min}(H_{k})|}{6}\|\eta_{k}^{E}\|^{2}}\leq
1-\rho_{TH},$

where the second inequality thanks to the second item of condition (b), and
then (12) is ture. The proof is complete. ∎

In the following lemma, we estimate the upper bound for the cubic
regularization parameters $\\{\sigma_{k}\\}$ before the algorithm terminates.
Here all iterations satisfying (12) are said to be successful and to be
unsuccessful otherwise.

######  Lemma 3.4.

Assume (A1)-(A4). Assume further that Algorithm 3.1 does not terminate at
$N$-th iteration, and the parameters satisfy
$\delta_{g}\leq\frac{9\delta_{H}^{2}}{4\sigma_{k}}$ and

$\delta_{H}\leq\min\\{\min\\{\frac{1}{18},\frac{1-\rho_{TH}}{9}\\}(\sqrt{K_{H}^{2}+4\sigma_{k}\varepsilon_{g}}-K_{H}),\min\\{\frac{1}{9},\frac{2(1-\rho_{TH})}{9}\\}\nu\varepsilon_{H}\\},$
(19)

then

$\sigma_{k}\leq\max\\{\sigma_{0},2\gamma L_{H}\\}\quad\forall k=1,2,\cdots,N.$
(20)

###### Proof.

By contradiction, we assume that the $k$-th iteration ($k\leq N$) is the first
unsuccessful iteration such that

$\sigma_{k+1}=\gamma\sigma_{k}\geq 2\gamma L_{H},$ (21)

which implies that $\sigma_{k}\geq 2L_{H}$. Since Algorithm 3.1 does terminate
at the $k$-th iteration, we have that either $\|G_{k}\|>\varepsilon_{g}$ or
$\lambda_{\min}(H_{k})<-\varepsilon_{H}$. Recalling that $\sigma_{k}\geq
2L_{H}$, $\delta_{g}\leq\frac{9\delta_{H}^{2}}{4\sigma_{k}}$ and (19), Lemma
3.3 is applicable to showing $\sigma_{k+1}=\frac{\sigma_{k}}{\gamma}$, which
contradicts with (21). The proof is complete. ∎

Without loss generality, we assume that $\sigma_{0}\leq 2\gamma L_{H}$.
Moreover, from now on, let $\mathcal{N}_{succ}$ and $\mathcal{N}_{fail}$ stand
the set of all the successful and unseccessful iterations of Algorithm 3.1,
respectively. The number of successful (unsuccessful) iterations is denoted by
$|\mathcal{N}_{succ}|$ $(|\mathcal{N}_{fail}|)$. The lemma below provides
estimation of the total number of successful iterations before the algorithm
terminates.

######  Lemma 3.5.

Assume that the assumptions made in Lemma 3.4 are satisfied. Then the
following estimate holds:

$|\mathcal{N}_{succ}|\leq\frac{f(x_{0})-f_{\min}}{\rho_{TH}\kappa_{\sigma}}\max\\{\varepsilon_{g}^{-2},\varepsilon_{H}^{-3}\\}+1,$
(22)

where $f_{\min}$ is the minimum of $f$ over $\mathcal{M}$ and

$\kappa_{\sigma}:=\\{\frac{\nu^{3}}{24\gamma^{2}L_{H}^{2}},\frac{1}{2\sqrt{3}}\min\\{\frac{1}{K_{H}},\frac{1}{\sqrt{2\gamma
L_{H}}}\\}\\}.$ (23)

###### Proof.

Noting that $\\{f(x_{k})\\}$ is monotonically decreasing, we have that

$\begin{array}[]{lll}f(x_{0})-f_{\min}&\geq\sum_{k\in\mathcal{N}_{succ}}\left(f(x_{k})-f(x_{k+1})\right)\\\
&\geq\sum_{k\in\mathcal{N}_{succ}}\rho_{TH}(-m_{k}(\eta_{k})),\end{array}$
(24)

where the second inequality is by Steps 4-5. Assume that $k\in N_{succ}$ and
Algorithm 3.1 does not terminate at the $k$-th iteration. Then we have either
$\|G_{k}\|>\varepsilon_{g}$ or $\lambda_{\min}(H_{k})<-\varepsilon_{H}$, and

$\sigma_{k}\leq 2\gamma L_{H}.$ (25)

by lemma 3.4(20) (noting $\sigma_{0}\leq 2\gamma L_{H}$). In the case when
$\|G_{k}\|>\varepsilon_{g}$, we get from (4) that

$\begin{array}[]{lll}-m_{k}(\eta_{k})&\geq\frac{\|G_{k}\|}{2\sqrt{3}}\min\\{\frac{\|G_{k}\|}{K_{H}},\sqrt{\frac{\|G_{k}\|}{\sigma_{k}}}\\}>\frac{\varepsilon_{g}}{2\sqrt{3}}\min\\{\frac{\varepsilon_{g}}{K_{H}},\sqrt{\frac{\varepsilon_{g}}{\sigma_{k}}}\\}\\\
&>\frac{\varepsilon_{g}^{2}}{2\sqrt{3}}\min\\{\frac{1}{K_{H}},\sqrt{\frac{1}{\sigma_{k}}}\\}>\frac{\varepsilon_{g}^{2}}{2\sqrt{3}}\min\\{\frac{1}{K_{H}},\sqrt{\frac{1}{2\gamma
L_{H}}}\\},\end{array}$ (26)

where the third inequality is by $\varepsilon_{g}\in(0,1)$, and the fourth
inequality is because of (25). Similarly, for the case when
$\lambda_{\min}(H_{k})<-\varepsilon_{H}$, we have from (5) that

$-m_{k}(\eta_{k})\geq\frac{\nu^{3}|\lambda_{\min}(H_{k})|^{3}}{6\sigma_{k}^{2}}\geq\frac{\nu^{3}\varepsilon_{H}^{3}}{24\gamma^{2}L_{H}^{2}}.$

This, together with (26), yields

$-m_{k}(\eta_{k})\geq\min\\{\frac{\varepsilon_{g}^{2}}{2\sqrt{3}}\min\\{\frac{1}{K_{H}},\frac{1}{\sqrt{2\gamma
L_{H}}}\\},\frac{\nu^{3}\varepsilon_{H}^{3}}{24\gamma^{2}L_{H}^{2}}\\}\geq\kappa_{\sigma}\min\\{\varepsilon_{g}^{2},\varepsilon_{H}^{3}\\},$
(27)

where the last inequality thanks to (23). Combing (24) and (27), we get that

$\displaystyle
f(x_{0})-f_{\min}\geq(|\mathcal{N}_{succ}|-1)\rho_{TH}\kappa_{\sigma}\min\\{\varepsilon_{g}^{2},\varepsilon_{H}^{3}\\}.$
(28)

Thus, (22) holds, completing the proof. ∎

We can now prove the main theorem, which indicates the iteration complexity of
Algorithm 3.1.

######  Theorem 3.6 (complexity of Algorithm 3.1).

Assume that the assumptions made in Lemma 3.4 are satisfied. Then, Algorithm
3.1 terminates after at most
$\mathcal{O}(\max\\{\varepsilon_{g}^{-2},\varepsilon_{H}^{-3}\\})$ iterations.

###### Proof.

Noting that the total number of iteration of Algorithm 3.1
$N=|\mathcal{N}_{succ}|+|\mathcal{N}_{fail}|$, we have by definition of
Algorithm 3.1 that

$\sigma_{N}=\sigma_{0}\gamma^{|\mathcal{N}_{fail}|-|\mathcal{N}_{succ}|}.$

Noting $\sigma_{k}\leq 2\gamma L_{H}$ for each $k=1,2,\cdots N$ by Lemma
3.4(20), it follows that

$|\mathcal{N}_{fail}|\leq|\mathcal{N}_{succ}|+\log_{\gamma}\frac{2\gamma
L_{H}}{\sigma_{0}}.$ (29)

This, together with Lemma 3.5(22), implies that

$\displaystyle N=|\mathcal{N}_{succ}|+|\mathcal{N}_{fail}|$ $\displaystyle\leq
2|\mathcal{N}_{succ}|+\log_{\gamma}{\frac{2\gamma L_{H}}{\sigma_{0}}}$
$\displaystyle\leq\frac{2(f(x_{0})-f_{\min})}{\rho_{TH}\kappa_{\sigma}}\max\\{\varepsilon_{g}^{-2},\varepsilon_{H}^{-3}\\}+1+\log_{\gamma}{\frac{2\gamma
L_{H}}{\sigma_{0}}},$

completing the proof. ∎

As corollaries of Theorem 3.6, we consider the (deterministic) Riemannian
adaptive cubic regularization algorithm and the inexact Riemannian adaptive
cubic regularization algorithm under the true gradient. The following
corollary shows the iteration complexity of the (deterministic) Riemannian
adaptive cubic regularization algorithm, which is immediate from Theorem 3.6
(noting that assumption (A3) is satisfied naturally).

######  Corollary 3.7.

Assume (A1)-(A2) and (A4). If Algorithm 3.1 employs

$G_{k}:=\mathrm{grad}f(x)\quad and\quad H_{k}:=\mathrm{Hess}f(x)\quad\forall
k\in\mathbb{N},$

then it terminates after at most
$\mathcal{O}(\max\\{\varepsilon_{g}^{-2},\varepsilon_{H}^{-3}\\})$ iterations.

Following the same line of proof of Theorem 3.6, we have the following
corollary about the iteration complexity of the inexact Riemannian adaptive
cubic regularization algorithm under the true gradient.

######  Corollary 3.8.

Assume (A1)-(A2), (A3)-b and (A4). Assume further that Algorithm 3.1 employs

$G_{k}:=\mathrm{grad}f(x)\quad\forall k\in\mathbb{N}$

and the parameter $\delta_{H}$ satisfies

$\delta_{H}\leq\min\\{\min\\{\frac{1}{12},\frac{1-\rho_{TH}}{6}\\}(\sqrt{K_{H}^{2}+4\sigma_{k}\varepsilon_{g}}-K_{H}),\min\\{\frac{1}{6},\frac{1-\rho_{TH}}{3}\\}\nu\varepsilon_{H}\\}.$

Then, the algorithm terminates after at most
$\mathcal{O}(\max\\{\varepsilon_{g}^{-2},\varepsilon_{H}^{-3}\\})$ iterations.

##  4 Application and numerical experiments

In this section, we shall apply Algorithm 3.1 to solve the joint
diagonalization (JD for short) problems on Stiefel manifolds, and provide some
numerical experiments to compare their performances. All numerical comparisons
are implemented in MATLAB on a 3.20 GHz Intel Core machine with 8 GB RAM.

We first introduce some notation and results about the Stiefel manifold and
the JD. The following notation and results about the Stiefel manifold can be
found in [12]. Let $r,d\in\mathbb{N}$, and let
$A=(a_{ij})\in\mathbb{R}^{d\times r}$ (we always assume that $r\leq d$). The
symmetric matrix of $A$ is denoted by $\mathrm{sym}(A)$, which is defined by
$\mathrm{sym}(A):=\frac{A^{T}+A}{2}$. The Frobenius norm for $A$ is defined as

$\|\mathrm{A}\|_{F}^{2}:=\sum_{i=1}^{d}\sum_{j=1}^{r}a_{ij}^{2}.$

We use ${\rm diag}(A)$ to denote the matrix formed by the diagonal elements of
$A$, and then have $\|\mathrm{diag(A)}\|_{F}^{2}=\sum_{j=1}^{r}a_{jj}^{2}$ .
The Stiefel manifold St$(r,d)$ is the set of all $r\times d$ orthogonal
matrices, i.e.,

$\mathrm{St}(r,d):=\\{X\in\mathbb{R}^{d\times r}:X^{T}X=I\\},$

where $I\in\mathbb{R}^{r\times r}$ is the identity matrix. Let
$U\in\mathrm{St}(r,d)$. The tangent space of St$(r,d)$ at $U$ is defined by

$T_{U}\mathrm{St}(r,d):={\\{\xi:\xi^{T}U+U^{T}\xi=0\\}}.$

The Riemannian inner product on $T_{U}\mathrm{St}(r,d)$ is given by

$\langle\eta,\xi\rangle_{U}:={\rm trace}(\eta^{T}\xi)\quad\forall\eta,\xi\in
T_{U}\mathrm{St}(r,d).$

In our practice, we adopt the following retraction defined by

$R_{U}(\xi):={\rm qf}(U+\xi)\quad\forall U\in\mathrm{St}(r,d),$

where ${\rm qf}(U+\xi)$ means the orthogonal factor based on the QR
decomposition of $U+\xi$.

The following introduction about the blind source separation (BSS for short)
and the JD can be found in the works in [9, 17, 25, 26] and references there
in. Let $x(t)$ be a given $d$-dimensional stochastic process and $t$ the time
index such that

$x(t)=As(t)$

where $A=(a_{ij})\in\mathbb{R}^{d\times r}$ is a mixing matrix and $s(t)$ is
the $r$-dimensional source vector, both of $A$ and $s(t)$ are unknown. The BSS
problem tries to recover the mixing matrix $A$ and the source vector$s(t)$. If
we pose additional constraints on $s(t)$ (see, e.g., [9]), the BSS problem can
be transformed into the JD problem on the Stiefel manifold, which is stated as
follows:

$\min_{{U}\in\mathrm{St}(r,d)}f_{\mathrm{{ica}}}{(U)}:=-\frac{1}{n}\sum_{i=1}^{n}\|\mathrm{diag}({U^{T}C}_{i}{U})\|_{F}^{2},$
(30)

where $\\{C_{i}\in\mathbb{R}^{d\times d}:i=1,2,\cdots,n\\}$ are matrices
estimated from the dada under some source conditions (see [9, 26] for more
details). The following expressions of the (Riemannian) gradient and Hessian
of $f{\rm{{}_{ica}}}$ are taken form [17]. Let ${U}\in\mathrm{St}(r,d)$. The
Riemannian gradient of $f\mathrm{{}_{ica}}(U)$ at $U$ is given by

$\mathrm{grad}f\mathrm{{}_{ica}}(U)=\mathrm{P}_{U}{\rm
egrad}f\mathrm{{}_{ica}}(U)=\mathrm{P}_{U}(-\frac{1}{n}\sum_{i=1}^{n}4{C}_{i}U\mathrm{diag}({U^{T}C}_{i}{U})),$

where $\mathrm{egrad}f\mathrm{{}_{ica}}(U)$ is the Euclidean gradient of
$f\mathrm{{}_{ica}}(U)$ and the orthogonal projection $\mathrm{P}_{U}$ is

$\mathrm{P}_{U}(W):={W-U{\rm sym}(U^{T}W)\quad\forall W}\in\mathbb{R}^{d\times
r}.$

The Riemannian Hessian of $f\mathrm{{}_{ica}}(U)$ is defined by

$\displaystyle\mathrm{Hess}f\mathrm{{}_{ica}}(U)[\xi]=$
$\displaystyle\mathrm{P}_{U}({\rm
Degrad}f_{\mathrm{{ica}}}{(U)[\xi]}-\xi\mathrm{sym}({U^{T}{\rm
egrad}}f\mathrm{{}_{ica}}(U))$ $\displaystyle-U\mathrm{sym}({\xi^{T}{\rm
egrad}}f_{\mathrm{{ica}}}{(U)})-U\mathrm{sym}({U^{T}{\rm
Degrad}}f_{\mathrm{{ica}}}{(U)[\xi]})\quad\forall\xi\in
T{{}_{U}}\mathrm{St}(r,d),$

where

$\mathrm{Degrad}f_{\mathrm{{ica}}}{(U)[\xi]}=-\frac{1}{n}\sum_{i=1}^{n}4{C}_{i}({\xi{\rm
diag}(U^{T}C}_{i}{U})+{U{\rm diag}}({\xi^{T}C}_{i}{U})+{U{\rm
diag}}({U^{T}C}_{i}{\xi})).$

For simplicity, Algorithm 3.1 with true gradient and true Hessain, with true
gradient and inexact Hessain, and with inexact gradient and inexact Hessain
are denoted by RACR, SRACR and SSRACR, respectively; and the inexact
Riemannian trust-region algorithm in [17, Algorithm 1] is denoted by SSRTR. We
shall apply RACR, SRACR, SSRACR, and SSRTR to solve the JD problem (30), and
provide some numerical comparisons with different $(n,d,r)$. In all numerical
experiments, we generate $\\{{C}_{i}:i=1,2,\cdots,n\\}$ randomly, and all
related inexact gradients and inexact Hessains are constructed by sampling
with sizes $|S_{g}|=\frac{n}{4}$ and $|S_{H}|=\frac{n}{40}$. We set the
parameters as $\sigma_{0}=0.001,\rho_{TH}=0.9$ and $\gamma=2$ in RACR, SRACR
and SSRACR, and set $\triangle_{0}=1,\rho_{TH}=0.9$ and $\gamma=2$ in SSRTR
(see [17, Algorithm 1] for details about the papmeters). The stopping
criterion of all executed algorithms is set as $\|G_{k}\|^{2}\leq 0.001$.

The numerical results are listed in Table 4, where $i$ and $T$ denote the
number of iterations and the CPUtime (in seconds), respectively. The
observations from Table 4 reveal that SSRACR costs much less CPUtime than
SRACR, RACR and SSRTR at the same accuracy.

Number of iterations and time for RACR, SRACR, SSRACR, and SSRTR.  RACR SRACR
SSRACR SSRTR $(n,d,r)$ $i$ $T$ $i$ $T$ $i$ $T$ $i$ $T$ (2015,5,5) 317 7.524
788 9.349 220 2.259 398 4.042 (2015,10,10) 1073 48.144 1247 29.432 718 14.048
1936 39.174 (2015,20,20) 2209 108.753 2114 87.839 2570 83.804 5140 172.033
(2015,30,30) 2693 395.710 2398 177.815 2603 202.081 12739 836.978 (2015,43,43)
2913 1280.422 3169 736.996 3563 691.293 12777 2622.408 (7200,43,43) 2873
5526.860 3019 2485.531 2909 2005.734 12837 8501.568 (60000,43,43) 298 4448.468
298 2395.300 303 2118.032 861 4579.865

## Funding

Research of the second author was supported in part by the National Natural
Science Foundation of China (grant 12161017) and the Guizhou Provincial
Science and Technology Projects (grant ZK[2022]110).

## References

  * [1] Shalev-Shwartz S, Ben-David S. Understanding Machine Learning: From Theory to Algorithms. Cambridge university press, 2014.
  * [2] Roosta-Khorasani F, Doel K, Ascher U. Data completion and stochastic algorithms for PDE inversion problems with many measurements. Electron. T. Numer. Ana., 42:177–196, 2014.
  * [3] Roosta-Khorasani F, Doel K, Ascher U. Stochastic algorithms for inverse problems involving pdes and many measurements. SIAM J. Sci. Comput., 36(5):S3–S22, 2014.
  * [4] Robbins H, Monro S. A stochastic approximation method. AMS, 22(3):400–407, 1951.
  * [5] Roux N, Schmidt M, Bach F. A stochastic gradient method with an exponential convergence rate for finite training sets. Advances in neural information processing systems, 25, 2012.
  * [6] Johnson R, Zhang T. Accelerating stochastic gradient descent using predictive variance reduction. In NIPS, 1(3):315–323, 2013.
  * [7] Liu DC, Jorge Nocedal J. On the limited memory bfgs method for large scale optimization. Math. Program., 45(1-3):503–528, 1989.
  * [8] Mishra B, Kasai H, Jawanpuria P, Saroop A. A Riemannian gossip approach to subspace learning on grassmann manifold. Mach. Learn., 108:1783–1803, 2019.
  * [9] Theis FJ, Cason TP, Absil PA. Soft dimension reduction for ICA by joint diagonalization on the Stiefel manifold. In ICA, 2019.
  * [10] Bonnabel S. Stochastic gradient descent on Riemannian manifolds. IEEE T. Automat. Contr., 58(9):2217–2229, 2013.
  * [11] Sato H, Kasai H, Mishra B. Riemannian stochastic variance reduced gradient. ArXiv preprint arXiv:1702.05594, 2017.
  * [12] Absil PA, Mahony R, Sepulchre R. Optimization Algorithms on Matrix Manifolds. Princeton University Press, 2008.
  * [13] Huang W, Gallivan KA, Absil PA. A broyden class of quasi-Newton methods for Riemannian optimization. SIAM J. Optimiz., 25(3):1660–1685, 2015.
  * [14] Gabay D. Minimizing a differentiable function over a differential manifold. J. Optimiz. Theory App., 37:177–219, 1982.
  * [15] Ring W, Wirth B. Optimization methods on Riemannian manifolds and their application to shape space. SIAM J. Optimiz., 22(2):596–627, 2012.
  * [16] Absil PA, Baker CG, Gallivan KA. Trust-Region Methods on Riemannian Manifolds. Found. Comput. Math., 2007;7(3): 303–330.
  * [17] Kasai H, Mishra B. Inexact trust-region algorithms on Riemannian manifolds. Advances in neural information processing systems, 31, 2018.
  * [18] DoCarmo MP. Riemannian Geometry. Boston MA: Birkhäuser Boston; 1992.
  * [19] Udriste C. Convex Functions and Optimization Methods on Riemannian Manifolds. In: Mathematics and Its Applications, Kluwer Academic, Dordrecht, 1994.
  * [20] Jorge N, Stephen JW. Numerical Optimization. Springer Series in Operations Research, Springer-Verlag, New York, 1999.
  * [21] Boumal N, Absil PA, Cartis C. Global rates of convergence for nonconvex optimization on manifolds. IMA J. Numer. Anal., 39(1):1–33, 2019.
  * [22] Conn AR, Gould NIM, Toint PL. Trust Region Methods. Society for Industrial and Applied Mathematics, 2000.
  * [23] Cartis C, Gould NIM, Toint PL. Adaptive cubic regularisation methods for unconstrained optimization. Part I: Motivation, convergence and numerical results. Math. Program., 127(2):245–295, 2011.
  * [24] Cartis C, Gould NIM, Toint PL. Adaptive cubic regularisation methods for unconstrained optimization. Part II: Worst-case function and derivative-evaluation complexity. Math. Program., 30(2):295–319, 2011.
  * [25] Hyvärinen A, Oja E. Independent component analysis: Algorithms and applications. Neural networks, 13(4-5): 411–430, 2000.
  * [26] Theis FJ, Inouye Y. On the use of joint diagonalization in blind signal processing. In: Proc. ISCAS 2006, Kos, Greece 2006.
  * [27] Bottou L, Curtis FE, Nocedal J. Optimization mehtods for large-scale machine learning. SIAM Rev., 60(2): 223–311, 2018.
  * [28] Byrd RH, Chin GM, Neveitt W, Nocedal J. On the use of stochastic Hessian information in optimization methods for machine learning. SIAM J. Optim., 21(3): 977–995, 2011.
  * [29] Erdogdu MA, Montanari A. Convergence rates of sub-sampled Newton methods. In NIPS, 2015.
  * [30] Kohler JM, Lucchi A. Sub-sampled cubic regularization for non-convex optimization. In ICML, 2017.
  * [31] Yao Z, Xu P, Roosta-Khorasani F, Mahoney MW. Inexact non-convex Newton-type methods. ArXiv preprint arXiv: 1802.06925, 2018.
  * [32] Xu P, Roosta F, Mahoney MW. Newton-type methods for non-convex optimization under inexact Hessian information. Math. Program., 184(1-2):35–70, 2020.
  * [33] Wang XM. Subgradient algorithms on Riemannian manifolds of lower bounded curvatures. Optimization, 67: 179–194, 2018.
  * [34] Wang XM, Li C, Wang JH, and Yao JC. Linear convergence of subgradient algorithm for convex feasibility on Reimannian manifolds. SIAM J. Optim., 25: 2334–2358, 2015.
  * [35] Ferreira OP, Louzeiro MS, Prudente LF. Gradient method for optimization on Riemannian manifolds with lower bounded curvature. SIAM J. Optim., 29: 2517–2541, 2019.
  * [36] Bento CG, Bitar SDB, Cruz Neto JX, Oliveira PR, Souza JCO. Computing Riemannian center of mass on Hadamard manifolds. J. Optim. Theory Appl., 183: 977–992, 2019.