A more accurate algorithm for computing the Christoffel transformation (Q2372952)

Consider a monic Jacobi matrix \(J\) associated with a real measure \(d\mu\) and let \(\alpha\) be a real number. If \(J-\alpha I =LU\) denotes the LU factorization without pivoting of \(J-\alpha I\), where \(L\) is unit lower triangular, then the matrix version of the basic Christoffel or Darboux transformation with shift \(\alpha\) is given by (1) \(J-\alpha I=LU\), \(\widetilde J=UL+\alpha I\), where \(\widetilde J\) is the monic Jacobi matrix associated with the measure \((x-\alpha)d\mu(x)\), and \(L\) is a monic bidiagonal matrix with first lower diagonal \(\ell_1,\ell_2,\ldots\), and \(U\) is a bidiagonal matrix with diagonal \(u_1,u_2,\ldots\) and unit first upper diagonal. Let \((M)_{n-1}\) denotes the leading principal submatrix of any matrix \(M\). For numerical analysis, a finite version of (1) is (2) \(J(B,G)-\alpha I=LU\), \(J(b,g)=(UL+\alpha I)_{n-1}\), where \(J(B,G)=(J)_n\), \(B=[B_1,\ldots,B_n]^T\), \(G=[G_1,\ldots,G_{n-1}]^T\), and \(J(b,g)=(\widetilde J)_{n-1}\). The transformation (2) is implemented in an inaccurate Algorithm 1. The authors derive a new stable and accurate Algorithm 2 from Algorithm 1 with the same cost of \(6n-9\) flops. A backward error analysis is performed and a tight first-order forward arror bound is computed for Algorithm 2.