This was already inaccurate in the original submission (Fig. S1 and former Fig. 4 considered sparse random matrices of size up to \(10^4\times10^4\) and \(10^6\times10^6\) respectively). The present version includes further timings and accuracy checks for dense, complex matrices without any special structure (other than being near-diagonal) of size up to \(10^5\times 10^5\), see comment 5 below. 
No complex matrices are exemplified in the illustrations of the D-PT.
The new tests of numerical accuracy of DPT in Fig. S4 are with non-Hermitian complex matrices. Generally speaking, moving from real to complex matrices is completely innocent, as DPT relies only on matrix multiplication.  
Stability of the eigensolutions in the D-PT critically depends on symmetry properties of M. This could impact detrimentally on stability of eigensolutions as incorporating good initial conditions becomes a notable problem.
This comment is unclear to us, and appears to conflate the following distinct statements, neither of which is specifically about DPT: (i) unlike general matrices, Hermitian matrices are always well-conditioned and are subject to strong stability results (realness of eigenvalues, Weyl's theorem on perturbations of eigenvalues) (ii) choosing good initial conditions is important for certain iterative eigenvalue algorithms such as Lanczos or Arnoldi. We answer these distinct comments as follows:
Besides, there is no particular intrinsic dependence of the stability on any symmetry of the matrix. For example, the condition \(\Vert \theta\Vert \Vert \Delta\Vert < (3-2\sqrt{2})\), which guarantees the convergence, does not rely on any symmetry considerations.

Comment 2

The Rayleigh-Schrödinger perturbation theory (RS-PT) fails for degenerate spectra. Nevertheless, the perturbation formalism of quantum mechanics has overcome this drawback by resorting to the brilliantly Padé-approximant-based resummation of the Brillouin-Wigner perturbation theory (BW-PT) for operators and/or matrices with degenerate spectra. 

Comment 3

The D-PT is not practical. The reason is in solving not only one, but many non-linear algebraic equations, say L in total (L = L1 , L2 , ..., Lmax ). Each of these L’s gives a set of eigensolutions.
As explained in the main text and SI, D-PT solves as many algebraic equations as eigenpairs are requested. For a complete set of eigenpairs, that is \(N\) algebraic equations in \(\mathbb{C}^N\). What makes this "not practical"? This is as it must be. Note that although each equation has a form \(z = F_n(z)\) and, indeed, generically has \(N\) different solutions, D-PT, when converges, nevertheless, picks only one of them for each equation, providing a single solution \(A\).