Scalable computation of Jordan chains

We present an algorithm to compute the Jordan chain of a nearly defective matrix with a $2\times2$ Jordan block. The algorithm is based on an inverse-iteration procedure and only needs information about the invariant subspace corresponding to the Jordan chain, making it suitable for use with large matrices arising in applications, in contrast with existing algorithms which rely on an SVD. The algorithm produces the eigenvector and Jordan vector with $O(\varepsilon)$ error, with $\varepsilon$ being the distance of the given matrix to an exactly defective matrix. As an example, we demonstrate the use of this algorithm in a problem arising from electromagnetism, in which the matrix has size $212^2\times 212^2$. An extension of this algorithm is also presented which can achieve higher order convergence [$O(\varepsilon^2)$] when the matrix derivative is known.