The Power of Bidiagonal Matrices

An upper bidiagonal matrix

\notag B = \begin{bmatrix} b_{11} & b_{12} & & \\ & b_{22} & \ddots & \\ & & \ddots & b_{n-1,n} \\ & & & b_{nn} \end{bmatrix} \in\mathbb{C}^{n\times n}

depends on just 2n-1 parameters, which appear on the main diagonal and the superdiagonal.

Such matrices arise commonly, for example as the U factor and the transpose of the L factor in the LU factorization of tridiagonal matrices, and as the intermediate matrix in the computation of the singular value decomposition by the Golub–Reinsch algorithm.

Bidiagonal matrices have many interesting properties, especially when we consider products of them. We describe some of their properties here.

Inverse

Consider the inverse of a nonsingular 4 \times 4 bidiagonal matrix:

\notag \begin{bmatrix} a & x & 0 & 0 \\ & b & y & 0 \\ & & c & z \\ & & & d \end{bmatrix}^{-1} = \begin{bmatrix} \frac{1}{a} & -\frac{x}{a\,b} & \frac{x\,y}{a\,b\,c} & -\frac{x\,y\,z}{a\,b\,c\,d} \\[3pt] & \frac{1}{b} & -\frac{y}{b\,c} & \frac{y\,z}{b\,c\,d} \\[3pt] & & \frac{1}{c} & -\frac{z}{c\,d} \\[3pt] & & & \frac{1}{d} \end{bmatrix}.

Notice that every element in the upper triangle is a product of off-diagonal elements of B and inverses of diagonal elements, that the superdiagonals have alternating signs attached, and that there are no additions. These properties hold for general n.

The consequences include:

  • the absolute values of the elements of B^{-1} depend on |B| = (|b_{ij}|), as there is no cancellation in the formulas for B^{-1} (unlike for general triangular matrices),
  • the singular values of B depend only on |B| (less obvious, and provable by a scaling argument).

Matrix Function

There is a simple formula for any function of B, not just the inverse. Here, the f[\dots] term is a divided difference.

Theorem 1. If B\in\mathbb{C}^{n\times n} is upper bidiagonal then F = f(B) is upper triangular with f_{ii} = f(t_{ii}) and

\notag f_{ij} = b_{i,i+1} b_{i+1,i+2} \dots b_{j-1,j}\, f[b_{ii},b_{i+1,i+1}, \dots, b_{jj}], \quad j > i.

A special case is the formula for a function of an m\times m Jordan block:

\notag f\left( \begin{bmatrix} \lambda & 1 & & \\ & \lambda & \ddots & \\ & & \ddots & 1 \\ & & & \lambda \end{bmatrix} \right) = \begin{bmatrix} f(\lambda) & f'(\lambda) & \dots & \displaystyle\frac{f^{(m-1)}(\lambda)}{(m-1)!} \\ & f(\lambda) & \ddots & \vdots \\ & & \ddots & f'(\lambda) \\ & & & f(\lambda) \end{bmatrix}.

Condition Number of Product of Bidiagonals

Perhaps the most interesting result is that if we have a factorization of a matrix A\in\mathbb{C}^{n\times n} into a product of 2n-1 bidiagonal matrices then we can compute the condition number \kappa_{\infty}(A) = \|A\|_{\infty} \|A^{-1}\|_{\infty} in O(n^2) flops and to high accuracy in floating-point arithmetic. Every matrix has such a factorization and in some cases the entries of the factors are known as explicit formulas. For example, we can compute the condition number of the Hilbert matrix H_n to high accuracy in O(n^2) flops.

n \kappa_{\infty}(H_n) Relative error for fast algorithm
4 2.84e4 1.28e-16
8 3.39e10 2.25e-16
16 5.06e22 3.67e-17
32 1.36e47 1.75e-15
64 1.10e96 1.77e-15

Any attempt to compute \kappa_{\infty}(H_n) by explicitly forming H_n is doomed to failure by the rounding errors incurred in the formation, no matter what algorithm is used for the computation.

All this analysis, and much more, is contained in

which is based on the Hans Schneider Prize talk that I gave at The 25th Conference of the International Linear Algebra Society, Madrid, Spain, June 12-16, 2023.

Leave a comment