NLMS algorithm for a MISO structure I am trying to implement an NLMS algorithm for a multi-input single-output(MISO) structure. 

We take a reference signal x, then we made a new set of P input signals from it as follows: x_op (k) = x(k)^p. k denotes the k-th sample of our reference signal x.

For the case where P = 1, our adaptive filter is just a vector of length N. But for P > 1, we have a matrix H (N x P) containing all the adaptive filters. (where N is the length of each filter).

Now I am confused as what would be dimenstion of the estimated output at each iteration (for each sample). For the case where P = 1, we can simply write:

y_hat = h_hat' * xk; where xk has the last N samples of x. 

and y_hat is then a scalar. What about this case when H is a matric and not a vector anymore? 

The way I understand it, y_hat for this case is no longer a scalar. But if it's not a scalar, then how I am supposed to define the error for each sample k, in order to write the adaption law equation?
I can provide you coding example, but it solves more generic problem than you posted. It is for general dynamic system that converts two inputs x(t), y(t) into single output z(t) and defined by Urysohn integral equation. Hammerstein is particular case of Urysohn and linear system is particular case of Hammerstein. The code and some explanation is here:
http://ezcodesample.com/NAF/MISO.html
http://ezcodesample.com/NAF/Wiener.html
http://ezcodesample.com/NAF/reallife2.html

And most generic case when dynamic system is an operator that converts fragment of input x(t) from t - T to t into scalar z(t) http://ezcodesample.com/NAF/kolmogorov.html