General-RHS Hermitian Toeplitz solve by reflection recursion and Gohberg–Semencul

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.