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