PeriodicSchurDecompositions

Given a series of NxN matrices A[j], j=1...p, a periodic Schur decomposition (PSD) is a factorization of the form:

Q[1]'*A[1]*Q[2] = T[1]
Q[2]'*A[2]*Q[3] = T[2]
...
Q[p]'*A[p]*Q[1] = T[p]

where the Q[j] are unitary (orthogonal) and the T[j] are upper triangular, except that one of the T[j] is quasi-triangular for real element types. It furnishes the eigenvalues and invariant subspaces of the matrix product prod(A).

The principal reason for using the PSD is that accuracy may be lost if one forms the product of the A_j before eigen-analysis. For some applications the intermediate Schur vectors are also useful.

This package provides a straightforward PSD for real and complex element types.

The basic API is as follows:

p = period_of_your_problem()
Aarg = [your_matrix(j) for j in 1:p]
pS = pschur!(Aarg, :R)
your_eigvals = pS.values

The result pS is a PeriodicSchur object. Computation of the Schur vectors is fairly expensive, so it is an option set by keyword argument; see the pschur! docstring for further details.

The :R argument above indicates that the product represented by pS is prod(Aarg) i.e., counting rightwards. In many applications it is more convenient to number the matrices leftwards (A[p]*...*A[2]*A[1]); this interpretation is available with an orientation argument :L.

Given a series of NxN matrices A[j], j=1...p, and a signature vector S where S[j] is 1 or -1, a generalized periodic Schur decomposition (GPSD) is a factorization of the formal product A[1]^(S[1]*A[2]^(S[2]*...*A[p]^(S[p]): Q[j]' * A[j] * Q[j+1] = T[j] if S[j] == 1 and Q[j+1]' * A[j] * Q[j] = T[j] if S[j] == -1.

The GPSD is an extension of the QZ decomposition used for generalized eigenvalue problems.

This package provides a GPSD for real and complex element types. The PSD is obviously a special case of the GPSD, but the implementations are separate for real eltypes.

The basic API is as follows:

p = period_of_your_problem()
Aarg = [your_complex_matrix(j) for j in 1:p]
signs = [sign_for_your_problem(j) for j in 1:p]
S = [x > 0 for x in signs] # boolean vector
gpS = pschur!(Aarg, S, :R)
your_eigvals = gpS.values

The result gpS is a GeneralizedPeriodicSchur object (see the docstring for further details; note that eigenvalues are stored in a decomposed form in this case).

Selected eigenvalues and their associated subspace can be moved to the top of a periodic Schur decomposition with ordschur! methods.

Reordering is available for both standard and generalized decompositions.

If only a few exterior eigenvalues (and corresponding Schur vectors) of a standard periodic system are needed, it may be more efficient to use the Krylov-Schur algorithm implemented as

    pps, hist = partial_pschur(Avec, nev, which; kw...)

where Avec is a vector of either matrices or linear maps. The result is a PartialPeriodicSchur object pps, with a summary hist of the iteration. pps usually includes the nev eigenvalues near the edge of the convex hull of the spectrum specified by which. The interface is derived from the ArnoldiMethod package, q.v. for additional details.

Note: partial_pschur is currently only implemented for the leftwards orientation.

This is a new package (2022) implementing complicated algorithms, so caveat emptor. Although tests to date indicate that the package is largely correct, there may be surprises; if you find any please file issues.

Little effort has gone into optimization so far, and the package takes a while to compile.

The API in v0.1 should be considered tentative; constructive suggestions for changes are welcome.

A. Bojanczyk, G. Golub, and P. Van Dooren, "The periodic Schur decomposition. Algorithms and applications," Proc. SPIE 1996.

D. Kressner, thesis and assorted articles.

Many functions in this package are translations of implementations in the SLICOT library.

Special thanks to Dr. A. Varga for making SLICOT available with a liberal license.

A few types and methods have been adapted from A. Noack's GenericLinearAlgebra package.