MRRR: Robust Eigenpair Computation

• Multiple Relatively Robust Representations (MRRR) is an algorithmic framework that efficiently solves symmetric tridiagonal eigenproblems using recursively generated robust representations.
• It leverages a spectrum-shifting and task-queue paradigm to break eigenvalue clusters, achieving O(kn) computational complexity and excellent parallel scalability.
• Mixed-precision implementations within MRRR enhance numerical stability and orthogonality, offering accuracy comparable to traditional methods while maintaining speed.

Multiple Relatively Robust Representations (MRRR) is an algorithmic framework for the efficient solution of the real symmetric tridiagonal eigenproblem (STEP), a key subproblem arising in the computation of eigenvalues and eigenvectors of dense Hermitian matrices after tridiagonal reduction. MRRR is distinguished by its use of recursively generated, carefully factored matrix representations—each being a "Relatively Robust Representation" (RRR)—to extract eigenpairs accurately and efficiently. Its primary contributions are a computational complexity of $O(kn)$ for $k$ eigenpairs of an $n\times n$ matrix (with $k \approx n$ in typical full-spectrum computations), modest memory requirements, and strong parallelizability. MRRR achieves these through a spectrum-shifting and task-queue paradigm that breaks clusters of close eigenvalues via recursively constructed RRRs. The method’s accuracy, performance, and scalability have been the subject of extensive theoretical and empirical investigation, especially in recent mixed-precision variants.

1. Mathematical Foundations and Classical Algorithm

The STEP requires the solution of

$$T z_i = \lambda_i z_i, \qquad \|z_i\|_2 = 1, \quad \lambda_1 \leq \dots \leq \lambda_n$$

for a real, symmetric, tridiagonal $T \in \mathbb{R}^{n \times n}$. The MRRR approach is built upon the following core concepts:

• Relatively Robust Representation (RRR): An RRR for a subset of eigenvalues indexed by $\mathcal{I} \subset \{1, \ldots, n\}$ is a factored form of a shifted tridiagonal,

$$M = T - \sigma I = L D L^T \text{ or } U \Omega U^T \text{ or \textit{twisted/block variants}},$$

where, for any small element-wise relative perturbation $|\widetilde x_i - x_i| \leq k\,\varepsilon\,|x_i|$, the associated eigenvalues and invariant subspaces in $\mathcal{I}$ remain well-conditioned. Specifically,

$k$0

for moderate $k$1 and small $k$2.

• Algorithmic Structure: MRRR constructs a queue of "tasks," each consisting of a representation $k$3 and index set $k$4. Singleton tasks (with $k$5) yield eigenpairs via Rayleigh-quotient iteration (RQI) and twisted factorizations; clustered tasks ($k$6) trigger a shift, produce a new RRR, and subdivide the cluster by recursively building child tasks.

Algorithmic steps:

1. Preprocessing: scale and split $k$7 at tiny off-diagonals.
2. Choose an initial root shift $k$8, form $k$9, perturb.
3. Compute eigenvalue estimates for $n\times n$0 to relative accuracy set by $n\times n$1.
4. Initialize a task queue with $n\times n$2.
5. Iteratively process queue:
   • Partition $n\times n$3 into clusters according to $n\times n$4.
   • For a singleton, execute RQI and back-shift.
   • For a cluster, select shift $n\times n$5, form new RRR $n\times n$6 and enqueue subtask.

The arithmetic complexity is $n\times n$7 for $n\times n$8 eigenpairs, with $n\times n$9 auxiliary storage, and the algorithm naturally supports both depth- and breadth-first traversals of the computation tree (Petschow et al., 2013, Petschow, 2014).

2. Theoretical Properties and Error Analysis

The quality and stability of MRRR depend on several analytical guarantees:

• Residual and Orthogonality Bounds: If all RRRs meet relative robustness and conditional element growth, for unit roundoff $k \approx n$0 and deepest cluster tree depth $k \approx n$1,

$k \approx n$2

$k \approx n$3

where $k \approx n$4 sets cluster splitting sensitivity, and $k \approx n$5 is the local residual (cf. (Petschow et al., 2013, Petschow, 2014)).

• Comparison to Other Methods: Classical MRRR yields worst-case orthogonality $k \approx n$6, as opposed to $k \approx n$7 for QR and Divide–Conquer methods. In empirical tests, MRRR may lose orthogonality up to $k \approx n$8, while alternatives remain below $k \approx n$9 (Petschow et al., 2013).
• Cluster Control via $$T z_i = \lambda_i z_i, \qquad \|z_i\|_2 = 1, \quad \lambda_1 \leq \dots \leq \lambda_n$$0: Fine tuning $$T z_i = \lambda_i z_i, \qquad \|z_i\|_2 = 1, \quad \lambda_1 \leq \dots \leq \lambda_n$$1 manages the trade-off between breaking clusters (improving performance and parallelism) and orthogonality (smaller $$T z_i = \lambda_i z_i, \qquad \|z_i\|_2 = 1, \quad \lambda_1 \leq \dots \leq \lambda_n$$2 generally favored).

3. Mixed-Precision MRRR and Improved Accuracy

Mixed-precision innovations address MRRR’s classical orthogonality limitations:

• Precision Levels: Input/output is in $$T z_i = \lambda_i z_i, \qquad \|z_i\|_2 = 1, \quad \lambda_1 \leq \dots \leq \lambda_n$$3-bit ($$T z_i = \lambda_i z_i, \qquad \|z_i\|_2 = 1, \quad \lambda_1 \leq \dots \leq \lambda_n$$4), internal sensitive steps in higher $$T z_i = \lambda_i z_i, \qquad \|z_i\|_2 = 1, \quad \lambda_1 \leq \dots \leq \lambda_n$$5-bit precision ($$T z_i = \lambda_i z_i, \qquad \|z_i\|_2 = 1, \quad \lambda_1 \leq \dots \leq \lambda_n$$6). Examples: single → double or double → quadruple.
• Critical Operations: All spectrum shifts, RRR constructions, and RQI are conducted in $$T z_i = \lambda_i z_i, \qquad \|z_i\|_2 = 1, \quad \lambda_1 \leq \dots \leq \lambda_n$$7-bit to make $$T z_i = \lambda_i z_i, \qquad \|z_i\|_2 = 1, \quad \lambda_1 \leq \dots \leq \lambda_n$$8. Non-critical routines (e.g., bisection for eigenvalue refinement) may remain in $$T z_i = \lambda_i z_i, \qquad \|z_i\|_2 = 1, \quad \lambda_1 \leq \dots \leq \lambda_n$$9-bit to minimize performance cost (Petschow et al., 2013).
• Enhanced Error Bounds: With appropriate $T \in \mathbb{R}^{n \times n}$0 and typically shallow $T \in \mathbb{R}^{n \times n}$1, mixed-precision MRRR achieves

$T \in \mathbb{R}^{n \times n}$2

and

$T \in \mathbb{R}^{n \times n}$3

This matches the orthogonality of QR/Divide–Conquer methods while preserving the speed and scalability of MRRR (Petschow et al., 2013).

4. Scalability, Parallel Performance, and Implementation

MRRR is designed for high scalability in both shared-memory and distributed environments:

• Multi-core (SMP) Parallelism: The algorithm is expressed in terms of independent tasks (singleton and cluster), enabling a shared FIFO queue processed by concurrent threads. Load balancing is enhanced via R-tasks which subdivide large clusters’ eigenvalue refinement (Petschow, 2014).
• Distributed-Memory (MPI) Parallelism: Eigenpairs are distributed over $T \in \mathbb{R}^{n \times n}$4 processes (approximately $T \in \mathbb{R}^{n \times n}$5 per process). Each process operates on its set of indices, with minimal communication except when clusters span processes. Nonblocking MPI is used for necessary broadcasts and to overlap computation with communication (Petschow, 2014).
• Memory and Data Layout: MRRR requires only $T \in \mathbb{R}^{n \times n}$6 workspace, substantially less than $T \in \mathbb{R}^{n \times n}$7 needed by Divide–Conquer for explicit eigenvector backtransforms. Mixed-precision adds a negligible $T \in \mathbb{R}^{n \times n}$8 overhead (Petschow, 2014).
• Performance Highlights:
  • On platforms like Intel Xeon X7550 (32 cores), single→double mixed-precision MRRR outperforms LAPACK’s single/double-precision MRRR ({\tt SSTEMR}, {\tt DSTEMR}) and is either faster or comparable to Divide–Conquer, with substantially improved orthogonality (Petschow et al., 2013).
  • In large dense eigensolvers (including the tridiagonal reduction stage), the mixed-precision tridiagonal solution is not rate limiting but produces vectors matching the quality of QR and DC (Petschow et al., 2013).
  • Parallel implementations such as PMRRR and MR3SMP achieve 80–90% efficiency scaling to thousands of cores (Petschow, 2014).

5. Algorithmic and Practical Considerations

Table 1 summarizes recommended parameter choices and practical guidelines for mixed-precision MRRR-based eigensolvers, as given in (Petschow et al., 2013):

Component	Recommendation / Notes
$T \in \mathbb{R}^{n \times n}$9	$\mathcal{I} \subset \{1, \ldots, n\}$0; push lower for parallelism
RRR constants	$\mathcal{I} \subset \{1, \ldots, n\}$1; $\mathcal{I} \subset \{1, \ldots, n\}$2
RQI stopping	$\mathcal{I} \subset \{1, \ldots, n\}$3, set $\mathcal{I} \subset \{1, \ldots, n\}$4
Orthogonality	Expect $\mathcal{I} \subset \{1, \ldots, n\}$5 with proper $\mathcal{I} \subset \{1, \ldots, n\}$6

Key workflow steps:

1. Input: Read $\mathcal{I} \subset \{1, \ldots, n\}$7 in $\mathcal{I} \subset \{1, \ldots, n\}$8-bit; convert to $\mathcal{I} \subset \{1, \ldots, n\}$9-bit for RRR work if $$M = T - \sigma I = L D L^T \text{ or } U \Omega U^T \text{ or \textit{twisted/block variants}},$$0-bit or mixed arithmetic not hardware-accelerated.
2. Preprocessing: Scale/split in $$M = T - \sigma I = L D L^T \text{ or } U \Omega U^T \text{ or \textit{twisted/block variants}},$$1-bit.
3. RRR Core: All sensitive operations in $$M = T - \sigma I = L D L^T \text{ or } U \Omega U^T \text{ or \textit{twisted/block variants}},$$2-bit; apply perturbation of magnitude $$M = T - \sigma I = L D L^T \text{ or } U \Omega U^T \text{ or \textit{twisted/block variants}},$$3.
4. Eigenvalue Refinement: Coarse and fine bisection in $$M = T - \sigma I = L D L^T \text{ or } U \Omega U^T \text{ or \textit{twisted/block variants}},$$4-bit.
5. Main Loop: Queue tasks; for singletons, twisted-factorization RQI in $$M = T - \sigma I = L D L^T \text{ or } U \Omega U^T \text{ or \textit{twisted/block variants}},$$5-bit to $$M = T - \sigma I = L D L^T \text{ or } U \Omega U^T \text{ or \textit{twisted/block variants}},$$6; for clusters, form child RRRs.
6. Backconvert: Final $$M = T - \sigma I = L D L^T \text{ or } U \Omega U^T \text{ or \textit{twisted/block variants}},$$7 to $$M = T - \sigma I = L D L^T \text{ or } U \Omega U^T \text{ or \textit{twisted/block variants}},$$8-bit; enforce orthonormality if desired.

6. Comparative Landscape and Applicability

MRRR is contrasted with alternatives as follows:

• Divide–Conquer (DC): $$M = T - \sigma I = L D L^T \text{ or } U \Omega U^T \text{ or \textit{twisted/block variants}},$$9 worst-case arithmetic cost, $|\widetilde x_i - x_i| \leq k\,\varepsilon\,|x_i|$0 storage, high BLAS-3 intensity; sometimes fastest for well-deflated, full-spectrum cases at moderate core counts ($|\widetilde x_i - x_i| \leq k\,\varepsilon\,|x_i|$116), but loses to MRRR for larger parallel jobs (Petschow, 2014).
• QR Iteration: Also $|\widetilde x_i - x_i| \leq k\,\varepsilon\,|x_i|$2, less favorable for large $|\widetilde x_i - x_i| \leq k\,\varepsilon\,|x_i|$3 or partial spectra.
• Bisection+Inverse Iteration: $|\widetilde x_i - x_i| \leq k\,\varepsilon\,|x_i|$4 cost, poor scalability in the presence of eigenvalue clusters.
• MRRR: $|\widetilde x_i - x_i| \leq k\,\varepsilon\,|x_i|$5, minimal memory, scalable, especially advantageous when $|\widetilde x_i - x_i| \leq k\,\varepsilon\,|x_i|$6 (partial spectrum), on massive or hybrid compute resources, or when orthogonality requirements demand mixed-precision (Petschow et al., 2013, Petschow, 2014).

Typical recommended scenarios for MRRR or mixed-precision MRRR:

• Large-scale $|\widetilde x_i - x_i| \leq k\,\varepsilon\,|x_i|$7, where MRRR is fastest, and high-quality eigenvectors are needed for spectral clustering or as Rayleigh–Ritz bases.
• Multi-threaded/distributed settings facing load imbalance from eigenvalue clustering; aggressively small $|\widetilde x_i - x_i| \leq k\,\varepsilon\,|x_i|$8 disaggregates clusters and greatly enhances parallelism.
• Emerging mixed-precision hardware where $|\widetilde x_i - x_i| \leq k\,\varepsilon\,|x_i|$9-bit operations are no longer performance limiting relative to $\mathcal{I}$0-bit (Petschow et al., 2013).

7. Implementation Notes and References

• Practical implementations utilize explicit factor storage (lower-bidiagonal $\mathcal{I}$1 or $\mathcal{I}$2-representation), with recommended block sizes for dense stages (96–128 for $\mathcal{I}$3 in Elemental or ScaLAPACK) (Petschow, 2014).
• Careful attention to stable spectrum shifting (typically via DSTQDS/DQDS) is required to maintain element-wise mixed stability.
• Depth-first recursion is particularly efficient in mixed-precision, as it collapses cluster depth to $\mathcal{I}$4 (Petschow, 2014).

Principal references include:

• Dhillon & Parlett, SIAM J. Matrix Anal. Appl. 2004 (original algorithm).
• Willems & Lang, BIT 2011 (perturbation theory).
• Bientinesi et al., PMRRR (parallel implementation).
• Vömel, LAPACK Working Note 2010 (ScaLAPACK’s MRRR).
• (Petschow et al., 2013) for mixed-precision methodology and extensive accuracy/performance data.
• (Petschow, 2014) for multi-core/distributed algorithms, tuning, and comparative studies.

MRRR, particularly when augmented with mixed-precision techniques, now provides eigensolvers that are simultaneously scalable, accurate, and memory-efficient, enabling robust large-scale spectral computations in both traditional and emerging computational environments.