How can I fit a polynomial given some prior knowledge but fewer observations than coefficients? I want to fit a second degree polynomial to data from a fractional factorial design in five dimensions. Because of the cost of trials, I am initially only going to do eight runs. So I have fewer observations than the full number of first and second degree terms, including interaction terms. 
My goal is to find a local optimal point or region  in the vicinity of my trials (if there is one)  or direction of greatest increase (if there isn’t).  I am saving budget for a few (I hope) confirmatory tests.
For some, but less than half, of the total number of polynomial coefficients I have a very rough idea of their magnitude from a literature review. For others I know the sign, but not the value. 
It is my belief that I can estimate the coefficients using some kind of penalized regression. But I don’t know how to incorporate my prior knowledge, and I don’t know whether I should alter my experimental design in some way based on that knowledge (e.g. by structuring is so that the aliased interaction terms are ones where I have some prior information).
I have done a lot of kinds of regression over the years, but I have not looked at anything Bayesian since I was an undergraduate, lo these many years ago. So I’d like a response that’s a little more of a how-to than “Oh, use a Bayesian LASSO” or whatever. Extra karma points for an answer with some pointers toward an R implimentation.
 A: In $d=5$ dimensions, a second degree polynomial has $\frac{(d+1)(d+2)}{2}=21$ coefficients, which you want to estimate with 8 runs. With such a tight budget, I definitely wouldn't go for a training/test set approach. With cross-validation, you will be able to use all your runs for fitting the final model. This way, you may be able to increase the number of runs a little bit (10, maybe?). In your situation, anything helps.
Having said that, fitting a polynomial with a number of coefficient larger than the size of your design, or in general fitting an underdetermined linear model, is really the perfect job for LASSO, ridge regression or elastic net (in case you don't know whether your "ideal" model is likely to be sparse or not). The glmnet function of the R package glmnetallows you to set lower.limits and upper.limits on each coefficient: this allows you to easily impose sign constraints. For example, suppose you have four coefficients in your model: let the first one be the intercept $\beta_0$, and assume that $\beta_1$ is known to be positive, $\beta_3$ is known to be negative while we have no prior knowledge on $\beta_2$. The intercept is never shrinked in shrinkage methods, nor it can be constrained (at least with the glmnet function). Thus you can only specify constraints for $\beta_i,\ i=1,\dots 3$. You will set
lower.limits <- c(0, -Inf, -Inf)
upper.limits <- c(Inf, Inf, 0)

For the order of magnitude constraints, you could introduce lower/upper bounds as before. Another option would be to specify different penalty factors for each coefficient with the parameter penalty.factor. A larger penalty factor will shrink the coefficient towards a  smaller absolute value (or even have it dropped out of the model, for LASSO). Suppose you have again 4 coefficients: as noted before, the intercept is not included in the penalty term and thus cannot be shrinked. This leaves us with three coefficients for which you can specify penalty.factor. Since the default value is 1, then if you want the last coefficient to be larger in module than the other ones, you can speficy
penalty.factor <- c(1, 1, 0.1)

There are two drawbacks with this approach: the first one is that, even assuming that all predictors have similar mean values and standard deviations, a penalty.factor of 0.1 doesn't mean that $\beta_3$ will be 10 times as larger as the other two coefficients. This depends on the data set and also on whether you're using ridge regression (which penalizes the squared coefficients), LASSO (which penalizes the absolute values of the coefficients) or elastic net (which is a combination of both). If we drop the assumption that all predictors have similar means and standard deviations, then the relationship between penalty.factor and the coefficient size is even more complicated, because by default glmnet standardizes all predictors before applying shrinkage. The reason is explained in paragraph 6.2.1 of An Introduction to Statistical Learning  by James, Witten, Hastie and Tibshirani. This means that the penalty is applied to the coefficients of the standardized predictors, but of course glmnet returns the coefficients in the original (unstandardized) scales. Thus, unless your predictors have similar ranges initially, a penalty.factor as in my example will not necessarily return a $\beta_3$ bigger than $\beta_1$ and $\beta_2$. 
Hope this helps: if you want a more "hands-on" answer, please provide a reproducible example.
