Estimate covariance matrix from a sample of linear combinations Suppose you have a sample of $p$-dimensional weight vectors $\boldsymbol{w_i} = (w_{i, 1}, w_{i, 2}, ..., w_{i, p})$ and observations $y_i$ where $n \in \{1, ..., N\}$ and $j \in \{1, ..., p\}$.
Let $y_i = \boldsymbol{w}_i^TX_i = \sum_{j=1}^p w_{i, j} X_{i, j}$ where $X_{i} = (X_{i, 1}, X_{i, 2}, ..., X_{i, p})$ is a $p$-dimensional unobserved latent vector sampled independently for each $i$ from a multivariate normal distribution $N(\mu, \Sigma)$.
Is there a way to estimate the covariance matrix $\Sigma$?
 A: Setting up the inference problem
As described in the question, $\{X_1, \dots, X_n\}$ are i.i.d. Gaussian random vectors. Each observation $Y_i$ is generated by projecting the corresponding $X_i$ onto its own weight vector $w_i$ (which is assumed to be known):
$$X_1, \dots, X_n \underset{\text{i.i.d.}}{\sim} \mathcal{N}(\mu, \Sigma)$$
$$Y_i = w_i^T X_i$$
Since we're linearly transforming Gaussian random vectors, each $Y_i$ is also Gaussian, with mean and variance given by standard rules:
$$Y_i \sim \mathcal{N}(w_i^T \mu, w_i^T \Sigma w_i)$$
Furthermore, the $Y_i$ are independent (but not identically distributed), so the joint distribution factorizes as:
$$p(y_1, \dots, y_n \mid \mu, \Sigma) =
\prod_{i=1}^n \mathcal{N}(y_i \mid w_i^T \mu, w_i^T \Sigma w_i)$$
The corresponding log likelihood function is:
$$\mathcal{L}(\mu, \Sigma) =
\sum_{i=1}^n \log \mathcal{N}(y_i \mid w_i^T \mu, w_i^T \Sigma w_i)$$
$$= -\frac{n}{2} \log (2 \pi)
- \frac{1}{2} \sum_{i=1}^n \log (w_i^T \Sigma w_i)
- \frac{1}{2} \sum_{i=1}^n \frac{(y_i - w_i^T \mu)^2}{w_i^T \Sigma w_i}$$
The parameters can estimated by maximum likelihood:
$$\max_{\mu, \Sigma} \ \mathcal{L}(\mu, \Sigma)$$
Alternatively, if you can specify a prior distribution representing a priori knowledge/assumptions about the parameters, you could use MAP estimation or compute a full Bayesian posterior distribution.
Maximum likelihood estimation
I haven't checked whether a closed-form maximum likelihood solution exists, but encourage you to look into this yourself (based on the expression for the log likelihood above). As a quick test, I numerically optimized the log likelihood for some toy examples with properties similar to those described in the comments: known/fixed mean $\mu = \vec{0}$, full rank covariance matrix $\Sigma$, and weight vectors sampled i.i.d. from $\mathcal{N}(\vec{0}, I)$. This produced seemingly unique solutions (multiple random initial guesses converged to the same solution) that were close to the true covariance matrices (given sufficiently many observations).
Existence of a unique solution
The joint distribution above is an $n$-dimensional Gaussian distribution. Let $\vec{y} = [y_1, \dots, y_n]^T$ contain the observations and $W \in \mathbb{R}^{p \times n}$ contain the weight vectors $\{w_1, \dots, w_n\}$ as its columns. Then the joint distribution can be written as:
$$p(\vec{y} \mid \mu, \Sigma) = \mathcal{N}(\vec{y} \mid W^T \mu, C)$$
where the covariance matrix $C$ is diagonal, with diagonal entries $C_{ii} = w_i^T \Sigma w_i$
Thus, the problem can be seen as finding a constrained $n$-dimensional covariance matrix $C$, which is parameterized by $\Sigma$. For a unique maximum likelihood solution to exist, every choice of $\Sigma$ must correspond to a unique $C$. This condition could be violated, for example, by linear dependence between the weight vectors, or by having too few observations.
