Reducing Internal State in Eigenvalue-Only Divide-and-Conquer Tridiagonal Eigensolvers Ruiyi Zhan
Shaoshuai Zhang∗
[email protected] University of Electronic Science and Technology of China Chengdu, China
[email protected] University of Electronic Science and Technology of China Chengdu, China
arXiv:2605.26599v1 [cs.DC] 26 May 2026
Abstract Divide and Conquer (D&C) is a widely used algorithmic strategy for symmetric eigenvalue decomposition. Its natural parallelism makes D&C attractive on modern multicore CPUs and GPUs, but existing eigenvalue-only routines often default to QR-based methods because conventional D&C still materializes or replays large transformation matrices during the conquer phase. This paper proposes a boundary-row D&C algorithm for eigenvalue-only computation. The key observation is that the conquer phase only needs selected boundary rows/columns rather than the full accumulated eigenvector matrix. By propagating these boundary rows directly through the recursion, the proposed algorithm reduces the memory requirement from quadratic to linear space while also eliminating unnecessary matrix-vector work in the conventional lazy-replay formulation. We provide the algorithm, its time and space complexity analysis, correctness and stability arguments, optimized CPU and GPU implementations, and an evaluation against QR and D&C routines in standard numerical libraries.
1
Introduction
Symmetric eigenvalue decomposition (EVD) is a core primitive in numerical linear algebra and scientific computing. Given a real symmetric matrix 𝐴 ∈ R𝑛×𝑛 , EVD computes 𝐴 = 𝑄Λ𝑄𝑇 , where Λ stores the eigenvalues and 𝑄 stores the corresponding orthonormal eigenvectors. Many applications use this spectral information in quantum chemistry [17], physics [8, 16], machine learning, and signal processing [1, 11, 18, 23]. In a substantial fraction of these workloads, however, the application needs only the eigenvalues, or needs the eigenvalues before deciding whether eigenvectors are necessary. This eigenvalue-only setting changes the optimization target: the solver should not pay quadratic memory cost for eigenvector data that will not be returned. Dense symmetric EVD is normally reduced to a tridiagonal eigenproblem before the final spectral solve. LAPACK describes the standard pipeline as first reducing 𝐴 to a real symmetric tridiagonal matrix 𝑇 , and then solving the tridiagonal problem for eigenvalues and, optionally, eigenvectors [2, 7]. This separation makes the tridiagonal eigensolver
a critical kernel in full EVD implementations. Once the problem reaches 𝑇 , classical QR/QL routines are attractive for eigenvalue-only computation because they update only the diagonal and off-diagonal arrays and require little auxiliary storage [2, 22]. Their limitation is performance: the computation is sequential in nature and exposes much less coarsegrained parallelism than D&C. Divide-and-conquer (D&C) is a more parallel alternative for the symmetric tridiagonal eigenproblem. Cuppen’s method recursively splits 𝑇 into two smaller tridiagonal matrices plus a rank-one correction, solves the child problems, and merges them by solving a secular equation [4, 9, 10, 14]. The merge has substantial parallel work across secular roots and, when eigenvectors are requested, uses matrix operations to combine child eigenvectors into the parent eigenvectors. This structure explains why production libraries include D&C drivers and why D&C can be much faster than QR for large eigensystems [2, 19]. Unfortunately, the same structure also explains why D&C is rarely the default choice for eigenvalueonly computation, because conventional D&C stores or constructs dense eigenvector information inside the recursion, leading to 𝑂 (𝑛 2 ) space complexity. LAPACK stedc routine 1 also comments that the default eigenvalue-only routine is set to QR algorithm (steqr), due to the heavy memory cost of D&C. The memory mismatch comes from the vector required by each rank-one merge. After splitting the tridiagonal matrix, the rank-one update is expressed in the eigenbasis of the two child subproblems. The update vector is obtained from the last row of the left child eigenvector matrix and the first row of the right child eigenvector matrix. Standard D&C therefore keeps enough child eigenvector information to produce these boundary rows at later levels. Even if the final output contains only eigenvalues, the internal algorithm still carries dense transformation data through the D&C tree. Existing lazy-replay formulations reduce unnecessary eager multiplication, but they still represent the accumulated orthogonal transformations and therefore retain a quadratic memory footprint in practical implementations. This work asks whether eigenvalue-only D&C actually needs the dense child eigenvector matrices, and our answer
1 https://netlib.org/lapack//explore-html/d3/d57/group__stedc_ ∗ Corresponding author.
gaec55368cca7558e3ac13e04c1347bc27.html
Zhan and Zhang
is no. For the purpose of future D&C merges, each child subproblem only needs to expose the boundary row or column that participates in the next rank-one update. We therefore introduce a boundary-row D&C algorithm that propagates exactly this information through the recursion. Instead of constructing or replaying the full eigenvector matrix, the algorithm updates the boundary rows needed by ancestor merges. This changes the memory behavior of eigenvalueonly D&C from storing dense intermediate transformations to storing linear-size spectral and boundary data, while preserving the same secular-equation structure as conventional D&C. The resulting algorithm makes D&C practical for memoryconstrained eigenvalue-only workloads on CPUs and GPUs. It retains the parallel merge structure that motivates D&C, but removes the main reason libraries fall back to QR-like routines when only eigenvalues are requested. The idea is especially important on GPUs, where a quadratic workspace can determine whether a large tridiagonal problem fits in device memory at all. More broadly, our work shows that the eigenvalue-only variant of D&C should not be treated as a degenerate eigenvector algorithm; it has a smaller data dependency that can be exploited directly. The rest of the paper follows this data dependency from library behavior to algorithm, implementation, and measurement. We first review how existing tridiagonal solvers expose the QR-versus-D&C trade-off, then formalize the boundaryrow state that replaces dense replay, implement the same state contract on CPU and GPU paths, and finally measure the resulting tridiagonal-stage prototypes against QR/QL, internal values-only D&C, and cuSOLVER D&C baselines. This paper makes the following contributions: • We identify the boundary-row dependency that causes conventional eigenvalue-only D&C to retain dense eigenvector state. • We propose a boundary-row D&C algorithm that computes all eigenvalues of a symmetric tridiagonal matrix using linear auxiliary memory. • We prove that the propagated boundary rows produce the same secular-equation inputs as conventional D&C and analyze the resulting time and space complexity. • We present CPU and GPU implementations that expose D&C parallelism while avoiding full eigenvector storage. • We evaluate CPU and GPU tridiagonal-stage prototypes against QR, internal values-only D&C, and cuSOLVER D&C routines across problem sizes and spectral structures.
2
Background and Related Work
2.1
Symmetric EVD and Tridiagonal Solvers
Dense symmetric EVD is conventionally organized as a reduction phase followed by a tridiagonal eigensolve. For a real symmetric matrix 𝐴, orthogonal transformations reduce the problem to 𝐴 = 𝑄 𝐴𝑇 𝑄𝑇𝐴 , where 𝑇 is real symmetric and tridiagonal [2, 7, 21]. The eigenvalues of 𝑇 are also the eigenvalues of 𝐴. If 𝑇 = 𝑄𝑇 Λ𝑄𝑇𝑇 is computed, then the eigenvectors of the original dense matrix are 𝑄 = 𝑄 𝐴𝑄𝑇 . This paper focuses on the tridiagonal stage, and specifically on the case where only Λ is required. LAPACK exposes several algorithmic choices for the symmetric tridiagonal eigensolvers. The QR/QL family includes routines such as xSTERF, which computes all eigenvalues without eigenvectors using a square-root-free QR/QL variant [2, 22]. Bisection computes selected or all eigenvalues by locating roots in intervals [2, 5]. Relatively robust representations (MRRR) compute eigenpairs through carefully chosen shifted factorizations and are designed for high accuracy and low workspace [3, 6]. D&C, implemented in routines such as xSTEDC, targets all eigenvalues and eigenvectors and can be much faster than QR for large problems, but it requires substantially more workspace [2, 10]. GPU libraries expose similar choices, but their public interfaces do not provide a low-memory eigenvalue-only D&C path. NVIDIA cuSOLVER provides symmetric eigensolver routines with options for eigenvalues only and eigenpairs, and its documentation states that eigenvectors are computed by D&C when requested [15]. However, cuSOLVER is distributed as a closed-source library, so the internal eigenvalueonly path cannot be audited or modified. MAGMA provides a source-available LAPACK-style baseline [12, 20]. Its magma_dsyevd and magma_dsyevd_gpu implementations show a sharp split. After reducing the dense matrix to tridiagonal form, MagmaNoVec calls LAPACK dsterf, whereas MagmaVec calls MAGMA’s D&C tridiagonal solver magma_dstedx and then applies the accumulated transformations [13]. The workspace requirements reflect the same design choice: the eigenvalue-only path uses linear workspace for the dense reduction and QR/QL tridiagonal solve, while the eigenvector path requires quadratic workspace for the D&C eigenvector state [12, 13]. MAGMA therefore confirms the gap addressed by this paper, that D&C is available in production software, but primarily as an eigenvector-producing algorithm rather than as a low-memory eigenvalue-only algorithm.
2.2
Divide-and-Conquer for Tridiagonal EVD
D&C solves the tridiagonal eigenproblem by recursively reducing it to diagonal-plus-rank-one problems. Given an irreducible tridiagonal matrix 𝑇 , the algorithm splits it at an
Reducing Internal State in Eigenvalue-Only Divide-and-Conquer Tridiagonal Eigensolvers
off-diagonal entry and writes 𝑇1 0 𝑇 = + 𝜌𝑢𝑢𝑇 , 0 𝑇2
(1)
where 𝑇1 and 𝑇2 are smaller tridiagonal matrices and 𝑢 has nonzeros only at the split boundary [4, 10]. If the child decompositions are 𝑇𝑖 = 𝑄𝑖 Λ𝑖 𝑄𝑇𝑖 ,
𝑖 ∈ {1, 2},
(2)
then the parent merge becomes 𝐷 = Λ1 ⊕ Λ2,
𝑇
𝐷 + 𝜌𝑣𝑣 ,
±𝑄𝑇1 𝑒𝑚 𝑣= , 𝑄𝑇2 𝑒 1
(3)
Λ1 0 where Λ1 ⊕ Λ2 is , 𝑒 1 = [1, 0, ..., 0]𝑇 and 𝑒𝑚 = 0 Λ2 [0, ..., 0, 1]𝑇 . Thus, the merge depends on the last row of 𝑄 1 and the first row of 𝑄 2 . The eigenvalues of the parent are then the roots of a secular equation associated with 𝐷 + 𝜌𝑣𝑣𝑇 [9, 14]. After the secular equation is solved, conventional D&C constructs the parent eigenvector matrix by multiplying the block-diagonal child eigenvector matrix with the eigenvectors of the diagonal-plus-rank-one problem. This step is the source of both D&C’s performance advantage and its memory cost. It uses dense matrix operations and exposes parallelism, especially near the top of the recursion tree, but it also requires dense eigenvector matrices at intermediate levels. LAPACK documentation explicitly notes that D&C can be many times faster than QR for large matrices while requiring quadratic workspace [2]. This trade-off is appropriate when eigenvectors are requested, but it is wasteful when the user asks only for eigenvalues.
2.3
Eigenvalue-Only D&C and Lazy Replay
Eigenvalue-only D&C tries to avoid forming final eigenvectors, but the merge vector 𝑣 still creates an implicit dependency on child eigenvectors. A straightforward implementation can compute 𝑣 by constructing child eigenvector matrices, using their boundary rows, and discarding the rest. This defeats the purpose of an eigenvalue-only routine because it stores and manipulates data whose only role is to recover a small number of rows. Lazy replay improves on this baseline by postponing the application of transformations and replaying them only when boundary information is needed. It reduces some unnecessary eager work, but the replay state still represents products of dense orthogonal transformations through the recursion. Thus lazy replay identifies the relevant dependency and reduces work relative to eagerly forming all intermediate eigenvectors, but it does not by itself give a low-memory eigenvalue-only D&C formulation.
The next section turns this dependency question into the algorithmic invariant used throughout the paper: every merge must preserve exactly the boundary information needed by future merges, and no full eigenvector block should be part of the persistent values-only state.
3
Method
The background reduces the eigenvalue-only D&C problem to a precise data-dependency question: which part of the child eigenvector basis is required by later secular merges? The key observation behind our method is that lazy replay keeps too much information for the eigenvalue-only problem. At each merge, the parent does not require the full child basis; it only requires the two boundary rows that define the secular vector. Once the secular equation has been solved, ancestors again require only boundary rows of the new parent basis. Therefore the data dependency is closed under boundary-row propagation: boundary rows of the parent can be computed from boundary rows and local merge information from the children. This observation separates eigenvalue-only D&C from full-eigenvector D&C and creates an opportunity to reduce memory from quadratic to linear. This section turns the observation into an algorithmic invariant. It defines boundary-row state, proves that this state is sufficient to reproduce the same secular problems as conventional D&C, and derives the resulting work and storage model. 3.1
Problem Setting
The input to the tridiagonal stage is a real symmetric tridiagonal matrix 𝑇 ∈ R𝑛×𝑛 . A conventional D&C merge splits 𝑇 into two child tridiagonal matrices and a rank-one coupling, 𝑇𝐿 0 𝑇 = + 𝜌𝑢𝑢𝑇 , (4) 0 𝑇𝑅 where 𝑇𝐿 ∈ R𝑛𝐿 ×𝑛𝐿 , 𝑇𝑅 ∈ R𝑛𝑅 ×𝑛𝑅 , and 𝑢 is nonzero only at the two split-boundary coordinates. If 𝑇𝐿 = 𝑄 𝐿 Λ𝐿 𝑄𝑇𝐿 ,
𝑇𝑅 = 𝑄 𝑅 Λ𝑅 𝑄𝑇𝑅 ,
(5)
then the parent merge is reduced to a diagonal-plus-rank-one eigenproblem 𝐴𝑣 = diag(𝐷) + 𝜌𝑧𝑧𝑇 , 𝐷 = Λ 𝐿 ⊕ Λ𝑅 ,
(6)
𝑧 = concat(bhi(𝑄 𝐿 ), blo(𝑄 𝑅 )). Here blo(𝑄) and bhi(𝑄) denote the first and last rows of 𝑄. The key point is that the parent secular equation needs the child spectra and only two child boundary rows; it does not need the complete child eigenvector matrices. 3.2
Boundary-Row State
BR replaces the dense eigenvector state in eigenvalue-only D&C with a boundary-row state. For each node 𝑣 in the
Zhan and Zhang merge vector z
left child
Algorithm 1 Boundary-row D&C for eigenvalue-only tridiagonal EVD
right child
secular merge
Figure 1. Data dependency of one D&C merge. The highlighted row of the left child represents bhi(𝑄 𝐿 ), and the highlighted row of the right child represents blo(𝑄 𝑅 ). Concatenating these two rows forms the secular vector 𝑧, so the merge needs only boundary rows rather than complete child eigenvector matrices. computation values-only D&C
r0
S1
S2
...
persistent state Sl
rl
split-boundary rows needed by the current parent; duplicate and unordered requests are allowed and are handled by metadata rather than by materializing full matrices.
rl = (((r0S1)S2) ... )Sl
BR D&C
BL
[BL; BR]Sv
Input: Symmetric tridiagonal subproblem 𝑇𝑣 and requested local rows 𝜎𝑣 Output: Eigenvalues Λ𝑣 and selected rows (𝑄 𝑣 )𝜎𝑣 1: if 𝑇𝑣 is a leaf problem then 2: Solve the leaf by a small tridiagonal eigensolver. 3: Return Λ𝑣 and the requested rows of the leaf eigenvector matrix. 4: end if 5: Split 𝑇𝑣 into left and right child problems. 6: Map 𝜎 𝑣 and split-boundary requests to child row lists 𝜎 𝐿 , 𝜎𝑅 . 7: Recursively compute (Λ𝐿 , (𝑄 𝐿 )𝜎𝐿 ) and (Λ𝑅 , (𝑄 𝑅 )𝜎𝑅 ). 8: Form 𝐷 = Λ𝐿 ⊕ Λ𝑅 and 𝑧 = concat(bhi(𝑄 𝐿 ), blo(𝑄 𝑅 )). 9: Solve the local secular equations for Λ 𝑣 and the local transform 𝑆 𝑣 needed by requested rows. 10: Compute (𝑄 𝑣 )𝜎𝑣 = (𝑄 𝐿 ⊕ 𝑄 𝑅 )𝜎𝑣 𝑆 𝑣 by selected-row updates. 11: return Λ 𝑣 and (𝑄 𝑣 )𝜎𝑣
3.2.1
Bv
Algebraic Correctness.
Bv = [BL; BR]Sv
Figure 2. State and computation carried between D&C levels. The internal values-only D&C path reconstructs a requested boundary row by replaying a chain of GEMV-like transformations, 𝑟 ℓ = (((𝑟 0𝑆 1 )𝑆 2 ) · · · )𝑆 ℓ , while retaining dense or replayable state with quadratic footprint. BR instead propagates boundary rows by the local recurrence 𝐵 𝑣 = [𝐵𝐿 ; 𝐵𝑅 ]𝑆 𝑣 , so only linear-size boundary state is persistent.
Lemma 3.1 (Boundary sufficiency). For any D&C merge with child eigenvector matrices 𝑄 𝐿 and 𝑄 𝑅 , the secular vector 𝑧 in Eq. (6) is fully determined by bhi(𝑄 𝐿 ) and blo(𝑄 𝑅 ). Proof. The rank-one coupling in Eq. (4) acts only on the last coordinate of the left child and the first coordinate of the right child. After transforming to the child eigenvector bases, these two coordinate vectors become 𝑄𝑇𝐿 𝑒𝑛𝐿 and 𝑄𝑇𝑅 𝑒 1 . Their entries are precisely the last row of 𝑄 𝐿 and the first row of 𝑄 𝑅 , giving Eq. (6).
merge tree, let 𝑄 𝑣 be the eigenvector matrix that conventional D&C would construct for the corresponding subproblem. BR stores
Lemma 3.2 (Selected-row multiplication). Let 𝑄 ∈ R𝑚×𝑘 , 𝑆 ∈ R𝑘 ×ℓ , and let 𝜎 be any row list, possibly unordered and with repetitions. Then
𝐵 𝑣 = (blo(𝑄 𝑣 ), bhi(𝑄 𝑣 )).
(𝑄𝑆)𝜎 = 𝑄 𝜎 𝑆.
(7)
For an internal node, the left component required by its parent is a selected row of the parent eigenvector block; the right component is another selected row. These selected rows can be computed directly from child selected rows and the local secular-vector block, rather than by forming all rows of 𝑄 𝑣 . At the root, only eigenvalues are returned. Boundary rows are nevertheless propagated through internal nodes because they are the information needed by ancestor secular equations. In implementation, the requested row list 𝜎𝑣 contains exactly the rows demanded by ancestors and the
(8)
Proof. For every requested row position 𝑖 and output column 𝑗, ((𝑄𝑆)𝜎 )𝑖 𝑗 =
𝑘 ∑︁ 𝑡 =1
𝑄 𝜎 (𝑖 ),𝑡 𝑆𝑡 𝑗 =
𝑘 ∑︁
(𝑄 𝜎 )𝑖𝑡 𝑆𝑡 𝑗 = (𝑄 𝜎 𝑆)𝑖 𝑗 .
𝑡 =1
Thus the two matrices are equal entrywise. Theorem 3.3 (BR computes the same eigenvalues as conventional D&C in exact arithmetic). Assume BR and conventional D&C use the same split tree, deflation decisions, column ordering, column signs, and secular-equation convention. In exact
Reducing Internal State in Eigenvalue-Only Divide-and-Conquer Tridiagonal Eigensolvers
arithmetic, every internal BR merge constructs the same secular problem as conventional D&C, and the root eigenvalues returned by BR are identical to those returned by conventional D&C. Proof. The proof is by induction over the merge tree. At leaves, BR solves the same leaf eigenproblem as conventional D&C and returns the requested rows of the same local eigenvector matrix. For an internal node, the induction hypothesis gives the same child spectra and the same requested child boundary rows. By Lemma 3.1, BR constructs the same 𝑧 and therefore the same diagonal-plus-rank-one secular problem as conventional D&C. The local eigenvalues are therefore the same. By Lemma 3.2, the parent rows requested by ancestors can be computed from selected child rows without forming unrequested rows. Hence the induction hypothesis also holds for the parent boundary state. Applying the argument up to the root proves the claim. 3.2.2 Conditioned Error Propagation. BR removes storage and arithmetic associated with unrequested rows, but it still relies on the same local numerical ingredients as D&C, including secular equation solves, deflation handling, and products used to form selected rows. Our stability statement is therefore conditioned on standard local stability interfaces rather than a fresh proof of the full LAPACK-style secular solver. b𝑣 be the computed boundary state at node 𝑣. If the Let 𝐵 two child boundary-row errors satisfy c 𝐿 − bhi𝐿 bhi
≤ 𝜂𝐿 , 2
c 𝑅 − blo𝑅 blo
≤ 𝜂𝑅 , 2
(9)
then the error in the constructed secular vector satisfies ∥b 𝑧 − 𝑧 ∥ 22 =
c 𝐿 − bhi𝐿 bhi
2 2
c 𝑅 − blo𝑅 + blo
2 2
𝑇lazy (𝐾) = (𝑐 sec + 𝑐 rep + 4)𝐾 2 + 𝑂 (𝐾), where the constant term represents the remaining quadratic passes for updated-vector reconstruction and dense local secular-vector materialization. BR removes the replay term and avoids writing a dense local secular-vector block, but it does not reduce the asymptotic cost of the secular root solves. At a non-root merge, BR streams each secular vector column through at most two selected boundary rows. In the same pass-count model, 𝑇BR,nonroot (𝐾) = (𝑐 sec + 4)𝐾 2 + 𝑂 (𝐾). At the root, no parent consumes boundary rows, so BR skips selected-row propagation and computes only the final secular roots: 𝑇BR,root (𝐾) = 𝑐 sec𝐾 2 + 𝑂 (𝐾). For a balanced D&C tree with constant leaf size, ∑︁ 𝐾 2 = 2𝑛 2 + 𝑂 (𝑛), internal nodes
∑︁
≤ 𝜂𝐿2 + 𝜂𝑅2 . (10)
Thus BR does not amplify the child boundary error when forming 𝑧 beyond the norm of the concatenation. The remaining local error is the error of the secular solve and the selected-row update. Section A gives the detailed conditioned proof, including compact delta reconstruction, pointwise-toFrobenius bounds, selected-row rounding-error composition, and tree-level propagation. 3.3
We use a pass-count model to make the constant-factor difference explicit. A length-𝐾 linear pass is counted as 𝐾, and 𝐾 such passes are counted as 𝐾 2 . For a merge with active rank 𝐾, both BR and internal values-only D&C solve the same 𝐾 secular equations. We write this shared rootsolving cost as 𝑐 sec 𝐾 2 , where 𝑐 sec absorbs the secular solver’s implementation-dependent iteration count and reduction cost. Lazy replay additionally reconstructs the current merge vector from stored transformation state; we denote this cost by 𝑐 rep 𝐾 2 . It also reconstructs updated secular quantities and materializes dense local 𝐾 × 𝐾 vector state for later replay. Thus a non-root internal values-only D&C merge has the leading form
Complexity
The relevant comparison is between values-only D&C paths, not between BR and the public QR/QL routine xSTERF. The public xSTEDC values-only interface may choose the lowworkspace QR/QL path, whereas LAPACK’s internal valuesonly D&C machinery, represented by DLAED0(ICOMPQ=0), retains enough replay state for DLAEDA to reconstruct future merge vectors. BR targets this second point in the design space: it keeps the D&C merge and secular-solve structure, but changes the vector-derived state that is carried across levels.
𝐾 2 = 𝑛 2 + 𝑂 (𝑛).
(11)
non-root internal nodes
Substituting Eqs. (3.3)–(3.3) gives the leading merge-path costs 𝑇lazy (𝑛) = (2𝑐 sec + 2𝑐 rep + 8)𝑛 2 + 𝑂 (𝑛 log 𝑛), and 𝑇BR (𝑛) = (2𝑐 sec + 4)𝑛 2 + 𝑂 (𝑛 log 𝑛). Thus the saved leading work is 𝑇lazy (𝑛) − 𝑇BR (𝑛) = (2𝑐 rep + 4)𝑛 2 + 𝑂 (𝑛 log 𝑛). This is a constant-factor reduction in the D&C merge path, not a reduction of eigenvalue computation below the 𝑂 (𝑛 2 ) secular root-solving work. The asymptotic gain is in storage. Conventional lazyreplay D&C stores dense local secular-vector blocks and replay metadata, yielding quadratic real workspace and 𝑂 (𝑛 log 𝑛) integer metadata in LAPACK’s internal values-only D&C path. BR stores pole arrays, compact secular data, local metadata, and only the selected boundary rows required by ancestors. The resulting auxiliary state is 𝑂 (𝑛), excluding the
Zhan and Zhang
input tridiagonal arrays and output eigenvalues. The implementation section uses this result as a contract: preserve the standard D&C merge semantics, but make boundary rows and compact secular data the only persistent eigenvectorderived state.
4
Implementation
The GPU and CPU code paths implement the same algorithmic modification to tridiagonal D&C. Both preserve the standard split tree, deflation, sorting, secular solves, and denominator reconstruction, but replace persistent full-eigenvector or replay state with the boundary-row state from Section 3.2. They differ only in platform realization: the GPU path maps the same state contract to resident CUDA kernels, while the CPU path exposes it through a LAPACK-style prototype for controlled comparison with existing routines. 4.1
GPU Implementation
The GPU values-only path is implemented as a resident tridiagonal D&C solver. It is used only when eigenvectors are not requested; full-vector tridiagonal D&C, dense back transformation, and other EVD stages use separate paths. The solver builds the split tree on the GPU, solves leaf subproblems, extracts only the first and last eigenvector rows from each leaf, and then processes the D&C tree bottom-up. Merge tasks at the same level are independent, so each level is handled by batched CUDA kernels over all active merges. The persistent eigenvector-derived state is the boundary-row state, while split metadata, active-column maps, compact secular-root data, and temporary reduction buffers are stored in GPU workspace. At each non-root merge, the GPU path constructs the rank-one update vector from child boundary rows, 𝑧 = [bhi(𝑄 𝐿 ), blo(𝑄 𝑅 )]. The implementation then performs the standard D&C merge steps on GPU, which include sorting child eigenvalues, detecting negligible 𝑧𝑖 , detecting close poles, applying the corresponding Givens rotations, compacting the active secular problem, and recording active/deflated column mappings. The values-only distinction is the scope of these updates. Permutations and rotations that would normally be applied to full eigenvector blocks are applied only to the selected boundary rows and the associated metadata. This keeps deflation consistent with the standard algorithm while avoiding full-vector state. After deflation, the active secular problem has rank 𝐾. The root solve evaluates the eigenvalues of 𝐷 + 𝜌𝑧𝑧𝑇 with CUDA kernels that parallelize both across roots and across the pole reductions inside each root. For large 𝐾, block-reduction kernels evaluate sums such as ∑︁ 𝑧 2 𝑖
𝑖
𝑑𝑖 − 𝜆
and the associated rescaling quantities in parallel within a CUDA block. This is important near the upper levels of the merge tree, where real dense-derived tridiagonal inputs often have limited deflation and therefore large active rank. For non-root merges, the parent boundary rows are produced by streaming secular eigenvector columns through the selected rows. For an active root 𝜆 𝑗 , the secular vector has entries 𝑧𝑖 /(𝑑𝑖 − 𝜆 𝑗 ) . 𝑦 𝑗 (𝑖) = 𝑧/(𝑑 − 𝜆 𝑗 ) 2 Instead of materializing the dense 𝐾 × 𝐾 secular eigenvector block 𝑌 , the kernel directly computes 𝑅parent (:, 𝑗) = 𝑅child 𝑦 𝑗 , where 𝑅child contains at most two selected rows. Thus each column update is reduced to two streamed dot products. At the root, no parent will consume boundary rows, so the implementation switches to a root-only mode that computes the final secular roots and skips boundary-row propagation entirely. The GPU implementation also avoids a dense 𝐾 × 𝐾 denominator matrix. Each secular root is stored in a compact representation consisting of an origin pole, an offset 𝜏, and the near-pole denominator entries most sensitive to cancellation. Conceptually, 𝜆 𝑗 = 𝑑 origin + 𝜏 𝑗 ,
𝛿𝑖 = 𝑑𝑖 − 𝑑 origin − 𝜏 𝑗 .
During selected-row propagation, the kernel reconstructs each secular vector column from this compact state and uses the cached near-pole denominators where direct subtraction would lose relative accuracy. This preserves the stable denominator representation needed for secular-vector reconstruction while keeping the temporary state linear in the active merge size. Overall, the GPU realization removes the persistent dense eigenvector matrix, dense secular eigenvector blocks, and replayable transformation history from the tridiagonal D&C stage. The algorithmic state remains the standard D&C secular problems and deflation metadata, but the vector-derived state is reduced to boundary rows. The CPU realization below follows the same state discipline in a LAPACK-compatible setting, which gives a controlled experimental vehicle for isolating the tridiagonal-stage behavior evaluated in Section 5. 4.2
CPU Implementation
The CPU prototype implements the same boundary-row algorithm as a private LAPACK-style path, DSTEDC_DC_BR. It deliberately does not change the public DSTEDC(’N’) behavior, which continues to use the standard low-memory QR/QL path. The private driver preserves the usual outer structure of a tridiagonal eigensolver: it validates arguments, splits at negligible off-diagonal entries, scales unreduced blocks, calls DSTERF for small blocks, and routes larger blocks to DLAED0_BR. This gives a controlled LAPACK-compatible
Reducing Internal State in Eigenvalue-Only Divide-and-Conquer Tridiagonal Eigensolvers
baseline for evaluating boundary-row D&C without changing public routine semantics. DLAED0_BR keeps the standard D&C tree and stores two arrays, BLO and BHI, for the first and last rows of each active local eigenvector block. Leaf blocks are solved with DSTEQR(’I’), but only these two rows are copied into persistent state. DLAED7_BR performs a single merge, DLAED8_BR mirrors LAPACK’s sorting, deflation, permutation, and closepole rotation logic on selected rows, and DLAED9_BR streams secular columns through those rows. DLAED4_BR provides the same compact denominator representation used by the GPU path, so the CPU implementation avoids storing full 𝐾 × 𝐾 denominator columns while retaining near-pole entries. The CPU path exposes OpenMP parallelism without adding thread-dependent workspace. DLAED0_BR parallelizes independent merge pairs at the same tree level; each merge receives a scratch slice based on its physical interval, so same-level merges use disjoint workspace. For the root-only merge, independent secular root solves are parallelized using existing 𝑂 (𝐾) scratch slices as per-thread denominator buffers. The large-block workspace query is 16𝑁 doubleprecision entries and 7𝑁 integer entries. In comparison, LAPACK’s internal values-only D&C path DLAED0(ICOMPQ=0) requires 1 + 3𝑁 + 2𝑁 ⌈log2 𝑁 ⌉ + 3𝑁 2 real workspace and 6 + 6𝑁 + 5𝑁 ⌈log2 𝑁 ⌉ integer workspace. Thus the CPU boundary-row path reduces the leading real workspace term from 𝑂 (𝑁 2 ) to 𝑂 (𝑁 ) and removes the 𝑂 (𝑁 log 𝑁 ) integer replay metadata. The CPU implementation therefore closes the loop from the method to measurement: it keeps the same D&C merge semantics as the internal values-only D&C baseline, exposes OpenMP parallelism without thread-dependent persistent storage, and reports an explicit linear workspace query that can be checked directly in experiments.
5
Experiments
The preceding sections make three claims about BR: it should preserve the D&C merge semantics, replace persistent valuesonly D&C replay state with linear-size boundary data, and expose more parallelism than the low-memory QR/QL path when the spectrum permits useful deflation. This section tests those claims as one evidence chain. We first measure workspace and runnable scale, then compare against CPU QR/QL and internal values-only D&C baselines, and finally test the same values-only D&C design point against cuSOLVER on an H100. Dense reduction and back transformation are outside the measured kernel because BR changes the tridiagonal D&C stage rather than the surrounding dense EVD pipeline.
5.1
Experimental Setup
All CPU experiments were run on one full compute node with two Intel Xeon Gold 6348 sockets, 28 physical cores per socket, and 56 cores in total. The jobs requested 56 CPU slots and 64 GB of memory from the scheduler. The code was compiled with gfortran -O2 -fopenmp -std=legacy and linked against Intel oneAPI MKL 2025.3 through mkl_rt. We set OMP_PROC_BIND=close, OMP_PLACES=cores, disable MKL dynamic threading, and use the LP64 GNU-threaded MKL runtime. For BR-only timing runs, MKL_NUM_THREADS=1 so that the measured parallelism comes from the private BR OpenMP implementation rather than from nested MKL calls. GPU experiments were run on one NVIDIA H100 PCIe GPU with 80 GB of device memory and CUDA 13.2. The scheduler job requested one h100_pcie GPU, eight CPU cores, and 80 GB of host memory. We use our GPU EVD prototype and build the tridiagonal D&C benchmark for SM90. The GPU benchmark runs the values-only tridiagonal D&C stage with –backend=2, –vectors=0, –run-self=1, and –run-cusolver=1; it therefore measures the same D&C stage targeted by BR rather than a full dense EVD pipeline. The synthetic and structured GPU rows use two warmup iterations and five timed iterations. The CPU experiments use two fixed-seed pseudo-random tridiagonal families and two deterministic structured stress cases. The uniform family uses 𝑑𝑖 ∼ 𝑈 [−1, 1] and 𝑒𝑖 ∼ 𝑈 [0.10, 0.30]. The normal family uses 𝑑𝑖 ∼ N (0, 1) and the same positive uniform off-diagonal distribution. Both pseudo-random families use a fixed xorshift seed determined by the distribution and 𝑁 , so every matrix is exactly reproducible. The Toeplitz family uses constant entries 𝑑𝑖 = 2 and 𝑒𝑖 = 0.25. The clustered family uses nearly equal diagonal entries 𝑑𝑖 = 1 + 10−12 (𝑖 − (𝑛 + 1)/2) and small off-diagonal entries 𝑒𝑖 = 10−4 (1 + 0.1 cos(0.33𝑖)). For the reduced-dense case, the H100 first reduces a dense symmetric input to tridiagonal form and writes only the compact 𝐷/𝐸 arrays; the CPU then reads that cache and measures only the tridiagonal eigensolver stage. Timing reports the best elapsed time over repeated runs for small sizes and a single run for large sizes, matching the benchmark scripts. Accuracy is measured against DSTERF whenever a DSTERF reference was computed; for sizes beyond the reference range we report successful completion and checksums rather than claiming a reference error. We use normalized accuracy metrics, ∥ 𝜆b − 𝜆ref ∥ ∞ , max(1, ∥𝜆ref ∥ ∞ ) ∥ 𝜆b − 𝜆ref ∥ ∞ 𝑒 bwd = . max(1, ∥𝑇 ∥ ∞ ) 𝑒 fwd =
(12) (13)
The second quantity is a backward-error certificate scaled by the input tridiagonal norm. Raw absolute eigenvalue differences are omitted because they depend directly on the
Zhan and Zhang
spectrum scale. All memory numbers use the workspace query returned by the BR routine, which is 16𝑁 double entries and 7𝑁 integer entries for large blocks. 5.2
Baselines
We compare BR with two CPU baselines. The first baseline is LAPACK/MKL DSTERF, the standard eigenvalue-only QR/QL routine. It has very small storage requirements, essentially the diagonal and off-diagonal arrays, but exposes limited parallelism. The second baseline is the internal values-only D&C path, represented by the LAPACK-style DLAED0(ICOMPQ=0) machinery in MKL. This path is the closest algorithmic comparison for BR because it preserves D&C secular merges and returns eigenvalues only, but it carries replay state for merge-vector reconstruction and therefore has much larger workspace. We use the term internal values-only D&C consistently for this baseline, rather than treating it as a public solver interface. The standard LP64 measurements cover this path through 𝑁 = 16,384; selected ILP64 measurements extend the comparison to larger rows whose queried LWORK exceeds the LP64 integer range. On the GPU we compare against cuSOLVER cusolverDnXstedc with compz=N, the public cuSOLVER tridiagonal D&C interface for eigenvalues only. 5.3
Table 1. Workspace design points at 𝑁 = 65,536 on the fixedseed uniform family. QR/QL is the lowest-memory baseline; BR uses more linear workspace to expose D&C parallelism, while internal values-only D&C is quadratic. Path
Mem. at 65,536 Time at 65,536
Space
DSTERF / QR-QL 𝑂 (𝑁 ) input only BR, 56T 𝑂 (𝑁 ) D&C state Internal values-only D&C 𝑂 (𝑁 2 ) replay state
1.00 MiB 9.75 MiB 96.0 GiB
56.52 s 10.56 ms OOM
Table 2. CPU comparison against DSTERF. Larger ratios mean BR is faster. The reduced-dense row uses 𝐷/𝐸 arrays produced once by GPU dense-to-tridiagonal reduction; timings exclude that shared reduction and measure only the tridiagonal eigensolver. Input
𝑁
DSTERF
BR 1T
BR 56T
Fixed-seed pseudo-random Uniform 16,384 4.17 s 79.2× 1330.8× Uniform 65,536 56.52 s 262.2× 5351.1× Normal 16,384 3.48 s 80.8× 1331.2× Normal 65,536 48.85 s 274.9× 4378.1× Structured stress cases Toeplitz 16,384 3.40 s Toeplitz 65,536 51.94 s Clustered 16,384 3.27 s Clustered 32,768 12.76 s
1.58× 1.58× 0.84× 0.83×
3.99× 4.03× 3.00× 6.75×
Reduced dense Reduced dense 49,152 29.85 s
0.95×
8.06×
Workspace and Scale
This workspace result is a design-space change, not a claim that BR uses less memory than QR/QL. DSTERF is the lowestmemory eigenvalue-only baseline: it stores only the tridiagonal arrays and therefore remains the right reference point when memory is the only objective. BR deliberately spends a larger, but still linear, workspace budget to recover the D&C merge parallelism without carrying dense eigenvector replay state. In production-style storage, the tridiagonal input plus BR workspace is about 172𝑁 bytes: 𝐷/𝐸, WORK(16N) in double precision, and IWORK(7N) in 32-bit integers. This is about 10.75× the QR/QL input storage, but it is still 𝑂 (𝑁 ); the internal values-only D&C baseline, in contrast, has a quadratic real-workspace term. Table 1 summarizes this trade-off at a shared reference size, 𝑁 = 65,536, on the fixed-seed uniform family. At this size, QR/QL needs only about 1.00 MiB for 𝐷/𝐸, but takes 56.52 s. BR uses 9.75 MiB of queried workspace and takes 10.56 ms on 56 threads. The internal values-only D&C formula would require about 96.0 GiB of real workspace at the same 𝑁 , so we report it as OOM in the standard LP64 experiment. Table 3 separately includes ILP64 internal-D&C rows at 𝑁 = 32,768 and on the reduced-dense input to show the cost once the interface limit is removed. The scale limit for BR in the current LP64 build is not node memory but the 32-bit LWORK count: BR successfully ran 𝑁 = 134,000,000, within 0.2% of the bound ⌊(231 − 1)/16⌋.
5.4
Performance Against QR/QL
BR is much faster than DSTERF on the fixed-seed pseudorandom families, and it still gives useful threaded speedups on the harder structured spectra and on a dense-derived tridiagonal input. Table 2 reports these CPU cases together so that the reader can compare the same baseline, metric, and thread count across input families. At 𝑁 = 65,536, 56-thread BR is 5351× faster than DSTERF on uniform inputs and 4378× faster on normal inputs. On Toeplitz and clustered inputs, single-thread BR is not the right comparison point because both algorithms are close to quadratic; the threaded BR path is the useful one, giving about 4.0× on Toeplitz and 6.7× on the largest updated clustered row. The reduced-dense row connects this tridiagonal study back to a dense EVD pipeline without charging either CPU solver for dense reduction. A GPU job first reduces a 49152 × 49152 dense symmetric matrix to tridiagonal form, stores the compact 𝐷/𝐸 arrays, and discards the dense intermediate. The CPU then solves exactly those cached arrays: singlethread BR is slightly slower than DSTERF, but the parallel BR path recovers an 8.06× speedup at 56 threads with the current root-solve scheduling. This result supports the intended hybrid usage: dense reduction can run on the GPU, while the compact tridiagonal stage remains transferable and can still benefit from BR on CPU.
Reducing Internal State in Eigenvalue-Only Divide-and-Conquer Tridiagonal Eigensolvers
The speedup comes from a change in empirical scaling on the pseudo-random families. For 𝑁 ≥ 65,536, uniform BR fits 𝑁 1.040 on one thread and 𝑁 0.973 on 56 threads, while normal BR fits 𝑁 1.033 and 𝑁 0.950 , respectively. Over the referenced range 𝑁 = 4096 to 65,536, DSTERF fits approximately 𝑁 1.916 on uniform inputs and 𝑁 1.907 on normal inputs. This should be read as an empirical property of the measured spectra, not as an asymptotic proof that BR makes all tridiagonal eigenvalue problems linear. 5.5
Comparison with Internal D&C
The closest CPU D&C baseline is MKL’s internal values-only path that calls DLAED0(ICOMPQ=0). This path has the same eigenvalue-only output as BR and uses the same secularmerge family, but it keeps the replay state needed to reconstruct merge vectors. Table 3 reports this comparison on the fixed-seed pseudo-random families, where BR has near-linear empirical scaling in Table 2. With one thread, BR is 3.1–3.8× faster than the internal values-only D&C path through 𝑁 = 16,384, while using 2.44 MiB rather than 6.0 GiB of workspace. At 𝑁 = 32,768, ILP64 internal-D&C runs require about 24.0 GiB, whereas BR still uses 4.88 MiB. The multithreaded rows show that the remaining CPU gap was an implementation bottleneck rather than a necessary cost of the BR formulation. The leaf eigenproblems are independent, but an earlier prototype initialized them serially with DSTEQR(’I’); parallelizing those leaf solves removes the dominant cost on pseudo-random inputs where later merge levels deflate strongly. At 𝑁 = 16,384, BR is 4.38× faster than internal values-only D&C on uniform inputs and 3.43× faster on normal inputs, while preserving the same linear workspace advantage. The structured and reduced-dense cases show that this advantage is not limited to pseudo-random spectra: on 56 threads, BR is 3.79× faster on the ILP64 Toeplitz row and 6.13× faster on the ILP64 clustered row at 𝑁 = 32,768, and 5.63× faster on the densederived tridiagonal input. The larger internal-D&C rows use ILP64 calls because their queried LWORK exceeds the LP64 integer range; their eigenvalues are checked against DSTERF with the same normalized forward- and backward-error metrics used for BR. Thus the internal-D&C evidence supports BR as both a state-reduction result and a practical faster values-only D&C path on the measured runnable cases. 5.6
GPU D&C Against cuSOLVER
The H100 experiment tests whether the same values-only state reduction is useful in a GPU D&C setting. Table 4 compares our GPU EVD prototype’s values-only D&C path with cuSOLVER Xstedc on the same tridiagonal inputs. The comparison is intentionally DC-only: both paths start from 𝐷/𝐸, compute eigenvalues only, and exclude dense reduction and back transformation. For the reduced-dense row, a dense symmetric input is first reduced to tridiagonal form; the
Table 3. BR versus MKL internal values-only D&C. The internal path calls DLAED0(ICOMPQ=0) and is output-equivalent to BR. Ratios are internal/BR, so values above one mean BR is faster. Rows at 𝑁 = 32,768 and the reduced-dense row use ILP64 internal-D&C calls when the LP64 LWORK count is insufficient. Mode
𝑁
Int. WS
BR WS Int. Time BR Time Int./BR
1-thread fixed-seed pseudo-random inputs Uniform 4,096 386 MiB 0.61 MiB 41.4 ms 13.2 ms Uniform 8,192 1.50 GiB 1.22 MiB 83.8 ms 27.0 ms Uniform 16,384 6.01 GiB 2.44 MiB 173.9 ms 54.4 ms Uniform 32,768 24.0 GiB 4.88 MiB 369.8 ms 105.9 ms
3.15× 3.10× 3.20× 3.49×
Normal Normal Normal Normal
4,096 386 MiB 0.61 MiB 40.0 ms 8,192 1.50 GiB 1.22 MiB 80.5 ms 16,384 6.01 GiB 2.44 MiB 166.3 ms 32,768 24.0 GiB 4.88 MiB 332.0 ms
10.5 ms 21.3 ms 43.4 ms 89.2 ms
3.81× 3.77× 3.84× 3.72×
56-thread ratio on runnable internal-D&C rows Uniform 16,384 6.01 GiB 2.44 MiB 15.6 ms Normal 16,384 6.01 GiB 2.44 MiB 11.0 ms Toeplitz 16,384 6.01 GiB 2.44 MiB 3.04 s Clustered 16,384 6.01 GiB 2.44 MiB 5.80 s Toeplitz 32,768 24.0 GiB 4.88 MiB 11.90 s Clustered 32,768 24.0 GiB 4.88 MiB 23.27 s Reduced dense 49,152 54.0 GiB 7.31 MiB 44.8 s
3.57 ms 3.20 ms 812 ms 972 ms 3.14 s 3.80 s 7.95 s
4.38× 3.43× 3.75× 5.97× 3.79× 6.13× 5.63×
shared preparation time is excluded, so the row still measures only the D&C solver stage. The results show two effects. First, cuSOLVER’s queried device workspace grows quadratically and reaches 54.0 GiB at 𝑁 = 49152, while our GPU EVD prototype uses a compact 𝐷/𝐸 cache and a 4096-column temporary slab, about 1.5 GiB for the largest case. This is the GPU form of the state-management issue addressed by BR: even before full dense EVD is considered, the public D&C path carries a large workspace. Second, after removing a serial split-metadata update from the GPU path, our GPU EVD prototype is faster than cuSOLVER on every measured 𝑁 = 49152 row in Table 4. On the fixed-seed pseudo-random inputs, the speedups are 1.86× for normal and 1.74× for uniform. The dense-derived and structured rows show larger gains: 2.60× on the reduceddense tridiagonal input, 5.73× on Toeplitz, and 3.57× on clustered inputs. Thus the GPU data supports the same interpretation as the CPU data. BR makes a much lower-memory D&C design point available, and the speedups are largest when conventional D&C moves substantially more state than the boundary-row path needs. Profiling explains why the pseudo-random speedups are smaller than the structured ones. Before the split-metadata fix, a serial split-update kernel consumed about 72 ms at 𝑁 = 49152 on normal and uniform inputs, larger than the secular root-solving time. Parallelizing this metadata update reduces the split-update cost to microsecond-scale overhead. The remaining dominant kernels are boundary-vector construction, about 42 ms, and local eigenvalue sorting, about
Zhan and Zhang
Table 4. H100 values-only tridiagonal D&C comparison with cuSOLVER Xstedc. Larger cuSOLVER/ours ratios mean our GPU EVD prototype is faster.
local secular-equation ingredients as D&C and constructs the same merge vectors in exact arithmetic.
6 Mode Normal Uniform Reduced dense Toeplitz Clustered
𝑁 Our GPU EVD cuSOLVER cuSOLVER/ours 49,152 49,152 49,152 49,152 49,152
78.73 ms 85.75 ms 683.18 ms 480.81 ms 697.22 ms
146.71 ms 149.09 ms 1.78 s 2.75 s 2.49 s
1.86× 1.74× 2.60× 5.73× 3.57×
15 ms. These costs are mostly per-level D&C bookkeeping rather than secular arithmetic, so they are most visible on pseudo-random spectra where deflation makes the root solves relatively cheap. 5.7
Effect of Spectrum Structure
The spectrum strongly affects the observed speedup, which is why the paper reports pseudo-random and structured cases together. The uniform and normal families have substantial deflation and produce near-linear BR timings in the measured large-𝑁 range. Toeplitz and clustered families are harder: for 𝑁 ≥ 4096, both DSTERF and BR fit close to quadratic growth. The Toeplitz fits are DSTERF ∝ 𝑁 1.972 , BR 1T ∝ 𝑁 1.975 , and BR 16T ∝ 𝑁 1.930 over the available 16-thread range. The clustered fits are similarly quadratic, with DSTERF ∝ 𝑁 1.953 , BR 1T ∝ 𝑁 1.979 , and BR 16T ∝ 𝑁 1.941 . In these harder cases, BR should be described as a lower-memory D&C formulation with useful constant-factor and threading benefits, not as a universally growing-speedup replacement for QR/QL. This caveat is important for interpreting the main result. BR removes dense eigenvector-derived state from valuesonly D&C and exposes more parallelism than DSTERF; it does not remove the secular root-solving work that remains when deflation is weak. The measurements therefore support a conditional performance claim: BR is especially effective for spectra where D&C deflation and selected-row propagation keep the merge work close to linear, while its most general guarantee is linear auxiliary storage. 5.8
Numerical Accuracy
For all benchmark rows with a DSTERF reference, BR produced eigenvalues close to the reference under the normalized metrics defined above. The fixed-seed pseudo-random cross-checks cover uniform and normal diagonal distributions through 𝑁 = 65,536, while the structured checks cover Toeplitz, clustered, and reduced-dense tridiagonal inputs. We report forward and backward errors rather than raw eigenvalue differences so that the reduced-dense case, whose spectrum has ∥𝜆∥ ∞ = 4.93×104 , is compared on the same scale as the synthetic cases. These normalized checks are consistent with the conditioned-error argument in Section 3.2.2: BR changes which rows are propagated, but it uses the same
Conclusion
This paper presents boundary-row D&C, a values-only tridiagonal eigensolver that keeps the parallel D&C merge structure while removing the dense eigenvector-derived state that makes conventional values-only D&C memory-heavy. The central idea is to propagate only the boundary rows needed by future rank-one merges. This state is sufficient to construct the same secular problems as conventional D&C in exact arithmetic, reduces persistent auxiliary storage from quadratic to linear, and maps naturally to both CPU and GPU implementations. The prototypes show the practical effect of this reduced state. On CPUs, BR reaches problem sizes near the LP64 workspace limit and reduces the internal values-only D&C state from quadratic replay storage to linear boundary-row storage. This is not a claim that BR uses less memory than QR/QL, which remains the minimal-storage eigenvalue-only baseline; rather, BR keeps D&C-style parallelism while avoiding the quadratic state of conventional values-only D&C. With parallel leaf initialization, the CPU BR prototype is faster than both DSTERF and MKL’s internal D&C path on the completed pseudo-random and structured rows. ILP64 internal-D&C runs confirm that larger cases can be executed when the interface and memory budget allow it, but the quadratic replay state remains expensive: 𝑁 = 32,768 already requires about 24.0 GiB, and BR remains faster on the completed ILP64 rows. On an H100, the values-only D&C path avoids the large quadratic cuSOLVER workspace and is faster than cuSOLVER on all measured 𝑁 = 49152 tridiagonal-stage cases, with 1.74–5.73× speedups on the synthetic and structured families after parallelizing split-metadata updates. These results make the boundary clear. BR is not a universal subquadratic eigensolver; it makes the D&C design point available to eigenvalue-only workloads without carrying full eigenvector state.
References [1] Hervé Abdi and Lynne J Williams. 2010. Principal component analysis. Wiley interdisciplinary reviews: computational statistics 2, 4 (2010), 433– 459. [2] Edward Anderson, Zhaojun Bai, Christian Bischof, L Susan Blackford, James Demmel, Jack Dongarra, Jeremy Du Croz, Anne Greenbaum, Sven Hammarling, Alan McKenney, et al. 1999. LAPACK Users’ guide. SIAM. [3] Paolo Bientinesi, Inderjit S Dhillon, and Robert A Van De Geijn. 2005. A parallel eigensolver for dense symmetric matrices based on multiple relatively robust representations. SIAM Journal on Scientific Computing 27, 1 (2005), 43–66. [4] J. J. M. Cuppen. 1981. A divide and conquer method for the symmetric tridiagonal eigenproblem. Numer. Math. 36 (1981), 177–195. [5] RZ Dautov, AD Lyashko, and SI Solov’ev. 1994. The bisection method for symmetric eigenvalue problems with a parameter entering nonlinearly. (1994).
Reducing Internal State in Eigenvalue-Only Divide-and-Conquer Tridiagonal Eigensolvers
[6] Inderjit S Dhillon and Beresford N Parlett. 2004. Multiple representations to compute orthogonal eigenvectors of symmetric tridiagonal matrices. Linear Algebra Appl. 387 (2004), 1–28. [7] Gene H Golub and Charles F Van Loan. 2013. Matrix computations. JHU press. [8] Roger Grimes, Henry Krakauer, John Lewis, Horst Simon, and SuHai Wei. 1987. The solution of large dense generalized eigenvalue problems on the Cray X-MP/24 with SSD. J. Comput. Phys. 69, 2 (1987), 471–481. [9] Ming Gu and Stanley C Eisenstat. 1994. A stable and efficient algorithm for the rank-one modification of the symmetric eigenproblem. SIAM journal on Matrix Analysis and Applications 15, 4 (1994), 1266–1276. [10] Ming Gu and Stanley C Eisenstat. 1995. A divide-and-conquer algorithm for the symmetric tridiagonal eigenproblem. SIAM J. Matrix Anal. Appl. 16, 1 (1995), 172–191. [11] Vineet Gupta, Tomer Koren, and Yoram Singer. 2018. Shampoo: Preconditioned stochastic tensor optimization. In International Conference on Machine Learning. PMLR, 1842–1850. [12] MAGMA Project. 2026. MAGMA 2.10.0 Documentation: Symmetric/Hermitian Eigenvalue Drivers. University of Tennessee, Knoxville. https://icl.utk.edu/projectsfiles/magma/doxygen/group_ _magma__heevd.html. MAGMA Source Code: dsyevd.cpp, [13] MAGMA Project. 2026. dsyevd_gpu.cpp, and dstedx.cpp. https://github.com/icl-utk-edu/ magma/tree/master/src. [14] A Melman. 1995. Numerical solution of a secular equation. Numer. Math. 69 (1995), 483–493. [15] NVIDIA Corporation. 2026. NVIDIA cuSOLVER Documentation. NVIDIA Corporation. https://docs.nvidia.com/cuda/cusolver/. [16] Matt Probert. 2011. Electronic Structure: Basic Theory and Practical Methods, by Richard M. Martin: Scope: graduate level textbook. Level: theoretical materials scientists/condensed matter physicists/computational chemists. [17] Notker Rösch, Sven Krüger, Vladimir A Nasluzov, and Alexei V Matveev. 2005. ParaGauss: The density functional program paragauss for complex systems in chemistry. In High Performance Computing in Science and Engineering, Garching 2004. Springer, 285–296. [18] Ruchi Shah, Shaoshuai Zhang, Ying Lin, and Panruo Wu. 2019. xSVM: Scalable distributed kernel support vector machine training. In 2019 IEEE International Conference on Big Data (Big Data). IEEE, 155–164. [19] Françoise Tisseur and Jack Dongarra. 1999. A parallel divide and conquer algorithm for the symmetric eigenvalue problem on distributed memory architectures. SIAM Journal on Scientific Computing 20, 6 (1999), 2223–2236. [20] Stanimire Tomov, Rajib Nath, Peng Du, and Jack Dongarra. 2011. MAGMA Users’ Guide. ICL, UTK (November 2009) (2011). [21] Hansheng Wang, Zhekai Duan, Zitian Zhao, Siqi Wu, Saiqi Zheng, Qiao Li, Xu Jiang, and Shaoshuai Zhang. 2025. Improving Tridiagonalization Performance on GPU Architectures. In Proceedings of the 30th ACM SIGPLAN Annual Symposium on Principles and Practice of Parallel Programming. 469–480. [22] David S Watkins. 1982. Understanding the QR algorithm. SIAM review 24, 4 (1982), 427–440. [23] Shaoshuai Zhang, Ruchi Shah, and Panruo Wu. 2020. TensorSVM: accelerating kernel machines with tensor engine. In Proceedings of the 34th ACM International Conference on Supercomputing. 1–11.
A
Proof Details
This appendix records the BR-specific algebra and conditioned error-propagation details used by Section 3.2.2. It does not reprove the stability of a complete secular-equation solver, nor does it formalize all IEEE 754 behavior. Those
components are treated as external numerical-analysis assumptions. A.1
Conventions
For a local eigenvector block 𝑄 ∈ R𝑚×𝑛 , define blo(𝑄) = 𝑄 1,:,
bhi(𝑄) = 𝑄𝑚,: .
The boundary state is 𝐵(𝑄) = (blo(𝑄), bhi(𝑄)). A secular problem is written as 𝑃 = (𝐷, 𝑧, 𝜌),
𝐴(𝑃) = diag(𝐷) + 𝜌𝑧𝑧𝑇 ,
where 𝐷 is the pole vector and 𝜌 carries the rank-one strength. For perturbation bookkeeping we use 1/2 2 b − 𝐷 + ∥b 𝑃b − 𝑃 = 𝐷 𝑧 − 𝑧 ∥ 22 + |b 𝜌 − 𝜌 |2 . 𝑃
2
This unweighted norm can be replaced by a weighted norm; the corresponding perturbation constants then absorb the scaling of 𝐷, 𝑧, and 𝜌. A.2
Row Selection and Metadata
Lemma A.1 (Column permutations commute with row selection). Let Π be a column permutation matrix and let 𝜎 be an arbitrary row list. Then (𝑄Π)𝜎 = 𝑄 𝜎 Π. Proof. For every requested row 𝑖 and column 𝑗, both sides equal 𝑄 𝜎 (𝑖 ),𝜋 ( 𝑗 ) . Lemma A.2 (Local rotations commute with row selection). Let 𝐺 be a Givens rotation or any other right-multiplied local column transformation. Then (𝑄𝐺)𝜎 = 𝑄 𝜎 𝐺 . Proof. This is the same entrywise argument as Lemma 3.2; a Givens rotation changes only two columns, and every affected entry is a fixed linear combination of the two selectedrow entries. In an implementation, column signs, deflation-induced permutations, and sorting metadata must be applied consistently to both conventional D&C and BR. Equivalently, one may regard the stored state as 𝐵(𝑄, 𝜋, 𝑠) = 𝑠 ⊙ 𝑄 1,𝜋 , 𝑠 ⊙ 𝑄𝑚,𝜋 , where 𝜋 is a column permutation and 𝑠 is a column-sign convention. Under this convention, the algebraic equality in Theorem 3.3 is exact. A.3
Compact Delta Reconstruction
For a secular root written as 𝜆 = 𝐷 org + 𝜏, define Δ𝑖 = 𝐷𝑖 − 𝐷 org − 𝜏 . Near-pole entries may be cached to avoid cancellation, while far entries are reconstructed by the formula above.
Zhan and Zhang
Lemma A.3 (Compact delta correctness). Assume the compact representation has legal near indices, non-near indices equal the complement of the cached near set, and cached near entries equal the corresponding exact Δ𝑖 . If every non-near entry is reconstructed as 𝐷𝑖 − 𝐷 org − 𝜏, then the reconstruction returns the exact vector Δ in exact arithmetic. Proof. For a near index, the returned value is the cached exact entry. For a non-near index, the returned value is exactly the defining formula for Δ𝑖 . Thus all entries match. Lemma A.4 (Pointwise delta error to vector error). If |b Δ𝑖 − Δ𝑖 | ≤ 𝜂 for all 𝑖, then √ b Δ − Δ ≤ 𝑛 𝜂. 2
Proof. Square the pointwise bound, sum over 𝑖, and take square roots. Lemma A.5 (Far-entry absolute rounding bound). Let a far entry be computed as b𝑖 − 𝐷 borg ) − 𝜏b), b Δ𝑖 = fl(fl(𝐷 and assume fl(𝑥 − 𝑦) = (𝑥 − 𝑦) (1 + 𝛿) with |𝛿 | ≤ 𝑢. If b𝑖 − 𝐷𝑖 | ≤ 𝜂𝐷 , |𝐷
borg − 𝐷 org | ≤ 𝜂𝐷 , |𝐷
|b 𝜏 − 𝜏 | ≤ 𝜂𝜏 ,
then b𝑖 − 𝐷 borg | + 𝑢 | fl(𝐷 b𝑖 − 𝐷 borg ) − 𝜏b| + 2𝜂𝐷 + 𝜂𝜏 . |b Δ𝑖 − Δ𝑖 | ≤ 𝑢 | 𝐷
then 𝑌b − 𝑌
Substituting 𝑎 and 𝑏 gives the result.
Proof. By the triangle inequality, every entry of 𝑌b − 𝑌 is bounded by 𝜂 round + 𝜂 input . The Frobenius bound follows by summing the squared entrywise bounds over the 𝑟 ℓ entries. The dot-product model gives 𝜂 round = 𝛾𝑘 𝑀abs whenever 𝑘 ∑︁
for a positive constant 𝛼.
for all actual dot products. This quantity must bound the sum of absolute products, not the cancellation-reduced exact dot product. A.5
Assumption A.1 (Dot-product error model). For a length-𝑘 dot product, the implementation satisfies 𝑘 ∑︁ 𝑡 =1
𝑥𝑡 𝑦𝑡 ≤ 𝛾𝑘
WPROD Product Chains
Some secular-vector entries are computed as products 𝑚 Ö 𝑓𝑖 . 𝑊 = 𝑊0 𝑖=1
Assume every effective factor satisfies 𝑓b𝑖 = 𝑓𝑖 (1 + 𝛿𝑖 ),
|𝛿𝑖 | ≤ 𝑢.
Then b = 𝑊 (1 + 𝜃 ), 𝑊
𝜃=
𝑚 Ö
(1 + 𝛿𝑖 ) − 1.
𝑖=1
If 0 ≤ 𝑢 ≤ 1, then This proves only composition of relative factor perturbations; stability of the individual factor evaluations must be supplied separately. A.6
Selected-Row Rounding Error
flDot(𝑥, 𝑦) −
b𝜎 )𝑖𝑡 𝑆b𝑡 𝑗 | ≤ 𝑀abs |(𝑄
𝑡 =1
|𝜃 | ≤ Γ𝑚± (𝑢) = max{(1 + 𝑢)𝑚 − 1, 1 − (1 − 𝑢)𝑚 }.
This is an absolute error bound. A relative bound for far entries requires a separation condition such as |Δ𝑖 | ≥ 𝛼 |𝐷𝑖 − 𝐷 org | + |𝜏 | ,
A.4
𝐹
where 𝑄 𝜎 ∈ R𝑟 ×𝑘 and 𝑆 ∈ R𝑘 ×ℓ .
b𝑖 − 𝐷 borg and 𝑏 = fl(𝑎) − 𝜏b. Then Proof. Let 𝑎 = 𝐷 |b Δ𝑖 − Δ𝑖 | ≤ | fl(𝑏) − 𝑏 | + |𝑏 − Δ𝑖 | ≤ 𝑢 |𝑏 | + 𝑢 |𝑎| + 2𝜂𝐷 + 𝜂𝜏 .
√ ≤ 𝑟 ℓ (𝜂 round + 𝜂 input ),
𝑘 ∑︁
|𝑥𝑡 𝑦𝑡 |.
𝑡 =1
For sequential summation with 𝑘𝑢 < 1, one may take 𝛾𝑘 = 𝑘𝑢/(1 − 𝑘𝑢). Proposition A.6 (Selected-row update with input perturb𝜎 , 𝑆). b If every bation). Let 𝑌 = 𝑄 𝜎 𝑆 and let 𝑌b = flMatMul(𝑄 output entry satisfies b𝜎 𝑆) b 𝑖 𝑗 | ≤ 𝜂 round, |𝑌b𝑖 𝑗 − (𝑄 and b𝜎 𝑆) b 𝑖 𝑗 − (𝑄 𝜎 𝑆)𝑖 𝑗 | ≤ 𝜂 input, |(𝑄
Tree-Level Propagation
Theorem A.7 (Conditional boundary-error propagation). bℓ − 𝐵 ℓ ≤ 𝐸 ℓ . For every internal Assume each leaf ℓ satisfies 𝐵 node 𝑣, assume the local merge routine has a declared output bound b𝑣 − 𝐵 𝑣 ≤ 𝐹 𝑣 (𝜂𝐿 , 𝜂𝑅 , 𝜂𝐷 , 𝜂 𝜌 , 𝜂 round, 𝜅 𝑣 ) ≡ 𝐸 𝑣 , 𝐵 whenever the left and right child boundary errors satisfy their input bounds. Then every node in the BR merge tree satisfies its declared bound. Moreover, every internal secular vector satisfies √︃ ∥b 𝑧 𝑣 − 𝑧 𝑣 ∥ 2 ≤ 𝜂𝐿2 (𝑣) + 𝜂𝑅2 (𝑣) ≤ 𝜂𝐿 (𝑣) + 𝜂𝑅 (𝑣) . Proof. The proof is induction on tree height. Leaf nodes satisfy the hypothesis. For an internal node, the induction hypothesis supplies the child boundary-error bounds. Equation (10) gives the secular-vector perturbation, and the asb𝑣 − 𝐵 𝑣 ≤ 𝐸 𝑣 . sumed local merge bound gives 𝐵
Reducing Internal State in Eigenvalue-Only Divide-and-Conquer Tridiagonal Eigensolvers
Theorem A.8 (Root values-only error decomposition). Let 𝑃 = (𝐷, 𝑧, 𝜌) be the exact root secular problem and 𝑃b the BR-computed perturbed problem. If 𝑃b − 𝑃
≤ 𝜂 prob, 𝑃
and the local eigenvalue map obeys b − 𝜆(𝑃) 𝜆(𝑃)
2
≤ 𝜅 pert𝜂 prob,
while the root solver obeys b 𝜆b − 𝜆(𝑃)
2
≤ 𝜂 solver,
then 𝜆b − 𝜆(𝑃)
2
≤ 𝜂 solver + 𝜅 pert𝜂 prob .
Proof. Apply the triangle inequality: 𝜆b − 𝜆(𝑃)
2
b ≤ 𝜆b − 𝜆(𝑃)
2
b − 𝜆(𝑃) + 𝜆(𝑃)
. 2
The two terms are bounded by the solver and perturbation assumptions. For a simplified balanced-tree model with per-level additive local error 𝜇, the loose recurrence 𝐸ℎ+1 = 2𝐸ℎ + 𝜇 gives 𝐸ℎ = 2ℎ 𝜂 0 + (2ℎ − 1)𝜇. √ Using the sharper concatenation bound gives 𝑈ℎ+1 = 2𝑈ℎ + 𝜇 and √ √ ℎ ( 2)ℎ − 1 𝜇. 𝑈ℎ = ( 2) 𝜂 0 + √ 2−1 These recurrences are illustrative upper-bound models; the actual local Lipschitz constants and secular conditioning are contained in 𝐹 𝑣 .