Exponential of a nonnormal matrix
This is an example of RKFIT being used for approximating , the action of the matrix exponential onto a vector . RKFIT is described in [1,2] and this code reproduces Example 3 in .
This example demonstrates that RKFIT can sometimes find sensible approximants even when is a nonnormal and all initial poles are chosen at infinity. However, we also demonstrate that the convergence can be sensitive to the intial guess, as is not surprising with nonlinear iterations. We show how real (or complex conjugate) poles can be enforced in RKFIT.
We first define the matrix , the vector , and the matrix corresponding to . Our aim is then to find a rational function of type such that is small.
N = 100; A = -5*gallery('grcar', N, 3); b = ones(N, 1); f = @(x) exp(x); fm = @(X) expm(full(X)); F = fm(A); exact = F*b; [fov, evs] = util_fovals(full(A), 100);
In order to run RKFIT we only need to specify the initial poles of . In this first test let's choose all 16 initial poles equal to infinity.
m = 16; init_xi = inf(1, m);
As all quantities , , and , as well as the initial poles are real or infinite, it is recommended to use the 'real' option. RKFIT will then attempt to produce a rational approximant with perfectly complex conjugate (or even real) poles. We set the tolerance for the relative misfit to :
maxit = 10; tol = 1e-12; xi = init_xi; [xi,ratfun,misfit,out] = rkfit(F,A,b,xi,maxit,tol,'real'); xi_rkfit = xi;
All computed poles appear in perfectly complex conjugate pairs.
xi_rkfit = Columns 1 through 2 5.0885e+00 - 2.9517e+01i 5.0885e+00 + 2.9517e+01i Columns 3 through 4 9.3537e+00 - 2.5243e+01i 9.3537e+00 + 2.5243e+01i Columns 5 through 6 1.2161e+01 - 2.1323e+01i 1.2161e+01 + 2.1323e+01i Columns 7 through 8 1.4113e+01 - 1.7513e+01i 1.4113e+01 + 1.7513e+01i Columns 9 through 10 1.5456e+01 - 1.3717e+01i 1.5456e+01 + 1.3717e+01i Columns 11 through 12 1.6329e+01 - 9.8811e+00i 1.6329e+01 + 9.8811e+00i Columns 13 through 14 1.6836e+01 - 5.9754e+00i 1.6836e+01 + 5.9754e+00i Columns 15 through 16 1.7062e+01 - 2.0015e+00i 1.7062e+01 + 2.0015e+00i
Here is a convergence history of RKFIT, showing the relative misfit defined as at each iteration. It turns out that only two iterations were required.
figure(2) semilogy(misfit, 'ro--') xlabel('iteration') title('relative 2-norm error')
The second output ratfun is an object that can be used to evaluate the computed rational approximant. This evaluation is implemented in two ways. The first option is to evaluate a matrix function by calling ratfun(A,b) with two input arguments. For example, here we are calculating the absolute misfit:
disp(norm(F*b - ratfun(A,b)))
Alternatively, we can evaluate pointwise by giving only one input argument. Let's plot the modulus of the scalar error function over a region in the complex plane:
[X,Y] = meshgrid(linspace(-18,18,500),linspace(-30,30,500)); Z = X + 1i*Y; E = f(Z) - ratfun(Z); figure(1) contourf(X,Y,log10(abs(E)),linspace(-16,8,25)) hold on plot(evs,'r.', 'MarkerSize', 6) plot(fov, 'm-') plot(xi_rkfit, 'gx') xlabel('real(z)'); ylabel('imag(z)') title('abs(exp(z) - r(z))') legend('error', 'evs', 'fov', 'poles', ... 'Location', 'NorthWest') colorbar
Choosing all initial poles equal to infinite seems to work fine. Let us try a finite initial guess, e.g., choosing all poles at :
init_xi = zeros(1, m);
Again, RKFIT requires only iterations:
[xi,ratfun,misfit,out] = rkfit(F,A,b,init_xi,maxit,tol,'real'); figure(2), hold on semilogy(misfit,'rs--') legend('RKFIT (init poles = inf)','RKFIT (init poles = 0)')
Finally, let's change the initial poles to . This turns out to be an unlucky initial guess and RKFIT fails to find a minimizer within iterations. Note that the matrix is highly nonnormal and RKFIT is a nonlinear iteration, which will probably make the convergence analysis of this example very involved.
init_xi = -10*ones(1, m); [xi,ratfun,misfit,out] = rkfit(F,A,b,init_xi,maxit,tol,'real'); semilogy(misfit,'r*--') legend('RKFIT (init poles = inf)','RKFIT (init poles = 0)',... 'RKFIT (init poles = -10)')
That's it. The following creates a thumbnail.
 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.
 M. Berljafa and S. Güttel. The RKFIT algorithm for nonlinear rational approximation, SIAM J. Sci. Comput., 39(5):A2049--A2071, 2017.