[go: up one dir, main page]

WO2008018188A1 - dispositif de décomposition de valeur propre et procédé de décomposition de valeur propre - Google Patents

dispositif de décomposition de valeur propre et procédé de décomposition de valeur propre Download PDF

Info

Publication number
WO2008018188A1
WO2008018188A1 PCT/JP2007/051575 JP2007051575W WO2008018188A1 WO 2008018188 A1 WO2008018188 A1 WO 2008018188A1 JP 2007051575 W JP2007051575 W JP 2007051575W WO 2008018188 A1 WO2008018188 A1 WO 2008018188A1
Authority
WO
WIPO (PCT)
Prior art keywords
matrix
eigenvalue
symmetric tridiagonal
storage unit
symmetric
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Ceased
Application number
PCT/JP2007/051575
Other languages
English (en)
French (fr)
Inventor
Yoshimasa Nakamura
Hiroaki Tsuboi
Taro Konda
Masashi Iwasaki
Masami Takata
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Kyoto University NUC
Original Assignee
Kyoto University NUC
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Kyoto University NUC filed Critical Kyoto University NUC
Priority to US12/376,444 priority Critical patent/US8255447B2/en
Priority to JP2008528725A priority patent/JP5017666B2/ja
Publication of WO2008018188A1 publication Critical patent/WO2008018188A1/ja
Anticipated expiration legal-status Critical
Ceased legal-status Critical Current

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization

Definitions

  • the present invention relates to an eigenvalue decomposition apparatus that performs eigenvalue decomposition.
  • 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, Numewitz Mathematik, Vol. 3 6, p. 177-195, 1981
  • 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.
  • an eigenvalue decomposition apparatus 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.
  • 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.
  • 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.
  • the eigenvector calculation unit may calculate an eigenvector by a qd-type twist decomposition method.
  • 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.
  • 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.
  • the eigenvector calculation unit may calculate an eigenvector by an LV type twist decomposition method.
  • the eigenvector can be calculated numerically stably.
  • 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 (+1) 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;
  • 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.
  • the Cholesky decomposition process can be performed in a short time.
  • 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.
  • the eigenvector calculation process can be performed in a short time.
  • 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.
  • the process of calculating the eigenvalue can be performed in a short time.
  • 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.
  • the eigenvalue decomposition process can be performed in a short time.
  • a matrix storage unit that stores a symmetric matrix A
  • the symmetric matrix A is read from the matrix storage unit
  • the symmetric matrix A is tridiagonalized.
  • a diagonalization unit that calculates a tridiagonal matrix T and accumulates it in the diagonal matrix storage unit.
  • an eigenvalue can be calculated for an arbitrary symmetric matrix A.
  • the matrix dividing unit may divide the symmetric tridiagonal matrix into two symmetric tridiagonal matrices of approximately half.
  • the eigenvalue calculation unit may calculate all eigenvalues of the symmetric tridiagonal matrix T.
  • the eigenvector calculation unit may calculate all eigenvectors of the symmetric tridiagonal matrix T.
  • the symmetric matrix A may be a covariance matrix in which an average force of a vector indicating face image data is calculated.
  • the eigenvalue decomposition apparatus can be used for facial image data recognition processing.
  • 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.
  • FIG. 1 is a block diagram showing the configuration of the eigenvalue decomposition apparatus 1 according to this embodiment.
  • the eigenvalue decomposition apparatus 1 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.
  • 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.
  • 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.
  • 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.
  • 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.
  • 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.
  • 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.
  • 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.
  • 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.
  • 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.
  • the eigenvalue decomposition of the symmetric tridiagonal matrix T is completed.
  • 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.
  • 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.
  • step 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.
  • step S102 of the flowchart of FIG. 2 will be described using the flowchart of FIG.
  • 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.
  • 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.
  • 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.
  • the symmetric matrix A can be tridiagonalized, as shown below, using, for example, No, Usholder transform, etc.
  • the symmetric tridiagonal matrix T is a square matrix in which the number of rows matches the number of columns.
  • 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).
  • ⁇ 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 ⁇ 2 is ( ⁇ -? 3 ⁇ 4) of X Real symmetric tridiagonal matrix,
  • ⁇ is a non-zero constant
  • ⁇ is an integer such that ⁇ ⁇ ⁇ .
  • the specific process of dividing an 8 x 8 symmetric tridiagonal matrix ⁇ into half the size is as follows.
  • each element specified in each matrix may be non-zero or zero.
  • 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
  • the matrix partitioning unit 14 is divided into two symmetric tridiagonal matrices ⁇ , T and two values
  • 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.
  • 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.
  • 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.
  • 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.
  • the matrix divider 14 may perform other processes in step S304.
  • 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.
  • 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). .
  • 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,
  • the matrix partitioning unit 14 has a symmetric tridiagonal matrix T that is also symmetric with the symmetric tridiagonal matrix T.
  • each symmetric tridiagonal matrix is a square matrix.
  • the eigenvalue decomposition of the matrix T is as follows.
  • 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.
  • 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).
  • 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.
  • the eigenvalue decomposition unit 15 performs each symmetric tri-diagonal matrix ⁇ , T,
  • Eigenvalue decomposition is performed on ⁇ , ⁇ ,.
  • 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.
  • 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.
  • the matrices ⁇ and ⁇ are both square matrices.
  • the symmetric tridiagonal matrix T of the division source is as follows.
  • 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 ⁇
  • D diag (, ⁇ ,..., ⁇ ).
  • Q is an eigenvalue of the matrix (D + P ZZ).
  • 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.
  • the eigenvector is the small matrix excluding the j-th row and j-th column of the matrix (D + p zz 1 ).
  • 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.
  • 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,
  • the eigenvalues are d 1, d 2,..., D and values obtained by solving the following eigen equations.
  • 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
  • ⁇ 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.
  • the matrix Q can be obtained from the matrix Q and the matrix Q, and the matrix D is also obtained.
  • 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).
  • F ⁇ is the first line of 3 ⁇ 4, and ⁇ 2 is the last line of ⁇ 2 .
  • f is the first row of the orthogonal matrix resulting from eigenvalue decomposition of the matrix ⁇ .
  • the first row and the last row of the orthogonal matrix are matrix elements.
  • 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.
  • the matrix T of the division source is constructed from the matrices, and T.
  • Equation 3 and Equation 4 are
  • 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.
  • 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.
  • 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.
  • 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.
  • Shall be rearranged. Specifically, the i-th component and j-th component of the diagonal matrix D are input.
  • 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.
  • 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.
  • 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
  • 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
  • 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.
  • the matrix Q which is not the matrix Q, is directly applied to f etc. I will calculate 1 You may do it.
  • FIG. 9 is a diagram for explaining processing for calculating eigenvalues.
  • the eigenvalue decomposition unit 15 performs eigenvalue decomposition on the matrices ⁇ , ⁇ , ...
  • the eigenvalue calculation unit 17 starts the first calculation of eigenvalues and the like (steps S401 and S402).
  • 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.
  • the matrix ⁇ and the matrix ⁇ are read from the eigenvalue decomposition storage unit 16.
  • the eigenvalue calculation unit 17 uses these values to calculate the eigenvalues of the matrix ⁇ .
  • the eigenvalue calculation unit 17 calculates an orthogonal matrix Q using the eigenvalues of the matrix T (step S502), and the orthogonal matrix
  • 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.
  • the eigenvalues and matrix elements of the matrix T are calculated (steps S402 and S403). Like this
  • 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.
  • 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.
  • 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.
  • 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.
  • 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.
  • the eigenvector calculation unit 19 includes a positive definite value conversion unit 21, a Cholesky decomposition unit 22, and a vector calculation unit 23.
  • 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 ⁇ .
  • 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.
  • 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.
  • 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 ⁇ .
  • 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.
  • 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 ( +1) and the lower 2 Cholesky decomposition into a bidiagonal matrix B ( _1) . Details of this process will be described later.
  • the Cholesky decomposition of a symmetric tridiagonal matrix T into an upper bidiagonal matrix B ( +1) 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 ( +1) and the lower bidiagonal matrix B (_1) May be included.
  • the vector calculation unit 23 includes each element of the upper double diagonal matrix B (+1) 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.
  • 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.
  • 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.
  • 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.
  • 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.
  • the Cholesky decomposition unit 22 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 ( +1) and the lower 2 Collessky decomposition into a multi-diagonal matrix B ( _ 1) .
  • Step S604 The vector calculation unit 23 calculates the symmetry of each element of the upper double diagonal matrix B ( +1) 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.
  • 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.
  • the Gram Schmitt processing unit 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.
  • matrix T is multiplied by a and shifted by si as follows.
  • I is a unit matrix of the same size as the matrix T, and a and s are real numbers.
  • the shift may be performed with the minimum eigenvalue, or (B) the shift may be performed with the maximum eigenvalue.
  • 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.
  • ⁇ 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.
  • the absolute value force of the minimum value of the eigenvalue is shifted to the absolute value of the minimum value of the eigenvalue.
  • the above shift (i) may be performed. By doing this, it is possible to reduce the amount of shift and reduce rounding errors.
  • 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.
  • 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 (+1) and lower double diagonal matrix B ( _1) . Details of this processing will be described later.
  • the Cholesky factorization of a symmetric tridiagonal matrix T into an upper bidiagonal matrix B ( +1) 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 ( +1) and a lower bidiagonal matrix B (_1 , Cholesky decomposition It may be included.
  • the vector calculation unit 32 includes each element of the upper double diagonal matrix B (+1) 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.
  • 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.
  • 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 31 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) .
  • Step S802 vector calculation unit 32, the eigenvalues for each element of the upper double diagonal matrix B (+1) and the lower double diagonal matrix B (_ 1>, symmetric tridiagonal matrix T And the eigenvector is calculated.
  • 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.
  • processing for calculating the eigenvector of the symmetric matrix A stored in the matrix storage unit 11 may be performed, but is omitted here.
  • 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.
  • step S601 and S602 the process of positively definite symmetric tridiagonal matrix T It is executed (steps S601 and S602).
  • the symmetric tridiagonal matrix T and the eigenvalues after the positive definite value ⁇ are used in the subsequent processing.
  • the diagonal matrix storage unit 13 stores the following symmetric tridiagonal matrix T of m X m, which is assumed to be:
  • the eigenvalue storage unit 18 stores an eigenvalue (l ⁇ i ⁇ m).
  • each element of the matrix B 1 ⁇ 2) can be calculated sequentially using each element of the matrix T.
  • [Equation 46] i ..
  • the eigenvector can be calculated by solving this simple equation.
  • T-lambda to Cholesky decomposition to I corresponds to ⁇ (0) ⁇ q (0 ), e (
  • each element of the diagonal matrix T is converted into each element of the matrix B (G) , and the matrix B
  • the force qd type conversion described for the standard qd type conversion may be a differen tial qd type conversion.
  • ⁇ ( ⁇ ) may be calculated again after resetting.
  • reverse Miura conversion may be called Miura conversion
  • Miura conversion may be called reverse Miura conversion
  • StdLVv conversion may be called stLVv conversion
  • RdLVv conversion may be called rLVv conversion.
  • 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 ⁇ (0) - ⁇ .
  • Step S903 The Cholesky decomposition portion 31 executes the processing of Procedure 1 described later.
  • 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.
  • 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
  • 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 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.
  • 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:
  • ⁇ ( ⁇ 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).
  • a constant for example, a real number such as 2, 3, 4 is preferable to select a value for reducing calculation error.
  • 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.
  • 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.
  • Step S1204 The Cholesky decomposition portion 31 sets ⁇ 2 to q (0) / (1 + ⁇ (0) u (0) ).
  • Step S 1205 The Cholesky decomposition portion 31 sets u ⁇ (°) u (°) to j81.
  • 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
  • ⁇ 1 becomes u (°), so when Miura transformation is performed on e (°) ⁇ ⁇ 8 1, a 1 is equal to ⁇ (0) u (°)
  • 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 ().
  • 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.
  • Step S1210 The Cholesky decomposition section 31 calculates ⁇ IX
  • a 1 is equal to ⁇ (0) u () and ⁇ 2 is q (0) / (1+ ⁇ () u
  • Step S1211 The Cholesky decomposition portion 31 sets
  • Step S1212 The Cholesky decomposition portion 31 sets the ⁇ to beta2-1 [delta] W.
  • ⁇ 82 is equal to q (°) ⁇ (1+ ⁇ ((°)), so
  • 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.
  • Step S1214 The Cholesky decomposition portion 31 calculates ⁇ 2 ⁇
  • 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.
  • 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 1 + ⁇ (+1) u.
  • Step S1302 The Cholesky decomposition portion 31 converts ⁇ 1 to u ⁇ (°) u ⁇ (+1
  • ⁇ 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
  • 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.
  • Step S1306 The Cholesky decomposition portion 31 sets j82 to u ⁇ () u ()) ⁇ 1
  • Step S 1307) The Cholesky decomposition portion 31 calculates ⁇ 1 X ⁇ 2. As mentioned above, a
  • 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.
  • j8 2 is equal to S (+1) u
  • 2 is 1 + ⁇ (+1) u (+1) . Therefore, to 0; 2
  • 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.
  • ⁇ 1 is equal to u and 2 is equal to S (+1) u, so the inverse Miura transform is applied to
  • the Cholesky decomposition unit 31 calculates ⁇ ⁇ ⁇ ⁇ 2.
  • Step S1311 The Cholesky decomposition portion 31 converts
  • 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.
  • 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 (+1) and e ends.
  • 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 (0) (1 + ⁇ (.) U (.))
  • Step S1402 The Cholesky decomposition portion 31 sets the counter k to “m ⁇ 1”.
  • Step S1403 The Cholesky decomposition portion 31 sets ⁇ 1 to
  • 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.
  • Step S1405 The Cholesky decomposition portion 31 sets j82 to u ⁇ (°) u (°)) ⁇ 1
  • 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 + ⁇ ( _1 ) u.
  • 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
  • Reskey decomposition unit 31 calculates ⁇ 1 ⁇ 2 and sets it to q.
  • Step S1409 The Cholesky decomposition portion 31 converts ⁇ 1 to u (0) (1 + ⁇ (0) u (0) ) / ⁇ 2
  • Step S 1410 The Cholesky decomposition portion 31 calculates j8 1 X j8 2. As mentioned above, ⁇
  • the Cholesky decomposition unit 31 calculates ⁇ ⁇ ⁇ ⁇ 2 k
  • 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.
  • 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
  • These values may be temporarily stored in a memory or the like (not shown) included in the Cholesky decomposition unit 31.
  • step 31 the process of the flowchart in FIG. 16 is repeated by the number of eigenvalues.
  • An array for the auxiliary variable u (n) is not necessarily prepared to reduce memory consumption k
  • 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.
  • 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.
  • the LV-type twist decomposition method has a parameter ⁇ (0) 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
  • This transformation is also a stqds transformation
  • the calculation amount is reduced, and the calculation speed can be improved.
  • 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
  • the orthogonal matrix Q can be calculated.
  • eigenvalue decomposition of the symmetric tridiagonal matrix T can be performed, eigenvalue decomposition of the symmetric matrix A is also performed as follows.
  • 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.
  • 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.
  • the case of using the LV type twist decomposition method will be described.
  • the Cholesky decomposition unit 31 first determines q (0) and e (0) 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 (0) , u (
  • step S1501 1 2 ⁇ , etc. are obtained sequentially (step S1501).
  • the Cholesky decomposition unit 31 uses the u ( ⁇ obtained by the Miura transform to convert the stdLVv transform to k
  • step S1502 u and the like are sequentially obtained.
  • the Cholesky decomposition unit 31 uses the u ((>) obtained by the Miura transform to execute the rdLVv transform to
  • step S1503 U ⁇ ⁇ , etc. are sequentially obtained.
  • 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
  • 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.
  • 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).
  • 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.
  • the matrix storage unit 11 may store the eigenvector of the symmetric matrix A! / ⁇ . The calculation method is as described above.
  • 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.
  • the output from the output unit 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.
  • 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.
  • 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.
  • the eigenvalue calculation unit 17 may calculate at least one eigenvalue of the symmetric tridiagonal matrix T.
  • the processing load can be reduced because it is not necessary to perform extra processing associated with calculating unnecessary eigenvalues.
  • the calculation time spent calculating the eigenvalues of the symmetric tridiagonal matrix T can be approximately (1 + kZm) Z2 times. sell.
  • m is the matrix size
  • k is the number of eigenvalues to be obtained.
  • the eigenvector calculation unit 19 calculates all eigenvectors corresponding to the calculated eigenvalue. 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.
  • 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.
  • 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.
  • the eigenvalue decomposition means 15a performs eigenvalue decomposition of the matrices ⁇ , ⁇ , ⁇ , ⁇ , and the eigenvalue decomposition means 15b
  • 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.
  • the eigenvalue calculating means 17a executes processing for calculating eigenvalues of the rows ⁇ , ⁇ , ⁇ , T, etc.
  • Processing for calculating eigenvalues and the like of the columns ⁇ 1, ⁇ 2, T may be executed.
  • 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.
  • the Cholesky decomposition means 22a may execute the process of performing Cholesky decomposition on half of the eigenvalues
  • the Cholesky decomposition means 22b may execute the process of performing Cholesky decomposition on the other half of the eigenvalues. ,.
  • 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.
  • the vector calculating means 23a executes processing for calculating eigenvectors for the eigenvalues subjected to Cholesky decomposition by the Cholesky decomposition means 22a
  • 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.
  • 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.
  • 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.
  • 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
  • the process of calculating the matrix elements may be executed in parallel. The process will be described below.
  • Equation 1 the processing for obtaining each eigenvalue can be executed in parallel.
  • 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.
  • each of the eigenvalue calculating means 17a and 17b has values of components of a matrix D and a vector z.
  • Equation 2 It can also be seen from Equation 2 that the processing for obtaining each eigenvector q can be executed in parallel.
  • the eigenvalue calculating unit 17a may calculate a half eigenvector q corresponding to the calculated eigenvalue
  • the eigenvalue calculating unit 17b may calculate the other half eigenvector q corresponding to the calculated eigenvalue.
  • 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.
  • 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.
  • the z component does not contain 0, and the diagonal component of D does not contain the same value.
  • 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.
  • 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.
  • each means may use a different memory.
  • the eigenvalue decomposition unit 15, the eigenvalue calculation unit 17 and the like have two hands using FIGS.
  • the eigenvalue decomposition unit 15, the eigenvalue calculation unit 17, and the like may include three or more means and execute parallel processing.
  • 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.
  • the eigenvalue decomposition unit 15 may not perform parallel processing.
  • the parallel processing may be executed by two or more devices, or may be executed by using two or more CPUs in one device.
  • 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.
  • 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.
  • 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.
  • 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.
  • 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.
  • the QR method DSTEQR provided by LAPACK was used as the QR method.
  • 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. .
  • 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.
  • 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.
  • 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.
  • the method of the eigenvalue decomposition apparatus 1 can perform eigenvalue decomposition at high speed, and has high eigenvalue and eigenvector accuracy.
  • parallelism is also excellent, it is possible to achieve higher speed by performing parallel processing.
  • Face image recognition processing by principal component analysis is performed by the following steps. (1) Creating a projection matrix
  • 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.
  • M face image data are handled.
  • the covariance matrix C is calculated as follows.
  • 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.
  • the eigenvalue is the maximum eigenvalue.
  • the projection matrix r? Is calculated as follows using d eigenvectors with the corresponding eigenvalues from the top.
  • the subscript k is a subscript for identifying an individual.
  • m d-dimensional vectors ⁇ are calculated as shown below, and the average value is calculated. Then, the average value / is registered.
  • the d-dimensional vector is calculated by multiplying the vector X of the face image data by the projection matrix 7 ?.
  • 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.
  • 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.
  • 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.
  • an optical reading device such as a scanner, a digital still camera, or a digital video camera.
  • the two-dimensional image to be handled is preferably a weak center projected image. At this time, the following equation holds.
  • step S5001 As can be seen from the shapes of M and S, the rank of D is 3.
  • step S5001 the rank of D is 3.
  • the matrix D is given.
  • the data M and shape S related to rotation are obtained below.
  • ⁇ is a singular value arranged in diagonal order in the order of magnitude
  • U and V are a left orthogonal matrix and a right orthogonal matrix, respectively.
  • 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.
  • the eigenbeta is equal to the right singular vector X that forms each column of the right orthogonal matrix V.
  • the left singular vector y can be obtained.
  • U (y, y, ..., y)
  • the left orthogonal matrix U can be calculated.
  • 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.
  • M ' L' ( ⁇ ') 1/2
  • S' ( ⁇ ') 1/2 R' T
  • the document 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.
  • the documents to be searched are d 1, d 2,.
  • Torr This is called a document vector.
  • 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.
  • the matrix D is called an index word document matrix.
  • each column of the index word document matrix is a document vector that represents information about the document.
  • 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.
  • FIG. 33 is a flowchart showing an example of a document search method using the eigenvalue decomposition apparatus 1 according to this embodiment.
  • step S6001,! / The question vector q is received.
  • 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.
  • D is used instead of D.
  • Latent semantic indexing 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.
  • 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.
  • 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 ".”
  • 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.
  • step S6002 select! / And k.
  • 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.
  • step S6003 singular value decomposition of the matrix D is performed.
  • 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.
  • 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.
  • the matrix D Singular value decomposition can be performed.
  • 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.
  • 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.
  • ⁇ is a k X k diagonal matrix where only the first k singular values are also forces.
  • step 6004 the similarity between the matrix D and the question vector q is calculated. V ⁇
  • the j-th document vector of D is represented by D e
  • T q is the coordinate (k-dimensional representation of the search question) under the base U of the search question vector q.
  • the eigenvalue decomposition apparatus 1 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.
  • the eigenvalue decomposition apparatus 1 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.
  • the process of calculating the eigenvalue force eigenbeta can be performed at high speed.
  • parallelization at the stage of calculating eigenvalues is difficult in the QR method.
  • each of the stage of calculating the eigenvalue and the stage of calculating the eigenvalue force eigenvector are performed. And essentially high parallelism.
  • the eigenvalue decomposition apparatus 1 according to the present embodiment can obtain higher accuracy than other methods.
  • the eigenvalue decomposition apparatus 1 may be a stand-alone apparatus or a server apparatus in a server client system, as described above. It may be a system.
  • 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.
  • 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.
  • 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.
  • each component may be configured by dedicated hardware, or a component that can be realized by software may be realized by executing a program.
  • 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.
  • the eigenvector may be calculated by the qd-type twist resolution method in the eigenvector calculation step.
  • 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.
  • the symmetric tridiagonal matrix T is converted into an upper bidiagonal matrix B (+1) and a lower bidiagonal matrix B (_1).
  • a Cholesky decomposition step for decomposing, and each element of the upper double diagonal matrix B (+1) and the lower double diagonal matrix B (_ 1), a positive value before Symbol symmetric tridiagonal matrix T
  • a vector calculation step of calculating an eigenvector using the eigenvalue and storing the eigenvector in the eigenvector storage unit.
  • the eigenvector may be calculated by the LV type twist decomposition method in the eigenvector calculation step.
  • 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,
  • the symmetric tridiagonal matrix T is transformed into an upper bidiagonal matrix B (+1) and lower 2 a Cholesky decomposing Cholesky decomposition step heavy diagonal matrix B (_1), and each element of the upper double diagonal matrix B (+1) and the lower double diagonal matrix B (_1), the symmetry 3
  • 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.
  • the step of accumulating information is performed only by hardware and does not include at least processing! /.
  • 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.
  • 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.
  • 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.
  • CD-ROM Compact Disk Read Only Memory
  • FD Flexible Disk
  • FIG. 35 is a diagram showing a computer system.
  • 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 computer 101 may include a network card (not shown) that provides connection to the LAN.
  • the function of the eigenvalue decomposition apparatus 1 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.
  • 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.
  • 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.
  • 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.
  • 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.
  • 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.
  • 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.
  • 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.
  • ⁇ 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.
  • 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

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Computing Systems (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Complex Calculations (AREA)

Description

明 細 書
固有値分解装置、及び固有値分解方法
技術分野
[0001] 本発明は、固有値分解を行う固有値分解装置等に関する。
背景技術
[0002] 固有値分解は、例えば、分子軌道法による化学計算、統計計算、情報検索等の非 常に広い分野で利用されている。また、固有値分解の方法も知られている（例えば、 非特許文献 1参照)。
非特許文献 1： J. J. M. Cuppen, 「A divide and conquer method for the symmetric tridiagonal eigenproblemj、 Numerische Mathematik, Vol. 3 6、 p. 177—195、 1981年
発明の開示
発明が解決しょうとする課題
[0003] 近年、これらの応用分野におけるデータ量の増大などに伴って、高速'高精度な固 有値分解が求められている。また、固有値分解を並列処理することができる固有値分 解装置等の開発が望まれていた。
[0004] 本発明は、上記状況の下になされたものであり、並列性に優れた高速、高精度な 固有値分解を実行可能な固有値分解装置等を提供することを目的とする。
課題を解決するための手段
[0005] 上記目的を達成するため、本発明による固有値分解装置は、対称 3重対角行列 T が記憶される対角行列記憶部と、前記対角行列記憶部から前記対称 3重対角行列 T を読み出し、当該対称 3重対角行列 Tを 2個の対称 3重対角行列に分割して前記対 角行列記憶部に蓄積し、その対称 3重対角行列を 2個の対称 3重対角行列に分割し て前記対角行列記憶部に蓄積する処理を、分割後の各対称 3重対角行列があらか じめ決められた大きさ以下となるまで繰り返す行列分割部と、前記あらかじめ決めら れた大きさ以下の各対称 3重対角行列の固有値と、当該各対称 3重対角行列の固有 ベクトル力 なる直交行列の一部の要素である行列要素とが記憶される固有値分解 記憶部と、前記あら力 め決められた大きさ以下の各対称 3重対角行列を前記対角 行列記憶部から読み出し、前記各対称 3重対角行列に対して固有値分解を行い、前 記各対称 3重対角行列の固有値と、前記各対称 3重対角行列の固有ベクトルからな る直交行列の行列要素とを少なくとも算出し、当該固有値と当該行列要素とを前記 固有値分解記憶部に蓄積する固有値分解部と、前記対称 3重対角行列 Tの固有値 が記憶される固有値記憶部と、各対称 3重対角行列の固有値と、行列要素とを前記 固有値分解記憶部から読み出し、前記固有値と、前記行列要素とから、分割元の対 称 3重対角行列の固有値と、分割元の対称 3重対角行列の行列要素とを算出して前 記固有値分解記憶部に蓄積し、その分割元の対称 3重対角行列の固有値と、行列 要素とを算出する処理を、対称 3重対角行列 Tの少なくとも 1個の固有値を算出する まで繰り返し、前記対称 3重対角行列 Tの少なくとも 1個の固有値を前記固有値記憶 部に蓄積する固有値算出部と、前記対称 3重対角行列 Tの固有ベクトルが記憶され る固有ベクトル記憶部と、前記対角行列記憶部から前記対称 3重対角行列 Tを読み 出し、前記固有値記憶部から前記対称 3重対角行列 Tの固有値を読み出し、前記対 称 3重対角行列 Tとその固有値とから、ツイスト分解法を用いて前記対称 3重対角行 列 Tの少なくとも 1個の固有ベクトルを算出して前記固有ベクトル記憶部に蓄積する 固有ベクトル算出部と、を備えたものである。
[0006] このような構成により、まず固有値を算出し、その固有値に基づいてツイスト分解法 を用いて固有べ外ルを算出することによって、高速で高精度な固有値分解を実現 することができうる。また、並列性にも優れている。さらに、全ての固有値や固有べタト ルを算出する必要がない場合には、必要な範囲で固有値や固有ベクトルを算出する ことができ、処理負荷を軽減することができる。
[0007] また、本発明による固有値分解装置では、前記固有ベクトル算出部は、 qd型ッイス ト分解法により固有ベクトルを算出してもよい。
このような構成により、高速で高精度な固有値分解を実現することができうる。また、 並列性にも優れている。
[0008] また、本発明による固有値分解装置では、前記固有ベクトル算出部は、前記固有 値記憶部から前記対称 3重対角行列 Tの固有値を読み出し、当該固有値のいずれ 力が負の値である場合に、前記対角行列記憶部力も前記対称 3重対角行列 Tを読 み出し、当該対称 3重対角行列 Τの固有ベクトルを変化させることなぐ当該対称 3重 対角行列 Τのすベての固有値が正の値となるように、当該対称 3重対角行列 Τを正 定値化し、正定値ィ匕した後の固有値を前記固有値記憶部に蓄積し、正定値化した 後の対称 3重対角行列 Τを前記対角行列記憶部に蓄積する正定値化部と、前記対 角行列記憶部からすべての固有値が正の値である対称 3重対角行列 Τを読み出し、 前記固有値記憶部から前記対称 3重対角行列 Τの正の値の固有値を読み出し、前 記対称 3重対角行列 Τの各要素に関して qd型変換を行うことによって、前記対称 3重 対角行列 Tを上 2重対角行列 B(+ 1)及び下 2重対角行列 B(_1)にコレスキー分解する コレスキー分解部と、前記上 2重対角行列 B(+1)及び下 2重対角行列 B(_1)の各要素 と、前記対称 3重対角行列 Tの正の値の固有値とを用いて固有ベクトルを算出して前 記固有ベクトル記憶部に蓄積するベクトル算出部と、を備えてもよい。
[0009] このような構成により、負の値の固有値が含まれる場合には、対称 3重対角行列 Tを 正定値ィ匕してから qd型ツイスト分解法を実行することができる。また、 LV型ツイスト分 解法よりも高速に計算することができうる。
[0010] また、本発明による固有値分解装置では、前記固有ベクトル算出部は、 LV型ツイ スト分解法により固有ベクトルを算出してもよい。
このような構成により、数値安定的に固有ベクトルの算出を行うことができうる。
[0011] また、本発明による固有値分解装置では、前記固有ベクトル算出部は、前記対角 行列記憶部から前記対称 3重対角行列 Tを読み出し、前記固有値記憶部から前記 対称 3重対角行列 Tの固有値を読み出し、前記対称 3重対角行列 Tの各要素に関し てミウラ変換、 dLVv型変換、逆ミウラ変換を行うことによって、前記対称 3重対角行列 Tを上 2重対角行列 B(+ 1)及び下 2重対角行列 B(_1)にコレスキー分解するコレスキー 分解部と、前記上 2重対角行列 B(+1)及び下 2重対角行列 B(_1)の各要素と、前記対 称 3重対角行列 Tの固有値とを用いて固有ベクトルを算出して前記固有ベクトル記憶 部に蓄積するベクトル算出部と、を備えてもよい。
このような構成により、コレスキー分解の処理において、複数の補助変数を用いるこ とによって、数値安定的に固有ベクトルの算出を行うことができうる。 [0012] また、本発明による固有値分解装置では、前記コレスキー分解部は、複数のコレス キー分解手段を備え、前記複数のコレスキー分解手段が、前記対称 3重対角行列 T をコレスキー分解する処理を並列実行してもよ ヽ。
このような構成により、コレスキー分解の処理を短時間で行うことができうる。
[0013] また、本発明による固有値分解装置では、前記ベクトル算出部は、複数のベクトル 算出手段を備え、前記複数のべ外ル算出手段が、固有べ外ルを算出する処理を 並列実行してもよい。
このような構成により、固有ベクトルを算出する処理を短時間で行うことができうる。
[0014] また、本発明による固有値分解装置では、前記固有値算出部は、複数の固有値算 出手段を備え、前記複数の固有値算出手段が、分割元の対称 3重対角行列の固有 値と、行列要素とを算出する処理を並列実行してもよい。
このような構成により、固有値を算出する処理を短時間で行うことができうる。
[0015] また、本発明による固有値分解装置では、前記固有値分解部は、複数の固有値分 解手段を備え、前記複数の固有値分解手段が、対称 3重対角行列に対して固有値 分解を行う処理を並列実行してもよ ヽ。
このような構成により、固有値分解する処理を短時間で行うことができうる。
[0016] また、本発明による固有値分解装置では、対称行列 Aが記憶される行列記憶部と、 前記対称行列 Aを前記行列記憶部から読み出し、前記対称行列 Aを 3重対角化した 前記対称 3重対角行列 Tを算出して前記対角行列記憶部に蓄積する対角化部と、を さらに備えてもよい。
[0017] このような構成により、任意の対称行列 Aについて固有値を算出することができる。
また、対称 3重対角行列 Tの固有ベクトルを用いることにより、対称行列 Aの固有べク トルを算出することもできうる。
[0018] また、本発明による固有値分解装置では、前記行列分割部は、対称 3重対角行列 を略半分の 2個の対称 3重対角行列に分割してもよい。
このような構成により、並列処理を適切に行うことができうる。
[0019] また、本発明による固有値分解装置では、前記固有値算出部は、対称 3重対角行 列 Tの全ての固有値を算出してもよ 、。 [0020] また、本発明による固有値分解装置では、前記固有ベクトル算出部は、対称 3重対 角行列 Tの全ての固有ベクトルを算出してもよい。
[0021] また、本発明による固有値分解装置では、前記対称行列 Aは、顔画像データを示 すベクトルの平均力 算出された共分散行列であってもよい。
このような構成により、固有値分解装置を顔画像データの認識の処理に用いること ができうる。
発明の効果
[0022] 本発明による固有値分解装置等によれば、高速で高精度な固有値分解を行うこと ができる。また、並列性にも優れている。
発明を実施するための最良の形態
[0023] 以下、本発明による固有値分解装置について、実施の形態を用いて説明する。な お、以下の実施の形態において、同じ符号を付した構成要素及びステップは同一ま たは相当するものであり、再度の説明を省略することがある。
[0024] (実施の形態 1)
本発明の実施の形態 1による固有値分解装置について、図面を参照しながら説明 する。
図 1は、本実施の形態による固有値分解装置 1の構成を示すブロック図である。図 1 において、本実施の形態による固有値分解装置 1は、行列記憶部 11と、対角化部 1 2と、対角行列記憶部 13と、行列分割部 14と、固有値分解部 15と、固有値分解記憶 部 16と、固有値算出部 17と、固有値記憶部 18と、固有ベクトル算出部 19と、固有べ タトル記憶部 20とを備える。
[0025] 行列記憶部 11では、任意の対称行列 Aが記憶される。この対称行列 Aは、各要素 が実数である実行列である。なお、対称行列 Aが記憶されているとは、対称行列 Aを 示すデータが記憶されている、という意味である。後述する記憶部においても同様で ある。行列記憶部 11は、所定の記録媒体 (例えば、半導体メモリや磁気ディスク、光 ディスクなど）によって実現されうる。行列記憶部 11での記憶は、 RAM等における一 時的な記憶でもよぐあるいは、長期的な記憶でもよい。行列記憶部 11に対称行列 A が記憶される過程は問わない。例えば、記録媒体を介して対称行列 Aが行列記憶部 11で記憶されるようになってもよぐ通信回線等を介して送信された対称行列 Aが行 列記憶部 11で記憶されるようになってもよぐあるいは、キーボードやマウス等の入力 デバイスを介して入力された対称行列 Aが行列記憶部 11で記憶されるようになって ちょい。
[0026] 対角化部 12は、対称行列 Aを行列記憶部 11から読み出し、その読み出した対称 行列 Aを 3重対角化した対称 3重対角行列 Tを算出する。そして、対角化部 12は、そ の算出した対称 3重対角行列 Tを対角行列記憶部 13に蓄積する。対角化部 12は、 例えば、ハウスホルダー（Householder)変換を必要なだけ繰り返し行う方法や、ラン チヨス (Lanczos)法を用いた方法、その他の 3重対角化法を用いて、対称行列 Aを 3 重対角化する。
[0027] 対角行列記憶部 13では、対称 3重対角行列 Tが記憶される。対角行列記憶部 13 は、所定の記録媒体 (例えば、半導体メモリや磁気ディスク、光ディスクなど）によって 実現されうる。対角行列記憶部 13での記憶は、 RAM等における一時的な記憶でも よぐあるいは、長期的な記憶でもよい。
[0028] 行列分割部 14は、対角行列記憶部 13から対称 3重対角行列 Tを読み出し、その 対称 3重対角行列 Tを 2個の対称 3重対角行列に分割して対角行列記憶部 13に蓄 積する。行列分割部 14は、その対称 3重対角行列を 2個の対称 3重対角行列に分割 して対角行列記憶部 13に蓄積する処理を、分割後の各対称 3重対角行列があらか じめ決められた大きさ以下となるまで再帰的に繰り返す。
[0029] 固有値分解部 15は、あらかじめ決められた大きさ以下の各対称 3重対角行列を対 角行列記憶部 13から読み出し、その各対称 3重対角行列に対して固有値分解を行 い、各対称 3重対角行列の固有値と、その各対称 3重対角行列の固有ベクトルを算 出する。固有値分解部 15は、例えば、 2分法と逆反復法を組み合わせた方法、 MR" 3法、 QR法等を用いて固有値分解を行ってもよい。ここで、 MR'3法による固有値分 解を行う場合には、 LAPACK3. 0 (http : //www. netlib. orgZlapackZ)で提 供されている DSTEGRを用いてもよい。また、 QR法による固有値分解を行う場合に は、 LAPACK3. 0で提供されている DSTEQRを用いてもよい。これらの固有値分 解の方法については、すでに公知であり、その詳細な説明を省略する。固有値分解 部 15は、算出した固有値を固有値分解記憶部 16に蓄積する。また、固有値分解部 15は、算出した固有ベクトル力 なる直交行列の一部の要素である行列要素も固有 値分解記憶部 16に蓄積する。固有ベクトル力もなる直交行列とは、固有ベクトルを各 列とする直交行列のことである。行列要素の詳細にっ 、ては後述する。
[0030] 固有値分解記憶部 16では、固有値分解部 15によって固有値分解された固有値と 、前述の行列要素とが記憶される。また、固有値分解記憶部 16では、固有値算出部 17によって算出された固有値や行列要素も記憶される。固有値分解記憶部 16は、 所定の記録媒体 (例えば、半導体メモリや磁気ディスク、光ディスクなど）によって実 現されうる。固有値分解記憶部 16での記憶は、 RAM等における一時的な記憶でも よぐあるいは、長期的な記憶でもよい。
[0031] 固有値算出部 17は、固有値分解部 15によって算出された固有値と、行列要素とを 固有値分解記憶部 16から読み出し、その固有値と行列要素とから、分割元の対称 3 重対角行列の固有値と、分割元の対称 3重対角行列の行列要素とを算出して固有 値分解記憶部 16に蓄積する。固有値算出部 17は、その分割元の対称 3重対角行 列の固有値と、行列要素とを算出する処理を、対称 3重対角行列 Tの固有値を算出 するまで再帰的に繰り返す。そして、固有値算出部 17は、算出した対称 3重対角行 列丁の固有値を固有値記憶部 18に蓄積する。
[0032] 固有値記憶部 18では、対称 3重対角行列 Tの固有値が記憶される。固有値記憶部 18は、所定の記録媒体 (例えば、半導体メモリや磁気ディスク、光ディスクなど）によ つて実現されうる。固有値記憶部 18での記憶は、 RAM等における一時的な記憶で もよぐあるいは、長期的な記憶でもよい。
[0033] 固有ベクトル算出部 19は、対角行列記憶部 13から対称 3重対角行列 Tを読み出し 、固有値記憶部 18から対称 3重対角行列 Tの固有値を読み出す。そして、固有べク トル算出部 19は、対称 3重対角行列 Tとその固有値とから、ツイスト分解法を用いて 対称 3重対角行列 Tの固有ベクトルを算出して固有ベクトル記憶部 20に蓄積する。 固有ベクトル算出部 19は、 LV型ツイスト分解法により固有ベクトルを算出してもよぐ あるいは、 qd型ツイスト分解法により固有ベクトルを算出してもよい。なお、具体的な ツイスト分解の手法にっ ヽては後述する。 [0034] 固有ベクトル記憶部 20では、対称 3重対角行列 Tの固有ベクトルが記憶される。固 有ベクトル記憶部 20は、所定の記録媒体 (例えば、半導体メモリや磁気ディスク、光 ディスクなど）によって実現されうる。固有ベクトル記憶部 20での記憶は、 RAM等に おける一時的な記憶でもよぐあるいは、長期的な記憶でもよい。
[0035] なお、行列記憶部 11、対角行列記憶部 13、固有値分解記憶部 16、固有値記憶部 18、固有ベクトル記憶部 20の任意の 2以上の記憶部は、同一の記録媒体によって 実現されてもよぐあるいは、別々の記録媒体によって実現されてもよい。前者の場合 には、例えば、対称行列 Aの記憶されている領域が行列記憶部 11となり、対称 3重 対角行列 T等の記憶されている領域が対角行列記憶部 13となる。
[0036] また、行列記憶部 11、対角行列記憶部 13、固有値分解記憶部 16、固有値記憶部 18、固有ベクトル記憶部 20の各記憶部は、 2以上の記録媒体から構成されてもよい
[0037] 次に、本実施の形態による固有値分解装置 1の動作について、図 2のフローチヤ一 トを用いて説明する。
(ステップ S101)対角化部 12は、行列記憶部 11で記憶されて!ヽる対称行列 Aを読 み出し、その行列 Aを 3重対角化して対称 3重対角行列 Tを算出して対角行列記憶 部 13に蓄積する。
[0038] (ステップ S102)行列分割部 14、固有値分解部 15、固有値算出部 17によって対 称 3重対角行列 Tの固有値が算出され、固有値記憶部 18に蓄積される。この処理の 詳細については後述する。
[0039] (ステップ S103)固有ベクトル算出部 19は、対角行列記憶部 13から対称 3重対角 行列 Tを読み出し、固有値記憶部 18から対称 3重対角行列 Tの固有値を読み出し、 対称 3重対角行列 Tの固有ベクトルを算出して固有ベクトル記憶部 20に蓄積する。こ の処理の詳細については後述する。
[0040] このようにして、対称 3重対角行列 Tの固有値分解が終了する。ここで、対称行列 A の固有値は対称 3重対角行列丁の固有値に等 ヽため、対称行列 Aの固有値も算出 されたことになる。また、後述するように、所定の変換を行うことによって、対称行列 A の固有ベクトルも対称 3重対角行列 Tの固有ベクトル力 容易に算出することができ る。
[0041] 次に、図 2のフローチャートのステップ S101, S102の処理、すなわち、対称 3重対 角行列 Tの固有値を算出するまでの処理について、より詳細に説明する。まず、図 2 のフローチャートのステップ S102の処理について、図 3のフローチャートを用いて説 明する。
[0042] (ステップ S201)行列分割部 14は、対角行列記憶部 13から対称 3重対角行列 Tを 読み出し、その対称 3重対角行列 Tを 2個の対称 3重対角行列に分割して対角行列 記憶部 13に蓄積する。行列分割部 14は、その対称 3重対角行列を 2個の対称 3重 対角行列に分割する処理を、分割後の各対称 3重対角行列があらかじめ決められた 大きさ以下となるまで繰り返す。この処理の詳細については後述する。
[0043] (ステップ S202)固有値分解部 15は、対角行列記憶部 13で記憶されているあらか じめ決められた大きさ以下の各対称 3重対角行列について、固有値分解を行う。固 有値分解部 15は、固有値分解の結果である固有値と、固有ベクトルカゝらなる直交行 列の一部の要素である行列要素とを固有値分解記憶部 16に蓄積する。
[0044] (ステップ S203)固有値算出部 17は、対称 3重対角行列の固有値と、行列要素と を固有値分解記憶部 16から読み出し、分割元の対称 3重対角行列の固有値と、分 割元の対称 3重対角行列の行列要素とを算出して固有値分解記憶部 16に蓄積する 。固有値算出部 17は、その分割元の対称 3重対角行列の固有値と、行列要素とを算 出する処理を、対称 3重対角行列 Tの固有値を算出するまで繰り返し、対称 3重対角 行列 Tの固有値を固有値記憶部 18に蓄積する。このようにして、固有値を算出する 処理が終了する。
[0045] [対称行列 Aの対角化]
対称行列 Aは、例えば、ノ、ウスホルダー変換等を用いて、以下に示されるように 3重 対角化することができることが知られている。ここで、対称 3重対角行列 Tは、行の数と 列の数とがー致する正方行列である。
[0046] [数 1]
Figure imgf000012_0001
ノ
ただし、 は《x«実対称行列、 Ρは直交行列である。
[0047] したがって、対角化部 12は、上記のようにして、行列記憶部 11から対称行列 Αを読 み出し、対称 3重対角行列 Tを算出することができる。その算出された対称 3重対角 行列 Tは、対角行列記憶部 13で記憶される (ステップ S101)。
[0048] [対称 3重対角行列 Tの分割]
対称 3重対角行列 Tを分割する処理について説明する。次式で示されるように、 n
Xnの対称 3重対角行列 Tが与えられた場合に、それらを 2個の対称 3重対角行列 T
, Tと、その他の要素（下記の Θと β)とに分けることができる。
[0049] [数 2]
Figure imgf000012_0002
Figure imgf000012_0003
ただし、 Γは《x« («は 2以上の整数） の実対称 3重対角行列、 7；は f^ x の 実対称 3重対角行列、 Γ2は (《 - ?¾) X の実対称 3重対角行列、
e; ≡(0,···,0,1)τ eR"\ ef ≡ (1,0,···,0)τ e R""\ β は Γの (?¾ +1,?¾)の要素、
Θ は非ゼロ定数である。
ここで、 ηは、 Κηく ηとなる整数である。例えば、 8 X 8の対称 3重対角行列 Τを半 分の大きさに分割する具体的な処理は、次のようになる。
[数 3]
k
h
t
T = ， in
^15ノ
Figure imgf000013_0001
[0051] ただし、 β =tであり、 Θは計算誤差が少なくなるように任意に設定することができ
8
る。また、上記各行列において、明記されている以外の行列の要素は 0である。なお 、上記各行列の明記されている各要素は 0以外であってもよぐ 0であってもよい。ここ で、並列処理を適切に実行するためには、 ηを、 ηΖ2を超えない最大の整数にとる 力 あるいは、 ηΖ2を下まわらない最小の整数にとることが好適である（この ηの値を 用いて行列を 2個の行列に分割する場合を、「行列を略半分の 2個の行列に分割す る」と呼ぶことにする)。
[0052] [数 4] ί¾ =L"/2」
本実施の形態では、 を上式のようにとるものとする。なお、 の値は、前述のよう に、 l〈nく nの範囲内で任意であることは言うまでもない。
[0053] したがって、行列分割部 14は、上述のようにして、対角行列記憶部 13が記憶して いる対称 3重対角行列 Tを読み出し、その対称 3重対角行列 Tを 2個の対称 3重対角 行列 Τ , Tと、 2個の値 |8、 Θに分割することができ、その分割の処理を繰り返すこと
1 2
ができる。行列分割部 14は、分割後の 2個の対称 3重対角行列 Τ , Tと、 2個の値 |8
1 2
、 Θとを、対角行列記憶部 13に蓄積する。なお、 2個の値 、 Θは、どの分割で発生 した値であるのかがわ力るように対角行列記憶部 13で記憶されることが好適である。 図 4は、図 3のフローチャートのステップ S201における行列分割部 14による行列を 分割する処理を示すフローチャートである。
[0054] (ステップ S301)行列分割部 14は、カウンタ Iを「1」に設定する。
(ステップ S302)行列分割部 14は、対角行列記憶部 13から I番目の分割をしてい ない対称 3重対角行列を読み出し、その対称 3重対角行列を 2個の対称 3重対角行 列と、その他の要素とに分割する。そして、行列分割部 14は、分割した 2個の対称 3 重対角行列と、その他の要素とを対角行列記憶部 13に蓄積する。
[0055] (ステップ S303)行列分割部 14は、 I番目の分割を行って 、な 、対称 3重対角行列 力 対角行列記憶部 13で記憶されているかどうか判断する。そして、 I番目の分割を 行って 、な 、対称 3重対角行列が、対角行列記憶部 13で記憶されて 、る場合には、 ステップ S302に戻り、そうでな ヽ場合に ίま、ステップ S304に進む。
[0056] (ステップ S304)行列分割部 14は、 I番目の分割を行った対称 3重対角行列の大き さがあら力じめ決められている大きさ以下力どうか判断する。行列分割部 14は、例え ば、目的とする行列の大きさ（例えば、 25 X 25など）の記憶されて ヽる図示しな 、記 録媒体から、その目的とする行列の大きさを読み出し、対角行列記憶部 13で記憶さ れている、 I番目の分割後の対称 3重対角行列がその大きさ以下であるかどうかを判 断してもよい。そして、 I番目の分割を行った対称 3重対角行列の大きさがあらかじめ 決められている大きさ以下である場合には、対称 3重対角行列を分割する処理は終 了となり、そうでない場合には、ステップ S305に進む。 (ステップ S305)行列分割部 14は、カウンタ Iを 1だけインクリメントする。そして、ス テツプ S302に戻る。
[0057] なお、この図 4のフローチャートでは、上述のように、行列分割部 14が各行列を、略 半分の 2個の行列に分割するため、 I番目の分割後の各対称 3重対角行列の大きさ がほとんど同じである場合について説明したが、行列分割部 14が各行列を、略半分 の 2個の行列に分割しな ヽ場合には、行列分割部 14は、分割後の各行列があらかじ め決められた大きさ以下となるように、分割を繰り返すものとする。
[0058] また、図 4のフローチャートでは、ステップ S304において、行列分割部 14が、対角 行列記憶部 13で記憶されている分割後の行列の大きさをあら力じめ決められた大き さと比較する処理を実行する場合について説明したが、これは一例であって、行列分 割部 14は、ステップ S304において、それ以外の処理を行ってもよい。例えば、上述 のように、行列分割部 14が各行列を、略半分の 2個の行列に分割する場合には、元 の対称 3重対角行列 Tの大きさを知ることができれば、何番目の分割で目的とする行 列の大きさとなるのかを知ることができる。したがって、 N番目（Nは 1以上の整数）の 分割で目的とする行列の大きさとなる場合には、ステップ S304において、 Iが Nであ るかどうかを比較し、 Nでなければステップ S305に進み、 Nであれば一連の処理を 終了するようにしてもよい。
[0059] 図 5は、行列の分割について説明するための図である。まず、行列分割部 14は、 1 番目の分割として、対称 3重対角行列 Tを、対称 3重対角行列 Tと、対称 3重対角行 列 Tとに分割する (ステップ S301, S302)。対称 3重対角行列 Tは、 1個しかないた
2
め、行列分割部 14は、 1番目の分割をして ヽな 、行列がな 、と判断する (ステップ S3 03)。また、対称 3重対角行列 T等はあらかじめ決められた大きさ以下の行列でない とすると (ステップ S304)、行列分割部 14は、 2番目の分割として、対称 3重対角行列 Tを、対称 3重対角行列 T と、対称 3重対角行列 T とに分割する (ステップ S305,
1 11 12
S302)。この場合には、 2番目の分割を行っていない対称 3重対角行列 Tが存在す
2 るため、行列分割部 14は、対称 3重対角行列 Tも、対称 3重対角行列 T と、対称 3
2 21 重対角行列 T とに分割する (ステップ S303, S302)。このようにして、分割後の各
22
対称 3重対角行列が目的とする大きさ以下になるまで、行列を分割する処理が繰り返 される。なお、図 5では、対称 3重対角行列以外の要素（上述の 13、 Θ )については省 略している。また、図 5において、各対称 3重対角行列は、正方行列である。
[0060] [分割された行列の固有値分解]
行列 Tが対称行列であるとすると、行列 Tの固有値分解は、次のようになる。
[数 5]
[0061] ここで、 Dは、行列 Tと同じ大きさの対角行列である。 Dの各対角成分は Tの固有 値である。また、 Qは、行列 Tの固有ベクトル力もなる直交行列である。
[0062] したがって、固有値分解部 15は、対角行列記憶部 13から分割後の各対称 3重対 角行列を読み出し、上記のようにして、固有値分解を行う (ステップ S 202)。固有値 分解の方法として、例えば、 2分法と逆反復法を組み合わせた方法、 MR" 3法、 QR 法等を用いてもよいことは前述の通りである。図 5で示されるように対称 3重対角行列 の分割が行われた場合には、固有値分解部 15は、各対称 3重対角行列 Τ , T ,
111 112
Τ , Τ , · ··, Τ に対して固有値分解を行うことになる。なお、固有値分解部 15
121 122 222
は、固有値分解の結果得られた各固有値と、固有値分解の結果得られた直交行列 の一部の要素である行列要素とを固有値分解記憶部 16に蓄積する。ここで、行列要 素力 直交行列のどの要素を含むのかについては後述する。なお、固有値分解部 1 5が固有値と行列要素とを固有値分解記憶部 16に蓄積するのではなぐ固有値と固 有ベクトル力もなる直交行列とを固有値分解記憶部 16に蓄積してもよい。その場合 でも、固有ベクトル力もなる直交行列を蓄積することによって、少なくとも、行列要素が 蓄積されたことになる。
[0063] なお、後述するように、固有値を算出するときに用いるために必要なのは、直交行 列 Qの全体ではなぐその直交行列 Qの第 1行と最終行である。したがって、この固 有値分解において、直交行列 Qの全体を求めるのではなぐ直交行列 Qの第 1行と 最終行を求めるようにしてもよい。そのためには、例えば、 LAPACK3. 0で提供され て！、る QR法の固有値分解ルーチンである DSTEQRにお!/、て、直交行列 Qの算出 における初期値として、単位行列 Iを用いる代わりに、次式で示される行列 Υを用いて ちょい。 [数 6]
Figure imgf000017_0001
ただし、 ρ,は η χ ηの直交行列であり、 ； Τは 2 χ ηの行列である。
[0064] このようにすることで、 DSTEQRで直交行列 Qの第 1行と最終行を求めるための計 算領域を削減することができる。直交行列 Qそのものを算出する場合には、 0 (η2)の 作業領域が必要であるが、上記の行列 Υを用いて直交行列 Qの第 1行と最終行を算 出する場合には、 Ο (η)の作業領域でよいことになる。その結果、計算量も少なくなり 、固有値分解の処理速度も向上する。このように、固有値分解部 15は、分割後の対 称 3重対角行列に対して、直交行列 Qそのものを算出する固有値分解を行ってもよく 、あるいは、直交行列の一部の成分 (ここでは、第 1行と最終行)を算出する処理を行 つてもよい。この後者の処理も固有値分解と呼ぶことにする。このように、固有値分解 部 15は、あらかじめ決められた大きさ以下の各対称 3重対角行列に対して固有値分 解を行い、各対称 3重対角行列の固有値と、各対称 3重対角行列の固有ベクトルか らなる直交行列の行列要素 (ここでは、第 1行と最終行)とを少なくとも算出し、その固 有値と、その行列要素とを固有値分解記憶部 16に蓄積してもよい。
[0065] [固有値の算出]
まず、図 6で示されるように、分割された 2個の対称 3重対角行列 Τ , Τから、分割
1 2
元の対称 3重対角行列 Τの固有値等を算出する処理について説明する。対称 3重
0
対角行列 Τ , Τは、次のように固有値分解されていたとする。ここで、対称 3重対角
1 2
行列 Τ , Τは、両者共に正方行列であるとする。
1 2
[数 7]
[0066] この場合に、分割元の対称 3重対角行列 Tは、次のようになる。
0
[数 8]
Figure imgf000018_0001
ただし、 は の最終行であり、 f2は ρ2の第 1行である。
上記の式において、 j8、 0は、対称 3重対角行列の分割の処理において説明した 値である。ここで、行列（D + p ζζτ)の固有値分解が得られれば、行列 Τの固有値
12 0 分解が得られたことになる。なぜなら、行列（D + ρ ζζΤ)
12 の固有値分解 D + ρ ζζΤ
12
= QDQTが求まると、次式のようになり、行列 Tの固有値分解が完了するからである
Figure imgf000018_0002
[0068] 次に、行列（D + p ZZ T)の固有値を算出する方法について説明する。行列（D
+ ρ ζζ )の固有値をえ、固有ベクトルを q.とすると、
[数 10]
( + pzz T)q, =へ q,
となる。
[0069] (l) zの成分に 0が含まれず、 D の対角成分に同じ値が含まれない場合
12
まず、 zの成分に 0が含まれず、 D の対角成分に同じ値が含まれない場合につい
12
て説明する。その場合には、上記の式を、次に示すように順次、変形することができ る。
[数 11]
Figure imgf000018_0003
1 + ρ zT(D12 — λ,·/)- = 0
[0070] ここで、
[数 12] D12 = dia ic/j ,d2,---,dn)
z =(ζ1,ζ2,···,ζπ)τ
とすると、次の固有方程式を得ることができる。固有値え ^=1, 2, ···, n)は、その固 有方程式 n個の解である。下記の f ( λ )は、単調増加であり、かつ f ( λ ) =0の解の範 囲が既知であるため、ニュートン法を用いて計算することができる。
[数 13] (λ) = ι+ Ρ =0 (式
[0071] なお、 D = diag( , λ , ···, λ )となる。行列 Qは、行列（D + P ZZ )の固有べ
1 2 n 12
クトルカ なる直交行列であるため、 Q= (q , q , ···, q )となる。したがって、固有べ
1 2 n
タトル qを求めると、行列 Qを求めることができる。なお、固有ベクトル q.は、次式で与 えられる。
[数 14]
q' (式 2
Figure imgf000019_0001
[0072] (2) zの成分に 0が含まれる場合
zの成分に 0が含まれる場合には、行列の次数を下げることができ、その次数を下げ た行列を用いて固有値、及び固有ベクトルを求めることができる。この行列の次数を 下げる操作をデフレーションと呼ぶ。
[0073] z=0の場合には、次式で示すように、 (d, e)がー組の固有値と固有ベクトルとなる 。ただし、 eは、 j番目の成分が 1であり、その他の成分が 0である単位ベクトルである。
[数 15]
(( + ρζζΊ )-dJI)eJ = (D12 -dJI)eJ + pz(zT ）
=0
[0074] また、行列（D + ρ ζζ1)の固有方程式は、次式のようになる c
12
[数 16] \pu + pzz T - In \ = ( j - λ; )|θ12 + pzz T - ただし、 f)u + pii Tは、 Du + pzz τの第 j行、 第 j列を除いた小行列である。 また、 /„は、 n次の単位行列である。
[0075] したがって、 z =0の場合の固有値は、 dと、行列 (D + p ζζτ)の第 j行、第 j列を除 j j 12
いた小行列について、上記式 1の固有方程式を解いた値である。すなわち、固有値 は、 dと、次の固有方程式を解いた値となる。なお、次の固有方程式を解いた固有値 を、固有値え (i= l, 2, · ··, j - 1, j + 1, · ··, n)とする。
[数 17]
： 0
Figure imgf000020_0001
固有ベクトルは、行列（D + p zz1)の第 j行、第 j列を除 、た小行列にっ 、て上記
12
式 2で算出した n—l次元の固有ベクトルの第 (j 1)行と、第 j行との間に 0を挿入し た次式の n次元のベクトルと、上述の eとなる。
[数 18] z1 ∑2
d, - λ, d ― Λ, —へ d -λ, ' d„ -λ,
Figure imgf000020_0002
[0077] このように、 zの成分に 0が含まれる場合には、行列（D + p ζζτ)よりも次元の小さ
12
い行列について固有値分解を行うことになるため、計算を高速ィ匕することが可能であ り、結果として、計算時間が短縮されうる。
[0078] (3) D の対角成分に同じ値が含まれる場合
12
d = d 二… = dの場合にも、上記説明と同様のデフレーションを行うことによつ j + 1 j + 2 k
て、固有値と固有ベクトルとを算出することができる。
[0079] [数 19] τ = ^j+i , zj+i 7…， zt )
に対して、次の 2個の直交行列を考える。 Q ζ = (*,0,0,···,0)τ
Figure imgf000021_0001
ただし、 Qe R ( ) )である c
また、 *は所定の実数である。
[0080] すると、 d :dより、
j
[数 21]
QD12QT=D12 となる。したがって、次式が得られる。
[数 22]
Figure imgf000021_0002
[0081] ここで、
[数 23]
D --— · X
12 + pzz
の固有値と固有ベクトルとは、上記の zの成分に 0が含まれるデフレーションの場合と 同様にして求めることができる。また、
[数 24]
Q
は、符号の任意性を除いて一意に定まる。ここで、その算出方法について説明する。 まず、次式を満たす直交行列について考える。
[0082] [数 25]
Q , '( +い…， ，…， ）— =(*,z J+2,… ; "O w,…， ）
この直交行列を次式のようにとると、上記式を満たすことになる。
[数 26]
Figure imgf000022_0001
[0083] ここで、上記の式の *の値は、次式から求まる
[数 27]
Q l++ι,ί■ j ,++ι ¾ )
Figure imgf000022_0002
[0084] したがって、
[数 28]
Q一 ^j ' Q
となる。その結果、
[数 29]
Q
も、符号の任意性を除いて求めることができる。このようにして、 d =d =〜 = d
j +1 j + 2 k の場合にも行列（D— + pzz"1)の固有値と固有ベクトルとを求めることができる。すな
12
わち、固有値は、 d , d , ·, dと、次の固有方程式を解いた値となる。なお、次 j+2 j+3
の固有方程式を解いた固有値を、固有値え (i -, 2,…， j + k- ·, n)とす る。
[数 30] (λ) = ι+Ρ ∑ -Γ ^ = °
m-1 —人 固有ベクトルは、次のようになる。
[数 31] ρτ + l≤i≤n)
Figure imgf000023_0001
Le u + 2≤i≤k
[0086] なお、上記（2) zの成分に 0が含まれる場合、（3) D の対角成分に同じ値が含まれ
12
る場合において、 zの成分力 SOであるかどうかの判断や、 D の対角成分に同じ値が
12
含まれているかどうかの判断は、微少な誤差を考慮して行ってもよい。例えば、次式 が満たされる場合には、 zの成分が 0であると判断してもよぐまた、 D の対角成分が
12
同じ値であると判断してもよ!、。
[数 32]
Figure imgf000023_0002
たたし、 rQLj
Figure imgf000023_0003
rQL ≡maxf<i1，(i2，- - -, <i.- ) x e ?iiD 。
[0087] なお、 εは、マシン'ィプシロンに定数 (例えば、 2, 4, 8などの実数であり、計算誤 差を少なくするための値を選択することが好適である）を掛けた値を用いてもよぐ他 の適切な値であってもよ ヽ。
[0088] このようにして、（1) ζの成分に 0が含まれず、 D の対角成分に同じ値が含まれない
12
場合、（2) ζの成分に 0が含まれる場合、 (3) D の対角成分に同じ値が含まれる場合
12
のそれぞれについて、固有値分解 D + p ζζτ二 QDQTが求まることになる。したがつ
12
て、行列 Qと行列 Q から行列 Qを求めることができ、行列 Dも求まっているため、分
12 0
割元の対称 3重対角行列 Tの固有値分解を行うことができる。このように、対象となる
0
行列を 2個の副行列に分割し、その分割した後の各副行列について固有値分解を 行 、、その後に各副行列の固有値分解を用 、て分割元の行列の固有値分解を行う 方法は、分割統治法（Divide and Conquer : D&C)と呼ばれている。
[0089] 分割統治法では、固有ベクトルまで算出することになり、その固有ベクトルを算出す る処理にお!、て行列の計算をしなければならな 、ため、非常に負荷の大き 、処理と なる。分割統治法において固有値分解をする場合には、例えば、計算時間の 95% 程度が固有ベクトルを算出するベクトル更新の処理 (例えば、上述のように、行列 Qと 行列 Q 12から行列 Q 0を算出するために行列を掛け合わせる処理）に費やされることも ありうる。
[0090] 一方、分割元の行列 Tの固有値のみを算出する場合には、上記処理を簡略化す
0
ることができる。次に、その方法について説明する。上記の結果より、次のようになる。
[0091] [数 33]
f = f122 (式 3 )
Ι = Ι12β (式 4 )
ただし、 fは ρ。の第 1行であり、 Iは ρ。の最終行である。
また、 f12は g12の第 1行であり、 f12 = ,(),·' ·, 0)となる。
また、 ι12は ρ12の最終行であり、 ι12 = (ο,- ·-,ο,ι2 )となる。
また、 f\は ¾の第 1行であり、 ι2は ρ2の最終行である。
[0092] ここで、上記のように、 fは行列 Τを固有値分解した結果の直交行列の最初の行で
0
あり、 1は行列 Tを固有値分解した結果の直交行列の最後の行である。この fと 1、すな
0
わち、直交行列の最初の行と、最後の行とが行列要素となる。なお、直交行列 Q 12に ついて、後述する列の入れ替えを行った場合には、 f 、1
12 12はそれぞれ、その列を入 れ替えた後の直交行列 Q の第 1行、最終行となる。
12
[0093] 上記の結果から、図 6で示されるように、行列 Τ , Tから分割元の行列 Tを構成す
1 2 0 る場合に、行列 τについて、
0
[数 34]
λ f, i
をそれぞれ算出することができる。このような処理を繰り返すことにより、最終的に、対 称行列 Aを 3重対角化した対称 3重対角行列 Tの固有値を算出することができる。
[0094] なお、上記の（3) D の対角成分に同じ値が含まれる場合には、式 3,式 4が次のよ
12
うになる。
[数 35] f = f12Q = f12Q TQ
Figure imgf000024_0001
ただし、 = QQである。
[0095] ここで、上記式を計算する場合に、 [数 36]
を先に計算して、その後に
[数 37]
Q
を掛けた方が、計算量が少なくなるため、その順番で計算するようにしてもよい。
[0096] したがって、固有値算出部 17は、上述のようにして、固有値分解記憶部 16から固 有値と行列要素とを読み出し、対角行列記憶部 13から行列の分割時に発生したそ の他の要素に関する β、 Θを読み出し、それらを用いることによって、分割元の対称 3重対角行列の固有値と行列要素とを算出する。そして、それらを算出する処理を繰 り返すことによって、最終的に対称 3重対角行列 Τの固有値を算出する。図 7は、図 3 のフローチャートのステップ S203における固有値算出部 17が固有値を算出する処 理を示すフローチャートである。
[0097] (ステップ S401)固有値算出部 17は、カウンタ Jを「1」に設定する。
(ステップ S402)固有値算出部 17は、 J番目の固有値の算出は最後の固有値の算 出であるかどうか判断する。ここで、最後の固有値の算出とは、対称 3重対角行列 T の固有値を算出することである。そして、最後の固有値の算出である場合には、ステ ップ S406に進み、そうでな ヽ場合に ίま、ステップ S403に進む。
[0098] (ステップ S403)固有値算出部 17は、分割元の対称 3重対角行列の固有値等を算 出する。この処理の詳細については後述する。
(ステップ S404)固有値算出部 17は、 J番目の固有値の算出において、分割元の 全ての対称 3重対角行列の固有値等を算出したかどうか判断する。そして、 J番目の 固有値の算出において、分割元の全ての対称 3重対角行列の固有値等を算出した 場合には、ステップ S405〖こ進み、そうでない場合には、ステップ S403〖こ戻る。
[0099] (ステップ S405)固有値算出部 17は、カウンタ Jを 1だけインクリメントする。そして、 ステップ S402に戻る。
(ステップ S406)固有値算出部 17は、対称 3重対角行列 Tの固有値を算出し、その 算出した固有値を固有値記憶部 18に蓄積する。このようにして、対称 3重対角行列 T の固有値を算出する一連の処理は終了となる。
[0100] 図 8は、図 7のフローチャートにおけるステップ S403の詳細な処理を示すフローチ ヤートである。なお、この図 8のフローチャートの処理を実行する前に、対角行列 D
12 の対角成分が昇順、あるいは降順となるように並べ替えておくことが好適である。なお 、その場合には、直交行列 Q の列、及び zの成分も、対角行列 D の並び替えと同
12 12
様に並び替えるものとする。具体的には、対角行列 D の第 i成分と、第 j成分とを入
12
れ替えた場合には、直交行列 Q の第洌と、第 j列とを入れ替え、 zの第珩と、第 j行
12
とを入れ替えればよい。
[0101] (ステップ S501)固有値算出部 17は、分割元の行列の固有値を、式 1等を用いて 算出する。そして、固有値算出部 17は、算出した固有値を固有値分解記憶部 16〖こ 蓄積する。
[0102] (ステップ S502)固有値算出部 17は、ステップ S501で算出した固有値を用いて、 式 2等力も q.を算出する。この q.を算出することにより、 q.を列ベクトルに有する行列 Q を算出したことになる。固有値算出部 17は、算出した行列 Qを図示しない記録媒体 にお 、て一次記憶してもよ 、。
[0103] (ステップ S503)固有値算出部 17は、ステップ S502で算出した行列 Qを用いて、 分割元の行列に関する行列要素 f、 1を算出する。この行列要素は、式 3,式 4で示さ れるものである。そして、固有値算出部 17は、算出した分割元の行列要素を固有値 分解記憶部 16に蓄積する。このようにしてステップ S403の処理は終了となる。
[0104] なお、図 8のフローチャートにおいて、 zの成分に 0が含まれず、 D の対角成分に
12
同じ値が含まれない場合には、上記（1)で説明したようにして、固有値、固有ベクトル 力 なる行列 Qを算出することができる。また、 zの成分に 0が含まれ、 D の対角成分
12
に同じ値が含まれない場合には、上記（2)で説明したようにして、固有値、固有べタト ルカもなる行列 Qを算出することができる。また、 D の対角成分に同じ値が含まれる
12
場合には、 zの成分に 0が含まれるかどうかに関わらず、上記（3)で説明したようにし て、固有値、固有ベクトル力もなる行列 Qを算出することができる。上記（3)で説明し たようにして固有値等を算出する場合には、前述のように、行列 Qを算出するのでは なぐ行列 Qを構成する行列を直接 f 等に作用させることによって、 f、 1を算出するよ うにしてもよい。
[0105] また、図 8のフローチャートでは、行列 Qから、行列要素 f、 1を一度に算出する場合 について示した力 行列要素 f、 1をその成分ごとに計算してもよい。上記の式 3,式 4 力 明らかなように、 f 、 1 に固有ベクトル qを掛けたもの力 行列要素 f、 1の i番目の
12 12 i
成分となる。このような方法で行列要素 f、 1を計算する場合には、固有ベクトル固有べ タトル qのみを iの値を順次、変えながら記憶することによってその計算を行うことがで きる。したがって、行列要素 f、 1の計算において行列 Q全体を記憶する必要がなくなり 、その計算で用いるメモリ上の作業領域を削減することが可能である。
[0106] 図 9は、固有値を算出する処理について説明するための図である。固有値分解記 憶部 16では、固有値分解部 15によって行列 Τ , Τ , · ··, Τ が固有値分解され
111 112 222
た結果、すなわち、それらの固有値と、行列要素とが記憶されているものとする。まず 、固有値算出部 17は、固有値等の 1番目の算出を開始する (ステップ S401, S402) 。ここで、固有値等の 1番目の算出とは、図 9の一番下の行の行列の固有値と、行列 要素とから、下から 2番目の行の行列の固有値と、行列要素とを算出することである。 固有値算出部 17は、行列 Τ と、行列 Τ との各固有値、及び行列 Τ と、行列 Τ
111 112 111 11 との各行列要素を固有値分解記憶部 16から読み出す。また、行列 Τ を行列 Τ と
2 11 111
、行列 Τ とに分解したときに発生した 2個の値 13、 Θを対角行列記憶部 13から読
112
み出す。そして、固有値算出部 17は、それらの値を用いて、行列 Τ の固有値を算
11
出して固有値分解記憶部 16に蓄積する (ステップ S501)。次に、固有値算出部 17 は、行列 T の固有値を用いて直交行列 Qを算出し (ステップ S502)、その直交行列
11
Qを用いて行列 T の行列要素を算出して固有値分解記憶部 16に蓄積する (ステツ
11
プ S503)。
[0107] 固有値等の 1番目の算出はまだ終わっていないため（ステップ S404)、固有値算出 部 17は、上記説明と同様にして、行列 T 、行列 T の固有値、行列要素等に基づ
121 122
いて、行列 T の固有値、行列要素を算出する（ステップ S402, S403)。このようにし
12
て、固有値等の 1番目の算出が終了すると、固有値算出部 17は、固有値等の 2番目 の算出、すなわち、行列 T 、行列 T の固有値、行列要素等に基づいて行列丁の
11 12 1 固有値、行列要素の算出を行う（ステップ S404, S405, S402, S403)。その後、固 有値算出部 17は、同様にして、行列 T 、行列 Τ の固有値、行列要素等に基づい
21 22
て行列 Τの固有値、行列要素の算出を行う（ステップ S404, S403)。
2
[0108] 固有値算出部 17が固有値等の 2番目の算出を終了すると (ステップ S404)、次の 固有値等の算出は、対称 3重対角行列 Tの固有値を求める最後の処理となるため（ ステップ S405, S402)、固有値算出部 17は、固有値の算出のみを行い、その算出 した固有値を固有値記憶部 18に蓄積する (ステップ S406)。このようにして、固有値 の算出が終了する。
[0109] なお、図 9の行列 T、行列 Tの固有値、行列要素等に基づいて行列 Tの固有値を
1 2
算出する場合には、行列要素 f 、 1、あるいは行列要素 f 、 1を用いるだけでよい。し
1 2 2 1
たがって、行列 T、行列 Tの行列要素を算出する場合に、行列要素 f 、 1、あるいは
1 2 1 2 行列要素 f 、 1
2 1を計算しなくてもよい。このようにすることで、計算量を削減することが でき、処理スピードを向上させることが可能である。
[0110] 次に、固有ベクトル算出部 19について説明する。本実施の形態では、固有ベクトル 算出部 19が（l) qd型ツイスト分解法により固有ベクトルを算出する場合と、（2) LV型 ツイスト分解法により固有ベクトルを算出する場合とについて説明する。
[0111] (l) qd型ツイスト分解法により固有ベクトルを算出する場合
図 10は、 qd型ツイスト分解法により固有ベクトルを算出する場合の固有ベクトル算 出部 19の構成を示す図である。図 10において、固有ベクトル算出部 19は、正定値 化部 21と、コレスキー分解部 22と、ベクトル算出部 23とを備える。
[0112] 正定値ィ匕部 21は、固有値記憶部 18から対称 3重対角行列 Tの固有値を読み出し 、その固有値のいずれかが負の値であるかどうか判断する。そして、その固有値のい ずれかが負の値である場合に、正定値ィ匕部 21は、対角行列記憶部 13から対称 3重 対角行列 Tを読み出し、その対称 3重対角行列 Tの固有ベクトルを変化させることなく 、対称 3重対角行列 Tのすベての固有値が正の値となるように、その対称 3重対角行 列 Tを正定値ィ匕し、正定値ィ匕した後の固有値を固有値記憶部 18に蓄積し、正定値 ィ匕した後の対称 3重対角行列 Tを対角行列記憶部 13に蓄積する。この正定値化の 処理を行うのは、 qd型ツイスト分解法により固有ベクトルを算出する場合には、対称 3 重対角行列 Tの固有値がすべて正の値でなければならないからである。なお、対称 3 重対角行列 Tを正定値ィ匕した行列も、対称 3重対角行列 Τと呼ぶことにする。また、こ の正定値ィ匕の処理が行われた場合には、コレスキー分解や固有ベクトルを算出する 処理は、正定値化後の対称 3重対角行列 Τ、固有値に対して行われることになる。こ の処理の詳細については後述する。ここで、正定値ィ匕の行われた場合であっても、 固有値記憶部 18には、正定値ィ匕の前の固有値が記憶されているものとする。その正 定値化の前の固有値が、固有値分解装置 1が最終的に求める固有値だカゝらである。 また、正定値ィ匕の行われた場合であって、 qd型ツイスト分解法による固有ベクトルの 算出のみを行う場合には、対角行列記憶部 13で記憶されている対称 3重対角行列 T を、正定値ィ匕の後の対称 3重対角行列 Tで上書きしてもよい。一方、 LV型ツイスト分 解法による固有ベクトルの算出も行う場合には、対角行列記憶部 13には、正定値ィ匕 の前の対称 3重対角行列 Tと、正定値化の後の対称 3重対角行列 Tとが記憶されるこ とになる。
[0113] コレスキー分解部 22は、対角行列記憶部 13から、すべての固有値が正の値である 対称 3重対角行列 Tを読み出し、固有値記憶部 18から、対称 3重対角行列 Tの正の 値の固有値を読み出す。そして、コレスキー分解部 22は、対称 3重対角行列 Tの各 要素に関して qd型変換を行うことによって、対称 3重対角行列 Tを上 2重対角行列 B( +1)及び下 2重対角行列 B(_1)にコレスキー分解する。この処理の詳細については後 述する。なお、「対称 3重対角行列 Tを上 2重対角行列 B(+1)及び下 2重対角行列 B(_ "にコレスキー分解する」ことには、対称 3重対角行列 Tの固有値に単位ベクトルを掛 けたものを対称 3重対角行列 Tから引いた行列を、上 2重対角行列 B(+1)及び下 2重 対角行列 B(_1)にコレスキー分解することが含まれてもよい。
[0114] ベクトル算出部 23は、コレスキー分解部 22が算出した上 2重対角行列 B(+1)及び下 2重対角行列 B(_1)の各要素と、対称 3重対角行列 Tの正の値の固有値とを用いて固 有ベクトルを算出して固有ベクトル記憶部 20に蓄積する。この固有ベクトルを各列に 有する行列が、固有値分解の直交行列となる。
[0115] 次に、図 2のフローチャートのステップ S103の処理、すなわち、対称 3重対角行列 Tの固有ベクトルを算出する処理について、より詳細に説明する。まず、図 2のフロー チャートのステップ S 103の処理につ!、て、図 11のフローチャートを用いて説明する [0116] (ステップ S601)正定値ィ匕部 21は、固有値記憶部 18から対称 3重対角行列 Tの固 有値を読み出し、その固有値のいずれかが負の値であるかどうか判断する。そして、 負の値がある場合には、ステップ S602〖こ進み、そうでない場合には、ステップ S603 に進む。
[0117] (ステップ S602)正定値ィ匕部 21は、対角行列記憶部 13から対称 3重対角行列 Tを 読み出し、その対称 3重対角行列 Tの固有ベクトルを変化させることなぐ対称 3重対 角行列 Tのすベての固有値が正の値となるように、その対称 3重対角行列 Tを正定値 化し、正定値ィ匕した後の対称 3重対角行列 Tを対角行列記憶部 13に蓄積し、正定値 化した後の固有値を固有値記憶部 18に蓄積する。
[0118] (ステップ S603)コレスキー分解部 22は、対角行列記憶部 13から、すべての固有 値が正の値である対称 3重対角行列 Tを読み出し、固有値記憶部 18から、対称 3重 対角行列 Tの正の値の固有値を読み出す。なお、正定値化が行われた場合には、 正定値ィ匕の後の対称 3重対角行列 T、及び固有値が読み出されることになる。そして 、コレスキー分解部 22は、対称 3重対角行列 Τの各要素に関して qd型変換を行うこと によって、対称 3重対角行列 Tを上 2重対角行列 B(+1)及び下 2重対角行列 B(_ 1)にコ レスキー分解する。
[0119] (ステップ S604)ベクトル算出部 23は、コレスキー分解部 22が算出した上 2重対角 行列 B(+1)及び下 2重対角行列 B(_1)の各要素と、対称 3重対角行列 Tの正の値の固 有値とを用いて固有ベクトルを算出する。
[0120] (ステップ S605)ベクトル算出部 23は、算出した固有ベクトルを正規ィ匕する。すな わち、ベクトル算出部 23は、算出した固有ベクトルのノルムを算出し、固有ベクトルを 算出したノルムで割ったものを最終的な固有ベクトルとして固有ベクトル記憶部 20に 蓄積する。このようにして、対称 3重対角行列 Tを固有値分解する処理は終了となる。 なお、この後に、行列記憶部 11が記憶している対称行列 Aの固有ベクトルを算出 する処理が行われてもよ 、が、ここでは省略して 、る。
[0121] また、固有ベクトル算出部 19が算出した固有ベクトルの精度を上げるために、図 12 で示されるように、固有ベクトルに関する処理を実行してもよぐあるいは、実行しなく てもよい。すなわち、図示しない逆反復法処理部は、固有ベクトル記憶部 20から固 有ベクトル算出部 19が算出した固有ベクトルを読み出し、その固有ベクトルに対して 逆反復法の処理を実行し、その結果の固有ベクトルを固有ベクトル記憶部 20に蓄積 する (ステップ S701)。このように、逆反復法の処理を実行することによって、固有べ タトルの精度を向上させることができうる。次に、図示しないグラムシュミット処理部は、 高精度が必要力どうか判断する (ステップ S 702)。この判断は、あらかじめ高精度が 必要かどうか設定されている記録媒体等から、その設定を読み出すことによって判断 してもよい。そして、図示しないグラムシュミット処理部は、高精度が必要な場合には、 固有ベクトル記憶部 20から逆反復法の処理の実行された固有ベクトルを読み出し、 その固有ベクトルに対してグラムシュミット法の処理を実行し、その結果の固有べタト ルを固有ベクトル記憶部 20に蓄積する (ステップ S 703)。このように、グラムシュミット 法の処理を実行することによって、固有ベクトルの直交性を高めることができうる。な お、逆反復法や、グラムシュミット法については、すでに公知であり、その詳細な説明 を省略する。
[0122] [対称 3重対角行列 Tの正定値化]
行列 Tを次のように固有値分解できたとする。
[数 38]
T = QDQT
[0123] その行列 Tを次のように a倍して siだけシフトしたとする。ここで、 Iは行列 Tと同じ大き さの単位行列であり、 a, sは実数である。
[数 39]
T→T = aT + sI
[0124] すると、シフト後の固有値分解は、次のようになる。
[数 40] = QDQT
ただし、 D = aD + siである。
[0125] このように、行列 Tを a倍して siだけシフトした場合には、各固有値が a倍されて sだけ シフトされる力 固有ベクトルは変化しない。したがって、この実数倍とシフトとによつ て、すべての固有値を正の値にしてから、ツイスト分解法を用いた固有ベクトルの算 出を行えばよいことになる。
[0126] なお、算出された固有値に負の値が含まれる場合には、例えば、（A)最小の固有 値によりシフトを行ってもよぐ（B)最大の固有値によりシフトを行ってもよい。ここで、 固有値には、負の値と正の値とが含まれているものとする。固有値がすべて正の値で あれば、シフトを行う必要はなぐ固有値がすべて負の値であれば、単に行列 Tを Tとすれば、正定値ィ匕することができる。
[0127] (A)最小の固有値によるシフト
最小の固有値をえ く 0とすると、 a = l、s =— λ + λ X εとする。ここで、 εと
mm min max
しては、マシン'ィプシロンに定数を掛けた値を用いてもよぐ他の適切な値であって もよい。このように、行列 Tをシフトすることによって、行列 Tのすベての固有値を正の 値にすることができる。
[0128] (B)最大の固有値によるシフト
最大の固有値をえ 〉0とすると、 a=— 1、 s = - λ Χ εとする。この場合に
max max mm
は、次のように行列 τが変換されることになる。
[数 41]
Τ→Τ = -Τ + λίΆ3Χ - λαήη χ εΙ = -(Τ - (kmaK - λ^ χ ε )/)
[0129] 上式の右辺の括弧内では、シフトによってすベての固有値が負の値になる力 全 体にマイナスを掛けることによって、正定値ィ匕することができている。ここで、 εは、（Α )の場合と同様にマシン'ィプシロンに定数 (例えば、 2, 3, 4などの実数であり、計算 誤差を少なくするための値を選択することが好適である）を掛けた値を用いてもよぐ 他の適切な値であってもよ ヽ。
[0130] 固有値の最小値の絶対値が、固有値の最大値の絶対値よりも小さい場合には、上 記 (Α)のシフトを行い、固有値の最小値の絶対値力 固有値の最大値の絶対値より も大きい場合には、上記（Β)のシフトを行うようにしてもよい。このようにすることで、シ フトの量を少なくすることができ、丸め誤差を削減することが可能だ力 である。
なお、正定値ィ匕の後に固有ベクトルを算出する処理の詳細については、次の LV型 ツイスト分解法の場合と一緒に説明する。 [0131] (2) LV型ツイスト分解法により固有ベクトルを算出する場合
図 13は、 LV型ツイスト分解法により固有ベクトルを算出する場合の固有ベクトル算 出部 19の構成を示す図である。図 13において、固有ベクトル算出部 19は、コレスキ 一分解部 31と、ベクトル算出部 32とを備える。
[0132] コレスキー分解部 31は、対角行列記憶部 13から、対称 3重対角行列 Tを読み出し 、固有値記憶部 18から、対称 3重対角行列 Tの固有値を読み出す。そして、コレスキ 一分解部 31は、対称 3重対角行列 Tの各要素に関してミウラ変換、 dLVv型変換、逆 ミウラ変換を行うことによって、対称 3重対角行列 Tを上 2重対角行列 B(+1)及び下 2重 対角行列 B(_1)にコレスキー分解する。この処理の詳細については後述する。なお、「 対称 3重対角行列 Tを上 2重対角行列 B(+ 1)及び下 2重対角行列 B(_1)にコレスキー 分解する」ことには、対称 3重対角行列 Tの固有値に単位ベクトルを掛けたものを対 称 3重対角行列 Tから引いた行列を、上 2重対角行列 B(+1)及び下 2重対角行列 B(_1 、にコレスキー分解することが含まれてもよ 、。
[0133] ベクトル算出部 32は、コレスキー分解部 31が算出した上 2重対角行列 B(+1)及び下 2重対角行列 B(_1)の各要素と、対称 3重対角行列 Tの固有値とを用いて固有べタト ルを算出して固有ベクトル記憶部 20に蓄積する。この固有ベクトルを各列に有する 行列が、固有値分解の直交行列となる。
[0134] 次に、図 2のフローチャートのステップ S103の処理、すなわち、対称 3重対角行列 Tの固有ベクトルを算出する処理について、より詳細に説明する。まず、図 2のフロー チャートのステップ S 103の処理につ!、て、図 14のフローチャートを用いて説明する
[0135] (ステップ S801)コレスキー分解部 31は、対角行列記憶部 13から対称 3重対角行 列 Tを読み出し、固有値記憶部 18から対称 3重対角行列 Tの固有値を読み出す。そ して、コレスキー分解部 31は、対称 3重対角行列 Tの各要素に関してミウラ変換、 dL Vv型変換、逆ミウラ変換を行うことによって、対称 3重対角行列 Tを上 2重対角行列 B 及び下 2重対角行列 B(_1)にコレスキー分解する。
[0136] (ステップ S802)ベクトル算出部 32は、上 2重対角行列 B(+1)及び下 2重対角行列 B (_1〉の各要素と、対称 3重対角行列 Tの固有値とを用いて固有ベクトルを算出する。 [0137] (ステップ S803)ベクトル算出部 32は、算出した固有ベクトルを正規ィ匕する。すな わち、ベクトル算出部 32は、算出した固有ベクトルのノルムを算出し、固有ベクトルを 算出したノルムで割ったものを最終的な固有ベクトルとして固有ベクトル記憶部 20に 蓄積する。このようにして、対称 3重対角行列 Tを固有値分解する処理は終了となる。
[0138] なお、この後に、行列記憶部 11が記憶している対称行列 Aの固有ベクトルを算出 する処理が行われてもよいが、ここでは省略している。また、固有ベクトル算出部 19 が算出した固有ベクトルの精度を上げるために、図 12で示されるように、固有べタト ルに関する処理を実行してもよぐあるいは、実行しなくてもよいことは、 qd型ツイスト 分解法の場合と同様である。
[0139] [固有値からの固有ベクトルの算出]
固有値から固有ベクトルを算出する処理については、 qd型ツイスト分解法を用いる 方法と、 LV型ツイスト分解法を用いる方法とについて説明する。図 11のフローチヤ ートで示されるように、 qd型ツイスト分解法を用いて固有ベクトルを算出する場合に、 負の固有値があれば、対称 3重対角行列 Tを正定値ィ匕する処理をあらかじめ実行し ておく（ステップ S601, S602)。 qd型ツイスト分解法を用いる場合には、この後の処 理において、その正定値ィ匕の後の対称 3重対角行列 T、及び固有値を用いるものと する。
[0140] 対角行列記憶部 13では、次のような m X mの対称 3重対角行列 Tが記憶されて 、 るちのとする。
[数 42]
Figure imgf000034_0001
[0141] また、固有値記憶部 18では、固有値え （l≤i≤m)が記憶されているものとする。
のような場合に、対称 3重対角行列 Tの固有ベクトルを求めることは、次の方程式の 固有ベクトル Xを求めることになる。
[数 43] (Γ— λ ) Xj =γΛ = l,2"..,m)
ただし、 行列 Γは mxm行列である。
は、 番目の要素力 で他の要素力^)のべクトルである。
( は、 έ単位行列/の番目の列である )
[0142] 本来であれば、上記式の右辺は 0になるはずである力 行列 Tの固有値を算出する ときに、いくらかの誤差を含むため、固有ベクトル Xが真値であれば、上記のように右 辺に残差項が存在する。
[0143] ここで、以下のようにコレスキー分解できたとする。
[数 44]
T . λ / =
(« = 0，+l).
l ("m）-2
Figure imgf000035_0001
u2 ("m）-l
T ·· ≡Β B
0
2m— 1ノ
, (
u2k-l - ir
¾ )y = i， ^" =1， （" = 0，+l，— 1)，
ξ =ξΓη。)， &^)=^v
[0144] ここで、行列 Tの各要素と、行列 BW)の各要素との対応は、次のようになる。
[数 45]
k = l，z，'，'，
l2k 1k-\ulk = 1,2,… "一 1
ただし、 =0 である。
[0145] なお、行列 Tの各要素と、行列 B(G)の各要素との対応を次のように書くこともできる。
下記の式力も明らかなように、行列 B½)の各要素を、行列 Tの各要素を用いて順次、 算出することができる。 [数 46]
Figure imgf000036_0001
=i ..
)= = 1,2,...，— 1 ただし、 も。 (°) = 0 である。
上記のコレスキー分解した結果を、次式のように書くことができる c
[数 47]
I LD+LT
1 {UD~UL
ただし、 ,ひ は以下の通りである。
Figure imgf000036_0002
そして、次のようになる。
[数 48]
T-%}I = NkDk{N
Figure imgf000036_0003
Dk≡ di D D; ,...,Dk +—い ,Dk+ …, Dm) ここで、行列 Nをツイスト行列と呼ぶ。また
k
[数 49] であるので、（T— λ ΐ)χ = γ eは、
j j k k
[数 50]
N
となる。この簡単な式を解くことによって固有ベクトルを算出することができる。具体的 には、
[数 51]
Figure imgf000037_0001
ただし、 (p)は の p番目の要素である。
のようにわずかな演算で固有ベクトルを算出することができる。なお、ある p 0につい て、 D + = 0あるいは D —ニ でぁったとして 行列 の ， 成分である b p 0 ρθ 2p
(G)と、行列 Bの ， p + 1)成分である b (G)を用いて、あるいは、行列 Tの成分を
-1 2p
用いて、
[数 52] + 2), ρθく k,
（
Figure imgf000037_0002
p0) =
(°) (°)
ρ0>^
^2 0-3ち^2 (°)0-
^2 0+2
(p0 + 2), pO < : ^2 0-4
（p0— 2), ρθ >
po- と固有ベクトルを算出することができる（このように固有ベクトルを算出する処理を例 外処理と呼ぶことにする）。なお、残差項のパラメータ γ及び kの値は、
k
[数 53]
ョ "+ —1)- ) _λ （式 5)
の絶対値が最小となるように決定する。したがって、上述のコレスキー分解を求め、ッ ィスト行列 Nを求めることができれば、固有ベクトルを求めることができることになる。
k
そこで、次にコレスキー分解について説明する。
[0149] 図 15で示すように、 T- λ Iをコレスキー分解することは、 Β(0)に対応する {q (0), e (
j k k
0)}から、 B(±1)に対応する {q (±1), e (±1)}を求めることである。
k k
[0150] (qd型ツイスト分解法）
まず、 qd型ツイスト分解法で用いる qd型変換について説明する（図 15参照)。
[数 54] q d型変換： Γ,^)}!^^^1),^1)}
この変換は、さらに stationary qd with shift (stqds)変換
[数 55]
s t q d s変換 : { ), 。)}4{^+1), +1)}
； ^+ 卜 k=l,2 ,m,
«+1) =«。), k = l,2,.,m-l,
ョ 0, ≡ 0
及び re verse— time progressive qd with shift (rpqds)変換
[数 56] r p q d s変換：
Figure imgf000038_0001
^+^^^+^Π- ,, =
M— ^r^, k=l,2,...,m-l, に分けられる。固有値えが既知であれば反復的な計算が不要なため、計算量を少 なく抑えることができるが、常に数値安定性と精度が高い方法ではない。それは、 stq ds変換、 rpqds変換は、共に減算による桁落ちが発生する可能性があるからである。 例えば、 stqds変換において、 q (0)+e (0)—え〜 e (+"であれば、 q を求め k k-1 j k-1 k
る際に、倍精度計算でもの有効数字がわず力 1桁になることもある。その場合には、 q (0)e (0)/q を計算すると誤差が生じる。つまり、 e が精度よく計算できないこ k k k k
とになる。また、 q を求めるのに e が必要であり、 e を求めるのに q
k+1 k k k が必要であると 、つたように逐次的な計算が要求されるので、 1箇所で発生した桁落 ちによる誤差が波及し、さらなる誤差増大の可能性も秘めている。その結果、理論上 は q (+ 1)≠0であるが、誤差蓄積により q (+ 1) = 0となり、 q (0)e (Vq の計算にお k k k k k
いてオーバーフローが起こるといった数値不安定な状況も想定される。 B(G)の成分 {b (0) , b (0) }が与えられる、すなわち {q (0) , e (0) }が与えられると、 λ及び e
2k- 1 2k k k j k- 1 がー意的に決まるので、この状況は避けることはできない。 rpqds変換も同様の性質 を持っため、実用的なレベルにまで達したとはいいがたい。 LAP ACKにおいて FO RTRANルーチン DSTEGRとして改良版が公開されているものの欠点は完全に解 決されてはいない。
[0151] なお、このように対角行列 Tの各要素を、行列 B(G)の各要素に変換して、その行列 B
の各要素について qd型変換を行う場合であっても、間接的に対角行列 Tの各要素 に関して qd型変換を行うことになるため、「対角行列 Tの各要素に関して qd型変換を 行う」と言うことにする。
また、ここではスタンダードな qd型変換について説明した力 qd型変換は differen tial qd型変換であってもよい。
[0152] (LV型ツイスト分解法）
次に、 LV型ツイスト分解法で用いるミウラ変換、 dLVv型変換、逆ミウラ変換につい て説明する（図 15参照)。
[数 57] ミウラ変換：
s t d L V v変換： °) ^ uk (+1)
r d L V v変換：
逆ミウラ変換 :
Figure imgf000039_0001
[0153] まず、ミウラ変換にっ 、て説明する。この変換は、次のように示される。
[数 58] k = \,2,...,m.
\,2,...,m-\,
Figure imgf000040_0001
≡ 0，
ただし、 Swは任意である。
[0154] 次に、 dLVv型変換について説明する。この変換は、さらに stationary discrete Lotka― Volterra with variable step— size (stdLVv)変換
[数 59] k = 1， …，2»ί— 1，
1 + δ(+1)0' 及び re verse— time discrete Lotka― Volterra with variable step― size ( rdLVv)変換
[数 60]
Figure imgf000040_0002
に分けられる。ただし、
[数 61] = λ (" = +1'— D を満たす範囲内で δ (±1)は任意である。
[0155] 最後に、逆ミウラ変換にっ 、て説明する。この変換は、次のように示される。
[数 62]
(l + ^u^2),(l + ^u^)
¾Jr _ g (+i) , κ ,Α'
Figure imgf000040_0003
≡ 0 [0156] このようにして、 qd型変換と同様に、コレスキー分解を行うことができる。 qd型変換 では見られな 、離散 Lotka— Volterra型変換の大きな特徴は、任意パラメータを持 つことである。すなわち、 λ = 1Ζ δ (°)— 1Z δ (±1)を満たす範囲で δ (η)の値を任意 に設定できる。 δ (η)を変動させると補助変数 u (η)の値も変化するが、桁落ちによる誤
k
差や数値不安定が発生するかどうかは事前に判定できる。この判定は、 if文によって 実装されてもよい。この場合は、 δ (η)を再設定後に再度計算すればよい。また、 u (±1
k
)が求まれば逆ミウラ変換によって q (±1)及び e (±1)が独立に計算されるので、誤差が
k k
伝播しないという性質を持つ。なお、逆ミウラ変換をミウラ変換、ミウラ変換を逆ミウラ 変換と呼んでもよぐ stdLVv変換を stLVv変換と呼んでもよぐ rdLVv変換を rLVv 変換と呼んでもよい。
[0157] なお、このように対角行列 Tの各要素を、行列 B(G)の各要素に変換して、その行列 B の各要素にっ 、て LV型変換を行う場合であっても、間接的に対角行列 Tの各要 素に関して LV型変換を行うことになるため、「対角行列 Tの各要素に関して LV型変 換を行う」と言うことにする。
[0158] ここで、 LV型ツイスト分解法による処理のより詳細な処理の一例について説明する 。図 16〜図 21は、 LV型ツイスト分解法による処理の一例を示すフローチャートであ る。
[0159] 図 16は、コレスキー分解の全体の処理の一例を示す図である。
(ステップ S901)コレスキー分解部 31は、ミウラ変換を行う。この処理の詳細につい ては後述する。
[0160] (ステップ S902)コレスキー分解部 31は、 1Ζ δ (±1)を lZ δ (0) - λとする。
(ステップ S903)コレスキー分解部 31は、後述する Procedure 1の処理を実行する
[0161] (ステップ S904)コレスキー分解部 31は、後述する Procedure2の処理を実行する
(ステップ S905)コレスキー分解部 31は、 q e がすでに算出されているか
k k
どうか判断する。そして、すでに算出されている場合には、コレスキー分解の一連の 処理は終了となり、一方、算出されていない場合には、ステップ S901に戻る。 [0162] 図 17は、図 16のフローチャートにおけるステップ S903の処理の詳細を示すフロー チャートである。
(ステップ S1001)コレスキー分解部 31は、 q e がすでに算出されている
k k
力どうか判断する。そして、すでに算出されている場合には、終了となり、算出されて ヽな ヽ場合に【ま、ステップ S1002に進む。
(ステップ S1002)コレスキー分解部 31は、 stdLVv変換等の処理を行う。この処理 の詳細については後述する。
[0163] 図 18は、図 16のフローチャートにおけるステップ S904の処理の詳細を示すフロー チャートである。
(ステップ S1101)コレスキー分解部 31は、 q e がすでに算出されている
k k
力どうか判断する。そして、すでに算出されている場合には、終了となり、算出されて ヽな ヽ場合に【ま、ステップ S1102に進む。
(ステップ S1102)コレスキー分解部 31は、 rdLVv変換等の処理を行う。この処理 の詳細については後述する。
[0164] 図 19は、図 16のフローチャートのステップ S901の処理の詳細を示すフローチヤ一 トである。
(ステップ S1201)コレスキー分解部 31は、 1Z δ (G)の値を決定する。この値は、前 述のように任意に決定することができるため、行列を正定値ィ匕する値に決定する。例 えば、コレスキー分解部 31は、 1Z δ (G)の値を次式のように決定する。
[数 63]
- ^ = | min| + e (°) χ |λ匪
[0165] なお、なお、 ε (ωは、マシン'ィプシロンに定数 (例えば、 2, 3, 4などの実数であり、 計算誤差を少なくするための値を選択することが好適である)を掛けた値を用いても よぐ他の適切な値であってもよい。その後、例えば、ステップ S 1203等で桁落ちの 発生する可能性があると判断された場合に、コレスキー分解部 31は、 1Z δ Wの値を 決める ε (G)の値を、例えば、マシン'ィプシロンに掛ける定数を 1ずつ大きくしながら 設定していってもよい。なお、 1Z δ の値を決定する方法は、これに限定されないこ とは言うまでもない。このように、 LV型ツイスト分解法では、この 1Z δ (ωの値を自由 に設定することができるため、 qd型ツイスト分解法の場合のように、行列 Tの正定値化 を行わなくてもよいことになる。なお、 qd型ツイスト分解法の場合と同様に、行列丁の 正定値ィ匕の処理を行ってから、 LV型ツイスト分解法を実行してもよ ヽことは言うまで もない。
(ステップ S1202)コレスキー分解部 31は、 β ΐを α ( )— lZ δ (ωに設定する。
[0166] (ステップ S1203)コレスキー分解部 31は、 β 1の絶対値が εよりも大きいかどうか 判断する。そして、大きい場合には、ステップ S1204に進み、そうでない場合には、 ステップ S1201に戻って、 1/ δ (G)の値を決定する処理を再度実行する。なお、この εも任意に定めることができる。例えば、 εとしてマシン'ィプシロンを用いてもよい。 εの値を大きくすると精度が向上し、 εの値を小さくすると精度が低下する。なお、こ のステップ S1203で行っている処理は、桁落ちが発生する可能性について判断する 処理である。 β 1の絶対値が εよりも大きくない場合には、桁落ちの発生する可能性 があると半 IJ断されることになる。
[0167] (ステップ S1204)コレスキー分解部 31は、 β 2を q (0)/ (1 + δ (0)u (0))に設定す
1 0
る。なお、前述のように u ((>) =0であるため、 |8 2は q (())に設定されたことになる。
0 1
[0168] (ステップ S1205)コレスキー分解部 31は、 u δ (°)u (°))を j8 1に設定する。
1 0
ここで、 13 1は、ステップ S 1202において q (0) - 1/ δ ( )に設定されている力 前述 のように u (。) = 0であるため、 β 1は、 q (。ソ（1 + δ (。)u (。）— 1Ζ δ (。)と等しぐこれ
0 1 0
にミウラ変換を行うと u (0) =u (0) I となる力もである。なお、 u (°) =0から（1 +
1 2k- 1 k=l 0
δ (% (°)) = 1である。
0
[0169] (ステップ S1206)コレスキー分解部 31は、カウンタ kを「1」に設定する。
(ステップ S1207)コレスキー分解部 31は、 a ^e (G)Z j8 1に設定する。前述のよ k
うに、 β 1は u (°)となるから、 e (°)Ζ ι8 1にミウラ変換を行うと、 a 1は δ (0)u (°)と等
2k- 1 k 2k し!ヽこと〖こなる。
[0170] (ステップ S1208)コレスキー分解部 31は、 α 2を α 1 + 1に設定する。ステップ S 12 07の説明からわ力るように、 ひ 2は 1 + S ( )u ( )と等しいことになる。
2k
[0171] (ステップ S1209)コレスキー分解部 31は、 α 2の絶対値が εよりも大きいかどうか 判断する。そして、大きい場合には、ステップ S1210に進み、そうでない場合には、 ステップ S1201に戻って、 1/ δ to)の値を決定する処理を再度実行する。なお、この ステップ S1209で行っている処理は、ステップ S1203の処理と同様に、桁落ちが発 生する可能性について判断する処理である。 α 2の絶対値が εよりも大きくない場合 には、桁落ちの発生する可能性があると判断されることになる。
[0172] (ステップ S1210)コレスキー分解部 31は、 《IX |82を算出して u (°)(1+ S (°)u
2k 2k
(0))に設定する。前述のように、 a 1は δ (0)u ( )と等しく、 β 2は q (0)/(1+ δ ( )u
-1 2k k 2
(。)）であるため、 αΙΧ β 2にミウラ変換を実行すると、 u (0) (1+ δ (o)u (。)）に k-2 2k 2k- 1 等しくなるからである。
[0173] (ステップ S1211)コレスキー分解部 31は、 |82を q (G)Z α 2に設定する。前述の k+1
ように、 《2は 1+ S (°)u (°)と等しいため、 |82は、 q (0)/(l+ 6 (0)u (°))となり、
2k k+ 1 2k
ι82における kの値を 1だけインクリメントしたことになる。
[0174] (ステップ S1212)コレスキー分解部 31は、 βΐを β2— 1 δ Wに設定する。前述 のように、 ι82は q (°)Ζ(1+ δ ( (°))と等しいため、 |81は、 q (°
k+1 2k k+1 ソ（1+ δ (
(0))-1/ δ (0)となり、ミウラ変換を実行すると、 β 1は、 u (0)となる。したがって、
2k 2k+l
ι81における kの値を 1だけインクリメントしたことになる。
[0175] (ステップ S1213)コレスキー分解部 31は、 131の絶対値が εよりも大きいかどうか 判断する。そして、大きい場合には、ステップ S1214に進み、そうでない場合には、 ステップ S1201に戻って、 1/ δ (G)の値を決定する処理を再度実行する。なお、この ステップ S1213で行っている処理は、ステップ S1203の処理と同様に、桁落ちが発 生する可能性について判断する処理である。 β 1の絶対値が εよりも大きくない場合 には、桁落ちの発生する可能性があると判断されることになる。
[0176] (ステップ S1214)コレスキー分解部 31は、 《2Χ |81を算出して u δ (°)
2k+l
u ( ))に設定する。前述のように、 《2は 1+ δ (% ( )と等しぐ j81は、 u ( )に等
2k 2k 2k+ 1 しいからである。
(ステップ S 1215)コレスキー分解部 31は、カウンタ kを 1だけインクリメントする。
[0177] (ステップ S1216)コレスキー分解部 31は、カウンタ kが mに等しいかどうか判断す る。そして、 mに等しい場合には、一連の処理は終了となり、そうでない場合には、ス テツプ SI 207に戻る。
[0178] このようにして、 u (0)(l+ 6 (0)u (0))と、 u (0)(l+ 6 (0)u (0))とが算出される
2k+l 2k 2k 2k- 1
ことになる。これらの値は、コレスキー分解部 31が有する図示しないメモリ等において 一時的に記憶されてもよい。
[0179] 図 20は、図 17のフローチャートのステップ S1002の処理の詳細を示すフローチヤ ートである。
(ステップ S1301)コレスキー分解部 31は、 ひ2を1+ δ (+1)u に設定する。なお
0
、前述のように u (+1)=0であるため、 α 2は 1に設定されたことになる。
0
[0180] (ステップ S1302)コレスキー分解部 31は、 β 1を u δ (°)u δ (+1
1 0
u に設定する。ここで、 u (0) (1+ δ (0)u (0))の値としては、ステップ S1205で算
0 1 0
出したものを用いる。なお、前述のように u (+1)=0である。また、この β 1の式に stdL
0
Vv変換を実行すると、 β 1は、 u に設定されたことになる。
[0181] (ステップ S1303)コレスキー分解部 31は、カウンタ kを「1」に設定する。
(ステップ S1304)コレスキー分解部 31は、 α 1を |81 + 1Z δ に設定する。
[0182] (ステップ S1305)コレスキー分解部 31は、 α 1の絶対値が εよりも大きいかどうか 判断する。そして、大きい場合には、ステップ S1306に進み、そうでない場合には、 図 16のフローチャートの Procedure2(ステップ S904)〖こ進む。なお、このステップ S 1305で行っている処理は、ステップ S1203の処理と同様に、桁落ちが発生する可 能性について判断する処理である。 a 1の絶対値が εよりも大きくない場合には、桁 落ちの発生する可能性があると判断されることになる。
[0183] (ステップ S1306)コレスキー分解部 31は、 j82を u δ ( )u ( ))Ζα1に
2k 2k- 1
設定する。ここで、 u (0)(1+ δ (0)u (0))の値としては、ステップ S1210で算出した
2k 2k- 1
ものを用いる。また、この 132の式に stdLVv変換を実行すると、 β 2は、 δ (+1)u
2k に設定されたことになる。
[0184] (ステップ S1307)コレスキー分解部 31は、 α 1 X α2を算出する。前述のように、 a
1は ι81 + 1/δ (+1)=ιι (+1) + 1/δ (+1)に等しぐ ひ 2は l+ S (+1)u に等
2k- 1 2k- 2 しいため、 α IXひ2に逆ミウラ変換を実行すると q に等しくなる。したがって、コレ k
スキー分解部 31は、 《IXひ2を算出して q に設定する。 [0185] (ステップ S1308)コレスキー分解部 31は、 0；2を1 + β 2に設定する。前述のように 、 j8 2は S (+1)u と等しいため、 ひ2は 1 + δ (+1)u (+1)となる。したがって、 0；2に
2k 2k
おける kの値を 1だけインクリメントしたことになる。
[0186] (ステップ S1309)コレスキー分解部 31は、 α 2の絶対値力 S εよりも大きい力どう力 判断する。そして、大きい場合には、ステップ S1310に進み、そうでない場合には、 図 16のフローチャートの Procedure2 (ステップ S904)〖こ進む。なお、このステップ S 1309で行っている処理は、ステップ S1203の処理と同様に、桁落ちが発生する可 能性について判断する処理である。 α 2の絶対値が εよりも大きくない場合には、桁 落ちの発生する可能性があると判断されることになる。
[0187] (ステップ S1310)コレスキー分解部 31は、 j8 1 X j8 2を算出する。前述のように、 β 1は u に等しく、 2は S (+1)u と等しいため、 |8 1 Χ |8 2に逆ミウラ変換を
2k- 1 2k
実行すると、 e に等しくなる。したがって、コレスキー分解部 31は、 β Ι Χ β 2を算
k
出して e に設定する。
k
[0188] (ステップ S1311)コレスキー分解部 31は、 |8 1を u (0) (1 + δ (0)u (0)) / α 2に
2k+l 2k
設定する。ここで、 u (0) (1 + δ (0)u の値としては、ステップ S1214で算出した
2k+l 2k
ものを用いる。また、この β 1の式に stdLVv変換を実行すると、 |8 1は、 u に
2k+l 設定されたことになる。したがって、 13 1における kの値を 1だけインクリメントしたことに なる。
[0189] (ステップ S1312)コレスキー分解部 31は、カウンタ kを 1だけインクリメントする。
(ステップ S1313)コレスキー分解部 31は、カウンタ kが mに等しいかどうか判断す る。そして、 mに等しい場合には、ステップ S1314に進み、そうでない場合には、ステ ップ S 1304に戻る。
[0190] (ステップ S1314)コレスキー分解部 31は、 α 2 X ( 1 + ΐΖ δ を算出する。こ れ ίま、ステップ S 1304で a 1を更新した後に、ステップ S1307で α Ι Χ α 2を計算す ることと等しい。したがって、コレスキー分解部 31は、 《2 Χ ( l + ΐΖ δ を算出 して q に設定する。このようにして、ミウラ変換された結果に stdLVv変換と、逆ミゥ ラ変換とを実行して、 q (+1)、 e を算出する処理が終了する。これらの値は、コレス
k k
キー分解部 31が有する図示しな 、メモリ等にぉ 、て一時的に記憶されてもょ 、。 [0191] 図 21は、図 18のフローチャートのステップ SI 102の処理の詳細を示すフローチヤ ートである。
(ステップ S1401)コレスキー分解部 31は、 131を u (0)(1+ δ (。)u (。)）
2m- 1 2m- 2 Z(1
+ δ (_1)u に設定する。ここで、 u (0)(l+ 6 (0)u (0))の値としては、ステ
2m 2m- 1 2m- 2
ップ S1214で算出したものを用いる。なお、前述のように u (_1)=0である。また、こ
2m
の 131の式に rdLVv変換を実行すると、 |81は、 u に設定されたことになる。
2m- 1
[0192] (ステップ S1402)コレスキー分解部 31は、カウンタ kを「m— 1」に設定する。
(ステップ S1403)コレスキー分解部 31は、 α 1を |81 + 1Z δ に設定する。
[0193] (ステップ S1404)コレスキー分解部 31は、 α 1の絶対値が εよりも大きいかどうか 判断する。そして、大きい場合には、ステップ S1405に進み、そうでない場合には、 図 16のフローチャートのミウラ変換 (ステップ S901)に戻る。なお、このステップ S 140 4で行っている処理は、ステップ S1203の処理と同様に、桁落ちが発生する可能性 について判断する処理である。 a 1の絶対値が εよりも大きくない場合には、桁落ち の発生する可能性があると判断されることになる。
[0194] (ステップ S1405)コレスキー分解部 31は、 j82を u δ (°)u (°))Ζα1に
2k 2k- 1
設定する。ここで、 u (0)(1+ δ (0)u (0))の値としては、ステップ S1210で算出した
2k 2k- 1
ものを用いる。また、この β 2の式に rdLVv変換を実行すると、 β 2は、 δ (_1)u
2k に設定されたことになる。
[0195] (ステップ S1406)コレスキー分解部 31は、 0；2を1+ β 2に設定する。前述のように 、 j82は S (_1)u と等しいため、 ひ2は 1+ δ (_1)u となる。
2k 2k
[0196] (ステップ S1407)コレスキー分解部 31は、 α 2の絶対値力 S εよりも大きい力どう力 判断する。そして、大きい場合には、ステップ S1408に進み、そうでない場合には、 図 16のフローチャートのミウラ変換 (ステップ S901)に戻る。なお、このステップ S 140 7で行っている処理は、ステップ S1203の処理と同様に、桁落ちが発生する可能性 について判断する処理である。 《2の絶対値が εよりも大きくない場合には、桁落ち の発生する可能性があると判断されることになる。
[0197] (ステップ S1408)コレスキー分解部 31は、 α 1 X α2を算出する。前述のように、 a
1は ι81 + 1/δ (_1)=ιι (_1) + 1/δ (_1)に等しぐ ひ 2は 1+ S (_1)u に等し いため、ひ 1 X ひ 2に逆ミウラ変換を実行すると q に等しくなる。したがって、コ k+1
レスキー分解部 31は、 α 1 Χ α 2を算出して q に設定する。
k+1
[0198] (ステップ S1409)コレスキー分解部 31は、 β 1を u (0) (1 + δ (0)u (0)) / α 2
2k- 1 2k- 2 に設定する。ここで、 u (0) (1 + δ (0)u (°))の値としては、ステップ S1205, S12
2k- 1 2k- 2
14で算出したものを用いる。また、この j8 1の式に stdLVv変換を実行すると、 1は 、 u に設定されたことになる。したがって、 13 1における kの値を 1だけデクリメン
2k- 1
卜したこと〖こなる。
[0199] (ステップ S1410)コレスキー分解部 31は、 j8 1 X j8 2を算出する。前述のように、 β
1は u 卜1)に等しく、 13 2は δ (_1)u と等しいため、 I X 2に逆ミウラ変換を
2k- 1 2k
実行すると、 e に等しくなる。したがって、コレスキー分解部 31は、 β Ι Χ β 2を算 k
出して e に設定する。
k
[0200] (ステップ S1411)コレスキー分解部 31は、カウンタ kを 1だけデクリメントする。
(ステップ S1412)コレスキー分解部 31は、カウンタ kが 0に等しいかどうか判断する 。そして、 0に等しい場合には、ステップ S 1413に進み、そうでない場合には、ステツ プ S 1403に戻る。
[0201] (ステップ S1413)コレスキー分解部 31は、 1 + ΐΖ δ を算出する。これは、ス テツプ S1403で a 1を更新した後に、ステップ S1408で α Ι Χ α 2を計算することと 等しい。この場合、 α 2は 1だ力もである。したがって、コレスキー分解部 31は、 |8 1 + 1/ δ を算出して q に設定する。このようにして、ミウラ変換された結果に rdL Vv変換と、逆ミウラ変換とを実行して、 q (_1)、 e を算出する処理が終了する。こ k k
れらの値は、コレスキー分解部 31が有する図示しないメモリ等において一時的に記 憶されてもよい。
[0202] なお、図 16のフローチャートの処理によって算出することができるのは 1個の固有 値に対応する固有ベクトルであるため、全ての固有値に対応する固有ベクトルを算出 する場合には、コレスキー分解部 31は、図 16のフローチャートの処理を固有値の数 だけ繰り返すことになる。
[0203] メモリ消費を抑えるために、補助変数 u (n)のための配列は必ずしも用意する必要は k
ない。一方、 1 + δ (G)u(G)のためのメモリ領域を確保し、ミウラ変換、 dLVv型変換、逆 ミウラ変換のステップにまたがつてこの値を利用することで、メモリ消費を抑え、計算 量を低減することができる。これにより、誤差も低減される。
[0204] ここで、誤差について説明する。人間が理想的な状況、すなわち、無限桁の計算を いくらでもできるとすると、 qd型ツイスト分解法であっても、 LV型ツイスト分解法であつ ても問題ない。しかしながら、コンピュータで計算を行う場合は注意が必要である。有 限桁の計算し力、行えないコンピュータ上では、数学的に正しい計算法を使ったとして も必ずしも正しい結果が得られる訳ではない。そればかりか、いつまでたっても計算 が終了しないといった思わぬ数値的な問題が発生する場合もある。
[0205] コンピュータ計算による誤差には、丸め誤差及び桁落ちによる誤差などが知られて いる。丸め誤差単独では、高々有効桁の最後の桁が真値と比べて異なる程度で大き な誤差にはならない。また、指数部が異なる 2つの実数の加算、乗算、除算の少なく とも 1つの演算を行えばやはり丸め誤差が生じるが、それ以上の誤差は発生しな 、。 さらに、このような丸め誤差を発生するような操作が繰り返されても、丸めモードが nea r (四捨五入)ならば、一方的に切り上げられたり、あるいは切り捨てられたりして誤差 が極端に蓄積することは少ない。よって、多くの数値計算法は加算、乗算、除算の少 なくとも 1つの演算によって発生する丸め誤差を特別注意することは少なぐ適切な 固有値計算でも結果的に丸め誤差は一様に増大しない。
[0206] 問題となるのは、同符号の実数の減算及び異符号の実数の加算により生じる、桁 落ちである。桁落ちによる誤差で値が 0となった後、その値による除算を行うと、 0が分 母にくるような不定形となり計算不可能となる。こうなるといつまでたつても計算が終了 しない。減算→除算と計算する部分が qd型ツイスト分解法、 LV型ツイスト分解法の 両方に存在するので、桁落ち誤差には十分に注意する必要がある。
[0207] LV型ツイスト分解法では、上述の桁落ちによる誤差を含んでいるかどうかは減算に よって得られた値が小さいかどうかで判断できる。 qd型ツイスト分解法の場合、桁落 ち誤差を含むことが分力つたとしても、それを回避することはできない。なぜならば、 初期値として {q (0), e (° が与えられると、 λは一意的に決定され、 {q e
k k j k k も一意的に導出されるためである。すなわち、任意パラメータを持たない自由度のな V、計算法であるためである。 [0208] それに対して、 LV型ツイスト分解法は、自由に設定できるパラメータ δ (0)を持った め、補助変数 u (η)の値を様々に変化させることができる（図 22参照)。すなわち、様 k
々な経路で {q (±1), e (±1)}を計算することができる。よって、桁落ちが発生する場合 k k
も回避できる。図 19〜図 21の条件判定によって桁落ちの影響をチェックし、減算に よって得られた値の絶対値力 、さな数 εより大きいという条件が満たされなければ、 ノ ラメータ δ (G)の設定に戻るというものである。この処理は、条件が満たされるまで繰 り返される。なお、精度よりも高速性を重視する場合は、数回条件が満たさなければ（ q (+1)=0あるいは q (_1)=0ならば）、例外処理を行ってもよい。
k k
[0209] なお、前述の行列 Tの各要素と、行列 B(G)の各要素との対応を示す式を用いると、 q d型変換を次のように書き換えることもできる。
[数 64]
q d型変換： - ^
[0210] この変換は、さらに stqds変換
[数 65]
s t q d s変換：
Figure imgf000050_0001
q( l +e^=tu_1- J, k = U
qi+1)ei+1)
Figure imgf000050_0002
= l,2"..,m_l,
E0
及び rpqds変換
[数 66]
r p q d s変換： {t2t i， }
+e x ~^j-> k = l,2,...,m,
に分けられる。このように、 {t , t }から {q (±1), e (±1)}を求める場合には、 {t
2k- 1 2k k k 2k- 1
, t }から {q (。)， e (。)}を求め、さらに q (。)， e (。)}から {q (土 e (±1)}を求める場合に
2k k k k k k k
比べて、計算量が少なくなることとなり、計算速度を向上させることも可能である。
[0211] なお、この場合には、対角行列 Tの各要素について直接 qd型変換が行われること になる。したがって、この変換が、「対角行列 Tの各要素に関して qd型変換を行う」場 合に含まれることは言うまでもな 、。
[0212] また、 qd型変換と同様に、前述の行列 Tの各要素と、行列 Β(ωの各要素との対応を 示す式を用いると、 LV型変換におけるミウラ変換を次のように書き換えることもできる
[数 67] ミウラ変換： {υ^ ί"^,^)}
s t d L V v変換： ） u[+1)
r d L V v変換： 。）
逆ミウラ変換：
Figure imgf000051_0001
[0213] この場合のミウラ変換は、次のように示される。
[数 68] t =^(1+5 (。) — 2 XI + δ ( ) + δ (。) 3 )— 2 k = l,2,.,m
(t2k)2
Figure imgf000051_0002
k = 1,2,… -1
≡ ο,
ただし、 Swは任意である。
このように、 {t , t in (0), u (。)}を求める場合には、 {t , 1; }から (
2k- 1 2k 2k- 1 2k 2k- 1 2k k o), e (。)}を求め、さらに q )， e (。)}から {u (0), u (。)}を求める場合に比べて、計 k k k 2k- 1 2k
算量が少なくなることとなり、計算速度を向上させることも可能である。
[0214] なお、この場合には、対角行列 Tの各要素について直接 LV型変換が行われること になる。したがって、この変換が、「対角行列 Tの各要素に関して LV型変換を行う」場 合に含まれることは言うまでもな 、。
[0215] qd型ツイスト分解法あるいは LV型ツイスト分解法によってコレスキー分解をすること ができると、上述のようにツイストされた行列 Nを算出することができ、その行列 Nの k k
N Tx =eを解くことによって、 χを算出することができる。ここで、
k j k j
[数 69]
INI
と Xを置き換えることにより、この Xを正規ィ匕する。このようにして、固有ベクトル Xを求 めること力できる。この固有べク卜ノレ Xを用！/、て、 Q = (x , X , ... , X )とすることにより j T 1 2 m
、直交行列 Qを算出することができる。
T
[0216] したがって、次式のように対称 3重対角行列 Tについて固有値分解が行われたこと になる。
[数 70]
T = QTDTQ
[0217] また、対称 3重対角行列 Tの固有値分解を行うことができれば、次のように、対称行 列 Aの固有値分解も行われる。
[数 71]
Α = ΡΤΡΊ
= P{QTDTQT T )PT
= (PQT )DT (PQT
Figure imgf000052_0001
[0218] 固有ベクトル算出部 19は、上述のようにして、固有値記憶部 18が記憶している対 称 3重対角行列 Tの固有値と、対角行列記憶部 13が記憶している対称 3重対角行列 Tとを用いて、対称 3重対角行列 Tの固有ベクトルを算出する。まず、コレスキー分解 部 31は、対角行列記憶部 13から対称 3重対角行列 Tを読み出し、固有値記憶部 18 から対称 3重対角行列 Tの固有値を読み出す。そして、コレスキー分解部 31は、図 2 3のフローチャートで示されるように各変換を行い、コレスキー分解の処理を行う。ここ では、 LV型ツイスト分解法を用いる場合にっ 、て説明する。
[0219] コレスキー分解部 31は、まず、対称 3重対角行列 Tの各要素の値から、 q (0), e (0) k k を求める。そして、コレスキー分解部 31は、ミウラ変換を実行することにより、 u (0), u (
1 2 ω等を順次求めていく（ステップ S1501)。
[0220] 次に、コレスキー分解部 31は、ミウラ変換で得られた u (ωを用いて、 stdLVv変換を k
実行することにより、 u 等を順次求めていく（ステップ S1502)。また、コレスキー 分解部 31は、ミウラ変換で得られた u ((>)を用いて、 rdLVv変換を実行することにより k
、 u 卜 υ等を順次求めていく（ステップ S1503)。
2m- 1
[0221] 最後に、コレスキー分解部 31は、 stdLVv変換で得られた u 及び rdLVv変換で 得られた u を用いて、逆ミウラ変換を実行することにより、 q (±1)、 e (±1)等を順次 k 1 1
求めていく（ステップ SI 504)。 q (±1)及び e (±1)が求まると、すなわち、上 2重対角行 k k
列 B(+1)と、下 2重対角行列 B(_1)が求まると、コレスキー分解部 31は、それらをべタト ル算出部 32に渡す。なお、コレスキー分解部 31は、各固有値について、それぞれコ レスキー分解を行うものとする。
[0222] なお、ここでは、説明の便宜上、図 23のフローチャートを用いてコレスキー分解の 処理を説明したが、図 16〜図 21のフローチャートで示されるように処理を行ってもよ いことは言うまでもない。
[0223] ベクトル算出部 32は、固有値記憶部 18から固有値を読み出し、式 5を用いて kの値 を決定する。ベクトル算出部 32は、その kの値を用いて、 q (±1)及び e (±1)からツイスト k k
行列 Nを求め、 N Tx =eを解くことによって、 Xを算出する（ステップ S802)。
k k j k j
[0224] 次に、ベクトル算出部 32は、算出した Xを正規ィ匕して、固有ベクトル Xを求めて固有 ベクトル記憶部 20に蓄積する (ステップ S803)。なお、ベクトル算出部 32は、各固有 値について、 Xを算出する処理と、正規ィ匕する処理とを行うことによって、全ての固有 ベクトルを算出することができる。このようにして、固有ベクトルの算出が終了する。 この後、行列記憶部 11が記憶して 、る対称行列 Aの固有ベクトルを算出してもよ!/ヽ 。その算出の方法は、前述の通りである。
[0225] また、固有値算出部 17が算出した固有値、固有ベクトル算出部 19が算出した固有 ベクトルを出力する図示しない出力部を固有値分解装置 1が備えてもよい。ここで、 図示しない出力部による出力は、例えば、固有値等の表示デバイス (例えば、 CRT や液晶ディスプレイなど)への表示でもよぐ固有値等の所定の機器への通信回線を 介した送信でもよぐ固有値等のプリンタによる印刷でもよぐ固有値等の記録媒体へ の蓄積でもよい。なお、その出力部は、出力を行うデバイス (例えば、表示デバイスや プリンタなど）を含んでもよぐあるいは含まなくてもよい。また、その出力部は、ハード ウェアによって実現されてもよぐあるいは、それらのデバイスを駆動するドライバ等の ソフトウェアによって実現されてもよい。
[0226] なお、ここでは、コレスキー分解部 31が LV型ツイスト分解法によってコレスキー分 解を行う場合について説明したが、コレスキー分解部 31は、 qd型ツイスト分解法によ つてコレスキー分解を行ってもよい。例えば、コレスキー分解部 31は、算出された固 有値の分布を見て、分布が密でない場合には、 qd型ツイスト分解法を用いるようにし てもよい。
[0227] また、上記説明では、行列 Bto)が上 2重対角行列である場合について説明したが、 行列 B(C>)は下 2重対角行列であっても固有値分解を同様に実行することができうる。 ただし、行列 Tの各要素と、行列 B(G)の各要素との対応を示す式が上述のものと異な ることになる。
[0228] また、上記説明では、固有値算出部 17が全ての固有値を算出する場合について 説明したが、一部の固有値のみを算出するようにしてもよい。例えば、図 9において、 固有値算出部 17は、行列 Τ , Tまでは、全ての固有値を算出する必要があるが、行
1 2
列 Τ , Tから対称 3重対角行列 Tの固有値を算出する際に、必要な範囲で固有値を
1 2
算出してもよい。したがって、固有値算出部 17は、対称 3重対角行列 Tの少なくとも 1 個の固有値を算出するものであってもよい。その場合には、不必要な固有値を算出 することに伴う余分な処理を実行しなくてよいため、処理負荷を軽減することができう る。このように、一部の固有値のみを算出する場合には、例えば、対称 3重対角行列 Tの固有値を算出する処理に費やす計算時間を約（1 +kZm) Z2倍にすることがで きうる。ここで、 mは行列サイズであり、 kは求める固有値の個数である。
[0229] また、上記説明では、固有ベクトル算出部 19が算出された固有値に対応する全て の固有ベクトルを算出する場合について説明したが、上記説明から明らかなように、 固有ベクトルを算出する処理は、固有値ごとに行うことができる。したがって、固有べ タトル算出部 19は、対称 3重対角行列 Tの少なくとも 1個の固有ベクトルを算出するも のであってもよい。このように、固有ベクトル算出部 19は、必要な範囲で固有ベクトル を算出することができ、不必要な固有ベクトルを算出することに伴う余分な処理を実 行しなくてよいため、処理負荷を軽減することができうる。
[0230] 次に、本実施の形態による固有値分解装置 1における並列処理について説明する 本実施の形態による固有値分解装置 1では、固有値の計算及び固有ベクトルの計 算において、並列的に処理を実行することができる。例えば、図 24〜図 26で示され るように、固有値分解装置 1において、固有値分解部 15、固有値算出部 17、コレス キー分解部 22、ベクトル算出部 23,コレスキー分解部 31、ベクトル算出部 32におい て、それぞれ並列処理を行ってもよい。
[0231] 固有値分解部 15は、複数の固有値分解手段 15a、 15bを備え、その複数の固有 値分解手段 15a、 15bが、分割後の対称 3重対角行列の固有値と、行列要素とを算 出する処理を並列実行してもよい。例えば、図 9において、固有値分解手段 15aが行 列 Τ , Τ , Τ , Τ の固有値分解を実行し、固有値分解手段 15bが行列 Τ ,
111 112 121 121 211 τ 。
212 , τ の固有値分解を実行してもよい
221 , τ 222
[0232] 固有値算出部 17は、複数の固有値算出手段 17a、 17bを備え、その複数の固有 値算出手段 17a、 17bが、分割元の対称 3重対角行列の固有値と、行列要素とを算 出する処理を並列実行してもよい。例えば、図 9において、固有値算出手段 17aが行 列 Τ , Τ , Τ , Tの固有値等を算出する処理を実行し、固有値算出手段 17bが行
11 12 1
列 τ , τ , Tの固有値等を算出する処理を実行してもよい。
21 22 2
[0233] コレスキー分解部 22は、複数のコレスキー分解手段 22a、 22bを備え、その複数の コレスキー分解手段 22a、 22bが、対称 3重対角行列 Tに関してコレスキー分解する 処理を並列実行してもよい。例えば、コレスキー分解手段 22aが半分の固有値につ いてコレスキー分解する処理を実行し、コレスキー分解手段 22bが残りの半分の固有 値にっ 、てコレスキー分解する処理を実行してもよ 、。
[0234] ベクトル算出部 23は、複数のベクトル算出手段 23a、 23bを備え、その複数のベタト ル算出手段 23a、 23bが、固有ベクトルを算出する処理を並列実行してもよい。例え ば、ベクトル算出手段 23aがコレスキー分解手段 22aによってコレスキー分解された 固有値に関して固有ベクトルを算出する処理を実行し、ベクトル算出手段 23bがコレ スキー分解手段 22bによってコレスキー分解された固有値に関して固有ベクトルを算 出する処理を実行してもよ ヽ。
[0235] コレスキー分解部 31が複数のコレスキー分解手段 31a、 31bを備えて、対称 3重対 角行列 Tに関してコレスキー分解する処理を並列実行してもよいことは、コレスキー分 解部 22と同様である。
[0236] ベクトル算出部 32が複数のベクトル算出手段 32a、 32bを備えて、固有ベクトルを 算出する処理を並列実行してもよいことは、ベクトル算出部 23と同様である。
[0237] なお、固有値算出部 17の複数の固有値算出手段 17a、 17bは、図 6で示される分 割された 2個の行列 Τ , Tの固有値と行列要素とから、分割元の行列 Tの固有値と
1 2 0
行列要素とを算出する処理を、並列実行してもよい。以下、その処理について説明 する。
[0238] まず、式 1から、各固有値を求める処理については、並列実行可能なことがわかる。
例えば、固有値算出手段 17aが半分の固有値を算出し、固有値算出手段 17bが残り の半分の固有値を算出してもよい。なお、固有値算出手段 17a、 17bのそれぞれは、 行列 D と、ベクトル zとの各成分の値を有するものとする。
12
[0239] また、式 2から、各固有ベクトル qを求める処理についても、並列実行可能なことが わかる。例えば、固有値算出手段 17aは、算出した固有値に対応する半分の固有べ タトル qを算出し、固有値算出手段 17bは、算出した固有値に対応する残りの半分の 固有ベクトル qを算出してもよい。
[0240] また、式 3,式 4から、 f、 1を求める処理にっ 、ても、並列実行可能なことがわかる。
例えば、固有値算出手段 17aは、算出した固有ベクトル qを用いて、 f、 1の一部の要 素を算出し、固有値算出手段 17bは、算出した固有ベクトル αを用いて、 f、 1の残りの 要素を算出してもよい。その後、固有値算出手段 17a、あるいは固有値算出手段 17 bが、固有値算出手段 17a、 17bの算出した f、 1の各要素を統合することによって、式 3,式 4で示される f、 1を算出することができる。
[0241] ここでは、 zの成分に 0が含まれず、 D の対角成分に同じ値が含まれない場合にお
12
ける固有値算出手段 17a、 17bの処理について説明した力 zの成分に 0が含まれる 場合、 D の対角成分に同じ値が含まれる場合であっても、固有方程式を解く処理、
12
固有ベクトル ^を算出する処理、固有ベクトル力 分割元の行列の行列要素を算出 する処理を同様に並列実行することができうる。
[0242] 上記のような各構成要素における並列処理において、各手段が固有値等の情報を 格納するメモリを用いる場合に、複数の手段が同一のメモリ、すなわち共有メモリを使 用してもよぐあるいは、各手段がそれぞれ別のメモリを使用してもよい。
[0243] なお、図 24〜図 26を用いて、固有値分解部 15や固有値算出部 17等が 2個の手 段によって並列処理を実行する場合について説明したが、固有値分解部 15や固有 値算出部 17等が 3以上の手段を備え、並列処理を実行してもよい。また、固有値分 解部 15、固有値算出部 17、コレスキー分解部 22、ベクトル算出部 23,コレスキー分 解部 31、ベクトル算出部 32のそれぞれにおいて並列処理が実行される場合につい て説明したが、いずれかの任意の 1以上の部において、並列処理が実行されなくても よい。例えば、固有値分解部 15において並列処理が行われなくてもよい。
[0244] また、その並列処理は、 1の装置において、 2以上の CPU等を用いて実行されても よぐ 2以上の装置において実行されてもよい。例えば、図 27で示されるように、装置 Aと、装置 Bとがそれぞれ固有値分解手段 15a、 15bを備え、各装置において、固有 値分解の処理が並列実行されてもよい。この場合には、装置 Aの固有値分解手段 1 5aを備える固有値分解部 15—1と、装置 Bの固有値分解手段 15bを備える固有値分 解部 15— 2とによって固有値分解部 15が構成されることになる。したがって、固有値 分解装置 1は、装置 Aと、装置 Bとからなるシステムを構成することになる。ここでは、 固有値分解部 15の並列処理にっ 、て説明したが、その他の固有値算出部 17ゃコ レスキー分解部 22等についても、 2以上の装置による並列処理を行ってもよい。
[0245] 次に、数値実験による性能の評価について説明する。ここでは、本実施の形態によ る固有値分解装置 1の方法 (DDC)と、標準的な分割統治法 (DC)と、 QR法 (QR)と 、二分法と逆反復法とを組み合わせた方法 (B&I)とを比較した。なお、本実施の形態 による固有値分解装置 1の方法では、コレスキー分解を qd型ツイスト分解法で行い、 逆反復（図 12のステップ S701参照）を 1回用いた。分割統治法としては、 LAPACK で提供されている DSTEDCを用いた。 QR法としては、 LAPACKで提供されている DSTEQRを用いた。二分法と逆反復法とを組み合わせた方法としては、 LAPACK で提供されて 、る DSTEVRのソースコードを一部変更したものを用いた。実験で用 いる LAPACKは BLAS (Basic Linear Algebra Subprograms)を呼び出す。 また、本実施の形態による固有値分解装置 1の固有値分解部 15では、 QR法を用い て固有値分解を行った。また、この数値実験では、ランダムに与えた対角行列を、ラ ンチヨス法を用いて対称 3重対角化した行列を用いた。数値実験結果のグラフにお いて、横軸は、数値実験で用いた行列の大きさである。 [0246] 図 28は、計算の実行時間を示す図である。実行時間が短いほど、高速に計算でき ることになる。図 28の数値実験結果のグラフからわ力るように、本実施の形態による 固有値分解装置 1の方法では、他の方法に比べて高速に固有値分解を行うことがで さることがゎカゝる。
[0247] 図 29は、算出された固有ベクトルの直交性を示す図である。値の小さい方力 直交 性が高いことになる。図 29の数値実験結果のグラフからわ力るように、本実施の形態 による固有値分解装置 1の方法では、他の方法に比べて、固有ベクトルの直交性は 高くない。本実施の形態による固有値分解装置 1の方法では、二分法と逆反復法と を組み合わせた方法と同程度の直交性を有することがわかる。これは、固有値をまず 算出し、その後に、固有ベクトルを算出していることに起因していると考えられる。
[0248] 図 30は、全固有値の相対誤差の総和を示す図である。値の小さい方が、誤差が小 さぐ精度が高いことになる。図 30の数値実験結果のグラフからわ力るように、本実施 の形態による固有値分解装置 1の方法では、他の方法に比べて固有値の精度が高 いことがわかる。分割統治法により算出した固有値と同程度の固有値の精度を有す る。本実施の形態による固有値分解装置 1の方法では、分割統治法によって固有値 を算出しているため、当然の結果であると考えられる。
[0249] 図 31は、前固有ベクトルの誤差の総和を示す図である。値の小さい方力 誤差が 小さぐ精度が高いことになる。図 31の数値実験結果のグラフからわ力るように、本実 施の形態による固有値分解装置 1の方法では、他の方法に比べて固有ベクトルの精 度が高いことがわかる。
[0250] このように、本実施の形態による固有値分解装置 1の方法では、高速に固有値分解 をすることができ、さらに、固有値、固有ベクトルの精度の高いことがわかる。また、前 述のように、並列性にも優れているため、並列処理を行うことによってさらに高速ィ匕を 図ることも可能である。
[0251] [主成分分析による顔画像認識]
次に、主成分分析による顔画像認識に固有値分解を応用する場合について説明 する。
[0252] 主成分分析による顔画像認識の処理は、以下の各ステップによって行われる。 (1)射影行列の作成
(2)個人の顔画像に関する情報の登録
(3)入力された顔画像と、登録されている個人の顔画像に関する情報とを比較するこ とによる認識
[0253] まず、上記（1)の射影行列の作成について説明する。
種々の人物の顔画像データを用意する。その顔画像データは、例えば、スキャナ やデジタルスチルカメラ、デジタルビデオカメラ等の光学的読み取り機器によって読 み取られた 2次元画像データであってもよい。各顔画像データは、各画素値を成分と するベクトルで表示することができる。そのベクトルを、次のようにする。ここでは、 M 個の顔画像データを扱うものとする。
[数 72]
Λ1 ' X2 ' ' M
[0254] 次に、その顔画像データのベクトルの平均値 を算出する。
[数 73]
Μ
[0255] その算出した; zを用いて、共分散行列 Cを次のように算出する。
[数 74]
M M
C = ∑∑( - - ) τ
[0256] その後、上述の固有値分解装置 1によって、その行列 Cに対して、固有値分解を行 う。すなわち、行列記憶部 11に行列 Cを蓄積し、その行列 Cに対して対角化や、行列 の分割等の処理を行うことによって、固有値分解を行う。その固有値分解された結果 である固有値えは、固有値記憶部 18に蓄積されることになる。また、固有値分解さ れた結果である固有ベクトル Vは、固有ベクトル記憶部 20で記憶されて 、る固有べク トルに対して、対角化で用いた直交行列を作用させることによって、前述のようにして 求めることができる。なお、添字 jは、固有値の降順になるように設定されているものと する。すなわち、固有値え が最大の固有値となる。 [0257] その求めた固有ベクトルのうち、対応する固有値が上位から d個の固有ベクトルを 用いて、射影行列 r?を次のように算出する。
[数 75]
= (v^ v, , - - - ^, )
[0258] 次に、上記（2)の個人の顔画像に関する情報の登録について説明する。
各個人について、 m個の顔画像データ
[数 76]
を用意する。ここで、添字 kは、個人を識別するための添字である。そのベクトルに射 影行列 7?をそれぞれ掛けることにより、次式のようにして m個の d次元ベクトル ξを算 出し、その平均値 算出する。そして、その平均値/ を登録しておく。
[数 77]
i = η i = l,2, " ', m
1
=丄∑
[0259] 次に、上記（3)の入力された顔画像と、登録されている個人の顔画像に関する情報 とを比較することによる認識について説明する。
ある人物の顔画像データのベクトル Xが入力されたとする。すると、その顔画像デー タのベクトル Xに射影行列 7?を掛けることにより、 d次元ベクトル を算出する。
[数 78]
[0260] 次に、その d次元ベクトル ξと、各個人の d次元ベクトル ξの平均値 kとの差のノ ルムを計算する。
[数 79]
I - <
[0261] そのノルムがあら力じめ設定されて ヽるしき 、値より小さ ヽ場合に、入力された顔画 像力 添字 kで識別される個人の顔画像であると判断してもよい。あるいは、すべての 添字 kについて上記のノルムを計算し、入力された顔画像が、最小のノルムに対応す る添字 kで識別される個人の顔画像であると判断してもよい。このようにして、固有値 分解装置 1を用いて、主成分分析による顔画像認識を行うことができうる。
[0262] [2次元画像から 3次元へ復元する画像処理への応用]
次に、物体の 2次元画像から 3次元へ復元する画像処理に固有値分解を応用する 場合について説明する。
[0263] 複数の 2次元画像から 3次元復元を行う処理は、以下の各ステップによって行われ る。
(1) 2次元画像力 特徴点を抽出するステップ。
(2)特徴点データより形状 (元の物体の特徴点の 3次元座標データ)及び回転（3次 元データ力 特徴点データへの変換）に関するデータを計算するステップ。
(3)形状及び回転のデータより可視化を行うステップ。
[0264] 以下、上記（1)、（2)の各ステップについて、図 32のフローチャートを用いて説明す る。ここで、 2次元画像は、例えば、スキャナやデジタルスチルカメラ、デジタルビデオ カメラ等の光学的読み取り機器によって読み取られた 2次元画像であってもよい。
[0265] ステップ S5001において、 2次元画像； j (j = l, · · · , m、 mは 3以上の整数）から特 徴点 i (i= l, · · · , n、nは 2以上の整数)の座標 (xj, yj)を抽出する。取り扱う 2次元 画像は、弱中心射影画像であることが好ましい。このとき、次の式が成り立つ。
[数 80]
Figure imgf000061_0001
、 'ノ
[0266] ここで、 sは物体のスケールに相対する j番目の画像のスケール、 r , r はそれぞ
j l, j 2, j
れ物体座標系に相対する j番目のカメラ座標系の回転行列の 1番目と 2番目の行べク トル、（Xi, Yi, Zi) 番目の点の 3次元座標である。物体のスケールは 1番目の画 像のスケールと同じにし (s = 1)、物体の座標系の姿勢は 1番目の画像のカメラ座標 系と同じにする (r = (1, 0, 0)T, r = (0, 1, 0)T)。 m枚全ての画像における η個 全ての点の座標を行列 Dの要素として並べると、次式のようになる。
[0267] [数 81] D = MS
{ 1
ただし、 D = y{
、
Figure imgf000062_0001
Figure imgf000062_0002
[0268] Mと Sの形から分かるように、 Dのランクは 3である。ここで、ステップ S5001において
、行列 Dが与えられている。以下、回転に関するデータ M及び形状 Sを求める。
[0269] そこで、行列 Dの特異値分解
D = U∑VT
を考える。ここで、∑は特異値を大小順に対角線上に並べたもので、 U及び Vはそれ ぞれ左直交行列、及び右直交行列である。この特異値分解において、前述の固有 値分解を用いることができる。すなわち、固有値分解装置 1において、行列記憶部 11 が記憶している対称行列 Aを、 DTDとすることにより、前述の説明のようにして、対称 行列 Aの固有値分解がなされることになる。なお、固有値分解装置 1において、固有 ベクトル記憶部 20に蓄積されるのは、対称行列 DTDを変換した対称 3重対角行列 T の固有ベクトルであるので、その固有ベクトルを、前述の説明のようにして、行列 DTD の固有ベクトルに変換する必要がある。また、その行列 DTDの各固有値の 1Z2乗（ 平方根)が各特異値となる。すなわち、特異値 σ = ( λ ) 1/2となる。また、固有べタト ルと、右直交行列 Vの各列を構成する右特異ベクトル Xとは等しいため、固有べタト
j
ルを各列に有する行列力 右直交行列 Vとなる。なお、特異値 σ力^でない場合に は、左直交行列 Uの各列を構成する左特異ベクトル yを次のようにして求めることが できる。 [0270] [数 82]
となる。一方、特異値 σ力^である場合には、
[数 83]
が = 0 を解くことにより、 yを算出することができる。ここで、次の置き換えを行うことによって、 左特異ベクトル yを正規化する。
[数 84]
Figure imgf000063_0001
[0271] このようにして、左特異ベクトル yを求めることができ、この左特異ベクトル yを用い j ] て、 U= (y , y , · ··, y )とすることにより、左直交行列 Uを算出することができる。
[0272] ここで、画像のデジタル誤差のため、ゼロでな!、特異値は 3つ以上出てくる。しかし 、 4番目以降の特異値はノイズによるもので、最初の 3つの特異値と比べて格段に小 さい。そこで、ステップ S5002では、最初の 3つの特異値に対して特異ベクトルを計 算する。すなわち、本実施の形態による固有値分解装置 1では、最初の 3個の固有 値に対して固有ベクトルを計算する。採用する 3個のベクトルをまとめると、次式となる
[0273] D' = L'∑'R'T=M'S'
ただし、 M' =L' (∑') 1/2、 S' = (∑') 1/2R'T、 D'は II D-D' IIを最小にするランク 3 の行列である。
[0274] 次に、 D'から M及び Sを求めたいが、その組合せは唯一ではない。なぜなら、任意 の正則行列じが
D' = (M'C) (C_1S，）
を満たす力もである。そこで、上式における Cを M = M'Cを満たすように決める。 ま 下記の式を満たす。
[数 85]
Figure imgf000064_0001
[0275] E = CC^すると、上式力 Eの 6つの要素に関する 2m+ 1個の線形方程式が得ら れる。 m≥ 3であるので、 Eの要素を一意に決めることができる。ステップ S5003にお いて、行列 Eを求める。
[0276] 次に、ステップ S5004にお!/、て、 Eから Cを求める。じの自由度（9)は Eの自由度（6 )より多い。そこで、条件 = (1, 0, 0)T, r = (0, 1, 0)τを加えれば、 Cを決めること lj ¾
ができる。このとき 2つの解（Necker Reversal)が出る。
次に、ステップ S5005において、 M = M'C及び S = C_1S'より、回転に関するデー タ M及び形状 Sが決まる。
[0277] [文書検索への応用]
次に、文書検索処理に固有値分解を応用する場合について説明する。文書中から その文書の内容に関連する索引語を抽出し、索引語の重みを計算する処理の後、 ベクトル空間モデルでは、この索引語の重みを要素とするベクトルで文書を表現する 。ここで、検索対象となる文書を d , d , · · · , dとし、これら文書集合全体を通して全
1 2 n
部で m個の索引語 w , w , · · · , wがあるとする。このとき、文書 dは、次のようなベタ
1 2 m j
トルで表現されることになる。これを文書ベクトルと呼ぶ。
[数 86]
Figure imgf000064_0002
[0278] ここで、 dは索引語 wの文書 dにおける重みである。また、文書集合全体は、次の ような m X n行列 Dによって表現することができる。
[数 87]
Figure imgf000065_0001
[0279] 行列 Dを索引語文書行列と呼ぶ。索引語文書行列の各列は文書に関する情報を 表している文書ベクトルである力 同様に、索引語文書行列の各行は索引語に関す る情報を表しているべ外ルであり、これを索引語べ外ルと呼ぶ。検索質問も、文書と 同様に、索引語の重みを要素とするベ外ルで表現することができる。検索質問文に 含まれる索引語 wの重みを qとすると、検索質問ベクトル qは次のように表されること になる。
[数 88]
Figure imgf000065_0002
[0280] 実際の文書検索においては、与えられた検索質問文と類似した文書を見つけ出す 必要があるが、検索質問ベクトル qと各文書べ外ル dの間の類似度を計算することに より行う。ベクトル間の類似度の定義としてはさまざまなものが考えられるが、文書検 索にお!、てよく用いられて 、るものはコサイン尺度（2つのベクトルのなす角度）また は内積である。
[0281] [数 89] d , q ∑ '
cos(d q) = (コサイン尺度
Figure imgf000065_0003
' q =∑ (内積 なお、ベクトルの長さが 1に正規ィ匕 (コサイン正規化）されている場合には、コサイン 尺度と内積とは一致する。
図 33は、本実施の形態による固有値分解装置 1を利用した文書検索方法の一例 を示すフローチャートである。 ステップ S6001にお!/、て、質問ベクトル qを受け取る。
[0283] ここでは、 Dの近似行列 Dを使った検索を考える。ベクトル空間モデルでは、検索
k
質問ベクトル qと索引語文書行列 D中の各文書ベクトル dの間の類似度を計算するこ とにより検索を行うが、ここでは Dの代わりに Dを使う。ベクトル空間モデルでは、文
k
書ベクトルの次元数は索引語の総数と等しい。したがって、検索対象となる文書の数 が増えるに従い、文書ベクトルの次元数も増加する傾向にある。しかし、次元数が増 カロしてくると、コンピュータのメモリによる制限や検索時間の増大などの問題が生じて くるばかりでなぐ文書中に含まれる不必要な索引語がノイズ的な影響を及ぼし、検 索精度を低下させてしまうという現象も起こってくる。潜在的意味インデキシング (late nt semantic indexing ; LSI)は、高次元の空間にある文書ベクトルを低次元の空 間へと射影することにより、検索精度の改善を図る技術である。高次元の空間では別 々に扱われていた索引語力 低次元の空間では相互に関連を持ったものとして扱わ れる可能性もあるため、索引語の持つ意味や概念に基づく検索を行うことができる。 たとえば、通常のベクトル空間モデルでは" car"という索引語と" automobile"という 索引語はまったく別物であり、一方の索引語による質問ではもう片方の索引語を含ん だ文書を検索することができない。しかし、低次元の空間ではこれらの意味的に関連 した索引語は 1つの次元に縮退することが期待できるため、 "car"という検索質問によ つて" car"を含む文書ば力りでなぐ automobile"を含む文書をも検索することが可能 となる。潜在的意味インデキシングでは、特異値分解により高次元ベクトルの次元削 減を行うが、これは基本的に多変量解析における主成分分析と等価である。
[0284] ステップ S6002にお!/、て、 kを選択する。ここでは、 k〈rなる kの値を選択する。 r= min (n, m)である。 kの値は、予め与えられていてもよいし、計算ごとに選択可能で あってもよい。
[0285] 次に、ステップ S6003では、行列 Dの特異値分解を行う。この特異値分解として、 前述の 2次元画像から 3次元へ復元する画像処理への応用で説明したように、固有 値分解による方法を用いることができる。すなわち、固有値分解装置 1において、行 列記憶部 11が記憶している対称行列 Aを、このたびの行列 DTDとすることにより、前 述の説明のようにして、対称行列 Aの固有値分解がなされ、その結果として、行列 D の特異値分解を行うことができうる。なお、この特異値分解では、計算された特異値 のうち、大きい順に 1番目力も k番目までの k個の特異値に対して Dの特異ベクトルを 算出する。すなわち、固有値分解において、計算された固有値のうち、大きい順に 1 番目から k番目までの k個の固有値に対して固有ベクトルを算出し、この固有ベクトル を用いて特異ベクトルを算出すればよい。 kは、ステップ S6002で選択された値であ る。
[0286] すなわち、
D =U∑V T
k k k
なる u及び Vを計算する。ここで、 Uは、最初の k個の左特異ベクトルのみから構成 k k k
される m X k行列であり、 Vは、最初の k個の右特異ベクトルのみから構成される n X
k
k行列であり、∑は、最初の k個の特異値のみ力も構成される k X k対角行列である。
[0287] 次に、ステップ 6004にお 、て、行列 Dと質問ベクトル qとの類似度を計算する。 Vヽ
k
ま、ベクトル eを n次元の単位ベクトルとすると、 Dの j番目の文書ベクトルは D eで表
j k k j すことができる。文書ベクトル D eと検索質問べ外ル qとの間の類似度計算は、
[数 90]
Figure imgf000067_0001
φ ejfo)
Figure imgf000067_0002
としてもよいが、別の定義を用いてもよい。上式では、 Dを U , ∑ , Vから再構成す
k k k k
る必要はなく特異値分解の結果から、直接、類似度を計算できることを示している。 上式の中に現われる∑ V Teは、
k k j
∑ V Te =U TD e
k k j k k j
と書き直すことができる。この式の右辺は、近似行列 Dにおける j番目の文書ベクトル
k
の基底 Uのもとでの座標（文書の k次元表現）を表している。同様に、上式の中の U
k k
Tqは、検索質問ベクトル qの基底 Uのもとでの座標 (検索質問の k次元表現)である。 [0288] ステップ S6005において、ステップ S6004において計算された類似度を基準に、 検索結果を出力する。
このように、本実施の形態による固有値分解装置 1は、画像処理や検索処理等の 種々の処理に対して用いることができうる。なお、上記以外の処理に用いてもよいこと は言うまでもない。
[0289] 以上のように、本実施の形態による固有値分解装置 1では、分割統治法を用いて 固有値のみを算出するため、標準的な分割統治法で実行されるベクトル更新が必要 なくなり、標準的な分割統治法よりも非常に高速である。また、固有値力 固有べタト ルを算出する処理も、高速に処理することができる。さらに、 QR法では固有値を算出 する段階の並列化が困難であるが、本実施の形態による固有値分解装置 1では、固 有値を算出する段階、及び固有値力 固有ベクトルを算出する段階のそれぞれにお いて、本質的に高い並列性を持つ。また、本実施の形態による固有値分解装置 1で は、他の方法と比べて高 、精度の得られることがわかる。
[0290] なお、本実施の形態による固有値分解装置 1は、スタンドアロンの装置であってもよ ぐサーバクライアントシステムにおけるサーバ装置であってもよぐ前述のように、複 数の装置カゝら構成されるシステムであってもよい。固有値分解装置 1がサーバクライ アントシステムにおけるサーバ装置である場合には、行列記憶部 11で記憶される対 称行列 Aや、固有値記憶部 18で記憶される固有値や固有ベクトル記憶部 20で記憶 される固有ベクトル等は、インターネットやイントラネット等の通信回線を介して送受信 されてちょい。
[0291] また、本実施の形態による固有値分解装置 1における固有ベクトルの算出において 、ある固有値を用いた固有ベクトルの算出では、 qd型ツイスト分解法を用い、他の固 有値を用いた固有ベクトルの算出では、 LV型ツイスト分解法を用いるようにしてもよ い。例えば、値の近接している固有値については、 LV型ツイスト分解法を用いて固 有ベクトルの算出を行い、値が近接していない固有値については、 qd型ツイスト分解 法を用いて固有ベクトルの算出を行ってもょ 、。固有値の値が近接して 、るかどうか は、所定の記録媒体等にぉ 、てあら力じめ設定されて 、るしき 、値と比較することに よって半 U断してもよ ヽ。 [0292] また、上記実施の形態にお!、て、各処理または各機能は、単一の装置または単一 のシステムによって集中処理されることによって実現されてもよぐあるいは、複数の 装置または複数のシステムによって分散処理されることによって実現されてもよい。
[0293] また、上記実施の形態において、各構成要素は専用のハードウェアにより構成され てもよく、あるいは、ソフトウェアにより実現可能な構成要素については、プログラムを 実行することによって実現されてもよい。例えば、ハードディスクや半導体メモリ等の 記録媒体に記録されたソフトウェア 'プログラムを CPU等のプログラム実行部が読み 出して実行することによって、各構成要素が実現され得る。なお、上記実施の形態に おける固有値分解装置を実現するソフトウェアは、以下のようなプログラムである。つ まり、このプログラムは、コンピュータに、対称 3重対角行列 Tが記憶される対角行列 記憶部から前記対称 3重対角行列 Tを読み出し、当該対称 3重対角行列 Tを 2個の 対称 3重対角行列に分割して前記対角行列記憶部に蓄積し、その対称 3重対角行 列を 2個の対称 3重対角行列に分割して前記対角行列記憶部に蓄積する処理を、分 割後の各対称 3重対角行列があらかじめ決められた大きさ以下となるまで繰り返す行 列分割ステップと、前記あらかじめ決められた大きさ以下の各対称 3重対角行列を前 記対角行列記憶部から読み出し、前記各対称 3重対角行列に対して固有値分解を 行い、前記各対称 3重対角行列の固有値と、前記各対称 3重対角行列の固有べタト ルカ なる直交行列の行列要素とを少なくとも算出し、当該固有値と当該行列要素と を固有値分解記憶部に蓄積する固有値分解ステップと、各対称 3重対角行列の固有 値と、行列要素とを前記固有値分解記憶部から読み出し、前記固有値と、前記行列 要素とから、分割元の対称 3重対角行列の固有値と、分割元の対称 3重対角行列の 行列要素とを算出して前記固有値分解記憶部に蓄積し、その分割元の対称 3重対 角行列の固有値と、行列要素とを算出する処理を、対称 3重対角行列 Tの少なくとも 1個の固有値を算出するまで繰り返し、前記対称 3重対角行列 Tの少なくとも 1個の固 有値を固有値記憶部に蓄積する固有値算出ステップと、前記対角行列記憶部から 前記対称 3重対角行列 Tを読み出し、前記固有値記憶部から前記対称 3重対角行 列 Tの固有値を読み出し、前記対称 3重対角行列 Tとその固有値とから、ツイスト分 解法を用いて前記対称 3重対角行列 Tの少なくとも 1個の固有ベクトルを算出して固 有ベクトル記憶部に蓄積する固有ベクトル算出ステップと、を実行させるためのもの である。
また、このプログラムでは、前記固有ベクトル算出ステップにおいて、 qd型ツイスト分 解法により固有ベクトルを算出してもよ 、。
[0294] また、このプログラムでは、前記固有ベクトル算出ステップ力 前記固有値記憶部か ら前記対称 3重対角行列 Tの固有値を読み出し、当該固有値のいずれかが負の値 である場合に、前記対角行列記憶部から前記対称 3重対角行列 Tを読み出し、当該 対称 3重対角行列 Tの固有ベクトルを変化させることなぐ当該対称 3重対角行列丁の すべての固有値が正の値となるように、当該対称 3重対角行列 Tを正定値化し、正定 値ィ匕した後の固有値を前記固有値記憶部に蓄積し、正定値化した後の対称 3重対 角行列 Tを前記対角行列記憶部に蓄積する正定値化ステップと、前記対角行列記 憶部からすべての固有値が正の値である対称 3重対角行列 Tを読み出し、前記固有 値記憶部から前記対称 3重対角行列 Tの正の値の固有値を読み出し、前記対称 3重 対角行列 Tの各要素に関して qd型変換を行うことによって、前記対称 3重対角行列 T を上 2重対角行列 B(+1)及び下 2重対角行列 B(_1)にコレスキー分解するコレスキー 分解ステップと、前記上 2重対角行列 B(+1)及び下 2重対角行列 B(_ 1)の各要素と、前 記対称 3重対角行列 Tの正の値の固有値とを用いて固有ベクトルを算出して前記固 有ベクトル記憶部に蓄積するベクトル算出ステップと、を備えてもよい。
また、このプログラムでは、前記固有ベクトル算出ステップにおいて、 LV型ツイスト 分解法により固有ベクトルを算出してもよ 、。
[0295] また、このプログラムでは、前記固有ベクトル算出ステップ力 前記対角行列記憶部 から前記対称 3重対角行列 Tを読み出し、前記固有値記憶部から前記対称 3重対角 行列 Tの固有値を読み出し、前記対称 3重対角行列 Tの各要素に関してミウラ変換、 dLVv型変換、逆ミウラ変換を行うことによって、前記対称 3重対角行列 Tを上 2重対 角行列 B(+1)及び下 2重対角行列 B(_1)にコレスキー分解するコレスキー分解ステップ と、前記上 2重対角行列 B(+1)及び下 2重対角行列 B(_1)の各要素と、前記対称 3重 対角行列 Tの固有値とを用いて固有ベクトルを算出して前記固有ベクトル記憶部に 蓄積するベクトル算出ステップと、を備えてもよい。 [0296] また、このプログラムでは、コンピュータに、対称行列 Aが記憶される行列記憶部か ら前記対称行列 Aを読み出し、前記対称行列 Aを 3重対角化した前記対称 3重対角 行列 Tを算出して前記対角行列記憶部に蓄積する対角ィ匕ステップをさらに実行させ てもよい。
なお、上記プログラムにおいて、情報を蓄積するステップなどでは、ハードウェアで しか行われな 、処理は少なくとも含まれな!/、。
[0297] また、このプログラムは、サーバなど力 ダウンロードされることによって実行されても よぐ所定の記録媒体 (例えば、 CD— ROMなどの光ディスクや磁気ディスク、半導 体メモリなど）に記録されたプログラムが読み出されることによって実行されてもよい。
[0298] また、このプログラムを実行するコンピュータは、単数であってもよぐ複数であって もよい。すなわち、集中処理を行ってもよぐあるいは分散処理を行ってもよい。
[0299] 図 34は、上記プログラムを実行して、上記実施の形態による固有値分解装置 1を実 現するコンピュータの外観の一例を示す模式図である。上記実施の形態は、コンビュ ータハードウェア及びその上で実行されるコンピュータプログラムによって実現されう る。
[0300] 図 34において、コンピュータシステム 100は、 CD— ROM (Compact Disk Rea d Only Memory)ドライブ 105、 FD (Flexible Disk)ドライブ 106を含むコンビュ ータ 101と、キーボード 102と、マウス 103と、モニタ 104とを備える。
[0301] 図 35は、コンピュータシステムを示す図である。図 35において、コンピュータ 101は 、 CD— ROMドライブ 105、 FDドライブ 106に加えて、 CPU (Central Processing Unit) 111と、ブートアッププログラム等のプログラムを記憶するための ROM (Rea d Only Memory) 112と、 CPU111に接続され、アプリケーションプログラムの命 令を一時的に記憶すると共に、一時記憶空間を提供する RAM (Random Access Memory) 113と、アプリケーションプログラム、システムプログラム、及びデータを 記憶するハードディスク 114と、 CPU111、 ROM112等を相互に接続するバス 115 とを備える。なお、コンピュータ 101は、 LANへの接続を提供する図示しないネットヮ ークカードを含んでもよい。
[0302] コンピュータシステム 100に、上記実施の形態による固有値分解装置 1の機能を実 行させるプログラムは、 CD— ROM121、または FD122に記憶されて、 CD— ROM ドライブ 105、または FDドライブ 106に挿入され、ハードディスク 114に転送されても よい。これに代えて、そのプログラムは、図示しないネットワークを介してコンピュータ 1 01に送信され、ハードディスク 114に記憶されてもよい。プログラムは実行の際に RA M113にロードされる。なお、プログラムは、 CD— ROM121や FD122、またはネット ワークから直接、ロードされてもよい。
[0303] また、行列記憶部 11、対角行列記憶部 13、固有値分解記憶部 16、固有値記憶部 18,固有ベクトル記憶部 20は、 RAMI 13やハードディスク 114によって実現されて ちょい。
[0304] プログラムは、コンピュータ 101に、上記実施の形態による固有値分解装置 1の機 能を実行させるオペレーティングシステム (OS)、またはサードパーティプログラム等 を必ずしも含まなくてもよい。プログラムは、制御された態様で適切な機能 (モジユー ル)を呼び出し、所望の結果が得られるようにする命令の部分のみを含んでいてもよ い。コンピュータシステム 100がどのように動作するのかについては周知であり、詳細 な説明を省略する。
また、本発明は、以上の実施の形態に限定されることなぐ種々の変更が可能であ り、それらも本発明の範囲内に包含されるものであることは言うまでもない。
また、本発明のほんのいくつかの典型的な実施例について上で詳細に説明したが 、その典型的な実施例において、発明の利益と新規な技術から実質的にはずれるこ となく多くの変更が可能であることを当業者は容易に認識することができるであろう。 したがって、そのようなすべての変更は、本発明の範囲に含まれるものである。
産業上の利用可能性
[0305] 以上のように、本発明による固有値分解装置等によれば、固有値分解を高速に処 理することができ、分子軌道法による化学計算処理や統計計算処理、情報検索処理 、その他の固有値分解を用いる処理を実行する装置等にぉ 、て有用である。
図面の簡単な説明
[0306] [図 1]本発明の実施の形態 1による固有値分解装置の構成を示すブロック図
[図 2]同実施の形態による固有値分解装置の動作を示すフローチャート 圆 3]同実施の形態による固有値分解装置の動作を示すフローチャート
圆 4]同実施の形態による固有値分解装置の動作を示すフローチャート
圆 5]同実施の形態における対称 3重対角行列 Tの分割について説明するための図 圆 6]同実施の形態における対称 3重対角行列 Tの固有値の算出について説明する ための図
圆 7]同実施の形態による固有値分解装置の動作を示すフローチャート
圆 8]同実施の形態による固有値分解装置の動作を示すフローチャート
圆 9]同実施の形態における対称 3重対角行列 Tの固有値の算出について説明する ための図
[図 10]同実施の形態における固有ベクトル算出部の構成を示すブロック図 圆 11]同実施の形態による固有値分解装置の動作を示すフローチャート 圆 12]同実施の形態による固有値分解装置の動作を示すフローチャート
[図 13]同実施の形態における固有ベクトル算出部の構成を示すブロック図 圆 14]同実施の形態による固有値分解装置の動作を示すフローチャート
[図 15]同実施の形態におけるコレスキー分解について説明するための図 圆 16]同実施の形態による固有値分解装置の動作を示すフローチャート 圆 17]同実施の形態による固有値分解装置の動作を示すフローチャート 圆 18]同実施の形態による固有値分解装置の動作を示すフローチャート 圆 19]同実施の形態による固有値分解装置の動作を示すフローチャート 圆 20]同実施の形態による固有値分解装置の動作を示すフローチャート 圆 21]同実施の形態による固有値分解装置の動作を示すフローチャート
[図 22]同実施の形態における qd型ツイスト分解法と LV型ツイスト分解法とについて 説明するための図
圆 23]同実施の形態による固有値分解装置の動作を示すフローチャート 圆 24]同実施の形態による固有値分解装置の構成の他の一例を示すブロック図
[図 25]同実施の形態における固有ベクトル算出部の構成の他の一例を示すブロック 図
[図 26]同実施の形態における固有ベクトル算出部の構成の他の一例を示すブロック 図
圆 27]同実施の形態による固有値分解装置の構成の他の一例を示すブロック図 圆 28]同実施の形態による固有値分解装置の数値実験の結果を示す図 圆 29]同実施の形態による固有値分解装置の数値実験の結果を示す図 圆 30]同実施の形態による固有値分解装置の数値実験の結果を示す図 圆 31]同実施の形態による固有値分解装置の数値実験の結果を示す図
[図 32]同実施の形態における画像処理の一例を示すフローチャート
[図 33]同実施の形態における文書検索処理の一例を示すフローチャート 圆 34]同実施の形態におけるコンピュータシステムの外観一例を示す模式図 [図 35]同実施の形態におけるコンピュータシステムの構成の一例を示す図

Claims

請求の範囲
対称 3重対角行列 Tが記憶される対角行列記憶部と、
前記対角行列記憶部から前記対称 3重対角行列 Τを読み出し、当該対称 3重対角 行列 Τを 2個の対称 3重対角行列に分割して前記対角行列記憶部に蓄積し、その対 称 3重対角行列を 2個の対称 3重対角行列に分割して前記対角行列記憶部に蓄積 する処理を、分割後の各対称 3重対角行列があらかじめ決められた大きさ以下となる まで繰り返す行列分割部と、
前記あらかじめ決められた大きさ以下の各対称 3重対角行列の固有値と、当該各対 称 3重対角行列の固有ベクトル力 なる直交行列の一部の要素である行列要素とが 記憶される固有値分解記憶部と、
前記あらかじめ決められた大きさ以下の各対称 3重対角行列を前記対角行列記憶部 から読み出し、前記各対称 3重対角行列に対して固有値分解を行い、前記各対称 3 重対角行列の固有値と、前記各対称 3重対角行列の固有ベクトルからなる直交行列 の行列要素とを少なくとも算出し、当該固有値と当該行列要素とを前記固有値分解 記憶部に蓄積する固有値分解部と、
前記対称 3重対角行列 Τの固有値が記憶される固有値記憶部と、
各対称 3重対角行列の固有値と、行列要素とを前記固有値分解記憶部から読み出 し、前記固有値と、前記行列要素とから、分割元の対称 3重対角行列の固有値と、分 割元の対称 3重対角行列の行列要素とを算出して前記固有値分解記憶部に蓄積し 、その分割元の対称 3重対角行列の固有値と、行列要素とを算出する処理を、対称 3 重対角行列 Τの少なくとも 1個の固有値を算出するまで繰り返し、前記対称 3重対角 行列 Τの少なくとも 1個の固有値を前記固有値記憶部に蓄積する固有値算出部と、 前記対称 3重対角行列 Τの固有ベクトルが記憶される固有ベクトル記憶部と、 前記対角行列記憶部から前記対称 3重対角行列 Τを読み出し、前記固有値記憶部 から前記対称 3重対角行列 Τの固有値を読み出し、前記対称 3重対角行列 Τとその 固有値とから、ツイスト分解法を用いて前記対称 3重対角行列 Τの少なくとも 1個の固 有ベクトルを算出して前記固有ベクトル記憶部に蓄積する固有ベクトル算出部と、を 備えた固有値分解装置。
[2] 前記固有ベクトル算出部は、 qd型ツイスト分解法により固有ベクトルを算出する、請 求項 1記載の固有値分解装置。
[3] 前記固有ベクトル算出部は、
前記固有値記憶部から前記対称 3重対角行列 Tの固有値を読み出し、当該固有値 のいずれかが負の値である場合に、前記対角行列記憶部から前記対称 3重対角行 列 Tを読み出し、当該対称 3重対角行列 Tの固有ベクトルを変化させることなぐ当該 対称 3重対角行列 Tのすベての固有値が正の値となるように、当該対称 3重対角行 列 Tを正定値ィ匕し、正定値ィ匕した後の固有値を前記固有値記憶部に蓄積し、正定値 ィ匕した後の対称 3重対角行列 Tを前記対角行列記憶部に蓄積する正定値ィ匕部と、 前記対角行列記憶部からすべての固有値が正の値である対称 3重対角行列 Tを読 み出し、前記固有値記憶部から前記対称 3重対角行列 Tの正の値の固有値を読み 出し、前記対称 3重対角行列 Tの各要素に関して qd型変換を行うことによって、前記 対称 3重対角行列 Tを上 2重対角行列 B(+ 1)及び下 2重対角行列 B(_1)にコレスキー 分解するコレスキー分解部と、
前記上 2重対角行列 B(+1)及び下 2重対角行列 B(_1)の各要素と、前記対称 3重対角 行列 Tの正の値の固有値とを用いて固有ベクトルを算出して前記固有ベクトル記憶 部に蓄積するベクトル算出部と、を備えた、請求項 2記載の固有値分解装置。
[4] 前記固有ベクトル算出部は、 LV型ツイスト分解法により固有ベクトルを算出する、請 求項 1記載の固有値分解装置。
[5] 前記固有ベクトル算出部は、
前記対角行列記憶部から前記対称 3重対角行列 Tを読み出し、前記固有値記憶部 から前記対称 3重対角行列 Tの固有値を読み出し、前記対称 3重対角行列 Tの各要 素に関してミウラ変換、 dLVv型変換、逆ミウラ変換を行うことによって、前記対称 3重 対角行列 Tを上 2重対角行列 B(+ 1)及び下 2重対角行列 B(_1)にコレスキー分解する コレスキー分解部と、
前記上 2重対角行列 B(+1)及び下 2重対角行列 B(_1)の各要素と、前記対称 3重対角 行列丁の固有値とを用 、て固有ベクトルを算出して前記固有ベクトル記憶部に蓄積 するベクトル算出部と、を備えた、請求項 4記載の固有値分解装置。
[6] 前記コレスキー分解部は、複数のコレスキー分解手段を備え、
前記複数のコレスキー分解手段が、前記対称 3重対角行列 Tをコレスキー分解する 処理を並列実行する、請求項 3または請求項 5記載の固有値分解装置。
[7] 前記ベクトル算出部は、複数のベクトル算出手段を備え、
前記複数のベクトル算出手段が、固有ベクトルを算出する処理を並列実行する、請 求項 3、請求項 5または請求項 6記載の固有値分解装置。
[8] 前記固有値算出部は、複数の固有値算出手段を備え、
前記複数の固有値算出手段が、分割元の対称 3重対角行列の固有値と、行列要素 とを算出する処理を並列実行する、請求項 1から請求項 7のいずれか記載の固有値 分解装置。
[9] 前記固有値分解部は、複数の固有値分解手段を備え、
前記複数の固有値分解手段が、対称 3重対角行列に対して固有値分解を行う処理 を並列実行する、請求項 1から請求項 8のいずれか記載の固有値分解装置。
[10] 対称行列 Aが記憶される行列記憶部と、
前記対称行列 Aを前記行列記憶部から読み出し、前記対称行列 Aを 3重対角化した 前記対称 3重対角行列 Tを算出して前記対角行列記憶部に蓄積する対角化部と、を さらに備えた請求項 1から請求項 9のいずれか記載の固有値分解装置。
[11] 前記行列分割部は、対称 3重対角行列を略半分の 2個の対称 3重対角行列に分割 する、請求項 1から請求項 10のいずれか記載の固有値分解装置。
[12] 前記固有値算出部は、対称 3重対角行列 Tの全ての固有値を算出する、請求項 1か ら請求項 11の 、ずれか記載の固有値分解装置。
[13] 前記固有ベクトル算出部は、対称 3重対角行列 Tの全ての固有ベクトルを算出する、 請求項 12記載の固有値分解装置。
[14] 前記対称行列 Aは、顔画像データを示すベクトルの平均力 算出された共分散行列 である、請求項 10記載の固有値分解装置。
[15] 対称 3重対角行列 Tが記憶される対角行列記憶部と、行列分割部と、固有値分解部 と、前記固有値分解部によって固有値分解された固有値と固有べ外ルカゝらなる直交 行列の一部の要素である行列要素とが記憶される固有値分解記憶部と、前記対称 3 重対角行列 Tの固有値が記憶される固有値記憶部と、固有値算出部と、前記対称 3 重対角行列 Τの固有ベクトルが記憶される固有ベクトル記憶部と、固有ベクトル算出 部とを備えた固有値分解装置で用いられる固有値分解方法であって、
前記行列分割部が、前記対角行列記憶部から前記対称 3重対角行列 Τを読み出し、 当該対称 3重対角行列 Τを 2個の対称 3重対角行列に分割して前記対角行列記憶 部に蓄積し、その対称 3重対角行列を 2個の対称 3重対角行列に分割して前記対角 行列記憶部に蓄積する処理を、分割後の各対称 3重対角行列があらかじめ決められ た大きさ以下となるまで繰り返す行列分割ステップと、
前記固有値分解部が、前記あらかじめ決められた大きさ以下の各対称 3重対角行列 を前記対角行列記憶部から読み出し、前記各対称 3重対角行列に対して固有値分 解を行い、前記各対称 3重対角行列の固有値と、前記各対称 3重対角行列の固有 ベクトル力 なる直交行列の行列要素とを少なくとも算出し、当該固有値と当該行列 要素とを前記固有値分解記憶部に蓄積する固有値分解ステップと、
前記固有値算出部が、各対称 3重対角行列の固有値と、行列要素とを前記固有値 分解記憶部から読み出し、前記固有値と、前記行列要素とから、分割元の対称 3重 対角行列の固有値と、分割元の対称 3重対角行列の行列要素とを算出して前記固 有値分解記憶部に蓄積し、その分割元の対称 3重対角行列の固有値と、行列要素と を算出する処理を、対称 3重対角行列 Τの少なくとも 1個の固有値を算出するまで繰 り返し、前記対称 3重対角行列 Τの少なくとも 1個の固有値を前記固有値記憶部に蓄 積する固有値算出ステップと、
前記固有ベクトル算出部が、前記対角行列記憶部から前記対称 3重対角行列 Τを読 み出し、前記固有値記憶部から前記対称 3重対角行列 Τの固有値を読み出し、前記 対称 3重対角行列 Τとその固有値とから、ツイスト分解法を用いて前記対称 3重対角 行列 Τの少なくとも 1個の固有ベクトルを算出して前記固有ベクトル記憶部に蓄積す る固有ベクトル算出ステップと、を備えた固有値分解方法。
コンピュータに、
対称 3重対角行列 Τが記憶される対角行列記憶部力 前記対称 3重対角行列 Τを読 み出し、当該対称 3重対角行列 Τを 2個の対称 3重対角行列に分割して前記対角行 列記憶部に蓄積し、その対称 3重対角行列を 2個の対称 3重対角行列に分割して前 記対角行列記憶部に蓄積する処理を、分割後の各対称 3重対角行列があらかじめ 決められた大きさ以下となるまで繰り返す行列分割ステップと、
前記あらかじめ決められた大きさ以下の各対称 3重対角行列を前記対角行列記憶部 から読み出し、前記各対称 3重対角行列に対して固有値分解を行い、前記各対称 3 重対角行列の固有値と、前記各対称 3重対角行列の固有ベクトルからなる直交行列 の行列要素とを少なくとも算出し、当該固有値と当該行列要素とを固有値分解記憶 部に蓄積する固有値分解ステップと、
各対称 3重対角行列の固有値と、行列要素とを前記固有値分解記憶部から読み出 し、前記固有値と、前記行列要素とから、分割元の対称 3重対角行列の固有値と、分 割元の対称 3重対角行列の行列要素とを算出して前記固有値分解記憶部に蓄積し 、その分割元の対称 3重対角行列の固有値と、行列要素とを算出する処理を、対称 3 重対角行列 Tの少なくとも 1個の固有値を算出するまで繰り返し、前記対称 3重対角 行列 Tの少なくとも 1個の固有値を固有値記憶部に蓄積する固有値算出ステップと、 前記対角行列記憶部から前記対称 3重対角行列 Tを読み出し、前記固有値記憶部 から前記対称 3重対角行列 Tの固有値を読み出し、前記対称 3重対角行列 Tとその 固有値とから、ツイスト分解法を用いて前記対称 3重対角行列 Tの少なくとも 1個の固 有ベクトルを算出して固有ベクトル記憶部に蓄積する固有ベクトル算出ステップと、を 実行させるためのプログラム。
PCT/JP2007/051575 2006-08-08 2007-01-31 dispositif de décomposition de valeur propre et procédé de décomposition de valeur propre Ceased WO2008018188A1 (fr)

Priority Applications (2)

Application Number Priority Date Filing Date Title
US12/376,444 US8255447B2 (en) 2006-08-08 2007-01-31 Eigenvalue decomposition apparatus and eigenvalue decomposition method
JP2008528725A JP5017666B2 (ja) 2006-08-08 2007-01-31 固有値分解装置、及び固有値分解方法

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2006-215660 2006-08-08
JP2006215660 2006-08-08

Publications (1)

Publication Number Publication Date
WO2008018188A1 true WO2008018188A1 (fr) 2008-02-14

Family

ID=39032728

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2007/051575 Ceased WO2008018188A1 (fr) 2006-08-08 2007-01-31 dispositif de décomposition de valeur propre et procédé de décomposition de valeur propre

Country Status (3)

Country Link
US (1) US8255447B2 (ja)
JP (1) JP5017666B2 (ja)
WO (1) WO2008018188A1 (ja)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110262479A (zh) * 2019-05-28 2019-09-20 南京天辰礼达电子科技有限公司 一种履带式拖拉机运动学估计及偏差校准方法

Families Citing this family (17)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP4325877B2 (ja) * 2004-06-03 2009-09-02 独立行政法人科学技術振興機構 行列の高速高精度特異値分解方法、プログラムおよび装置
US8612445B2 (en) * 2009-05-13 2013-12-17 Hamid Hatami-Hanza System and method for a unified semantic ranking of compositions of ontological subjects and the applications thereof
JP2012234258A (ja) * 2011-04-28 2012-11-29 Sony Corp 画像処理装置と画像処理方法およびプログラム
US20120296637A1 (en) * 2011-05-20 2012-11-22 Smiley Edwin Lee Method and apparatus for calculating topical categorization of electronic documents in a collection
US8776250B2 (en) * 2011-07-08 2014-07-08 Research Foundation Of The City University Of New York Method of comparing private data without revealing the data
US11587172B1 (en) 2011-11-14 2023-02-21 Economic Alchemy Inc. Methods and systems to quantify and index sentiment risk in financial markets and risk management contracts thereon
US9495357B1 (en) * 2013-05-02 2016-11-15 Athena Ann Smyros Text extraction
US9262799B2 (en) * 2013-11-08 2016-02-16 Silicon Graphics International Corp. Shared memory eigensolver
JP6087848B2 (ja) * 2014-01-16 2017-03-01 日本電信電話株式会社 行列演算装置、行列演算方法、およびプログラム
EP3195146A4 (en) * 2014-09-18 2018-04-18 Empire Technology Development LLC Three-dimensional latent semantic analysis
US9965318B2 (en) * 2015-03-16 2018-05-08 Tata Consultancy Services Limited Concurrent principal component analysis computation
US9984041B2 (en) 2016-06-30 2018-05-29 International Business Machines Corporation System, method, and recording medium for mirroring matrices for batched cholesky decomposition on a graphic processing unit
JP2020502695A (ja) * 2016-12-22 2020-01-23 ライブランプ インコーポレーテッド 主成分分析による混合データ・フィンガープリンティング
EP3577781B1 (en) 2017-02-02 2025-06-18 Nokia Technologies Oy Adaptive explicit csi feedback and overhead reduction
KR101943518B1 (ko) * 2017-08-28 2019-01-30 한국과학기술연구원 소재의 전자 구조를 예측하는 방법 및 전자 장치
US10867008B2 (en) 2017-09-08 2020-12-15 Nvidia Corporation Hierarchical Jacobi methods and systems implementing a dense symmetric eigenvalue solver
US11366876B2 (en) * 2020-06-24 2022-06-21 International Business Machines Corporation Eigenvalue decomposition with stochastic optimization

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20030182518A1 (en) * 2002-03-22 2003-09-25 Fujitsu Limited Parallel processing method for inverse matrix for shared memory type scalar parallel computer
JP4325877B2 (ja) * 2004-06-03 2009-09-02 独立行政法人科学技術振興機構 行列の高速高精度特異値分解方法、プログラムおよび装置
US7895254B2 (en) * 2004-11-15 2011-02-22 Qualcomm Incorporated Eigenvalue decomposition and singular value decomposition of matrices using Jacobi rotation
US7711762B2 (en) * 2004-11-15 2010-05-04 Qualcomm Incorporated Efficient computation for eigenvalue decomposition and singular value decomposition of matrices
JP5011545B2 (ja) * 2005-12-05 2012-08-29 国立大学法人京都大学 特異値分解装置、及び特異値分解方法
JP4953239B2 (ja) * 2006-12-11 2012-06-13 インターナショナル・ビジネス・マシーンズ・コーポレーション 観測対象の異常を検出する技術

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
KATAGIRI T.: "Heiretsu Koyu Vector Keisan ni okeru Kyosei Taikaku Block-ka no Koka", (EFFECT ON ENDIAGNOAL-BLOCK METHOD FOR PARALLEL EIGENVECTOR COMPUTATIONS) INFORMATION PROCESSING SOCIETY OF JAPAN KENKYU HOKOKU, vol. 2002, no. 8, 31 August 2005 (2005-08-31), pages 43 - 48, XP003020985 *
YAMAMOTO Y. ET AL.: "Hitaisho Sanju Taikaku Gyoretsu Muke no Heiretsu Renritsu Ichiji Hoteishiki Kaiho", (A PARALLEL LINEAR EQUATION SOLVER FOR NONSYMMETRIC TRIDIAGONAL MATRICES) TRANSACTIONS OF INFORMATION PROCESSING SOCIETY OF JAPAN, vol. 43, no. SIG 9, 15 August 2001 (2001-08-15), pages 19 - 27, XP003020986 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110262479A (zh) * 2019-05-28 2019-09-20 南京天辰礼达电子科技有限公司 一种履带式拖拉机运动学估计及偏差校准方法
WO2020238011A1 (zh) * 2019-05-28 2020-12-03 南京天辰礼达电子科技有限公司 一种履带式拖拉机运动学估计及偏差校准方法

Also Published As

Publication number Publication date
US8255447B2 (en) 2012-08-28
US20100185716A1 (en) 2010-07-22
JPWO2008018188A1 (ja) 2009-12-24
JP5017666B2 (ja) 2012-09-05

Similar Documents

Publication Publication Date Title
WO2008018188A1 (fr) dispositif de décomposition de valeur propre et procédé de décomposition de valeur propre
JP5011545B2 (ja) 特異値分解装置、及び特異値分解方法
Hitzer Quaternion and Clifford Fourier Transforms
Stonyakin et al. Gradient methods for problems with inexact model of the objective
Wei Wavelets generated by using discrete singular convolutionkernels
Ceruti et al. Rank-adaptive time integration of tree tensor networks
US11971906B2 (en) Clustering apparatus, clustering method, program and data structure
Larsson et al. Dynamical pruning of the multiconfiguration time-dependent Hartree (DP-MCTDH) method: An efficient approach for multidimensional quantum dynamics
CN111033495B (zh) 用于快速相似性搜索的多尺度量化
Barvinsky et al. Heat kernel expansion for higher order minimal and nonminimal operators
Nazih et al. Using the proximal gradient and the accelerated proximal gradient as a canonical polyadic tensor decomposition algorithms in difficult situations
Alhaidari Exponentially confining potential well
Hopper Harwell subroutine library. a catalogue of subroutines (1973)
Canepa et al. Corner structure of four-dimensional general relativity in the coframe formalism
Kukelova et al. Singly-bordered block-diagonal form for minimal problem solvers
Mau Nam et al. Smoothing techniques and difference of convex functions algorithms for image reconstructions
Eriksson et al. Kähler quantization of vortex moduli
WO2018131105A1 (en) Information processing device, information processing method and storage medium
EP4616338A2 (en) Determining logical stabilizer instrument for stabilizer circuit
Dreier Hardware-oriented Krylov methods for high-performance computing
Stewart et al. Updating a generalized URV decomposition
Moreno et al. Sparse representation for white matter fiber compression and calculation of inter-fiber similarity
Santos Symmetric Central Configurations and the Inverse Problem
Thangavelu A Scalar-Valued Fourier Transform for the Heisenberg Group
Neretin On the correspondence between boson Fock space and the L2 space with respect to Poisson measure

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 07707777

Country of ref document: EP

Kind code of ref document: A1

WWE Wipo information: entry into national phase

Ref document number: 2008528725

Country of ref document: JP

WWE Wipo information: entry into national phase

Ref document number: 12376444

Country of ref document: US

NENP Non-entry into the national phase

Ref country code: DE

NENP Non-entry into the national phase

Ref country code: RU

122 Ep: pct application non-entry in european phase

Ref document number: 07707777

Country of ref document: EP

Kind code of ref document: A1