1. Introduction

Canonical correlation analysis plays an important role for assessing the relationship between sets of variables. In order to find the linear combination of the set of variables, the maximal correlation problem (MCP) is proposed by Hotelling ^[1,2]. When we consider the optimal solution to the MCP, which then corresponds to the global maximizer of the following equality constrained optimization problem:

$\begin{equation} \label{eq1} \begin{cases} \text{maximize}&r(x): = x^T Ax, \\ \text{subject to}&\left\| {x_i } \right\|_2 = 1, ~i = 1, 2, \cdots, m, \\ \end{cases} \end{equation}$

(1.1)

where $x = \left[{x_1^T, \cdots, x_m^T } \right] \in \mathbb{R}^n$, $x_i \in \mathbb{R}^{n_i}$, $A_{ij} \in \mathbb{R}^{n_i\times n_j}$ and

$A = \left( {\begin{array}{*{20}c} {A_{11} }&{A_{12} }& \cdots &{A_{1m} } \\ {A_{21} }&{A_{22} }& \cdots &{A_{2m} } \\ \vdots & \vdots & \ddots & \vdots \\ {A_{m1} }&{A_{m2} }& \cdots &{A_{mm} } \\ \end{array}} \right) \in \mathbb{R}^{n \times n}$

is a symmetric and positive definite matrix.

Using the method of Lagrange multipliers, we will get the following equations:

$\begin{equation} \label{eq2} \begin{cases} Ax = \Lambda x, \\ \left\| {x_i } \right\|_2 = 1, ~i = 1, 2, \cdots, m. \\ \end{cases} \end{equation}$

(1.2)

where

$\Lambda : = \text{diag}\left\{ {\begin{array}{*{20}c} {\lambda _1 I^{\left[{n_1 } \right]}, }&{\lambda _2 I^{\left[{n_2 } \right]}, }&{ \cdots, }&{\lambda _m I^{\left[{n_m } \right]} } \\ \end{array}} \right\}$

with $\lambda _1, \ldots, \lambda _m$ as the Lagrange multipliers and $I^{\left[{n_i } \right]} \in \mathbb{R}^{n_i \times n_i}$ denoting the $n_i \times n_i$ identity matrix. The multivariate eigenvalue problem (MEP) is cast exactly as the system of equations (1.2), which serves only as a necessary condition for the global maximum of the MCP ^[4]. $\lambda _1, \ldots, \lambda _m$ are usually called the multivariate eigenvalues.

The Horst-Jacobi algorithm as a generalization of the so called power method is introduced to solve MEP in ^[5], its convergence theory was established much later in ^[3]. The Gauss-Seidel algorithm has also been suggested in ^[3] and its monotone convergence is recently established in ^[7]. Based on a core engine in seeking global maximum of the MCP, the authors proposed an alternating variable method (AVM) in ^[9]. This algorithm is proved to enjoy the global and monotone convergence and it is also shown that for a nonnegative irreducible matrix $A$, the algorithm converges to the global maximizer from any nonnegative starting point ^[8,9].

In this paper, we propose the SOR-like AVM and the monotone convergence of the SOR-like AVM is proved.

---

2. Main results

If we express the objective function $r_i(x)$ as

$r_i(x) = \left( {x_i^T A_{ii} x_i + 2x_i^T \sum\limits_{j \ne i} {A_{ij} x_j } } \right) + \sum\limits_{k \ne i} {\sum\limits_{j \ne i} {x_k^T A_{kj} x_j } },$

It is very clear that if $x^* = \left[{(x_1^* )^T, \cdots, (x_m^* )^T } \right]^T \in \mathbb{R}^n$, is a global maximizer of the MCP (1.1), $x_i^* \in \mathbb{R}^{n_i}$ must be a global maximizer for the following subproblem

$\begin{equation} \label{eq4} \mathop {\max }\limits_{\left\| {x_i } \right\|_2 = 1} r_i(x_1^*, \cdots, x_{i - 1}^*, x_i, x_{i + 1}^*, \cdots, x_m^* ). \end{equation}$

(2.1)

Then we can get the following AVM algorithm, which have shown its efficient in ^[9].

Algorithm 1 (The framework of the alternating variable method (AVM)).

Select $x^{(0)} \in \mathcal {M}$, and set $k: = 0$.

while the stoping criterion is not met do

for $i = 1, \cdots, m$ do

$x_i^{(k + 1)} : = \arg \max _{\left\| {x_i } \right\|_2 = 1} r_i(x_1^{(k + 1)}, \cdots, x_{i - 1}^{(k + 1)}, x_i^{(k)}, x_{i + 1}^{(k)}, \cdots, x_m^{(k)} )$

end for

$k: = k+1$

end while

Denote

$\begin{align} x^{(k+1)}: & = \left[{(x_1^{(k + 1)} )^T, \cdots, (x_{i-1}^{(k + 1)} )^T, (x_i^{(k+1)} )^T, \cdots, (x_m^{(k+1)} )^T } \right]^T \in \mathbb{R}^n, \notag \\ g_i^{(k)} : & = \sum\limits_{j = 1}^{i - 1} {A_{ij} } x_j^{(k + 1)} + \sum\limits_{j = i + 1}^m {A_{ij} } x_j^{(k)} \in \mathbb{R}^{n_i }, \notag \\ m^{(k)} (x_i ): & = r_i\left( {x_1^{(k + 1)}, \cdots x_{i - 1}^{(k + 1)}, x_i^{(k)}, x_{i + 1}^{(k)}, \cdots, x_m^{(k)} } \right) \in \mathbb{R}, \notag \\ c_i^{(k)} : & = m^{(k)} (x_i ) - (x_i^{(k)})^T A_{ii} x_i^{(k)} - 2(x_i^{(k)})^T g_i^{(k)} \in \mathbb{R}. \end{align}$

(2.2)

Then, the $i$th inner-loop iteration of Algorithm 1 at the $k$th outer-loop iteration is then to solve the following subproblem:

$x_i^{(k + 1)} : = \arg \mathop {\max }\limits_{\left\| {x_i } \right\|_2 = 1} m^{(k)} (x_i ) = (x_i^{(k)})^T A_{ii} x_i^{(k)} + 2(x_i^{(k)})^T g_i^{(k)} + c_i^{(k)}.$

Let $\lambda _i^{(k + 1)} = (x_i^{(k + 1)})^TA_ix^{(k+1)}$, we can get the necessary and sufficient condition for the global solution $x_i^{(k + 1)}$:

$\label{eq6}\left( { - A_{ii} + \lambda _i^{(k + 1)} I^{[n_i]} } \right)x_i^{(k + 1)} = g_i^{(k)},$

(2.3)

$\label{eq7}\left\| {x_i^{(k + 1)} } \right\|_2 = 1,$

(2.4)

$\label{eq8}\left( { - A_{ii} + \lambda _i^{(k + 1)} I^{[n_i]} } \right) \ \ \ \ \text{is positive semidefinite}.$

(2.5)

If we split $A$ as

$\begin{equation} \label{eq9} A = D+U+U^T, \end{equation}$

(2.6)

where $U$ is the strictly block triangular matrix of $A$, and

$\begin{equation} D = \text{diag}\left\{A_{11}, \cdots, A_{mm} \in \mathbb{R}^{n \times n}\right\} \end{equation}$

(2.7)

is the block diagonal matrix of $A$. Then Algorithm 1 may be written in matrix form:

$\begin{equation} \Lambda ^{(k+1)}x^{(k+1)} = (D+U^T)x^{(k+1)}+Ux^{(k)}. \end{equation}$

(2.8)

We give the following SOR-like AVM.

Algorithm 2 (SOR-like AVM).

Select $x^{(0)} \in \mathcal {M}$, and set $k: = 0$.

while the stoping criterion is not met do

for $i = 1, \cdots, m$ do

$y_i^{(k + 1)} : = \arg \max _{\left\| {x_i } \right\|_2 = 1} r_i(x_1^{(k + 1)}, \cdots, x_{i - 1}^{(k + 1)}, x_i^{(k)}, x_{i + 1}^{(k)}, \cdots, x_m^{(k)} )$,

$\bar x_i^{(k + 1)} : = \omega _i y_i^{(k + 1)} + (1 - \omega _i )x_i^{(k)}$,

$x_i^{(k + 1)} : = \frac{{\bar x_i^{(k + 1)} }}{{\left\| {\bar x_i^{(k + 1)} } \right\|}}$.

end for

$k: = k+1$

end while

From the Algorithm 2, let $\xi_i^{(k)} = \left\| {\bar x_i^{(k + 1)} } \right\|$, we have

$\lambda _i^{(k + 1)} y_i^{(k + 1)} = \sum\limits_{j = 1}^i {A_{ij} x_j^{(k + 1)} } + \sum\limits_{j = i + 1}^m {A_{ij} x_j^{(k)} },$

and

$\begin{equation} \label{eq12} \xi_i^{(k)} x_i^{(k + 1)} = (1 - \omega _i^{(k)} )x_i^{(k)} + \omega _i^{(k)} y_i^{(k + 1)}. \end{equation}$

(2.9)

Then, we get

$\begin{equation} \label{eq13} \lambda _i^{(k + 1)} \xi_i^{(k)} x_i^{(k + 1)} = \lambda _i^{(k + 1)} (1 - \omega _i^{(k)} )x_i^{(k)} + \omega _i^{(k)} \sum\limits_{j = 1}^i {A_{ij} x_j^{(k + 1)} } + \omega _i^{(k)} \sum\limits_{j = i + 1}^m {A_{ij} x_j^{(k)} }. \end{equation}$

(2.10)

If we define

$\begin{equation} \label{eq14} \Xi ^{(k)} : = \text{diag}\left\{\xi ^{(k)}_1I^{\left[ n_1 \right]}, \cdots, \xi ^{(k)}_mI^{\left[n_m \right]}\right\}, \end{equation}$

(2.11)

and

$\begin{equation} \label{eq15} \Omega ^{(k)} : = \text{diag}\left\{\omega ^{(k)}_1I^{\left[n_1 \right]}, \cdots, \omega ^{(k)}_mI^{\left[n_m \right]}\right\}, \end{equation}$

(2.12)

then Algorithm 2 may be written in matrix form:

$\begin{equation} \label{eq16} \left[\Lambda ^{(k+1)}\Xi^{(k)}-\Omega^{(k)}(U^T+D)\right]x^{(k+1)} = \left[(I-\Omega^{(k)})\Lambda^{(k+1)}+\Omega^{(k)}U\right]x^{(k)}. \end{equation}$

(2.13)

Then, we can get

$\begin{align} r(x^{(k+1)}) &= (x^{(k + 1)} )^T Ax^{(k + 1)}\notag \\ & = (x^{(k + 1)} )^T (D+U+U^T)x^{(k + 1)} \notag \\ & = (x^{(k + 1)} )^T \left[{U + (\Omega ^{(k)} )^{-1} \Lambda ^{(k + 1)} \Xi ^{(k)}-(\Omega ^{(k)} )^{-1} \\\left( {\Lambda ^{(k + 1)} \Xi ^{(k)} - \Omega ^{(k)} (U^T + D)} \right)} \right]x^{(k + 1)}\notag \\ & = (x^{(k + 1)} )^T \left[{U + (\Omega ^{(k)} )^{-1} \Lambda ^{(k + 1)} \Xi ^{(k)} } \right]x^{(k + 1)} \\ &\quad - (x^{(k + 1)} )^T (\Omega ^{(k)} )^{ - 1} \left[{(I- \Omega ^{(k)} )\Lambda ^{(k + 1)} + \Omega ^{(k)} U} \right]x^{(k)}, \end{align}$

(2.14)

and

$\begin{align} r(x^{(k)}) & = (x^{(k)} )^T Ax^{(k)}\notag \\ & = (x^{(k)} )^T (D+U+U^T)x^{(k)} \notag \\ & = (x^{(k)} )^T \left[{\left( {(I-\Omega ^{(k)} )\Lambda ^{(k + 1)} + \Omega ^{(k)} U} \right)(\Omega ^{(k)} )^{-1}-\\(I - \Omega ^{(k)} )\Lambda ^{(k + 1)} (\Omega ^{(k)} )^{ - 1} + U^T + D} \right]x^{(k)} \notag \\ & = (x^{(k + 1)} )^T \left[{\Lambda ^{(k + 1)} \Xi ^{(k)}-\Omega ^{(k)} (U + D)} \right](\Omega ^{(k)} )^{ - 1} x^{(k)} \\ &\quad-(x^{(k)} )^T \left[{-(I-\Omega ^{(k)} )\Lambda ^{(k + 1)} (\Omega ^{(k)} )^{-1} + U^T + D} \right]x^{(k)}. \end{align}$

(2.15)

Hence, by (2.14) and (2.15),

$\begin{align} \Delta r(x^{(k)} ) & = r(x^{(k + 1)} ) - r(x^{(k)} )\notag \\ & = (x^{(k + 1)} )^T Ux^{(k + 1)} - (x^{(k)} )^T Ux^{(k)} + (x^{(k + 1)} )^T Dx^{(k)} - (x^{(k)} )^T Dx^{(k)} \\ &\quad + \sum\limits_{i = 1}^m {\frac{{\lambda _i^{(k + 1)} }}{{\omega _i^{(k)} }}} (\xi _i^{(k)} + 1 - \omega _i^{(k)} )(1 - (x_i^{(k + 1)} )^T x_i^{(k)} ) \end{align}$

(2.16)

$\begin{align} & = \frac{1}{2}\left[{r(x^{(k + 1)} )-r(x^{(k)} )} \right] - \frac{1}{2}(s^{(k)} )^T Ds^{(k)} \\ &\quad + \frac{1}{2}(s^{(k)} )^T \left[{(\Omega ^{(k)} )^{-1} \Lambda ^{(k + 1)} (\Xi ^{(k)} + I-\Omega ^{(k)} )} \right]s^{(k)}, \end{align}$

(2.17)

where $s^{(k)} : = x^{(k + 1)} - x^{(k)}$. Then we have

$\begin{equation} \label{eq21} \Delta r(x^{(k)} ) = (s^{(k)} )^T \left[{-D+(\Omega ^{(k)} )^{-1} \Lambda ^{(k + 1)} (\Xi ^{(k)} + I-\Omega ^{(k)} )} \right]s^{(k)}. \end{equation}$

(2.18)

Theorem 2.1. Let $\left\{x^{(k)}\right\}$ be generated from Algorithm 2, with the corresponding $\left\{\Lambda^{(k)}\right\}$ given in (2.8). If $\left\{{-D+(\Omega ^{(k)} )^{ - 1} \Lambda ^{(k + 1)} (\Xi ^{(k)} + I - \Omega ^{(k)} )}\right\}$ is uniformly positive definite. Then

$\begin{equation} \label{eq22} \mathop {\lim }\limits_{k \to + \infty } s^{(k)} = 0. \end{equation}$

(2.19)

Proof. Under these assumptions, and by (2.18), we know that $\left\{r(x^{(k)})\right\}$ is nondecreasing but bounded. This directly implies (2.19).

Let $a_i \in \sigma(A_{ii})$, $a_i$ is an eigenvalue of the matrix $A_{ii}$. Obviously, $a_i > 0$. If

$\left\{{-D+(\Omega ^{(k)} )^{ - 1} \Lambda ^{(k + 1)} (\Xi ^{(k)} + I - \Omega ^{(k)} )}\right\}$

is uniformly positive definite, then we have

$- a_i + \frac{{\lambda _i^{(k + 1)} (\xi _i^k - \omega _i^{(k)} + 1)}}{{\omega _i^{(k)} }} > 0.$

By direct computation, we have

$\begin{equation} \label{eq23} \omega _i^{(k)} < \frac{{\lambda _i^{(k + 1)} (\xi _i^k + 1)}}{{\lambda _i^{(k + 1)} + a_i }} < \frac{{\lambda _i^{(k + 1)} (\xi _i^k + 1)}}{{\lambda _i^{(k + 1)} }} = \xi _i^k + 1. \end{equation}$

(2.20)

By the equation (2.9), if $\omega_i^{(k)}\geq1$, we get

$\xi _i^k \left\| {x_i^{(k + 1)} } \right\|_2 = \left\| {(1 - \omega _i^{(k)} )x_i^{(k)} + \omega _i^{(k)} y_i^{(k + 1)} } \right\| \ge \omega _i^{(k)} - (\omega _i^{(k)} - 1) = 1.$

Then, we have

$\omega _i^{(k)} < 2.$

And we find that, if $\omega _i^{(k)}$ is very small, then $x_i^{(k + 1)} \to x_i^{(k)}$, so we usually let $\omega _i^{(k)}>0.5$ in our numerical experiments.

Lemma 2.2. ([1, Lemma 4.4]). Let $a_k$ be a bounded sequence of real numbers with the property $\left| {a_k - a_{k + 1} } \right| \to 0$ as $k \to 0$. If there are only finitely many limit points for the sequence, then $\left\{a_k\right\}$ converges to a unique limit point.

Theorem 2.3. Under the assumptions of Theorem 2.1, and assume further that A has n distinct eigenvalues, $\omega_i^{(k)} = \omega$ is a constant. Then the sequence $\left\{\Lambda^{(k)}\right\}$ converges as $k \to +\infty$, the sequence $\left\{x^{(k)}\right\}$ from Algorithm 2 converges monotonically to a solution of the MEP.

Proof. From the compactness, we know that $\left\{x^{(j_k)}\right\}$ has a convergent subsequence. Without loss of generality, we may assume $\left\{x^{(j_k)}\right\}$ converges to $\bar x$ and $\left\{\xi^{(j_k)}\right\}$ converges to $\bar \xi$. By (2.9), we get

$\bar \xi \bar x = (1 - \omega )\bar x + \omega \mathop {\lim }\limits_{k \to + \infty } y_i^{(k + 1)},$

and

$\left\| {\mathop {\lim }\limits_{k \to + \infty } y_i^{(k + 1)} } \right\|_2 = \left\| {\bar x} \right\|_2 = 1.$

Then we get

$\mathop {\lim }\limits_{k \to + \infty } y_i^{(k + 1)} = \bar x.$

Therefore,

$\begin{align} \mathop {\lim }\limits_{k \to + \infty } \delta x^{(k)} & = \mathop {\lim }\limits_{k \to + \infty } (Ax^{(k)} - \Lambda x^{(k)} )\notag \\ & = \mathop {\lim }\limits_{k \to + \infty } (Ax^{(k)} - \Lambda y^{(k)} )\notag \\ & = \mathop {\lim }\limits_{k \to + \infty } \left[{(U + U^T + D)x^{(k)}-(U^T + D)x^{(k)}-Ux^{(k-1)} } \right] \notag \\ & = \mathop {\lim }\limits_{k \to + \infty } Us^{(k - 1)} = 0. \end{align}$

(2.21)

So, subsequences $\left\{x^{(j_k)}, \Lambda^{(j_k)}\right\}$ converges to $\left\{\bar x, \bar \Lambda\right\}$, and $A\bar x = \bar \Lambda \bar x$.

From the compactness and the result that the MEP has only finitely many solutions ^[3], we know that $\left\{x^{(k)}, \Lambda^{(k)}\right\}$ has finitely many limit points. By (2.9), we have

$\left| {\lambda _i^{(k + 1)} - \lambda _i^{(k)} } \right| \le \sum\limits_{j = 1}^i {\left\| {A_{ij} } \right\|_2 } \left\| {x_i^{(k + 1)} - x_i^{(k)} } \right\|_2 + \sum\limits_{j = i + 1}^m {\left\| {A_{ij} } \right\|_2 } \left\| {x_i^{(k)} - x_i^{(k - 1)} } \right\|_2.$

By Theorem 2.1, we see that $\left| {\lambda _i^{(k + 1)} - \lambda _i^{(k)} } \right| \to 0$ as $k \to +\infty$. Together with Lemma 2.2, we can get that, the sequence $\left\{\Lambda^{(k)}\right\}$ converges as $k \to +\infty$, and the sequence $\left\{x^{(k)}\right\}$ from Algorithm 2 converges monotonically to a solution of the MEP.

We shall show that the SOR-like AVM is able to converge globally to the global maximizer of the MCP when $A$ is nonnegative irreducible.

Theorem 2.4. Suppose that A is a nonnegative irreducible matrix. For any $x^{(0)} \geq_e 0$, let $\left\{x^{(k)}\right\}$ be generated from Algorithm 2, with the corresponding $\left\{\Lambda^{(k)}\right\}$ given in (2.8). If $\left\{{-D+(\Omega ^{(k)} )^{ - 1} \Lambda ^{(k + 1)} (\Xi ^{(k)} + I - \Omega ^{(k)} )}\right\}$ is uniformly positive definite, then $\left\{x^{(k)}\right\}$ converges to the positive global maximizer of the MCP.

Proof. Similar to the proof of the Theorem 3.6 in ^[9].

---

3. Numerical experiments

In this section, we present our numerical experiments of the SOR-like AVM to show its efficiency, and most importantly, the effectiveness in seeking the global maximum of the MCP. All of our tests were conducted in MATLAB 7.0. We run the algorithm starting from $10^3$ random initial points, and use $\left\| {\delta x^{(k)} } \right\|_2 \le 10^{ - 6}$ as the stopping criterion.

Example 1. We choose the symmetric and positive definite matrix $BCSSTK04\in R^{132\times132}$ in the set BCSSTRUC1 from the Harwell-Boeing collection in Matrix Market as the original matrices $A$'s for the MCP, with $m = 2$, and $\left\{ {P = 50, 82} \right\}$, each of which consequently corresponds to a particular MCP problem.

Example 2. The matrix $A$ is given by

$A = \left( {\begin{array}{*{20}c} I&{A_{12} }&{A_{13} } \\ {A_{12}^T }&I&{A_{23} } \\ {A_{13}^T }&{A_{23}^T }&I \\ \end{array}} \right), A_{12} = \left( {\begin{array}{*{20}c} {0.636}&{0.126}&{0.059} \\ { - 0.021}&{0.633}&{0.049} \\ {0.016}&{0.157}&{0.521} \\ \end{array}} \right),$

$A_{13} = \left( {\begin{array}{*{20}c} {0.626}&{0.195}&{0.059} \\ { 0.035}&{0.459}&{0.129} \\ {0.048}&{0.238}&{0.426} \\ \end{array}} \right), A_{23} = \left( {\begin{array}{*{20}c} {0.709}&{0.050}&{-0.002} \\ { 0.039}&{0.532}&{0.190} \\ {0.067}&{0.258}&{0.299} \\ \end{array}} \right),$

with $m = 3$, and $\left\{ {P = 2, 3, 4} \right\}$. This example is from ^[5].

In Table 1, under the columns "Avg. Iter. #" are the average numbers of iterations needed to meet the stopping criterion. Under columns "% to Global" are the sample probabilities, out of the $10^3$ random tests, of convergence to a global maximizer. We observe, from the results, that for almost all randomly chosen starting points, the SOR-like AVM is able to reach a global maximizer as the standard AVM, for different choice of the $\omega$, the SOR-like AVM may iterate faster than the standard AVM. This is very attractive and has consequential effect on applications.

Furthermore, to demonstrate more clearly the performance of each algorithm, we plot the history of $r(x^{(k)})$, the residual $\left\| {\delta x^{(k)} } \right\|_2$, and the multivariate eigenvalues $\lambda_i^{(k)}$ in the following figures. The iteration $\left\{x^{(k)}\right\}$ of Example 1 in Figure 1, Figure 2 and Figure 3 starts from

$x^{(0)} = \left[{\begin{array}{*{20}c} {\frac{1}{{\sqrt {50} }}}, \ldots, {\frac{1}{{\sqrt {50} }}}, {\frac{1}{{\sqrt {82} }}}, \ldots, {\frac{1}{{\sqrt {82} }}} \\ \end{array}} \right]^T,$

while the sequence $\left\{x^{(k)}\right\}$ of Example 2 in Figure 4, Figure 5 and Figure 6 starts from

$x^{(0)} = \left[{\begin{array}{*{20}c} 1, 0, 1, 0, 0, 1, 0, 0, 0 \\ \end{array}} \right]^T.$

---

4. Conclusion

In this paper, we prove the monotone convergence of the SOR-like alternating variable method. The SOR-like AVM shows its computational advantage, for different choice of the $\omega$, the SOR-like AVM may have a better performance than the standard AVM.

---

Acknowledgments

This research is supported by Science and technology Foundation of Guizhou Province (Qian Ke He Ji Chu [2016]1161); Guizhou Province natural science foundation in China (Qian Jiao He KY [2016]255); The doctoral scientific research foundation of Zunyi Normal College(BS[2015]09);Highlevel innovative talents of Guizhou Province(Zun Ke He Ren Cai[2017]8).

---

Conflict of interest

The author declares no conflicts of interest in this paper.

---