%% Barycentric interpolation via the AAA algorithm % Steven Elsworth \and Stefan Guettel % % November 2017 % % Tags: AAA, RKFUN, NEPs %% Introduction % Since version 2.7 the RKToolbox provides two new utility functions % |util_bary2rkfun| and |util_aaa| for working with rational functions % in barycentric representation % % $$% r(z) = \frac{ \displaystyle \sum_{j=0}^m \frac{w_j f_j}{z - % z_j}}{\displaystyle \sum_{j=0}^m \frac{w_j}{z - z_j}}. %$$ % % The function |util_bary2rkfun| takes as inputs the interpolation % points (vector |z|), function values (vector |f|), and barycentric % weights (vector |w|) of the interpolant. It then outputs the same % rational function converted into RKFUN format. % See [2] for details on this conversion. %% % The |util_aaa| is a wrapper for the AAA algorithm developed in [3], % which appears to be a very robust method for scalar rational % approximation on arbitrary discrete point sets. The only differences of % |util_aaa| compared to the original AAA implementation in [3] are that % % * the output is returned in RKFUN or RKFUNM format, using the % conversion described in [2], % * it also works for a matrix-valued function F by using a scalar % surrogate function $f(z) = u^T F(z) v$ for the AAA sampling. %% A simple scalar example % In an example taken from [3, page 27] the AAA algorithm is used % to compute a rational approximant $r(z)$ to the Riemann zeta function % $\zeta(z)$ on the interval $[4 - 40i, 4 + 40i]$. % Here we do the same, using the |util_aaa| wrapper: zeta = @(z) sum(bsxfun(@power,(1e5:-1:1)',-z)); Z = linspace(4-40i, 4+40i); rat = util_aaa(zeta,Z); %% % The computed rational approximant |rat| is of degree 29. As it is % an object of RKFUN type we can, for example, evaluate it for matrix % arguments. Following [2], we compute $r(A)\mathbf{b}$ for a shifted % skew-symmetric matrix $A$ having eigenvalues in $[4 - 40i; 4 + 40i]$ % and a vector $\mathbf{b}$. This evaluation uses the efficient % rerunning algorithm described in [1, Section 5.1] and requires % no diagonalization of $A$. We then compute and display the relative % error of the approximant: A = 10*gallery('tridiag',10); S = 4*speye(10); A = [ S , A ; -A , S ]; b = ones(20,1); f = rat(A, b); % approximates zeta(A)*b [V,D] = eig(full(A)); ex = V*(zeta(diag(D).').'.*(V\b)); norm(ex - f)/norm(ex) %% Solving a nonlinear eigenproblem % While the AAA algorithm is designed for scalar-valued functions $f$, % the computed interpolants can easily be modified to interpolate % matrix-valued functions $F$. % The |util_aaa| implements a surrogate % approximation approach where the matrix-valued function $F$ is first reduced % to $f(z) = u^T F(z) v$ with random unit vectors and then sampled via AAA. % The scalar rational interpolant $r$ is then recast as a matrix-valued % interpolant by replacing the function values $f_j$ in the barycentric % formula with matrices $F(z_j)$. % % We demonstrate this procedure with the help of a simple nonlinear eigenvalue % problem (NEP) $F(z) = A - z^{1/2} I$, with the matrix $A$ defined above. % We choose a disk of radius 50 centered at $10 + 50i$ as the target set, % inside of which there are three eigenvalues of $F$: F = @(z) A - sqrt(z)*speye(20); evs = eig(full(A)).^2; % exact eigenvalues of F figure(1), plot(evs,'ro'); hold on; shg Z = 10 + 50i + 50*exp(1i*linspace(0,2*pi,100)); plot(Z,'k:'); axis equal; axis([-100,60,-60,100]) legend('exact evs','target set','Location','SouthWest') %% % The boundary circle is discretized by 100 equispaced points, % providing the set $Z$ of candidate points. We can now % run |util_aaa| to sample this function, say with an error tolerance % of $10^{-12}$. ratm = util_aaa(F,Z,1e-12) %% % The output is an object of class RKFUNM, a $20\times 20$ matrix-valued % rational function of degree $14$. Note that the error tolerance we % used when calling |util_aaa| corresponds to the approximation error % for the scalar surrogate problem, not necessarily for the matrix-valued NEP. % Anyway, let's plot the poles of |ratm|, which appear to align on a % (nonstandard) branch cut for the square root: plot(poles(ratm),'bs') legend('exact evs','target set','poles',... 'Location','SouthWest') %% % The solution of the rational NEP is now straightforward: We % simply need to linearize the RKFUNM and compute the eigenvalues of the % resulting matrix pencil: AB = linearize(ratm); [A,B] = AB.get_matrices(); evs = eig(full(A),full(B)); plot(evs,'m+') legend('exact evs','target set','poles','eigenvalues',... 'Location','SouthWest') %% % Note how the eigenvalues of the linearization are very good % approximations to the eigenvalues of $F$ inside the target set. % Of course, solving the linearized problem using |eig| as above is only % viable for small problems. However, this is not a problem as the pencil % structure |AB| can also be used as an input to the |rat_krylov| method. % The following lines compute Ritz approximations to the eigenvalues of the % linearization from a rational Krylov space of dimension $15$, having all % shifts identically placed at the center of the target disk: shifts = repmat(10 + 50i, 1, 15); [m,n] = type(ratm); dimlin = m*size(ratm,1); % dimension of linearization v = randn(dimlin, 1); % starting vector of Krylov space shifts = repmat(10 + 50i, 1, 14); [V, K, H] = rat_krylov(AB, v, shifts); [X, D] = eig(H(1:end-1, :), K(1:end-1, :)); ritzval = diag(D); plot(ritzval, 'g*') legend('exact evs','target set','poles','eigenvalues',... 'Ritz values','Location','SouthWest') %% % As one can see, the two Ritz values close to the center of the disk are % already quite close to the desired eigenvalues. The accuracy can be % increased by computing Ritz values of order higher than $15$. %% References % % [1] M. Berljafa and S. Guettel. % _The RKFIT algorithm for nonlinear rational approximation,_ % SIAM J. Sci. Comput., 39(5):A2049--A2071, 2017. % % RKT_BIGBREAK % % [2] S. Elsworth and S. Guettel. _Conversions between barycentric, RKFUN, % and Newton representations of rational interpolants,_ % MIMS Eprint 2017.40 (), % to appear in Linear Algebra Appl., 2019. % % RKT_BIGBREAK % % [3] Y. Nakatsukasa, O. Sete, and L. N. Trefethen. % _The AAA algorithm for rational approximation,_ % SIAM J. Sci. Comput. 40(3):A1494--A1522, 2018.