Toolbox Guide
Mario Berljafa and Stefan Güttel (The University of Manchester, UK)
Overview
This guide explains the main functionalities of the Rational Krylov Toolbox [1]. To run the MATLAB codes it is required to download the toolbox and have it added to your MATLAB path. For details about the download we refer to http://guettel.com/rktoolbox/.
Click here to view the PDF version of this guide.
Rational Krylov spaces
A rational Krylov space is a linear vector space of rational functions in a matrix times a vector. Let be a square matrix of size
,
an
starting vector, and let
be a sequence of complex or infinite poles all distinct from the eigenvalues of
. Then the rational Krylov space of order
associated with
is defined as
where is the common denominator of the rational functions associated with the rational Krylov space. The rational Krylov sequence method by Ruhe [5] computes an orthonormal basis
of
. The basis matrix
satisfies a rational Arnoldi decomposition of the form
where is an (unreduced) upper Hessenberg pencil of size
.
Rational Arnoldi decompositions are useful for several purposes. For example, the eigenvalues of the upper part of the pencil
can be excellent approximations to some of
's eigenvalues [5]. Other applications include matrix function approximation and rational quadrature, model order reduction, matrix equations, and rational least squares fitting (see below).
Computing rational Krylov bases
Relevant functions: rat_krylov, util_cplxpair
Let us compute ,
, and
using the rat_krylov function, and verify that the outputs satisfy the rational Arnoldi decomposition by computing the relative residual norm
. For
we take the tridiag matrix of size
from MATLAB's gallery, and
. The
poles
are, in order,
.
N = 100; % matrix size A = gallery('tridiag', N); b = eye(N, 1); % starting vector xi = [-1, inf, -1i, 0, 1i]; % m = 5 poles [V, K, H] = rat_krylov(A, b, xi); resnorm = norm(A*V*K - V*H)/norm(H) % residual check
resnorm = 4.1663e-15
As some of the poles in this example are complex, the matrices
,
, and
are complex, too:
[isreal(V), isreal(K), isreal(H)]
ans = 0 0 0
However, the poles can be reordered to appear in complex conjugate pairs using the function util_cplxpair. After reordering the poles, we can call the function rat_krylov with the 'real' option, thereby computing a real-valued rational Arnoldi decomposition [4].
% Group together poles appearing in complex-conjugate pairs. xi = util_cplxpair(xi); [V, K, H] = rat_krylov(A, b, xi, 'real'); resnorm = norm(A*V*K - V*H)/norm(H) [isreal(V), isreal(K), isreal(H)]
resnorm = 6.4057e-15 ans = 1 1 1
Our implementation rat_krylov supports many features not shown in the basic description above.
- It is possible the use matrix pencils
instead of a single matrix
. This leads to decompositions of the form
.
- Both the matrix
and the pencil
can be passed either explicitly, or implicitly by providing function handles to perform matrix-vector products and to solve shifted linear systems.
- Non-standard inner products for constructing the orthonormal bases are supported. Further, one can choose between CGS and MGS with or without reorthogonalization.
- Support for iterative refinement of linear system solutions.
For more details type help rat_krylov.
Moving poles of a rational Krylov space
Relevant functions: move_poles_expl, move_poles_impl
There is a direct link between the starting vector and the poles
of a rational Krylov space
. A change of the poles
to
can be interpreted as a change of the starting vector from
to
, and vice versa. Algorithms for moving the poles of a rational Krylov space are described in [2] and implemented in the functions move_poles_expl and move_poles_impl.
Example: Let us move the poles
and
into
,
.
N = 100;
A = gallery('tridiag', N);
b = eye(N, 1);
xi = [-1, inf, -1i, 0, 1i];
[V, K, H] = rat_krylov(A, b, xi);
xi_new = -1:-1:-5;
[KT, HT, QT, ZT] = move_poles_expl(K, H, xi_new);
The poles of a rational Krylov space are the eigenvalues of the lower part of the pencil
in a rational Arnoldi decomposition
associated with that space [2]. By transforming a rational Arnoldi decomposition we are therefore effectively moving the poles:
VT = V*QT'; resnorm = norm(A*VT*KT - VT*HT)/norm(HT) moved_poles = util_pencil_poles(HT, KT).'
resnorm = 6.8685e-15 moved_poles = -1.0000e+00 + 1.0476e-16i -5.0000e-01 - 2.1919e-16i -3.3333e-01 - 1.7310e-16i -2.5000e-01 - 2.7182e-16i -2.0000e-01 + 3.3675e-16i
Rational Krylov fitting (RKFIT)
Relevant function: rkfit
RKFIT [2, 3] is an iterative Krylov-based algorithm for nonlinear rational approximation. Given two families of matrices
and
, a
block of vectors
, and a
matrix
, the algorithm seeks a family of rational functions
of type
, all sharing a common denominator
, such that the relative misfit
is minimal. The matrices are optional, and if not provides
is assumed. The algorithm takes an initial guess for
and iteratively tries to improve it by relocating the poles of a rational Krylov space.
We now show how to use the rkfit function on a simple example. Consider again the tridiagonal matrix and the vector
from above and let
.
N = 100;
A = gallery('tridiag', N);
b = eye(N, 1);
F = sqrtm(full(A));
exact = F*b;
Now let us find a rational function of type
with
such that
is small. The function rkfit requires an input vector of
initial poles and then tries to return an improved set of poles. If we had no clue about where to place the initial poles we can easily set them all to infinity. In the following we run RKFIT at most
iterations of RKFIT and aim at relative misfit
below
. We display the error after each iteration.
[xi, ratfun, misfit] = rkfit(F, A, b, ... repmat(inf, 1, 10), ... 15, 1e-10, 'real'); misfit
misfit = 7.8110e-07 1.4769e-10 4.6371e-11
The rational function of type
approximates
to about
decimal places. A useful output of rkfit is the RKFUN object ratfun representing the rational function
. It can be used,for example, to evaluate
:
- ratfun(A,v) evaluates
as a matrix function times a vector, or
- ratfun(z) evaluates
as a scalar function in the complex plane.
For example, here is a plot of the error over the spectral interval of
(approximately
), together with the values at the eigenvalues of
:
figure ee = eig(full(A)).'; xx = sort([logspace(-4.3, 1, 500) , ee]); loglog(xx,abs(sqrt(xx) - ratfun(xx))); hold on loglog(ee,abs(sqrt(ee) - ratfun(ee)), 'r.') axis([4e-4, 8, 1e-14, 1e-3]); xlabel('x'); grid on title('| x^{1/2} - r_m(x) |','interpreter','tex')

As expected the rational function is a good approximation of the square root over
. It is, however, not a uniform approximation because we are approximately minimizing the 2-norm error on the eigenvalues of
, and moreover we are implicitly using a weight function given by the components of
in
's eigenvector basis.
Additional features of RKFIT are listed below.
- An automated degree reduction procedure is implemented; it takes place if a relative misfit below tolerance is achieved, unless deactivated.
- Nondiagonal rational approximants are supported; can be specified via an additional param structure.
- Utility functions are provided for transforming scalar data appearing in complex-conjugate pairs into real-valued data, as explained in [3, Section 3.5].
For more details type help rkfit.
Some of the capabilities of RKFUN are shown in the following section.
The RKFUN class
The rkfun class is the fundamental data type to represent and work with rational functions. It has already been described above how to evaluate rkfun object for scalar and matrix arguments by calling ratfun(z) or ratfun(A,v), respectively. There are more than 20 other methods implemented for rkfun, and a list of all these can be obtained by typing help rkfun:
basis - Orthonormal rational basis functions of a rkfun. coeffs - Expansion coefficients of an rkfun. contfrac - Convert rkfun into continued fraction form. diff - Differentiate an rkfun. disp - Display information about an rkfun. double - Convert rkfun into double precision (undo vpa or mp). ezplot - Easy-to-use function plotter. feval - Evaluate rkfun at scalar or matrix arguments. isreal - Returns true if a rkfun is real. minus - Scalar subtraction. mp - Convert rkfun into Advanpix Multiple Precision format. mrdivide - Scalar division. mtimes - Scalar multiplication. plus - Scalar addition. poles - Return the poles of an rkfun. poly - Convert rkfun into a quotient of two polynomials. residue - Convert a rkfun into partial fraction form. roots - Compute the roots of an rkfun. size - Returns the size of an rkfun. subsref - Evaluate an rkfun (calls feval). type - Return the type (m+k,m) of an rkfun. uminus - Unary minus. uplus - Unary plus. vpa - Convert rkfun into variable precision format.
The names of these methods should be self-explanatory. For example, roots(ratfun) will return the roots of a ratfun, and residue will compute the partial fraction form. Most methods support the use of MATLAB's Variable Precision Arithmetic (VPA) or the Advanpix Multiple Precision toolbox (MP). So, for example, contfrac(mp(ratfun)) will compute a continued fraction expanion of ratfun using multiple precision arithmetic. For more details on each of the methods, type help [name of method].
References
[1] M. Berljafa and S. Güttel. A Rational Krylov Toolbox for MATLAB, MIMS EPrint 2014.56 (http://eprints.ma.man.ac.uk/2199/), Manchester Institute for Mathematical Sciences, The University of Manchester, UK, 2014.
[2] M. Berljafa and S. Güttel. Generalized rational Krylov decompositions with an application to rational approximation, SIAM J. Matrix Anal. Appl., 2015. To appear. Available also as MIMS EPrint 2014.59 (http://eprints.ma.man.ac.uk/2278/).
[3] M. Berljafa and S. Güttel. The RKFIT algorithm for nonlinear rational approximation, MIMS EPrint 2015.38 (http://eprints.ma.man.ac.uk/2309/), Manchester Institute for Mathematical Sciences, The University of Manchester, UK, 2015.
[4] A. Ruhe. Rational Krylov: A practical algorithm for large sparse nonsymmetric matrix pencils, SIAM J. Sci. Comput., 19(5):1535--1551, 1998.
[5] A. Ruhe. The rational Krylov algorithm for nonsymmetric eigenvalue problems. III: Complex shifts for real matrices, BIT, 34:165--176, 1994.