In a program I'm writing, I have a sub-routine that does some vectorized linear operations (specifically differentiation). Say for convenience I have defined the following inline function, which approximates a partial derivative from a 3D array:

$$ F = @(X)(X(3:\text{end},:,:)-X(1:\text{end}-2,:,:))/(2h) $$

This is a linear transformation which would clearly have a matrix representation. Is there a simple way to derive the index form of its transpose, $F^{T}$? I would like to avoid building matrices if possible because the array formulation is much cleaner (and I've already implemented it...). The "interior" grid points seem obvious - the transpose is also a difference operator. But is there an easy way to see what the boundary terms will be?



I don't know of a general way to do this, but it is certainly not generally possible to write the transpose of a given linear operator with such a simple vectorized expression. Consider your example, but applied to a simple 1d vector x instead of a 3d array, I will call this function $G(x)$:
G = @(x) ( x(3:end)-x(1:end-2) )/(2*h)

  • $G$ operates on vectors of length $n$ and returns vectors of length $n-2$
  • If I write $G$ as a matrix it is $$ G_{(n-2)\times n} = \dfrac{1}{2h}\left[\begin{array}{c} -1 & 0 & 1 & 0 & \dots & 0\\ 0 & -1 & 0 & 1 & \dots & 0\\ & & \ddots & \ddots & \ddots & \\ 0 & \dots & 0 & -1 & 0 & 1\\ \end{array}\right] $$
  • The transpose of $G$ is $$ G^\intercal_{n\times(n-2)} = \dfrac{1}{2h}\left[\begin{array}{c} -1 & 0 & & 0\\ 0 & -1 & & \vdots\\ 1 & 0 & \ddots & 0 \\ 0 & 1 & \ddots & -1 \\ \vdots & \vdots & \ddots & 0\\ 0 & 0 & & 1 \end{array}\right] $$
  • There are a couple obvious ways to write a function for $G^\intercal x$, but none of them quite as simple as the expression for $G$. For example:
    Gtranspose = @(x) [-x(1:2); x(1:end-2)-x(3:end); x(end-1:end)]/(2*h);

A primary reason for the difference is the fact that G is non-square. What you requested can be done very easily if the operator is square and symmetric or anti-symmetric. As you probably know, a symmetric operator is its own transpose and the negative of an anti-symmetric operator is its transpose. For example, your derivative operator $F$ is very nearly anti-symmetric. If you consider periodic boundary conditions then the operator is symmetric and can be written as

H = @(X) [ X(end,:,:) - X(2,:,:); ...
           X(3:end,:,:) - X(1:end-2,:,:); ...
           X(end-1,:,:)-X(1,:,:) ]/(2*h);

the transpose is then given by $H^\intercal = -H$ so

Htranspose = @(X) [ X(2,:,:) - X(end,:,:); ...
                    X(1:end-2,:,:) - X(3:end,:,:); ...
                    X(1,:,:) - X(end-1,:,:) ]/(2*h);

This approach should work for finite differences based on central differences which are always either symmetric (even derivatives) or anti-symmetric (odd derivatives).

  • 1
    $\begingroup$ Thanks for the response! I decided to stop being a wimp and learn how to use Kronecker products for constructing difference matrices. Sparse matrix multiply ends up being just as fast as indexing anyway, so overall a net benefit. $\endgroup$
    – icurays1
    Jan 29 '15 at 16:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.