We have been investigating whether one-time Householder tridiagonalization before Jacobi can improve cusolverDnDsyevj performance.
The combination is usually considered unhelpful because tridiagonalization costs O(n³), while subsequent Jacobi rotations do not preserve the tridiagonal zeros. But that does not necessarily imply that the tridiagonal starting form cannot accelerate Jacobi convergence.
We compared
A → cusolverDnDsyevj
with
A → cusolverDnDsytrd → T → cusolverDnDsyevj,
where T = QᵀAQ.
We tested six dense symmetric matrix families for n = 16, 32, 64, 128, 256, 512, 1024.
For five of six families, preprocessing reduced the subsequent Jacobi sweep count, by 23.5% on average. Beginning at n = 32, the Jacobi-time reduction recovered the cost of sytrd and construction of T for all five favorable families; across those family/size combinations, the complete pipeline was 28.0% faster on average than standalone syevj.
The result was not universal. A deliberately adverse unequal-variance covariance family lost substantial diagonal Frobenius energy after Householder reduction, and both sweep count and runtime became worse.
Our current interpretation is that two effects compete.
First, after Householder tridiagonalization, all remaining off-diagonal energy is concentrated in the first band. For even n, one natural parallel first batch is
(1,2), (3,4), (5,6), …
For approximately evenly distributed off-diagonal Frobenius energy, these pivots contain about 1/(n−1) of the energy in a dense matrix, but about one half after tridiagonalization.
Second, Householder similarity can redistribute energy between the diagonal and off-diagonal parts. Favorable families either gained diagonal energy or changed little; the adverse covariance family lost substantial diagonal energy.
Our questions are:
- Does cusolverDnDsyevj use this parallel pivot ordering internally, or something structurally similar?
- If the behavior remains consistent across architectures and matrix families, could conditional tridiagonalization be useful as preprocessing for syevj in the small/medium matrix regime?
- We would welcome comments on the algebraic interpretation, and can provide working CUDA code or run additional tests ourselves.
Full derivation, tables, accuracy checks, and results: