The Bunch-Kaufman algorithm and Aasen's algorithm are two of the most widely used methods for solving symmetric indefinite linear systems, yet they both are known to suffer from occasional numerical instability due to potentially exponential element growth or unbounded entries in the matrix factorization. In this work, we develop a randomized complete pivoting (RCP) algorithm for solving symmetric indefinite linear systems. RCP is comparable to the Bunch-Kaufman algorithm and Aasen's algorithm in computational efficiency, yet enjoys theoretical element growth and bounded entries in the factorization comparable to that of complete-pivoting, up to a theoretical failure probability that exponentially decays with an oversampling parameter. Our finite precision analysis shows that RCP is as numerically stable as Gaussian elimination with complete pivoting, and RCP has been observed to be numerically stable in our extensive numerical experiments. 1 sparse symmetric indefinite linear systems. All but Aasen's algorithm use diagonal pivoting to factorize A into block LDL T , where L is unit lower triangular and D is symmetric block diagonal with each block of order 1 or 2. Aasen's algorithm factorizes A into LT L T , where T is a symmetric tridiagonal matrix.The block LDL T factorization is a generalization of the Cholesky factorization, which requires the input matrix A to be positive semi-definite. While the Cholesky factorization is numerically stable with or without diagonal pivoting [21], block LDL T factorization with partial pivoting can have numerical instability issues [5,24]. There are a number of strategies to choose permutation matrices for numerical stability. Bunch and Parlett [14] proposed a complete pivoting strategy, which requires searching the whole Schur complement at each stage of the block LDL T factorization and therefore requires up to O(n 3 ) comparisons. Bunch [12] proved that Bunch-Parlett algorithm satisfies a backward error bound almost as good as that for Gaussian elimination with complete pivoting (GECP), and its element growth factor is within a factor 3.07(n − 1) 0.446 of Wilkinson's element growth factor bound [38] for GECP. To reduce the number of comparisons, Bunch and Kaufman devised a partial pivoting strategy, which searches at most two columns at each stage and the number of comparisons is reduced from O(n 3 ) to O(n 2 ). However, the multipliers can't be controlled with Bunch-Kaufman partial pivoting strategy. To overcome this instability problem, Ashcraft, Grimes, and Lewis [5] proposed a bounded Bunch-Kaufman algorithm, where multipliers can be bounded near one, at a potentially high cost in comparisons. In the worst-case scenario, bounded Bunch-Kaufman is no better than the Bunch-Parlett algorithm, while in the best-case scenario bounded Bunch-Kaufman costs no more than the Bunch-Kaufman algorithm. LAPACK [2] routines SYSV and SYSVX, and LINPACK [17] routines SIFA and SISL are all based on the Bunch-Kaufman algorithm, whereas LAPACK routines SYSV rook and SYSV aa are based on ...

We present the Flip-Flop Spectrum-Revealing QR (Flip-Flop SRQR) factorization, a significantly faster and more reliable variant of the QLP factorization of Stewart for low-rank matrix approximations. Flip-Flop SRQR uses SRQR factorization to initialize a partial column-pivoted QR factorization and then computes a partial LQ factorization. As observed by Stewart in his original QLP work, Flip-Flop SRQR tracks the exact singular values with "considerable fidelity". We develop singular value lower bounds and residual error upper bounds for the Flip-Flop SRQR factorization. In situations where singular values of the input matrix decay relatively quickly, the low-rank approximation computed by Flip-Flop SRQR is guaranteed to be as accurate as the truncated SVD. We also perform a complexity analysis to show that Flip-Flop SRQR is faster than the randomized subspace iteration for approximating the SVD, the standard method used in the Matlab tensor toolbox. We additionally compare Flip-Flop SRQR with alternatives on two applications, a tensor approximation and a nuclear norm minimization, to demonstrate its efficiency and effectiveness.

No abstract

Aasen's algorithm factorizes a symmetric indefinite matrix A as A = P T LT L T P , where P is a permutation matrix, L is unit lower triangular with its first column being the first column of the identity matrix, and T is tridiagonal. In this note, we provide a growth factor upper bound for Aasen's algorithm which is much smaller than that given by Higham. We also show that the upper bound we have given is not tight when the matrix dimension is greater than or equal to 6.

No abstract