Transition between field representation
Solution 1:
I try to describe the passage between $GF(16)$ and $GF(4^2)$. The same principles can be applied in going from $GF(256)$ to $GF(16^2)$ and back.
Let us begin with $GF(4)$. It's elements are $0,1$, an element called $\alpha$ (aka the generator of the field) and $\alpha+1$. The arithmetic follows from the rule that $\alpha$ is a zero of the polynomial $x^2+x+1$, in other words $\alpha^2=\alpha+1$. In a program we represent the element $a_1\alpha+a_0$ with the pair of bits $(a_1,a_0)$. So the individual bits are viewed as coefficients of a (at most linear) polynomial that we imagine to be evaluated at $\alpha$.
Similarly elements of $GF(16)$ are represented as sequences of bits $(a_3,a_2,a_1,a_0)$. This sequence represents the element $a_3\beta^3+a_2\beta^2+a_1\beta+a_0$, where $\beta$ is the (possibly, but hopefully not mysterious) generator. This time the arithmetic follows from the fact that $\beta$ is a zero of the polynomial $x^4+x+1$, in other words $\beta^4=\beta+1$. For example, we then see that $$ \beta^5=\beta\cdot\beta^4=\beta(\beta+1)=\beta^2+\beta, $$ and $$ \beta^{10}=(\beta^5)^2=(\beta^2+\beta)^2=\beta^4+2\beta^2\beta+\beta^2 =\beta^4+\beta^2=\beta^2+\beta+1, $$ using the fact that in all fields of characteristic two we have $2=1+1=0$. This piece of information is systematically used in a lot of software and hardware in the sense that field addition is implemented as bitwise XOR of the sequences representing the elements of the field.
We make the important observation that $\beta^2+\beta$ is a root of the same equation $x^2=x+1$ as $\alpha$. This allows us to identify $GF(4)$ as a subset of $GF(16)$. We can choose to equate $\alpha=\beta^2+\beta$.
Transformation from $GF(16)$ to $GF(4^2)$: Consider an arbitrary element $z=a_3\beta^3+a_2\beta^2+a_1\beta+a_0$ of $GF(16)$. We can get rid of the high powers of $\beta$ here using the equation $\beta^2+\beta=\alpha$, or rather the equivalent form $\beta^2=\beta+\alpha$, and its consequence $\beta^3=\beta\cdot\beta^2=\beta(\beta+\alpha)=\beta^2+\alpha\beta=\alpha\beta+\beta+\alpha.$ Substituting these allows us to rewrite $$ \begin{aligned} z&=a_3(\alpha\beta+\beta+\alpha)+a_2(\beta+\alpha)+a_1\beta+a_0 =(a_3\alpha+[a_3+a_2+a_1])\beta+([a_3+a_2]\alpha+a_0)\\ &=a_h\beta+a_\ell, \end{aligned} $$ where $a_h=(a_3\alpha+[a_3+a_2+a_1])$ and $a_\ell=([a_3+a_2]\alpha+a_0)$ are elements of $GF(4)$. Internally we represent elements of $GF(4)$ with the sequence of coefficient of $\alpha^j,j=0,1$, e.g. $a_h=a_{h1}\alpha+a_{h0}$, so here $$ a_{h1}=a_3,\quad a_{h0}=a_3+a_2+a_1,\quad a_{\ell1}=a_3+a_2,\quad a_{\ell0}=a_0. $$
Transformation from $GF(4^2)$ to $GF(16)$: The other direction is also straightforward with the identification $\alpha=\beta^2+\beta$ in place. If we are given an element of $GF(4^2)$ in the form $$ z=a_h\beta+a_\ell, $$ where $a_h=a_{h1}\alpha+a_{h0}$ and $a_\ell=a_{\ell1}\alpha+a_{\ell0}$, then $$ \begin{aligned} z&=(a_{h1}[\beta^2+\beta]+a_{h0})\beta+(a_{\ell1}[\beta^2+\beta]+a_{\ell0})\\ &=a_{h1}\beta^3+(a_{\ell1}+a_{h1})\beta^2+(a_{h0}+a_{\ell1})\beta+a_{\ell0}\\ &=z_3\beta^3+z_2\beta^2+z_1\beta+z_0, \end{aligned} $$ where $z_3=z_{h1}$, $z_2=(a_{\ell1}+a_{h1})$, $z_1=(a_{h0}+a_{\ell1})$ and $z_0=a_{\ell0}$.
Inversion in $GF(4^2)$: Here the idea is to take advantage of the fact that inversion is easier in $GF(4)$ than it would be in $GF(16)$. Namely $1^{-1}, \alpha^{-1}=\alpha+1$ and $(\alpha+1)^{-1}=\alpha$. Equivalently $x^{-1}=x^2$ for any non-zero $x\in GF(4)$. The idea is similar to inverting a complex number $z=x+yi$ where $x,y$ are real, at least one non-zero. There we make the observation that $$ (x+yi)(x-yi)=x^2+y^2=Q(x,y), $$ where we can deduce that the real number $Q(x,y)\neq0$. As we know how to compute the reciprocal of $Q(x,y)$ we can take advantage and compute $$ \frac1{x+yi}=\frac1{Q(x,y)}(x-yi). $$ The key to the success of such calculation is (this is bread and butter in field theory, but does require a more extensive background in abstract algebra) that $i$ and $-i$ are roots of the same polynomial $x^2+1$ with real coefficients.
We can do something similar here. Remember that $\beta$ is a root of the polynomial $x^2+x+\alpha$ with coefficients in $GF(4)$. The other root of that polynomial is $\beta'=\beta+1$, because $$ \beta'^2+\beta'=(\beta+1)^2+(\beta+1)=\beta^2+1+\beta+1=\beta^2+\beta=\alpha. $$ So given an element $z=a_h\beta+a_\ell$ of $GF(4^2)$ let's calculate $$ z(a_h\beta'+a_\ell)=a_h^2\beta\beta'+a_ha_\ell(\beta+\beta')+a_\ell^2. $$ Here $\beta\beta'=\beta(\beta+1)=\beta^2+\beta=\alpha$ and $\beta+\beta'=1$, so we get $$ z(a_h\beta'+a_\ell)=Q(a_h,a_\ll), $$ where $Q(x,y)=\alpha x^2+xy+y^2$. Here we can prove that $Q(x,y)\neq0$ if at least one of $x,y$ is non-zero. If $x=0$, then $Q(x,y)=Q(0,y)=y^2\neq0$, because then $y\neq0$. OTOH, if $x\neq0$, then $$ Q(x,y)=\frac1{x^2}\left(\alpha+\frac{y}{x}+(\frac{y}{x})^2\right). $$ Here the factor in round parens is $P(y/x)$, where $P(T)=T^2+T+\alpha$. We have seen that $P(T)=0$, when $T=\beta$ or $T=\beta'$. Because $P(T)$ is a quadratic polynomial it can have at most two zeros in the field $GF(16)$. Because neither of those zeros was in the smaller field $GF(4)$ we can deduce that here $P(y/x)\neq0$, because $y/x\in GF(4)$.
Putting all this together gives us, at long last, $$ z^{-1}=\frac{1}{Q(a_h,a_\ell)}(a_h\beta'+a_\ell) =\frac{1}{Q(a_h,a_\ell)}(a_h\beta+[a_h+a_\ell]). $$ In other words, using the representation of $GF(16)$ as $GF(4^2)$ allows us to calculate the inverse in $GF(16)$ as a single inversion, a single addition, and two multiplications in $GF(4)$. I guess that may be worth the engineers' while in this setting.
What is needed to extend this to conversion between $GF(256)$ and $GF(16^2)$: AES specifies that the field $GF(256)$ be given by using a generator $\theta$ that is a zero of the polynomial $x^8+x^4+x^3+x+1$ [EDIT Earlier I had another polynomial here. Corrected.]. IOW the arithmetic of $GF(256)$ follows from the rule $$ \theta^8=\theta^4+\theta^3+\theta+1. $$ We need to identify a copy of $GF(16)$ inside that field. A little bit of tinkering allowed me to find that $$ (\theta^{16}+\theta)=\theta^2+\theta^3+\theta^4+\theta^6 $$ is a solution of $x^4+x+1=0$. So let us declare $$\beta=\theta^2+\theta^3+\theta^4+\theta^6.$$ Then we can calculate that $${\bf\{e\}}=(1,1,1,0)_2=\beta^3+\beta^2+\beta=\theta^2+\theta^3+\theta^5+\theta^6+\theta^7.$$
It is possible to prove that the polynomial $q(x)=x^2+x+{\bf\{e\}}$ has no solutions in $GF(16)$, and hence two in $GF(256)$. The simplest proof of this relies on the properties of the trace mapping.
A little bit of searching reveals that the roots of $q(x)$ in $GF(256)$ are $$ \gamma=\theta+\theta^2+\theta^3+\theta^4+\theta^5+\theta^6+\theta^7 $$ and $$\gamma+1=\gamma^{16}.$$ We can then represent elements $z\in GF(256)$ in "$GF(16^2)$" format like $$ z=a_h\gamma+a_\ell, $$ where this time $a_h,a_\ell\in GF(16)$. The formula for converting the 8-tuple of bits $(z_0,z_1,\ldots,z_7)$ to the pair $(a_h,a_\ell)\in GF(16)^2$ such that $$ z_0+z_1\theta+z_2\theta^2+\cdots+z_7\theta^7=z=a_h\gamma+a_\ell $$ is quite messy, and I haven't calculated it.
Solution 2:
Late answer, but this answer may be useful for others that find this question.
The two links in the question mention the method of sub-byte transformation also known as sub-field mapping, composite field mapping, field reduction, ... , which has been used for Galois Field inversion (1/z) to reduce gate count in hardware implementations for Reed Solomon codes going back some time around 1990, perhaps prior to that. Link to a report from professor E J Weldon Jr to a tape company I worked for back then. The mapping in this report doesn't minimize gate count, but does reduce it compared to a hardware lookup table, and there was only one instance of an inverter, unlike an S box which could have 10 or more encoders + decoders in the same chip.
https://github.com/jeffareid/finite-field/blob/master/wldnrpt.pdf
In software, unless a field is very large, lookup tables can be used for inversion. For a software based GF(2^8), a 256 byte table could be used, index by z with values 1/z in each entry of the table.
Getting back to the mapping, all Galois (finite) fields of the same size are isomorphic, but require mapping of elements between fields in order for the fields to be isomorphic in addition and multiplication: map(a + b) = map(a) + map(b) and map(a b) = map(a) map(b). Typically, to map from GF(2^8) to GF((2^4)^2) or to GF(((2^2)^2)^2), and 8 row by 8 bit matrix is used to multiply an element of GF(2^8) treated as an 8 row by 1 bit matrix. The inverse matrix is used to map back to GF(2^8).
I have yet to find an article that explains the derivation and how such a mapping matrix is created, so I explained this in the answer linked to below. Short version: the parameters (polynomials, primitive elements) for GF((2^4)^2) or GF(((2^2)^2)^2) are chosen to minimize gate count. Generally GF(2^8) polynomial is fixed, so the only variable is which of 128 possible primitive elements to use for the mapping. This is found by brute force trial and error search for a primitive element that will produce a isomorphic mapping matrix (there are ways to optimize this for larger fields, but with only 128 possible candidates, there is no point in optimizing a one time search). The matrix is based on the field polynomials and the primitive elements of the two fields. Link to a full explanation:
https://math.stackexchange.com/questions/3739707#3756361