The helical coordinate

Cholesky decomposition

Conceptually, the simplest computational method of spectral factorization might be ``Cholesky decomposition.'' For example, the matrix of (15) could have been found by Cholesky factorization of (14). The Cholesky algorithm takes a positive-definite matrix $\bold Q$ and factors it into a triangular matrix times its transpose, say $\bold Q = \bold T\T \, \bold T$ .

It is easy to reinvent the Cholesky factorization algorithm. To do so, simply write all the components of a $3\times 3$ triangular matrix $\bold T$ and then explicitly multiply these elements times the transpose matrix $\bold T\T$ . You then find you have everything you need to recursively build the elements of $\bold T$ from the elements of $\bold Q$ . Likewise, for a $4\times 4$ matrix, etc.

The $1\times 1$ case shows that the Cholesky algorithm requires square roots. Matrix elements are not always numbers. Sometimes, matrix elements are polynomials, such as $Z$ -transforms. To avoid square roots, there is a variation of the Cholesky method. In this variation, we factor $\bold Q$ into $\bold Q=\bold T\T\bold D\bold T$ , where $\bold D$ is a diagonal matrix.

Once a matrix has been factored into upper and lower triangles, solving simultaneous equations is simply a matter of two back substitutions: (We looked at a special case of back substitution with Equation ().) For example, we often encounter simultaneous equations of the form $\bold B\T \, \bold B\bold m=\bold B\T\, \bold d$ . Suppose the positive-definite matrix $\bold B\T\, \bold B$ has been factored into triangle form $\bold T\T\, \bold T\bold m=\bold B\T\, \bold d$ . To find $\bold m$ , first backsolve $\bold T\T\, \bold x=\bold B\T\, \bold d$ for the vector $\bold x$ . Then, we backsolve $\bold T\bold m=\bold x$ . When $\bold T$ happens to be a band matrix, then the first back substitution is filtering down a helix, and the second is filtering back up it. Polynomial division is a special case of back substitution.

Poisson's equation $\nabla^2 \bold p = -\bold q$ requires boundary conditions, that we can honor when we filter starting from both ends. We cannot simply solve Poisson's equation as an initial-value problem. We could insert the Laplace operator into the polynomial division program, but the solution would diverge.

The helical coordinate