Neural-Network-Based Approach for Extracting Eigenvectors and Eigenvalues of Real Normal Matrices and Some Extension to Real Matrices

Abstract

This paper introduces a novel neural-network-based approach for extracting some eigenpairs of real normal matrices of order n. Based on the proposed algorithm, the eigenvalues that have the largest and smallest modulus, real parts, or absolute values of imaginary parts can be extracted, respectively, as well as the corresponding eigenvectors. Although the ordinary differential equation on which our proposed algorithm is built is only n-dimensional, it can succeed to extract n-dimensional complex eigenvectors that are indeed 2n-dimensional real vectors. Moreover, we show that extracting eigen-pairs of general real matrices can be reduced to those of real normal matrices by employing the norm-reducing skill. Numerical experiments verified the computational capability of the proposed algorithm.

1. Introduction

The problem of extracting special eigenpairs of real matrices has attracted much attention both in theory [1–4] and in many engineering fields such as real-time signal processing [5–8] and principal or minor component analysis [9–12]. For example, we may wish to get eigenvectors and the corresponding eigenvalues that (1) have the largest or smallest modulus; (2) have the largest or smallest real parts; (3) have the largest or smallest imaginary parts in absolute value. Two most popular methods for this problem are the power method and the Rayleigh quotient method in their direct forms or in the context of inverse iteration [13]. Recently, many neural-network-based methods have also been proposed to solve this problem [14–23]. However, most of those neural network based methods focused on computing eigenpairs of real symmetric matrices. The following two ordinary differential equations (ODEs): d𝑥(𝑡)(𝑡)(𝑡)d𝑡=𝐴𝑥−𝑥^𝑇(𝑡)𝑥(𝑡),(𝑡)𝐴𝑥d𝑥(𝑡)d𝑡=𝑥^𝑇𝑥(𝑡)𝐴𝑥(𝑡)−𝑥(𝑡)^𝑇𝐴𝑥(𝑡)𝑥(𝑡) (1)(2) were proposed by [19, 23], respectively, where 𝐴 is a real symmetric matrix. Both (1) and (2) are efficient to compute the largest eigenvalue of 𝐴, as well as the corresponding eigenvector. In addition, they can succeed to compute the smallest eigenvalue of 𝐴 and the corresponding eigenvector by simply replacing 𝐴 with −𝐴, for example, d𝑥(𝑡)(𝑡)d𝑡=−𝑥^𝑇𝑥(𝑡)𝐴𝑥(𝑡)+𝑥(𝑡)^𝑇𝐴𝑥(𝑡)𝑥(𝑡). (3)

The following ODE for solving the generalized eigenvalue problem was proposed by [22] d𝑥(𝑡)(𝑡)(𝑥(𝑡(𝑡),d𝑡=𝐴𝑥−𝑓))𝐵𝑥 (4) where 𝐴 and 𝐵 are two real symmetric matrices, and 𝑓 can be a general form to some degree. Particularly, if 𝐵 is the identity matrix, (4) can be used to solve the standard eigenvalue problem as (1) and (2).

References [16–18] extended those neural network based approaches to the case of real antisymmetric or special real matrices of order 𝑛, where the proposed neural networks can be summarized by 2𝑛-dimensional ODEs since eigenvectors in those cases may be 𝑛-dimensional complex vectors, that is, 2𝑛-dimensional real vectors.

In this paper, we propose an approach for extracting six types of eigenvalues of 𝑛-by-𝑛 real normal matrices and the corresponding eigenvectors based on (2) or (3). Although eigenvectors of real normal matrices may be 𝑛-dimensional complex vectors, the computation of our proposed method can be achieved in 𝑛-dimensional real vector space, which can reduce the scale of networks a lot. Then, we show that any real matrix can be made arbitrarily close to a normal matrix by a series of similarity transformations, based on which our proposed algorithm can be extended to the case of arbitrary real matrices.

2. Main Results

Let 𝑖=√−1 be the imaginary unit, 𝑥 the conjugate of 𝑥, and diag[𝐴₁,…,𝐴_𝑝] a block diagonal matrix, where 𝐴_𝑖, 𝑖=1,…,𝑝, is a square matrix at the 𝑖th diagonal block. Unless specially stated, 𝐴 is a real normal matrix in this paper.

Lemma 1. In [13], 𝐴 is a real normal matrix of order 𝑛 if and only if there exists an orthogonal matrix 𝑈 such that 𝑈^𝑇𝐴𝑈=diag𝜆₁,…,𝜆_𝑟𝑟,𝐴_𝑟+1,…,𝐴_{𝑟+𝑘𝑘},𝑟+2𝑘=𝑛,0≤𝑟≤𝑛, (5) where 𝑈=𝑢₁,…,𝑢_{𝑟𝑟realeigenvectors},𝑢^𝑅_𝑟+1,𝑢^𝐼_𝑟+1,…,𝑢^𝑅_𝑟+𝑘,𝑢^𝐼_{𝑟+𝑘𝑘pairsofcomplexeigenvectors},𝐴_𝑟+𝑗=𝑎_𝑟+𝑗𝑏_𝑟+𝑗−𝑏_𝑟+𝑗𝑎_𝑟+𝑗,𝑗=1,…,𝑘. (6)
Here, 𝜆_𝑗, 𝑗=1,…,𝑟, are 𝑟 real eigenvalues of 𝐴 corresponding to the 𝑟 real eigenvectors 𝑢_𝑗, and 𝑎_𝑟+𝑠±𝑖𝑏_𝑟+𝑠,  𝑠=1,…,𝑘, are 𝑘 pairs of complex eigenvalues of 𝐴 corresponding to the 𝑘 pairs of complex eigenvectors 𝑢^𝑅_𝑟+𝑠±𝑖𝑢^𝐼_𝑟+𝑠.

For simplicity, let 𝑎_𝑗=𝜆_𝑗, 𝑏_𝑗=0 for 𝑗=1,…,𝑟. Based on (5) and (6), it is straightforward to verify 𝑈^𝑇𝐴−𝐴^𝑇𝑇𝐴−𝐴^𝑇𝑈=4diag𝑏²₁,…,𝑏²_𝑟,𝑏²_𝑟+1,𝑏²_𝑟+1,…,𝑏²_𝑟+𝑘,𝑏²_𝑟+𝑘𝑈^𝑇𝐴+𝐴^𝑇𝑈=2diag𝑎₁,…,𝑎_𝑟,𝑎_𝑟+1,𝑎_𝑟+1,…,𝑎_𝑟+𝑘,𝑎_𝑟+𝑘𝑈^𝑇𝐴𝐴^𝑇𝑈=diag𝑎²₁,…,𝑎²_𝑟,𝑎²_𝑟+1+𝑏²_𝑟+1,𝑎²_𝑟+1+𝑏²_𝑟+1,…,𝑎²_𝑟+𝑘+𝑏²_𝑟+𝑘,𝑎²_𝑟+𝑘+𝑏²_𝑟+𝑘. (7)(8)(9)

Then, the following six definitions and two lemmas are presented, which will be involved much in the sequel. (C1) Let 𝒥₁={𝑠∣|𝑏_𝑠|=max{|𝑏₁|,…,|𝑏_𝑟+𝑘|}. Then for all 𝑗,𝑘∈𝒥₁,  𝑎_𝑗=𝑎_𝑘. (C2) Let 𝒥₂={𝑠∣|𝑏_𝑠|=min{|𝑏₁|,…,|𝑏_𝑟+𝑘|}. Then for all 𝑗,𝑘∈𝒥₂,  𝑎_𝑗=𝑎_𝑘. (C3) Let 𝒥₃={𝑠∣𝑎_𝑠=max{𝑎₁,…,𝑎_𝑟+𝑘}. Then for all 𝑗,𝑘∈𝒥₃,  |𝑏_𝑗|=|𝑏_𝑘|. (C4) Let 𝒥₄={𝑠∣𝑎_𝑠=min{𝑎₁,…,𝑎_𝑟+𝑘}. Then for all 𝑗,𝑘∈𝒥₄,  |𝑏_𝑗|=|𝑏_𝑘|. (C5) Let 𝒥₅={𝑠∣𝑎²_𝑠+𝑏²_𝑠=max{𝑎²₁+𝑏²₁,…,𝑎²_𝑟+𝑘+𝑏²_𝑟+𝑘}. Then for all 𝑗,𝑘∈𝒥₅,  𝑎_𝑗=𝑎_𝑘. (C6) Let 𝒥₆={𝑠∣𝑎²_𝑠+𝑏²_𝑠=min{𝑎²₁+𝑏²₁,…,𝑎²_𝑟+𝑘+𝑏²_𝑟+𝑘}. Then for all 𝑗,𝑘∈𝒥₆,  𝑎_𝑗=𝑎_𝑘.

Lemma 2 (Theorem  4 in [23]). Assume that nonzero 𝑥(0)∈𝑅^𝑛 is not orthogonal to the eigensubspace corresponding to the largest eigenvalue of 𝐴. Then, the solution of (2) starting from 𝑥(0) converges to an eigenvector corresponding to the largest eigenvalue of 𝐴 that is equal to lim_𝑡→+∞(𝑥(𝑡)^𝑇𝐴𝑥(𝑡)/𝑥(𝑡)^𝑇𝑥(𝑡)).

Lemma 3 (Theorem  5 in [23]). Assume that nonzero 𝑥(0)∈𝑅^𝑛 is not orthogonal to the eigensubspace corresponding to the smallest eigenvalue of 𝐴. Then, the solution of (3) starting from 𝑥(0) converges to an eigenvector corresponding to the smallest eigenvalue of 𝐴 that is equal to lim_𝑡→+∞(𝑥(𝑡)^𝑇𝐴𝑥(𝑡)/𝑥(𝑡)^𝑇𝑥(𝑡)).

Remark 4. If we randomly choose 𝑥(0), the projection of 𝑥(0) on the eigensubspace corresponding to the largest or smallest eigenvalue of 𝐴 will be nonzero with high probability. Hence, (2) and (3) can almost work well with randomly generated 𝑥(0).

2.1. Computing the Eigenvalues with the Largest or Smallest Imaginary Parts in Absolute Value, as well as the Corresponding Eigenvectors

Without loss of generality, in this subsection we assume that 𝑏₁=⋯=𝑏_𝑟=0<𝑏²_𝑟+1≤⋯≤𝑏²_𝑟+𝑘 in (7). Note that (𝐴−𝐴^𝑇)^𝑇(𝐴−𝐴^𝑇)=−(𝐴−𝐴^𝑇)².

Based on (7), we know that 0 (if any) is the eigenvalue of −(𝐴−𝐴^𝑇)² corresponding to the eigenvector 𝑢_𝑗, 𝑗=1,…,𝑟, and that 4𝑏²_𝑟+𝑠, 𝑠=1,…,𝑘, is the eigenvalue of −(𝐴−𝐴^𝑇)² corresponding to the eigenvectors 𝑢^𝑅_𝑟+𝑠 or 𝑢^𝐼_𝑟+𝑠. If replacing 𝐴 with the symmetric matrix −(𝐴−𝐴^𝑇)² in (2) as d𝑥(𝑡)(𝑡)d𝑡=−𝑥^𝑇𝑥(𝑡)𝐴−𝐴^𝑇2𝑥(𝑡)+𝑥(𝑡)^𝑇𝐴−𝐴^𝑇2𝑥(𝑡)𝑥(𝑡) (10) based on Lemma 2, we know that 𝛼=lim_𝑡→∞𝑥(𝑡) is an eigenvector of −(𝐴−𝐴^𝑇)² corresponding to 4𝑏²_𝑟+𝑘 and 4𝑏²_𝑟+𝑘=lim_𝑡→∞(−𝑥(𝑡)^𝑇(𝐴−𝐴^𝑇)²𝑥(𝑡)/𝑥(𝑡)^𝑇𝑥(𝑡)) (thus, |𝑏_𝑟+𝑘| has been gotten), where 𝑥(𝑡) is a solution of (10). If 𝑏_𝑟+𝑘=0, 𝑏_𝑗’s are all zero, meaning that 𝐴 is symmetric. In this case, we can directly use (2) or (3) to get the largest or smallest eigenvalue of 𝐴 and the corresponding eigenvector, respectively. Hence, we assume 𝑏_𝑟+𝑘≠0.

The following lemma introduces an approach for computing 𝑎_𝑟+𝑘 under the condition (C1).

Lemma 5. Assume that (C1) holds. Let 𝛼=lim_𝑡→∞𝑥(𝑡), where 𝑥(𝑡) is a solution of (10). Then, 𝑎_𝑟+𝑘=𝛼^𝑇(𝐴+𝐴^𝑇)𝛼/2𝛼^𝑇𝛼. In addition, there exists 𝑎²_𝑟+𝑘+𝑏²_𝑟+𝑘=𝛼^𝑇(𝐴𝐴^𝑇)𝛼/𝛼^𝑇𝛼.

Proof. Assume that 4𝑏²_𝑟+𝑗,  𝑗=𝑠,𝑠+1,…,𝑘,  𝑠≥1, are the largest eigenvalues of −(𝐴−𝐴^𝑇)², that is; |𝑏_𝑟+𝑠|=|𝑏_𝑟+𝑠+1|=⋯=|𝑏_𝑟+𝑘|. Therefore, 𝒥₁={𝑟+𝑠,…,𝑟+𝑘}. (11)
Based on Lemma 2 and (7), we know that 𝛼 should be a linear combination of 𝑢^𝑅_𝑗 and 𝑢^𝐼_𝑗,  𝑗∈𝒥₁. Let 𝛼=_𝑗∈𝒥1𝛾_𝑗𝑢^𝑅_𝑗+𝜏_𝑗𝑢^𝐼_𝑗,𝛾_𝑗,𝜏_𝑗∈𝑅. (12)
In addition, by (8) we have 𝐴+𝐴^𝑇𝑢^𝑅_𝑗=2𝑎_𝑗𝑢^𝑅_𝑗,𝐴+𝐴^𝑇𝑢^𝐼_𝑗=2𝑎_𝑗𝑢^𝐼_𝑗,𝑗∈𝒥₁. (13)
And by (9), we have 𝐴𝐴^𝑇𝑢^𝑅_𝑗=𝑎²_𝑗+𝑏²_𝑗𝑢^𝑅_𝑗,𝐴𝐴^𝑇𝑢^𝐼_𝑗=𝑎²_𝑗+𝑏²_𝑗𝑢^𝐼_𝑗,𝑗∈𝒥₁. (14)
Because 𝑈 is an orthogonal matrix and for all 𝑗∈𝒥₁, 𝑎_𝑗=𝑎_𝑟+𝑘 holds due to (C1), it is straightforward to verify 𝛼^𝑇𝐴+𝐴^𝑇𝛼2𝛼^𝑇𝛼=2𝑎_𝑟+𝑘∑_𝑗∈𝒥1𝛾²_𝑗+𝜏²_𝑗2∑_𝑗∈𝒥1𝛾²_𝑗+𝜏²_𝑗=𝑎_𝑟+𝑘,𝛼^𝑇𝐴𝐴^𝑇𝛼𝛼^𝑇𝛼=𝑎²_𝑟+𝑘+𝑏²_𝑟+𝑘∑_𝑗∈𝒥1𝛾²_𝑗+𝜏²_𝑗∑_𝑗∈𝒥1𝛾²_𝑗+𝜏²_𝑗=𝑎²_𝑟+𝑘+𝑏²_𝑟+𝑘, (15) thus proving the lemma.

Remark 6. Based on Lemma 5, if 𝑎²_𝑟+𝑘+𝑏²_𝑟+𝑘≠𝛼^𝑇(𝐴𝐴^𝑇)𝛼/𝛼^𝑇𝛼, (C1) does not hold surely, which can be used to check whether (C1) holds or not.
The following lemma introduces an approach for computing a pair of conjugated eigenvectors of 𝐴 corresponding to the eigenvalues 𝑎_𝑟+𝑘±𝑖|𝑏_𝑟+𝑘| under the condition (C1).

Lemma 7. Assume that (C1) holds. Given any nonzero 4𝑏²_𝑟+𝑘, the largest eigenvalue of −(𝐴−𝐴^𝑇)², and the corresponding eigenvector 𝛼 obtained by (10), let 𝑎_𝑟+𝑘=𝛼^𝑇(𝐴+𝐴^𝑇)𝛼/2𝛼^𝑇𝛼 and 𝛽=(𝐴−𝐴^𝑇)𝛼/2|𝑏_𝑟+𝑘|. Then, 𝑎_𝑟+𝑘±𝑖|𝑏_𝑟+𝑘| are two eigenvalues of 𝐴 corresponding to the eigenvectors 𝛽±𝑖𝛼, respectively.

Proof. Let 𝒥₁ and 𝛼 take the form as (11) and (12), respectively. Based on (13), we can write 𝐴+𝐴^𝑇2𝛼=𝑎_𝑟+𝑘𝛼. (16)
Following the decomposition of 𝐴 as (5) and the definition of 𝛽, we have 𝑏_𝑟+𝑘𝛽=𝐴−𝐴^𝑇𝛼2=𝑈^𝑇diag0,…,0_𝑟,𝐵_𝑟+1,…,𝐵_𝑟+𝑘𝑈𝛼2, (17) where 𝐵_𝑗=02𝑏_𝑗−2𝑏_𝑗0,𝑗∈{𝑟+1,…,𝑟+𝑘}⧵𝒥₁02𝑏_𝑟+𝑘−2𝑏_𝑟+𝑘0,𝑗∈𝒥₁. (18)
Taking (12) into (17), we get 𝛽=_𝑗∈𝒥1𝜏_𝑗𝑢^𝑅_𝑗−𝛾_𝑗𝑢^𝐼_𝑗. (19)
Based on (13) and (19), it is straightforward to verify 𝐴+𝐴^𝑇2𝛽=𝑎_𝑟+𝑘𝛽. (20)
In addition, since 4𝑏²_𝑟+𝑘 is the eigenvalue of −(𝐴−𝐴^𝑇)² corresponding to the eigenvector 𝛼, we have −(𝐴−𝐴^𝑇)²𝛼=4𝑏²_𝑟+𝑘𝛼. Hence, −(𝐴−𝐴^𝑇)|𝑏_𝑟+𝑘|𝛽=2𝑏²_𝑟+𝑘𝛼. Since 𝑏_𝑟+𝑘≠0, we have 𝐴−𝐴^𝑇2𝛽=−𝑏_𝑟+𝑘𝛼. (21)
By (20) and (21), we have 𝐴𝛽=𝐴+𝐴^𝑇2𝛽+𝐴−𝐴^𝑇2𝛽=𝑎_𝑟+𝑘𝛽−𝑏_𝑟+𝑘𝛼. (22)
By (16) and (17), we have 𝐴𝛼=𝐴+𝐴^𝑇2𝛼+𝐴−𝐴^𝑇2𝛼=𝑎_𝑟+𝑘𝛼+𝑏_𝑟+𝑘𝛽. (23)
Then, it is straightforward to verify 𝐴𝛽+𝑖𝛼=𝑎_𝑟+𝑘+𝑖𝑏_𝑟+𝑘,𝐴=𝑎𝛽+𝑖𝛼𝛽−𝑖𝛼_𝑟+𝑘−𝑖𝑏_𝑟+𝑘𝛽−𝑖𝛼, (24) thus proving the lemma.

To get |𝑏₁|, we can in advance get 4𝑏²₁, the smallest eigenvalue of −(𝐴−𝐴^𝑇)² by replacing 𝐴 with −(𝐴−𝐴^𝑇)² in (3) as follows: d𝑥(𝑡)(𝑡)d𝑡=𝑥^𝑇𝑥(𝑡)𝐴−𝐴^𝑇2𝑥(𝑡)−𝑥(𝑡)^𝑇𝐴−𝐴^𝑇2𝑥(𝑡)𝑥(𝑡). (25)

Let 𝛼^∗=lim_𝑡→∞𝑥(𝑡), where 𝑥(𝑡) is a solution of (25). From Lemma 3, we know 4𝑏²₁=−(𝛼^∗)^𝑇(𝐴−𝐴^𝑇)²𝛼^∗/(𝛼^∗)^𝑇𝛼^∗. Then, the following lemma similar to Lemma 5 can be used to compute 𝑎₁.

Lemma 8. Assume that (C2) holds. Then, 𝑎₁=(𝛼^∗)^𝑇(𝐴+𝐴^𝑇)𝛼^∗/2(𝛼^∗)^𝑇𝛼^∗. In addition, there exists 𝑎²₁+𝑏²₁=(𝛼^∗)^𝑇(𝐴𝐴^𝑇)𝛼^∗/(𝛼^∗)^𝑇𝛼^∗.

Proof. The proof is almost the same to that in Lemma 5.

Note that 𝑏₁ may be zero; that is, −(𝐴−𝐴^𝑇)² has real eigenvalues. In this case, we have the following lemma.

Lemma 9. Assume that (C2) holds and 0 is the smallest eigenvalue of −(𝐴−𝐴^𝑇)² corresponding to the eigenvector 𝛼^∗. Then, 𝑎₁ is the eigenvalue of 𝐴 corresponding to the eigenvector 𝛼^∗, where 𝑎₁=(𝛼^∗)^𝑇(𝐴+𝐴^𝑇)𝛼^∗/2(𝛼^∗)^𝑇𝛼^∗.

Proof. Following the conditions, we have −(𝐴−𝐴^𝑇)²𝛼^∗=0. Hence, (𝛼^∗)^𝑇(𝐴−𝐴^𝑇)^𝑇(𝐴−𝐴^𝑇)𝛼^∗=‖(𝐴−𝐴^𝑇)𝛼^∗‖²₂=0; that is, 𝐴−𝐴^𝑇𝛼^∗=0. (26)
Note that 𝑏₁=⋯=𝑏_𝑟=0,  𝑟≥1, because 0 is the smallest eigenvalue of −(𝐴−𝐴^𝑇)². Based on the definition of 𝒥₂, we have 𝑎₁=⋯=𝑎_𝑟,𝒥₂={1,…,𝑟},𝑟≥1. (27)
Applying Lemma 3 to (25), we know that 𝛼^∗ should be a linear combination of 𝑢₁,…,𝑢_𝑟. Let 𝛼^∗=_𝑗∈𝒥2𝛾_𝑗𝑢_𝑗,𝛾_𝑗∈𝑅. (28)
By (8), we have 𝐴+𝐴^𝑇𝑢_𝑗=2𝑎_𝑗𝑢_𝑗=2𝑎₁𝑢_𝑗,𝑗∈𝒥₂. (29)
Therefore, 𝐴+𝐴^𝑇2𝛼^∗=12_𝑗∈𝒥2𝛾_𝑗𝐴+𝐴^𝑇𝑢_𝑗=𝑎₁𝛼^∗. (30)
Then, by (26) and (30), it is straightforward to verify 𝐴𝛼^∗=𝐴+𝐴^𝑇2𝛼^∗+𝐴−𝐴^𝑇2𝛼^∗=𝑎₁𝛼^∗, (31)
thus proving the lemma.

In the case of 𝑏₁≠0, that is, all the eigenvalues of −(𝐴−𝐴^𝑇)² are complex numbers, we have the following lemma similar to Lemma 7.

Lemma 10. Assume that (C2) holds. Given any nonzero 4𝑏²₁, the smallest eigenvalue of −(𝐴−𝐴^𝑇)², and the corresponding eigenvector 𝛼^∗ obtained by (25), let 𝑎₁=(𝛼^∗)^𝑇(𝐴+𝐴^𝑇)𝛼^∗/2(𝛼^∗)^𝑇𝛼^∗ and 𝛽^∗=(𝐴−𝐴^𝑇)𝛼^∗/2|𝑏₁|. Then, 𝑎₁±𝑖|𝑏₁| are two eigenvalues of 𝐴 corresponding to the eigenvectors 𝛽^∗±𝑖𝛼^∗, respectively.

Proof. The proof is almost the same to that in Lemma 7.

Remark 11. Among the following four real normal matrices ,,,,70000050000050000062000−2670000051000−150000062000−267000005−2000250000052000−2570000052000−250000062000−26 (32) only the first three meet (C1), but the last one does not. And only the last three matrices meet (C2), but the first one does not.

2.2. Computing the Eigenvalues with the Largest or Smallest Real Parts, as well as the Corresponding Eigenvectors

As shown in (8), 2𝑎_𝑗,  𝑗=1,…,𝑟+𝑘, are the eigenvalues of the symmetric matrix (𝐴+𝐴^𝑇). In this subsection, we assume that 𝑎₁≤⋯≤𝑎_𝑟 and 𝑎_𝑟+1≤⋯≤𝑎_𝑟+𝑘, which can be achieved by reordering 𝑎_𝑗 and the corresponding columns of 𝑈.

If replacing 𝐴 with (𝐴+𝐴^𝑇) in (2), we get d𝑥(𝑡)(𝑡)d𝑡=𝑥^𝑇𝑥(𝑡)𝐴+𝐴^𝑇𝑥(𝑡)−𝑥(𝑡)^𝑇𝐴+𝐴^𝑇𝑥(𝑡)𝑥(𝑡). (33)

Without loss of generality, assume that 𝑎_𝑟+𝑘 is the largest real part of the eigenvalues of 𝐴 (it may be 𝑎_𝑟. However, Lemmas 12 and 13 have no difference in that case). Applying Lemma 2 to (33), we can get 2𝑎_𝑟+𝑘=lim_𝑡→∞(𝑥(𝑡)^𝑇(𝐴+𝐴^𝑇)𝑥(𝑡)/𝑥(𝑡)^𝑇𝑥(𝑡)), the largest eigenvalue of (𝐴+𝐴^𝑇), and the corresponding eigenvector 𝜁=lim_𝑡→∞𝑥(𝑡), where 𝑥(𝑡) is a solution of (33).

The following lemma introduces an approach for computing |𝑏_𝑟+𝑘| under the condition (C3).

Lemma 12. Assume that (C3) holds. Let 𝜁=lim_𝑡→∞𝑥(𝑡), where 𝑥(𝑡) is a solution of (33). Then, |𝑏_𝑟+𝑘|=(𝜁^𝑇(𝐴𝐴^𝑇)𝜁/𝜁^𝑇𝜁)−𝑎²_𝑟+𝑘.

Proof. By Lemma 2, we know that 𝜁 should be a linear combination of 𝑢^𝑅_𝑗 and 𝑢^𝐼_𝑗,  𝑗∈𝒥₃. Let 𝜁=_𝑗∈𝒥3𝛾_𝑗𝑢^𝑅_𝑗+𝜏_𝑗𝑢^𝐼_𝑗,𝛾_𝑗,𝜏_𝑗∈𝑅. (34)
Then, based on (9), we get 𝐴𝐴^𝑇𝑢^𝑅_𝑗=𝑎²_𝑗+𝑏²_𝑗𝑢^𝑅_𝑗,𝐴𝐴^𝑇𝑢^𝐼_𝑗=𝑎²_𝑗+𝑏²_𝑗𝑢^𝐼_𝑗,𝑗∈𝒥₃. (35)
Because 𝑈 is an orthogonal matrix and for all 𝑗∈𝒥₃,  𝑎_𝑗=𝑎_𝑟+𝑘 and |𝑏_𝑗|=|𝑏_𝑟+𝑘| hold, that is, 𝑎²_𝑗+𝑏²_𝑗=𝑎²_𝑟+𝑘+𝑏²_𝑟+𝑘 for all 𝑗∈𝒥₃, we have 𝜁^𝑇𝐴𝐴^𝑇𝜁𝜁^𝑇𝜁=𝑎²_𝑟+𝑘+𝑏²_𝑟+𝑘∑_𝑗∈𝒥3𝛾²_𝑗+𝜏²_𝑗∑_𝑗∈𝒥3𝛾²_𝑗+𝜏²_𝑗=𝑎²_𝑟+𝑘+𝑏²_𝑟+𝑘, (36) thus proving the lemma.

The following lemma introduces an approach for computing a pair of conjugated eigenvectors of 𝐴 corresponding to the eigenvalues 𝑎_𝑟+𝑘±𝑖|𝑏_𝑟+𝑘| under the condition (C3).

Lemma 13. Assume that (C3) holds. Given any nonzero 2𝑎_𝑟+𝑘, the largest eigenvalue of (𝐴+𝐴^𝑇), and the corresponding eigenvector 𝜁 obtained by (33), let |𝑏_𝑟+𝑘|=(𝜁^𝑇(𝐴𝐴^𝑇)𝜁/𝜁^𝑇𝜁)−𝑎²_𝑟+𝑘. If 𝑏_𝑟+𝑘=0, 𝑎_𝑟+𝑘 is an eigenvalue of 𝐴 corresponding to the eigenvector 𝜁. If 𝑏_𝑟+𝑘≠0, let 𝜂=(𝐴−𝐴^𝑇)𝜁/2|𝑏_𝑟+𝑘|. Then, 𝑎_𝑟+𝑘±𝑖|𝑏_𝑟+𝑘| are two eigenvalues of 𝐴 corresponding to the eigenvectors 𝜂±𝑖𝜁, respectively.

Proof. Combining the proofs of Lemmas 9 and 7, we can prove this lemma.

If replacing 𝐴 with (𝐴+𝐴^𝑇) in (3), we get d𝑥(𝑡)(𝑡)d𝑡=−𝑥^𝑇𝑥(𝑡)𝐴+𝐴^𝑇𝑥(𝑡)+𝑥(𝑡)^𝑇𝐴+𝐴^𝑇𝑥(𝑡)𝑥(𝑡). (37)

Without loss of generality, assume that 𝑎_𝑟+1 is the smallest real part of the eigenvalues of 𝐴 (it may be 𝑎₁. However, Lemma 14 has no difference in that case). Applying Lemma 3 to (37), we can obtain 2𝑎_𝑟+1, the smallest eigenvalue of (𝐴+𝐴^𝑇), as well as the corresponding eigenvector, denoted by 𝜁^∗. Then, we have the following lemma.

Lemma 14. Assume that (C4) holds. Given any nonzero 2𝑎_𝑟+1, the smallest eigenvalue of (𝐴+𝐴^𝑇), and the corresponding eigenvector 𝜁^∗ obtained by (37), let |𝑏_𝑟+1|=((𝜁^∗)^𝑇(𝐴𝐴^𝑇)𝜁^∗/(𝜁^∗)^𝑇𝜁^∗)−𝑎²_𝑟+1. If 𝑏_𝑟+1=0, 𝑎_𝑟+1 is the eigenvalue of 𝐴 corresponding to the eigenvector 𝜁^∗. If 𝑏_𝑟+1≠0, let 𝜂^∗=(𝐴−𝐴^𝑇)𝜁^∗/2|𝑏_𝑟+1|. Then, 𝑎_𝑟+1±𝑖|𝑏_𝑟+1| are two eigenvalues of 𝐴 corresponding to the eigenvectors 𝜂^∗±𝑖𝜁^∗, respectively.

Proof. Combining the proofs of Lemmas 12, 9, and 7, we can prove this lemma.

Remark 15. Among the following four real normal matrices ,,,,10000040000050000012000−2130000040000040000012000−2130000042000−240000042000−2430000040000050000052000−25 (38) only the first three meet (C3), but the last one does not. And only the last three matrices meet (C4), but the first one does not.

2.3. Computing the Eigenvalues with the Largest or Smallest Modulus, as well as the Corresponding Eigenvectors

Reorder the eigenvalues of the symmetric matrix 𝐴𝐴^𝑇 in (9) and the corresponding columns of 𝑈 such that 𝑎²₁≤⋯≤𝑎²_𝑟 and 𝑎²_𝑟+1+𝑏²_𝑟+1≤⋯≤𝑎²_𝑟+𝑘+𝑏²_𝑟+𝑘. Without loss of generality, assume that 𝑎_𝑟+𝑘±𝑖|𝑏_𝑟+𝑘| are the eigenvalues of 𝐴 that have the largest modulus (it may be 𝑎_𝑟. However, Lemma 16 has no difference in that case), and that 𝑎_𝑟+1±𝑖|𝑏_𝑟+1| are the eigenvalues of 𝐴 that have the smallest modulus (it may be 𝑎₁. However, Lemma 17 has no difference in that case).

Replacing 𝐴 with 𝐴𝐴^𝑇 in (2), we get d𝑥(𝑡)(𝑡)d𝑡=𝑥^𝑇𝑥(𝑡)𝐴𝐴^𝑇𝑥(𝑡)−𝑥(𝑡)^𝑇𝐴𝐴^𝑇𝑥(𝑡)𝑥(𝑡). (39)

Applying Lemma 2 to (39), we can obtain 𝑎²_𝑟+𝑘+𝑏²_𝑟+𝑘=lim_𝑡→∞(𝑥(𝑡)^𝑇𝐴^𝑇𝐴𝑥(𝑡)/𝑥(𝑡)^𝑇𝑥(𝑡)), the largest eigenvalue of 𝐴𝐴^𝑇, and the corresponding eigenvector 𝜃=lim_𝑡→∞𝑥(𝑡), where 𝑥(𝑡) is the solution of (39). Then, we have the following lemma.

Lemma 16. Assume that (C5) holds. Given any 𝑎²_𝑟+𝑘+𝑏²_𝑟+𝑘, the largest eigenvalue of 𝐴𝐴^𝑇, and the corresponding eigenvector 𝜃 obtained by (39). Then, 𝑎_𝑟+𝑘=𝜃^𝑇(𝐴+𝐴^𝑇)𝜃/2𝜃^𝑇𝜃. Thus, |𝑏_𝑟+𝑘| can be gotten. If 𝑏_𝑟+𝑘=0, 𝑎_𝑟+𝑘 is an eigenvalue of 𝐴 corresponding to the eigenvector 𝜃. If 𝑏_𝑟+𝑘≠0, let 𝜔=(𝐴−𝐴^𝑇)𝜃/|2𝑏_𝑟+𝑘|. Then, 𝑎_𝑟+𝑘±𝑖|𝑏_𝑟+𝑘| are two eigenvalues of 𝐴 corresponding to the eigenvectors 𝜔±𝑖𝜃, respectively.

Proof. Combining the proofs of Lemmas 5, 7, and 9, we can prove this lemma.

Replacing 𝐴 with 𝐴𝐴^𝑇 in (3), we get d𝑥(𝑡)(𝑡)d𝑡=−𝑥^𝑇𝑥(𝑡)𝐴𝐴^𝑇𝑥(𝑡)+𝑥(𝑡)^𝑇𝐴𝐴^𝑇𝑥(𝑡)𝑥(𝑡). (40)

Applying Lemma 3 to (40), we can obtain 𝑎²_𝑟+1+𝑏²_𝑟+1, the smallest eigenvalue of 𝐴𝐴^𝑇, and the corresponding eigenvector, denoted by 𝜃^∗. Then, we have the following lemma.

Lemma 17. Assume (C6) holds. Given any 𝑎²_𝑟+1+𝑏²_𝑟+1, the smallest eigenvalue of 𝐴𝐴^𝑇, and the corresponding eigenvector 𝜃^∗ obtained by (40). Then, 𝑎_𝑟+1=(𝜃^∗)^𝑇(𝐴+𝐴^𝑇)𝜃^∗/2(𝜃^∗)^𝑇𝜃^∗. So one can get that |𝑏_𝑟+1|. If 𝑏_𝑟+1=0, 𝑎_𝑟+1 is an eigenvalue of 𝐴 corresponding to the eigenvector 𝜃^∗. If 𝑏_𝑟+1≠0, let 𝜔^∗=(𝐴−𝐴^𝑇)𝜃^∗/2|𝑏_𝑟+1|. Then, 𝑎_𝑟+1±𝑖|𝑏_𝑟+1| are two eigenvalues of 𝐴 corresponding to the eigenvectors 𝜔^∗±𝑖𝜃^∗, respectively.

Proof. Combining the proofs of Lemmas 5, 7, and 9, we can prove this lemma.

Remark 18. Among the following four real normal matrices ,,,,200000−20000060000012000−2130000060000060000012000−21300003−4004300003000−430000040000050000034000−43 (41) only the first three meet (C5), but the last one does not. And only the last three meet (C6), but the first one does not.
However, there exists some specially constructed 𝐴 that meet none of (C1) to (C6), for example, 𝐴=diag[𝐴₁,𝐴₂,𝐴₃,𝐴₄], where 𝐴₁=−41−1−4,𝐴₂=−43−3−4,𝐴₃=,𝐴43−34₄=.41−14 (42)
Nevertheless, a randomly generated real normal matrix can meet (C1) to (C6) with high probability.

2.4. Extension to Arbitrary Real Matrices

In this subsection, 𝐴 is an arbitrary real matrix of order 𝑛. Let ‖𝐴‖ be the Frobenius norm of 𝐴, and, 𝜆_𝑗,  𝑗=1,…,𝑛, be the eigenvalue of 𝐴. Denote the set of all complex nonsingular matrices by 𝒯.

By the Schur inequality [13], we know ^𝑛_𝑗=1𝜆²_𝑗≤‖𝐴‖² (43) with equality if and only if 𝐴 is a normal matrix. Since the spectrum of 𝐴 does not change by a similarity transformation, the inequality ^𝑛_𝑗=1𝜆²_𝑗≤𝑃⁻¹𝐴𝑃²,𝑃∈𝒯 (44) holds with equality if and only if 𝑃⁻¹𝐴𝑃 is a normal matrix. In addition, [24] proved ^𝑛_𝑖=1𝜆²_𝑖=inf_𝑃∈𝒯𝑃⁻¹𝐴𝑃². (45)

Based on (45), if we can find a sequence 𝐴_𝑗 as follows: 𝐴₁=𝐴,𝐴_𝑗+1=𝑃_𝑗⁻¹𝐴_𝑗𝑃_𝑗,𝐴_𝑗+1<𝐴_𝑗,𝑃_𝑗∈𝒯,𝑗=1,2,… (46) such that ^𝑛_𝑖=1𝜆²_𝑖=lim_𝑗→∞𝐴_𝑗², (47)

where 𝐴_∞=lim_𝑗→∞𝐴_𝑗 is to be a normal matrix with the same eigenvalues as 𝐴. Such skill, termed as the norm-reducing technique, has been proposed by [25–27]. Moreover, following the idea presented by [26], it is easy to find that when 𝐴 is real, 𝐴_∞ can be chosen to be a real normal matrix.

In a word, any real matrix 𝐴 can be translated into a real normal matrix 𝐴_∞=𝑃_∞⁻¹𝐴𝑃_∞ by a similarity transformation 𝑃_∞=𝑃₁𝑃₂⋯. Typical approaches for constructing 𝑃_𝑗 can be found in [26]. Note that if 𝜆 is the eigenvalue of 𝐴_∞ corresponding to the eigenvector 𝑢, 𝜆 is the eigenvalue of 𝐴 corresponding to the eigenvector 𝑃_∞𝑢. Hence, our proposed algorithm can be extended to extract eigenpairs of arbitrary real matrices by employing the norm-reducing technique.

Without loss of generality, we use the following random matrix 𝐴 as an example to describe the norm-reducing technique: 𝐴=1.05410.4052−1.0191−0.57710.4158−1.9797−0.7025−1.3852−0.83640.0430−1.86741.49900.95490.8530−0.9489−1.83240.1378−0.60110.47730.54160.8486−1.5868−1.17190.3023−0.8211. (48)

The Frobenius norm of 𝐴 is ‖𝐴‖=√27.7629, and the eigenvalues of matrix 𝐴 are 𝜆₁=2.8797, 𝜆₂=𝜆₃=−0.6900+1.8770𝑖, 𝜆₄=−1.5144, and 𝜆₅=0.9774; obviously ‖𝐴‖>∑⁵_𝑖=1𝜆²_𝑖=19.5401, so 𝐴 is a nonnormal real matrix. According to the approach that presented in [26], we can construct 𝑃_𝑗, 𝑗=1,2,…. If the condition ‖𝐴_𝑖−𝐴_𝑖+1‖<10⁻¹⁰ can be satisfied, we break the iteration. After 129 iterations, we can obtain a approximate normal matrix 𝐴₁₃₀: 𝐴₁₃₀=−0.74410.00071.63340.82100.0010−0.00352.8797−0.00380.0000−0.0000−1.4162−0.0048−0.91320.82230.0034−1.1560−0.0019−0.1181−1.2370−0.00140.0011−0.00000−0.00360.9774, (49) from which we have ‖𝐴^∗₁₃₀𝐴₁₃₀−𝐴₁₃₀𝐴^∗₁₃₀‖=√0.00121231723027107 and ‖𝐴₁₃₀‖=√19.540133042214 that is very close to ∑⁵_𝑖=1𝜆²_𝑖=√19.5400526824864; so the matrix 𝐴₁₃₀ can be regarded as a normal matrix in practical application, the corresponding 𝑃_∞=𝑃₁𝑃₂𝑃₃⋯𝑃₁₂₉ are as follows: 𝑃_∞=−0.3499−0.5290−0.00900.2228−0.0151−0.4809−0.0259−0.6935−0.40790.1559−0.29140.72110.67430.03950.3437−0.44840.1475−0.36680.2444−0.87040.0611−0.32650.4294−1.1028−0.5166, (50) which satisfy 𝐴₁₃₀=𝑃_∞⁻¹𝐴𝑃_∞. We also can see that the eigenvalues of 𝐴₁₃₀ are 𝜆₁=−1.5144,  𝜆₂=𝜆₃=−0.6900+1.8770𝑖, 𝜆₄=2.8797, and 𝜆₅=0.9774, which are just the eigenvalues of original nonnormal matrix 𝐴. We presented here the transient behavior of ‖𝐴_𝑖‖² in Figure 1, where 𝑖=1,2,3,… is the iterations and 𝐴₁=𝐴.

Figure 1 

Trajectory of ‖𝐴𝑖‖2.

In the following, we use 𝑁=1000 random matrices to verify the average performance of the norm-reducing technique. Let MeanΔ²_𝑖=1𝑁^𝑁_𝑗=1Δ²_𝑗𝑖 (51) denotes the average measure of a large number of nonnormal matrices in a statistical sense at 𝑖th iteration, where Δ²_𝑗𝑖=‖𝐴^𝑇_𝑗𝑖𝐴_𝑗𝑖−𝐴_𝑗𝑖𝐴^𝑇_𝑗𝑖‖² is the measure of the nonnormal matrix 𝐴_𝑗 at 𝑖th iteration, Δ²_𝑗𝑖=0 if and only if 𝐴_𝑗 is normal matrix at 𝑖th iteration. We also presented the dynamic behavior trajectory of MeanΔ² in Figure 2, from which we can see that for most of nonnormal matrices, after 350 iterations, the MeanΔ² is very close to zero.

Figure 2 

Average performance trajectory of norm-reducing technique.

3. Neural Implementation Description

In the presented paper, we mainly focus on the classical neural network differential equation as shown in (2), where 𝐴=(𝑎_𝑖𝑗), 𝑖,𝑗=1,2,…,𝑛 are symmetric matrices that need to calculate eigenvalues and the corresponding eigenvectors, 𝑥(𝑡)=[𝑥₁(𝑡),𝑥₂(𝑡),…,𝑥_𝑛(𝑡)]^𝑇 is a column vector which denotes the states of neurons in the neural network dynamic system, and the elements of symmetric matrix 𝐴 denote the connection weights between those neurons. We presented the schematic diagram of the neural network in Figure 3, from which we can see that it is a recurrent neural network since the input is just the output of the system.

Figure 3 

The schematic diagram of the neural network equation (2) for solving the eigenvector corresponding to the modulus largest eigenvalues of real symmetric matrix 𝐴.

In the practical applications, we often only need a nonzero column vector 𝑥(0)=[𝑥₁(0),𝑥₂(0),…,𝑥_𝑛(0)]^𝑇 to start the neural network system by the following update rule: 𝑥𝑗+1=𝑥𝑗+𝜏𝑥𝑗^𝑇𝑥𝑗𝑗𝑗𝐴𝑥−𝑥^𝑇𝐴𝑥𝑗𝑥𝑗, (52) where 𝑗 denote the 𝑗th iteration and 𝜏 is a small time step. The iteration stops once ‖𝑥(𝑗+1)−𝑥(𝑗)‖<𝜀, where 𝜀 is a small constraint error that can be set in advance. If ‖𝑥(𝑗+1)−𝑥(𝑗)‖<𝜀, we could regard that ‖𝑥(𝑗+1)−𝑥(𝑗)‖=0, that is, 𝑥(𝑗)^𝑇𝑥(𝑗)𝐴𝑥(𝑗)−𝑥(𝑗)^𝑇𝐴𝑥(𝑗)𝑥(𝑗)=0; so we have 𝐴𝑥(𝑗)=(𝑥(𝑗)^𝑇𝐴𝑥(𝑗)/𝑥(𝑗)^𝑇𝑥(𝑗))𝑥(𝑗), according to the theory in [23], 𝑥(𝑗) is the eigenvector corresponding to the modulus largest eigenvalue which can be denoted as 𝑥(𝑗)^𝑇𝐴𝑥(𝑗)/𝑥(𝑗)^𝑇𝑥(𝑗).

4. Examples and Discussion

Three experiments are presented to verify our results. The following real normal matrix 𝐴 (randomly generated) was used in those three experiments: 𝐴=1.8818−2.1563−0.5979−0.47680.72700.4133−1.42670.52581.49460.71250.47002.53840.42730.77842.86450.85871.3792−1.0153−0.19291.3273−1.31253.16290.25730.63291.5966−0.6593−0.37990.20582.98450.64041.06271.9973−0.0249−1.57270.92892.3259. (53)

Using the [𝑉,𝐷]=eig(𝐴) function in Matlab, we got the eigenvalues of 𝐴 as 𝜆₁=−2.3420, 𝜆₂=1.005, 𝜆₃=𝜆₄=3.1967+1.7860𝑖, and 𝜆₅=𝜆₆=4.3445+1.2532𝑖, as well as the corresponding eigenvectors 𝑢₁,…,𝑢₆ as follows: 𝑢₁=,𝑢0.37540.7510−0.2146−0.16640.0280−0.4696₂=,𝑢0.6657−0.02420.29940.2594−0.58880.2295₃=,𝑢−0.0294+0.1297𝑖−0.0243−0.1385𝑖0.6165−0.1917+0.5472𝑖0.0913+0.2785𝑖−0.2708−0.2950𝑖₅=,𝑢0.3207+0.2956𝑖0.0472−0.4424𝑖0.0711−0.2167𝑖−0.1146−0.3213𝑖0.49030.3692−0.2583𝑖₄=𝑢₃,𝑢₆=𝑢₅. (54)

We can see that all of (C1) to (C6) hold except (C2). For simplicity, denote lim_𝑡→∞𝑥(𝑡) by 𝑥(∞).

Example 19 (for Section 2.1). We used (10) with the following initial condition (randomly generalized) 𝑥(0)=[]0.2329,−0.0023,−0.8601,−0.2452,−0.1761,−0.1577^𝑇 (55) to get |𝑏₃| (the largest absolute value of the imaginary part among 𝜆_𝑖) that −𝑥(𝑡)^𝑇(𝐴−𝐴^𝑇)²𝑥(𝑡)/4𝑥(𝑡)^𝑇𝑥(𝑡) should converge to. From Lemma 5, 𝑥(𝑡)^𝑇(𝐴+𝐴^𝑇)𝑥(𝑡)/2𝑥(𝑡)^𝑇𝑥(𝑡)→𝑎₃ should hold. By Lemma 7, 𝑥(𝑡) and (𝐴−𝐴^𝑇)𝑥(𝑡)/2|𝑏₃| should converge to the imaginary and real parts of an eigenvector corresponding to 𝜆₃, respectively. The transient behaviors of the above four variables are shown in Figures 4, 5, and 6, respectively. After convergence, we saw 𝐴−𝐴^𝑇𝑥(∞)2𝑏₃(∞)𝑖==()𝑢+𝑥0.1795−0.0005𝑖−0.1747+0.0737𝑖−0.1857−0.8107𝑖0.7772+0.0872𝑖0.3386−0.2039𝑖−0.3064+0.4450𝑖−0.3012−1.3149𝑖₃. (56)
Thus, the estimated complex vector is an eigenvector of 𝐴 corresponding to 𝜆₃.
Although we can use (25) to get the smallest absolute value of the imaginary part among 𝜆_𝑖 (it is zero in this experiment), neither the corresponding real part nor the eigenvector can be obtained from Lemmas 8 or 9 since (C2) does not hold.

Figure 4 

Trajectories of 𝜙(𝑡)=−𝑥(𝑡)𝑇(𝐴−𝐴𝑇)2𝑥(𝑡)/4𝑥(𝑡)𝑇𝑥(𝑡) and 𝜓(𝑡)=𝑥(𝑡)𝑇(𝐴+𝐴𝑇)𝑥(𝑡)/2𝑥(𝑡)𝑇𝑥(𝑡) based on (10) with initial 𝑥(0) as (55), which should converge to |𝑏3| and 𝑎3, respectively.

Figure 5 

Trajectories of 𝑥(𝑡), the solution of (10) with initial 𝑥(0) as (55), which should converge to the imaginary part of an eigenvector of 𝐴 corresponding to the eigenvalue 𝜆3.

Figure 6 

Trajectories of 𝑦(𝑡)=(𝐴−𝐴𝑇)𝑥(𝑡)/2𝑏3 based on (10) with initial 𝑥(0) as (55), which should converge to the real part of an eigenvector of 𝐴 corresponding to the eigenvalue 𝜆3.

Example 20 (for Section 2.2). We used (33) with the same 𝑥(0) as (55) to get 𝑎₅ (the largest real part among 𝜆_𝑖) that 𝑥(𝑡)^𝑇(𝐴+𝐴^𝑇)𝑥(𝑡)/2𝑥(𝑡)^𝑇𝑥(𝑡) should converge to. Based on Lemma 12, (𝑥(𝑡)^𝑇(𝐴𝐴^𝑇)𝑥(𝑡)/𝑥(𝑡)^𝑇𝑥(𝑡))−𝑎²₅→|𝑏₅| should hold. The transient behaviors of such two variables are shown in Figure 7. Since 𝑏₅≠0, from Lemma 13, 𝑥(𝑡) and (𝐴−𝐴^𝑇)𝑥(𝑡)/2|𝑏₅| should converge to the imaginary and real parts of an eigenvector corresponding to 𝜆₅, as shown in Figures 8 and 9. After convergence, we saw 𝐴−𝐴^𝑇𝑥(∞)2𝑏₅(∞)𝑖==()𝑢+𝑥−0.5695−0.1474𝑖0.2254+0.5561𝑖0.0534+0.3029𝑖0.3407+0.3091𝑖−0.5830+0.3120𝑖−0.2746+0.5420𝑖−1.1890+0.6363𝑖₅. (57)
Hence, the estimated complex vector is an eigenvector of 𝐴 corresponding to 𝜆₅.
Based on (37), we got 𝑎₁ (the smallest real part among 𝜆_𝑖). After convergence, we saw that 𝑥(∞)^𝑇(𝐴+𝐴^𝑇)𝑥(∞)/2𝑥(∞)^𝑇𝑥(∞)=𝑎₁=𝜆₁=−2.3420, (𝑥(∞)^𝑇(𝐴𝐴^𝑇)𝑥(∞)/𝑥(∞)^𝑇𝑥(∞))−𝑎²₁=|𝑏₁|=0, and 𝑥(∞) was equal to 𝑥(∞)=−0.3581−0.71650.20480.1588−0.02670.4480=−0.9541𝑢₁, (58) just as expected from Lemma 14.

Figure 7 

Trajectories of 𝜓(𝑡)=𝑥(𝑡)𝑇(𝐴+𝐴𝑇)𝑥(𝑡)/2𝑥(𝑡)𝑇𝑥(𝑡) and 𝜙(𝑡)=(𝑥(𝑡)𝑇(𝐴𝐴𝑇)𝑥(𝑡)/𝑥(𝑡)𝑇𝑥(𝑡))−𝑎25 based on (33) with initial 𝑥(0) as (55), which should converge to 𝑎5 and |𝑏5|, respectively.

Figure 8 

Trajectories of 𝑥(𝑡), the solution of (33) with initial 𝑥(0) as (55), which should converge to the imaginary part of an eigenvector of 𝐴 corresponding to the eigenvalue 𝜆5.

Figure 9 

Trajectories of 𝑦(𝑡)=(𝐴−𝐴𝑇)𝑥(𝑡)/2|𝑏5| based on (33) with initial 𝑥(0) as (55), which should converge to the real part of an eigenvector of 𝐴 corresponding to the eigenvalue 𝜆5.

Example 21 (for Section 2.3). Based on (39) and Lemma 16, we can get 𝜆₅ and one corresponding eigenvector again because |𝜆₅| is the largest modulus of 𝜆_𝑖. In addition, we used (40) with the same 𝑥(0) as (55) to get |𝜆₂|, the smallest modulus of 𝜆_𝑖. By Lemma 17, 𝑎₂=𝑥(∞)^𝑇(𝐴+𝐴^𝑇)𝑥(∞)/2𝑥(∞)^𝑇𝑥(∞) and |𝑏₂|=(𝑥(∞)^𝑇(𝐴𝐴^𝑇)𝑥(∞)/𝑥(∞)^𝑇𝑥(∞))−𝑎²₂ should hold. The transient behaviors of such two variables are shown in Figure 10. After convergence, we saw that 𝑏₂=0. Hence, from Lemma 17, 𝑥(𝑡) should converge to a constant multiple of 𝑢₂, which is shown in Figure 11. As expected, 𝑥(∞) was equal to 𝑥(∞)=−0.63370.0230−0.2850−0.24690.5605−0.2185=−0.9519𝑢₂. (59)

Figure 10 

Trajectories of 𝜓(𝑡)=𝑥(𝑡)𝑇(𝐴+𝐴𝑇)𝑥(𝑡)/2𝑥(𝑡)𝑇𝑥(𝑡) and 𝜙(𝑡)=(𝑥(𝑡)𝑇(𝐴𝐴𝑇)𝑥(𝑡)/𝑥(𝑡)𝑇𝑥(𝑡))−𝑎22 based on (40) with initial 𝑥(0) as (55), which should converge to 𝑎2 and |𝑏2|, respectively.

Figure 11 

Trajectories of 𝑥(𝑡), the solution of (40), with initial 𝑥(0) as (55), which should converge to an eigenvector of 𝐴 corresponding to the eigenvalue 𝜆2.

Example 22 (for extension to arbitrary real matrices). According to the theory in Section 2.1, we present here an experiment for arbitrary real matrices to verify the effectiveness of the norm-reducing technique in Section 2.4. Considering the following nonnormal matrix 𝐴: 𝐴=3.4633−3.74872.9281−0.5106−0.04554.48603.4259−3.20105.06251.21200.5268−1.03111.5386−0.4766−1.4478−2.47681.25852.05252.3539−0.84870.0498−3.1267−6.5297−1.51062.2083, (60) the Frobenius norm of 𝐴 is ‖𝐴‖=√187.8869>∑⁵_𝑖=1𝜆²_𝑖=97.1780, where 𝜆_𝑖,  𝑖=1,2,…,5 are the eigenvalues of matrix 𝐴, so matrix 𝐴 is a nonnormal matrix. According to the theory in Section 2.4, the eigenpairs problem of nonnormal real matrix can be converted into the eigenpairs problem of the corresponding normal matrix 𝐴_𝑝, which can be calculated by norm-reducing technique: 𝐴_𝑝=2.34380.79151.10711.8624−2.6867−1.59182.10161.72492.20671.98430.8659−1.74344.5103−1.27080.04081.7777−2.3015−0.98983.58150.7201−2.4688−2.3136−0.47460.05450.4528, (61) and the corresponding 𝑃_∞ as follows: 𝑃_∞=28.2974−28.8442115.669816.476513.5575−35.968948.54754.346558.13849.0642−10.2373−14.297614.116823.6416−40.2882−3.4360−25.4122−114.241715.019115.1720−79.1532−35.346575.4347−24.8455−48.2893, (62) which satisfies the relationship 𝐴_𝑝=𝑃_∞⁻¹𝐴𝑃_∞. In order to solve the eigenpairs of matrix 𝐴, we have to calculate the eigenpairs of the matrix 𝐴_𝑝 at first. Direct calculations of the eigenvalues of matrix 𝐴_𝑝 are 𝜆₁=−1.8342, 𝜆₂=𝜆₃=2.1114+3.786𝑖,  𝜆₄=𝜆₅=5.3007+0.1328𝑖 and those are just the eigenvalues of nonnormal matrix 𝐴 and the corresponding eigenvectors as follows: 𝑢₁=0.60790.0886−0.1195−0.28050.7278,𝑢₂==0.0035−0.2240𝑖0.70200.0189+0.3280𝑖0.0502+0.4231𝑖−0.0660+0.4040𝑖𝑢₃,𝑢₄==0.1769−0.4835𝑖−0.0129+0.0561𝑖0.6204−0.2522−0.4643𝑖−0.1415+0.2181𝑖𝑢₅. (63)
The eigenvectors of nonnormal matrix 𝐴, which corresponding to the eigenvalues with largest imaginary parts in absolute are 𝑢^𝐴₂==−0.2353−0.5331𝑖−0.64160.0449−0.0581𝑖0.3856+0.0270𝑖0.0747−0.2973𝑖𝑢^𝐴₃. (64)
According to the results above, the eigenvalues with largest imaginary in absolute are 𝜆₂ and 𝜆₃, and let 𝑎₂±|𝑏₂|𝑖 denote them.
We used (10) with the following initial condition (randomly generalized): 𝑥(0)=[0.3891,0.0856,0.2215,0.3151,0.4929]^𝑇 (65) to get |𝑏₂| (the largest absolute value of the imaginary part among 𝜆_𝑖) that −𝑥(𝑡)^𝑇(𝐴_𝑝−𝐴^𝑇_𝑝)²𝑥(𝑡)/4𝑥(𝑡)^𝑇𝑥(𝑡) should converge to. From Lemma 5, 𝑥(𝑡)^𝑇(𝐴_𝑝+𝐴^𝑇_𝑝)𝑥(𝑡)/2𝑥(𝑡)^𝑇𝑥(𝑡)→𝑎₂ should hold. By Lemma 7, 𝑥(𝑡) and (𝐴_𝑝−𝐴^𝑇_𝑝)𝑥(𝑡)/2|𝑏₂| should converge to the imaginary and real parts of an eigenvector corresponding to 𝜆₂, respectively. The transient behaviors of the above four variables are shown in Figures 12, 13, and 14, respectively.
After convergence, we saw 𝐴_𝑝−𝐴^𝑇_𝑝𝑥(∞)2𝑏₂(∞)𝑖==()𝑢+𝑥0.0394−0.2317𝑖0.7277+0.1121𝑖−0.0327+0.3430𝑖−0.0155+0.4466𝑖−0.1329+0.4083𝑖1.0367+0.1597i₂. (66)
Thus, the estimated complex vector is an eigenvector of 𝐴_𝑝 corresponding to 𝜆₂. According to the theory in Section 2.4, the estimated complex eigenvector of nonnormal matrix 𝐴 corresponding to the eigenvalue 𝜆₂ should be 𝑃_∞𝐴_𝑝−𝐴^𝑇_𝑝𝑥(∞)2𝑏₂(∞)𝑖==()𝑢+𝑥−25.7189+42.7862𝑖31.6606+44.9313𝑖−6.2833−0.2807𝑖−17.1372−28.3402𝑖−24.5113+9.4386𝑖−49.3476−70.0318𝑖^𝐴₂, (67) from which we can see that the estimated complex vector is an eigenvector of 𝐴 corresponding to 𝜆₂.

Figure 12 

Trajectories of 𝜙(𝑡)=−𝑥(𝑡)𝑇(𝐴𝑝−𝐴𝑇𝑝)2𝑥(𝑡)/4𝑥(𝑡)𝑇𝑥(𝑡) and 𝜓(𝑡)=𝑥(𝑡)𝑇(𝐴𝑝+𝐴𝑇𝑝)𝑥(𝑡)/2𝑥(𝑡)𝑇𝑥(𝑡) based on (10) with initial 𝑥(0) as (65), which should converge to |𝑏2| and 𝑎2, respectively.

Figure 13 

Trajectories of 𝑥(𝑡), the solution of (10) with initial 𝑥(0) as (65), which should converge to the imaginary part of an eigenvector of 𝐴𝑝 corresponding to the eigenvalue 𝜆2.

Figure 14 

Trajectories of 𝑦(𝑡)=(𝐴𝑝−𝐴𝑇𝑝)𝑥(𝑡)/2𝑏2 based on (10) with initial 𝑥(0) as (65), which should converge to the real part of an eigenvector of 𝐴𝑝 corresponding to the eigenvalue 𝜆2.

5. Conclusion

This paper introduces a neural network based approach for computing eigenvectors of real normal matrices and the corresponding eigenvalues that have the largest or smallest modulus, have the largest or smallest real part, and have the largest or smallest imaginary part in absolute value. All the computation can be carried out in real vector space although eigenpairs may be complex, which can reduce the scale of networks a lot. We also shed light on extending this method to the case of general real matrices by employing the norm-reducing technique proposed in other literatures. Four simulation examples verified the validity of our proposed algorithm.