Here's an approach for solving this type of problem using latent variable models. It's not a specific model, but a general way to formulate a model by breaking the description of the system into two parts: the relationship between individual inputs and (unobserved) individual outputs, and the relationship between individual outputs and (observed) aggregate group outputs. This gives a natural way to think about the problem that (hopefully somewhat) mirrors the data generating process, and makes assumptions explicit. Linear or nonlinear relationships can be accommodated, as well as various types of noise model. There's well-developed, general-purpose machinery for performing inference in latent variable models (mentioned below). Finally, explicitly including individual outputs in the model gives a principled way to make predictions about them. But, of course there's no free lunch--aggregating data destroys information.
General approach
The central idea is to treat the individual outputs as latent variables, since they're not directly observed.
Suppose the individual inputs are $\{x_1, \dots, x_n\}$, where each $x_i \in \mathbb{R}^d$ contains both individual and group-level features for the $i$th individual (group-level features would be duplicated across individuals). Inputs are stored on the rows of matrix $X \in \mathbb{R}^{n \times d}$. The corresponding individual outputs are represented by $y = [y_1, \dots, y_n]^T$ where $y_i \in \mathbb{R}$.
The first step is to postulate a relationship between the individual inputs and outputs, even though the individual outputs are not directly observed in the training data. This takes the form of a joint conditional distribution $p(y \mid X, \theta)$ where $\theta$ is a parameter vector. Of course, it factorizes as $\prod_{i=1}^n p(y_i \mid x_i, \theta)$ if the outputs are conditionally independent, given the inputs (e.g. if error terms are independent).
Next, we relate the unobserved individual outputs to the observed aggregate group outputs $\bar{y} = [\bar{y}_1, \dots, \bar{y}_k]^T$ (for $k$ groups). In general, this takes the form of another conditional distribution $p(\bar{y} \mid y, \phi)$, since the observed group outputs may be a noisy function of the individual outputs (with parameters $\phi$). Note that $\bar{y}$ is conditionally independent of $X$, given $y$. If group outputs are a deterministic function of the individual outputs, then $p(\bar{y} \mid y)$ takes the form of a delta function.
The joint likelihood of the individual and group outputs can then be written as:
$$p(y, \bar{y} \mid X, \theta, \phi) = p(\bar{y} \mid y, \phi) p(y \mid X, \theta)$$
Since the individual outputs are latent variables, they must be integrated out of the joint likelihood to obtain the marginal likelihood for the observed group outputs:
$$p(\bar{y} \mid X, \theta, \phi) =
\int p(\bar{y} \mid y, \phi) p(y \mid X, \theta) dy$$
If group outputs are a known, deterministic function of the individual outputs, the marginal likelihood can be written directly without having to think about this integral (and $\phi$ can be ignored).
Maximum likelihood estimation
Maximum likelihood estimation of the parameters proceeds by maximizing the marginal likelihood:
$$\theta_{ML}, \phi_{ML} \ = \
\arg \max_{\theta,\phi} \ p(\bar{y} \mid X, \theta, \phi)$$
If the above integral can be solved analytically, it's possible to directly optimize the resulting marginal likelihood (either analytically or numerically). However, the integral may be intractable, in which case the expectation maximization algorithm can be used.
The maximum likelihood parameters $\theta_{ML}$ could be studied to learn about the data generating process, or used to predict individual outputs for out-of-sample data. For example, given a new individual input $x_*$, we have the predictive distribution $p(y_* \mid x_*, \theta_{ML})$ (whose form we already chose in the first step above). Note that this distribution doesn't account for uncertainty in estimating the parameters, unlike the Bayesian version below. But, one could construct frequentist prediction intervals (e.g. by bootstrapping).
Care may be needed when making inferences about individuals based on aggregated data (e.g. see various forms of ecological fallacy). It's possible that these issues may be mitigated to some extent here, since individual inputs are known, and only the outputs are aggregated (and parameters are assumed to be common to all individuals). But, I don't want to make any strong statements about this without thinking about it more carefully.
Bayesian inference
Alternatively, we may be interested in the posterior distribution over parameters:
$$p(\theta, \phi \mid \bar{y}, X) =
\frac{1}{Z} p(\bar{y} \mid X, \theta, \phi) p(\theta, \phi)$$
where $Z$ is a normalizing constant. Note that this is based on the marginal likelihood, as above. It also requires that we specify a prior distribution over parameters $p(\theta, \phi)$. In some cases, it may be possible to find a closed form expression for the posterior. This requires an analytical solution to the integral in the marginal likelihood, as well as the integral in the normalizing constant. Otherwise, the posterior can be approximated, e.g. by sampling (as in MCMC) or variational methods.
Given a new individual input $x_*$, we can make predictions about the output $y_*$ using the posterior predictive distribution. This is obtained by averaging the predictive distributions for each possible choice of parameters, weighted by the posterior probability of these parameters given the training data:
$$p(y_* \mid x_*, X, \bar{y}) =
\iint p(y_* \mid x_*, \theta) p(\theta, \phi \mid \bar{y}, X) d\theta d\phi$$
As above, approximations may be necessary.
Example
Here's an example showing how to apply the above approach with a simple, linear model, similar to that described in the question. One could naturally apply the same techniques using nonlinear functions, more complicated noise models, etc.
Generating individual outputs
Let's suppose the unobserved individual outputs are generated as a linear function of the inputs, plus i.i.d. Gaussian noise. Assume the inputs include a constant feature (i.e. $X$ contains a column of ones), so we don't need to worry about an extra intercept term.
$$y_i = \beta \cdot x_i + \epsilon_i \quad \quad
\epsilon_i \sim \mathcal{N}(0, \sigma^2)$$
Therefore, $y = [y_1, \dots, y_n]^T$ has a Gaussian conditional distribution:
$$p(y \mid X, \beta, \sigma^2) = \mathcal{N}(y \mid X \beta, \sigma^2 I)$$
Generating aggregate group outputs
Suppose there are $k$ non-overlapping groups, and the $i$th group contains $n_i$ known points. For simplicity, assume we observe the mean output for each group:
$$\bar{y} = W y$$
where $W$ is a $k \times n$ weight matrix that performs averaging over individuals in each group. $W_{ij} = \frac{1}{n_i}$ if group $i$ contains point $j$, otherwise $0$. Alternatively, we might have assumed that observed group outputs are contaminated with additional noise (which would lead to a different expression for the marginal likelihood below).
Marginal likelihood
Note that $\bar{y}$ is a deterministic, linear transformation $y$, and $y$ has a Gaussian conditional distribution. Therefore, the conditional distribution of $\bar{y}$ (i.e. the marginal likelihood) is also Gaussian, with mean $W X \beta$ and covariance matrix $\sigma^2 W W^T$. Note that $W W^T = \text{diag}(\frac{1}{n_1}, \dots, \frac{1}{n_k})$, which follows from the structure of $W$ above. Let $\bar{X} = W X$ be a matrix whose $i$th row contains the mean of the inputs in the $i$th group. Then, the marginal likelihood can be written as:
$$p(\bar{y} \mid X, \beta, \sigma^2) = \mathcal{N} \left(
\bar{y} \ \Big| \
\bar{X} \beta, \
\sigma^2 \text{diag} \big( \frac{1}{n_1}, \dots, \frac{1}{n_k} \big)
\right)$$
The covariance matrix is diagonal, so the observed outputs are conditionally independent. But, they're not identically distributed; the variances are scaled by the reciprocal of the number of points in each group. This reflects the fact that larger groups average out the noise to a greater extent.
Maximum likelihood estimation
Maximizing the likelihood is equivalent to minimizing the following loss function, which was obtained by writing out the negative log marginal likelihood and then discarding constant terms:
$$\mathcal{L}(\beta, \sigma^2) =
k \log(\sigma^2)
+ \frac{1}{\sigma^2}
(\bar{y} - \bar{X} \beta)^T N (\bar{y} - \bar{X} \beta)$$
where $N = \text{diag}(n_1, \dots, n_k)$. From the loss function, it can be seen that the maximum likelihood weights $\beta_{ML}$ are equivalent to those obtained by a form of weighted least squares. Specifically, by regressing the group-average outputs $\bar{y}$ against the group-average inputs $\bar{X}$, with each group weighted by the number of points it contains.
$$\beta_{ML} = (\bar{X}^T N \bar{X})^{-1} \bar{X}^T N \bar{y}$$
The estimated variance is given by a weighted sum of the squared residuals:
$$\sigma^2_{ML} = \frac{1}{k}
(\bar{y} - \bar{X} \beta_{ML})^T N (\bar{y} - \bar{X} \beta_{ML})$$
Prediction
Given a new input $x_*$, the conditional distribution for the corresponding individual output $y_*$ is:
$$p(y_* \mid x_*, \beta_{ML}, \sigma^2_{ML}) = \mathcal{N}(y_* \mid \beta_{ML} \cdot x_*, \sigma^2_{ML})$$
The conditional mean $\beta_{ML} \cdot x_*$ could be used as a point prediction.
References
Machine learning: A probabilistic perspective (Murphy 2012). I don't recall that it speaks specifically about aggregated data, but, it covers concepts related to latent variable models quite well.
