Toeplitz circulant and displacement structured solvers

All vector and matrix indices below start at $0$. A rectangular Toeplitz matrix $T\in {\mathbb{C}}^{m\times n}$, where $m,n\ge 1$, is specified by lags $t_r\in \mathbb{C}$ for $-(n-1)\le r\le m-1$ and has entries $T_{ij}:=t_{i-j}(0\le i<m,0\le j<n)\mathrm{.}$ For a square matrix $T\in {\mathbb{C}}^{n\times n}$, Hermitian Toeplitz means $t_{-r}=\overline{t_r}$ for $0\le r<n$. It is positive definite when $x^{\ast }Tx>0$ for every nonzero $x\in {\mathbb{C}}^n$, where $x^{\ast }:={\bar{x}}^T$ and the product is the matrix product of Rectangular matrix multiplication and the identity matrix $I_n$, including zero-sized shapes.

For $N\ge 1$, a vector $c=(c_0,\dots ,c_{N-1}{)}^T\in {\mathbb{C}}^N$ generates the circulant matrix $C(c{)}_{ij}:=c_{(i-j)<mrow><mi\; mathvariant="normal">m</mi><mi\; mathvariant="normal">o</mi><mi\; mathvariant="normal">d</mi></mrow>N}(0\le i,j<N)\mathrm{.}$ For $n\ge 1$ and $0\le j<n$, let $e_j\in {\mathbb{C}}^n$ be the coordinate vector with $(e_j{)}_i=1$ when $i=j$ and $0$ otherwise. For $a\in {\mathbb{C}}^n$, $J_{n}a$ is the reversed vector, $(J_{n}a{)}_i=a_{n-1-i}$; $J_0$ is the unique map on the empty vector. The lower-triangular Toeplitz matrix of $a$ is $L(a)\in {\mathbb{C}}^{n\times n}$ with $L(a{)}_{ij}=a_{i-j}$ when $i\ge j$, and $0$ otherwise. The canonical unitary DFT is the map ${\mathcal{F}}_N$ of The unitary discrete Fourier transform on $\mathbb{Z}\mathrm{/}N\mathbb{Z}$; the unnormalised DFT used for operation counts is $D_N:=\sqrt{N}{\mathcal{F}}_N$.

These conventions specify matrix orientation, conjugation, normalization, and boundary indexing; in particular, circulant eigenvalues below use the negative-sign canonical DFT.

Let $s\ge 1$, $A\in {\mathbb{R}}^{s\times s}$, $x,b\in {\mathbb{R}}^s$, and let each row dot product be evaluated left to right: round each product, then accumulate the rounded products from the first term onward. After computing $\widehat{p}=\mathrm{fl}(Ax)$, compute $\widehat{r}=\mathrm{fl}(b-\widehat{p})$ componentwise. Assume the standard relative model with unit roundoff $u$, all products, intermediate sums and subtractions are zero or in the normal range, and $(s+1)u<1$. Put ${\gamma }_q:=qu\mathrm{/}(1-qu)$.

For the exact residual $r=b-Ax$, the specified operation DAG satisfies $\mathrm{\mid }r-\widehat{r}\mathrm{\mid }\le {\gamma }_{s+1}(\mathrm{\mid }b\mathrm{\mid }+\mathrm{\mid }A\mathrm{\mid }\mathrm{\mid }x\mathrm{\mid })$ componentwise, and therefore $\mathrm{\parallel }r{\mathrm{\parallel }}_{\mathrm{\infty }}\le \mathrm{\parallel }\widehat{r}{\mathrm{\parallel }}_{\mathrm{\infty }}+{\gamma }_{s+1}(\mathrm{\parallel }b{\mathrm{\parallel }}_{\mathrm{\infty }}+\mathrm{\parallel }A{\mathrm{\parallel }}_{\mathrm{\infty }}\mathrm{\parallel }x{\mathrm{\parallel }}_{\mathrm{\infty }}),$ where $\mathrm{\parallel }A{\mathrm{\parallel }}_{\mathrm{\infty }}:={\max }_i{\sum }_j\mathrm{\mid }{A}_{ij}\mathrm{\mid }$ and $\mathrm{\parallel }x{\mathrm{\parallel }}_{\mathrm{\infty }}:={\max }_i\mathrm{\mid }{x}_i\mathrm{\mid }$. If $d:=\mathrm{\parallel }A{\mathrm{\parallel }}_{\mathrm{\infty }}\mathrm{\parallel }x{\mathrm{\parallel }}_{\mathrm{\infty }}+\mathrm{\parallel }b{\mathrm{\parallel }}_{\mathrm{\infty }}>0$, the normwise backward-error ratio ${\eta }_{\mathrm{\infty }}:=\mathrm{\parallel }r{\mathrm{\parallel }}_{\mathrm{\infty }}\mathrm{/}d$ is at most ${\gamma }_{s+1}+\mathrm{\parallel }\widehat{r}{\mathrm{\parallel }}_{\mathrm{\infty }}\mathrm{/}d$; if $d=0$, then $r=0$.

For $n\ge 1$, complex $A\in {\mathbb{C}}^{n\times n}$ and $x,b\in {\mathbb{C}}^n$, define $R_{n}z:=(\mathrm{Re}z,\mathrm{Im}z{)}^T\in {\mathbb{R}}^{2n}$ and $\mathcal{R}(A)=\left(\begin{array}{cc}\mathrm{Re}A& -\mathrm{Im}A\\ \mathrm{Im}A& \mathrm{Re}A\end{array}\right)$; put $\tilde{x}=R_{n}x$ and $\tilde{b}=R_{n}b$. The following square expresses the realification identity proved by direct multiplication in step 1.1.

Applying the same real DAG to $\mathcal{R}(A),\tilde{x},\tilde{b}$ gives the same componentwise bound with $s=2n$ and the realified residual. This is a certificate for that explicit realified residual check; it is not a bound for an unspecified native complex-operation DAG or a solver-stability theorem.

The complex branch uses the explicit $2n$-dimensional realification DAG, so its parameter is ${\gamma }_{2n+1}$. It does not treat one complex multiply-add as one real operation. A different residual evaluator needs its own operation-DAG analysis.

Let $N\ge 1$, $c\in {\mathbb{C}}^N$, and let $C(c)$ be the circulant of Toeplitz, circulant, and finite-index conventions for structured solvers. Put ${\lambda }_k:={\sum }_{j=0}^{N-1}{c}_j{e}^{-2\pi ikj\mathrm{/}N}=\sqrt{N}({\mathcal{F}}_{N}c)(k),0\le k<N\mathrm{.}$ Then $C(c)={\mathcal{F}}_N^{-1}\mathrm{diag}({\lambda }_0,\dots ,{\lambda }_{N-1}){\mathcal{F}}_N\mathrm{.}$ Consequently, $C(c)$ is invertible exactly when every ${\lambda }_k\mathrm{\ne }0$. For $b\in {\mathbb{C}}^N$, $C(c)x=b$ is solvable exactly when $({\mathcal{F}}_{N}b)(k)=0$ for every zero mode ${\lambda }_k=0$; on nonzero modes every solution has $({\mathcal{F}}_{N}x)(k)=({\mathcal{F}}_{N}b)(k)\mathrm{/}{\lambda }_k$, while its zero-mode coordinates are free.

The square records the identity ${\mathcal{F}}_{N}C(c)=\mathrm{diag}(\lambda ){\mathcal{F}}_N$; the two paths send each input vector to the same Fourier-coordinate vector. The proof below establishes this identity from the finite convolution sums.

diagram failed to render:
\begin{tikzcd}
\mathbb C^N \arrow[r,"C(c)"] \arrow[d,"\mathcal F_N"'] & \mathbb C^N \arrow[d,"\mathcal F_N"] \\
\mathbb C^N \arrow[r,"\operatorname{diag}(\lambda)"] & \mathbb C^N
\end{tikzcd}

Let $m,n\ge 1$, let $T\in {\mathbb{C}}^{m\times n}$ be Toeplitz with $T_{ij}=t_{i-j}$ for $0\le i<m$, $0\le j<n$, and let $x\in {\mathbb{C}}^n$. Set $L=m+n-1$ and let $N$ be the least power of two with $N\ge L$. Define $c\in {\mathbb{C}}^N$ by $c_{\mathrm{\ell }}=t_{\mathrm{\ell }}(0\le \mathrm{\ell }<m),c_{N-d}=t_{-d}(1\le d<n),$ and set every other coordinate to zero. Let $\tilde{x}=(x_0,\dots ,x_{n-1},0,\dots ,0{)}^T\in {\mathbb{C}}^N$. In exact arithmetic, the first $m$ coordinates of the cyclic convolution $c{\ast }_N\tilde{x}$ equal $Tx$.

The convolution can be computed with three exact length-$N$ DFT calls, pointwise products, and scalar scaling. Its arithmetic count is $O(N{\log }_{2}N)$ for $N\ge 2$ and constant when $N=1$. This is an exact-arithmetic count for the formal exact-constant FFT model; a finite-precision execution is not claimed to be exact.

The computed convolution has a separate conditional error theorem. Under its successful-operation, unit-roundoff, and twiddle-certificate hypotheses, Computed FFT convolution error separates transform, product, inverse, and final-scale terms bounds the cyclic-convolution error; restricting to the first $m$ coordinates cannot increase its Euclidean norm. No finite-precision exactness is asserted here.

Let $n\ge 1$ and let $T_n=(t_{i-j}{)}_{0\le i,j<n}$ be Hermitian positive definite, so $t_{-r}=\overline{t_r}$. For $0\le k\le n$, set $T_k=(t_{i-j}{)}_{0\le i,j<k}$, with $T_0$ empty, and for $0\le k<n$ set $g_k=(t_1,\dots ,t_k{)}^T$, with $g_0$ empty. Define ${\varphi }_0$ empty and ${\nu }_0=t_0$. For each $0\le k\le n-2$, let $J_k$ reverse $k$ coordinates, put $c_k=J_k\overline{g_k}$, and set ${\alpha }_{k+1}:=\frac{t_{k+1}-c_k^{\ast }{\varphi }_k}{{\nu }_k},{\varphi }_{k+1}:=\left(\begin{array}{c}{\varphi }_k-J_k\overline{{\varphi }_k}{\alpha }_{k+1}\\ {\alpha }_{k+1}\end{array}\right),{\nu }_{k+1}:={\nu }_k(1-\mathrm{\mid }{\alpha }_{k+1}{\mathrm{\mid }}^2)\mathrm{.}$ Every denominator used is strictly positive; $T_k{\varphi }_k=g_k$ and ${\nu }_k=t_0-g_k^{\ast }{\varphi }_k>0$ for $0\le k<n$. Computing through order $n-1$ takes $O(n^2)$ complex arithmetic operations and $O(n)$ working storage. This is the complex matrix-bordering extension; the real ACVF forecast procedure remains Durbin–Levinson forecast recursion.

For $A\in {\mathbb{C}}^{n\times n}$, $n\ge 1$, let $Z=Z_n$ be the lower shift with $Z_{i+1,i}=1$ for $0\le i<n-1$ and all other entries zero. Its lower-shift (Stein) displacement is ${\mathrm{\nabla }}_Z(A):=A-ZA{Z}^{\ast }\mathrm{.}$ The displacement rank of $A$ with respect to this fixed shift is $\mathrm{rank}({\mathrm{\nabla }}_Z(A))$. A pair of matrices $G,H\in {\mathbb{C}}^{n\times r}$ is a generator of width $r$ when ${\mathrm{\nabla }}_Z(A)=G{H}^{\ast }\mathrm{.}$ Such a factorization is a representation of the displacement, not a factorization of $A$ itself. The term Toeplitz-like refers to matrices whose displacement rank is small relative to $n$; the rank bound $r$ must be specified rather than unrestricted. For a family of growing dimensions, a uniform bound means that $r$ is independent of $n$. This definition makes no conditioning, pivot-growth, or numerical-stability assertion.

Input: An integer $n\ge 1$, a Hermitian Toeplitz lag list $(t_{-(n-1)},\dots ,t_{n-1})$ satisfying $t_{-r}=\overline{t_r}$, and a right-hand side $b\in {\mathbb{C}}^n$. All arithmetic in this definition is exact. The real covariance forecast procedure remains Durbin–Levinson forecast recursion; this procedure solves a general Hermitian matrix equation with arbitrary $b$.

Output: Either a vector $x\in {\mathbb{C}}^n$ or one of the statuses invalid-input, nonpositive-pivot(k), or zero-gohberg-semencul-denominator.

1. Reject invalid-input if the lag list or right-hand side has the wrong length, or if the Hermitian lag equalities fail. Set $T_{ij}=t_{i-j}$.
2. If $t_0\le 0$, return nonpositive-pivot(0). Initialize ${\varphi }_0$ empty and ${\nu }_0=t_0$. For $k=0,\dots ,n-2$, form $g_k=(t_1,\dots ,t_k{)}^T$, $c_k=J_k\overline{g_k}$, and ${\alpha }_{k+1}=(t_{k+1}-c_k^{\ast }{\varphi }_k)\mathrm{/}{\nu }_k,{\varphi }_{k+1}=({\varphi }_k-J_k\overline{{\varphi }_k}{\alpha }_{k+1},{\alpha }_{k+1}{)}^T,{\nu }_{k+1}={\nu }_k(1-\mathrm{\mid }{\alpha }_{k+1}{\mathrm{\mid }}^2)\mathrm{.}$ If a required ${\nu }_k$ is not positive, return nonpositive-pivot(k) before dividing; if a newly computed ${\nu }_{k+1}\le 0$, return nonpositive-pivot(k+1).
3. With $\varphi ={\varphi }_{n-1}$ and $\nu ={\nu }_{n-1}$, form $v=(-J_{n-1}\bar{\varphi }\mathrm{/}\nu ,1\mathrm{/}\nu {)}^T,u=J_n\bar{v}\mathrm{.}$ If $u_0=0$, return zero-gohberg-semencul-denominator.
4. Do not form dense triangular matrices. For any $a,y\in {\mathbb{C}}^n$, compute $(L(a)y{)}_i={\sum }_{j=0}^i{a}_{i-j}{y}_j,(L(a{)}^{\ast }y{)}_i={\sum }_{j=i}^{n-1}\overline{a_{j-i}}{y}_j\mathrm{.}$ Return $x=(L(u)L(u{)}^{\ast }b-L(Z_{n}v)L(Z_{n}v{)}^{\ast }b)\mathrm{/}{u}_0,$ evaluating the four triangular products successively with vector workspace.

For a Hermitian positive-definite input, Complex Hermitian Toeplitz bordering yields reflection coefficients and positive pivots gives positive reflection pivots, so $u_0=1\mathrm{/}{\nu }_{n-1}>0$. The formula is the complex Hermitian extension; it makes no finite-precision stability claim.

Let $T=(t_{i-j}{)}_{0\le i,j<n}$ be any square Toeplitz matrix, and let $Z=Z_n$ be the lower shift of Lower-shift displacement rank and generator pairs. With $e_0=(1,0,\dots ,0{)}^T$ and $t_0=T_{00}$, $T-ZT{Z}^{\ast }=e_0(e_0^{\ast }T-t_0{e}_0^{\ast })+(T{e}_0)e_0^{\ast }\mathrm{.}$ Consequently, its displacement rank is at most two. A width-two generator is $G=[e_{0}T{e}_0],H^{\ast }=\left[\begin{array}{c}{e}_0^{\ast }T-t_0{e}_0^{\ast }\\ e_0^{\ast }\end{array}\right]\mathrm{.}$

Let $n\ge 2$, let $A=\left(\begin{array}{cc}d& r^{\ast }\\ \mathrm{\ell }& M\end{array}\right)\in {\mathbb{C}}^{n\times n}$ with $d\mathrm{\ne }0$, and let $s\ge 1$ be an integer. Let $Z_n=\left(\begin{array}{cc}0& 0\\ e_0& Z_{n-1}\end{array}\right)$ be the lower shift partitioned after its first row and column. Suppose $A-Z_{n}A{Z}_n^{\ast }=G{H}^{\ast },G=\left(\begin{array}{c}\widehat{G}\\ G_b\end{array}\right),H=\left(\begin{array}{c}\widehat{H}\\ H_b\end{array}\right),$ where $G,H\in {\mathbb{C}}^{n\times s}$, $\widehat{G},\widehat{H}\in {\mathbb{C}}^{1\times s}$, and $G_b,H_b\in {\mathbb{C}}^{(n-1)\times s}$. Let $S:=M-\mathrm{\ell }{d}^{-1}{r}^{\ast }$ be the scalar-pivot Schur complement. Define $p:=e_0+(Z_{n-1}\mathrm{\ell }-\mathrm{\ell })\mathrm{/}d,q:=e_0+(Z_{n-1}r-r)\mathrm{/}\bar{d},G^{\mathrm{\prime }}:=G_b+p\widehat{G},H^{\mathrm{\prime }}:=H_b+q\widehat{H}\mathrm{.}$ Then $S-Z_{n-1}S{Z}_{n-1}^{\ast }=G^{\mathrm{\prime }}{H}^{\mathrm{\prime }\ast }\mathrm{.}$ Thus the updated displacement rank is at most $s$, although it can decrease. This is an exact identity conditional on the nonzero pivot $d$; it makes no conditioning or growth claim.