%% 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.