How to factorise a covariance matrix with one observation per row-column combination?

Question

I have $N$ correlated random variables. I assume that these random variables are given by the following expression:

$ \tilde{x}_i = \alpha_i + \beta_i \cdot \tilde{m} + \gamma_i \cdot \tilde{\varepsilon_i}, $

where $\tilde{m}$ is a "global" random variable and $\tilde{\varepsilon_i}$ are "variable specific" random variables (as can be seen from absence and presence of the index $i$, respectively). The mean and sigma of both $\tilde{m}$ and $\tilde{\varepsilon_i}$ are assumed to be zero and one, respectively. The $\tilde{\varepsilon_i}$ are also assumed to be independent. As a consequence, the covariance matrix should be given by the following expression:

$ C_{ij} = \beta_i \cdot \beta_j + \delta_{ij} \cdot \gamma_i \cdot \gamma_j, $

where $\delta_{ij}$ is Kronecker delta.

Now I say that each random variable comes with one number (feature $f_i$) that determines values of $\alpha_i$, $\beta_i$ and $\gamma_i$:

$ \alpha_i = \alpha (f_i), $

$ \beta_i = \beta (f_i), $

$ \gamma_i = \gamma (f_i), $

where $\alpha$, $\beta$ and $\gamma$ are some "universal" functions (the same for all N random variables).

Using the available observations of $x_i$ I can calculate the covariance matrix $C_{ij}$ and try to find such functions $\beta$ and $\gamma$ that approximate it well:

$ C_{ij} = C(f_i, f_j) = \beta(f_i) \cdot \beta(f_j) + \delta_{ij} \cdot \gamma(f_i) \cdot \gamma(f_j). $

So far no problems. The problem comes from the fact that features $f_i$ are not constants as well as the number of random variables.

For example, on the first time step I might have 3 random variables with the following values of features: $f_1 = 1.3, f_2 = 4.5, f_3 = 0.3$ and I also have the corresponding observations of the random variables: $x_1 = 1.0, x_2 = -0.5, x_3 = 4.0$. On the second step I might have 5 random variables coming with some new 5 values of features $f_i$ and 5 new observations $x_i$. How can I find functions $\beta(f)$ and $\gamma(f)$ in this case? Or, in other words, I can assume one pair of functions ($\beta_1(f)$, $\gamma_1(f)$) and another pair ($\beta_2(f)$, $\gamma_2(f)$). How can I determine which pair of functions approximate my data set better?

ADDED (to cover questions from the comments):

1. What is the difference between factor analysis and my problem? In the factor analysis we have a (covariance) matrix that we want to factorise. In my case I do not have a matrix. I would have a covariance matrix if I have a constant number of random variables and if the statistical properties of these variables (i.e. correlation between them) is constant.
2. What do I mean by a "pair of functions". I pair of functions is my hypothesis about how $\beta$ and $\gamma$ depend on feature $f$. Given a set of observations, I would like to check what hypothesis is more plausible (accurate).

Once again, my set up is as follow:

1. On each time step $t$ I have $n_t$ observations ($n_t$ random numbers): $y_1, y_2, \dots , y_{t_{n}}$
2. On each time step $t$, for each random number, I have a corresponding feature: $f_1, f_2, \dots , f_{t_{n}}$
3. I assume that $\beta$ and $\gamma$ are functions of features and I want to find out what functions describe my data in the best way.

What can also say, that my random variables instead of being indexed by an integer $i$ are "indexed" by a real valued feature $f$.

ADDED 2:

Here is an example of my data set:

   time  feature    y
0     1      1.0 -4.0
1     1     -0.5  2.0
2     1     -3.7  3.2
3     2      2.2  5.6
4     2      1.3  0.3
5     2      0.2  0.7
6     2     -4.5  2.2
7     3      7.2  4.5
8     3      0.3  5.9

Answer 1

If your model is like

$$\tilde{x}_{it} = \alpha_i + \beta_i \cdot \tilde{m}_{t} + \gamma_i \cdot \tilde{\varepsilon}_{it}$$

Where $\tilde{m}_{j}, \tilde{\varepsilon}_{it} \sim \mathcal{N}(0,1)$ then we can rewrite is as a multivariate normal distribution

$$\textbf{x} \sim N(\boldsymbol{\alpha}, \boldsymbol{\Sigma})$$

• Where $\textbf{x}$ is the vector of all observations $\lbrace x_{it} \rbrace$.

  For instance the vector $\lbrace x_{1,1}, x_{2,1}, x_{3,1}, x_{1,2}, x_{2,2}, x_{3,2}, x_{4,2}, x_{5,2} \rbrace$ corresponds to measurements three measurements in the first time step and five measurements in the second time step. The index $i$ is repeated, and so the parameter $\beta_i$ will be the same for all these measurements at different times $j$ but with the same $i$ (I am not sure whether this is what you want?).

• Where $\boldsymbol{\alpha}$ is the vector of the corresponding means.

• Where $\boldsymbol{\Sigma}$ is the covariance matrix which will have a block form

  $$\Sigma = \begin{bmatrix} C_{1ij} & 0 & \dots &0 \\ 0 & C_{2ij} & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & C_{nij} \\ \end{bmatrix}$$

  with $n$ blocks equal to the number of time steps and each block $C_{tij}$ is like your original $C_{ij}$

This is similar to mixed-effects models explained here: Intuition about parameter estimation in mixed models (variance parameters vs. conditional modes) A code example to manually construct those blocks (instead of using build in mixed model functions) is here: https://stats.stackexchange.com/a/337348

So for given features $f_{it}$ (one for each $x_{it}$?) and a given model to compute the $\alpha_i,\beta_i,\gamma_i$ the model is fully specified and this allows you to compute the likelihood and make comparisons of models based on likelihood. Or if $\tilde{m}_{j}, \tilde{\varepsilon}_{ij}$ are not really normal distributed, then the covariance matrix still holds and you may see it as an approximation of the true likelihood, the result is a quasi likelihood.

(Or potentially you want to optimize the parameters of the model and optimize the likelihood? I am not sure whether that is what you want because you explicitly ask for comparing two models. Doing that might be possible but it is not easy to fit a non-linear mixed model where the variance also depends on the mean. You could try to just put it into some optimizer, but maybe, depending on the problem, simplifications can be made to made the convergence easier. Finding those simplifications is a bit of an art and there is not a general straigtforward method.)