In ordinary least squares regression, for outcome vector $y$ and design matrix $X$ (full rank), the estimated coefficient values are $\hat{\beta} = (X^TX)^{-1}X^TY$. Given a new set of covariates $X_{new}$, the predicted values are $y_{new} = X_{new}\hat{\beta}.$
For smoothing splines, 
My question: is $\hat{\theta}$ analogous to $\hat{\beta}_{OLS}$? That is, are the predicted values given by a smoothing spline $y_{new} = X_{new}\hat{\theta}$?
Smoothing splines are a basis expansion method and for all of these methods prediction goes the same way. We have a model like $$ f(x) = \sum_{j=1}^m \beta_j h_j(x) $$ where the $h_j$ are completely known and the only thing to be estimated is $\beta$ which determines the relative weights of the basis functions. For a point $x_0 \in \mathbb R^p$, we first represent it in terms of the basis functions as $$ (h_1(x_0), \dots, h_m(x_0))\in\mathbb R^m $$ and then it's like a linear regression on this new representation of $x_0$ so our prediction is $$ \langle (h_1(x_0), \dots, h_m(x_0)), \beta\rangle = \sum_{j=1}^m \beta_j h_j(x_0) = f(x_0). $$ The mapping from $\mathbb R^p$ to $\mathbb R^m$ given by $$ x_0 \mapsto (h_1(x_0), \dots, h_m(x_0)) $$ is completely known (when we choose to use a cubic regression spline, say, we've implicitly chosen our basis functions [up to changes of basis]), so once $\hat\beta$ is obtained from the training data we have everything we need to predict for any new point $x_0$.
In your case the form of $\hat\theta$ represents the fact that a smoothing spline leads to a generalized ridge regression, so $\theta$ is estimated with shrinkage, but for prediction we just do $$ x_0 \mapsto \sum_j \hat\theta_j N_j(x_0) $$ which is exactly because this is a basis expansion method.
