The non-linear equations necessitate the starting values to initiate the iteration process. The boundary conditions to each of the algebraic equations in the D-PT are unknown. Guessing the initial values is subjective, to say the least.
The initial (not "boundary") conditions are not "unknown": given a partition of \(M\) as  \(D+\lambda\ \Delta\), DPT uses a non-ambiguous set of initial conditions for fixed-point iteration, namely \(A^{(0)}=I\). The only "subjective" element in this procedure is the partitioning itself—but this is a feature of any perturbation theory, by definition. The need for starting values has nothing to do with the "non-linear" nature of the equations. 
Some surmised/estimated, insufficiently good starters to nonlinear equations can lead to the wrong results, irrespective of the status of the convergence issue (convergence to the wrong result is not infrequent for implicit non-linear algebraic equations).
Once again, we prove convergence to the right results under the condition \(\Vert \theta\Vert \Vert \Delta\Vert < (3-2\sqrt{2})\) with the explicit (not "surmised/estimated") initial condition  \(A^{(0)}=I\). Here too, the referee appears to be referring to issues of other iterative methods for the eigenvalue problems (e.g. the Lanczos or Arnoldi algorithms), and not to our proposed algorithm. 

Comment 4

In practice, the elements of M are not ideal entries, i.e. they enter the analysis with their intrinsic inaccuracies. Matrix elements can come from some elaborate computations with finite arithmetics (computational round-off errors) or they can be some empirical data stemming from experimental measurements (inevitably contaminated with noise–systematic, random, etc).
This is true, but since we never require or mention "ideal entries", we do not see the relevance of this comment. 
For such matrices routinely encountered in practice, the D-PT would give some non-unique eigensolutions.
This claim is unsubstantiated and incorrect. Once again, we prove convergence of DPT to a unique attracting fixed point containing a complete set of eigenvectors of \(M\) whenever \(\Vert \theta\Vert \Vert \Delta\Vert < (3-2\sqrt{2})\). Unless errors in \(M\) violate this bound, uniqueness of eigensolutions is guaranteed. 
Stated equivalently, there is no guarantee that the eigensolutions would not vary from one to another set of the eigensolutions. Some of the eigensolutions might be found in one or more subsets of all the L sets of the eigensolutions. However, there could also be some spurious eigensolutions that would differ from one to another set in {L1,L2,...,Lmax}. A key procedure is lacking in the D-PT for separating the unphysical from physical eigensolutions. 
Unlike other algorithms, DPT does not generate "spurious" or "unphysical solutions". Our method provably converges to a complete set of eigenpairs under the condition recalled above; other solutions of \(A = F(A)\) are not seen by DPT, as explained in SI A. No "key procedure" is lacking. 

Comment 5

Some of the existing generic well-known eigenvalue solvers can expediently extract millions of eigen- values from general non-Hermitean complex matrices. Any newly proposed eigenvalue solver (perturbative or nonperturbative alike) should be benchmarked on at least modestly sized non-Hermitean complex matrices (e.g. 1000 × 1000 or so) with elements corrupted by ∼ 5-20% random Gaussian-distributed noise (to mimic the experimental data) in a fixed interval. 
As already noted, Fig. S1 and 4 examine sparse matrices of size up to \(10^4\times10^4\) and \(10^6\times10^6\) respectively. The new figures address performance and numerical stability (two distinct issues) as follows: