construct_krylov

Contructs a Krylov subspace in an iterative manner following the methods of Gugercin [1].

The construction of the Krylov space is focused on Pade and partial realisation cases for the purposes of model reduction. I.e. the partial realisation form of the Krylov space is used if approx_type = 'partial_realisation'

\[\text{range}(\textbf{V}) = \mathcal{K}_r(\mathbf{A}, \mathbf{b})\]

Else, it is replaced by the Pade approximation form:

\[\text{range}(\textbf{V}) = \mathcal{K}_r((\sigma\mathbf{I}_n - \mathbf{A})^{-1}, (\sigma\mathbf{I}_n - \mathbf{A})^{-1}\mathbf{b})\]

Note that no inverses are actually computed but rather a single LU decomposition is performed at the beginning of the algorithm. Forward and backward substitution is used thereinafter to calculate the required vectors.

The algorithm also builds the Krylov space for the \(\mathbf{C}^T\) matrix. It should simply replace B and side should be side = 'c'.

Examples

Partial Realisation:

>>> V = construct_krylov(r, A, B, 'partial_realisation', 'b')
>>> W = construct_krylov(r, A, C.T, 'partial_realisation', 'c')

Pade Approximation:

>>> V = construct_krylov(r, (sigma * np.eye(nx) - A), B, 'Pade', 'b')
>>> W = construct_krylov(r, (sigma * np.eye(nx) - A), C.T, 'Pade', 'c')

References

[1]. Gugercin, S. - Projection Methods for Model Reduction of Large-Scale Dynamical Systems. PhD Thesis. Rice University. 2003.

param r

Krylov space order

type r

int

param lu_A

For Pade approximations it should be the LU decomposition of \((\sigma I - \mathbf{A})\) in tuple form, as output from the scipy.linalg.lu_factor(). For partial realisations it is simply \(\mathbf{A}\).

type lu_A

np.ndarray

param B

If doing the B side it should be \(\mathbf{B}\), else \(\mathbf{C}^T\).

type B

np.ndarray

param approx_type

Type of approximation: partial_realisation or Pade.

type approx_type

str

param side

Side of the projection b or c.

returns

Projection matrix

rtype

np.ndarray