%% Moving the poles of a rational Krylov space
%  Mario Berljafa \and Stefan Guettel
%
%  June 2015
%
%  Tags: RAT_KRYLOV, poles


%% Rational Krylov spaces
% A rational Krylov space is a linear vector space of rational functions in a
% matrix times a vector [5]. Let $A$ be a square matrix of size $N\times N$,
% $\mathbf{b}$ an $N\times 1$ nonzero starting vector, and let
% $\xi_1,\xi_2,\ldots,\xi_m$ be a sequence of
% complex or infinite _poles_ all distinct from the eigenvalues of $A$.
% Then the rational Krylov space of order $m+1$ associated with
% $A,\mathbf{b},\xi_j$ is defined as
%
% $$\displaystyle \mathcal{Q}_{m+1}\equiv\mathcal{Q}_{m+1}(A,\mathbf{b}, q_m) = q_m(A)^{-1} \mathrm{span} \{ \mathbf{b},A\mathbf{b},\ldots,A^m \mathbf{b}\},$$
%
% where $q_m(z) = \prod_{{j=1,\xi_j\neq \infty}}^m (z - \xi_j)$ is the
% common denominator of the rational functions associated with $\mathcal{Q}_{m+1}$. 
% The rational Krylov sequence method by Ruhe [5] computes an orthonormal
% basis $V_{m+1}$ of $\mathcal{Q}_{m+1}$. The first column
% of $V_{m+1}$ can be chosen as $V_{m+1}\mathbf e_1 = \mathbf{b}/\|\mathbf{b}\|_2$. 
% The basis matrix $V_{m+1}$ satisfies a rational Arnoldi decomposition of the form
%
% $$\displaystyle  A V_{m+1} \underline{K_m} = V_{m+1} \underline{H_m}, $$
%
% where $(\underline{H_m},\underline{K_m})$ is an (unreduced) upper
% Hessenberg pencil of size $(m+1)\times m$.

%% The poles of a rational Krylov space
% Given a rational Arnoldi decomposition of the above form, it can be 
% shown [1] that the poles $\xi_j$ of the associated rational Krylov space
% are the generalized eigenvalues of the lower $m\times m$ subpencil 
% of $(\underline{H_m},\underline{K_m})$. Let us verify this at a simple
% example by first constructing a rational Krylov space associated with 
% the $m=5$ poles $-1,\infty,-\mathrm{i},0,\mathrm{i}$. 
% The matrix $A$ is of size $N=100$ and chosen as the |tridiag| matrix
% from MATLAB's |gallery|, and $\mathbf{b}$ is the first canonical unit
% vector. The |rat_krylov| command is used to compute the quantities in the
% rational Arnoldi decomposition:

N  = 100;
A  = gallery('tridiag', N);
b  = eye(N, 1);
m  = 5;
xi = [-1, inf, -1i, 0, 1i];
[V, K, H] = rat_krylov(A, b, xi);

%%
% Indeed, the rational Arnoldi decomposition is satisfied with a residual 
% norm close to machine precision:

format shorte
disp(norm(A*V*K - V*H) / norm(H))

%%
% And the chosen poles $\xi_j$ are the eigenvalues of the lower
% $m\times m$ subpencil:

disp(eig(H(2:m+1,1:m),K(2:m+1,1:m)))

%% Moving the poles explicitly
% There is a direct link between the starting vector $\mathbf{b}$ and the
% poles $\xi_j$ of a rational Krylov space $\mathcal{Q}_{m+1}$.
% A change of the poles $\xi_j$ to $\breve \xi_j$ can be interpreted as
% a change of the starting vector from $\mathbf{b}$ to $\mathbf{\breve b}$,
% and vice versa. Algorithms for moving the poles of a rational Krylov
% space are described in [1] and implemented in the functions
% |move_poles_expl| and |move_poles_impl|.
%
% RKT_BIGBREAK
%
% For example, let us move the poles of the above 
% rational Krylov space $\mathcal{Q}_{m+1}$ to 
% the points $-1,-2,\ldots,-5$:

xi_new = -1:-1:-5;
[KT, HT, QT, ZT] = move_poles_expl(K, H, xi_new);

%%
% The output of |move_poles_expl| are unitary matrices $Q$ and $Z$, 
% and transformed upper Hessenberg
% matrices $\underline{\breve K_m} = Q\underline{K_m} Z$ and  
% $\underline{\breve H_m} = Q\underline{H_m} Z$,
% so that the lower $m\times m$ part of the pencil $(\underline{\breve 
% H_m},\underline{\breve K_m})$ has as generalized eigenvalues 
% the new poles $\breve \xi_j$:

disp(eig(HT(2:m+1,1:m),KT(2:m+1,1:m)))

%%
% Defining $\breve V_{m+1} = V_{m+1} Q^*$, the transformed rational Arnoldi decomposition is 
%
% $$\displaystyle  A \breve V_{m+1} \underline{\breve K_m} = \breve V_{m+1} \underline{\breve H_m}. $$
%
% This can be verified numerically by looking at the residual norm:

VT = V*QT';
disp(norm(A*VT*KT - VT*HT) / norm(HT))

%%
% It should be noted that the function |move_poles_expl| can be used to
% move the $m$ poles to arbitrary locations, including to infinity,
% and even to the eigenvalues of
% $A$. In latter case, the transformed space $\breve V_{m+1}$ does not
% correspond to a rational Krylov space generated with starting vector 
% $\breve V_{m+1}\mathbf e_1$ and poles $\breve \xi_j$, but must be
% interpreted as a _filtered_ rational Krylov space. Indeed, the pole
% relocation problem is very similar to that of applying an implicit filter
% to the rational Krylov space [3,4]. See also [1] for more details.

%% Moving the poles implicitly
% Assume we are given a nonzero vector $\mathbf{\breve b}\in\mathcal{Q}_{m+1}$ with
% coefficient representation $\mathbf{\breve b} = V_{m+1}\mathbf{c}$,
% where $\mathbf{c}$ is a vector with $m+1$ entries. The function |move_poles_impl| can be
% used to obtain a transformed rational Arnoldi decomposition with starting
% vector $\mathbf{\breve b}$. 
%
% RKT_BIGBREAK
%
% As an example, let us take $\mathbf{c} = [0,\ldots,0,1]^T$ and hence transform
% the rational Arnoldi decomposition so that $\breve V_{m+1}\mathbf{e}_1 = \mathbf{v}_{m+1}$,
% the last basis vector in $V_{m+1}$: 

c = zeros(m+1,1); c(m+1) = 1;
[KT, HT, QT, ZT] = move_poles_impl(K, H, c);
VT = V*QT';

%%
% The poles of the rational Krylov space with the modified starting vector
% can again be read off as the generalized eigenvalues of the lower
% $m\times m$ part of $(\underline{\breve H_m},\underline{\breve K_m})$:

disp(eig(HT(2:m+1,1:m),KT(2:m+1,1:m)))

%%
% This implicit pole relocation procedure is key element of the RKFIT
% algorithm described in [1,2].

%% Some fun with moving poles
% To conclude this example, let us consider a $10\times 10$ random 
% matrix $A$, a random vector $\mathbf{b}$, and the corresponding 
% $6$-dimensional rational Krylov space with poles at $-2,-1,0,1,2$:

A  = (randn(10) + 1i*randn(10))*.5;
b  = randn(10,1) + 1i*randn(10,1);
m  = 5;
xi = -2:2; 
[V, K, H] = rat_krylov(A, b, xi);

%%
% Here are the eigenvalues of $A$:

figure
plot(eig(A),'ko','MarkerFaceColor','y')
axis([-2.5,2.5,-2.5,2.5]), grid on, hold on

%%
% We now consider a $t$-dependent coefficient vector $\mathbf{c}(t)$ such
% that $V_{m+1}\mathbf{c}(t)$ is continuously "morphed" from $\mathbf{v}_1$ 
% to $\mathbf{v}_2$. The poles of the
% rational Krylov space with the transformed starting vector $V_{m+1}\mathbf{c}(t)$
% are then plotted as a function of $t$.

for t = linspace(1,2,51),
  c = zeros(m+1,1); 
  c(floor(t))   = cos(pi*(t-floor(t))/2);
  c(floor(t)+1) = sin(pi*(t-floor(t))/2);
  [KT, HT, QT] = move_poles_impl(K, H, c);% transformed pencil
  xi_new = sort(eig(HT(2:m+1,1:m),KT(2:m+1,1:m))); % new poles
  plot(real(xi_new), imag(xi_new), 'b+')
end

%%
% As one can see, only one of the five poles starts moving away from $-2$, with the
% remaining four poles staying at their positions. This is because "morphing"
% the starting vector from $\mathbf{v}_1$ to $\mathbf{v}_2$ only affects a
% two-dimensional subspace of $\mathcal{Q}_{m+1}$ which includes the vector
% $\mathbf{b}$ and is itself a rational Krylov space, and this
% space is parameterized by one pole only. 
%
% RKT_BIGBREAK
%
% As we now continue morphing from $\mathbf{v}_2$ to $\mathbf{v}_3$,
% another pole starts moving:

for t = linspace(2,3,51),
  c = zeros(m+1,1); 
  c(floor(t))   = cos(pi*(t-floor(t))/2);
  c(floor(t)+1) = sin(pi*(t-floor(t))/2);
  [KT, HT, QT, ZT] = move_poles_impl(K, H, c);
  xi_new = sort(eig(HT(2:m+1,1:m),KT(2:m+1,1:m)));
  plot(xi_new, 'r+')
end

%%
% Morphing from $\mathbf{v}_3$ to $\mathbf{v}_4$, then to $\mathbf{v}_5$, and
% finally to $\mathbf{v}_6$ will eventually affect all five poles of the rational
% Krylov space:

for t = linspace(3, 5.99, 150)
  c = zeros(m+1,1); 
  c(floor(t))   = cos(pi*(t-floor(t))/2);
  c(floor(t)+1) = sin(pi*(t-floor(t))/2);
  [KT, HT, QT, ZT] = move_poles_impl(K, H, c);
  xi_new = sort(eig(HT(2:m+1, 1:m), KT(2:m+1, 1:m)));
  switch floor(t)
   case 3, plot(xi_new', 'g+')
   case 4, plot(xi_new', 'm+')
   case 5, plot(xi_new', 'c+')
  end
end

%% References
%
% [1] M. Berljafa and S. Guettel. _Generalized rational Krylov
%     decompositions with an application to rational
%     approximation,_ SIAM J. Matrix Anal. Appl., 36(2):894--916,
%     2015.
%
% RKT_BIGBREAK
%
% [2] M. Berljafa and S. Guettel.
%     _The RKFIT algorithm for nonlinear rational approximation,_ 
%     SIAM J. Sci. Comput., 39(5):A2049--A2071, 2017.
%
% RKT_BIGBREAK
%
% [3] G. De Samblanx and A. Bultheel. _Using implicitly filtered RKS for 
%     generalised eigenvalue problems_, J. Comput. Appl. Math., 
%     107(2):195--218, 1999.
%
% RKT_BIGBREAK
%
% [4] G. De Samblanx, K. Meerbergen, and A. Bultheel. _The implicit 
%     application of a rational filter in the RKS method_, 
%     BIT, 37(4):925--947, 1997.
%
% RKT_BIGBREAK
%
% [5] A. Ruhe. _Rational Krylov: A practical algorithm for large sparse
%     nonsymmetric matrix pencils_, SIAM J. Sci. Comput., 19(5):1535--1551,
%     1998.
