Suppose that we had two channels encoded in some form in two signals according to a matrix. Recovering the original signals would be straightforward using the matrix inverse (or impossible, if the inverse did not exist). The situation is no different here, except that the matrix is not square, so you cannot use the conventional matrix inverse. You can, however, use a pseudoinverse or generalized inverse to reconstruct the original signals.
Other matrices which are not equal to an inverse may work also succeed, and may in fact sound better. I don't know the psychoacoustic situation, but it may more sense to do some transform based on the position of speakers for example. Certainly, the inverse seems like a good place to start.
I rewrote your encoder table as a 2x5 input matrix $ E $, and used Wolfram Alpha to determine the inverse, $ D = E^{-1} $. Here is the 5x2 result it provided, using $ i $ as the imaginary unit, with $ j = i $ and $ k = -i $.
$$
D = \frac{1}{18(46+6\sqrt{2})} \begin{bmatrix}
405 & -3(27 - 36\sqrt{2}) \\
-3(27 - 36\sqrt{2}) & 405 \\
6(18+27\sqrt{2}) & 6(18+27\sqrt{2}) \\
-i(27\sqrt{3} + 99\sqrt{6}) & i(63\sqrt{3} + 27\sqrt{6}) \\
-i(63\sqrt{3} + 27\sqrt{6}) & i(27\sqrt{3} + 99\sqrt{6})
\end{bmatrix}
$$
As one would expect, the decoding matrix shows left/right symmetry.
Now we can express the encoding/decoding operations as:
$$
\begin{bmatrix}
L_{encoded} \\
R_{encoded}
\end{bmatrix} = E \space{} \begin{bmatrix}
L \\
R \\
C \\
L_{rear} \\
R_{rear} \\
\end{bmatrix}
$$
$$
\begin{bmatrix}
L' \\
R' \\
C' \\
L_{rear}' \\
R_{rear}' \\
\end{bmatrix} = D \space{} \begin{bmatrix}
L_{encoded} \\
R_{encoded}
\end{bmatrix}
$$
Now you can easily obtain formulas to express the three front components as functions of the encoded left and right signals. For the two rear components, you have an equation in the time domain with a complex component.
Looking at the left rear component, and converting constants from the matrix $ D $ to numeric form:
$$
L_{rear}'[n] = -i0.2949 L_{encoded}[n] + i0.1787 R_{encoded}[n]
$$
The Hilbert transform can be used to remove the imaginary unit, because this has the same effect on the frequency domain of performing a ±90º phase shift. The Hilbert transform only needs to be applied once, since convolution ($ * $) is distributive and commutative and all other good things.
$$
L_{rear}'[n] = (-0.2949 L_{encoded}[n] + 0.1787 R_{encoded}[n]) * H[n]
$$
The operating principle, then, is that the encoding technique functions in part by encoding the left/right rear components in quadrature.
The Hilbert transform is infinitely long, so it must be truncated and expressed as either a FIR filter, or depending on your circumstances, a FFT-IFFT operation. A causal Hilbert transform will delay the rear channels, so your other filters will need to be delayed by some amount to compensate.
For example, if $ k $ is the group delay of a causal approximation of a Hilbert filter $ H_k[n] $ then you delay must delay the three front channels by $ k $ as well:
$$
C'[n] = (0.3437L_{encoded}[n] + 0.3437R_{encoded}[n]) * \delta[n-k]
$$
In practice, all of the output signals depend on the two input signals. That is, you take two input channels and produce one output channel five times. Your JACK pipeline will need to accommodate this.