I'm working on implementation of a Fast surace interpolation using hierarchial basis functions (Szeliski et al) algorithm.
The idea is: given a discrete function measurements of its values (depths)
$\{p_i\} = \{(u_i, v_i, d_i)\} = \{(u_i, v_i, \hat{f}(u_i, v_i))\}$
we are looking for such a function $f(u, v)$ so that it is very close to the original function, $\hat{f}(u, v)$ on its whole domain.
Note: we are dealing with discrete values. Thus, for convenience, let's denote a set of interpolated values as $\mathrm{x}$ and a set of measurement values as $\mathrm{d}$.
To do that, we try to minimize the energy function, $E(\mathrm{x}) = E_s(\mathrm{x}) + E_d(\mathrm{x}, \mathrm{d})$.
The data cost term, $E_d(\mathrm{x}, \mathrm{d})$, determines how close the interpolated surface (or function) is in reference to the original function.
The smoothness term, $E_s(\mathrm{x})$ encodes the curvature of the surface interpolated.
Now, the source problem is discretized using finite element method, which uses the stiffness matrix, $A$. To minimize the energy function, described above, the stiffness matrix is used too. It is said, to find matrix $A$ the authors solve this equation: $A\mathrm{x}=\mathrm{b}$ by setting:
$a_{(i, j), (k, l)} = \frac{\partial^2E(\mathrm{x}) }{\partial x_{i, j}\partial x_{k, l}} \Bigg|_{\mathrm{x}=0}$
and
$b_{i, j} = - \frac{\partial E(\mathrm{x}) }{\partial x_{i, j}} \Bigg|_{\mathrm{x}=0}$
My question is: why does the element of $A$ have those strange pair indices $a_{(i, j), (k, l)}$? Or maybe what is the stiffness matrix - is it 4-dimensional? Or what it looks like?
Would appreciate any help!