We are obviously not limited to the above polynomial expansions. We could replace the various powers of \( x \) with elements of Fourier series or instead of \( x_i^j \) we could have \( \cos{(j x_i)} \) or \( \sin{(j x_i)} \), or time series or other orthogonal functions. For every set of values \( y_i,x_i \) we can then generalize the equations to

$$ \begin{align*} y_0&=\theta_0x_{00}+\theta_1x_{01}+\theta_2x_{02}+\dots+\theta_{n-1}x_{0n-1}+\epsilon_0\\ y_1&=\theta_0x_{10}+\theta_1x_{11}+\theta_2x_{12}+\dots+\theta_{n-1}x_{1n-1}+\epsilon_1\\ y_2&=\theta_0x_{20}+\theta_1x_{21}+\theta_2x_{22}+\dots+\theta_{n-1}x_{2n-1}+\epsilon_2\\ \dots & \dots \\ y_{i}&=\theta_0x_{i0}+\theta_1x_{i1}+\theta_2x_{i2}+\dots+\theta_{n-1}x_{in-1}+\epsilon_i\\ \dots & \dots \\ y_{n-1}&=\theta_0x_{n-1,0}+\theta_1x_{n-1,2}+\theta_2x_{n-1,2}+\dots+\theta_{n-1}x_{n-1,n-1}+\epsilon_{n-1}.\\ \end{align*} $$ Note that we have \( p=n \) here. The matrix is symmetric. This is generally not the case!