Under what conditions is the least-squares approximation bounded with the same Lipschitz gradient constants?

Question

Let $f(x):\mathbb{R}^K\Longrightarrow \mathbb{R}^L$ denote a multivariate continuously differentiable function. All the partial derivatives of $f$ (all its Jacobian elements) are bounded from above and below $C_\text{min}\leq \nabla f(x) \leq C_\text{max}$ for any $x\in\mathbb{R}^k$ ($f$ is bi-lipschitz). Let $X\in\mathbb{R}^{N\times K}$ denote an input matrix, such that its rows represent the input points ($N$ points) and the columns the features ($K$ features). We run the following least squares in order to find the linear approximation of $f$ on $X$ with Jacobian $\theta^*\in\mathbb{R}^{L\times K}$ and intercept $b^*\in\mathbb{R}^{L}$: \begin{equation} \theta^*,b^*=\operatorname{argmin}_{\theta, b} \frac{1}{N}\sum_{i=1}^N (X\theta +b - f(X_i))^2. \end{equation} The approximated Jacobian $\theta^*$ solution is given by: $\theta^*=((X-\bar{X})^T(X-\bar{X}))^{-1}(X-\bar{X})^T(f(X)-\bar{f}(X))$, where $\bar{X}, \bar{f}(X)$ are the means over the points. Now, the question is, under what conditions (preferably tight) on the data $X$ and on the function $f$, can it be assured that all the elements of $\theta^*$ are also bounded with the same Jacobian bounds, possibly with error $\epsilon$: $\forall i,j$ $\theta^*_{i,j}\in [C_\text{min}-\epsilon,C_\text{max}+\epsilon]$. Note that it is not generally true that $\theta^*_{i,j}\in[C_\text{min},C_\text{max}]$ when $K>1$.