# RKToolbox Guide

Mario Berljafa, Steven Elsworth, and Stefan Güttel (The University of Manchester, UK)

## Overview

Thank you for your interest in the Rational Krylov Toolbox (RKToolbox). The RKToolbox is a collection of scientific computing tools based on rational Krylov techniques. The development started in 2013 and the current version 2.6 provides

- an implementation of Ruhe's (block) rational Krylov sequence method [6, 7], allowing to control various options, including user-defined inner products, exploitation of complex-conjugate shifts, orthogonalization, rerunning [4], and parallelism [5],
- algorithms for the implicit and explicit relocation of the poles of a rational Krylov space [3],
- a collection of utility functions, e.g., for solving nonlinear eigenvalue problems,
- an implementation of RKFIT [3, 4], a robust algorithm for rational L2 approximation, including automated degree reduction, and
- the RKFUN class [4] allowing for numerical computations with rational functions, including support for MATLAB Variable Precision Arithmetic and the Advanpix Multiple Precision toolbox [1].

This guide explains the main functionalities of the toolbox. To run the embedded MATLAB codes the RKToolbox needs to be in MATLAB's search path. For details about the installation we refer to the **Download** section on http://rktoolbox.org/.

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 method by Ruhe [6, 7] 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 [6, 7]. Other applications include matrix function approximation and rational quadrature, model order reduction, matrix equations, nonlinear eigenproblems, 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.7969e-15

As some of the poles in this example are complex, the matrices , , and are complex, too:

disp([isreal(V), isreal(K), isreal(H)])

0 0 0

However, the poles can be reordered so that complex-conjugate pairs appear next to each other 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 [6].

% 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) disp([isreal(V), isreal(K), isreal(H)])

resnorm = 5.7460e-15 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.
- One can choose between CGS and MGS with or without reorthogonalization.
- Iterative refinement for the linear system solves is supported.

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 [3] 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 [3]. 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 = 8.0477e-15 moved_poles = -1.0000e+00 + 9.3728e-17i -5.0000e-01 + 6.2977e-16i -3.3333e-01 + 2.3253e-17i -2.5000e-01 + 3.5200e-16i -2.0000e-01 + 2.7114e-16i

## Rational Krylov fitting (RKFIT)

**Relevant function:** `rkfit`

RKFIT [3, 4] is an iterative Krylov-based algorithm for nonlinear rational approximation. Given two families of matrices and , an block of vectors , and an 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 provided 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 on a simple example how to use the `rkfit` function. 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 for at most iterations 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'); disp(misfit)

7.8082e-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,`ratfun(A,V)`evaluates as a matrix function times a matrix, e.g., setting as the identity matrix will return the full matrix function , or`ratfun(z)`evaluates as a scalar function in the complex plane.

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.', 'markers', 15) 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 [4, Section 4] 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 [4, 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 an `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 `methods rkfun`. Here we provide a complete list with brief descriptions.

basis - Orthonormal rational basis functions of an rkfun. coeffs - Expansion coefficients of an rkfun. contfrac - Convert an rkfun into continued fraction form. diff - Differentiate an rkfun. disp - Display information about an rkfun. double - Convert an rkfun into double precision (undo vpa or mp). ezplot - Easy-to-use function plotter. feval - Evaluate an rkfun at scalar or matrix arguments. hess - Convert an rkfun pencil to (strict) upper-Hessenberg form. inv - Invert an rkfun corresponding to a Moebius transform. isreal - Returns true if an rkfun is real-valued. minus - Scalar subtraction. mp - Convert an rkfun into Advanpix Multiple Precision format. mrdivide - Scalar division. mtimes - Scalar multiplication. plus - Scalar addition. poles - Return the poles of an rkfun. poly - Convert an rkfun into a quotient of two polynomials. power - Integer exponentiation of an rkfun. rdivide - Division of two rkfuns. residue - Convert an rkfun into partial fraction form. rkfun - The rkfun constructor. roots - Compute the roots of an rkfun. size - Returns the size of an rkfun. subsref - Evaluate an rkfun (calls feval). times - Multiplication of two rkfuns. type - Return the type (m+k,m) of an rkfun. uminus - Unary minus. uplus - Unary plus. vpa - Convert rkfun into MATLAB's 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 expansion of `ratfun` using multiple precision arithmetic. For more details on each of the methods, type `help rkfun.<method name>`.

The RKFUN gallery provides some predefined rational functions that may be useful. A list of the options can be accessed as follows:

```
help rkfun.gallery
```

GALLERY Collection of rational functions. obj = rkfun.gallery(funname, param1, param2, ...) takes funname, a case-insensitive string that is the name of a rational function family, and the family's input parameters. See the listing below for available function families. constant Constant function of value param1. cheby Chebyshev polynomial (first kind) of degree param1. cayley Cayley transformation (1-z)/(1+z). moebius Moebius transformation (az+b)/(cz+d) with param1 = [a,b,c,d]. sqrt Zolotarev sqrt approximation of degree param1 on the positive interval [1,param2]. invsqrt Zolotarev invsqrt approximation of degree param1 on the positive interval [1,param2]. sqrt0h balanced Remez approximation to sqrt(x+(h*x/2)^2) of degree param3 on [param1,param2], where param1 <= 0 <= param2 and h = param4. sqrt2h balanced Zolotarev approximation to sqrt(x+(hx/2)^2) of degree param5 on [param1,param2]U[param3,param4], param1 < param2 < 0 < param3 < param4, h = param6. invsqrt2h balanced Zolotarev approximation to 1/sqrt(x+(hx/2)^2) of degree param5 on [param1,param2]U[param3,param4], param1 < param2 < 0 < param3 < param4, h = param6. sign Zolotarev sign approximation of degree 2*param1 on the union of [1,param2] and [-param2,-1]. step Unit step function approximation for [-1,1] of degree 2*param1 with steepness param2.

Another way to create an RKFUN is to make use of MATLAB's symbolic engine. For example, `r = rkfun('(x+1)*(x-2)/(x-3)^2')` will create a rational function as expected. Alternatively, one can specify a rational function by its roots and poles (and an optional scaling factor) using the `rkfun.nodes2rkfun` function. For example, `r = rkfun.nodes2rkfun([-1,2],[3,3])` will create the same rational function as above.

## References

[1] Advanpix LLC., *Multiprecision Computing Toolbox for MATLAB,* version 4.3.3.12213, Tokyo, Japan, 2017. http://www.advanpix.com/.

[2] M. Berljafa and S. Güttel. *A Rational Krylov Toolbox for MATLAB,* MIMS EPrint 2014.56 (http://eprints.ma.man.ac.uk/2390/), Manchester Institute for Mathematical Sciences, The University of Manchester, UK, 2014.

[3] M. Berljafa and S. Güttel. *Generalized rational Krylov decompositions with an application to rational approximation,* SIAM J. Matrix Anal. Appl., 36(2):894--916, 2015.

[4] M. Berljafa and S. Güttel. *The RKFIT algorithm for nonlinear rational approximation,* MIMS EPrint 2015.38 (http://eprints.ma.man.ac.uk/2530/), Manchester Institute for Mathematical Sciences, The University of Manchester, UK, 2015.

[5] M. Berljafa and S. Güttel, *Parallelization of the rational Arnoldi algorithm,* MIMS EPrint 2016.32 (http://eprints.ma.man.ac.uk/2503/), Manchester Institute for Mathematical Sciences, The University of Manchester, UK, 2016. Accepted for publication in SIAM J. Sci. Comput., 2017.

[6] A. Ruhe. *Rational Krylov: A practical algorithm for large sparse nonsymmetric matrix pencils,* SIAM J. Sci. Comput., 19(5):1535--1551, 1998.

[7] A. Ruhe. *The rational Krylov algorithm for nonsymmetric eigenvalue problems. III: Complex shifts for real matrices,* BIT, 34:165--176, 1994.