DCT in 2D: Symmetric Functions

This is an extension of a previous post that sought to approximate a 1D function \(f(x), ~ x \in [0,1)\), that is symmetric about \(x = 1/2\). That is, the shape of the function above the line of symmetry is a mirror image of the function below it.

As as example, we considered the function:  \[f(x) = (1/2-x)^2 + \cos(2 \pi x).^2.\]

We sample the function at \(n\) equispaced points \(x_i = i/(n+1)\), where \(i = 0, 1, 2, ..., n\). Given the bunch of points \(\{x_i, f_i\}\), we approximate the function using a sum of discrete cosine functions, \[f(x) \approx \sum_{j=0}^{N-1} a_j v_j(x_i),\] where,  \(v_j(x) = \cos(2 \pi j x)\) is a typical orthogonal basis function, and \(N\) is the number of basis functions used.

In this post, we add a small wrinkle to the original problem.

Suppose we had a 2D function \(G(x,y) = f(x) f(y)\). This is certainly a special kind of 2D function. It has mirror symmetry about the centerline which it inherits from \(f(x)\). In addition, it has a very special symmetry along the \(x\) and \(y\) directions. For example, \(G(x,y) = G(y,x)\).

For concreteness, let us use the particular \(f(x)\) that we considered above, and build a \(G(x,y)\) out of it. It looks like:

clear all
clf

% number of collocation points is (n+1)
n  = 20;  
% collocation points for discrete cosine series
xi = [0:n]'/(n+1);
% grid is equivalent in the x and y-directions
yi = xi;

%
% Sample Function
%
f = @(x) (1/2-x).^2 + cos(2*pi*x).^2;

[X, Y] = meshgrid(xi);
G      = f(X).*f(Y);

surf(X,Y,G)
view(40,40)
xlabel('x')
ylabel('y')
zlabel('G(x,y)')

Since the "x" and "y" grid-points are at identical locations, we can assert\[G(x_i, y_j) =  \sum_{j=0}^{N-1} a_j v_j(x_i) \sum_{k=0}^{N-1} a_k v_k(y_j.)\] In the matrix language developed in the original post, this is equivalent to \[\mathbf{G} = \mathbf{V a}  (\mathbf{V a})^T = \mathbf{V A V^T},\] where \(\mathbf{A} = \mathbf{a a^T}\).

Using the property of orthogonality, we hit both sides with the appropriate factors to obtain \[\mathbf{V^T G V} =   \frac{(n+1)^2}{4} \begin{bmatrix} 2 & 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 1 & 0& 0\\ 0 & 0& 0& 1 & 0\\ 0 & 0 & 0& 0& 1\end{bmatrix} \mathbf{A}  \begin{bmatrix} 2 & 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 1 & 0& 0\\ 0 & 0& 0& 1 & 0\\ 0 & 0 & 0& 0& 1\end{bmatrix}.\]  If \(\mathbf{a} = (a_0, a_1, ..., a_{N-1})\) then this implies that the RHS is equal to, \[\frac{(n+1)^2}{4} \begin{bmatrix} 4 a_0^2 & 2 a_0 a_1 & 2 a_0 a_2 & ... & 2 a_0 a_{N-1}\\ 2 a_1 a_0 & a_1^2 & a_1 a_2 & ... & a_1 a_{N-1}\\ 2 a_2 a_0 & a_2 a_1 & a_2^2 & ...&  a_1 a_{N-1}\\  \vdots & \vdots & \vdots & \vdots & \vdots \\ 2 a_{N-1} a_0 & a_{N-1} a_1 & a_{N-1} a_2 & ... & a_{N-1}^2\end{bmatrix},\]

In the above equation, the LHS is known. The diagonal of the the computed matrix yields the vector \(\mathbf{a}\) that we seek.

Programmatically,

nmodes   = 7;

%
% Build the matrix of basis vector V = [v_0 v_1 ...  v_N-1]
%

V     = zeros(n+1, nmodes);
V(:,1) = ones(n+1,1);

for j = 2:nmodes
  V(:,j) = cos(2*pi*(j-1)*xi);
end

%
% Compose the matrix: The key step
% Replace small negative values with 0 to prevent imaginary roots
%

M      = diag(V'*G*V);
M(M<0) = 0;                
M      = sqrt(M);

a    = zeros(nmodes,1);
a(1) = M(1)/(n+1);

a(2:nmodes) = 2*M(2:nmodes)/(n+1);

%
% Compute Norm on the Same Grid
%
g       = V*a;
Gapprox = g*g';
disp('Norm of Error');
norm(Gapprox - G)/(n+1)^2

This yields the coefficients for the basis vectors which can be used to reconstruct the function. Here are contour plots of the original data on the grid, and the smooth fitted function.