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'.


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')


[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:


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:


param B:

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

type B:


param approx_type:

Type of approximation: partial_realisation or Pade.

type approx_type:


param side:

Side of the projection b or c.


Projection matrix