Consider a finite field $F$ and suppose we have a system of equations
$$h_1(\alpha,\beta)=0,h_2(\alpha,\beta)=0,...,h_t(\alpha,\beta)=0$$
where $\alpha=(\alpha_1,...,\alpha_s)$ and $\beta=(\beta_1,...,\beta_k)$ with the property that: For all $i\in\{1,...,t\}:$ the function $h_i(\alpha,\beta)$ is linear in the variables $\beta$ for all $\alpha\in F^{s}$ and linear in the variables $\alpha$ for all $\beta\in F^{k}$ (i.e. if we fix the variables $\alpha$ or the variables $\beta$ we get a linear system).
This property implies that we can write each $h_i$ in the form
$$h_i(\alpha,\beta)=\sum_{j=1}^{k} h_i^j(\alpha)\beta_j$$
where $h_i^j(\alpha)$ is a linear function. My question is: Can we build the matrix $t\times k$ where the $(i,j)$ entry is the function $h_i^j(\alpha)$, and then apply Gaussian Reduction (or any other method) to solve matrix in terms of the variables $\alpha$ (and the free variables from $\beta$)?
I will write down a little example to make things clearer. Suppose we have the functions $h_1(\alpha_1,\alpha_2,\beta_1,\beta_2)=2\alpha_1\beta_2+\alpha_1\beta_1$ and $h_2(\alpha_1,\alpha_2,\beta_1,\beta_2)=\alpha_2\beta_2-\alpha_1\beta_2+\alpha_1\beta_1$, then we can write each of them as
$$h_1=h_1^1(\alpha_1,\alpha_2)\beta_1+h_1^2(\alpha_1,\alpha_2)\beta_2=\alpha_1\beta_1+2\alpha_1\beta_2$$ $$h_2=h_2^1(\alpha_1,\alpha_2)\beta_1+h_2^2(\alpha_1,\alpha_2)\beta_2=\alpha_1\beta_1+(\alpha_2-\alpha_1)\alpha_1\beta_2$$
We can build the matrix
$$\left[\begin{array}{cc} \alpha_1 & 2\alpha_1 \\ \alpha_1 & \alpha_2-\alpha_1\end{array}\right]$$
and apply Gaussian reduction with the variables $\alpha$, for instance if we sum $-1$ times the first row to the second one, we obtain:
$$\left[\begin{array}{cc} \alpha_1 & 2\alpha_1 \\ 0 & \alpha_2-3\alpha_1\end{array}\right]$$
and now we make reverse substitution to obtain the solution in terms of the $\alpha$ (and the possible free variables from $\beta$, but we didn't have any in this example).
One problem seemed to be the fact that we can't multiply by inverses since we don't know a priori who wil be zero and who not, but I think that the Smith Canonical Form can take out that issue. I will precise my main question a little bit: Is this process affordable in practice, e.g., for the following set of values?: $n=55$, $F=GF(7^{n})$, $t\approx n^2$, $k\approx n^2$, $s\approx n$.
SOME BACKGROUND FOR THE PROBLEM
We are developing a Public Key Cryptosystem based on multivariate polynomials over a finite field. One of the steps of the cryptosystem requires solving a system just like the one I posted. We can choose random values for the variables $\alpha$, and what we're doing is assigning those values BEFORE reducing the matrix and after that apply Gaussian Elimination to the pure linear system that results.
We thought that maybe we can obtain a solution in terms of the $\alpha$ and $\beta$ and AFTER that just evaluate the random values of $\alpha$ in that solution to obtain a solution in terms of the free variables of $\beta$ (we can even assign random values to some of the $\beta$ variables in order to limit a little bit the range). Note that this process would be more efficient since reduction is performed only ONE TIME so it doesn't matter is that process takes long time, anyway, we're concerned with the size of the possible solution that we find.
Thanks a lot for your time. Any doubt or any suggestion, I'll be happy to answer it or hear it.