[Paper Review] Accelerating eigenvector and pseudospectra computation using blocked multi-shift triangular solves

Tim Moon, Jack Poulson|arXiv (Cornell University)|Jul 6, 2016

Parallel Computing and Optimization Techniques17 references3 citations

TL;DR

This paper proposes blocked multi-shift triangular solve algorithms that leverage Level 3 BLAS to accelerate eigenvector and pseudospectra computations. By reformulating triangular solves with multiple shifts and incorporating safeguards for ill-conditioned matrices, the method achieves up to 60× speedup over LAPACK for triangular eigenvectors and 9× speedup for pseudospectra relative to EigTool, significantly improving performance in eigenvalue analysis workloads.

ABSTRACT

Multi-shift triangular solves are basic linear algebra calculations with applications in eigenvector and pseudospectra computation. We propose blocked algorithms that efficiently exploit Level 3 BLAS to perform multi-shift triangular solves and safe multi-shift triangular solves. Numerical experiments indicate that computing triangular eigenvectors with a safe multi-shift triangular solve achieves speedups by a factor of 60 relative to LAPACK. This algorithm accelerates the calculation of general eigenvectors threefold. When using multi-shift triangular solves to compute pseudospectra, we report ninefold speedups relative to EigTool.

Motivation & Objective

• To address performance bottlenecks in eigenvector and pseudospectra computation caused by reliance on Level 2 BLAS in existing LAPACK routines.
• To enable efficient computation of triangular and general eigenvectors by replacing unblocked back substitution with blocked multi-shift triangular solves.
• To develop a safe, numerically stable variant of the multi-shift triangular solve that prevents overflow and division by zero in ill-conditioned cases.
• To accelerate pseudospectra computation by enabling Level 3 BLAS-optimized multi-shift solves across multiple shift points.
• To provide a foundation for future integration into LAPACK through Fortran implementations and support for real arithmetic.

Proposed method

• Proposes a blocked algorithm for multi-shift triangular solves that replaces Level 2 BLAS operations with Level 3 BLAS routines like xGEMM to improve data locality and cache efficiency.
• Adapts the standard blocked back substitution algorithm to handle multiple shifts simultaneously by applying the shift to the diagonal blocks of the triangular matrix.
• Introduces a safe multi-shift triangular solve using dynamic scaling and error monitoring to prevent numerical overflow and maintain backward stability.
• Employs a safeguarded iterative scaling strategy that monitors the infinity norm of intermediate vectors and adjusts scaling factors to stay within safe bounds.
• Generalizes the method to the generalized eigenvalue problem by solving (U - λV)x = b using modified blocked algorithms with xTRSV and xGEMM operations.
• Uses the generalized safe multi-shift solve as a core subroutine in a robust eigensolver pipeline for both triangular and general matrices.

Experimental results

Research questions

• RQ1 Can blocked multi-shift triangular solves using Level 3 BLAS significantly reduce computational time for eigenvector computation compared to LAPACK’s Level 2 BLAS-based xTREVC?
• RQ2 Can a safe multi-shift triangular solve algorithm prevent numerical overflow and division by zero while maintaining backward stability in ill-conditioned cases?
• RQ3 To what extent can multi-shift triangular solves accelerate pseudospectra computation, especially in large-scale problems?
• RQ4 How does the performance of the proposed algorithm scale with matrix size and number of shifts?
• RQ5 Can the proposed algorithms be generalized to the generalized eigenvalue problem and integrated into standard linear algebra libraries?

Key findings

• Computing triangular eigenvectors using the safe multi-shift triangular solve achieves a 60× speedup over LAPACK’s xTREVC routine.
• General eigenvector computation is accelerated by a factor of three using the proposed method compared to the LAPACK baseline.
• Pseudospectra computation using the new algorithm achieves a ninefold speedup over EigTool when computing the resolvent norm for a 3200×3200 matrix at 10,000 grid points.
• The performance gains are attributed to efficient use of Level 3 BLAS, particularly xGEMM, which improves data reuse and cache utilization.
• The safe multi-shift triangular solve maintains numerical stability by dynamically scaling solution vectors to prevent overflow, even for ill-conditioned matrices.
• The generalized multi-shift triangular solve extends the method to the generalized eigenvalue problem, enabling robust and efficient computation of generalized eigenvectors.