FACTORING THE LAPLACIAN ON A HELIX

Discretize the (x,y)-plane to an $N\times M$ array and pack the array into a vector of $N\times M$ components. Likewise pack the Laplacian operator $\partial_{xx}+\partial_{yy}$into a matrix. For a $4\times 3$ plane, that matrix is shown in equation (7).  

$$-\ \nabla^2 \eq \begin{array} {\vert cccc\vert cccc\vert ccc... ... \cdot & \cdot & 1 & -4 \\ &&&& &&&& &&& \\ \hline\end{array}$$ (7)

The two-dimensional matrix of coefficients for the Laplacian operator is shown in (7), where, on a cartesian space, h=0, and in the helix geometry, h=1. (A similar partitioned matrix arises from packing a cylindrical surface into a $4\times4$ array.) Notice that the partitioning becomes transparent for the helix, h=1.

With the partitioning thus invisible, the matrix simply represents one-dimensional convolution and we have an alternative analytical approach, one-dimensional Fourier Transform. We often need to solve sets of simultaneous equations with a matrix similar to (7). A costly method is to factor the matrix into upper and lower triangular form that can be ``backsolved'' which in this case amounts to recursive filtering.

The Fourier approach of the next section is roughly equivalent to factoring $-\nabla^2$ into upper and lower triangular matrices, but it is much faster owing to the convolutional nature of the matrix. The (negative of the) Laplacian operator is regarded as an autocorrelation  

$$\bold r = (\cdots,-1,\cdots,-1,4,-1,\cdots,-1,\cdots)$$ (8)

Using the Kolmogoroff spectral-factorization method, a wavelet is found which has this autocorrelation. We'll see how to compute this wavelet  

$$\bold a \eq (1.791,-.651,-.044,-.024,\cdots ,-.044,-.087,-.200,-.558 )$$ (9)

whose autocorrelation is (8). Displaying these signals on a helix, the 2-D autocorrelation is:

$$\bold r \eq \left[ \begin{array} {rrr} & -1 & \\ -1 & 4 & -1 \\ & -1 & \end{array}\right]$$ (10)

and the 2-D wavelet with this autocorrelation is  

$$\bold a \eq \left[ \begin{array} {rrrrrrrrr} & & & & 1.79... ...ts &-.044 & -.087 & -.200 & -.558 & & & & \end{array} \right]$$ (11)

where now I am wrapping the top row around to a second row. In the representation (11) we see the coefficients diminishing rapidly away from maximum value 1.791. In a more abstract notation, we might write

$$-\nabla^2 \eq \bold A' \bold A$$ (12)

Where the triangular matrix $\bold A$ is constructed from the row (9) as follows: Each row of the matrix $\bold A'$contains the row (9) shifted along itself so that the 1.791 lies along the main diagonal of the matrix.

Unfortunately, we see that the factored operator $\bold A$has a great number of nonzero terms, but fortunately they seem to be converging rapidly so truncation seems reasonable. We can use this feedback operator to solve Poisson's equation very rapidly. Polynomial division and its adjoint gives us $\bold p =(\bold q/\bold A)/\bold A'$which means that we have solved the PDE $\nabla^2 \bold p = -\bold q$by using polynomial division on a helix. Using the seven coefficients shown, the cost is fourteen multiplications (because we need to run both ways) per mesh point. An example is shown in Figure 7.

 

lapfac90
Figure 7 Deconvolution by a filter whose autocorrelation is the two-dimensional Laplacian operator. Amounts to solving the Poisson equation. Left is $\bold q$; Middle is $\bold q/\bold A$; Right is $(\bold q/\bold A)/\bold A'$.

We are all aware of the factorization of $\nabla^2$into a divergence dotted into a gradient, where the divergence is the adjoint of the gradient. In two-dimensional physical space, the gradient maps one field to two fields. The factorization of $-\nabla^2$ with the helix gives us an operator that maps one field to one field.

Any fact this basic should be well known in some arcane field of mathematics or theoretical physics. Meanwhile, being ignorant of any pre-existant name, I have chosen the name ``helix derivative'' or ``helical derivative'' for the operator $\bold A$.(Admittedly, the concept exists on an infinite cartesian plane without a helix, but all my codes in a finite space involve the helix, and the helix concept led me to it.)

Our construction makes $\bold A$ have the energy spectrum k_x²+k_y², so the magnitude of the Fourier transform is $\sqrt{\vert k_x^2+k_y^2\vert}$.This rotationally invariant function in the Fourier domain contrasts sharply with the nonrotationally invariant function shape in (x,y)-space. The difference must arise from the phase spectrum. The factorization (11) is nonunique in that causality associated with the helix mapping can be defined along either x or y axes; thus the operator (11) can be rotated or reflected.