With inspiration from this question, I'm wondering if biorthogonal bases for arbitrary dimensions are possible to construct with the same mechanism. I am thinking a subsampling of a factor of $N$ in each of $D$ dimensions.
In particular we want (as in the one dimensional case) perfect reconstruction. In the one dimensional case that is given by from Wavelet Tutorial source, page 19: $$\begin{align}H(z)H_i(z) + G(z)G_i(z) = 2&\\ H(-z)H_i(z) + G(-z)G_i(z) = 0&\end{align}$$ Or in matrix terms: $$\left[\begin{array}{c}H_i(z)\\G_i(z)\end{array}\right] = \frac{2}{\det(H_m)} \left[\begin{array}{r}G(-z)\\-H(-z)\end{array}\right]$$
Wikipedia the following well known relation for the one dimensional case of biorthogonality:
$$\sum_{n\in\mathbb{N}}a_n{\tilde a}_{n+2m} = 2 \delta_{m,0}$$
This basically means the convolution of $a$ and $\tilde a$ is the identity.
How to expand this to more dimensions? I will give some of my thoughts, but I hope to be challenged and get counter proposals.
Own work, or suspicions: Out of practicality let us replace tilde with an indexing vector $l$. Let us denote the sequence of filters $a(l)_n$, $l\in {{\mathbb Z}_N}^D$.
One expansion we can do is
$$a([0,0,\cdots]) * a(l) = \delta, \forall l \neq [0,\cdots,0]$$
Where this $a([0,0,\cdots])$ is a mean-preserving filter.
I think this would be compatible with separable transforms, but allow for increased freedom to design the high pass filters. Does it make sense?
Observation 2:
Note that the above are polynomial equations or systems of polynomial equations (the z-transforms $H,H_i,G,G_i$ of the filters). In one dimension these are one variable polynomials, but as the dimensionality increases, so does the number of variables of the polynomials. So if there exist similar constraints for the Z-transforms of filters for multi dimensional wavelets, they would lead to solving systems of multivariate polynomial equations.
If you want to use dyadic dilation for your wavelets, then for $n=2,$ the corresponding matrix equations are $$\begin{bmatrix} H(z_{1},z_{2})&G_{1}(z_{1},z_{2})&G_{2}(z_{1},z_{2})&G_{3}(z_{1},z_{2})\\ H(-z_{1},z_{2})&G_{1}(-z_{1},z_{2})&G_{2}(-z_{1},z_{2})&G_{3}(-z_{1},z_{2})\\ H(z_{1},-z_{2})&G_{1}(z_{1},-z_{2})&G_{2}(z_{1},-z_{2})&G_{3}(z_{1},-z_{2})\\ H(-z_{1},-z_{2})&G_{1}(-z_{1},-z_{2})&G_{2}(-z_{1},-z_{2})&G_{3}(-z_{1},-z_{2}) \end{bmatrix} \times\begin{bmatrix} \tilde{H}(z_{1},z_{2})&\tilde{H}(-z_{1},z_{2})&\tilde{H}(z_{1},-z_{2})&\tilde{H}(-z_{1},-z_{2})\\ \tilde{G}_{1}(z_{1},z_{2})&\tilde{G}_{1}(-z_{1},z_{2})&\tilde{G}_{1}(z_{1},-z_{2})&\tilde{G}_{1}(-z_{1},-z_{2})\\ \tilde{G}_{2}(z_{1},z_{2})&\tilde{G}_{2}(-z_{1},z_{2})&\tilde{G}_{2}(z_{1},-z_{2})&\tilde{G}_{2}(-z_{1},-z_{2})\\ \tilde{G}_{3}(z_{1},z_{2})&\tilde{G}_{3}(-z_{1},z_{2})&\tilde{G}_{3}(z_{1},-z_{2})&\tilde{G}_{3}(-z_{1},-z_{2}) \end{bmatrix}=4I, $$ these can be extended to higher dimensions than $2$, but usually this is done using trigonometric polynomials rather than the $z$-transform, since this lets us write things like $H(-z_{1},z_{2},z_{3},-z_{4},-z_{5})=\tau(\omega+(\pi,0,0,\pi,\pi)),$ the latter of which is easier to write in a general way.
If $\tau,\tilde{\tau}$ are multidimensional trigonometric polynomials with $\tau(0)=\tilde{\tau}(0)=1$, then the biorthogonality condition is $$\sum_{\gamma\in\{0,\pi\}^{n}}\tau(\omega+\pi)\overline{\tilde{\tau}(\omega+\pi)}=1\text{ for all }\omega\in[-\pi,\pi]^{n}.$$
As you suspected, we need to find multivariate polynomials satisfying several conditions when we go to higher dimensions. This is complicated further when we add in accuracy number/ flatness number/ vanishing moments conditions, symmetry conditions, and much harder than these, smoothness conditions on the function $\phi$ satisfying $\hat{\phi}(\omega)=\tau(\omega/2)\hat{\phi}(\omega/2)$ (the refinable function).