WO2008018188A1

WO2008018188A1 - Eigen value decomposing device and eigen value decomposing method

Info

Publication number: WO2008018188A1
Application number: PCT/JP2007/051575
Authority: WO
Inventors: Yoshimasa Nakamura; Hiroaki Tsuboi; Taro Konda; Masashi Iwasaki; Masami Takata
Original assignee: Kyoto University
Priority date: 2006-08-08
Filing date: 2007-01-31
Publication date: 2008-02-14
Also published as: JPWO2008018188A1; JP5017666B2; US20100185716A1; US8255447B2

Abstract

[PROBLEMS] An eigenvalue decomposing device capable of performing parallel processing at high speed with high accuracy. [MEANS FOR SOLVING PROBLEMS] An eigenvalue decomposing device comprises a matrix decomposing section (14) for repeating division of a symmetric tridiagonal matrix (T) into two symmetric tridiagonal matrices, an eigenvalue decomposing section (15) for performing eigenvalue decomposition of the divided symmetric tridiagonal matrices, an eigenvalue calculating section (17) for repeating calculation of the eigenvalue of the symmetric tridiagonal matrix before the division and the matrix elements of the symmetric tridiagonal matrix before the division from the eigenvalue subjected to the eigenvalue decomposition by the eigenvalue decomposition section (15) and matrix elements which are a part of the elements of an orthogonal matrix composed of eigenvectors until the eigenvalue of the symmetric tridiagonal matrix (T) is calculated, and an eigenvector calculating section (19) for calculating the eigenvectors of the symmetric tridiagonal matrix (T) by the twist decomposition form the symmetric tridiagonal matrix T and its eigenvalue.

Description

Specification

Eigenvalue decomposition apparatus and eigenvalue decomposition method

Technical field

The present invention relates to an eigenvalue decomposition apparatus that performs eigenvalue decomposition.

Background art

[0002] Eigenvalue decomposition is used in a very wide range of fields such as chemical calculation by molecular orbital method, statistical calculation, and information retrieval. Also, a method of eigenvalue decomposition is known (for example, see Non-Patent Document 1).

Non-Patent Document 1: J. J. M. Cuppen, "A divide and conquer method for the symmetric tridiagonal eigenproblemj, Numerische Mathematik, Vol. 3 6, p. 177-195, 1981

Disclosure of the invention

Problems to be solved by the invention

[0003] In recent years, with an increase in the amount of data in these application fields, fast and highly accurate eigenvalue decomposition is required. In addition, it has been desired to develop an eigenvalue decomposition apparatus that can process eigenvalue decomposition in parallel.

[0004] The present invention has been made under the above circumstances, and an object thereof is to provide an eigenvalue decomposition apparatus and the like capable of executing eigenvalue decomposition with high parallelism and high speed and high accuracy.

Means for solving the problem

In order to achieve the above object, an eigenvalue decomposition apparatus according to the present invention includes a diagonal matrix storage unit that stores a symmetric tridiagonal matrix T, and the symmetric tridiagonal matrix from the diagonal matrix storage unit. T is read out, the symmetric tridiagonal matrix T is divided into two symmetric tridiagonal matrices and stored in the diagonal matrix storage unit, and the symmetric tridiagonal matrix is divided into two symmetric tridiagonal matrices. A matrix dividing unit that repeats the process of dividing into diagonal matrices and storing in the diagonal matrix storage unit until each symmetric tridiagonal matrix after division is less than or equal to a predetermined size; and Eigenvalues that store the eigenvalues of each symmetric tridiagonal matrix that is less than or equal to the predetermined size and the matrix elements that are part of the orthogonal matrix that is the eigenvector force of each symmetric tridiagonal matrix Disassembly The storage unit and each symmetric tridiagonal matrix having a size less than or equal to the predetermined size are read from the diagonal matrix storage unit, and eigenvalue decomposition is performed on each symmetric tridiagonal matrix. Calculate at least eigenvalues of each symmetric tridiagonal matrix and matrix elements of an orthogonal matrix consisting of eigenvectors of each symmetric tridiagonal matrix, and store the eigenvalues and the matrix elements in the eigenvalue decomposition storage unit The eigenvalue decomposition unit, the eigenvalue storage unit storing the eigenvalues of the symmetric tridiagonal matrix T, the eigenvalues of each symmetric tridiagonal matrix, and the matrix elements are read from the eigenvalue decomposition storage unit, and the eigenvalues And the matrix element, the eigenvalues of the symmetric tridiagonal matrix of the split source and the matrix elements of the symmetric tridiagonal matrix of the split source are calculated and stored in the eigenvalue decomposition storage unit. The eigenvalues of the symmetric tridiagonal matrix of the split source, the matrix elements and The calculation process is repeated until at least one eigenvalue of the symmetric tridiagonal matrix T is calculated, and at least one eigenvalue of the symmetric tridiagonal matrix T is stored in the eigenvalue storage unit. , The eigenvector storage unit storing the eigenvectors of the symmetric tridiagonal matrix T, the symmetric tridiagonal matrix T is read from the diagonal matrix storage unit, and the symmetric tridiagonal matrix is read from the eigenvalue storage unit. Read the eigenvalues of the matrix T, calculate at least one eigenvector of the symmetric tridiagonal matrix T from the symmetric tridiagonal matrix T and its eigenvalues using the twist decomposition method, and store the eigenvectors And an eigenvector calculation unit that accumulates in the unit.

With such a configuration, first, eigenvalues are calculated, and eigenvalues are calculated using the twist decomposition method based on the eigenvalues, whereby high-speed and high-precision eigenvalue decomposition can be realized. . It also has excellent parallelism. Furthermore, when it is not necessary to calculate all eigenvalues and eigenvalues, eigenvalues and eigenvectors can be calculated within a necessary range, and the processing load can be reduced.

[0007] In the eigenvalue decomposition apparatus according to the present invention, the eigenvector calculation unit may calculate an eigenvector by a qd-type twist decomposition method.

With such a configuration, high-speed and high-precision eigenvalue decomposition can be realized. It also has excellent parallelism.

[0008] In the eigenvalue decomposition apparatus according to the present invention, the eigenvector calculation unit reads the eigenvalues of the symmetric tridiagonal matrix T from the eigenvalue storage unit, and When the force is a negative value, the diagonal matrix storage unit also reads the symmetric tridiagonal matrix T and changes the eigenvector of the symmetric tridiagonal matrix Τ. The symmetric tridiagonal matrix Τ is converted to a positive definite value so that all eigenvalues of the angular matrix Τ are positive values, and the eigenvalues after the positive definite value 匕 are stored in the eigenvalue storage unit. A symmetric tridiagonal matrix Τ stored in the diagonal matrix storage unit, and a symmetric tridiagonal matrix ある in which all eigenvalues are positive values from the diagonal matrix storage unit Τ And reading out the positive eigenvalues of the symmetric tridiagonal matrix から from the eigenvalue storage unit and performing qd conversion on each element of the symmetric tridiagonal matrix Τ to obtain the symmetric tridiagonal matrix Τ. over double diagonal matrix B diagonal matrix T ^{(+ 1)} and Cholesky decompose Kore under double diagonal matrix B ^(_1) A key decomposition unit, and each element of the upper double diagonal matrix B ⁽⁺¹⁾ and the lower double diagonal matrix B ^(_1), the eigenvalues of the positive values of the symmetric tridiagonal matrix T using A vector calculation unit that calculates an eigenvector and stores it in the eigenvector storage unit.

[0009] With such a configuration, when a negative eigenvalue is included, the qd-type twist decomposition method can be executed after making the symmetric tridiagonal matrix T positive definite. In addition, it can be calculated faster than the LV twist method.

[0010] In the eigenvalue decomposition apparatus according to the present invention, the eigenvector calculation unit may calculate an eigenvector by an LV type twist decomposition method.

With such a configuration, the eigenvector can be calculated numerically stably.

In the eigenvalue decomposition apparatus according to the present invention, the eigenvector calculation unit reads the symmetric tridiagonal matrix T from the diagonal matrix storage unit, and the symmetric tridiagonal matrix T from the eigenvalue storage unit. And by performing Miura transform, dLVv type transform, and inverse Miura transform on each element of the symmetric tridiagonal matrix T, the symmetric tridiagonal matrix T is transformed into an upper bidiagonal matrix B ( ^{+ 1)} and the Cholesky decomposing Cholesky decomposition portion under double diagonal matrix B ^(_1), each of the upper double diagonal matrix B ⁽⁺¹⁾ and the lower double diagonal matrix B ^(_1) A vector calculation unit that calculates an eigenvector using the elements and eigenvalues of the symmetric tridiagonal matrix T and stores the eigenvector in the eigenvector storage unit;

With such a configuration, it is possible to calculate the eigenvector numerically stably by using a plurality of auxiliary variables in the Cholesky decomposition process. In the eigenvalue decomposition apparatus according to the present invention, the Cholesky decomposition unit includes a plurality of Cholesky decomposition units, and the plurality of Cholesky decomposition units perform Cholesky decomposition on the symmetric tridiagonal matrix T. You can also execute parallel processes.

With such a configuration, the Cholesky decomposition process can be performed in a short time.

In the eigenvalue decomposition apparatus according to the present invention, the vector calculation unit includes a plurality of vector calculation units, and the plurality of vector calculation units execute a process of calculating the eigenvectors in parallel. Also good.

With this configuration, the eigenvector calculation process can be performed in a short time.

[0014] Further, in the eigenvalue decomposition apparatus according to the present invention, the eigenvalue calculating unit includes a plurality of eigenvalue calculating means, and the plurality of eigenvalue calculating means includes eigenvalues of a symmetric tridiagonal matrix of a division source, Processing for calculating matrix elements may be executed in parallel.

With such a configuration, the process of calculating the eigenvalue can be performed in a short time.

[0015] In the eigenvalue decomposition apparatus according to the present invention, the eigenvalue decomposition unit includes a plurality of eigenvalue decomposition means, and the plurality of eigenvalue decomposition means performs eigenvalue decomposition on a symmetric tridiagonal matrix. May be executed in parallel.

With this configuration, the eigenvalue decomposition process can be performed in a short time.

[0016] Further, in the eigenvalue decomposition apparatus according to the present invention, a matrix storage unit that stores a symmetric matrix A, the symmetric matrix A is read from the matrix storage unit, and the symmetric matrix A is tridiagonalized. A diagonalization unit that calculates a tridiagonal matrix T and accumulates it in the diagonal matrix storage unit.

With such a configuration, an eigenvalue can be calculated for an arbitrary symmetric matrix A.

It is also possible to calculate the eigenvector of the symmetric matrix A by using the eigenvector of the symmetric tridiagonal matrix T.

[0018] Further, in the eigenvalue decomposition apparatus according to the present invention, the matrix dividing unit may divide the symmetric tridiagonal matrix into two symmetric tridiagonal matrices of approximately half.

With such a configuration, parallel processing can be appropriately performed.

In the eigenvalue decomposition apparatus according to the present invention, the eigenvalue calculation unit may calculate all eigenvalues of the symmetric tridiagonal matrix T. In the eigenvalue decomposition apparatus according to the present invention, the eigenvector calculation unit may calculate all eigenvectors of the symmetric tridiagonal matrix T.

In the eigenvalue decomposition apparatus according to the present invention, the symmetric matrix A may be a covariance matrix in which an average force of a vector indicating face image data is calculated.

With such a configuration, the eigenvalue decomposition apparatus can be used for facial image data recognition processing.

The invention's effect

[0022] The eigenvalue decomposition apparatus and the like according to the present invention can perform eigenvalue decomposition at high speed and with high accuracy. It also has excellent parallelism.

BEST MODE FOR CARRYING OUT THE INVENTION

Hereinafter, an eigenvalue decomposition apparatus according to the present invention will be described using embodiments. In the following embodiments, the components and steps denoted by the same reference numerals are the same or equivalent, and the description thereof may be omitted.

[Embodiment 1]

An eigenvalue decomposition apparatus according to Embodiment 1 of the present invention will be described with reference to the drawings.

FIG. 1 is a block diagram showing the configuration of the eigenvalue decomposition apparatus 1 according to this embodiment. In FIG. 1, the eigenvalue decomposition apparatus 1 according to the present embodiment includes a matrix storage unit 11, a diagonalization unit 1 2, a diagonal matrix storage unit 13, a matrix division unit 14, an eigenvalue decomposition unit 15, and an eigenvalue. The decomposition storage unit 16, the eigenvalue calculation unit 17, the eigenvalue storage unit 18, the eigenvector calculation unit 19, and the eigenvector storage unit 20 are provided.

The matrix storage unit 11 stores an arbitrary symmetric matrix A. This symmetric matrix A is a real matrix where each element is a real number. Note that the fact that symmetric matrix A is stored means that data indicating symmetric matrix A is stored. The same applies to the storage unit described later. The matrix storage unit 11 can be realized by a predetermined recording medium (for example, a semiconductor memory, a magnetic disk, an optical disk, etc.). The storage in the matrix storage unit 11 may be a temporary storage in a RAM or the like, or a long-term storage. The process in which the symmetric matrix A is stored in the matrix storage unit 11 does not matter. For example, a symmetric matrix A is stored in a matrix storage unit via a recording medium. The symmetric matrix A transmitted via the communication line etc. can be stored in the matrix storage unit 11 or input devices such as a keyboard and a mouse. The symmetric matrix A input via is stored in the matrix storage unit 11.

The diagonalization unit 12 reads the symmetric matrix A from the matrix storage unit 11 and calculates a symmetric tridiagonal matrix T obtained by tridiagonalizing the read symmetric matrix A. Then, the diagonalization unit 12 stores the calculated symmetric tridiagonal matrix T in the diagonal matrix storage unit 13. For example, the diagonalization unit 12 performs a symmetric matrix A using a method of repeatedly performing householder transformation as necessary, a method using the Lanczos method, and other tridiagonalization methods. Is tridiagonalized.

The diagonal matrix storage unit 13 stores a symmetric tridiagonal matrix T. The diagonal matrix storage unit 13 can be realized by a predetermined recording medium (for example, a semiconductor memory, a magnetic disk, an optical disk, etc.). Storage in the diagonal matrix storage unit 13 may be temporary storage in RAM or the like, or may be long-term storage.

[0028] The matrix dividing unit 14 reads the symmetric tridiagonal matrix T from the diagonal matrix storage unit 13, divides the symmetric tridiagonal matrix T into two symmetric tridiagonal matrices, and performs diagonal processing. Stores in the matrix storage unit 13. The matrix dividing unit 14 divides the symmetric tridiagonal matrix into two symmetric tridiagonal matrices and stores them in the diagonal matrix storage unit 13. It repeats recursively until it becomes less than the predetermined size.

[0029] The eigenvalue decomposition unit 15 reads each symmetric tridiagonal matrix of a predetermined size or less from the diagonal matrix storage unit 13, and performs eigenvalue decomposition on each symmetric tridiagonal matrix. The eigenvalues of each symmetric tridiagonal matrix and the eigenvectors of each symmetric tridiagonal matrix are calculated. The eigenvalue decomposition unit 15 may perform eigenvalue decomposition using, for example, a method combining the bisection method and the inverse iteration method, MR "3 method, QR method, etc. Here, eigenvalue decomposition using the MR'3 method is performed. DSTEGR provided in LAPACK3.0 (http: //www.netlib.orgZlapackZ) may be used when performing eigenvalue decomposition with the LAPACK3.0. The provided DSTEQR may be used, and the method of eigenvalue decomposition is already known and will not be described in detail. The unit 15 stores the calculated eigenvalues in the eigenvalue decomposition storage unit 16. Further, the eigenvalue decomposition unit 15 also accumulates in the eigenvalue decomposition storage unit 16 matrix elements that are part of the orthogonal matrix that is the calculated eigenvector force. An orthogonal matrix that also has eigenvector forces is an orthogonal matrix that has eigenvectors as columns. Details of the matrix elements will be described later.

[0030] The eigenvalue decomposition storage unit 16 stores the eigenvalues subjected to the eigenvalue decomposition by the eigenvalue decomposition unit 15 and the matrix elements described above. The eigenvalue decomposition storage unit 16 also stores eigenvalues and matrix elements calculated by the eigenvalue calculation unit 17. The eigenvalue decomposition storage unit 16 can be realized by a predetermined recording medium (for example, a semiconductor memory, a magnetic disk, an optical disk, etc.). The storage in the eigenvalue decomposition storage unit 16 may be temporary storage in RAM or the like, or may be long-term storage.

[0031] The eigenvalue calculation unit 17 reads the eigenvalues calculated by the eigenvalue decomposition unit 15 and the matrix elements from the eigenvalue decomposition storage unit 16, and uses the eigenvalues and matrix elements to calculate the symmetric tridiagonal matrix of the division source. The eigenvalue and the matrix element of the symmetric tridiagonal matrix of the division source are calculated and stored in the eigenvalue decomposition storage unit 16. The eigenvalue calculation unit 17 recursively repeats the process of calculating the eigenvalues and matrix elements of the symmetric tridiagonal matrix of the division source until the eigenvalues of the symmetric tridiagonal matrix T are calculated. Then, the eigenvalue calculation unit 17 stores the calculated eigenvalues of the symmetric tridiagonal rows in the eigenvalue storage unit 18.

[0032] In the eigenvalue storage unit 18, eigenvalues of the symmetric tridiagonal matrix T are stored. The eigenvalue storage unit 18 can be realized by a predetermined recording medium (for example, a semiconductor memory, a magnetic disk, an optical disk, etc.). Storage in the eigenvalue storage unit 18 may be temporary storage in RAM or the like, or may be long-term storage.

The eigenvector calculation unit 19 reads the symmetric tridiagonal matrix T from the diagonal matrix storage unit 13 and reads the eigenvalues of the symmetric tridiagonal matrix T from the eigenvalue storage unit 18. Then, the eigenvector calculation unit 19 calculates the eigenvectors of the symmetric tridiagonal matrix T from the symmetric tridiagonal matrix T and its eigenvalues using the twist decomposition method, and stores them in the eigenvector storage unit 20. . The eigenvector calculation unit 19 may calculate the eigenvector by the LV type twist decomposition method, or may calculate the eigenvector by the qd type twist decomposition method. A specific twist decomposition method will be described later. The eigenvector storage unit 20 stores eigenvectors of the symmetric tridiagonal matrix T. The eigenvector storage unit 20 can be realized by a predetermined recording medium (for example, a semiconductor memory, a magnetic disk, an optical disk, etc.). Storage in the eigenvector storage unit 20 may be temporary storage in RAM or the like, or may be long-term storage.

Note that any two or more storage units of the matrix storage unit 11, the diagonal matrix storage unit 13, the eigenvalue decomposition storage unit 16, the eigenvalue storage unit 18, and the eigenvector storage unit 20 are realized by the same recording medium. Alternatively, it may be realized by a separate recording medium. In the former case, for example, the area where the symmetric matrix A is stored becomes the matrix storage unit 11, and the area where the symmetric tridiagonal matrix T is stored becomes the diagonal matrix storage unit 13.

[0036] Each storage unit of the matrix storage unit 11, the diagonal matrix storage unit 13, the eigenvalue decomposition storage unit 16, the eigenvalue storage unit 18, and the eigenvector storage unit 20 may be composed of two or more recording media.

[0037] Next, the operation of the eigenvalue decomposition apparatus 1 according to the present embodiment will be described using the flowchart of FIG.

(Step S101) The diagonalization unit 12 reads the symmetric matrix A stored in the matrix storage unit 11! And calculates the symmetric tridiagonal matrix T by tridiagonalizing the matrix A. And stored in the diagonal matrix storage unit 13.

(Step S102) The eigenvalues of the symmetric tridiagonal matrix T are calculated by the matrix dividing unit 14, the eigenvalue decomposing unit 15, and the eigenvalue calculating unit 17 and stored in the eigenvalue storage unit 18. Details of this processing will be described later.

(Step S103) The eigenvector calculation unit 19 reads the symmetric tridiagonal matrix T from the diagonal matrix storage unit 13, reads the eigenvalues of the symmetric tridiagonal matrix T from the eigenvalue storage unit 18, The eigenvectors of the diagonal matrix T are calculated and stored in the eigenvector storage unit 20. Details of this processing will be described later.

[0040] In this way, the eigenvalue decomposition of the symmetric tridiagonal matrix T is completed. Here, since the eigenvalues of the symmetric matrix A are equal to those of the symmetric tridiagonal matrix, the eigenvalues of the symmetric matrix A are also calculated. In addition, as will be described later, the eigenvector of the symmetric matrix A can be easily calculated by performing a predetermined transformation, as well as the eigenvector force of the symmetric tridiagonal matrix T. The

Next, the processing in steps S101 and S102 in the flowchart of FIG. 2, that is, the processing until the eigenvalue of the symmetric tridiagonal matrix T is calculated will be described in more detail. First, the processing in step S102 of the flowchart of FIG. 2 will be described using the flowchart of FIG.

[0042] (Step S201) The matrix dividing unit 14 reads the symmetric tridiagonal matrix T from the diagonal matrix storage unit 13, and divides the symmetric tridiagonal matrix T into two symmetric tridiagonal matrices. And stored in the diagonal matrix storage unit 13. The matrix partitioning unit 14 performs the process of dividing the symmetric tridiagonal matrix into two symmetric tridiagonal matrices until each symmetric tridiagonal matrix after the division becomes a predetermined size or less. repeat. Details of this processing will be described later.

(Step S 202) The eigenvalue decomposition unit 15 performs eigenvalue decomposition on each symmetric tridiagonal matrix having a size equal to or less than a predetermined size stored in the diagonal matrix storage unit 13. The eigenvalue decomposition unit 15 accumulates eigenvalues as a result of eigenvalue decomposition and matrix elements that are partial elements of an orthogonal matrix such as eigenvectors in the eigenvalue decomposition storage unit 16.

[0044] (Step S203) The eigenvalue calculation unit 17 reads the eigenvalues and matrix elements of the symmetric tridiagonal matrix from the eigenvalue decomposition storage unit 16, and the eigenvalues and divides of the symmetric tridiagonal matrix of the division source. The matrix elements of the original symmetric tridiagonal matrix are calculated and stored in the eigenvalue decomposition storage unit 16. The eigenvalue calculation unit 17 repeats the process of calculating the eigenvalues and matrix elements of the symmetric tridiagonal matrix of the division source until the eigenvalues of the symmetric tridiagonal matrix T are calculated, and the symmetric tridiagonal matrix is calculated. The eigenvalues of the matrix T are stored in the eigenvalue storage unit 18. In this way, the process for calculating the eigenvalue is completed.

[0045] [Diagonalization of symmetric matrix A]

It is known that the symmetric matrix A can be tridiagonalized, as shown below, using, for example, No, Usholder transform, etc. Here, the symmetric tridiagonal matrix T is a square matrix in which the number of rows matches the number of columns.

[0046] [Equation 1]

No

Where «x« real symmetric matrix and Ρ is an orthogonal matrix.

Accordingly, the diagonalization unit 12 can read the symmetric matrix から from the matrix storage unit 11 and calculate the symmetric tridiagonal matrix T as described above. The calculated symmetric tridiagonal matrix T is stored in the diagonal matrix storage unit 13 (step S101).

[0048] [Partition of symmetric tridiagonal matrix T]

The process of dividing the symmetric tridiagonal matrix T will be described. N as shown in the following equation

Given a symmetric tridiagonal matrix T of Xn, they are divided into two symmetric tridiagonal matrices T

, T and other elements (Θ and β below).

[0049] [Equation 2]

Where Γ is a real symmetric tridiagonal matrix of «x« («is an integer greater than 2), 7; is a real symmetric tridiagonal matrix of f ^ x, and Γ ₂ is (<<-? ¾) of X Real symmetric tridiagonal matrix,

e _; ≡ (0, ..., 0,1) ^τ eR "\ e _f ≡ (1,0, ..., 0) ^τ e R""\ β is (? ¾ +1,? ¾ of Γ ) Elements,

Θ is a non-zero constant.

Here, η is an integer such that Κη く η. For example, the specific process of dividing an 8 x 8 symmetric tridiagonal matrix Τ into half the size is as follows.

[Equation 3]

k

h

t

T =, in

^ 15

[0051] However, β = t, and Θ can be arbitrarily set to reduce the calculation error.

8

The In each of the above matrices, the elements of the matrix other than those specified are 0. It should be noted that each element specified in each matrix may be non-zero or zero. Here, in order to appropriately execute parallel processing, it is preferable to set η to the maximum integer that does not exceed ηΖ2, or the minimum integer that does not fall below ηΖ2 (the value of η is Using this method to divide a matrix into two matrices is called “dividing a matrix into two substantially half matrices”).

[0052] [Equation 4] ί¾ = L "/ 2"

In the present embodiment, is given by the above equation. Needless to say, the value of is arbitrary within the range of l <n <n, as described above.

Therefore, the matrix dividing unit 14 reads the symmetric tridiagonal matrix T stored in the diagonal matrix storage unit 13 as described above, and converts the symmetric tridiagonal matrix T into two pieces. Can be divided into symmetric tridiagonal matrices と, T and two values | 8, Θ, and repeat the division process

1 2

Can do. The matrix partitioning unit 14 is divided into two symmetric tridiagonal matrices Τ, T and two values | 8

1 2

, Θ are stored in the diagonal matrix storage unit 13. It is preferable that the two values Θ are stored in the diagonal matrix storage unit 13 so that it can be determined which division the value is generated in. FIG. 4 is a flowchart showing the process of dividing the matrix by the matrix dividing unit 14 in step S201 of the flowchart of FIG.

(Step S301) The matrix dividing unit 14 sets the counter I to “1”.

(Step S302) The matrix partitioning unit 14 reads a symmetric tridiagonal matrix that has not been subjected to the I-th partitioning from the diagonal matrix storage unit 13, and converts the symmetric tridiagonal matrix into two symmetric tridiagonals. Divide into rows and other elements. The matrix dividing unit 14 stores the two divided tridiagonal matrices and other elements in the diagonal matrix storage unit 13.

(Step S 303) The matrix division unit 14 performs the I-th division and determines whether or not the symmetric tridiagonal matrix force diagonal matrix storage unit 13 stores them. Then, if the symmetric tridiagonal matrix is stored in the diagonal matrix storage unit 13 after performing the I-th division, the process returns to step S302. Proceed to step S304.

[0056] (Step S304) The matrix partitioning unit 14 determines whether the symmetric tridiagonal matrix that has undergone the I-th partition is a force that is less than or equal to a predetermined size. For example, the matrix dividing unit 14 reads the size of the target matrix from a recording medium (not shown) that stores the size of the target matrix (for example, 25 X 25). It may be determined whether or not the I-th divided symmetric tridiagonal matrix stored in the diagonal matrix storage unit 13 is less than or equal to the size thereof. If the size of the symmetric tridiagonal matrix that has undergone the I-th division is less than or equal to a predetermined size, the process of dividing the symmetric tridiagonal matrix ends, otherwise it is not. If yes, go to Step S305. (Step S305) The matrix dividing unit 14 increments the counter I by 1. Then, the process returns to step S302.

In the flowchart of FIG. 4, as described above, the matrix dividing unit 14 divides each matrix into two substantially half matrices, and thus each symmetric tridiagonal after the I-th division. The case where the matrix sizes are almost the same has been described. However, if the matrix dividing unit 14 does not divide each matrix into two substantially half matrices, the matrix dividing unit 14 The division is repeated so that the matrix is less than or equal to the predetermined size.

In the flowchart of FIG. 4, in step S304, the matrix dividing unit 14 compares the size of the matrix after division stored in the diagonal matrix storage unit 13 with a predetermined size. However, this is only an example, and the matrix divider 14 may perform other processes in step S304. For example, as described above, when the matrix dividing unit 14 divides each matrix into two substantially half matrices, as long as the size of the original symmetric tridiagonal matrix T can be known, what number It is possible to know whether the size of the target matrix is obtained by dividing. Therefore, if the target matrix size is obtained in the Nth division (N is an integer equal to or greater than 1), whether or not I is N is compared in step S304. If it is N, a series of processes may be terminated.

FIG. 5 is a diagram for explaining matrix division. First, the matrix dividing unit 14 divides a symmetric tridiagonal matrix T into a symmetric tridiagonal matrix T and a symmetric tridiagonal matrix T as the first division (steps S301 and S302). . There is only one symmetric tridiagonal matrix T

2

Therefore, the matrix dividing unit 14 determines that there is no matrix after the first division (step S303). Also, assuming that the symmetric tridiagonal matrix T etc. is not a matrix smaller than a predetermined size (step S304), the matrix dividing unit 14 uses the symmetric tridiagonal matrix T as the second division. Divide into tridiagonal matrix T and symmetric tridiagonal matrix T (step S305,

1 11 12

S302). In this case, there is a symmetric tridiagonal matrix T that does not perform the second partition.

Therefore, the matrix partitioning unit 14 has a symmetric tridiagonal matrix T that is also symmetric with the symmetric tridiagonal matrix T.

2. Divide into 21-diagonal matrix T (steps S303, S302). In this way, each divided

twenty two

The process of dividing the matrix is repeated until the symmetric tridiagonal matrix is less than or equal to the target size. Is done. In FIG. 5, elements other than the symmetric tridiagonal matrix (13, Θ described above) are omitted. In FIG. 5, each symmetric tridiagonal matrix is a square matrix.

[0060] [Eigenvalue decomposition of divided matrix]

Assuming that the matrix T is a symmetric matrix, the eigenvalue decomposition of the matrix T is as follows.

[Equation 5]

[0061] Here, D is a diagonal matrix having the same size as the matrix T. Each diagonal component of D is an eigenvalue of T. Q is an orthogonal matrix that also contains the eigenvector force of the matrix T.

Therefore, the eigenvalue decomposition unit 15 reads each symmetric tridiagonal matrix after the division from the diagonal matrix storage unit 13 and performs eigenvalue decomposition as described above (step S 202). As described above, the eigenvalue decomposition method, for example, the method combining the bisection method and the inverse iteration method, the MR "3 method, the QR method, etc., may be used. When the multi-diagonal matrix is divided, the eigenvalue decomposition unit 15 performs each symmetric tri-diagonal matrix Τ, T,

111 112

Eigenvalue decomposition is performed on Τ, Τ,. The eigenvalue decomposition unit 15

121 122 222

Stores the eigenvalues obtained as a result of the eigenvalue decomposition and matrix elements which are part of the orthogonal matrix obtained as a result of the eigenvalue decomposition in the eigenvalue decomposition storage unit 16. Here, which elements of the matrix element orthogonal matrix are included will be described later. The eigenvalue decomposition unit 15 may store eigenvalues and matrix elements in the eigenvalue decomposition storage unit 16 instead of storing the eigenvalues and matrix elements in the eigenvalue decomposition storage unit 16. Even in this case, at least matrix elements are accumulated by accumulating orthogonal matrices that also have eigenvector forces.

[0063] As will be described later, what is necessary to use when calculating eigenvalues is the first row and the last row of the orthogonal matrix Q, not the entire orthogonal matrix Q. Therefore, in this eigenvalue decomposition, the first row and the last row of the orthogonal matrix Q may be obtained instead of obtaining the entire orthogonal matrix Q. To that end, for example, provided in LAPACK3.0! Instead of using the unit matrix I as the initial value in the calculation of the orthogonal matrix Q, use the matrix 示 shown in the following equation in DSTEQR, the eigenvalue decomposition routine of the QR method. [Equation 6]

Where ρ, is an orthogonal matrix of η χ η; Τ is a matrix of 2 χ η.

[0064] By doing this, it is possible to reduce the calculation area for obtaining the first and last rows of the orthogonal matrix Q by DSTEQR. When calculating the orthogonal matrix Q itself, a work area of 0 (η ² ) is required, but when calculating the first and last rows of the orthogonal matrix Q using the above matrix Υ, A work area of Ο (η) is sufficient. As a result, the amount of calculation is reduced and the processing speed of eigenvalue decomposition is improved. As described above, the eigenvalue decomposition unit 15 may perform eigenvalue decomposition on the divided tridiagonal matrix after the division to calculate the orthogonal matrix Q itself, or a part of the orthogonal matrix (here, Then, the process of calculating the first and last lines) may be performed. This latter process is also called eigenvalue decomposition. In this way, the eigenvalue decomposition unit 15 performs eigenvalue decomposition on each symmetric tridiagonal matrix with a predetermined size or less, and the eigenvalues of each symmetric tridiagonal matrix and each symmetric tridiagonal pair. Even if at least the matrix elements (here, the first row and the last row) of the orthogonal matrix composed of the eigenvectors of the angular matrix are calculated and the unique values and the matrix elements are stored in the eigenvalue decomposition storage unit 16 Good.

[0065] [Calculation of eigenvalues]

First, as shown in Fig. 6, from two divided tridiagonal matrices Τ and Τ,

1 2

The process of calculating the eigenvalues of the original symmetric tridiagonal matrix Τ will be described. Symmetrical triple

0

Suppose that the diagonal matrices Τ and Τ are eigenvalue decomposed as follows. Where symmetric tridiagonal

1 2

The matrices Τ and Τ are both square matrices.

1 2

[Equation 7]

[0066] In this case, the symmetric tridiagonal matrix T of the division source is as follows.

0

[Equation 8]

Where is the last row of, and f ₂ is the first row of ρ ₂ .

In the above equation, j8, 0 are the values described in the process of dividing the symmetric tridiagonal matrix. If the eigenvalue decomposition of the matrix (D + p ζζ ^τ ) is obtained, the eigenvalues of the matrix Τ

12 0 Decomposition was obtained. Because the matrix (D + ρ ζζ ^Τ )

Eigenvalue decomposition of 12 D + ρ ζζ ^Τ

12

= When QDQ ^T is obtained, the following equation is obtained, and eigenvalue decomposition of the matrix T is completed.

Next, a method for calculating the eigenvalues of the matrix (D + p _ZZ ^T ) will be described. Matrix (D

+ ρ ζζ) and the eigenvector is q.

[Equation 10]

(+ pzz ^T ) q, = to q,

It becomes.

[0069] (l) When the component of z does not contain 0 and the diagonal component of D does not contain the same value

12

First, if the z component does not contain 0 and the diagonal component of D does not contain the same value,

12

I will explain. In that case, the above equation can be modified sequentially as shown below.

[Equation 11]

1 + ρ z ^T (D ₁₂ — λ, · /)-= 0

[0070] where

[Equation 12] D ₁₂ = dia ic / j, d ₂ , ---, d _n )

z = (ζ ₁ , ζ ₂ , ..., ζ _π ) ^τ

Then, the following eigen equation can be obtained. The eigenvalues ^ = 1, 2,..., N) are the solutions of the eigen equation n. The following f (λ) is monotonically increasing, and since the range of solutions of f (λ) = 0 is known, it can be calculated using the Newton method.

[Equation 13] (λ) = ι _{+ Ρ} = 0 (Formula

Note that D = diag (, λ,..., Λ). The matrix Q is an eigenvalue of the matrix (D + P ZZ).

1 2 n 12

Q = (q, q,. Therefore,

1 2 n

Finding the tuttle q gives the matrix Q. The eigenvector q. Is given by the following equation.

[Equation 14]

q '(Formula 2

[0072] (2) When 0 is included in the component of z

When 0 is included in the component of z, the order of the matrix can be lowered, and the eigenvalue and the eigenvector can be obtained using the matrix with the lowered order. The operation of lowering the order of this matrix is called deflation.

[0073] When z = 0, (d, e) is a set of eigenvalues and eigenvectors as shown in the following equation. However, e is a unit vector in which the j-th component is 1 and the other components are 0.

[Equation 15]

^{= ((+ Ρζζ Ί) -d} J I) e J (D 12 -d J I) e J + pz (z T)

= 0

[0074] Further, specific equation of the matrix (D + ρ ζζ ¹⁾ is expressed by the following equation _c

12

[Equation 16] \ p _u + pzz ^T -I _n \ = (j-λ _; ) | θ ₁₂ + pzz T-where f) _u + pii ^T excludes the j th row and j th column of D _u + pzz ^τ It is a small matrix. / „Is an nth-order unit matrix.

[0075] Therefore, the eigenvalue when z = 0 is d and the j-th row and j-th column of the matrix (D + p ζζ ^τ ) jj 12

This is the value obtained by solving the eigen equation of Equation 1 above for the submatrix. That is, the eigenvalue is d and the value obtained by solving the following eigen equation. The eigenvalues obtained by solving the following eigenequation are assumed to be eigenvalues (i = l, 2, ···, j-1, j + 1, ···, n).

[Equation 17]

: 0

The eigenvector is the small matrix excluding the j-th row and j-th column of the matrix (D + p zz ¹ ).

12

The above-mentioned e is the n-dimensional vector of the following formula with 0 inserted between the (j 1) -th row of the n−l-dimensional eigenvector calculated by Equation 2 and the j-th row.

[Equation 18] z ₁ ^∑ 2

d,-λ, d ― Λ, — to d -λ, 'd „-λ,

[0077] Thus, when the component of z includes 0, the dimension is smaller than that of the matrix (D + p ζζ ^τ ).

12

Since eigenvalue decomposition is performed for a large matrix, the calculation can be performed at high speed, and as a result, the calculation time can be shortened.

[0078] (3) When the diagonal values of D contain the same value

12

If d = d 2 ... = d, the same deflation as described above is performed. j + 1 j + 2 k

Thus, eigenvalues and eigenvectors can be calculated.

[0079] [Equation 19] τ = ^ j + i, ^z j + i 7…, ^z t)

In contrast, consider the following two orthogonal matrices. Q ζ = (*, 0,0, ..., 0) ^τ

Where Qe R ()) _c

* Is a predetermined real number.

[0080] Then, from d: d,

j

[Number 21]

QD ₁₂ Q ^T = D ₁₂ Therefore, the following equation is obtained.

[Number 22]

[0081] where

[Equation 23]

D --— · X

1 ₂ + pzz

The eigenvalues and eigenvectors of can be obtained in the same manner as in the case of deflation where the z component contains 0. Also,

[Number 24]

Q

Is uniquely determined except for the arbitraryness of the code. Here, the calculation method will be described. First, consider an orthogonal matrix that satisfies the following equation.

[0082] [Equation 25]

Q, '(+ I…,,…,) — ⁼ (*, ^z _{J + 2} ,… _; “O w,…,)

If this orthogonal matrix is taken as the following equation, the above equation is satisfied.

[Equation 26]

[0083] Here, the value of * in the above equation is obtained from the following equation:

[Equation 27]

Q _{l +} + ι, ί ■ j, ++ ι ¾)

[0084] Therefore,

[Equation 28]

Q i ^ j 'Q

It becomes. as a result,

[Equation 29]

Q

Can also be obtained by removing the arbitraryness of the code. In this way, d = d = ~ = d

Even in the case of j + 1j + 2k, the eigenvalues and eigenvectors of the matrix (D— + pzz " ¹ ) can be obtained.

12

In other words, the eigenvalues are d 1, d 2,..., D and values obtained by solving the following eigen equations. Next j + 2 j + 3

Let the eigenvalues obtained by solving the eigenequation of Eigenvalues (i-, 2, ..., j + k- ·, n).

[Equation 30] (λ) = ι + Ρ ∑ -Γ ^ = °

m-1 —People The eigenvector is

[Equation 31] (ρ ^τ + l≤i≤n)

Le u + 2≤i≤k

[0086] When (2) z component includes 0, (3) D diagonal component includes the same value.

12

The same value for the diagonal component of D

12

The determination of whether or not it is included may be made in consideration of a minute error. For example, if the following equation is satisfied, it may be determined that the component of z is 0, and the diagonal component of D is

12

You can judge that they are the same value!

[Equation 32]

RQLj

rQL ≡ maxf <i ₁ , (i ₂ ,---, <i.-) xe? iiD.

[0087] Note that ε is a value obtained by multiplying a machine epsilon by a constant (for example, a real number such as 2, 4, 8, and the like, which is preferably selected to reduce a calculation error). Other suitable values may be used.

[0088] In this way, (1) 0 is not included in the component of ζ, and the same value is not included in the diagonal component of D

12

(2) When the component of ζ contains 0, (3) When the diagonal component of D contains the same value

12

For each of these, the eigenvalue decomposition D + p ζζ ^τ 2 QDQ ^T is obtained. Gatsutsu

12

Since the matrix Q can be obtained from the matrix Q and the matrix Q, and the matrix D is also obtained,

12 0

Eigenvalue decomposition of the symmetric tridiagonal matrix T of the quotient can be performed. Thus, the target

0

The method of dividing a matrix into two sub-matrices, performing eigenvalue decomposition on each sub-matrix after the division, and then performing eigenvalue decomposition on the original matrix using eigen-value decomposition of each sub-matrix is as follows: It is called Divide and Conquer (D & C).

In the divide-and-conquer method, even eigenvectors are calculated, and in the process of calculating the eigenvectors, a matrix must be calculated. Therefore, the processing is very heavy. When performing eigenvalue decomposition in the divide-and-conquer method, for example, about 95% of the computation time is a vector update process that calculates the eigenvector (for example, as described above, the matrix Q and The process of multiplying the matrix to calculate the matrix Q 0 from the matrix Q 12 may be spent.

[0090] On the other hand, when calculating only the eigenvalues of the matrix T to be divided, the above process is simplified.

0

Can. Next, the method will be described. From the above results, it is as follows.

[0091] [Numerical 33]

f = f ₁₂ 2 (Equation 3)

Ι = Ι ₁₂ β (Equation 4)

Where f is ρ. The first row of I and I is ρ. Is the last line.

F ₁₂ is the first line of g ₁₂ and f ₁₂ =, (), ···, 0).

Ι ₁₂ is the last line of ρ ₁₂ , and ι ₁₂ = (ο,-·-, ο, ι ₂ ).

F \ is the first line of ¾, and ι ₂ is the last line of ρ ₂ .

[0092] Here, as described above, f is the first row of the orthogonal matrix resulting from eigenvalue decomposition of the matrix Τ.

0

Yes, 1 is the last row of the orthogonal matrix resulting from eigenvalue decomposition of the matrix T. This f and 1, sand

0

That is, the first row and the last row of the orthogonal matrix are matrix elements. For the orthogonal matrix Q 12, f, 1

12 12 are the first row and the last row of the orthogonal matrix Q after the columns are replaced.

12

[0093] From the above results, as shown in Fig. 6, the matrix T of the division source is constructed from the matrices, and T.

1 2 0, for matrix τ,

0

[Equation 34]

λ f, i

Can be calculated respectively. By repeating such processing, the eigenvalues of the symmetric tridiagonal matrix T obtained by tridiagonalizing the symmetric matrix A can be finally calculated.

[0094] When the same value is included in the diagonal component of (3) D above, Equation 3 and Equation 4 are

12

I will become.

[Equation 35] f = f ₁₂ Q = f ₁₂ Q ^T Q

However, = QQ.

[0095] Here, when calculating the above formula, [Equation 36]

Is calculated first, then

[Equation 37]

Q

Since the amount of calculation is reduced by multiplying, it may be calculated in that order.

Therefore, the eigenvalue calculation unit 17 reads the eigenvalues and matrix elements from the eigenvalue decomposition storage unit 16 as described above, and the other values generated when the matrix is divided from the diagonal matrix storage unit 13. By reading β and Θ related to the elements and using them, the eigenvalues and matrix elements of the symmetric tridiagonal matrix of the division source are calculated. Then, by repeating the process of calculating them, the eigenvalues of the symmetric tridiagonal matrix Τ are finally calculated. FIG. 7 is a flowchart showing a process in which the eigenvalue calculation unit 17 calculates the eigenvalue in step S203 of the flowchart of FIG.

(Step S401) The eigenvalue calculation unit 17 sets the counter J to “1”.

(Step S402) The eigenvalue calculation unit 17 determines whether the calculation of the Jth eigenvalue is the calculation of the last eigenvalue. Here, the calculation of the last eigenvalue is to calculate the eigenvalue of the symmetric tridiagonal matrix T. If the last eigenvalue is calculated, the process proceeds to step S406. If not, the process proceeds to step S403.

[0098] (Step S403) The eigenvalue calculation unit 17 calculates eigenvalues and the like of the symmetric tridiagonal matrix of the division source. Details of this processing will be described later.

(Step S404) The eigenvalue calculation unit 17 determines whether eigenvalues and the like of all the symmetric tridiagonal matrices of the division source have been calculated in the calculation of the Jth eigenvalue. Then, in calculating the J-th eigenvalue, if eigenvalues and the like of all the symmetric tridiagonal matrices of the division source are calculated, the process proceeds to step S405, and if not, the process returns to step S403.

(Step S 405) The eigenvalue calculation unit 17 increments the counter J by 1. Then, the process returns to step S402.

(Step S406) The eigenvalue calculation unit 17 calculates the eigenvalues of the symmetric tridiagonal matrix T and stores the calculated eigenvalues in the eigenvalue storage unit 18. In this way, a symmetric tridiagonal matrix T The series of processes for calculating the eigenvalues of is finished.

[0100] FIG. 8 is a flowchart showing the detailed processing of step S403 in the flowchart of FIG. Before executing the processing of the flowchart of FIG.

It is preferable to rearrange the 12 diagonal components in ascending order or descending order. In this case, the columns of the orthogonal matrix Q and the components of z are the same as the rearrangement of the diagonal matrix D.

12 12

Shall be rearranged. Specifically, the i-th component and j-th component of the diagonal matrix D are input.

12

In this case, the 洌 and j-th columns of the orthogonal matrix Q are swapped, and the z-th and z-th rows of z

12

Can be replaced.

[0101] (Step S501) The eigenvalue calculation unit 17 calculates the eigenvalue of the matrix of the division source using Equation 1 or the like. Then, the eigenvalue calculation unit 17 accumulates the calculated eigenvalues by 16 eigenvalue decomposition storage units.

[0102] (Step S502) The eigenvalue calculating unit 17 calculates q. Of the equal force of equation 2 using the eigenvalue calculated in step S501. By calculating q., The matrix Q having q. In the column vector is calculated. The eigenvalue calculation unit 17 may primarily store the calculated matrix Q on a recording medium (not shown).

(Step S503) The eigenvalue calculation unit 17 calculates the matrix elements f and 1 related to the matrix of the division source, using the matrix Q calculated in step S502. These matrix elements are shown in Equation 3 and Equation 4. Then, the eigenvalue calculation unit 17 accumulates the calculated matrix elements of the division source in the eigenvalue decomposition storage unit 16. In this way, the process of step S403 ends.

[0104] In the flowchart of Fig. 8, the z component does not include 0, and the diagonal component of D

12

If the same value is not included, the matrix Q consisting of eigenvalues and eigenvector forces can be calculated as described in (1) above. Also, the z component contains 0, and the diagonal component of D

12

If the same value is not included in, the matrix Q with eigenvalues and eigenbeta tokens can be calculated as described in (2) above. Also, the same value is included in the diagonal component of D

12

In this case, the matrix Q including eigenvalues and eigenvector forces can be calculated as described in (3) above regardless of whether or not 0 is included in the z component. When calculating eigenvalues, etc. as described in (3) above, as described above, the matrix Q, which is not the matrix Q, is directly applied to f etc. I will calculate 1 You may do it.

In the flowchart of FIG. 8, the force matrix elements f and 1 shown for the case where the matrix elements f and 1 are calculated at once from the matrix Q may be calculated for each component. Eq. 3, Eq. 4 Force As can be seen, f, 1 multiplied by eigenvector q Force Matrix element f, i-th of 1

12 12 i

Become an ingredient. When the matrix elements f and 1 are calculated in this way, the calculation can be performed by storing only the eigenvector eigenvector q while changing the value of i sequentially. Therefore, it is not necessary to store the entire matrix Q in the calculation of the matrix elements f and 1, and the work area on the memory used for the calculation can be reduced.

[0106] FIG. 9 is a diagram for explaining processing for calculating eigenvalues. In the eigenvalue decomposition storage unit 16, the eigenvalue decomposition unit 15 performs eigenvalue decomposition on the matrices Τ, Τ, ...

111 112 222

As a result, that is, those eigenvalues and matrix elements are stored. First, the eigenvalue calculation unit 17 starts the first calculation of eigenvalues and the like (steps S401 and S402). Here, the first calculation of eigenvalues, etc. is to calculate the eigenvalues and matrix elements of the matrix in the second row from the bottom from the eigenvalues and matrix elements of the bottom row of FIG. That is. The eigenvalue calculation unit 17 calculates the eigenvalues of the matrix Τ and the matrix Τ, the matrix Τ, and the matrix Τ

The matrix elements 111 112 111 11 are read from the eigenvalue decomposition storage unit 16. In addition, the matrix を and the matrix Τ

2 11 111

2 values 13 and Θ generated by decomposing into matrix Τ are read from diagonal matrix storage unit 13.

112

Stick out. Then, the eigenvalue calculation unit 17 uses these values to calculate the eigenvalues of the matrix Τ.

11

And stored in the eigenvalue decomposition storage unit 16 (step S501). Next, the eigenvalue calculation unit 17 calculates an orthogonal matrix Q using the eigenvalues of the matrix T (step S502), and the orthogonal matrix

11

Calculate the matrix elements of the matrix T using Q and store them in the eigenvalue decomposition storage unit 16 (steps

11

S503).

[0107] Since the first calculation of eigenvalues and the like is not yet completed (step S404), eigenvalue calculation unit 17 is based on eigenvalues, matrix elements, and the like of matrix T and matrix T in the same manner as described above.

121 122

The eigenvalues and matrix elements of the matrix T are calculated (steps S402 and S403). Like this

12

When the first calculation of eigenvalues, etc. is completed, the eigenvalue calculation unit 17 performs the second calculation of the eigenvalues, i.e., the matrix T, the eigenvalues of the matrix T, the matrix elements, etc.

11 12 1 Eigenvalues and matrix elements are calculated (steps S404, S405, S402, S403). Then fix Similarly, the value calculation unit 17 is based on the eigenvalues, matrix elements, etc. of the matrix T and the matrix Τ.

21 22

Then, eigenvalues and matrix elements of the matrix Τ are calculated (steps S404 and S403).

2

[0108] When the eigenvalue calculation unit 17 finishes the second calculation of eigenvalues, etc. (step S404), the next eigenvalue calculation is the final process for obtaining eigenvalues of the symmetric tridiagonal matrix T (step S404). (S405, S402), the eigenvalue calculation unit 17 only calculates the eigenvalue, and stores the calculated eigenvalue in the eigenvalue storage unit 18 (step S406). In this way, the calculation of the eigenvalue is completed.

Note that the eigenvalues of the matrix T are calculated based on the matrix T, the eigenvalues of the matrix T, the matrix elements, etc. in FIG.

1 2

When calculating, it is only necessary to use the matrix element f 1, or the matrix element f 1. Shi

1 2 2 1

Therefore, when calculating the matrix element of matrix T, matrix T, matrix element f, 1, or

1 2 1 2 matrix element f, 1

2 1 need not be calculated. By doing so, the amount of calculation can be reduced and the processing speed can be improved.

Next, the eigenvector calculation unit 19 will be described. In the present embodiment, a case where the eigenvector calculation unit 19 calculates (1) an eigenvector by the qd type twist decomposition method and (2) a case where an eigenvector is calculated by the LV type twist decomposition method will be described.

[0111] (l) When calculating eigenvectors by qd-type twist decomposition

FIG. 10 is a diagram showing a configuration of the eigenvector calculation unit 19 when the eigenvector is calculated by the qd-type twist decomposition method. In FIG. 10, the eigenvector calculation unit 19 includes a positive definite value conversion unit 21, a Cholesky decomposition unit 22, and a vector calculation unit 23.

[0112] The positive definite value i section 21 reads the eigenvalues of the symmetric tridiagonal matrix T from the eigenvalue storage section 18, and determines whether any of the eigenvalues is a negative value. Then, when any of the eigenvalues is negative, the positive definite part 21 reads the symmetric tridiagonal matrix T from the diagonal matrix storage unit 13, and the symmetric tridiagonal matrix T Without changing the eigenvector of, the symmetric tridiagonal matrix T is positive definite and the positive definite Then, the eigenvalues are stored in the eigenvalue storage unit 18, and the symmetric tridiagonal matrix T after the positive definite value is stored in the diagonal matrix storage unit 13. This positive definite processing is performed because all eigenvalues of the symmetric tridiagonal matrix T must be positive when eigenvectors are calculated by the qd-type twist decomposition method. Symmetry 3 A matrix with a positive definite matrix を is also called a symmetric tridiagonal matrix Τ. In addition, when this positive definite value 匕 processing is performed, the Cholesky decomposition and eigenvector calculation processing are performed on the symmetric tridiagonal matrix Τ and eigenvalue after the positive definite value. Become. Details of this processing will be described later. Here, even when positive definite value 匕 is performed, the eigenvalue storage unit 18 stores the eigenvalue before the positive definite value 匕. The eigenvalues before the definiteization are eigenvalues finally obtained by the eigenvalue decomposition apparatus 1. When the positive definite value 匕 is performed and only the eigenvector is calculated by the qd-type twist decomposition method, the symmetric tridiagonal matrix T stored in the diagonal matrix storage unit 13 is used. You can overwrite it with a symmetric tridiagonal matrix T after the positive definite 匕. On the other hand, when eigenvectors are also calculated by the LV-type twist decomposition method, the diagonal matrix storage unit 13 stores the symmetric tridiagonal matrix T before positive definite 匕 and the symmetric 3 after positive definite conversion. The multi-diagonal matrix T is stored.

[0113] The Cholesky decomposition unit 22 reads a symmetric tridiagonal matrix T in which all eigenvalues are positive values from the diagonal matrix storage unit 13, and the symmetric tridiagonal matrix T from the eigenvalue storage unit 18. Read the eigenvalue of the positive value of. The Cholesky decomposition unit 22 performs a qd-type transformation on each element of the symmetric tridiagonal matrix T, thereby transforming the symmetric tridiagonal matrix T into the upper bidiagonal matrix B ⁽ ⁺¹⁾ and the lower 2 Cholesky decomposition into a bidiagonal matrix B ( ^_1) . Details of this process will be described later. Note that “the Cholesky decomposition of a symmetric tridiagonal matrix T into an upper bidiagonal matrix B ( ⁺¹⁾ and a lower bidiagonal matrix B ^(_ ” ⁾ means that the symmetric tridiagonal matrix T Cholesky decomposition of the matrix obtained by multiplying the eigenvalue of x by the unit vector from the symmetric tridiagonal matrix T into the upper bidiagonal matrix B ( ⁺¹⁾ and the lower bidiagonal matrix B ^(_1) May be included.

[0114] The vector calculation unit 23 includes each element of the upper double diagonal matrix B ⁽⁺¹⁾ and the lower double diagonal matrix B ( ^_1) calculated by the Cholesky decomposition unit 22, and a symmetric tridiagonal matrix. The eigenvector is calculated using the positive eigenvalue of T and stored in the eigenvector storage unit 20. A matrix having this eigenvector in each column is an orthogonal matrix for eigenvalue decomposition.

Next, the process of step S103 in the flowchart of FIG. 2, that is, the process of calculating eigenvectors of the symmetric tridiagonal matrix T will be described in more detail. First, the processing in step S 103 of the flowchart of FIG. 2 will be described with reference to the flowchart of FIG. [0116] (Step S601) The positive definite key section 21 reads the eigenvalues of the symmetric tridiagonal matrix T from the eigenvalue storage section 18, and determines whether any of the eigenvalues is a negative value. If there is a negative value, the process proceeds to step S602. If not, the process proceeds to step S603.

[0117] (Step S602) The positive definite part 21 reads the symmetric tridiagonal matrix T from the diagonal matrix storage unit 13 and changes the eigenvector of the symmetric tridiagonal matrix T without changing the eigenvector. The symmetric tridiagonal matrix T after the symmetric tridiagonal matrix T is converted to positive definite so that all eigenvalues of the diagonal matrix T are positive, and the symmetric tridiagonal matrix T is converted to a diagonal matrix. The eigenvalues accumulated in the storage unit 13 and converted into positive definite values are accumulated in the eigenvalue storage unit 18.

[0118] (Step S603) The Cholesky decomposition unit 22 reads from the diagonal matrix storage unit 13 a symmetric tridiagonal matrix T in which all eigenvalues are positive values, and from the eigenvalue storage unit 18 the symmetric 3 Read the eigenvalues of the positive values of the superdiagonal matrix T. In addition, when positive definite conversion is performed, the symmetric tridiagonal matrix T and the eigenvalue after the positive definite value 匕 are read out. Then, the Cholesky decomposition unit ₂₂ performs a qd-type transformation on each element of the symmetric tridiagonal matrix Τ to convert the symmetric tridiagonal matrix T into the upper bidiagonal matrix B ( ⁺¹⁾ and the lower 2 Collessky decomposition into a multi-diagonal matrix B ( ^{_ 1)} .

[0119] (Step S604) The vector calculation unit 23 calculates the symmetry of each element of the upper double diagonal matrix B ( ⁺¹⁾ and the lower double diagonal matrix B ( ^_1) calculated by the Cholesky decomposition unit 22. The eigenvector is calculated using the eigenvalue of the positive value of the multidiagonal matrix T.

(Step S 605) Vector calculation unit 23 normalizes the calculated eigenvector. That is, the vector calculation unit 23 calculates the norm of the calculated eigenvector and accumulates the eigenvector divided by the calculated norm in the eigenvector storage unit 20 as the final eigenvector. In this way, the process of eigenvalue decomposition of the symmetric tridiagonal matrix T is completed. After this, processing for calculating the eigenvector of the symmetric matrix A stored in the matrix storage unit 11 may be performed, but is omitted here.

[0121] Further, in order to increase the accuracy of the eigenvector calculated by the eigenvector calculation unit 19, as shown in Fig. 12, the process related to the eigenvector may or may not be executed. May be. That is, the inverse iteration method processing unit (not shown) reads the eigenvector calculated by the eigenvector calculation unit 19 from the eigenvector storage unit 20, executes the inverse iteration process for the eigenvector, and stores the resulting eigenvector in the eigenvector storage. Store in part 20 (step S701). In this way, the accuracy of the eigenvector can be improved by executing the inverse iteration process. Next, the Gram Schmitt processing unit (not shown) determines whether high accuracy is necessary (step S702). This determination may be made by reading the setting from a recording medium or the like that has been previously set to determine whether high accuracy is required. Then, when high accuracy is required, the Gram Schmitt processing unit (not shown) reads the eigenvector that has been subjected to the inverse iteration process from the eigenvector storage unit 20, and executes the Gram Schmitt process on the eigenvector. Then, the resulting eigenbeta is stored in the eigenvector storage unit 20 (step S703). In this way, the orthogonality of eigenvectors can be improved by executing the Gramschmitt process. Note that the inverse iteration method and the Gramschmitt method are already known and will not be described in detail.

[0122] [Make symmetric tridiagonal matrix T positive definite]

Suppose that the matrix T can be decomposed into eigenvalues as follows.

[Equation 38]

T = QDQ ^T

[0123] Suppose that matrix T is multiplied by a and shifted by si as follows. Where I is a unit matrix of the same size as the matrix T, and a and s are real numbers.

[Equation 39]

T → T = aT + sI

[0124] Then, the eigenvalue decomposition after the shift is as follows.

[Equation 40] = QDQ ^T

However, D = aD + si.

[0125] Thus, when the matrix T is multiplied by a and shifted by si, each eigenvalue is multiplied by a and shifted by s. The eigenvector does not change. Therefore, the real number multiple and the shift Thus, after all eigenvalues have been set to positive values, eigenvectors can be calculated using the twist decomposition method.

[0126] When the calculated eigenvalue includes a negative value, for example, (A) the shift may be performed with the minimum eigenvalue, or (B) the shift may be performed with the maximum eigenvalue. . Here, it is assumed that the eigenvalue includes a negative value and a positive value. If the eigenvalues are all positive, it is not necessary to shift them. If all the eigenvalues are negative, if the matrix T is simply T, it can be positive definite.

[0127] (A) Shift by the smallest eigenvalue

If the smallest eigenvalue is 0, a = l and s = — λ + λ X ε. Where ε and

mm min max

As such, it may be any other appropriate value, which may be a value obtained by multiplying the machine epsilon by a constant. Thus, by shifting the matrix T, all eigenvalues of the matrix T can be made positive.

[0128] (B) Shift by maximum eigenvalue

If the maximum eigenvalue is> 0, then a = —1 and s = -λ Χ ε. In this case

max max mm

Will transform the matrix τ as follows.

[Equation 41]

Τ → Τ = -Τ + λ _ίΆ3Χ -λ _αήη χ εΙ =-( _Τ- (k _maK -λ ^ χ ε) /)

[0129] In the parentheses on the right side of the above equation, a positive definite value can be obtained by multiplying the whole force that makes all eigenvalues negative by shifting. Here, ε is a constant in the machine epsilon as in the case of (Α) (for example, it is a real number such as 2, 3, 4 and it is preferable to select a value for reducing the calculation error) You can use the value multiplied by or any other appropriate value.

[0130] If the absolute value of the minimum value of the eigenvalue is smaller than the absolute value of the maximum value of the eigenvalue, the absolute value force of the minimum value of the eigenvalue is shifted to the absolute value of the minimum value of the eigenvalue. When the value is larger than the value, the above shift (i) may be performed. By doing this, it is possible to reduce the amount of shift and reduce rounding errors.

The details of the process of calculating the eigenvector after the positive definite value 匕 will be explained together with the case of the following LV type twist decomposition method. [0131] (2) When calculating eigenvectors using the LV twist decomposition method

FIG. 13 is a diagram showing a configuration of the eigenvector calculation unit 19 when eigenvectors are calculated by the LV type twist decomposition method. In FIG. 13, the eigenvector calculation unit 19 includes a Cholesky decomposition unit 31 and a vector calculation unit 32.

The Cholesky decomposition unit 31 reads the symmetric tridiagonal matrix T from the diagonal matrix storage unit 13 and reads the eigenvalues of the symmetric tridiagonal matrix T from the eigenvalue storage unit 18. The Cholesky decomposition unit 31 performs the Miura transformation, dLVv type transformation, and inverse Miura transformation on each element of the symmetric tridiagonal matrix T, thereby transforming the symmetric tridiagonal matrix T into the upper bidiagonal matrix B. Cholesky decomposition into ⁽⁺¹⁾ and lower double diagonal matrix B ( ^_1) . Details of this processing will be described later. Note that “the Cholesky factorization of a symmetric tridiagonal matrix T into an upper bidiagonal matrix B ( ⁺¹⁾ and a lower bidiagonal matrix B ( ^_1) ” means that a symmetric tridiagonal matrix T A matrix obtained by multiplying the eigenvalue of x by a unit vector and subtracting it from the tridiagonal matrix T is converted into an upper bidiagonal matrix B ( ⁺¹⁾ and a lower bidiagonal matrix B ^(_1 , Cholesky decomposition It may be included.

[0133] The vector calculation unit 32 includes each element of the upper double diagonal matrix B ⁽⁺¹⁾ and the lower double diagonal matrix B ( ^_1 ) calculated by the Cholesky decomposition unit 31, and a symmetric tridiagonal matrix. The eigenvector is calculated using the eigenvalue of T and stored in the eigenvector storage unit 20. A matrix having this eigenvector in each column is an orthogonal matrix for eigenvalue decomposition.

Next, the process of step S103 in the flowchart of FIG. 2, that is, the process of calculating eigenvectors of the symmetric tridiagonal matrix T will be described in more detail. First, the processing of step S 103 in the flowchart of FIG. 2 will be described with reference to the flowchart of FIG.

(Step S801) The Cholesky decomposition unit 31 reads the symmetric tridiagonal matrix T from the diagonal matrix storage unit 13, and reads the eigenvalues of the symmetric tridiagonal matrix T from the eigenvalue storage unit 18. Then, the Cholesky decomposition unit ₃₁ performs the Miura transformation, dL Vv type transformation, and inverse Miura transformation on each element of the symmetric tridiagonal matrix T, so that the symmetric tridiagonal matrix T is Cholesky decomposition into angle matrix B and lower bidiagonal matrix B ^(_1) .

[0136] (Step S802) vector calculation unit 32, the eigenvalues for each element of the upper double diagonal matrix B ⁽⁺¹⁾ and the lower double diagonal matrix B ^(_ ^1>, symmetric tridiagonal matrix T And the eigenvector is calculated. [0137] (Step S803) The vector calculation unit 32 normalizes the calculated eigenvector. That is, the vector calculation unit 32 calculates the norm of the calculated eigenvector and accumulates the eigenvector divided by the calculated norm in the eigenvector storage unit 20 as the final eigenvector. In this way, the process of eigenvalue decomposition of the symmetric tridiagonal matrix T is completed.

[0138] Note that, after this, processing for calculating the eigenvector of the symmetric matrix A stored in the matrix storage unit 11 may be performed, but is omitted here. In addition, as shown in FIG. 12, in order to increase the accuracy of the eigenvector calculated by the eigenvector calculation unit 19, it may or may not be necessary to execute the process related to the eigenvector. The same as in the case of the type twist decomposition method.

[0139] [Calculation of eigenvectors from eigenvalues]

Regarding the process of calculating eigenvectors from eigenvalues, the method using the qd-type twist decomposition method and the method using the LV-type twist decomposition method are explained. As shown in the flow chart of Fig. 11, if eigenvectors are calculated using the qd-type twist decomposition method, if there are negative eigenvalues, the process of positively definite symmetric tridiagonal matrix T It is executed (steps S601 and S602). When the qd-type twist decomposition method is used, the symmetric tridiagonal matrix T and the eigenvalues after the positive definite value 匕 are used in the subsequent processing.

[0140] The diagonal matrix storage unit 13 stores the following symmetric tridiagonal matrix T of m X m, which is assumed to be:

[Number 42]

[0141] Also, it is assumed that the eigenvalue storage unit 18 stores an eigenvalue (l≤i≤m).

In such a case, finding the eigenvector of the symmetric tridiagonal matrix T results in the eigenvector X of the following equation.

[Equation 43] (Γ— λ) _Xj = γ _Λ = l, 2 ".., _m )

Where matrix Γ is an mxm matrix.

Is the vector of the other element force ^) with the element force.

(Is the 列 th column of the identity matrix /)

[0142] Originally, the right side of the above equation should be 0. When calculating the eigenvalue of the force matrix T, it contains some error, so if the eigenvector X is a true value, There is a residual term on the right side.

[0143] Here, it is assumed that the Cholesky decomposition can be performed as follows.

[Number 44]

T. Λ / =

(«= 0, + l).

l ("m) -2

^u 2 ("m) -l

T ·· ≡Β B

0

2m—1

, (

u2k-l-ir

¾) y = i, ^ "= 1, (" = 0, + l, — 1),

ξ = ξΓη. ), _& ^ _{) =} ^ _V

[0144] Here, the correspondence between each element of the matrix T and each element of the matrix B ^W) is as follows.

[Equation 45]

k = l, z, ',',

l2k 1k- \ ^u lk = 1,2,… "One 1

However, = 0.

[0145] The correspondence between each element of the matrix T and each element of the matrix B ^(G) can also be written as follows.

As is clear from the following formula force, each element of the matrix B ^½) can be calculated sequentially using each element of the matrix T. [Equation 46]

= i ..

) = = 1,2, ..., — 1 where. (°) = 0.

The results of the Cholesky decomposition of the above, can be written as: _c

[Equation 47]

I LD ⁺ L ^T

1 {UD ~ U ^L

Where, and are as follows.

And it becomes as follows.

[Number 48]

T-% _} I = N _k D _k {N

D _k ≡ di DD;, ..., D _k ⁺ —I, D _{k +} …, D _m ) Here, the matrix N is called a twisted matrix. Also

k

[Number 49] (T— λ ΐ) χ = γ e is

j j k k

[Number 50]

N

It becomes. The eigenvector can be calculated by solving this simple equation. In particular,

[Equation 51]

Where (p) is the pth element of.

Thus, the eigenvector can be calculated with a few operations. Note that for some p 0, D + = 0 or D — D, and the component of the matrix is b p 0 ρθ 2p

^(G) and b ^(G) , which is the, p + 1) component of matrix B, or the component of matrix T

-1 2p

make use of,

[Equation 52] + 2), ρθ く k,

(

p0) =

(°) (°)

ρ0> ^

^ 2 0-3 ^ 2 (°) 0-

^ 2 0 + 2

(p0 + 2), pO < _: ^ 2 0-4

(P0— 2), ρθ>

po- and eigenvectors can be calculated (this eigenvector calculation process is called exception processing). The values of the residual term parameters γ and k are

k

[Equation 53]

"" + — ¹ )-) _λ (Formula 5)

Is determined so that the absolute value of becomes the minimum. Therefore, the above-mentioned Cholesky decomposition is obtained and If it is possible to obtain the error matrix N, the eigenvector can be obtained.

k

Therefore, Cholesky decomposition will be described next.

[0149] As shown in Figure 15, T-lambda to Cholesky decomposition to I corresponds to ^{^{Β (0) {q (0}} ), e (

j k k

⁰⁾ } from {q ^{(± 1)} , e ^{(± 1} )} corresponding to B ^{(± 1)} .

k k

[0150] (qd-type twist decomposition method)

First, the qd-type transformation used in the qd-type twist decomposition method will be described (see Fig. 15).

[Equation 54] qd conversion: Γ, ^)}! ^^^ ¹ ), ^ ¹ )}

This transformation is also a stationary qd with shift (stqds) transformation

[Equation 55]

stqds conversion: {),. )} 4 {^ ⁺¹ ), ⁺¹ )}

^ + 卜 k = l, 2, m,

« ⁺¹⁾ =«. ), k = l, 2,., ml,

0, ≡ 0

And re verse—time progressive qd with shift (rpqds) conversion

[Equation 56] rpqds conversion:

^ + ^^^ + ^ Π- ,, =

M—Divided into ^ r ^, k = l, 2, ..., ml. If the eigenvalue is known, iterative calculation is not necessary, so the amount of calculation can be reduced. However, this method is not always highly numerically stable and accurate. This is because the stq ds conversion and the rpqds conversion may cause a loss of digits due to subtraction. For example, in the stqds transformation, if q ⁽⁰⁾ + e ⁽⁰⁾ — e ~ e (+ ", then find q and k k-1 j k-1 k

In this case, the number of significant digits in double-precision calculations may be too small. In that case, an error occurs if q ( ⁰⁾ e ⁽⁰⁾ / q is calculated. In other words, e cannot be calculated accurately kkkk

It becomes. Also, e is needed to find q, and q to find e

If k + 1 kkk is required, it is necessary to perform sequential calculation as described above. The error caused by this has spread, and there is a possibility of further error increase. As a result, q ⁽⁺¹⁾ ≠ 0 in theory, but q ⁽⁺¹⁾ = 0 due to error accumulation, and q ⁽⁰⁾ e (Vkk

In addition, numerical instability such as overflow may occur. Given the component {b ( ⁰⁾ , b ⁽⁰⁾ } of B ^(G) , that is, given {q ⁽⁰⁾ , e ⁽⁰⁾ }, λ and e

This situation is unavoidable because 2k- 1 2k k k j k- 1 is determined arbitrarily. Since the rpqds transform has similar properties, it cannot be said that it has reached a practical level. Although the improved version of the LAP ACK as the FO RTRAN routine DSTEGR has been released, the drawbacks have not been completely resolved.

[0151] In this way, each element of the diagonal matrix T is converted into each element of the matrix B ^(G) , and the matrix B

Even if qd-type conversion is performed for each element of q, since qd-type conversion is indirectly performed for each element of diagonal matrix T, “qd-type conversion is performed for each element of diagonal matrix T I will say.

Also, here, the force qd type conversion described for the standard qd type conversion may be a differen tial qd type conversion.

[0152] (LV type twist decomposition method)

Next, the Miura transform, dLVv transform, and inverse Miura transform used in the LV twist decomposition method are described (see Fig. 15).

[Equation 57] Miura transformation:

std LV v conversion: °) ^ u _k ⁽⁺¹⁾

r d L V v conversion:

Reverse Miura transform:

First, the Miura transform will be described. This conversion is shown as follows.

[Equation 58] k = \, 2, ..., m.

\, 2, ..., m- \,

≡ 0,

However, S ^w is arbitrary.

Next, dLVv type conversion will be described. This conversion is further converted to stationary discrete Lotka-Volterra with variable step-size (stdLVv)

[Equation 59] k = 1,…, 2 »ί— 1,

1 + δ ⁽⁺¹⁾ 0 'and re verse— time discrete Lotka― Volterra with variable step― size (rdLVv) conversion

[Equation 60]

It is divided into. However,

[Equation 61] = ^λ ("= ⁺¹ '— In the range that satisfies D, δ ^{(± 1)} is arbitrary.

[0155] Finally, the inverse Miura transform will be described. This conversion is shown as follows.

[Numerical 62]

(l ₊ ^ u ^ ₂ ), (l ₊ ^ u ^)

¾Jr _ g (+ i), ^κ , Α '

≡ 0 [0156] In this manner, Cholesky decomposition can be performed in the same manner as the qd conversion. A major feature of discrete Lotka-Volterra type transformation that is not seen in qd type transformation is that it has arbitrary parameters. That is, the value of δ ( ^η ) can be arbitrarily set within a range satisfying λ = 1Ζδ (°) −1Z δ ( ^{± 1} ). When δ ^(η) is changed, the value of the auxiliary variable u ^(η) also changes.

k

Whether differences or numerical instability will occur can be determined in advance. This determination may be implemented by an if statement. In this case, δ ^(η) may be calculated again after resetting. U ^{(± 1}

k

), Q ^{(± 1)} and e ^{(± 1)} are calculated independently by inverse Miura transform, so the error is

k k

It does not propagate. Note that reverse Miura conversion may be called Miura conversion, and Miura conversion may be called reverse Miura conversion. StdLVv conversion may be called stLVv conversion. RdLVv conversion may be called rLVv conversion.

[0157] Even if each element of the diagonal matrix T is converted into each element of the matrix B ^(G) and LV type conversion is performed on each element of the matrix B in this way, Indirectly, LV type conversion is performed for each element of the diagonal matrix T, so “LV type conversion is performed for each element of the diagonal matrix T”.

[0158] Here, an example of more detailed processing of processing by the LV type twist decomposition method will be described. 16 to 21 are flowcharts showing an example of processing by the LV type twist decomposition method.

FIG. 16 is a diagram showing an example of the entire process of Cholesky decomposition.

(Step S901) The Cholesky decomposition unit 31 performs Miura conversion. Details of this processing will be described later.

(Step S902) The Cholesky decomposition portion 31 sets 1Ζδ ^{(± 1} ) as lZδ ⁽⁰⁾ -λ.

(Step S903) The Cholesky decomposition portion 31 executes the processing of Procedure 1 described later.

[0161] (Step S904) The Cholesky decomposition portion 31 performs the process of Procedure 2 described later.

(Step S905) The Cholesky decomposition unit 31 determines whether q e has already been calculated.

k k

Judge whether. If it has already been calculated, the series of processes of Cholesky decomposition ends, whereas if it has not been calculated, the process returns to step S901. FIG. 17 is a flowchart showing details of the process of step S903 in the flowchart of FIG.

(Step S1001) Cholesky decomposition unit 31 has already calculated q e

k k

Judge whether it is power. If it has already been calculated, the process ends. If it has been calculated, the process proceeds to step S1002.

(Step S1002) The Cholesky decomposition unit 31 performs processing such as stdLVv conversion. Details of this process will be described later.

FIG. 18 is a flowchart showing details of the process of step S904 in the flowchart of FIG.

(Step S1101) Cholesky decomposition unit 31 has already calculated q e

k k

Judge whether it is power. If it has already been calculated, the process ends. If it has been calculated, the process proceeds to step S1102.

(Step S1102) The Cholesky decomposition unit 31 performs processing such as rdLVv conversion. Details of this process will be described later.

FIG. 19 is a flowchart showing details of the process in step S901 of the flowchart of FIG.

(Step S1201) The Cholesky decomposition portion 31 determines the value of 1Z δ ^(G) . Since this value can be arbitrarily determined as described above, the value is determined to be a positive definite value. For example, the Cholesky decomposition unit 31 determines the value of 1Z δ ( ^G) as follows:

[Equation 63]

-^ = | _min | + e (°) χ | λ 匪

[0165] It should be noted that ε ^(ω is a machine's epsilon multiplied by a constant (for example, a real number such as 2, 3, 4 is preferable to select a value for reducing calculation error). After that, for example, when it is determined in step S 1203 or the like that a digit loss may occur, the Cholesky decomposition unit 31 The value of ε ^(G) that determines the value of 1Z δ ^W may be set, for example, by increasing the constant applied to the machine epsilon by 1. The method of determining the value of 1Z δ is as follows: Not limited to this Needless to say. Thus, in the LV-type twist decomposition method, the value of 1Z δ ^(ω can be set freely, so that the matrix T does not have to be positively definite as in the case of the qd-type twist decomposition method. Of course, as with the qd-type twist decomposition method, the LV-type twist decomposition method can be executed after processing the positive definite value 匕 of the matrix string. .

(Step S1202) The Cholesky decomposition portion 31 sets β α to α () —lZ δ ( ^ω .

(Step S1203) The Cholesky decomposition portion 31 determines whether the absolute value of β 1 is greater than ε. If it is larger, the process proceeds to step S1204. If not, the process returns to step S1201, and the process of determining the value of 1 / δ ( ^G) is executed again. This ε can also be determined arbitrarily. For example, a machine epsilon may be used as ε. Increasing the value of ε improves the accuracy, and decreasing the value of ε decreases the accuracy. Note that the processing performed in step S1203 is processing for determining the possibility of a digit loss. If the absolute value of β1 is not larger than ε, half-JI will be cut off if there is a possibility of a digit loss.

[0167] (Step S1204) The Cholesky decomposition portion 31 sets β 2 to q ⁽⁰⁾ / (1 + δ ⁽⁰⁾ u ⁽⁰⁾ ).

Ten

The Since u ^((>) = 0 as described above, | 82 is set to q ⁽⁽⁾⁾ .

0 1

(Step S 1205) The Cholesky decomposition portion 31 sets u δ (°) u (°) to j81.

Ten

Here, 13 1 is the force set to q ⁽⁰⁾ -1 / δ () in step S 1202 Since u (.) = 0 as described above, β 1 is q (. 1 + δ (.) U (.) — This is equal to 1Ζ δ (.)

0 1 0

When Miura transformation is applied to, u ⁽⁰⁾ = u ⁽⁰⁾ I is also the force. From u (°) = 0 to (1 +

1 2k- 1 k = l 0

δ (% (°)) = 1.

0

(Step S 1206) The Cholesky decomposition portion 31 sets the counter k to “1”.

(Step S1207) The Cholesky decomposition portion 31 sets a ^ e ( ^G) Z j8 1. Said k

Thus, β 1 becomes u (°), so when Miura transformation is performed on e (°) Ζ ι8 1, a 1 is equal to δ ⁽⁰⁾ u (°)

2k- 1k 2k!

(Step S1208) The Cholesky decomposition portion 31 sets α 2 to α 1 +1. As can be seen from the explanation of step S 12 07, 22 is equal to 1 + S () u ().

2k

[0171] (Step S1209) The Cholesky decomposition portion 31 determines whether the absolute value of α 2 is greater than ε. to decide. If it is larger, the process proceeds to step S1210. If not, the process returns to step S1201, and the process of determining the value of 1 / Δto ⁾ is executed again. Note that the process performed in step S1209 is a process for determining the possibility of a digit loss, similar to the process in step S1203. If the absolute value of α 2 is not greater than ε, it is determined that there is a possibility of a digit loss.

[0172] (Step S1210) The Cholesky decomposition section 31 calculates <IX | 82 and u (°) (1+ S (°) u

2k 2k

Set to ⁽⁰⁾ ). As mentioned above, a 1 is equal to δ ⁽⁰⁾ u () and β 2 is q ⁽⁰⁾ / (1+ δ () u

-1 2k k 2

(.)), When Miura transform is performed on αΙΧβ2, k-2 2k 2k- 1 is equal to u (0) (1+ δ (o) u (.)).

(Step S1211) The Cholesky decomposition portion 31 sets | 82 to q ( ^G) Z α 2. K + 1 above

Since << 2 is equal to 1+ S (°) u (°), | 82 becomes q ⁽⁰⁾ / (l + 6 ⁽⁰⁾ u (°)),

2k k + 1 2k

This means that the value of k in ι82 has been incremented by 1.

[0174] (Step S1212) The Cholesky decomposition portion 31 sets the βΐ to beta2-1 [delta] ^W. As mentioned above, ι82 is equal to q (°) Ζ (1+ δ ((°)), so | 81 becomes q (°

k + 1 2k k + 1 so (1+ δ (

⁽⁰⁾ ) -1 / δ ^{(0) When} Miura transform is executed, β 1 becomes u ( ⁰⁾ . Therefore,

2k 2k + l

This means that the value of k in ι81 has been incremented by 1.

(Step S1213) The Cholesky decomposition portion 31 determines whether or not the absolute value of 131 is greater than ε. If it is larger, the process proceeds to step S1214. If not, the process returns to step S1201, and the process of determining the value of 1 / δ ( ^G) is executed again. Note that the process performed in step S1213 is a process for determining the possibility of a digit loss, similar to the process in step S1203. If the absolute value of β 1 is not larger than ε, it is judged that there is a possibility that a digit loss will occur.

[0176] (Step S1214) The Cholesky decomposition portion 31 calculates <2Χ | 81 and u δ (°)

2k + l

Set to u ()). As mentioned above, << 2 is equal to 1+ δ (% () j81 is equal to u ()

2k 2k 2k + 1 1

(Step S 1215) The Cholesky decomposition portion 31 increments the counter k by 1.

(Step S 1216) The Cholesky decomposition portion 31 determines whether or not the counter k is equal to m. If it is equal to m, the series of processing ends. Return to step SI 207.

[0178] In this way, u ⁽⁰⁾ (l + 6 ⁽⁰⁾ u ( ⁰⁾ ) and u ⁽⁰⁾ (l + 6 ⁽⁰⁾ u ⁽⁰⁾ ) are calculated.

2k + l 2k 2k 2k- 1

It will be. These values may be temporarily stored in a memory or the like (not shown) included in the Cholesky decomposition unit 31.

FIG. 20 is a flowchart showing details of the process of step S1002 of the flowchart of FIG.

(Step S1301) The Cholesky decomposition portion 31 sets ひ 2 to ¹ + δ ⁽⁺¹⁾ u. In addition

0

As described above, since u ⁽⁺¹⁾ = 0, α 2 is set to 1.

0

[0180] (Step S1302) The Cholesky decomposition portion 31 converts β 1 to u δ (°) u δ ⁽⁺¹

Ten

Set to u. Here, the value of u ⁽⁰⁾ (1+ δ ⁽⁰⁾ u ⁽⁰⁾ ) is calculated in step S1205.

0 1 0

Use the ones that were taken out. Note that u ⁽⁺¹⁾ = 0 as described above. In addition, stdL

0

When the Vv conversion is performed, β 1 is set to u.

(Step S1303) The Cholesky decomposition portion 31 sets the counter k to “1”.

(Step S1304) The Cholesky decomposition portion 31 sets α 1 to | 81 + 1Zδ.

[0182] (Step S1305) The Cholesky decomposition portion 31 determines whether the absolute value of α 1 is greater than ε. If it is larger, the process proceeds to step S1306; otherwise, the process proceeds to Procedure 2 (step S904) in the flowchart of FIG. Note that the process performed in step S 1305 is a process for determining the possibility of a digit loss, similar to the process in step S 1203. If the absolute value of a 1 is not greater than ε, it is determined that there is a possibility of a digit loss.

[0183] (Step S1306) The Cholesky decomposition portion 31 sets j82 to u δ () u ()) Ζα1

2k 2k- 1

Set. Here, the value of u ⁽⁰⁾ (1+ δ ⁽⁰⁾ u ⁽⁰⁾ ) was calculated in step S1210.

2k 2k- 1

Use things. Also, when stdLVv conversion is performed on this 132 formula, β 2 becomes δ ⁽⁺¹⁾ u

It is set to 2k.

(Step S 1307) The Cholesky decomposition portion 31 calculates α 1 X α2. As mentioned above, a

1 equals ι81 + 1 / δ ⁽⁺¹⁾ = ιι ( ⁺¹⁾ + 1 / δ ⁽⁺¹⁾ひ 2 equals l + S ⁽⁺¹⁾ u

2k- 1 2k- 2 Therefore, when the inverse Miura transform is performed on α IX, it becomes equal to q. Therefore, this k

The ski disassembling unit 31 calculates << IX 2 and sets it to q. (Step S 1308) The Cholesky decomposition portion 31 sets 0; 2 to 1 + β 2. As mentioned above, since j8 2 is equal to S ⁽⁺¹⁾ u, 2 is 1 + δ ⁽⁺¹⁾ u ⁽⁺¹⁾ . Therefore, to 0; 2

2k 2k

It means that the value of k has been incremented by 1.

(Step S 1309) The Cholesky decomposition portion 31 determines whether the force is greater than the absolute value force S _{ε of} α 2. If it is larger, the process proceeds to step S1310. Otherwise, the process proceeds to Procedure 2 (step S904) in the flowchart of FIG. Note that the process performed in step S1309 is a process for determining the possibility of a digit loss, similar to the process in step S1203. If the absolute value of α2 is not greater than ε, it is determined that there is a possibility of a digit loss.

(Step S1310) The Cholesky decomposition portion 31 calculates j8 1 X j8 2. As mentioned above, β 1 is equal to u and 2 is equal to S ⁽⁺¹⁾ u, so the inverse Miura transform is applied to | 8 1 Χ | 8 2

2k- 1 2k

When executed, it is equal to e. Therefore, the Cholesky decomposition unit 31 calculates β Ι Χ β 2.

k

And set to e.

k

[0188] (Step S1311) The Cholesky decomposition portion 31 converts | 8 1 into u ⁽⁰⁾ (1 + δ ⁽⁰⁾ u ⁽⁰⁾ ) / α 2

2k + l 2k

Set. Here, the value of u ⁽⁰⁾ (1 + δ ⁽⁰⁾ u is calculated in step S1214.

2k + l 2k

Use things. Also, when stdLVv transformation is performed on this β 1 expression, | 8 1 becomes u

2k + l is set. Therefore, the value of k at 13 1 is incremented by 1.

(Step S 1312) The Cholesky decomposition portion 31 increments the counter k by 1.

(Step S1313) The Cholesky decomposition portion 31 determines whether the counter k is equal to m. If it is equal to m, the process proceeds to step S1314; otherwise, the process returns to step S1304.

[0190] (Step S1314) The Cholesky decomposition portion 31 calculates α 2 X (1 + ΐΖ δ. Until this time, after updating a 1 in step S 1304, α Ι Χ α 2 in step S1307 Therefore, the Cholesky decomposition unit 31 calculates << 2 Χ (l + ΐΖ δ and sets it to q. In this way, the result of the Miura transform is converted to the stdLVv transform and the inverse. The Muller transformation is executed, and the process of calculating q ⁽⁺¹⁾ and e ends.

k k

Although not shown in the key disassembly unit 31, it may be temporarily stored in a memory or the like. FIG. 21 is a flowchart showing details of the processing of step SI 102 in the flowchart of FIG.

(Step S1401) The Cholesky decomposition unit 31 converts 131 to u ⁽⁰⁾ (1 + δ (.) U (.))

2m- 1 2m- 2 Z (1

+ Set to δ ^(_1) u. Here, the value of u ⁽⁰⁾ (l + 6 ⁽⁰⁾ u ( ⁰⁾ ) is

2m 2m- 1 2m- 2

Use the value calculated in S1214. Note that u ( ^_1) = 0 as described above. Also this

2m

If the rdLVv conversion is performed on the 131 expression, | 81 is set to u.

2m- 1

(Step S1402) The Cholesky decomposition portion 31 sets the counter k to “m−1”.

(Step S1403) The Cholesky decomposition portion 31 sets α 1 to | 81 + 1Z δ.

[0193] (Step S1404) The Cholesky decomposition portion 31 determines whether or not the absolute value of α 1 is larger than ε. If larger, the process proceeds to step S1405. Otherwise, the process returns to the Miura conversion (step S901) in the flowchart of FIG. Note that the processing performed in step S 1404 is processing for determining the possibility of a digit loss, similar to the processing in step S1203. If the absolute value of a 1 is not greater than ε, it is determined that there is a possibility of a digit loss.

[0194] (Step S1405) The Cholesky decomposition portion 31 sets j82 to u δ (°) u (°)) Ζα1

2k 2k- 1

Use things. Also, when rdLVv transformation is performed on this β 2 equation, β 2 becomes δ ^(_1) u

It is set to 2k.

(Step S1406) The Cholesky decomposition portion 31 sets 0; 2 to 1 + β2. As described above, since j82 is equal to S ( ^_1 ) u, one 2 is ¹ + δ ( ^_1 ) u.

2k 2k

(Step S 1407) The Cholesky decomposition portion 31 determines whether the force is larger than the absolute value force S _{ε of} α 2. If it is larger, the process proceeds to step S1408; otherwise, the process returns to the Miura conversion (step S901) in the flowchart of FIG. Note that the process performed in step S 1407 is a process for determining the possibility of a digit loss, similar to the process in step S1203. << If the absolute value of 2 is not larger than ε, it is judged that there is a possibility of a digit loss.

(Step S 1408) The Cholesky decomposition portion 31 calculates α 1 X α2. As mentioned above, a

1 equals ι81 + 1 / δ ^(_1) = ιι ( ^_1) + 1 / δ ^(_1) 2 equals 1+ S ^(_1) u Therefore, when the inverse Miura transform is performed on 1 1 X 2 2, it becomes equal to q. Therefore, k + 1

Reskey decomposition unit 31 calculates α 1 Χα 2 and sets it to q.

k + 1

[0198] (Step S1409) The Cholesky decomposition portion 31 converts β 1 to u ⁽⁰⁾ (1 + δ ⁽⁰⁾ u ⁽⁰⁾ ) / α 2

2k-1 Set to 2k-2. Here, the value of u ⁽⁰⁾ (1 + δ ⁽⁰⁾ u (°)) is the value of steps S1205, S12

2k- 1 2k- 2

Use the value calculated in 14. Also, when stdLVv conversion is performed on this j8 1 expression, 1 is set to u. Therefore, the value of k in 13 1 is decremented by 1.

2k- 1

That's what happened.

(Step S 1410) The Cholesky decomposition portion 31 calculates j8 1 X j8 2. As mentioned above, β

Since 1 is equal to u ¹ ) and 13 2 is equal to δ ^(_1) u, the inverse Miura transform is applied to IX 2.

2k- 1 2k

When executed, it is equal to e. Therefore, the Cholesky decomposition unit 31 calculates β Ι Χ β 2 k

And set to e.

k

[0200] (Step S1411) The Cholesky decomposition portion 31 decrements the counter k by 1.

(Step S1412) The Cholesky decomposition portion 31 determines whether or not the counter k is equal to zero. If the value is equal to 0, the process proceeds to step S 1413; otherwise, the process returns to step S 1403.

[0201] (Step S1413) The Cholesky decomposition portion 31 calculates 1 + δ δ. This is equivalent to calculating α Ι Χ α 2 in step S1408 after updating a 1 in step S1403. In this case, α 2 is 1 force. Therefore, the Cholesky decomposition unit 31 calculates | 8 1 + 1 / δ and sets it to q. In this way, the rdL Vv conversion and the inverse Miura conversion are executed on the result of the Miura conversion, and the process of calculating q ^(_1) and e is completed. This kk

These values may be temporarily stored in a memory or the like (not shown) included in the Cholesky decomposition unit 31.

[0202] Since the eigenvector corresponding to one eigenvalue can be calculated by the processing of the flowchart of Fig. 16, the Cholesky decomposition unit is used when eigenvectors corresponding to all eigenvalues are calculated. In step 31, the process of the flowchart in FIG. 16 is repeated by the number of eigenvalues.

[0203] An array for the auxiliary variable u ⁽ⁿ⁾ is not necessarily prepared to reduce memory consumption k

Absent. On the other hand, memory area for 1 + δ ^(G) u ^(G) is secured, Miura conversion, dLVv type conversion, reverse By using this value across the Miura conversion steps, memory consumption can be reduced and the amount of computation can be reduced. Thereby, an error is also reduced.

[0204] Here, the error will be described. Assuming that humans can perform any number of calculations in an infinite digit, there is no problem with either the qd-type twist decomposition method or the LV-type twist decomposition method. However, care must be taken when computing with a computer. On computers that cannot calculate and use limited digits, the correct results are not always obtained even if mathematically correct calculation methods are used. Not only that, there may be unexpected numerical problems that the calculation will not end indefinitely.

[0205] Errors due to computer calculations are known to include rounding errors and errors caused by dropping digits. Rounding error alone does not cause a large error to the extent that the last significant digit differs from the true value at most. In addition, rounding errors will still occur if at least one of the addition, multiplication, and division of two real numbers with different exponents is performed, but no further error will occur. In addition, even if operations that generate such rounding errors are repeated, if the rounding mode is near (rounded off), errors will not accumulate excessively due to unilateral rounding up or down. Few. Therefore, many numerical computation methods pay special attention to the rounding error caused by at least one operation of addition, multiplication, and division. Even with proper eigenvalue calculation, the rounding error does not increase uniformly.

[0206] The problem is the loss of digits caused by subtraction of real numbers with the same sign and addition of real numbers with different signs. If the value becomes 0 due to an error due to a digit loss and then division by that value is performed, it becomes indefinite so that 0 comes to the denominator and cannot be calculated. When this happens, the calculation will not end indefinitely. Since the subtraction-> division calculation part exists in both the qd-type twist decomposition method and the LV-type twist decomposition method, it is necessary to pay close attention to the precision error.

[0207] In the LV type twist decomposition method, whether or not the error due to the above-mentioned digit loss is included can be determined by whether or not the value obtained by subtraction is small. In the case of the qd-type twist decomposition method, even if it includes a fractional error, it cannot be avoided. Because given {q ⁽⁰⁾ , e ⁽ ° as an initial value, λ is uniquely determined and (qe

This is because kkjkk is also derived uniquely. In other words, it is a free V, calculation method that does not have arbitrary parameters. [0208] On the other hand, the LV-type twist decomposition method has a parameter δ ⁽⁰⁾ that can be set freely, so the value of the auxiliary variable u ^(η) can be changed variously (see Fig. 22). That is, like k

{Q ^{(± 1)} , e ^{(± 1)} } can be calculated by various routes. Therefore, when a digit loss occurs kk

Can also be avoided. The condition judgment of Fig. 19 to Fig. 21 is used to check the influence of the digit loss, and if the condition that the absolute value power of the value obtained by subtraction is greater than the number ε is not satisfied, the parameter δ ^(G) Return to setting. This process is repeated until the condition is met. If high speed is more important than accuracy, exception processing may be performed if the condition is not satisfied several times (if q ( ⁺¹ ) = 0 or q ( ^_1 ) = 0).

k k

[0209] It should be noted that the qd-type transformation can be rewritten as follows by using an expression indicating the correspondence between each element of the matrix T and each element of the matrix B ^(G) .

[Equation 64]

q d-type conversion:-^

[0210] This transformation is also a stqds transformation

[Equation 65]

stqds conversion:

q ^(l + e ^ = t _u _ ₁ - _J , k = U

qi ⁺¹⁾ ei ⁺¹⁾

= l, 2 ".., _m _l,

E0

And rpqds conversion

[Equation 66]

rpqds conversion: {t _2t i,}

+ e ^x ~ ^ j-> k = l, 2, ..., m,

It is divided into. Thus, when {q ^{(± 1)} , e ^{(± 1)} } is obtained from {t, t}, {t

2k- 1 2k k k 2k- 1

, t} to {q (.), e (.)} and {q (soil e ^{(± 1)} } from q (.), e (.)}

2k k k k k k k

In comparison, the amount of calculation is reduced, and the calculation speed can be improved.

[0211] In this case, qd type conversion is directly performed on each element of the diagonal matrix T. Therefore, if this transformation is `` do qd transformation on each element of diagonal matrix T '' Needless to say, it is included in the event.

[0212] Also, like the qd type transformation, using the formulas that show the correspondence between each element of the matrix T and the matrix Β ⁽ each element of the ^ω , the Miura transformation in the LV type transformation is rewritten as Can also

[Equation 67] Miura transformation: {υ ^ ί "^, ^)}

std LV v conversion:) u [ ⁺¹⁾

r d L V v conversion:. )

Reverse Miura transformation:

[0213] The Miura transform in this case is shown as follows.

[Equation 68] t = ^ (1 + 5 (.) — ₂ XI + δ () + δ ⁽ .) ₃ ) — ₂ k = l, 2,., M

(t _2k ) ²

k = 1,2,… -1

≡ ο,

However, S ^w is arbitrary.

Thus, to find {t, t in ⁽⁰⁾ , u (.)}, {T, 1;

2k- 1 2k 2k- 1 2k 2k- 1 2k ko), e (.)} And then {u ⁽⁰⁾ , u (.)} From q), e (.)} , Total kkk 2k- 1 2k

The calculation amount is reduced, and the calculation speed can be improved.

[0214] In this case, LV conversion is directly performed on each element of the diagonal matrix T. Therefore, it goes without saying that this conversion is included in the case of “performing LV type conversion for each element of the diagonal matrix T”.

[0215] If the Cholesky decomposition can be performed by the qd-type twist decomposition method or the LV-type twist decomposition method, the twisted matrix N can be calculated as described above, and k k of the matrix N

By solving N ^T x = e, χ can be calculated. here,

k j k j

[Equation 69]

INI

By substituting and X, this X is normalized. In this way, the eigenvector X is obtained. I can power. Use this unique vector X! / T, j T 1 2 m by setting Q = (x, X, ..., X)

The orthogonal matrix Q can be calculated.

T

[0216] Therefore, the eigenvalue decomposition is performed for the symmetric tridiagonal matrix T as shown in the following equation.

[Equation 70]

T = Q _T D _T Q

[0217] If eigenvalue decomposition of the symmetric tridiagonal matrix T can be performed, eigenvalue decomposition of the symmetric matrix A is also performed as follows.

[Equation 71]

Α = ΡΤΡ ^Ί

= P (Q _T D _T Q _T ^T ) P ^T

= (PQ _T ) D _T (PQ _T

[0218] As described above, the eigenvector calculation unit 19 performs the eigenvalues of the symmetric tridiagonal matrix T stored in the eigenvalue storage unit 18 and the symmetric triples stored in the diagonal matrix storage unit 13. The eigenvector of the symmetric tridiagonal matrix T is calculated using the diagonal matrix T. First, the Cholesky decomposition unit 31 reads the symmetric tridiagonal matrix T from the diagonal matrix storage unit 13 and reads the eigenvalues of the symmetric tridiagonal matrix T from the eigenvalue storage unit 18. Then, the Cholesky decomposition unit 31 performs each conversion as shown in the flowchart of FIG. 23 to perform the Cholesky decomposition process. Here, the case of using the LV type twist decomposition method will be described.

[0219] The Cholesky decomposition unit 31 first determines q ⁽⁰⁾ and e ⁽⁰⁾ kk from the values of the elements of the symmetric tridiagonal matrix T. Then, the Cholesky decomposition unit 31 performs the Miura transformation to obtain u ⁽⁰⁾ , u ⁽

1 2 ω, etc. are obtained sequentially (step S1501).

[0220] Next, the Cholesky decomposition unit 31 uses the u ^(ω obtained by the Miura transform to convert the stdLVv transform to k

As a result, u and the like are sequentially obtained (step S1502). In addition, the Cholesky decomposition unit 31 uses the u ^((>) obtained by the Miura transform to execute the rdLVv transform to

, U 卜^υ, etc. are sequentially obtained (step S1503).

2m- 1

[0221] Finally, the Cholesky decomposition unit 31 performs the u and rdLVv conversion obtained by the stdLVv conversion. Q ^{(± 1)} , e ^{(± 1),} etc. are sequentially converted into k 1 1 by performing inverse Miura transform using the obtained u.

(Step SI 504). When q ^{(± 1)} and e ^{(± 1)} are obtained, that is, upper double diagonal kk

When the column B ⁽⁺¹⁾ and the lower double diagonal matrix B ^(_1) are obtained, the Cholesky decomposition unit 31 passes them to the beta calculation unit 32. Note that the Cholesky decomposition unit 31 performs the Kolesky decomposition for each eigenvalue.

[0222] Here, for convenience of explanation, the Cholesky decomposition process has been described using the flowchart of FIG. 23. However, the process may be performed as shown in the flowcharts of FIGS. Needless to say.

[0223] The vector calculation unit 32 reads the eigenvalue from the eigenvalue storage unit 18 and determines the value of k using Equation 5. The vector calculation unit 32 uses the k value to calculate the twist kk from q ^{(± 1)} and e ^{(± 1).}

A matrix N is obtained and X is calculated by solving N ^T x = e (step S802).

k k j k j

[0224] Next, the vector calculation unit 32 obtains the eigenvector X by normalizing the calculated X and stores it in the eigenvector storage unit 20 (step S803). The vector calculation unit 32 can calculate all eigenvectors by performing the process of calculating X and the process of normalizing each eigenvalue. In this way, the calculation of the eigenvector ends. Thereafter, the matrix storage unit 11 may store the eigenvector of the symmetric matrix A! / ヽ. The calculation method is as described above.

In addition, the eigenvalue decomposition apparatus 1 may include an output unit (not shown) that outputs the eigenvalue calculated by the eigenvalue calculation unit 17 and the eigenvector calculated by the eigenvector calculation unit 19. Here, the output from the output unit (not shown) may be, for example, display on a display device such as an eigenvalue (for example, a CRT or a liquid crystal display), or transmission via a communication line to a predetermined device such as an eigenvalue. It is also possible to store eigenvalues on a recording medium, such as printing by eigenvalues. The output unit may or may not include an output device (for example, a display device or a printer). The output unit may be realized by hardware, or may be realized by software such as a driver for driving these devices.

[0226] Here, the case where the Cholesky decomposition unit 31 performs the Cholesky decomposition by the LV type twist decomposition method has been described, but the Cholesky decomposition unit 31 uses the qd type twist decomposition method. Cholesky decomposition may be performed. For example, the Cholesky decomposition unit 31 may look at the distribution of the calculated eigenvalues and use the qd-type twist decomposition method when the distribution is not dense.

[0227] Also, in the above description, the case where the matrix B ^to) is an upper bidiagonal matrix has been described, but the eigenvalue decomposition is similarly performed even if the matrix B ^(C>) is a lower bidiagonal matrix. Could be executed. However, the expression indicating the correspondence between each element of matrix T and each element of matrix B ^(G) is different from the above.

[0228] In the above description, the case where the eigenvalue calculation unit 17 calculates all eigenvalues has been described, but only some eigenvalues may be calculated. For example, in FIG. 9, the eigenvalue calculation unit 17 needs to calculate all eigenvalues up to the matrices Τ and T.

1 2

When calculating eigenvalues of symmetric tridiagonal matrix T from columns 列 and T, eigenvalues are calculated within the required range.

1 2

It may be calculated. Therefore, the eigenvalue calculation unit 17 may calculate at least one eigenvalue of the symmetric tridiagonal matrix T. In this case, the processing load can be reduced because it is not necessary to perform extra processing associated with calculating unnecessary eigenvalues. Thus, when calculating only some eigenvalues, for example, the calculation time spent calculating the eigenvalues of the symmetric tridiagonal matrix T can be approximately (1 + kZm) Z2 times. sell. Here, m is the matrix size, and k is the number of eigenvalues to be obtained.

[0229] In the above description, the case where the eigenvector calculation unit 19 calculates all eigenvectors corresponding to the calculated eigenvalue has been described. However, as is clear from the above description, the process of calculating the eigenvector is performed for each eigenvalue. Can be done. Therefore, the eigenvector calculation unit 19 may calculate at least one eigenvector of the symmetric tridiagonal matrix T. As described above, the eigenvector calculation unit 19 can calculate the eigenvector within a necessary range, and does not need to perform an extra process associated with calculating an unnecessary eigenvector, thereby reducing the processing load. Can be done.

Next, parallel processing in the eigenvalue decomposition apparatus 1 according to the present embodiment will be described. In the eigenvalue decomposition apparatus 1 according to the present embodiment, processing is executed in parallel in calculation of eigenvalues and eigenvectors. Can do. For example, as shown in FIGS. As described above, in the eigenvalue decomposition apparatus 1, the eigenvalue decomposition unit 15, the eigenvalue calculation unit 17, the correspondence key decomposition unit 22, the vector calculation unit 23, the Cholesky decomposition unit 31, and the vector calculation unit 32 perform parallel processing, respectively. Also good.

[0231] The eigenvalue decomposition unit 15 includes a plurality of eigenvalue decomposition means 15a and 15b, and the plurality of eigenvalue decomposition means 15a and 15b calculate eigenvalues and matrix elements of the symmetric tridiagonal matrix after division. The processing to be performed may be executed in parallel. For example, in FIG. 9, the eigenvalue decomposition means 15a performs eigenvalue decomposition of the matrices Τ, Τ, Τ, Τ, and the eigenvalue decomposition means 15b

111 112 121 121 211 τ.

May perform eigenvalue decomposition of 212, τ

221, τ 222

[0232] The eigenvalue calculation unit 17 includes a plurality of eigenvalue calculation means 17a and 17b, and the plurality of eigenvalue calculation means 17a and 17b calculate the eigenvalues and matrix elements of the symmetric tridiagonal matrix of the division source. The processing to be performed may be executed in parallel. For example, in FIG. 9, the eigenvalue calculating means 17a executes processing for calculating eigenvalues of the rows Τ, Τ, Τ, T, etc., and the eigenvalue calculating means 17b

11 12 1

Processing for calculating eigenvalues and the like of the columns τ 1, τ 2, T may be executed.

21 22 2

[0233] The Cholesky decomposition unit 22 includes a plurality of Cholesky decomposition units 22a and 22b, and the plurality of Cholesky decomposition units 22a and 22b execute a Cholesky decomposition process on a symmetric tridiagonal matrix T in parallel. May be. For example, the Cholesky decomposition means 22a may execute the process of performing Cholesky decomposition on half of the eigenvalues, and the Cholesky decomposition means 22b may execute the process of performing Cholesky decomposition on the other half of the eigenvalues. ,.

[0234] The vector calculation unit 23 may include a plurality of vector calculation units 23a and 23b, and the plurality of vector calculation units 23a and 23b may execute processes for calculating eigenvectors in parallel. For example, the vector calculating means 23a executes processing for calculating eigenvectors for the eigenvalues subjected to Cholesky decomposition by the Cholesky decomposition means 22a, and the vector calculating means 23b calculates eigenvectors for the eigenvalues subjected to Cholesky decomposition by the Cholesky decomposition means 22b. You can execute the calculation.

[0235] The Cholesky decomposition unit 31 includes a plurality of Cholesky decomposition units 31a and 31b, and may perform parallel Cholesky decomposition processing on the symmetric tridiagonal matrix T. Same as 22.

[0236] The vector calculation unit 32 includes a plurality of vector calculation means 32a and 32b, and the eigenvector is calculated. Similar to the vector calculation unit 23, the calculation processes may be executed in parallel.

[0237] The multiple eigenvalue calculating means 17a, 17b of the eigenvalue calculating unit 17 uses the divided eigenvalues and matrix elements of the two matrices Τ, T shown in FIG. Eigenvalues and

1 2 0

The process of calculating the matrix elements may be executed in parallel. The process will be described below.

[0238] First, it can be seen from Equation 1 that the processing for obtaining each eigenvalue can be executed in parallel.

For example, the eigenvalue calculating means 17a may calculate half of the eigenvalues, and the eigenvalue calculating means 17b may calculate the other half of the eigenvalues. Note that each of the eigenvalue calculating means 17a and 17b has values of components of a matrix D and a vector z.

12

[0239] It can also be seen from Equation 2 that the processing for obtaining each eigenvector q can be executed in parallel. For example, the eigenvalue calculating unit 17a may calculate a half eigenvector q corresponding to the calculated eigenvalue, and the eigenvalue calculating unit 17b may calculate the other half eigenvector q corresponding to the calculated eigenvalue.

[0240] From Equations 3 and 4, it can be seen that parallel processing is possible even with the processing for obtaining f and 1.

For example, the eigenvalue calculating unit 17a calculates some elements of f and 1 using the calculated eigenvector q, and the eigenvalue calculating unit 17b uses the calculated eigenvector α to calculate the remaining elements of f and 1. May be calculated. After that, the eigenvalue calculating means 17a or the eigenvalue calculating means 17b calculates f and 1 represented by the expressions 3 and 4 by integrating the elements f and 1 calculated by the eigenvalue calculating means 17a and 17b. be able to.

[0241] Here, the z component does not contain 0, and the diagonal component of D does not contain the same value.

12

The eigenvalue calculation means 17a, 17b explained in the above processing, if the component of the force z contains 0, even if the diagonal component of D contains the same value,

12

The process of calculating the eigenvector ^ and the process of calculating the matrix elements of the eigenvector force source matrix can be executed in parallel as well.

[0242] In parallel processing in each component as described above, when each means uses a memory for storing information such as eigenvalues, a plurality of means may use the same memory, that is, a shared memory. Alternatively, each means may use a different memory.

[0243] Note that the eigenvalue decomposition unit 15, the eigenvalue calculation unit 17 and the like have two hands using FIGS. Although the case where parallel processing is executed by stages has been described, the eigenvalue decomposition unit 15, the eigenvalue calculation unit 17, and the like may include three or more means and execute parallel processing. Also, the case where parallel processing is executed in each of the eigenvalue decomposition unit 15, the eigenvalue calculation unit 17, the Cholesky decomposition unit 22, the vector calculation unit 23, the Cholesky decomposition unit 31, and the vector calculation unit 32 has been described. In any one or more parts, parallel processing may not be executed. For example, the eigenvalue decomposition unit 15 may not perform parallel processing.

[0244] Further, the parallel processing may be executed by two or more devices, or may be executed by using two or more CPUs in one device. For example, as shown in FIG. 27, apparatus A and apparatus B may include eigenvalue decomposition means 15a and 15b, respectively, and eigenvalue decomposition processing may be executed in parallel in each apparatus. In this case, the eigenvalue decomposition unit 15-1 including the eigenvalue decomposition unit 15a of the device A and the eigenvalue decomposition unit 15-2 including the eigenvalue decomposition unit 15b of the device B constitute the eigenvalue decomposition unit 15. become. Therefore, the eigenvalue decomposition apparatus 1 constitutes a system including the apparatus A and the apparatus B. Here, the parallel processing of the eigenvalue decomposition unit 15 has been described, but the parallel processing by two or more devices may be performed for the other eigenvalue calculation units 17 and the Koresky decomposition unit 22 and the like.

[0245] Next, performance evaluation by numerical experiments will be described. Here, the method (DDC) of the eigenvalue decomposition apparatus 1 according to the present embodiment, the standard divide-and-conquer method (DC), the QR method (QR), and the method combining the bisection method and the inverse iteration method (B & I) was compared. In the method of the eigenvalue decomposition apparatus 1 according to the present embodiment, the Cholesky decomposition is performed by the qd-type twist decomposition method, and the inverse iteration (see step S701 in FIG. 12) is used once. DSTEDC provided by LAPACK was used as the divide-and-conquer method. As the QR method, DSTEQR provided by LAPACK was used. As a method of combining the bisection method and the inverse iteration method, a part of the DSTEVR source code provided by LAPACK was used. LAPACK used in the experiment calls BLAS (Basic Linear Algebra Subprograms). Further, the eigenvalue decomposition unit 15 of the eigenvalue decomposition apparatus 1 according to the present embodiment performs eigenvalue decomposition using the QR method. In this numerical experiment, a diagonal matrix given at random was symmetric and tridiagonalized using the Lunchos method. In the graph of the numerical experiment results, the horizontal axis represents the size of the matrix used in the numerical experiment. [0246] FIG. 28 is a diagram showing calculation execution time. The shorter the execution time, the faster the calculation. As can be seen from the graph of the numerical experiment results in FIG. 28, the method of the eigenvalue decomposition apparatus 1 according to the present embodiment can perform eigenvalue decomposition at a higher speed than other methods. .

[0247] Fig. 29 is a diagram showing the orthogonality of the calculated eigenvectors. Smaller value of force is high. As can be seen from the graph of the numerical experiment results in FIG. 29, in the method of the eigenvalue decomposition apparatus 1 according to the present embodiment, the orthogonality of the eigenvectors is not high compared to the other methods. It can be seen that the method of the eigenvalue decomposition apparatus 1 according to the present embodiment has the same degree of orthogonality as the method combining the bisection method and the inverse iteration method. This is probably because the eigenvalues are calculated first and then the eigenvectors are calculated.

FIG. 30 is a diagram showing the sum of relative errors of all eigenvalues. The smaller the value, the more accurate the error is. As can be seen from the graph of the numerical experiment results in FIG. 30, it can be seen that the method of the eigenvalue decomposition apparatus 1 according to the present embodiment has higher eigenvalue accuracy than the other methods. It has the same eigenvalue accuracy as the eigenvalue calculated by the divide-and-conquer method. In the method of the eigenvalue decomposition apparatus 1 according to the present embodiment, the eigenvalue is calculated by the divide-and-conquer method.

FIG. 31 is a diagram showing the sum of errors of the previous eigenvectors. Small value of the force error is small and accuracy is high. As can be seen from the graph of the numerical experiment results in FIG. 31, it can be seen that the method of the eigenvalue decomposition apparatus 1 according to this embodiment has higher eigenvector accuracy than the other methods.

[0250] Thus, it can be seen that the method of the eigenvalue decomposition apparatus 1 according to the present embodiment can perform eigenvalue decomposition at high speed, and has high eigenvalue and eigenvector accuracy. In addition, as described above, since parallelism is also excellent, it is possible to achieve higher speed by performing parallel processing.

[0251] [Face image recognition by principal component analysis]

Next, the case where eigenvalue decomposition is applied to face image recognition by principal component analysis will be described.

[0252] Face image recognition processing by principal component analysis is performed by the following steps. (1) Creating a projection matrix

(2) Registration of information on personal face images

(3) Recognition by comparing the input face image with information related to registered individual face images

[0253] First, the creation of the projection matrix (1) will be described.

Prepare face image data of various people. The face image data may be two-dimensional image data read by an optical reading device such as a scanner, a digital still camera, or a digital video camera. Each face image data can be displayed as a vector having each pixel value as a component. The vector is as follows. Here, M face image data are handled.

[Equation 72]

^Λ 1 ' ^X 2''M

Next, the average value of the vector of the face image data is calculated.

[Equation 73]

Μ

[0255] Using the calculated z, the covariance matrix C is calculated as follows.

[Equation 74]

M M

^C = ∑∑ (--) ^τ

[0256] After that, the eigenvalue decomposition apparatus 1 performs eigenvalue decomposition on the matrix C. That is, the matrix C is stored in the matrix storage unit 11, and eigenvalue decomposition is performed by performing processing such as diagonalization and matrix division on the matrix C. The eigenvalues resulting from the eigenvalue decomposition are accumulated in the eigenvalue storage unit 18. The eigenvector V, which is the result of eigenvalue decomposition, is stored in the eigenvector storage unit 20, and the orthogonal matrix used for diagonalization is applied to the eigenvector as described above. Can be requested. Note that the subscript j is set to be in descending order of the eigenvalues. In other words, the eigenvalue is the maximum eigenvalue. [0257] Among the obtained eigenvectors, the projection matrix r? Is calculated as follows using d eigenvectors with the corresponding eigenvalues from the top.

[Equation 75]

= (v ^ v,,---^,)

[0258] Next, registration of information related to the individual face image in (2) will be described.

M face image data for each individual

[Equation 76]

Prepare. Here, the subscript k is a subscript for identifying an individual. By multiplying each vector by the projection matrix 7 ?, m d-dimensional vectors ξ are calculated as shown below, and the average value is calculated. Then, the average value / is registered.

[Equation 77]

i = η i = l, 2, "', m

1

= 丄 ∑

[0259] Next, recognition by comparing the input face image in (3) above with information related to registered individual face images will be described.

Assume that a vector X of face image data of a person is input. Then, the d-dimensional vector is calculated by multiplying the vector X of the face image data by the projection matrix 7 ?.

[Equation 78]

[0260] Next, the norm of the difference between the d-dimensional vector ξ and the average value ^{k of} each individual's d-dimensional vector ξ is calculated.

[Equation 79]

I-<

[0261] If the norm is preliminarily set and is smaller than the value, it may be determined that the face image of the individual is identified by the input facial image power index k. Alternatively, the above norm is calculated for all subscripts k, and the input face image corresponds to the minimum norm. It may be determined that the face image of the individual is identified by the subscript k. In this way, it is possible to perform face image recognition by principal component analysis using the eigenvalue decomposition apparatus 1.

[0262] [Application to image processing to restore 2D image to 3D]

Next, the case where eigenvalue decomposition is applied to image processing to restore an object from a 2D image to a 3D image is explained.

[0263] The process of performing a 3D restoration from a plurality of 2D images is performed by the following steps.

(1) Two-dimensional image force A step of extracting feature points.

(2) A step of calculating data on the shape (three-dimensional coordinate data of the feature point of the original object) and rotation (conversion to the three-dimensional data force feature point data) from the feature point data.

(3) A step of visualizing from the shape and rotation data.

[0264] The steps (1) and (2) will be described below with reference to the flowchart of FIG. Here, the two-dimensional image may be a two-dimensional image read by an optical reading device such as a scanner, a digital still camera, or a digital video camera.

[0265] In step S5001, a two-dimensional image; j (j = l, ···, m, m is an integer greater than or equal to 3), and feature point i (i = l, ···, n, n is greater than 2) The coordinates (x ^j , y ^j ). The two-dimensional image to be handled is preferably a weak center projected image. At this time, the following equation holds.

[Equation 80]

, 'No

[0266] where s is the scale of the jth image relative to the scale of the object, and r and r are respectively

j l, j 2, j

This is the first and second row vectors of the rotation matrix of the jth camera coordinate system relative to the object coordinate system, and the three-dimensional coordinates of the (Xi, Yi, Zi) th point. The scale of the object is the same as the scale of the first image (s = 1), and the orientation of the object's coordinate system is the same as the camera coordinate system of the first image (r = (1, 0, 0) ^T , r = (0, 1, 0) ^T ). If the coordinates of all η points in all m images are arranged as elements of matrix D, the following equation is obtained.

[0267] [Equation 81] D = MS

{1

Where D = y {

,

[0268] As can be seen from the shapes of M and S, the rank of D is 3. Here, in step S5001

The matrix D is given. The data M and shape S related to rotation are obtained below.

[0269] So, singular value decomposition of matrix D

D = U∑V ^T

think of. Here, ∑ is a singular value arranged in diagonal order in the order of magnitude, and U and V are a left orthogonal matrix and a right orthogonal matrix, respectively. In this singular value decomposition, the above eigenvalue decomposition can be used. That is, in the eigenvalue decomposition apparatus 1, the eigenvalue decomposition of the symmetric matrix A is performed as described above by setting the symmetric matrix A stored in the matrix storage unit 11 to D ^T D. . In the eigenvalue decomposition apparatus 1, what is stored in the eigenvector storage unit 20 is the eigenvector of the symmetric tridiagonal matrix T obtained by transforming the symmetric matrix D ^T D. Therefore, the eigenvector is expressed as described above. Thus, it is necessary to convert to the eigenvector of the matrix D ^T D. Each singular value is the 1Z square (square root) of each eigenvalue of the matrix D ^T D. That is, the singular value σ = (λ) ^1/2 . In addition, the eigenbeta is equal to the right singular vector X that forms each column of the right orthogonal matrix V.

j

Matrix force with right column in each column. If the singular value is not σ force ^, the left singular vector y composing each column of the left orthogonal matrix U can be obtained as follows. [0270] [Equation 82]

It becomes. On the other hand, if the singular value σ force ^

[Number 83]

By solving for = 0, y can be calculated. Here, the left singular vector y is normalized by performing the following replacement.

[Number 84]

[0271] In this way, the left singular vector y can be obtained. By using this left singular vector y j] and U = (y, y, ..., y), the left orthogonal matrix U can be calculated.

[0272] Here, because of the digital error of the image, it is not zero! There are 3 or more singular values. However, the fourth and subsequent singular values are due to noise and are much smaller than the first three singular values. Therefore, in step S5002, singular vectors are calculated for the first three singular values. That is, the eigenvalue decomposition apparatus 1 according to the present embodiment calculates eigenvectors for the first three eigenvalues. Summing up the three vectors to be adopted, the following equation is obtained.

[0273] D '= L'∑'R' ^T = M'S '

Where M '= L' (∑ ') ^1/2 , S' = (∑ ') ^1/2 R' ^T , and D 'are rank 3 matrices that minimize II DD' II.

[0274] Next, we want to find M and S from D ', but the combination is not unique. Because any regular matrix is

D '= (M'C) (C ^_1 S,)

It is also the power to satisfy. Therefore, C in the above equation is determined to satisfy M = M'C. Also satisfies the following formula.

[Equation 85]

[0275] E = CC ^, we get 2m + 1 linear equations for the six elements of the above equation E. Since m≥3, the elements of E can be uniquely determined. In step S5003, a matrix E is obtained.

[0276] Next, in step S5004,! The same degree of freedom (9) is greater than the degree of freedom of E (6). Therefore, if condition = (1, 0, 0) ^T , r = (0, 1, 0) ^τ is added, C is determined lj ¾

Can do. At this time, two solutions (Necker Reversal) come out.

Next, in step S5005, M = from M'C and S = C ^_1 S ', data M and the shape S is determined relating to the rotation.

[0277] [Application to document search]

Next, a case where eigenvalue decomposition is applied to document search processing will be described. After the process of extracting the index word related to the content of the document from the document and calculating the weight of the index word, in the vector space model, the document is expressed by a vector having the index word weight as an element. Here, the documents to be searched are d 1, d 2,.

1 2 n

Suppose that there are m index words w 1, w 2,. At this time, the document d is

1 2 m j

It will be expressed in Torr. This is called a document vector.

[Equation 86]

[0278] Here, d is the weight of the index word w in the document d. The entire document set can be expressed by the m X n matrix D as follows.

[Equation 87]

[0279] The matrix D is called an index word document matrix. Similarly, each column of the index word document matrix is a document vector that represents information about the document. Similarly, each row of the index word document matrix is a vector that represents information about the index word. This is called the outer ring. Search queries can also be expressed in the same way as documents, with the index word weight as an element. If the weight of the index word w included in the search question sentence is q, the search question vector q is expressed as follows.

[Equation 88]

[0280] In an actual document search, it is necessary to find a document that is similar to a given search query sentence, but this is done by calculating the similarity between the search query vector q and each document rule d. . There are various definitions of the similarity between vectors, but what is often used in document search is the cosine measure (the angle between two vectors) or the inner product.

[0281] [Equation 89] d, q ∑ '

cos (d q) = (cosine measure

'q = ∑ (Inner product If the vector length is normal 正規 (cosine normalized), the cosine scale and inner product match.

FIG. 33 is a flowchart showing an example of a document search method using the eigenvalue decomposition apparatus 1 according to this embodiment. In step S6001,! /, The question vector q is received.

[0283] Here, a search using an approximate matrix D of D is considered. Search in vector space model

k

The search is performed by calculating the similarity between the query vector q and each document vector d in the index word document matrix D. Here, D is used instead of D. In the vector space model, the sentence

k

The number of dimensions of the book vector is equal to the total number of index words. Therefore, the number of document vectors tends to increase as the number of documents to be searched increases. However, as the number of dimensions increases, problems such as computer memory limitations and search time increase will arise, and unnecessary index words included in documents will have a noisy effect. There is also a phenomenon that decreases the search accuracy. Latent semantic indexing (LSI) is a technology that improves search accuracy by projecting document vectors in a high-dimensional space into a low-dimensional space. Index word power treated separately in high dimensional space Since it may be treated as interrelated in low dimensional space, it is possible to search based on the meaning and concept of index words. it can. For example, in an ordinary vector space model, the index word "car" and the index word "automobile" are completely different, and a query with one index word cannot retrieve documents that contain the other index word. However, in a low-dimensional space, these semantically related index terms can be expected to degenerate into one dimension, so that the search query “car” does not rely on the power of documents containing “car”. It is possible to search for documents that contain "." In latent semantic indexing, dimensionality reduction of high-dimensional vectors is performed by singular value decomposition, which is basically equivalent to principal component analysis in multivariate analysis. is there.

[0284] In step S6002, select! / And k. Here, k <r is selected as k <r. r = min (n, m). The value of k may be given in advance or may be selectable for each calculation.

[0285] Next, in step S6003, singular value decomposition of the matrix D is performed. As this singular value decomposition, a method based on eigenvalue decomposition can be used as described in the application to image processing for restoring from a two-dimensional image to a three-dimensional image. In other words, in the eigenvalue decomposition apparatus 1, the eigenvalue decomposition of the symmetric matrix A is performed as described above by setting the symmetric matrix A stored in the matrix storage unit 11 as the matrix D ^T D. And, as a result, the matrix D Singular value decomposition can be performed. In this singular value decomposition, the singular vector of D is calculated for k singular values of the calculated singular values in descending order, up to the kth singular value. That is, in the eigenvalue decomposition, eigenvectors may be calculated for k eigenvalues from the first to the k-th in descending order among the calculated eigenvalues, and a singular vector may be calculated using these eigenvectors. k is the value selected in step S6002.

[0286] That is,

D = U∑V ^T

k k k

Calculate u and V. Where U is composed of only the first k left singular vectors k k k

M X k matrix, where V is only the first k right singular vectors n X

k

k matrix, and ∑ is a k X k diagonal matrix where only the first k singular values are also forces.

[0287] Next, in step 6004, the similarity between the matrix D and the question vector q is calculated. V ヽ

k

If the vector e is an n-dimensional unit vector, the j-th document vector of D is represented by D e

j k k j The similarity calculation between the document vector D e and the search query exterior q is

[Number 90]

φ ejfo)

However, another definition may be used. In the above equation, D is reconstructed from U,, and V

k k k k

This shows that the similarity can be calculated directly from the result of singular value decomposition. ∑ V ^T e that appears in the above equation is

k k j

∑ V ^T e = U ^T D e

k k j k k j

Can be rewritten. The right side of this expression is the jth document vector in the approximation matrix D

k

Represents the coordinates (k-dimensional representation of the document) under the base U of. Similarly, U in the above equation

k k

^T q is the coordinate (k-dimensional representation of the search question) under the base U of the search question vector q. [0288] In step S6005, a search result is output based on the similarity calculated in step S6004.

Thus, the eigenvalue decomposition apparatus 1 according to the present embodiment can be used for various processes such as image processing and search processing. Needless to say, it may be used for processes other than those described above.

[0289] As described above, since the eigenvalue decomposition apparatus 1 according to the present embodiment calculates only the eigenvalues using the divide-and-conquer method, the vector update executed in the standard divide-and-conquer method is not necessary, and the standard It is much faster than a simple divide-and-conquer method. In addition, the process of calculating the eigenvalue force eigenbeta can be performed at high speed. Furthermore, parallelization at the stage of calculating eigenvalues is difficult in the QR method. However, in the eigenvalue decomposition apparatus 1 according to the present embodiment, each of the stage of calculating the eigenvalue and the stage of calculating the eigenvalue force eigenvector are performed. And essentially high parallelism. In addition, it can be seen that the eigenvalue decomposition apparatus 1 according to the present embodiment can obtain higher accuracy than other methods.

[0290] The eigenvalue decomposition apparatus 1 according to the present embodiment may be a stand-alone apparatus or a server apparatus in a server client system, as described above. It may be a system. When the eigenvalue decomposition device 1 is a server device in the server client system, it is stored in the symmetric matrix A stored in the matrix storage unit 11 or the eigenvalue or eigenvector storage unit 20 stored in the eigenvalue storage unit 18. Eigenvectors are sent and received via communication lines such as the Internet and Intranet.

[0291] Also, in the calculation of the eigenvector in the eigenvalue decomposition apparatus 1 according to the present embodiment, the calculation of the eigenvector using a certain eigenvalue uses the qd-type twist decomposition method, and the calculation of the eigenvector using another eigenvalue. The LV type twist decomposition method may be used. For example, for eigenvalues close to each other, the eigenvector is calculated using the LV type twist decomposition method, and for eigenvalues where the values are not close, the eigenvector is calculated using the qd type twist decomposition method. Let's go. Whether or not the values of the eigenvalues are close to each other can be set in a predetermined recording medium or the like, and can be cut in half by comparing with the values. [0292] Also, in the above embodiment, each processing or each function may be realized by centralized processing by a single device or a single system, or may be realized by a plurality of devices or It may be realized by being distributedly processed by a plurality of systems.

[0293] In the above-described embodiment, each component may be configured by dedicated hardware, or a component that can be realized by software may be realized by executing a program. For example, each component can be realized by a program execution unit such as a CPU reading and executing a software program recorded on a recording medium such as a hard disk or a semiconductor memory. The software that realizes the eigenvalue decomposition apparatus in the above embodiment is the following program. In other words, this program reads out the symmetric tridiagonal matrix T from the diagonal matrix storage unit in which the symmetric tridiagonal matrix T is stored in the computer, and stores the symmetric tridiagonal matrix T in two pieces. Dividing into a symmetric tridiagonal matrix and storing it in the diagonal matrix storage unit, dividing the symmetric tridiagonal matrix into two symmetric tridiagonal matrices and storing it in the diagonal matrix storage unit This process is repeated until each divided tridiagonal matrix after division is less than or equal to a predetermined size, and each symmetric tridiagonal matrix that is less than or equal to the predetermined size. Read from the diagonal matrix storage unit, perform eigenvalue decomposition on each symmetric tridiagonal matrix, and determine the eigenvalues of each symmetric tridiagonal matrix and the eigenvalues of each symmetric tridiagonal matrix. Luka computes at least the matrix element of the orthogonal matrix, and the eigenvalue and the matrix element Is stored in the eigenvalue decomposition storage unit, the eigenvalues and matrix elements of each symmetric tridiagonal matrix are read from the eigenvalue decomposition storage unit, and the eigenvalues and the matrix elements are read from the eigenvalue decomposition storage unit. The eigenvalues of the symmetric tridiagonal matrix and the matrix elements of the symmetric tridiagonal matrix of the division source are calculated and stored in the eigenvalue decomposition storage unit, the eigenvalues of the symmetric tridiagonal matrix of the division source, The process of calculating matrix elements is repeated until at least one eigenvalue of the symmetric tridiagonal matrix T is calculated, and at least one eigenvalue of the symmetric tridiagonal matrix T is stored in the eigenvalue storage unit. The symmetric tridiagonal matrix T is read from the diagonal matrix storage unit, the eigenvalues of the symmetric tridiagonal matrix T are read from the eigenvalue storage unit, and the symmetric tridiagonal matrix is read out. Twist resolution from T and its eigenvalue Method to calculate at least one eigenvector of the symmetric tridiagonal matrix T and And an eigenvector calculation step stored in the vector storage unit.

In this program, the eigenvector may be calculated by the qd-type twist resolution method in the eigenvector calculation step.

[0294] Also, in this program, the eigenvector calculation step force The eigenvalues of the symmetric tridiagonal matrix T are read from the eigenvalue storage unit, and when any of the eigenvalues is a negative value, the diagonal Read out the symmetric tridiagonal matrix T from the matrix storage unit so that all eigenvalues of the symmetric tridiagonal matrix T are positive values without changing the eigenvectors of the symmetric tridiagonal matrix T. Then, the symmetric tridiagonal matrix T is converted to a positive definite value, the eigenvalues after the positive definite value 匕 is stored in the eigenvalue storage unit, and the symmetric tridiagonal matrix T after being converted to a positive definite value is stored in the diagonal matrix. A positive definite value accumulating step and a symmetric tridiagonal matrix T in which all eigenvalues are positive values are read from the diagonal matrix storage unit, and the symmetric tridiagonal matrix is read from the eigenvalue storage unit. Read the eigenvalue of the positive value of T. By performing qd-type transformation on each element of the diagonal matrix T, the symmetric tridiagonal matrix T is converted into an upper bidiagonal matrix B ⁽⁺¹⁾ and a lower bidiagonal matrix B ^(_1). a Cholesky decomposition step for decomposing, and each element of the upper double diagonal matrix B ⁽⁺¹⁾ and the lower double diagonal matrix B ^{(_ 1),} a positive value before Symbol symmetric tridiagonal matrix T And a vector calculation step of calculating an eigenvector using the eigenvalue and storing the eigenvector in the eigenvector storage unit.

In this program, the eigenvector may be calculated by the LV type twist decomposition method in the eigenvector calculation step.

[0295] Also, in this program, the eigenvector calculation step force reads the symmetric tridiagonal matrix T from the diagonal matrix storage unit, reads the eigenvalues of the symmetric tridiagonal matrix T from the eigenvalue storage unit, By performing Miura transformation, dLVv type transformation, and inverse Miura transformation for each element of the symmetric tridiagonal matrix T, the symmetric tridiagonal matrix T is transformed into an upper bidiagonal matrix B ⁽⁺¹⁾ and lower 2 a Cholesky decomposing Cholesky decomposition step heavy diagonal matrix B ^(_1), and each element of the upper double diagonal matrix B ⁽⁺¹⁾ and the lower double diagonal matrix B ^(_1), the symmetry 3 A vector calculation step of calculating an eigenvector using the eigenvalues of the superdiagonal matrix T and storing the eigenvector in the eigenvector storage unit. [0296] In this program, the computer reads the symmetric matrix A from the matrix storage unit in which the symmetric matrix A is stored, and the symmetric tridiagonal matrix T is obtained by tridiagonalizing the symmetric matrix A. Further, a diagonal step of calculating and storing in the diagonal matrix storage unit may be executed.

In the above program, the step of accumulating information is performed only by hardware and does not include at least processing! /.

[0297] In addition, this program is recorded on a predetermined recording medium (for example, an optical disk such as a CD-ROM, a magnetic disk, a semiconductor memory, etc.) that can be executed by being downloaded by force such as a server. The program may be executed by being read.

[0298] Further, the computer that executes this program may be singular or plural. That is, centralized processing or distributed processing may be performed.

FIG. 34 is a schematic diagram showing an example of the appearance of a computer that executes the program and realizes the eigenvalue decomposition apparatus 1 according to the embodiment. The above embodiment can be realized by computer hardware and a computer program executed on the computer hardware.

In FIG. 34, the computer system 100 includes a computer 101 including a CD-ROM (Compact Disk Read Only Memory) drive 105, an FD (Flexible Disk) drive 106, a keyboard 102, a mouse 103, and a monitor. 104.

[0301] FIG. 35 is a diagram showing a computer system. In FIG. 35, a computer 101 includes a CD-ROM drive 105 and an FD drive 106, a CPU (Central Processing Unit) 111, and a ROM (Read Only Memory) 112 for storing a program such as a bootup program. The CPU 111 is connected to the CPU 111 to temporarily store the instructions of the application program, and to provide a temporary storage space. The RAM (Random Access Memory) 113, the hard disk 114 to store the application program, the system program, and the data, and the CPU 111 And a bus 115 for mutually connecting the ROM 112 and the like. The computer 101 may include a network card (not shown) that provides connection to the LAN.

[0302] The function of the eigenvalue decomposition apparatus 1 according to the above embodiment is implemented in the computer system 100. The program to be executed may be stored in the CD-ROM 121 or FD 122, inserted into the CD-ROM drive 105 or FD drive 106, and transferred to the hard disk 114. Alternatively, the program may be transmitted to the computer 101 via a network (not shown) and stored in the hard disk 114. The program is loaded into RAM 113 when executed. Note that the program may be loaded directly from the CD-ROM 121, the FD 122, or the network.

[0303] Further, the matrix storage unit 11, the diagonal matrix storage unit 13, the eigenvalue decomposition storage unit 16, the eigenvalue storage unit 18, and the eigenvector storage unit 20 may be realized by the RAMI 13 or the hard disk 114.

[0304] The program does not necessarily include an operating system (OS) or a third-party program that causes the computer 101 to execute the functions of the eigenvalue decomposition apparatus 1 according to the above-described embodiment. The program may contain only the part of the instruction that calls the appropriate function (module) in a controlled manner and achieves the desired result. How the computer system 100 operates is well known and will not be described in detail.

Further, the present invention can be variously modified without being limited to the above embodiment, and it goes without saying that these are also included in the scope of the present invention.

Also, while only a few exemplary embodiments of the present invention have been described in detail above, many changes may be made in the exemplary embodiments without substantially departing from the benefits of the invention and the new technology. One skilled in the art will readily recognize that this is possible. Accordingly, all such modifications are intended to be included within the scope of this invention.

Industrial applicability

[0305] As described above, according to the eigenvalue decomposition apparatus and the like according to the present invention, eigenvalue decomposition can be processed at high speed, and chemical calculation processing, statistical calculation processing, information retrieval processing, and other eigenvalues by the molecular orbital method can be performed. It is useful for an apparatus that performs processing using decomposition.

Brief Description of Drawings

FIG. 1 is a block diagram showing a configuration of an eigenvalue decomposition apparatus according to Embodiment 1 of the present invention.

FIG. 2 is a flowchart showing the operation of the eigenvalue decomposition apparatus according to the embodiment. [3] A flowchart showing the operation of the eigenvalue decomposition apparatus according to the embodiment.

[4] A flowchart showing the operation of the eigenvalue decomposition apparatus according to the embodiment.

圆 5] Diagram for explaining division of symmetric tridiagonal matrix T in the same embodiment 圆 6] Diagram for explaining calculation of eigenvalues of symmetric tridiagonal matrix T in the same embodiment

[7] A flowchart showing the operation of the eigenvalue decomposition apparatus according to the embodiment.

[8] A flowchart showing the operation of the eigenvalue decomposition apparatus according to the embodiment.

[9] Diagram for explaining calculation of eigenvalues of symmetric tridiagonal matrix T in the same embodiment

FIG. 10 is a block diagram showing the configuration of the eigenvector calculation unit in the embodiment. 圆 11] A flowchart showing the operation of the eigenvalue decomposition apparatus according to the embodiment. 圆 12] A flowchart showing the operation of the eigenvalue decomposition apparatus according to the embodiment.

FIG. 13 is a block diagram showing the configuration of the eigenvector calculation unit in the embodiment. 圆 14] A flowchart showing the operation of the eigenvalue decomposition apparatus according to the embodiment.

FIG. 15 is a diagram for explaining Cholesky decomposition in the same embodiment. 圆 16] A flowchart showing the operation of the eigenvalue decomposition apparatus according to the embodiment. 圆 17] A flowchart showing the operation of the eigenvalue decomposition apparatus according to the embodiment.圆 18] Flowchart showing the operation of the eigenvalue decomposition apparatus according to the embodiment 圆 19] Flowchart showing the operation of the eigenvalue decomposition apparatus according to the embodiment] 20] Flowchart showing the operation of the eigenvalue decomposition apparatus according to the embodiment 圆 21 Flowchart showing the operation of the eigenvalue decomposition apparatus according to the embodiment

FIG. 22 is a diagram for explaining the qd-type twist decomposition method and the LV-type twist decomposition method in the same embodiment.

圆 23] Flow chart showing the operation of the eigenvalue decomposition apparatus according to the embodiment 圆 24] Block diagram showing another example of the configuration of the eigenvalue decomposition apparatus according to the embodiment

FIG. 25 is a block diagram showing another example of the configuration of the eigenvector calculation unit in the embodiment.

FIG. 26 is a block diagram showing another example of the configuration of the eigenvector calculation unit in the embodiment. Figure

圆 27] A block diagram showing another example of the configuration of the eigenvalue decomposition apparatus according to the embodiment. 圆 28] A diagram showing the result of a numerical experiment of the eigenvalue decomposition apparatus according to the embodiment. 圆 29] The eigenvalue decomposition according to the embodiment. Fig. 30 shows the results of the numerical experiment of the device. 圆 30] Fig. 31 shows the results of the numerical experiment of the eigenvalue decomposition device according to the embodiment. 圆 31] The result of the numerical experiment of the eigenvalue decomposition device according to the embodiment.

FIG. 32 is a flowchart showing an example of image processing in the embodiment.

FIG. 33 is a flowchart showing an example of a document search process in the embodiment. 圆 34] A schematic diagram showing an example of the external appearance of the computer system in the embodiment. FIG. 35 shows an example of the configuration of the computer system in the embodiment. Illustration

Claims

The scope of the claims

A diagonal matrix storage for storing a symmetric tridiagonal matrix T;

Reading the symmetric tridiagonal matrix から from the diagonal matrix storage unit, dividing the symmetric tridiagonal matrix Τ into two symmetric tridiagonal matrices, and storing them in the diagonal matrix storage unit, The symmetrical tridiagonal matrix is divided into two symmetric tridiagonal matrices and stored in the diagonal matrix storage unit. A matrix partitioning unit that repeats until

Eigenvalues that store the eigenvalues of each symmetric tridiagonal matrix that is less than or equal to the predetermined size and the matrix elements that are part of the orthogonal matrix that is the eigenvector force of each symmetric tridiagonal matrix A decomposition storage unit;

Each symmetric tridiagonal matrix having a predetermined size or less is read from the diagonal matrix storage unit, eigenvalue decomposition is performed on each symmetric tridiagonal matrix, and each symmetric tridiagonal matrix is obtained. And an eigenvalue decomposition unit that stores at least the eigenvalue and the matrix element in the eigenvalue decomposition storage unit.

An eigenvalue storage for storing eigenvalues of the symmetric tridiagonal matrix Τ,

The eigenvalues and matrix elements of each symmetric tridiagonal matrix are read from the eigenvalue decomposition storage unit, and the eigenvalues and the eigenvalues of the symmetric tridiagonal matrix of the division source are divided from the eigenvalues and the matrix elements. The matrix element of the original symmetric tridiagonal matrix is calculated and stored in the eigenvalue decomposition storage unit, and the process of calculating the eigenvalues and matrix elements of the symmetric tridiagonal matrix of the division source is performed as symmetric 3 An eigenvalue calculator that stores at least one eigenvalue of the symmetric tridiagonal matrix Τ in the eigenvalue storage unit, and repeats until at least one eigenvalue of the multidiagonal matrix Τ is calculated; and the symmetric tridiagonal An eigenvector storage unit storing eigenvectors of the matrix Τ, reading the symmetric tridiagonal matrix から from the diagonal matrix storage unit, and reading eigenvalues of the symmetric tridiagonal matrix から from the eigenvalue storage unit, Symmetric tridiagonal matrix Τ and its eigenvalues From eigenvalue decomposition apparatus and a eigenvector computing section for storing the eigenvector storage unit calculates at least one of the eigenvectors of Τ the symmetric tridiagonal matrix using twisted factorization.

[2] The eigenvalue decomposition apparatus according to claim 1, wherein the eigenvector calculation unit calculates an eigenvector by a qd-type twist decomposition method.

[3] The eigenvector calculation unit

Read eigenvalues of the symmetric tridiagonal matrix T from the eigenvalue storage unit, and read the symmetric tridiagonal matrix T from the diagonal matrix storage unit when any of the eigenvalues is a negative value. The symmetric tridiagonal matrix T is positively definite so that all eigenvalues of the symmetric tridiagonal matrix T without changing the eigenvectors of the symmetric tridiagonal matrix T are positive. The eigenvalue is stored in the eigenvalue storage unit after the value is definite and positive definite, and the definite tridiagonal matrix T is stored in the diagonal matrix storage unit after being positive definite. And a symmetric tridiagonal matrix T in which all eigenvalues are positive values from the diagonal matrix storage unit, and a positive eigenvalue of the symmetric tridiagonal matrix T from the eigenvalue storage unit. By performing qd-type transformation on each element of the symmetric tridiagonal matrix T. A Cholesky decomposition unit that performs Cholesky decomposition of a tridiagonal matrix T into an upper double diagonal matrix B ( ⁺¹⁾ and a lower double diagonal matrix B ( ^_1) ;

An eigenvector is calculated using each element of the upper bidiagonal matrix B ( ⁺¹⁾ and lower bidiagonal matrix B ( ^_1) and the eigenvalue of the positive value of the symmetric tridiagonal matrix T. The eigenvalue decomposition apparatus according to claim 2, further comprising a vector calculation unit that accumulates in the eigenvector storage unit.

[4] The eigenvalue decomposition apparatus according to claim 1, wherein the eigenvector calculation unit calculates an eigenvector by an LV-type twist decomposition method.

[5] The eigenvector calculation unit includes:

The symmetric tridiagonal matrix T is read out from the diagonal matrix storage unit, the eigenvalues of the symmetric tridiagonal matrix T are read out from the eigenvalue storage unit, and Miura for each element of the symmetric tridiagonal matrix T is read out. Cholesky decomposition of the symmetric tridiagonal matrix T into upper bidiagonal matrix B ( ⁺¹⁾ and lower bidiagonal matrix B ( ^_1) by performing transformation, dLVv type transformation, and inverse Miura transformation Cholesky decomposition part,

The eigenvector is calculated by using each element of the upper bidiagonal matrix B ( ⁺¹⁾ and lower bidiagonal matrix B ( ^_1) and the eigenvalues of the symmetric tridiagonal matrix. 5. The eigenvalue decomposition apparatus according to claim 4, further comprising a vector calculation unit that accumulates in the storage unit.

[6] The Cholesky decomposition unit includes a plurality of Cholesky decomposition means,

6. The eigenvalue decomposition apparatus according to claim 3, wherein the plurality of Cholesky decomposition units execute a process of performing Cholesky decomposition on the symmetric tridiagonal matrix T in parallel.

[7] The vector calculation unit includes a plurality of vector calculation means,

7. The eigenvalue decomposition apparatus according to claim 3, wherein the plurality of vector calculating means execute processing for calculating eigenvectors in parallel.

[8] The eigenvalue calculation unit includes a plurality of eigenvalue calculation means,

The eigenvalue decomposition apparatus according to any one of claims 1 to 7, wherein the plurality of eigenvalue calculation means execute in parallel processing of calculating eigenvalues of a symmetric tridiagonal matrix and matrix elements.

[9] The eigenvalue decomposition unit includes a plurality of eigenvalue decomposition means,

The eigenvalue decomposition apparatus according to any one of claims 1 to 8, wherein the plurality of eigenvalue decomposition means execute in parallel a process of performing eigenvalue decomposition on a symmetric tridiagonal matrix.

[10] a matrix storage unit for storing the symmetric matrix A;

A diagonalization unit that reads the symmetric matrix A from the matrix storage unit, calculates the symmetric tridiagonal matrix T obtained by tridiagonalizing the symmetric matrix A, and stores the symmetric matrix A in the diagonal matrix storage unit; The eigenvalue decomposition apparatus according to any one of claims 1 to 9, further comprising:

11. The eigenvalue decomposition apparatus according to any one of claims 1 to 10, wherein the matrix dividing unit divides a symmetric tridiagonal matrix into two substantially half symmetric tridiagonal matrices.

[12] The eigenvalue decomposition apparatus according to any one of claims 1 to 11, wherein the eigenvalue calculation unit calculates all eigenvalues of the symmetric tridiagonal matrix T.

13. The eigenvalue decomposition apparatus according to claim 12, wherein the eigenvector calculation unit calculates all eigenvectors of the symmetric tridiagonal matrix T.

14. The eigenvalue decomposition apparatus according to claim 10, wherein the symmetric matrix A is a covariance matrix calculated as an average force of vectors representing face image data.

[15] A diagonal matrix storage unit storing a symmetric tridiagonal matrix T, a matrix partitioning unit, an eigenvalue decomposition unit, and an eigenvalue decomposed by the eigenvalue decomposition unit and an orthogonality such as an eigenvalue Luke An eigenvalue decomposition storage unit that stores matrix elements that are part of the matrix, and the symmetric 3 An eigenvalue decomposition apparatus comprising an eigenvalue storage unit that stores eigenvalues of a multidiagonal matrix T, an eigenvalue calculation unit, an eigenvector storage unit that stores eigenvectors of the symmetric tridiagonal matrix Τ, and an eigenvector calculation unit The eigenvalue decomposition method used in

The matrix division unit reads the symmetric tridiagonal matrix から from the diagonal matrix storage unit, divides the symmetric tridiagonal matrix Τ into two symmetric tridiagonal matrices, and The process of accumulating in the storage unit, dividing the symmetric tridiagonal matrix into two symmetric tridiagonal matrices and storing in the diagonal matrix storage unit is performed by each symmetric tridiagonal matrix after division. A matrix partitioning step that repeats until it is less than or equal to a predetermined size;

The eigenvalue decomposition unit reads each symmetric tridiagonal matrix having a predetermined size or less from the diagonal matrix storage unit, performs eigenvalue decomposition on each symmetric tridiagonal matrix, and Calculate at least the eigenvalues of each symmetric tridiagonal matrix and the matrix elements of the orthogonal matrix that is the eigenvector force of each symmetric tridiagonal matrix, and store the eigenvalues and the matrix elements in the eigenvalue decomposition storage unit Eigenvalue decomposition step,

The eigenvalue calculation unit reads the eigenvalues and matrix elements of each symmetric tridiagonal matrix from the eigenvalue decomposition storage unit, and from the eigenvalues and the matrix elements, the eigenvalues of the symmetric tridiagonal matrix of the division source And the matrix element of the symmetric tridiagonal matrix of the division source are stored in the eigenvalue decomposition storage unit, and the eigenvalue and matrix element of the symmetric tridiagonal matrix of the division source are calculated. The process is repeated until at least one eigenvalue of the symmetric tridiagonal matrix Τ is calculated, and at least one eigenvalue of the symmetric tridiagonal matrix Τ is accumulated in the eigenvalue storage unit. When,

The eigenvector calculation unit reads the symmetric tridiagonal matrix から from the diagonal matrix storage unit, reads the eigenvalues of the symmetric tridiagonal matrix から from the eigenvalue storage unit, and the symmetric tridiagonal An eigenvalue decomposition step comprising: calculating at least one eigenvector of the symmetric tridiagonal matrix から from the matrix Τ and its eigenvalue using a twist decomposition method and storing the eigenvector in the eigenvector storage unit. Method.

On the computer,

Symmetric tridiagonal matrix Diagonal matrix storage power storing される Read the symmetric tridiagonal matrix Τ and divide the symmetric tridiagonal matrix Τ into two symmetric tridiagonal matrices The diagonal line The process of accumulating in the column storage unit and dividing the symmetric tridiagonal matrix into two symmetric tridiagonal matrices and storing them in the diagonal matrix storage unit is performed for each symmetric tridiagonal after division. A matrix partitioning step that repeats until the matrix is less than or equal to a predetermined size;

Each symmetric tridiagonal matrix having a predetermined size or less is read from the diagonal matrix storage unit, eigenvalue decomposition is performed on each symmetric tridiagonal matrix, and each symmetric tridiagonal matrix is obtained. An eigenvalue decomposition step of calculating at least an eigenvalue of the matrix and a matrix element of an orthogonal matrix composed of eigenvectors of each symmetric tridiagonal matrix and storing the eigenvalue and the matrix element in an eigenvalue decomposition storage unit;

The eigenvalues and matrix elements of each symmetric tridiagonal matrix are read from the eigenvalue decomposition storage unit, and the eigenvalues and the eigenvalues of the symmetric tridiagonal matrix of the division source are divided from the eigenvalues and the matrix elements. The matrix element of the original symmetric tridiagonal matrix is calculated and stored in the eigenvalue decomposition storage unit, and the process of calculating the eigenvalues and matrix elements of the symmetric tridiagonal matrix of the division source is performed as symmetric 3 An eigenvalue calculation step of repeatedly storing at least one eigenvalue of the symmetric tridiagonal matrix T in an eigenvalue storage unit until at least one eigenvalue of the multidiagonal matrix T is calculated; from the diagonal matrix storage unit; The symmetric tridiagonal matrix T is read out, the eigenvalues of the symmetric tridiagonal matrix T are read out from the eigenvalue storage unit, and the symmetric tridiagonal matrix T and its eigenvalues are used to extract the eigenvalues using the twist decomposition method. At least one fixed tridiagonal matrix T Program for executing the eigenvector computing step of storing the eigenvector storage part to calculate a vector, a.