The scalar target $z$ is modeled as $$f(x,y) = \underline c^T \underline b, \qquad \underline b=\begin{bmatrix} 1 \\ x \\ ln(y+d) \\ x \cdot ln(y+d) \end{bmatrix},$$ with unknown parameter $d$ and unknown coefficients $\underline c \in \mathbb R^4$. A small set $S=\{ (x_i, y_i, z_i \}_{i=1}^6$ is available to fit $\underline c$. A larger dataset L ($\lvert L \rvert = 378$) is reserved exclusively for validation and must not be used to estimate $d$ or $\underline c$. The objective is to choose $d$ and $\underline c$ so that a chosen metric (e.g. MSE or RMSE) between model predctions and observed $z$ on the dataset of interest is minimized.
Linear Regression can recover $\underline c$ for any fixed $d$ $$\hat{\underline c}(d) = \displaystyle \arg \min_{\underline c} \displaystyle \sum_{i=1}^{\lvert S \rvert} (\underline b(x_i, y_i; d)^T \underline c - z_i)^2$$ but $d$ itself is unknown and with only $\lvert S \rvert = 6$ points, jointly estimating $d$ and $\underline c$ is numerically unstable. In practice it was found that manually tuning $d$ produced a better fit (lower MSE) than the automated procedure. For two candidate values $d_h$ and $d_l$ this was observed as $$MSE(d_h) < MSE(d_l)$$
Which estimation strategies, regularization choices and validation workflows do you recommend to robustly determine $d$ and $\underline c$ from only $S$ given the numerical instability caused by $\lvert S \rvert = 6$?
