Final Report Summary - MATLAN (High relative accuracy matrix techniques for linear and non-linear structured eigenvalue problems and applications)
The visitor Prof. Dr Ivan Slapnicar spent a scientifically highly interesting year at the Technische Univesitaet Berlin. Most of the research was done on the first objective of the MATLAN project, which was to develop highly accurate algorithms for Hamiltonian and skew-Hamiltonian eigenvalue problems (EVPs), focussing in particular on the Hamiltonian case. The problem proved to be harder than expected and took almost the entire project time. The work continued after the visitor's return to his home institution. Since the Hamiltonian eigenvalue problem had all features of general non-symmetric EVPs, the visitor first thoroughly studied existing Jacoby-type algorithms for these problems. The summary of the obtained results is as follows:
1. the Jacobi algorithm for normal matrices converged, and actually converged quadratically. It was sufficient to apply the standard Jacobi algorithm to the Hermitian, i.e. symmetric, part or Paaredkooper algorithm on the skew-symmetric part (Loizou, 1970). The algorithm was run implicitly on a given matrix, while each transformation was to be computed from the symmetric part of the current two by two pivot submatrix. This algorithm was simpler than the known norm-reducing method, as described by Goldstine and Hurwitz, 1959.
2. for skew-symmetric matrices the algorithm of choice was the well known Paardekooper (1970) method which used four by four orthogonal transformations and reduced the starting matrix to the Murnaghan form. The method converged and was quadratically convergent.
3. for complex symmetric matrices the norm-reducing Jacobi type algorithm which used plane rotations by complex orthogonal matrices (Seaton, 1968; Anderson and Loizou, 1974) was convergent (Eberlein, 1971) and converged quadratically (Anderson and Loizou, 1973).
4. the best convergent algorithm for the general real matrix was the norm reducing method which diagonalised the given matrix, as mentioned by Eberlein, 1970. This method, which employed complex arithmetic, had better convergence properties than the methodology which used only real arithmetic (Eberlein and Boothroyd, 1968). The method converged quadratically (Ruhe, 1968) and was also more reliable than other related techniques (Veselic, 1975 and 1979; Shroff, 1989). The drawback of the method was that nearly defective eigenvalues resulted in a highly ill-conditioned eigenvector matrix, up to the inverse of the square root of machine precision.
5. for the general real matrix the Jacobi algorithm (Mehl, 2008), which was characterised by bottom-up-left-right pivoting, converged quadratically towards the complex Schur form, but did not converge in general. The algorithm utilised unitary transformations and was not norm-reducing. The adapted method for Hamiltonian matrices used unitary symplectic transformations to reduce the matrix to the Hamiltonian Schur form. This variant was also quadratically convergent, but did not need to converge. This adaptation was an improvement over the method from Byers, 1990, with respect to the pivoting scheme.
6. in case the Hamiltonian or skew-Hamiltonian matrix was additionally symmetric or skew-symmetric, the Jacobi-type methods based on unitary quaternion transformations (Mackey, 1995; Fassbender, Mackey nd Mackey, 2001) were convergent and quadratically convergent.
7. for the real Hamiltonian matrix the state of the art at the beginning of the project was the algorithm from Mehrmann, Schroeder and Watkins, 2009. The algorithm was developed from series of results (Benner, Mehrmann and Xu, 1998; Chu, Liu and Mehrmann, 2004; Watkins, 2006) and used orthogonal symplectic transformations that were applied on the implicitly squared matrix and reductions to real Schur forms of the parts of the current matrix. The symplectic transformations depended on their property of working on isotropic subspaces. However, the algorithm might not perform well in the presence of purely imaginary eigenvalues or nearly purely imaginary eigenvalues. The visitor's task was to improve this algorithm in the sense that it always computed the Hamiltonian Schur form. The visitor successfully fulfilled this task by developing three new algorithms.
The first development was the norm-reducing algorithm A1 which incorporated the Jacobi type transformations as in the method of Eberlein for general real matrices with special Hamiltonian pivoting derived by Mehl. The algorithm was convergent and converged quadratically. However, since the algorithm diagonalised the given matrix in the presence of defective or nearly defective eigenvalues, the condition of the eigenvalue matrix could grow up to the inverse of the square root of machine precision. In addition, the algorithm utilised complex transformations even in case the starting matrix was real.
The second developed algorithm A2 was the modification of the algorithm by Mehrmann, Schroeder and Watkins (MSW). The MSW algorithm carefully advanced in blocks as small as possible, trying to preserve all the necessary numerical properties of the current matrices. The new proposal was that the blocks of real, complex and purely imaginary eigenvalues were determined at the beginning from the squared starting matrix, after which the computation proceeded through these three blocks. The proposal was somewhat slower than the MSW approach, but more stable. The final block B was the block with purely imaginary eigenvalues and was then reduced by A1. Here the complex arithmetic was necessary, yet it was confined to the last block.
The final development was the algorithm A3, a modification of A2 in which the final block B was reduced to Hamiltonian Schur form by unitary symplectic transformations. The proposal was essentially the method from MSW applied to the matrix i times B which had only real eigenvalues. In addition, the subspaces that the method worked on were isotropic. The only difference to the real case was that, in this case, all transformation matrices were unitary and symplectic of the block-form (U1 U2; -conj(U2) conj(U1)). After the algorithm was finished, and if the lower triangular part of the computed Schur form was not sufficiently small, the Hamiltonian Jacobi method by Mehl could be used to reduce it quadratically. The algorithm was numerically stable since only orthogonal and unitary transformations were used.
The socioeconomic impact of the preformed research was that the algorithms would be used in subsequent projects, e.g. the European Research Council advanced grant 'Modelling, simulation and regulation of multi-physics systems' for computations with Hamiltonian matrices. The appropriate publications were pending by the time of the project completion.
1. the Jacobi algorithm for normal matrices converged, and actually converged quadratically. It was sufficient to apply the standard Jacobi algorithm to the Hermitian, i.e. symmetric, part or Paaredkooper algorithm on the skew-symmetric part (Loizou, 1970). The algorithm was run implicitly on a given matrix, while each transformation was to be computed from the symmetric part of the current two by two pivot submatrix. This algorithm was simpler than the known norm-reducing method, as described by Goldstine and Hurwitz, 1959.
2. for skew-symmetric matrices the algorithm of choice was the well known Paardekooper (1970) method which used four by four orthogonal transformations and reduced the starting matrix to the Murnaghan form. The method converged and was quadratically convergent.
3. for complex symmetric matrices the norm-reducing Jacobi type algorithm which used plane rotations by complex orthogonal matrices (Seaton, 1968; Anderson and Loizou, 1974) was convergent (Eberlein, 1971) and converged quadratically (Anderson and Loizou, 1973).
4. the best convergent algorithm for the general real matrix was the norm reducing method which diagonalised the given matrix, as mentioned by Eberlein, 1970. This method, which employed complex arithmetic, had better convergence properties than the methodology which used only real arithmetic (Eberlein and Boothroyd, 1968). The method converged quadratically (Ruhe, 1968) and was also more reliable than other related techniques (Veselic, 1975 and 1979; Shroff, 1989). The drawback of the method was that nearly defective eigenvalues resulted in a highly ill-conditioned eigenvector matrix, up to the inverse of the square root of machine precision.
5. for the general real matrix the Jacobi algorithm (Mehl, 2008), which was characterised by bottom-up-left-right pivoting, converged quadratically towards the complex Schur form, but did not converge in general. The algorithm utilised unitary transformations and was not norm-reducing. The adapted method for Hamiltonian matrices used unitary symplectic transformations to reduce the matrix to the Hamiltonian Schur form. This variant was also quadratically convergent, but did not need to converge. This adaptation was an improvement over the method from Byers, 1990, with respect to the pivoting scheme.
6. in case the Hamiltonian or skew-Hamiltonian matrix was additionally symmetric or skew-symmetric, the Jacobi-type methods based on unitary quaternion transformations (Mackey, 1995; Fassbender, Mackey nd Mackey, 2001) were convergent and quadratically convergent.
7. for the real Hamiltonian matrix the state of the art at the beginning of the project was the algorithm from Mehrmann, Schroeder and Watkins, 2009. The algorithm was developed from series of results (Benner, Mehrmann and Xu, 1998; Chu, Liu and Mehrmann, 2004; Watkins, 2006) and used orthogonal symplectic transformations that were applied on the implicitly squared matrix and reductions to real Schur forms of the parts of the current matrix. The symplectic transformations depended on their property of working on isotropic subspaces. However, the algorithm might not perform well in the presence of purely imaginary eigenvalues or nearly purely imaginary eigenvalues. The visitor's task was to improve this algorithm in the sense that it always computed the Hamiltonian Schur form. The visitor successfully fulfilled this task by developing three new algorithms.
The first development was the norm-reducing algorithm A1 which incorporated the Jacobi type transformations as in the method of Eberlein for general real matrices with special Hamiltonian pivoting derived by Mehl. The algorithm was convergent and converged quadratically. However, since the algorithm diagonalised the given matrix in the presence of defective or nearly defective eigenvalues, the condition of the eigenvalue matrix could grow up to the inverse of the square root of machine precision. In addition, the algorithm utilised complex transformations even in case the starting matrix was real.
The second developed algorithm A2 was the modification of the algorithm by Mehrmann, Schroeder and Watkins (MSW). The MSW algorithm carefully advanced in blocks as small as possible, trying to preserve all the necessary numerical properties of the current matrices. The new proposal was that the blocks of real, complex and purely imaginary eigenvalues were determined at the beginning from the squared starting matrix, after which the computation proceeded through these three blocks. The proposal was somewhat slower than the MSW approach, but more stable. The final block B was the block with purely imaginary eigenvalues and was then reduced by A1. Here the complex arithmetic was necessary, yet it was confined to the last block.
The final development was the algorithm A3, a modification of A2 in which the final block B was reduced to Hamiltonian Schur form by unitary symplectic transformations. The proposal was essentially the method from MSW applied to the matrix i times B which had only real eigenvalues. In addition, the subspaces that the method worked on were isotropic. The only difference to the real case was that, in this case, all transformation matrices were unitary and symplectic of the block-form (U1 U2; -conj(U2) conj(U1)). After the algorithm was finished, and if the lower triangular part of the computed Schur form was not sufficiently small, the Hamiltonian Jacobi method by Mehl could be used to reduce it quadratically. The algorithm was numerically stable since only orthogonal and unitary transformations were used.
The socioeconomic impact of the preformed research was that the algorithms would be used in subsequent projects, e.g. the European Research Council advanced grant 'Modelling, simulation and regulation of multi-physics systems' for computations with Hamiltonian matrices. The appropriate publications were pending by the time of the project completion.