\section{Introduction}

Automatic synthesis of visually pleasant texture that resembles
exemplary texture finds applications in computer graphics. Texture
synthesis has been studied for several decades since it is also of
theoretical interest in texture analysis and modeling. Texture can be
synthesized pixel-by-pixel \cite{de1997multiresolution,
efros1999texture, wei2000fast} or patch-by-patch~\cite{efros2001image,
liang2001real, cohen2003wang,kwatra2003graphcut,wu2004feature} based on
an exemplary pattern. For the pixel-based synthesis, a pixel conditioned
on its squared neighbor was synthesized using the conditional
probability and estimated by a statistical method in
\cite{efros1999texture}. Generally, patch-based texture synthesis yields
higher quality than pixel-based texture synthesis. Yet, to search the
whole image for patch-based synthesis is extremely slow
\cite{efros1999texture, liang2001real}. To speed up, one can stitch
small patches of the exemplary texture to form a larger region
~\cite{efros2001image, cohen2003wang, kwatra2005texture}. Although these
methods can produce texture of higher quality, its diversity is limited.
Besides texture synthesis in the spatial domain, one can transform
texture images from the spatial domain to the spectral domain with
certain filters (or kernels) and exploit the statistical correlation of
filter responses for texture synthesis. Commonly used kernels include
the Gabor filters~\cite{heeger1995pyramid} and the steerable pyramid
filter banks~\cite{portilla2000parametric}.

We have witnessed amazing quality improvement of synthesized texture in
the last 5-6 years due to the resurgence of neural networks. Texture
synthesis based on deep learning (DL), such as convolutional neural
networks (CNNs) and generative adversarial networks (GANs), yield
visually pleasant results. DL-based methods learn transform kernels from
numerous training data~\cite{gatys2015texture, liu2016texture,
risser2017stable, li2017diversified, li2017universal,
ustyuzhaninov2017does, shi2020fast} through end-to-end optimization.
However, these methods have two main shortcomings: 1) lack of
mathematical transparency and 2) higher training and inference
complexity.  To address them, we investigate a non-parametric and
interpretable texture synthesis method, called NITES, in this work.

NITES consists of three steps. First, it analyzes the exemplary texture
to obtain its joint spatial-spectral representations.  Second, the
probabilistic distributions of training samples in the joint
spatial-spectral spaces are characterized. Finally, new texture images
are generated by mimicking probabilities of source texture images. In
particular, we adopt a data-driven transform, known as the channel-wise
(c/w) Saab transform~\cite{chen2020pixelhopplus}, which provides a
powerful representation in the joint spatial-spectral space.  The c/w
Saab transform is derived from the successive subspace learning (SSL)
theory~\cite{kuo2016understanding, kuo2017cnn,
kuo2019interpretable,chen2020pixelhop}. We will show that NITES can
generate high quality texture at lower complexity.

The rest of the paper is organized as follows. The framework of
successive subspace embedding and generation is described in Sec.
\ref{sec:SSL}. The NITES method is proposed in Sec.~\ref{sec:method}.
Experimental results are shown in Sec.~\ref{sec:experiments}. Concluding
remarks are given in Sec.~\ref{sec:conclusion}.

\section{Successive Subspace Embedding and Generation}\label{sec:SSL}

In this section, we would like to explain a high-level idea behind the
NITES method, which is called the successive subspace embedding and
generation principle. Consider an input signal space denoted by
$\tilde{S}_0$, and a sequence of subspaces denoted by $\tilde{S}_1,
\cdots, \tilde{S}_n$. Their dimensions are denoted by $\tilde{D}_0$,
$\tilde{D}_1$, $\cdots$, $\tilde{D}_n$.  They are related with each
other by the constraint -- any element in $\tilde{S}_{i+1}$ is formed by
an affine combination of elements in $\tilde{S}_{i}$, where $i=0,
\cdots, n-1$.

An affine transform can be converted to a linear transform by augmenting
vector $\tilde{{\bf a}}$ in $\tilde{S}_{i}$ via ${\bf a}=(\tilde{{\bf
a}}^T, 1)^T$. We use $S_i$ to denote the augmented space of
$\tilde{S}_{i}$ and $D_i=\tilde{D}_i+1$. Then, we have the following
embedding relationship
\begin{equation}\label{eq:embedding}
S_n \subset S_{n-1} \subset \cdots \subset S_1 \subset S_0,
\end{equation}
and
\begin{equation}\label{eq:dim}
D_n < D_{n-1} < \cdots <  D_1 < D_0.
\end{equation}
This concept is illustrated in Fig. \ref{fig:subspace} with $n=2$.

\begin{figure*}[t]
\centering
\includegraphics[width=0.9\linewidth]{fig/subspace.pdf}
\caption{Illustration of the concept of successive subspace embedding and
generation, where blue arrows indicates a sequence of subspace embedding
processes and red arrows indicates a sequence of subspace generation
processes.}\label{fig:subspace}
\end{figure*}

We use texture embedding and generation as an example to explain this
concept. To generate homogeneous texture, we collect a collection of
texture patches cropped out of exemplary texture as the input set.
Suppose that each texture patch has RGB three color channels, and a
spatial resolution $N\times N$. Then, the input set has a dimension of $3
N^2$ and its augmented space $S_0$ has a dimension of $D_0=3 N^2 + 1$.
If $N=32$, we have $D_0=3073$ which is too high to find an effective
generation model directly.

To address this challenge, we build a sequence of embedded subspaces
$S_0$, $S_1$, $\cdots$, $S_n$ with decreasing dimensions. We call $S_0$
and $S_n$ the "source" space and the "core" subspace, respectively.
We need to find an effective embedded subspace $S_{i+1}$ from $S_i$, and
such an embedding model is denote by $E_i^{i+1}$. Proper subspace
embedding is important since it determines how to decompose an input
space into the direct sum of two subspaces in the forward embedding path.
Although we choose one of the two for further processing, we also need
to record the probabilistic relationship of the two decomposed
subspaces so that samples of diversity and fidelity can be generated in
the reverse generation path.

In the reverse generation path, we begin with the generation of samples
in $S_n$ by studying its own statistics. This is accomplished by
generation model $G_n$. Then, conditioned on a generated sample in
$S_{i+1}$, we generate a new sample in $S_{i}$ through a generation
model denoted by $G_{i+1}^i$.  In Fig. \ref{fig:subspace}, we use blue
and red arrows to indicate a sequence of subspace embedding and
generation processes, respectively. This idea can be implemented as a
non-parametric method since we can choose subspaces $S_1$, $\cdots$,
$S_n$, flexibly in a feedforward manner. One specific design is
elaborated in the next section.

\section{Proposed NITES Method}\label{sec:method}

The proposed NITES method is presented in this section. To begin with,
we provide a system overview of the NITES method in
Sec.~\ref{subsec:system}. Next, we discuss the successive subspace
embedding (SSE) scheme based on the c/w Saab transform in
Sec.~\ref{subsec:sse}. After that, we examine the problem of core
subspace generation (CSG) in Sec.~\ref{subsec:csg}.  Finally, we
describe the successive subspace generation (SSG) process in
Sec.~\ref{subsec:ssg}.

\subsection{System Overview}\label{subsec:system}

An overview of the NITES method is given in Fig.~\ref{fig:framework}.
The exemplary color texture image has a spatial resolution of $256
\times 256$ and RGB three channels. We are interested in generating
multiple texture images that are visually similar to it. By randomly
cropping texture patches of size $32 \times 32$ out of the source image,
we obtain a collection of texture samples, which serves as the input to
the NITES system. The dimension of these patches is $32 \times 32
\times 3=3072$. Their augmented vectors form source space $S_0$. The
NITES system is designed to generate texture patches of the same size
that are visually similar to samples in $S_0$.  This is feasible if we
can capture both global and local patterns of these samples. There are
two paths in Fig.~\ref{fig:framework}. The blue arrows go from left to
right, denoting the successive subspace embedding process. The red
arrows go from right to left, denoting the successive subspace
generation process.  We can generate as many texture patches as desired
using this procedure. In order to generate a texture image of a larger
size, we perform image quilting~\cite{efros2001image} based on
synthesized patches.

\begin{figure*}[tb]
\centering
\includegraphics[width=1.0\linewidth]{fig/framework.pdf}
\caption{An overview of the proposed NITES method.}\label{fig:framework}
\end{figure*}

\subsection{Successive Subspace Embedding (SSE)}\label{subsec:sse}

The global structure of an image (or an image patch) can be well
characterized by spectral analysis. Yet, it is weak in capturing local
detail such as boundaries between regions. Joint spatial-spectral
representations offer an ideal solution to the description of both
global shape and local detail information.  Embedding model $E_0^1$ find
a proper subspace, $S_1$, in $S_0$ while embedding model $E_1^2$ find a
proper subspace, $S_2$, in $S_1$. As shown in Fig.~\ref{fig:framework}
NITES applies two-stage transforms.  They correspond to $E_0^1$ and
$E_1^2$, respectively.  Specifically, we can apply the c/w Saab
transform in each stage to achieve the successive subspace embedding
(SSE) task.  In the following, we provide a brief review on the Saab
transform~\cite{kuo2019interpretable} and the c/w Saab
transform~\cite{chen2020pixelhopplus}.

We partition each input patch into non-overlapping blocks, each of which
has a spatial resolution of $I_0 \times I_0$ with $K_0$ channels.  We
flatten 3D tensors into 1D vectors, and decompose each vector into the
sum of one DC and multiple AC spectral components. The DC filter is a
all-ones filter weighted by a constant. AC filters are obtained by
applying the principal component analysis (PCA) to DC-removed residual
tensor.  By setting $I_0=2$ and $K_0=3$, we have a tensor block of
dimension $2 \times 2 \times 3=12$.  Filter responses of PCA can be
positive or negative. There is a sign confusion problem
~\cite{kuo2016understanding, kuo2017cnn} if both of them are allowed to
enter the transform in the next stage. To avoid sign confusion, a
constant bias term is added to all filter responses to ensure that all
responses become positive - leading to the name of the "subspace
approximation with adjusted bias (Saab)" transform.  The Saab transform
is a date-driven transform, which is significantly different from
traditional transforms (e.g. Fourier and wavelet transforms) which are
data independent.  We partition AC channels into low- and high-frequency
two bands. The energy of high-frequency channels (shaded by gray color
in Fig.~\ref{fig:framework}) is low and they are discarded for dimension
reduction without affecting the performance much. The energy of
low-frequency channels (shaded by blue color in
Fig.~\ref{fig:framework}) is higher.  For a tensor of dimension 12, we
have one DC and 11 AC components. Typically, we select $K_1=6$ to 10
leading AC components and discard the rest. Thus, after $E_0^1$, one 12D
tensor becomes $K_1$-D vector, which is illustrated by dots in subspace
$S_1$.  The $K_1$-D response vectors are fed into the next stage for
another transform.

The channel-wise (c/w) Saab transform \cite{chen2020pixelhopplus} exploits
the weak correlation property between channels so that the Saab
transform can be applied to each channel separately (see the middle part
of Fig.~\ref{fig:framework}). The c/w Saab transform offers an improved
version of the standard Saab transform with a smaller model size.

One typical setting used in our experiments is shown below.
\begin{itemize}
    \item Dimension of the input patch ($\tilde{D}_0$): $32\times32\times3=3072$;
    \item Dimension of subspace $\tilde{S}_1$ ($\tilde{D}_1$):
    $16\times16\times10=2560$ (by keeping 10 channels in Hop-1);
    \item Dimension of subspace $\tilde{S}_2$ ($\tilde{D}_2$):
    $8\times8\times27=1728$ (by keeping 27 channels in Hop-2).
\end{itemize}
Note that the ratio between $\tilde{D}_1$ and $\tilde{D}_0$ is 83.3\%
while that between $\tilde{D}_2$ and $\tilde{D}_1$ is 67.5\%. We are
able to reduce the dimension of the source space to that of the core
subspace by a factor of 56.3\%. In the reverse path indicated by red
arrows, we need to develop a multi-stage generation process. It should
also be emphasized that users can choose channel numbers in Hop-1 and
Hop-2 flexibly. Thus, NITES is a non-parametric method.

The first-stage Saab transform provides the spectral information of the
nearest neighborhood, which is the first hop of the center pixel. By
generalizing from one to multiple hops, we can capture the information
in the short-, mid- and long-range neighborhoods. This is analogous to
increasingly larger receptive fields in deeper layers of CNNs. However,
filter weights in CNNs are learned from end-to-end optimization via
backpropagation while weights of the Saab filters in different hops are
determined by a sequence of PCAs in a feedforward unsupervised manner.

\subsection{Core Subspace Generation (CSG)}\label{subsec:csg}

As stated before, we begin with the generation of samples in $S_n$,
which is accomplished by generation model $G_n$, in the reverse
generation path. In the current case, $n=2$. We need to characterize the
sample statistics in core subspace $S_2$ for the purpose of synthesis.
The statistics of a large texture image can be complicated.  We leverage
the semi-periodic property of texture so as to focus on patches of size
$32 \times 32$, whose statistics can be analyzed more easily.

After the two-stage c/w Saab transform, the dimension of $S_2$ is
typically less than 2000. They are in form of c/w Saab coefficients. We
flatten these coefficients into a 1D vector, denoted by $\mathbf{z}$,
which is a random vector in $S_2$. To simplify the distribution
characterization of a high-dimensional random vector, we cluster
training samples into clusters and transform random vectors in each
cluster to a set of independent random variables. This is inspired by
the divide-and-conquer principle; namely, breaking down a difficult
problem into several easier sub-problems and solving them individually.

We adopt a hierarchical K-Means clustering algorithm
\cite{nister2006scalable} to cluster the training samples of
$\mathbf{z}$ into $N$ clusters, which are denoted by $\{C_i\}$,
$i=0,\cdots,N-1$. Rather than modeling probability $P(\mathbf{z})$
directly, we model condition probability $P(\mathbf{z}\mid
\mathbf{z}\in C_i)$ with a fixed cluster index. The probability,
$P(\mathbf{z})$, can be written as
\begin{equation}\label{eq:prob}
P(\mathbf{z}) = \sum\limits_{i=0}^{N-1} P(\mathbf{z}\mid \mathbf{z}\in
C_i) \cdot P(\mathbf{z}\in C_i),
\end{equation}
where $P(\mathbf{z}\in C_i)$ is the percentage of data points in cluster
$C_i$. It is abbreviated as $p_i$, $i=0,\dots,N-1$ (see the right part of
Fig.~\ref{fig:framework}).

Typically, a set of independent Gaussian random variables is used for
image generation. To do the same, we need to convert a collection of
correlated random vectors into a set of independent Gaussian random
variables. To achieve this objective, we transform random vector
$\mathbf{z}$ in cluster $C_i$ to a set of independent random variables
through independent component analysis (ICA), where non-Gaussianity
serves as an indicator for statistical independence. ICA finds
applications in noise reduction~\cite{hyvarinen1999sparse}, face
recognition~\cite{bartlett2002face}, and image
infusion~\cite{mitianoudis2007pixel}.

Our implementation is detailed below.
\begin{enumerate}
    \item Apply PCA to $\mathbf{z}$ in cluster $C_i$ for dimension reduction and data whitening.
    \item Apply FastICA~\cite{hyvarinen2000independent}, which is conceptually simple, computational efficient and robust to outliers, to the PCA output.
    \item Compute the cumulative density function (CDF) of each ICA
    component of random vector $\mathbf{z}$ in each cluster based on its
    histogram of training samples.
    \item Match the CDF in Step 3 with the CDF of a Gaussian random variable
    (see the right part of Fig.  \ref{fig:framework}), where the inverse CDF
    is obtained by resampling between bins with linear interpolation.  To
    reduce the model size, we quantize N-dimensional CDFs, which have $N$
    bins, with vector quantization (VQ) and store the codebook of quantized CDFs.
\end{enumerate}

We encode $P(\mathbf{z}\in C_i)$ in Eq. (\ref{eq:prob}) to be the length
of a segment in the unit interval, $[0,1]$. All segments are
concatenated in order to build the unit interval. The segment index is
the cluster index. These segments is called the interval representation.
To draw a sample from subspace $S_2$, we use the uniform random number
generator to select a random number from interval $[0,1]$ to decide the
cluster index using the interval representation.

To generate a new sample in $S_2$, we perform the following steps:
\begin{enumerate}
    \item Select a random number from the uniform random number generator to decide the cluster index.
    \item Draw a set of sample independently from the Gaussian distribution.
    \item Conduct histogram matching on the generated Gaussian samples using the inverse CDFs in the chosen cluster.
    \item Repeat Steps 1-3 if the generated sample of Step 3 has no value larger than a pre-set threshold.
    \item Perform the inverse transform of ICA and the inverse transform of PCA.
    \item Reshape the 1D vector into a 3D tensor and this tensor is the generated sample in $S_2$.
\end{enumerate}

\subsection{Successive Subspace Generation (SSG)}\label{subsec:ssg}

In this section, we examine generation model $G_{i+1}^i$, whose role is
to generate a sample in $S_i$ given a generated sample in $S_{i+1}$.
The core subspace generation process can generate a sample that captures
the global structure of an image patch but lacks in local detail such as
boundaries between regions. Embedding model, $E_i^{i+1}$, provides a
good description of local detail, which is obtained by the c/w Saab
transform. Thus, in the generation process, we conduct the inverse c/w
Saab transform on the generated sample in $S_{i+1}$.

It is worthwhile to point out that both Saab and c/w Saab transforms
have the forward and the inverse transforms. If no high-frequency
channels are removed, they are lossless transforms. If some
high-frequency channels are removed, they become lossy transforms.  In
the current case, we use lossy transforms for dimension reduction in the
successive subspace embedding process and lossy inverse transforms in
the successive subspace generation process.  The inverse c/w Saab
transform is a deterministic one once the forward transform is
determined. It can be used for signal reconstruction.

In the following, we take generation model $G_2^1$ from $S_2$ and to
$S_1$ as an example to explain the generation process.  A generated
sample in $S_2$ can be partitioned into $K_1$ groups as shown in the
left part of Fig.~\ref{fig:ssg}. Each group of channels is composed by
one DC channel and several low-frequency AC channels. The $k$th group of
channels in $S_2$, whose number is denoted by $K_2^{(k)}$, is derived
from the $k$th channel in $S_1$. We apply the inverse c/w Saab transform
to each group individually. The inverse c/w Saab transform converts the
tensor at the same spatial location across $K_2^{(k)}$ channels
(represented by white dots in Fig.~\ref{fig:ssg}) in $S_2$ into a block
of size $I_i \times I_i$ (represented by the white parallelogram in
Fig.~\ref{fig:ssg}) in $S_1$, using the DC and AC components obtained in
the successive subspace embedding process.  After the inverse c/w Saab
transform, the Saab coefficients in $S_1$ form a generated sample in
$S_1$. The same procedure is repeated between $S_1$ and $S_0$.

\begin{figure}[tb]
\begin{center}
\includegraphics[width=0.7\linewidth]{fig/ssg.pdf}
\end{center}
\caption{Illustration of generation model $G_{i+1}^i$, which converts the
same spatial location of $K_{i+1}$ channels in subspace $S_{i+1}$ to a block of
size $I_i \times I_i$ of $K_i$ channels in space/subspace $S_i$.}\label{fig:ssg}
\end{figure}

Examples of several generated textures in core subspace $S_0$,
intermediate subspace $S_1$ and source space $S_0$ are shown in
Fig.~\ref{fig:between}, where we display the generated DC channel of
group 1 (the first channel in Fig.~\ref{fig:ssg}) in $S_2$ and $S_1$ and
the generated texture patch in $S_0$. These DC channels offer gray-scale
low-resolution patterns of a generated sample. NITES can generate
samples containing more local detail gradually.

\begin{figure*}[tb]
\begin{center}
\includegraphics[width=0.7\linewidth]{fig/between.pdf}
\end{center}
\caption{Examples of several generated textures in core subspace $S_0$,
intermediate subspace $S_1$ and source space $S_0$.}\label{fig:between}
\end{figure*}

\section{Experiments}\label{sec:experiments}

\subsection{Experimental Setup}

The following hyper parameters (see Fig. \ref{fig:framework}) are used in our experiments.
\begin{itemize}
    \item Input filter window size to Hop-1: $I_0=2$,
    \item Input filter window size to Hop-2: $I_1=2$,
    \item Selected channel numbers in Hop-1 ($K_1$): $6\sim10$,
    \item Selected channel numbers in Hop-2 ($K_2$): $20\sim30$.
\end{itemize}
The analysis filter window size is the same as the generation window
size. All windows are non-overlapping with each other. The actual
channel numbers $K_1$ and $K_2$ are texture dependent. That is, we
examine the energy distribution based on the PCA eigenvalue plot and
choose the knee point where the energy becomes nearly flat.

\subsection{An Example: Brick Wall Texture Generation}

We show generated \emph{brick\_wall} texture patches of size
$32\times32$ and $64\times64$ in Figs.~\ref{fig:synpatch}(a) and (c).
We performed two-stage c/w Saab transforms on $32\times32$ patches and
three-stage c/w Saab transforms on $64\times64$ patches, whose core
subspace dimensions are equal to 1728 and 4032, respectively. Patches in
these figures were synthesized by running the NITES method in one
hundred rounds. Randomness in each round primarily comes from two
factors: 1) random cluster selection, and 2) random seed vector
generation.

Generated patches retain the basic shape of bricks and the diversity of
brick texture. We observe some unseen patterns generated by NITES, which
are highlighted by red squared boxes in Figs.~\ref{fig:synpatch} (a) and
(c). As compared with generated $32\times32$ patches, generated
$64\times64$ patches could be blurry sometimes (e.g., the one in the
upper right corner) due to a higher source dimension.

As a non-parametric generation model, NITES has freedom in choosing
multiple settings under the same pipeline. For example, it can select
different channel numbers in $\tilde{S}_1$ and $\tilde{S}_2$ to derive
different generation results. Four settings are listed in
Table~\ref{table:settings}. The corresponding generation results are
shown in Fig.~\ref{fig:dimanalysis}. Dimensions decrease faster from
(a) to (d). The quality of generated results becomes poorer due to
smaller dimensions of the core subspace, $\tilde{S}_2$, and the
intermediate subspace, $\tilde{S}_1$.

\begin{table}[h]
\begin{center}
\caption{The settings of our four generation model.} \label{table:settings}
\begin{tabular}{cccc} \\ \hline
Setting & $\tilde{D}_0$ & $\tilde{D}_1$ & $\tilde{D}_2$ \\ \hline
a       & 3072          & 2560          & 2048 \\
b       & 3072          & 1536          & 768 \\
c       & 3072          & 1280          & 512 \\
d       & 3072          & 768           & 192 \\ \hline
\end{tabular}
\end{center}
\end{table}

To generate larger texture images, we first generate 5,000 texture
patches and perform image quilting~\cite{efros2001image} with them. The
results after quilting are shown in Figs.~\ref{fig:synpatch} (b) and
(d).  All eight images are of the same size, $256\times256$. They are
obtained using different initial patches for the image quilting process.
By comparing the two sets of stitched images, the global structure of
the brick wall is better preserved using larger patches (i.e. of size
$64 \times 64$) while its local detail is a little bit blurry sometimes.

\begin{figure*}[htb]
\begin{center}
\includegraphics[width=0.9\linewidth]{fig/synbrick.pdf}
\end{center}
\caption{Examples of generated $brick\_wall$ texture patches and
stitched images of larger sizes, where the image in the bottom-left
corner is the exemplary texture image and the patches highlighted by red
squared boxes are unseen patterns.}\label{fig:synpatch}
\end{figure*}

\begin{figure*}[htb]
\begin{center}
\includegraphics[width=0.9\linewidth]{fig/dimanalysis.pdf}
\end{center}
\caption{An illustration of how dimension growth affects the generation
results. The numbers above the figure indicates the dimensions of $S_0$,
$S_1$ and $S_2$, respectively.}\label{fig:dimanalysis}
\end{figure*}

\subsection{Performance Benchmarking with CNN-based Methods}

{\bf Visual Quality Comparison.} The quality of synthesized texture is
usually evaluated by human eyes. A diversity loss function was proposed
to measure texture diversity for CNN-based methods~\cite{shi2020fast, li2017diversified}. Yet, NITES dose not have a loss function.  Thus, we
show synthesis results of two CNN methods and NITES side by side in
Fig.~\ref{fig:synimg}.  Exemplary texture images are collected from
those in~\cite{gatys2015texture, ustyuzhaninov2017does, portilla2000parametric} or the Internet for the illustration purpose. The two benchmarking CNN methods were
proposed by Gatys {\em et al.}~\cite{gatys2015texture} and Ustyuzhaninov
{\em et al.} \cite{ustyuzhaninov2017does}. We run the codes provided by
them, and show their results in the second and third columns of
Fig.~\ref{fig:synimg}, respectively. We used the default setting of the iteration number, which is 2000 in~\cite{gatys2015texture} and 4000 in~\cite{gatys2015texture} for visualization. Two results generated by NITES are
shown in the last two columns. They are obtained in two different runs.
For texture \emph{meshed}, we see the brown fog artifact in
~\cite{gatys2015texture, ustyuzhaninov2017does}, which is visually
apparent. However, it does not show up in our generation.  As
demonstrated by these examples, NITES can generate high quality and
visually pleasant texture images.

\begin{figure*}[htb]
\begin{center}
\includegraphics[width=\textwidth, height=0.92\textheight, keepaspectratio]{fig/synimg.pdf}
\end{center}
\caption{Comparison of synthesized texture images using two CNN-based
methods and NITES (from left to right): exemplary texture images,
texture images generated by \cite{gatys2015texture} texture images
generated by \cite{ ustyuzhaninov2017does} and two texture images
generated by NITES.} \label{fig:synimg}
\end{figure*}

{\bf Comparison of Texture Generation Time.} It is worthwhile to compare
texture synthesis time of different texture image generation methods.
The results are shown in Table~\ref{table:time}. All experiments were
conducted in the same machine composed by 12 CPUs (Intel Core i7-5930K
CPU at 3.50GHz) and 1 GPU (GeForce GTX TITAN X).  No GPU was needed in
our model but it was needed in~\cite{gatys2015texture, ustyuzhaninov2017does}. We set the iteration number to 1000
for~\cite{gatys2015texture} and 100 for~\cite{ustyuzhaninov2017does}.
NITES generated 5,000 $32\times32$ texture patches, 81 of which were
stitched into the final texture image.  For all three methods, we
measure the time needed in generating one image of size $256\times 256$.
As shown in Table~\ref{table:time}, NITES generates one texture image
with 213.04 seconds while Gatys' method and Ustyuzhaninov's method
demand 513.98 and 949.64 seconds, respectively.  NITES is significantly
faster.

\begin{table}[htb]
\begin{center}
\caption{Comparison of time needed to generate one texture image.} \label{table:time}
\begin{tabular}{cc} \\ \hline
Methods                                  &  Time (seconds) \\ \hline
Gatys et al.~\cite{gatys2015texture}     &  513.98  \\
Multi-Scale~\cite{ustyuzhaninov2017does} &  949.64  \\
NITES (Ours)                             &  213.04  \\ \hline
\end{tabular}
\end{center}
\end{table}

{\bf Breakdown of NITES' Generation Time.} We can break down the texture
generation time of NITES into three parts: 1) successive subspace
embedding (i.e., the forward embedding path), 2) core and successive
subspace generations (i.e., the reverse generation path) and 3) the
quilting process. Time required by each part is shown in
Table~\ref{table:processtime}. They demand 24.24, 108.03 and 8.08
seconds, respectively. To generate multiple texture images from the same
exemplary texture, we need to run the first part once but the second and
third parts multiple times (one run for one new synthesis).  Thus, it is
fair to focus on the last two parts only for single texture image
generation, which is equal to 116 seconds.  In contrast, the two
benchmarking CNN methods do not have such a breakdown.  They need to go
through the whole pipeline to generate one new texture image.

\begin{table}[htb]
\begin{center}
\caption{The time of three processes in our method.}\label{table:processtime}
\begin{tabular}{cc} \\ \hline
Process           &   Time (seconds)  \\ \hline
Forware Embedding &   24.24  \\
Reverse Generation&  108.03 \\
Quilting          &    8.08 \\  \hline
\end{tabular}
\end{center}
\end{table}

\section{Conclusion and Future Work}\label{sec:conclusion}

A non-parametric interpretable texture synthesis (NITES) method was
proposed based on a new texture analysis and synthesis framework in this
work. Texture can be analyzed and represented effectively using the
multi-stage c/w Saab transforms that offer a sequence of joint
spatial-spectral representations. The sample distribution in the core
subspace was carefully studied, which allows us to build a core subspace
generation model. Furthermore, a successive subspace generation model
was developed to build a higher-dimensional subspace based on a
lower-dimensional subspace. As a result, new texture samples can be
generated by mimicking probabilities and/or conditional probabilities of
source texture patches. Extensive experimental results were conducted to
demonstrate the power of the proposed NITES method.  It can generate
visually pleasant texture images effectively, including some unseen
patterns.

As to future research extension, there are several directions worth
exploration. It appears to be an important topic in controlling the
growth of subspace dimensions in the generation process. Is it
beneficial to introduce more intermediate subspaces between the source
and the core? Can we apply the same model for other image generation
such as human faces, digits, scenes and objects? Is it possible to
generalize the framework to image inpainting? How to compare our
generation model and GAN? These are all open and interesting questions
for further investigation.

\section*{Acknowledgment}

This work was supported by MediaTek. Computation for the work was supported by the University of Southern
California's Center for High Performance Computing (hpc.usc.edu).

\bibliographystyle{unsrt}


\section{Introduction}
In Euclidean space, nonlinear conjugate gradient methods aim to solve unconstrained optimization problems.
Conjugate gradient methods have been developed by Hestenes and Stiefel \cite{hestenes1952methods} for solving linear systems whose coefficient matrix is symmetric positive-definite.
Fletcher and Reeves \cite{fletcher1964function} extended the conjugate gradient method to unconstrained nonlinear optimization problems.
Theirs is the first nonlinear conjugate gradient method in Euclidean space.
Since then, various nonlinear conjugate gradient methods have been proposed (see \cite{al1985descent, polak1969note, polyak1969conjugate, dai1999nonlinear});
these have summarized by Hager and Zhang in \cite{hager2006survey}.
A sufficient descent condition is used to analyze the global convergence of conjugate gradient methods with inexact line searches.
Hager and Zhang \cite{hager2005new} proposed a conjugate gradient method whose search direction satisfies the sufficient descent condition regardless of whether a line search is used or not.
In addition, Dai \cite{dai2010nonlinear} proposed a nonlinear conjugate gradient methods that are the generalizations of the Hager-Zhang method.
His method also satisfies the sufficient descent condition regardless of whether a line search is used or not.
A nonlinear conjugate gradient method that satisfies the sufficient descent condition is called a sufficient descent nonlinear conjugate gradient method.
Narushima and Yabe summarized the sufficient descent nonlinear conjugate gradient methods in \cite{narushima2014survey}.

The conjugate gradient method in Euclidean space is applicable to a Riemannian manifold.
In \cite{smith1994optimization}, Smith introduced the notion of Riemannian optimization. He used the exponential map and parallel transport to generalize the optimization method on a Riemannian manifold.
However, in general, using the exponential map or parallel transport on a Riemannian manifold is not computationally efficient.
Absil, Mahony, and Sepulchre \cite{absil2008optimization} proposed to use a mapping called a retraction that approximates the exponential map.
Moreover, they introduced the notion of vector transport, which approximates parallel transport.
Retractions and vector transports on Stiefel manifolds have been summarized and numerically compared by Zhu \cite{zhu2017riemannian}.

Ring and Wirth \cite{ring2012optimization} proposed the Fletcher-Reeves type of nonlinear conjugate gradient method on Riemannian manifolds with retraction and vector transport.
They indicated that the Fletcher-Reeves methods converges globally when each step size satisfies the strong Wolfe conditions \cite{wolfe1969convergence, wolfe1971convergence}.
However, their convergence analysis assumed that vector transport satisfies the Ring-Wirth nonexpansive condition (see \eqref{eq:nonexpansive} for the definition of the Ring-Wirth nonexpansive condition).
Vector tansports that do not satisfy this condition have also been used (see \cite[Section 5]{sato2013new}).
In \cite{sato2013new}, Sato and Iwai introduced the notion of scaled vector transport \cite[Definition 2.2]{sato2013new} to remove this impractical assumption from the convergence analysis.
They proved that by using scaled vector transport, the Fletcher-Reeves method on a Riemannian manifold generates a descent direction at every iteration
and converges globally without the Ring-Wirth nonexpansive condition. Similarly, Sato \cite{sato2015dai} used scaled vector transport in a convergence analysis.
He indicated that the Dai-Yuan-type Riemannian conjugate gradient method generates a descent direction at every iteration and converges globally under the Wolfe conditions.
In \cite{sakai2020hybrid}, Sakai and Iiduka proposed the hybrid Riemannian conjugate gradient method, which combines the Hestenes-Stiefel and Dai-Yuan methods.
They proved that by using scaled vector transport, this hybrid method generates a descent direction at every iteration and converges globally under the strong Wolfe conditions.

In this paper, we focus on the sufficient descent condition and sufficient descent conjugate gradient method on Riemannian manifolds.
We propose two kinds of sufficient descent nonlinear conjugate methods to Riemannian manifold.
One is a hybrid formula combining the Fletcher-Reeves method with the Polak-Ribi\`ere-Polyak method, and we prove that using scaled vector transport,
this hybrid method has the global convergence property under the strong Wolfe conditions.
The other is a formula that satisfies the sufficient descent condition regardless of whether a line search is used or not, and we prove that this method has this property even on Riemannian manifolds.
This formula is a generalization of the Hager-Zhang method defined on Euclidean space.
Moreover, we generalize two typical line search algorithms in Euclidean space, i.e., the backtracking line search and line search algorithm with a \emph{zoom} phase.
In numerical experiments, we compare the sufficient descent Riemannian conjugate gradient methods with the above two line search algorithms.
The numerical results show that the proposed hybrid method should use the step sizes satisfying the strong Wolfe conditions, which guarantee its convergence (Theorem \ref{thm:Hybrid2convergence}).
Moreover, the results show that the benefit of the Hager-Zhang-type method is its fast convergence property regardless of the type of line search used,
as promised by its sufficient descent property (Theorem \ref{thm:SDsufficient}).
The intellectual contribution of this paper is to show the fast convergence property of the sufficient descent Riemannian conjugate gradient method regardless of the type of line search used.

This paper is organized as follows. Section \ref{sec:RCGM} reviews the Riemannian conjugate gradient methods and some useful concepts.
Moreover, two Riemannian conjugate gradient methods are proposed in this section. Section \ref{sec:sufficient} proves that several Riemannian conjugate gradient methods satisfy the sufficient descent condition.
Section \ref{sec:LineSearch} generalizes two kinds of typical line search algorithms on Riemannian manifolds.
Section \ref{sec:NE} provides the numerical experiments on several Riemannian optimization problems.
Section \ref{sec:conclusion} concludes the paper.

\section{Riemannian Conjugate Gradient Methods}~\label{sec:RCGM}
Let $M$ be a Riemannian manifold and $T_xM$ be a tangent space at a point $x \in M$. $\ip{\cdot}{\cdot}_x : T_xM \times T_xM \rightarrow \mathbb{R}$ denotes a Riemannian metric at a point $x \in M$. The Riemannian gradient of a smooth function $f:M \rightarrow \mathbb{R}$ at $x \in M$ is denoted by $\grad f(x)$. Let $TM$ be the tangent bundle of $M$ and $\oplus$ be the Whitney sum. For a smooth mapping $F:M \rightarrow N$ between two manifolds $M$ and $N$, $\mathrm{D}F(x) : T_xM \rightarrow T_{F(x)}N$ denotes the differential of $F$ at $x \in M$ (see \cite[Section 3]{absil2008optimization}). An unconstrained optimization problem on a Riemannian manifold $M$ is expressed as follows:
\begin{problem}~\label{pbl:main}
Let $f:M\rightarrow\mathbb{R}$ be smooth. Then, we would like to
\begin{align*}
\textrm{minimize}& \quad f(x), \\
\textrm{subject to}& \quad x \in M.
\end{align*}
\end{problem}
In order to generalize line search optimization algorithms to Riemannian manifolds, we will use the notions of a retraction \cite[Chapter 4, Definition 4.1.1]{absil2008optimization} and vector transport \cite[Chapter 8, Definition 8.1.1]{absil2008optimization}, which are defined as follows:

\begin{definition}
[Retraction]~\label{def:retraction}
Any smooth map $R:TM{\rightarrow}M$ is called a retraction on $M$ if it has the following properties.
\begin{itemize}
\item $R_x(0_x)=x$, where $0_x$ denotes the zero element of $T_xM$;
\item With the canonical identification $T_{0_x}T_xM \simeq T_xM$, $R_x$ satisfies
\begin{align*}
\mathrm{D}R_x(0_x)[\xi]=\xi
\end{align*}
for all $\xi \in T_xM$,
\end{itemize}
where $R_x$ denotes the restriction of $R$ to $T_xM$.
\end{definition}

\begin{definition}
[Vector transport]~\label{def:vector transport}
Any smooth map $\mathcal{T}:TM \oplus TM \rightarrow TM:(\eta ,\xi) \mapsto \mathcal{T}_\eta(\xi)$ is called a vector transport on $M$ if it has the following properties.
\begin{itemize}
\item There exists a retraction $R$, called the retraction associated with $\mathcal{T}$, such that $\mathcal{T}_{\eta}(\xi) \in T_{R_x(\eta)}M$ for all $x \in M$, and for all $\eta ,\xi \in T_xM$;
\item $\mathcal{T}_{0_x}(\xi)=\xi$ for all $\xi \in T_xM$;
\item $\mathcal{T}_{\eta}(a\xi + b\zeta)=a\mathcal{T}_{\eta}(\xi)+b\mathcal{T}_{\eta}(\zeta)$ for all $a,b \in \mathbb{R}$, and for all $\eta ,\xi ,\zeta \in T_xM$.
\end{itemize}
\end{definition}
Retraction and vector transport are generalizations of the exponential map and parallel transport, respectively. We will use the Ring-Wirth nonexpansive condition \cite{ring2012optimization}, that is a vector transport $\mathcal{T}$ satisfying
\begin{align}~\label{eq:nonexpansive}
\norm{\mathcal{T}_\eta(\xi)}_{R_x(\eta)} \leq \norm{\xi}_x,
\end{align}
to establish global convergence for the Fletcher-Reeves type Riemannian conjugate gradient method. In this paper, we will focus on the differentiated retraction $\mathcal{T}^R$ of $R$ as a vector transport, defined by
\begin{align*}
\mathcal{T}^{R}_{\eta}(\xi):=\mathrm{D}R_{x}(\eta)[\xi],
\end{align*}
where $x \in M$ and $\eta , \xi \in T_xM$. Then, the retraction $R$ is associated with $\mathcal{T}^R$. However, the differentiated retraction $\mathcal{T}^R$ does not always satisfy the Ring-Wirth nonexpansive condition \eqref{eq:nonexpansive}. To overcome this difficulty, Sato and Iwai \cite{sato2013new} introduced the notion of scaled vector transport. Scaled vector transport $\mathcal{T}^S$ respect to a retraction $R$ is defined as
\begin{align}~\label{eq:scaled}
\mathcal{T}^{S}_{\eta}(\xi):=
\begin{cases}
\mathcal{T}^R_{\eta}(\xi), & \textrm{if}\, \norm{\mathcal{T}_{\eta}^R(\xi)}_{R_{x}(\xi)}\leq\norm{\eta}_{x}, \\
\dfrac{\norm{\eta}_{x}}{\norm{\mathcal{T}^R_{\eta}(\xi)}_{R_{x}(\xi)}}\mathcal{T}^R_{\eta}(\xi), & \textrm{otherwise}.
\end{cases}
\end{align}

The general framework of Riemannian conjugate gradient methods is described in Algorithm \ref{alg:conjugate}.
\begin{algorithm}
\caption{General framework of Riemannian conjugate gradient method with scaled vector transport for solving Problem \ref{pbl:main} \cite{absil2008optimization, ring2012optimization, sato2013new, sato2015dai}. \label{alg:conjugate}}
\begin{algorithmic}[1]
\REQUIRE A Riemann manifold $M$, a retraction $R$, a smooth function $f:M\rightarrow\mathbb{R}$, an initial point $x_0 \in M$, convergence tolerance $\epsilon > 0$.
\ENSURE Sequence $\{x_k\}_{k=0,1,\cdots} \subset M$.
\STATE Set $\eta_0 = -g_0 := -\grad f(x_0)$
\STATE $k \leftarrow 0.$
\WHILE{$\norm{g_k}_{x_k} > \epsilon$}
\STATE Determine the positive step size $\alpha_k > 0$ and set \begin{align}~\label{eq:update} x_{k+1}=R_{x_k}(\alpha_k\eta_k). \end{align}
\STATE Compute $g_{k+1} = -\grad f(x_{k+1})$.
\STATE Compute the parameter $\beta_{k+1}$ (e.g., Hestenes-Stiefel formula \eqref{eq:HS}).
\STATE Set the search direction
\begin{align}~\label{eq:direction} \eta_{k+1} = -g_{k+1} + \beta_{k+1}\mathcal{T}^S_{\alpha_k\eta_k}(\eta_k), \end{align}
where $\mathcal{T}^S$ is the scaled vector transport \eqref{eq:scaled} respect to $R$.
\STATE $k \leftarrow k + 1.$
\ENDWHILE
\end{algorithmic}
\end{algorithm}

In this paper, we say that the search direction $\eta_k \in T_{x_k}M$ is a descent direction if $\ip{g_k}{\eta_k} < 0$ holds. In addition, $\eta_k$ is a sufficient descent direction (see \cite{narushima2014survey}) if the sufficient descent condition,
\begin{align}~\label{eq:sufficient}
\ip{g_k}{\eta_k} \leq -\kappa\norm{g_k}^2,
\end{align}
holds for some constant $\kappa > 0$. In \eqref{eq:update}, for a given descent direction $\eta_k \in T_xM$ at $x \in M$, one often chooses a step size $\alpha_k > 0$ to satisfy the Armijo condition \cite[Definition 2.3]{hosseini2018line}, \cite[(1a)]{ring2012optimization}, namely,
\begin{align}~\label{eq:Armijo}
f(R_{x_k}(\alpha_k\eta_k)) \leq f(x_k)+c_1\alpha_k\ip{\grad f(x_k)}{\eta_k}_{x_k},
\end{align}
where $0 < c_1 < 1$. The following condition is called the curvature condition \cite[Definition 2.5]{hosseini2018line}:
\begin{align}
\label{eq:Wolfe}
\ip{\grad f(R_{x_k}(\alpha_k\eta_k))} {\mathcal{T}^R_{\alpha_k\eta_k}(\eta_k)}_{R_{x_k}(\alpha_k\eta_k)} \geq c_2 \ip{\grad f(x_k)}{\eta_k}_{x_k},
\end{align}
where $0<c_1<c_2<1$. Conditions \eqref{eq:Armijo} and \eqref{eq:Wolfe} are called the Wolfe conditions \cite[Definition 2.7]{hosseini2018line}, \cite[(1a), (1b)]{ring2012optimization}. If condition \eqref{eq:Wolfe} is replaced by
\begin{align}
\label{eq:sWolfe}
\abs{ \ip{\grad f(R_{x_k}(\alpha_k\eta_k))} {\mathcal{T}^R_{\alpha_k\eta_k}(\eta_k)}_{R_{x_k}(\alpha_k\eta_k)}} \leq c_2\abs{\ip{\grad f(x_k)}{\eta_k}_{x_k}},
\end{align}
then \eqref{eq:Armijo} and \eqref{eq:sWolfe} are called the \emph{strong} Wolfe conditions \cite[(1a), (2)]{ring2012optimization}.

In \eqref{eq:direction}, $\beta_{k+1}$ is given by generalizations of the formulas in Euclidean space (see \cite{hestenes1952methods, fletcher1964function, polak1969note, polyak1969conjugate, dai1999nonlinear}), e.g.,
\begin{align}
\label{eq:HS}
\beta_{k+1}^\mathrm{HS}&=\frac{\ip{g_{k+1}}{y_{k+1}}_{x_{k+1}}}{\ip{g_{k+1}}{\mathcal{T}^S_{\alpha_{k}\eta_{k}}(\eta_{k})}_{x_{k+1}}-\ip{g_{k}}{\eta_{k}}_{x_{k}}}, \\
\label{eq:FR}
\beta_{k+1}^\mathrm{FR}&=\frac{\norm{g_{k+1}}^2_{x_{k+1}}}{\norm{g_k}^2_{x_{k}}}, \\
\label{eq:PRP}
\beta_{k+1}^\mathrm{PRP}&=\frac{\ip{g_{k+1}}{y_{k+1}}_{x_{k+1}}}{\norm{g_k}^2_{x_{k}}}, \\
\label{eq:DY}
\beta_{k+1}^\mathrm{DY}&=\frac{\norm{g_{k+1}}^2_{x_{k+1}}}{\ip{g_{k+1}}{\mathcal{T}^S_{\alpha_{k}\eta_{k}}(\eta_{k})}_{x_{k+1}}-\ip{g_{k}}{\eta_{k}}_{x_{k}}},
\end{align}
where $y_{k+1} := g_{k+1} - \mathcal{T}^S_{\alpha_k\eta_k}(g_k)$. Formulas \eqref{eq:HS}, \eqref{eq:FR}, \eqref{eq:PRP}, and \eqref{eq:DY} are called the Hestenes-Stiefel (HS), Fletcher-Reeves (FR), Polak-Ribi\`ere-Polyak (PRP), and Dai-Yuan (DY) formulas, respectively. In \cite{sato2013new}, Sato and Iwai indicated that by using scaled vector transport, the FR method converges globally under the \emph{strong} Wolfe conditions \eqref{eq:Armijo} and \eqref{eq:sWolfe}. In \cite{sato2015dai}, Sato proved that the DY method converges globally under the Wolfe conditions \eqref{eq:Armijo} and \eqref{eq:Wolfe}. The HS and PRP methods have good numerical performance; however, no useful convergence analyses have been presented for them. To make up for these shortcomings, hybrid-type formulas, such as
\begin{align}
\label{eq:Hybrid1}
\beta_{k+1}^\mathrm{Hyb1} &= \max\{0,\min\{\beta_{k+1}^\mathrm{HS},\beta_{k+1}^\mathrm{DY}\}\}, \\
\label{eq:Hybrid2}
\beta_{k+1}^\mathrm{Hyb2} &= \max\{0,\min\{\beta_{k+1}^\mathrm{FR},\beta_{k+1}^\mathrm{PRP}\}\},
\end{align}
have been developed in Euclidean space (see \cite{hu1991global, dai2001efficient}). Below, we call the hybrid methods using \eqref{eq:Hybrid1} and \eqref{eq:Hybrid2}, Hybrid1 and Hybrid2, respectively. The Hybrid1 method was proposed by Dai and Yuan \cite{dai2001efficient}, and the Hybrid2 method was suggested by Hu and Storey \cite{hu1991global}. In \cite{sakai2020hybrid}, Sakai and Iiduka generalized the Hybrid1 method on Riemannian manifolds and proved that it converges globally under the strong Wolfe conditions \eqref{eq:Armijo} and \eqref{eq:sWolfe}. They also showed that the numerical performance of the Hybrid1 method is better than that of the PRP method \cite[Section 4]{sakai2020hybrid}. In the next section (Theorem \ref{thm:Hybrid2sufficient}), we generalize the Hybrid2 method on a Riemannian manifold and prove that it satisfies the sufficient descent condition under the strong Wolfe conditions. Moreover, we give its convergence analysis (Theorem \ref{thm:Hybrid2convergence}).

We consider the nonlinear conjugate gradient methods that can guarantee the sufficient descent condition \eqref{eq:sufficient} regardless of the type of line search used. We generalize the Hager-Zhang method \cite{hager2005new} to Riemannian manifolds, as
\begin{align}~\label{eq:HZ}
\beta_{k+1}^\mathrm{HZ} = \beta^\mathrm{HS}_{k+1} - \mu\frac{\norm{y_{k+1}}_{x_{k+1}}^2\ip{g_{k+1}}{\mathcal{T}^S_{\alpha_k\eta_k}(\eta_k)}_{x_{k+1}}}{\left(\ip{g_{k+1}}{\mathcal{T}^S_{\alpha_{k}\eta_{k}}(\eta_{k})}_{x_{k+1}}-\ip{g_{k}}{\eta_{k}}_{x_{k}}\right)^2},
\end{align}
where $y_{k+1} := g_{k+1} - \mathcal{T}^S_{\alpha_k\eta_k}(g_k)$ and $\mu > 1/4$. Moreover, we modify $\beta_{k+1}$ of the form $\beta_{k+1} = \ip{g_{k+1}}{\xi_{k+1}}_{x_{k+1}}$ by
\begin{align}~\label{eq:SD}
\beta_{k+1}^\mathrm{SD} = \beta_{k+1} - \mu \norm{\xi_{k+1}}_{x_{k+1}}^2\ip{g_{k+1}}{\mathcal{T}^S_{\alpha_{k}\eta_{k}}(\eta_{k})}_{x_{k+1}},
\end{align}
where $\xi_{k+1} \in T_{x_{k+1}}M$ is any tangent vector and $\mu > 1/4$ (see \cite{dai2010nonlinear, narushima2014survey}). For instance, if we set $\xi_{k+1} = y_{k+1} / \left(\ip{g_{k+1}}{\mathcal{T}^S_{\alpha_{k}\eta_{k}}(\eta_{k})}_{x_{k+1}}-\ip{g_{k}}{\eta_{k}}_{x_{k}}\right)$, we have $\beta_{k+1}^\mathrm{SD} = \beta_{k+1}^\mathrm{HZ}$. We will show that the SD method always satisfies the sufficient descent condition \eqref{eq:sufficient} with $\kappa = 1 - (1/4\mu)$ (Theorem \ref{thm:SDsufficient}).

\section{Sufficient Descent Properties of the Riemannian Conjugate Gradient Methods}~\label{sec:sufficient}
In this section, we recall the properties of the FR \eqref{eq:FR}, DY \eqref{eq:DY} and Hybrid1 \eqref{eq:Hybrid1} methods (see \cite{sato2013new, sato2015dai, sakai2020hybrid}).

\begin{proposition}~\label{prop:sufficient}
The following statements hold:
\begin{description}
\item[(P1)] If $\beta_{k+1} = \beta_{k+1}^\mathrm{FR}$ and $\alpha_k$ satisfies the strong Wolfe conditions \eqref{eq:Armijo} and \eqref{eq:sWolfe} with $0<c_1<c_2<1/2$, then
\begin{align*}
-\frac{1}{1-c_2}\norm{g_k}_{x_k}^2 \leq \ip{g_k}{\eta_k}_{x_k} \leq -\frac{1-2c_2}{1-c_2}\norm{g_k}_{x_k}^2,
\end{align*}
for all $k=0,1,\cdots$. Thus, the FR method satisfies the sufficient descent condition \eqref{eq:sufficient} with $\kappa = (1-2c_2)/(1-c_2) > 0$.
\item[(P2)] If $\beta_{k+1} = \beta_{k+1}^\mathrm{DY}$ and $\alpha_k$ satisfies the Wolfe conditions \eqref{eq:Armijo} and \eqref{eq:Wolfe}, then
\begin{align*}
-\frac{1}{1-c_2}\norm{g_k}_{x_k}^2 \leq \ip{g_k}{\eta_k}_{x_k} \leq -\frac{1}{1+c_2}\norm{g_k}_{x_k}^2,
\end{align*}
for all $k=0,1,\cdots$. Thus, the DY method satisfies the sufficient descent condition \eqref{eq:sufficient} with $\kappa = 1/(1+c_2) > 0$.
\item[(P3)] If $\beta_{k+1} = \beta_{k+1}^\mathrm{Hyb1}$ and $\alpha_k$ satisfies the strong Wolfe conditions \eqref{eq:Armijo} and \eqref{eq:sWolfe}, then
\begin{align*}
-\frac{1 + c_2}{1-c_2}\norm{g_k}_{x_k}^2 \leq \ip{g_k}{\eta_k}_{x_k} \leq -\frac{1 - c_2}{1+c_2}\norm{g_k}_{x_k}^2,
\end{align*}
for all $k=0,1,\cdots$. Thus, the Hybrid1 method satisfies the sufficient descent condition \eqref{eq:sufficient} with $\kappa = (1-c_2)/(1+c_2) > 0$.
\end{description}
\end{proposition}
Proposition \ref{prop:sufficient} implies that the FR, DY and Hybrid1 methods satisfy the sufficient descent condition \eqref{eq:sufficient}, dependent of the type of line search used. Here, (P1) is the result in \cite[Lemma 4.1]{sato2013new}, and (P2) and (P3) are easily shown from \cite[(35)]{sakai2020hybrid}.

\subsection{A Sufficient Descent Property of the Hybrid2 method}
In this section, we show that the Hybrid2 method generates a sufficient descent direction \eqref{eq:sufficient} at every iteration. This result is a simple extension of the result in Proposition \ref{prop:sufficient} (P1).

\begin{theorem}~\label{thm:Hybrid2sufficient}
Let $f:M \rightarrow \mathbb{R}$ be a smooth function. If each $\alpha_k>0$ satisfies the strong Wolfe conditions \eqref{eq:Armijo} and \eqref{eq:sWolfe}, with $0 < c_1 < c_2 < 1/2$, and $\beta_{k+1}$ satisfies\footnote{The formulas defined by \eqref{eq:FR} and \eqref{eq:Hybrid2} satisfy $\abs{\beta_{k+1}} \leq \beta_{k+1}^\mathrm{FR}$.} $\abs{\beta_{k+1}} \leq \beta_{k+1}^\mathrm{FR}$, then any sequence $\{x_k\}_{k=0,1\cdots}$ generated by Algorithm \ref{alg:conjugate} satisfies
\begin{align}~\label{eq:Hybrid2descent}
-\frac{1}{1-c_2}\norm{g_k}_{x_k}^2 \leq \ip{g_k}{\eta_k}_{x_k} \leq -\frac{1-2c_2}{1-c_2}\norm{g_k}_{x_k}^2,
\end{align}
for all $k=0,1, \cdots$.
\end{theorem}
\begin{proof}
The proof is by induction. If $k=0$, \eqref{eq:Hybrid2descent} clearly holds. Assume that \eqref{eq:Hybrid2descent} holds for some $k \geq 0$. By $c_2 < 1/2$, we obtain $\ip{g_k}{\eta_k}_{x_k} < 0$. From the search direction \eqref{eq:direction}, we have
\begin{align*}
\frac{\ip{g_{k+1}}{\eta_{k+1}}_{x_{k+1}}}{\norm{g_{k+1}}_{x_{k+1}}^2} = -1 + \beta_{k+1}\frac{\ip{g_{k+1}}{\mathcal{T}^S_{\alpha_k\eta_k}(\eta_k)}_{x_{k+1}}}{\norm{g_{k+1}}_{x_{k+1}}^2},
\end{align*}
which implies
\begin{align}~\label{eq:Hybrid2induction}
\frac{\ip{g_{k+1}}{\eta_{k+1}}_{x_{k+1}}}{\norm{g_{k+1}}_{x_{k+1}}^2} = -1 + \frac{\beta_{k+1}}{\beta^\mathrm{FR}_{k+1}}\frac{s_k\ip{g_{k+1}}{\mathcal{T}^R_{\alpha_k\eta_k}(\eta_k)}_{x_{k+1}}}{\norm{g_{k}}_{x_{k}}^2},
\end{align}
where
\begin{align*}
s_k:=\min\left\{1,\frac{\norm{\eta_k}_{x_k}}{\norm{\mathcal{T}^R_{\alpha_k\eta_k}(\eta_k)}_{x_{k+1}}}\right\} \in [0,1].
\end{align*}
From the second condition of the strong Wolfe conditions \eqref{eq:sWolfe} and $\ip{g_k}{\eta_k}_{x_k}<0$, we obtain
\begin{align*}
\abs{\beta_{k+1}\ip{g_{k+1}}{\mathcal{T}^R_{\alpha_k\eta_k}(\eta_k)}_{x_{k+1}}} \leq -c_2\abs{\beta_{k+1}}\ip{g_k}{\eta_k}_{x_k},
\end{align*}
which together with \eqref{eq:Hybrid2induction} implies
\begin{align*}
-1 + c_2s_k \frac{\abs{\beta_{k+1}}}{\beta_{k+1}^\mathrm{FR}}\frac{\ip{g_k}{\eta_k}_{x_k}}{\norm{g_k}_{x_k}^2}
\leq \frac{\ip{g_{k+1}}{\eta_{k+1}}_{x_{k+1}}}{\norm{g_{k+1}}_{x_{k+1}}^2}
\leq -1 - c_2s_k \frac{\abs{\beta_{k+1}}}{\beta_{k+1}^\mathrm{FR}}\frac{\ip{g_k}{\eta_k}_{x_k}}{\norm{g_k}_{x_k}^2}.
\end{align*}
From the left-hand side of the induction hypothesis \eqref{eq:Hybrid2descent}, we have
\begin{align*}
-1 - c_2s_k \frac{\abs{\beta_{k+1}}}{\beta_{k+1}^\mathrm{FR}}\frac{1}{1-c_2}
\leq \frac{\ip{g_{k+1}}{\eta_{k+1}}_{x_{k+1}}}{\norm{g_{k+1}}_{x_{k+1}}^2}
\leq -1 + c_2s_k \frac{\abs{\beta_{k+1}}}{\beta_{k+1}^\mathrm{FR}}\frac{1}{1-c_2}.
\end{align*}
Utilizing the assumption $\abs{\beta_{k+1}} \leq \beta_{k+1}^\mathrm{FR}$ and $0 \leq s_k \leq 1$, we obtain
\begin{align*}
-1 - \frac{c_2}{1-c_2}
\leq \frac{\ip{g_{k+1}}{\eta_{k+1}}_{x_{k+1}}}{\norm{g_{k+1}}_{x_{k+1}}^2}
\leq -1 + \frac{c_2}{1-c_2}.
\end{align*}
This implies that \eqref{eq:Hybrid2descent} holds for $k+1$. \qed
\end{proof}

Moreover, we prove the global convergence of the Hybrid2 method under the strong Wolfe conditions and the following assumption.
\begin{assumption}~\label{asm:main}
Let $M$ be a Riemannian manifold and R be a retraction on M. Let $f:M\rightarrow\mathbb{R}$ be a smooth, bounded below function. Then, we assume that there exists $L>0$ such that
\begin{align*}
|\mathrm{D}(f \circ R_x)(t\eta)[\eta]-\mathrm{D}(f \circ R_x)(0_x)[\eta]| \leq Lt,
\end{align*}
where $x \in M$, $\eta \in T_xM$, $\norm{\eta}_x=1$ and $t \geq 0$.
\end{assumption}
This is the assumption for Zoutendijk's theorem (Theorem \ref{thm:Zoutendijk}) on Riemannian manifolds. Zoutendijk's theorem is described on Riemannian manifolds as follows:
\begin{theorem}[Zoutendijk \cite{ring2012optimization}]~\label{thm:Zoutendijk}
Let $(M,g)$ be a Riemannian manifold and R be a retraction on M. Suppose $f:M\rightarrow\mathbb{R}$ satisfies Assumption \ref{asm:main}. Suppose further that in Algorithm \ref{alg:conjugate}, each step size $\alpha_k>0$ satisfies the strong Wolfe conditions \eqref{eq:Armijo} and \eqref{eq:sWolfe}. Then the following series converges:
\begin{align}
\label{eq:Zoutendijk}
\sum_{k=0}^\infty \frac{\ip{g_k}{\eta_k}_{x_k}^2}{\norm{\eta_k}^2_{x_k}} < \infty .
\end{align}
\end{theorem}

The proof of this theorem is along the lines of Zoutendijk's theorem in Euclidean space (see \cite[Theorem 3.3]{ring2012optimization}). Global convergence proofs for Riemannian conjugate gradient methods are often based on Zoutendijk's theorem. Theorem \ref{thm:Hybrid2convergence} guarantees global convergence of the Hybrid2 method \eqref{eq:Hybrid2}. It is a generalization of the convergence theorem of the Hybrid2 method in Euclidean space \cite{gilbert1992global}.

\begin{theorem}~\label{thm:Hybrid2convergence}
Let $f:M \rightarrow \mathbb{R}$ be a function satisfying Assumption \ref{asm:main}. If each $\alpha_k>0$ satisfies the strong Wolfe conditions \eqref{eq:Armijo} and \eqref{eq:sWolfe}, with $0 < c_1 < c_2 < 1/2$, and $\beta_{k+1}$ satisfies $\abs{\beta_{k+1}} \leq \beta_{k+1}^\mathrm{FR}$, then any sequence $\{x_k\}_{k=0,1\cdots}$ generated by Algorithm \ref{alg:conjugate} satisfies
\begin{align}~\label{eq:Hybrid2converge}
\liminf_{k \to \infty}\norm{g_k}_{x_k}=0.
\end{align}
\end{theorem}
\begin{proof}
We prove \eqref{eq:Hybrid2converge} by contradiction. If $g_{k_0} = 0$ for some $k_0$, then \eqref{eq:Hybrid2converge} follows. Assume that
\begin{align*}
\liminf_{k \to \infty}\norm{g_k}_{x_k} > 0.
\end{align*}
Then, noting $\norm{g_k}_{x_k}\neq 0$ for all $k$, there exists $\gamma > 0$ such that
\begin{align*}
\norm{g_k}_{x_k} \geq \gamma > 0,
\end{align*}
for all $k$. From \eqref{eq:sWolfe} and \eqref{eq:Hybrid2descent}, we have
\begin{align*}
\abs{\ip{g_k}{\mathcal{T}^R_{\alpha_{k-1}\eta_{k-1}}(\eta_{k-1})}_{x_k}} &\leq -c_2\ip{g_{k-1}}{\eta_{k-1}}_{x_{k-1}} \\
&\leq \frac{c_2}{1-c_2}\norm{g_{k-1}}_{x_{k-1}}^2.
\end{align*}
Thus, from \eqref{eq:direction} and \eqref{eq:Hybrid2descent}, and using the condition $\abs{\beta_k} \leq \beta^\mathrm{FR}_k \leq \norm{g_k}_{x_k}^2 / \norm{g_{k-1}}_{x_{k-1}}^2$, we have
\begin{align*}
\norm{\eta_k}_{x_k}^2 &\leq \norm{g_k}_{x_k}^2 + 2s_k\abs{\beta_k\ip{g_k}{\mathcal{T}^R_{\alpha_{k-1}\eta_{k-1}}(\eta_{k-1})}_{x_k}} + \norm{\beta_k\mathcal{T}^S_{\alpha_{k-1}\eta_{k-1}}(\eta_{k-1})}^2_{x_k} \\
&\leq \norm{g_k}_{x_k}^2 + \frac{2c_2}{1-c_2}\abs{\beta_k}\norm{g_{k-1}}^2_{x_{k-1}} + \beta_k^2\norm{\eta_{k-1}}_{x_{k-1}}^2 \\
&\leq \hat{c}\norm{g_k}_{x_k}^2 + \beta_k^2\norm{\eta_{k-1}}_{x_{k-1}}^2,
\end{align*}
where $\hat{c} := (1 + c_2) / (1 - c_2) > 1$. Applying this equation repeatedly, we obtain
\begin{align*}
\norm{\eta_k}_{x_k}^2 &\leq \hat{c}\norm{g_k}_{x_k}^2 + \beta_k^2\left(\hat{c}\norm{g_{k-1}}_{x_{k-1}}^2 + \beta_{k-1}^2\norm{\eta_{k-2}}_{x_{k-2}}^2\right) \\
&\leq \hat{c}\left(\norm{g_k}_{x_k}^2 + \beta_k^2\norm{g_{k-1}}_{x_{k-1}}^2 + \cdots + \beta_k^2\beta_{k-1}^2\cdots\beta_2^2\norm{g_1}_{x_1}^2\right) + \beta_k^2\beta_{k-1}^2\cdots\beta_1^2\norm{\eta_0}_{x_0}^2 \\
&\leq \hat{c}\norm{g_k}_{x_k}^4\left(\frac{1}{\norm{g_k}_{x_k}^2} + \frac{1}{\norm{g_{k-1}}_{x_{k-1}}^2} + \cdots + \frac{1}{\norm{g_1}_{x_1}^2}\right) + \frac{\norm{g_k}_{x_k}^4}{\norm{g_0}_{x_0}^2} \\
&< \hat{c}\norm{g_k}_{x_k}^4\sum_{j=0}^k\frac{1}{\norm{g_j}_{x_j}^2} \\
&\leq \frac{\hat{c}}{\gamma^2}\norm{g_k}_{x_k}^4(k+1).
\end{align*}
This implies that
\begin{align*}
\frac{\norm{g_k}_{x_k}^4}{\norm{\eta_k}_{x_k}^2} \geq \frac{\gamma^2}{\hat{c}(k+1)},
\end{align*}
which together with \eqref{eq:Hybrid2descent}, gives
\begin{align*}
\sum_{k=0}^\infty \frac{\ip{g_k}{\eta_k}_{x_k}^2}{\norm{\eta_k}^2_{x_k}} &= \sum_{k=0}^\infty \frac{\norm{g_k}_{x_k}^4}{\norm{\eta_k}_{x_k}^2}\frac{\ip{g_k}{\eta_k}_{x_k}^2}{\norm{g_k}_{x_k}^4} \\
&\geq \left( \frac{2c_2 - 1}{1 -c_1} \right)^2 \sum_{k=0}^\infty \frac{\gamma^2}{\hat{c}(k+1)} \\
& = \infty.
\end{align*}
This contradicts \eqref{eq:Zoutendijk} in Zoutendijk's theorem (Theorem \ref{thm:Zoutendijk}) and completes the proof. \qed
\end{proof}

\subsection{Sufficient Descent Property of the SD method}
Theorem \ref{thm:SDsufficient} asserts that the SD method \eqref{eq:SD} produces sufficient descent directions \eqref{eq:sufficient} regardless of the choice of the step size $\alpha_k$.
\begin{theorem}~\label{thm:SDsufficient}
Let $f:M \rightarrow \mathbb{R}$ be a smooth function. If $\beta_{k+1} = \beta_{k+1}^\mathrm{SD}$, then any sequence $\{x_k\}_{k=0,1\cdots}$ generated by Algorithm \ref{alg:conjugate} satisfies
\begin{align}~\label{eq:SDsufficient}
\ip{g_k}{\eta_k}_{x_k} \leq -\left(1-\frac{1}{4\mu}\right)\norm{g_{k}}_{x_k}^2.
\end{align}
\end{theorem}
\begin{proof}
From \eqref{eq:SD}, we obtain
\begin{align}~\label{eq:HZtemp}
\begin{split}
\ip{g_k}{\eta_k}_{x_k} = &-\norm{g_k}_{x_k}^2 + \ip{g_k}{\xi_k}_{x_k}\ip{g_k}{\mathcal{T}^S_{\alpha_{k-1}\eta_{k-1}}(\eta_{k-1})}_{x_k} \\
& \quad - \mu\norm{\xi_k}_{x_k}^2\ip{g_k}{\mathcal{T}^S_{\alpha_{k-1}\eta_{k-1}}(\eta_{k-1})}_{x_k}^2.
\end{split}
\end{align}
An upper bound for the middle term in \eqref{eq:HZtemp} is obtained using the inequality,
\begin{align*}
\ip{u_k}{v_k}_{x_k} \leq \frac{\norm{u_k}_{x_k}^2 + \norm{v_k}_{x_k}^2}{2}
\end{align*}
with the choice
\begin{align*}
u_k := \frac{1}{\sqrt{2\mu}}g_k \quad and \quad v_k := \sqrt{2\mu}\ip{g_k}{\mathcal{T}^S_{\alpha_{k-1}\eta_{k-1}}(\eta_{k-1})}\xi_k.
\end{align*}
Then, we have
\begin{align*}
&\ip{g_k}{\xi_k}_{x_k}\ip{g_k}{\mathcal{T}^S_{\alpha_{k-1}\eta_{k-1}}(\eta_{k-1})}_{x_k} \\
&\qquad \leq \frac{1}{4\mu}\norm{g_k}_{x_k}^2 + \mu\norm{\xi_k}_{x_k}^2\ip{g_k}{\mathcal{T}^S_{\alpha_{k-1}\eta_{k-1}}(\eta_{k-1})}_{x_k}^2.
\end{align*}
Combining this with \eqref{eq:HZtemp}, we obtain \eqref{eq:SDsufficient}. \qed
\end{proof}

\section{Line Search Algorithm on Riemannian Manifolds}~\label{sec:LineSearch}
In Algorithm \ref{alg:conjugate}, we need to use a line search algorithm to determine the step size $\alpha_k$.
A backtracking line search algorithm (see \cite[Chapter 3, Algorithm 3.1]{nocedal2006numerical}) is widely used in optimization algorithms in Euclidean space to find a step size that satisfies the Armijo condition \eqref{eq:Armijo}. Algorithm \ref{alg:BacktrackArmijo} is a backtracking line search on Riemannian manifold. This algorithm multiplies a positive constant $\rho > 0$ until a step size $\alpha$ satisfying the Armijo condition is found.

\begin{algorithm}
\caption{Backtracking line search on Riemannian manifold $M$. \label{alg:BacktrackArmijo}}
\begin{algorithmic}[1]
\REQUIRE A smooth function $f:M\rightarrow\mathbb{R}$, a point $x \in M$, a descent direction $\eta \in T_xM$, scalars $0 < \alpha_{\mathrm{hi}}$, $\rho \in (0,1)$.
\ENSURE A positive step size $\alpha > 0$ satisfying the Armijo condition \eqref{eq:Armijo}.
\STATE $\alpha \leftarrow \alpha_{\mathrm{hi}}$
\WHILE{$f(R_{x}(\alpha\eta)) > f(x) + c_1\alpha \ip{\grad f(x)}{\eta}_{x}$}
\STATE $\alpha \leftarrow \rho\alpha$
\ENDWHILE
\RETURN $\alpha$
\end{algorithmic}
\end{algorithm}

However, a backtracking line search algorithm cannot be used for the Wolfe or the strong Wolf conditions. To find a step size satisfying the \emph{strong} Wolfe conditions, we present Algorithm \ref{alg:LineSearch}, a generalization of the algorithm in \cite[Chapter 3, Algorithm 3.5]{nocedal2006numerical} for strong Wolfe conditions in Euclidean space. Algorithm \ref{alg:LineSearch} calls the \emph{zoom} function (Algorithm \ref{alg:Zoom}), which successively decreases the size of the interval until an acceptable step size is found (see \cite[Chapter 3, Algorithm 3.6]{nocedal2006numerical}). The parameter $\alpha_\mathrm{hi}$ is a user-supplied bound on the maximum step size. Algorithm \ref{alg:LineSearch} returns a positive step size, $\alpha_\star > 0$, that satisfies the strong Wolfe conditions.

\begin{algorithm}
\caption{Line search algorithm on Riemannian manifold $M$. \label{alg:LineSearch}}
\begin{algorithmic}[1]
\REQUIRE A smooth function $f:M\rightarrow\mathbb{R}$, a point $x \in M$, a descent direction $\eta \in T_xM$, scalars $0<c_1<c_2<1$, $0 < \alpha_{\mathrm{hi}}$ and $\alpha_0 \in (0,\alpha_\mathrm{hi})$.
\ENSURE A positive step size $\alpha > 0$ satisfying the \emph{strong} Wolfe conditions \eqref{eq:Armijo} and \eqref{eq:sWolfe}.
\STATE Set $\phi(\alpha) = f(R_x(\alpha\eta))$.
\STATE $i \leftarrow 0$.
\LOOP
\IF{$\phi(\alpha_i) > \alpha(0) + \alpha c_1\phi^\prime(0)$ \OR $[\phi(\alpha_i) \geq \phi(\alpha_{i-1})$ \AND $i \geq 1]$ }
\STATE Set $\alpha_\star = \mathrm{Zoom}(\alpha_{i-1}, \alpha_i)$ and stop.
\ELSIF{$\abs{\phi^\prime(\alpha_i)} \leq -c_2\phi^\prime(0)$}
\STATE Set $\alpha_\star = \alpha_i$ and stop.
\ELSIF{$\phi^\prime(0) \geq 0$}
\STATE Set $\alpha_\star = \mathrm{Zoom}(\alpha_i, \alpha_{i-1})$ and stop.
\ENDIF
\STATE Choose $\alpha_{i+1} \in (\alpha_i, \alpha_\mathrm{hi})$.
\STATE $i \leftarrow i+1$.
\ENDLOOP
\RETURN $\alpha_\star$
\end{algorithmic}
\end{algorithm}

\begin{algorithm}
\caption{Zoom. \label{alg:Zoom} \cite[Chapter 3, Algorithm 3.6]{nocedal2006numerical}}
\begin{algorithmic}[1]
\REQUIRE Scalars $\alpha_\mathrm{min}, \alpha_\mathrm{max} > 0$, and $\phi(\alpha) = f(R_x(\alpha\eta))$.
\ENSURE $\alpha = \mathrm{Zoom}(\alpha_\mathrm{min}, \alpha_\mathrm{max})$.
\LOOP
\STATE Interpolate (using quadratic, cubic, or bisection) to find a trial step length $\alpha_j \in (\alpha_\mathrm{min}, \alpha_\mathrm{max})$.
\IF{$\phi(\alpha_j) > \phi(0) + c_1\alpha_j\phi^\prime(0)$ \OR $\phi(\alpha_j) \geq \phi(\alpha_\mathrm{min})$}
\STATE $\alpha_\mathrm{max} \leftarrow \alpha_j$
\ELSE
\IF{$\abs{\phi^\prime(\alpha_j)} \leq -c_2\phi^\prime(0)$}
\STATE Set $\alpha_\star = \alpha_j$ and stop.
\ELSIF{$\phi^\prime(\alpha_j) (\alpha_\mathrm{max} - \alpha_\mathrm{min}) \geq 0$}
\STATE $\alpha_\mathrm{max} \leftarrow \alpha_\mathrm{min}$.
\ENDIF
\STATE $\alpha_\mathrm{min} \leftarrow \alpha_j$.
\ENDIF
\ENDLOOP
\RETURN $\alpha_\star$
\end{algorithmic}
\end{algorithm}

\section{Numerical Experiments}~\label{sec:NE}
Our experiments used the source code based on \texttt{pymanopt}\footnote{\url{https://www.pymanopt.org/}} (see \cite{townsend2016pymanopt}). In addition, Algorithm \ref{alg:LineSearch} was based on an implementation by SciPy\footnote{\url{https://docs.scipy.org/doc/scipy/reference/}} in Euclidean space. We solved four kinds of Riemannian optimization problems (Problem \ref{pbl:rayleigh}--\ref{pbl:off-diag}) on several Riemannian manifolds and objective functions.

Problem \ref{pbl:rayleigh} is the Rayleigh-quotient minimization problem on the unit sphere (see \cite[Chapter 4.6]{absil2008optimization}).
\begin{problem}~\label{pbl:rayleigh}
For $A \in \mathcal{S}^n_{++}$,
\begin{align*}
\textrm{minimize}& \quad f(x)=x^\top Ax, \\
\textrm{subject to}& \quad x \in \mathbb{S}^{n-1} := \{ x \in \mathbb{R}^{n} : \norm{x} = 1\},
\end{align*}
where $\norm{\cdot}$ denotes the Euclidean norm and $\mathcal{S}^n_{++}$ denotes the set of all $n \times n$ symmetric positive-definite matrices.
\end{problem}
In the experiments, we set $n = 100$ and generated a matrix $A \in \mathcal{S}^n_{++}$ with randomly chosen elements by using \texttt{sklearn.datasets.make\_spd\_matrix}.

Problem \ref{pbl:brockett} is the Brockett-cost-function minimization problem on a Stiefel manifold (see \cite[Chapter 4.8]{absil2008optimization}).
\begin{problem}~\label{pbl:brockett}
For $A \in \mathcal{S}^n_{++}$ and $N = \mathrm{diag}(\mu_0, \cdots ,\mu_p)$ $(0 \leq \mu_0\leq \cdots \leq\mu_p)$,
\begin{align*}
\textrm{minimize}& \quad f(X)=\mathrm{tr}(X^\top AXN)\\
\textrm{subject to}& \quad X \in \mathrm{St}(p,n) := \{X \in \mathbb{R}^{n \times p} : X^\top X=I_p\}.
\end{align*}
\end{problem}
In the experiments, we set $p=5$, $n=20$ and $N := \mathrm{diag}(1, 2, 3, 4, 5)$ and generated a matrix $A \in \mathcal{S}^n_{++}$ with randomly chosen elements by using \texttt{sklearn.datasets.make\_spd\_matrix}.

In \cite{vandereycken2013low}, Vandereycken discussed the following robust matrix completion problem (Problem \ref{pbl:completion}).
\begin{problem}~\label{pbl:completion}
For $A \in \mathbb{R}^{m \times n}$, and a subset $\Omega$ of the complete set of entries $\{1, \cdots ,m\}\times\{1, \cdots ,n\}$,
\begin{align*}
\textrm{minimize}& \quad f(X)=\norm{P_\Omega(X - A)}^2_F, \\
\textrm{subject to}& \quad X \in M_k := \{ X \in \mathbb{R}^{m \times n} : \mathrm{rank}(X) = k\},
\end{align*}
where $\norm{\cdot}_F$ denotes the Frobenius norm and
\begin{align*}
P_\Omega : \mathbb{R}^{m \times n} \rightarrow \mathbb{R}^{m \times n}, X_{ij} \mapsto
\begin{cases}
X_{ij} & (i,j) \in \Omega \\
0 & (i,j) \not\in \Omega
\end{cases}.
\end{align*}
\end{problem}
In the experiments, we set $m=n=100$ and $k=4$, and $\Omega$ contained each pair $(i,j) \in \{1, \cdots ,m\}\times\{1, \cdots ,n\}$ with probability $1/2$. Moreover, we used a matrix $A \in \mathbb{R}^{m \times n}$ that was generated with randomly chosen elements by using \texttt{numpy.random.randn}.

In \cite{absil2006joint}, Absil and Gallivan introduced the following off-diagonal cost function minimization problem on oblique manifolds (Problem \ref{pbl:off-diag}).
\begin{problem}~\label{pbl:off-diag}
For $C_i \in \mathcal{S}^{n}$ $(i=1, \cdots ,N)$,
\begin{align*}
\textrm{minimize}& \quad f(X)=\sum_{i=1}^N\norm{X^\top C_iX - \mathrm{ddiag}(X^\top C_iX)}^2_F \\
\textrm{subject to}& \quad X \in \mathcal{OB}(n,p):=\{X \in \mathbb{R}^{n \times p} : \mathrm{ddiag}(X^\mathrm{T}X)=I_p\},
\end{align*}
where $\mathcal{S}^n$ denotes the set of all $n \times n$ symmetric matrices and $\mathrm{ddiag}(X)$ denotes a diagonal matrix whose diagonal elements are those of $X$.
\end{problem}
In the experiments, we set $N=10$, $n=100$ and $p=5$ and generated ten matrices $B_i \in \mathbb{R}^{n \times n}$ $(i=1,2, \cdots ,10)$ with randomly chosen elements by using \texttt{numpy.random.randn}. Then, we set symmetric matrices $C_i \in \mathcal{S}^n$ as $C_i := (B_i + B_i^\top) / 2$ $(i=1,2, \cdots ,10)$.

The experiments used a MacBook Air (2017) with a 1.8 GHz Intel Core i5, 8 GB 1600 MHz DDR3 memory, and the macOS Mojave version 10.14.5 operating system. The algorithms were written in Python 3.7.6 with the NumPy 1.19.0 package and the Matplotlib 3.2.2 package. We solved the above four problems 100 times with each algorithm, that is, 400 times in total. If the stopping condition,
\begin{align*}
\norm{{\grad}f(x_k)}_{x_k}<10^{-6}
\end{align*}
was satisfied, we determined that a sequence had converged to an optimal solution. We compared seven Riemannian conjugate gradient methods, i.e., FR, DY, PRP, HS, HZ, Hybrid1, and Hybrid2 methods, and two line search algorithms, i.e., Algorithm \ref{alg:BacktrackArmijo} and \ref{alg:LineSearch}. In the HZ method, we set $\mu = 2$. In the Armijo condition \eqref{eq:Armijo} and the second condition of the strong Wolfe conditions \eqref{eq:sWolfe}, we set $c_1 = 10^{-4}$ and $c_2 = 0.9$. In Algorithm \ref{alg:BacktrackArmijo}, we set the scalars as $\alpha_\mathrm{hi} = 1$ and $\rho = 0.5$. In Algorithm \ref{alg:LineSearch}, we set the scalar as $\alpha_0 = 1$, and in step 11, we set $\alpha_i = 2\alpha_{i-1}$ (see \texttt{scipy.optimize.line\_search}).

For comparison, we calculated the performance profile \cite{dolan2002benchmarking}. The performance profile $P_s : \mathbb{R} \rightarrow [0, 1]$ is defined as follows: let $\mathcal{P}$ and $\mathcal{S}$ be the set of problems and solvers, respectively. For each $p \in \mathcal{P}$ and $s \in \mathcal{S}$, we define
\begin{align*}
t_{p,s} := (\text{iterations or time required to solve problem }p\text{ by solver }s).
\end{align*}
Furthermore, we defined the performance ratio $r_{p,s}$ as
\begin{align*}
r_{p,s} := \frac{t_{p,s}}{\min_{s^\prime \in \mathcal{S}}t_{p,s^\prime}}
\end{align*}
and defined the performance profile, for all $\tau \in \mathbb{R}$, as
\begin{align*}
P_s(\tau) := \frac{\#\{p \in \mathcal{P} : r_{p,s} \leq \tau \}}{\#\mathcal{P}},
\end{align*}
where $\# S$ denotes the number of elements of a set $S$.

Figure \ref{fig1} plots the performance profiles of each algorithm by using Algorithm \ref{alg:BacktrackArmijo} to determine the step size. In particular, Figure \ref{fig1} (a) and (b) plot the performance profiles versus the number of iterations and the elapsed time, respectively. They show that the HZ method solved the most problems, which is about the same number as the Hybrid1 method solved. The HZ and Hybrid1 methods substantially outperformed the other methods. In particular, Hybrid1 boasted outstanding performance in terms of the number of solved problems, number of iterations, and elapsed time. It can also be seen that Hybrid2 is not compatible with Algorithm \ref{alg:BacktrackArmijo}.

\begin{figure}[htbp]
\centering
\subfigure[iteration]{\includegraphics[width=50mm]{Figure1.eps}}
\subfigure[elapsed time]{\includegraphics[width=50mm]{Figure2.eps}}
\caption{Performance profiles of each algorithm versus the number of iterations (a) and the elapsed time (b) by using Algorithm \ref{alg:BacktrackArmijo} to determine the step size. \label{fig1}}
\end{figure}

Figure \ref{fig2} plots the performance profiles of each algorithm by using Algorithm \ref{alg:LineSearch} to determine the step size. In particular, Figure \ref{fig2} (a) and (b) plot the performance profiles versus the number of iterations and the elapsed time, respectively. It shows that Hybrid1 solved the most problems, and Hybrid2 solved the second-most problems. Unlike the case of using Algorithm \ref{alg:BacktrackArmijo}, Hybrid2 performed well when using Algorithm \ref{alg:LineSearch}. It can be seen that the PRP and HS methods have about the same performance. Similarly, the FR and DY methods have about the same performance.

\begin{figure}[htbp]
\centering
\subfigure[iteration]{\includegraphics[width=50mm]{Figure3.eps}}
\subfigure[elapsed time]{\includegraphics[width=50mm]{Figure4.eps}}
\caption{Performance profiles of each algorithm versus the number of iterations (a) and the elapsed time (b) by using Algorithm \ref{alg:LineSearch} to determine the step size. \label{fig2}}
\end{figure}

Figure \ref{fig3} plots the performance profiles of HZ, Hybrid1 and Hybrid2 by using Algorithm \ref{alg:BacktrackArmijo} and \ref{alg:LineSearch} to determine the step size. In particular, Figure \ref{fig3} (a) and (b) plots the performance profile versus the number of iterations and elapsed time, respectively. Figure \ref{fig3} (a) shows that when Algorithm \ref{alg:LineSearch} is used, all methods solve the problem in fewer iterations than in the case of using Algorithm \ref{alg:BacktrackArmijo}. It can be seen from Figure \ref{fig3} (b) that Algorithm \ref{alg:LineSearch} often takes a long time to execute.

\begin{figure}[htbp]
\centering
\subfigure[iteration]{\includegraphics[width=50mm]{Figure5.eps}}
\subfigure[elapsed time]{\includegraphics[width=50mm]{Figure6.eps}}
\caption{Performance profiles of the HZ, Hybrid1 and Hybrid2 methods versus the number of iterations (a) and the elapsed time (b) by using Algorithm \ref{alg:BacktrackArmijo} and \ref{alg:LineSearch} to determine the step size. \label{fig3}}
\end{figure}

It can be seen from Figure \ref{fig1}--\ref{fig3} that Hybrid1 performed the best in all cases, but the HZ method also performed well. It can also be seen that the performances of the Riemannian conjugate gradient methods depend greatly the type of line search used. In particular, Hybrid2 may perform poorly when certain line searches are used.

\section{Conclusion}~\label{sec:conclusion}
We generalized two nonlinear conjugate gradient methods, i.e., the HZ and Hybrid2 method. We proved that the Hybrid2 method \eqref{eq:Hybrid2} satisfies the sufficient descent condition and converges globally under the strong Wolfe conditions. In addition, we proved that the HZ method \eqref{eq:HZ} satisfies the sufficient descent condition regardless of the type of line search used. We also generalized two kinds of line search algorithms, i.e., Algorithm \ref{alg:BacktrackArmijo} and \ref{alg:LineSearch}. In numerical experiments, we showed that the HZ and Hybrid1 methods perform well. Moreover, we showed that the performance of the Riemannian conjugate gradient method depends on the type of line search used. In particular, if we use Algorithm \ref{alg:BacktrackArmijo} to determine the step size, Hybrid2 performs poorly. On the other hand, it performs better with a step size computed by Algorithm \ref{alg:LineSearch}, as the convergence analysis guarantees. Meanwhile, the numerical results showed that the HZ method converges quickly without depending on line search conditions. Hence, the HZ method is good for solving Riemannian optimization problems from the viewpoints of both theory and practice.

\section*{acknowledgements}
This work was supported by a JSPS KAKENHI Grant, Number JP18K11184.

\bibliographystyle{abbrv}