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:
- (i) Although important in other contexts, the issue of condition numbers for non-Hermitian matrices is of little relevance to the physics audience of PRL. In physics, \(M\) is almost always the Hamiltonian operator of a quantum system, and is therefore Hermitian by definition. Given the space constraint of a Letter, we do not deem it appropriate to address this topic here.
- (ii) Unlike other iterative eigenvalue algorithms which leave the starting vector unspecified, our method provides a definite initial condition: the identity matrix \(A^{(0)}=I\), corresponding to approximating the eigenvalues of \(M\) by the diagonal elements of \(D\). There is no ambiguity regarding "good initial conditions" in our method other than the partitioning of \(M\) itself.
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.
- "Brilliantly Padé-approximant-based resummation of the Brillouin-Wigner perturbation theory" is not how degenerate eigenvalues are commonly dealt with in quantum mechanics. Instead, the standard procedure is to diagonalize the perturbation \(\Delta\) in the degenerate eigenspace and use the corresponding eigenvectors as starting points for RS-PT.
- While it is true that BW-PT does not require \(D\) to have distinct diagonal elements (unlike RS-PT), the comparison with D-PT is not pertinent: BW-PT is an implicit method which does not by itself provide numerical approximations for the eigenpairs of \(M\). Our method, by contrast, is explicit: it provides formulas approximating the eigenpairs of \(M\) to arbitrary order in \(\lambda\).
- Finally, the introduction recalls that Padé approximants are useful to improve the convergence of (RS- or BW-) PT—as are many other techniques. We are unsure why the referee singles out this particular method here, or what makes it more "brilliant" than other techniques in his/her eyes. However, we have added a reference to Brillouin-Wigner perturbation theory in the introduction.
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\).