# Learning parameters of noise and filter coefficients from data where data and noise both have Gaussian distributions

Assume $X$ and $N$ are two sets of vectors (observations) from two different normal distributions, where $X$ represents clean data and $N$ represents noise; and $A$ a projection matrix of a filter. the scenario is that our clean data was corrupted by a multiplicative noise via matrix $A$ and an additive noise of $N$. then:

$Y=A \times X + N$

where $Y$ are a set of projected vectors from $X$ using $A$, what are solutions to learn this projection matrix and $N$ from training data? Does the Gaussian assumption of $A, N$ and $X$ help to have a better estimation or guide to use a specific solution?

Here is matlab code for the training data, noise and a simple projection:

    dataVariance = .10;
noiseVariance = .05;
mixtureCenters=randn(13,1);
X=randn(13, 1000)*sqrt(dataVariance ) + repmat(mixtureCenters,1,1000);

%N and A are unknown and we want to estimate them.
N=randn(13, 1000)*sqrt(noiseVariance ) + repmat(mixtureCenters,1,1000);
A=2*eye(13);

Y=A*X+N;

for iter=1:1000
A_hat,N_hat = training(X_hat,X,Y);
end


Note: if necessary, for each estimation of $A$, an error can be calculated for an estimation of $N$ using a current $A$.

For example:

for iterate=1:1000
initiate A
estimate N using current A (N=Y-A*X)
calculate error of estimation (err=Y-A*X+N)
update A


But I would prefer not to go for gradient descent approaches.

I should clarify that the observations of $X$ and $Y$ are time independent and in $X_i$ and $Y_i$, i is not the time. They are just different observations sampled from two normal distributions.

• Have you considered trying to use a Kalman filter to estimate these parameters? – spektr Jan 22 '16 at 20:11
• @choward hi, no. not really. I thought of multivar. linear regression though. Do you think keeping the Gaussian assumptions in mind, Kalman filter can estimate $A$ and $N$? – PickleRick Jan 22 '16 at 20:30
• I think it is possible, yes. Kalman filters are based on an assumption of Gaussian noise, so it should be a decent fit for your problem. – spektr Jan 22 '16 at 20:44
• @choward there is also something else: not only my noise $N$ are samples from a normal distribution, but also the input observations $X$ are samples from a Gaussian distribution. is this in the assumption of Kalman filtering? – PickleRick Jan 22 '16 at 20:48
• Also posted on CS.SE. Please do not post the same question on multiple sites. Each community should have an honest shot at answering without anybody's time being wasted. – D.W. Jan 26 '16 at 7:59

So the way I went about formulating the problem was to essentially write the following equations:

The state that will be estimated, which is defined as a column vector, is the following: $$w = [vec(A)^{T},\bar{N}^{T}]^{T}$$ where $\bar{N}$ is the unknown average $N$ vector and $vec(A)$ is the vectorization operator on the unknown matrix A.

Based on the definition above, one can try to estimate $w$ by the following equations:

$$w_{i+1} = w_{i} + \eta_1$$ $$y_{i} = h(w_{i},X_{k}) + \eta_2$$

where $\eta_j$ are additive white noise vectors based on some variance values, $v_j$, and $X_{k}$ represents some $k^{th}$ data point from your set of values for $X$. You can select the $X_k$ value however you want (I did it randomly). One last equation you need is the following:

$$h(w,X) = [\hat{A},\hat{N}][X^{T},1]^{T}$$

where you extract $\hat{A}$ & $\hat{N}$ from the input $w$. Based on these equations, you can use a Kalman filter to iteratively solve this problem.

I implemented this earlier using an Unscented Kalman filter and validated it was converging to $A$ and $\bar{N}$. I will note that in the event $N$ is constant throughout all the samples, the algorithm I put together converged to $A$ and $N$ exactly (to numerical precision) within a reasonable number of iterations, but it didn't occur that way when $N$ varied randomly for each $X_k$, like you have in your problem.

Here is a sample image showing the convergence of the algorithm (which can be shown by the covariance trace moving towards 0):

Edit: I thought I would also post a plot of the $L_{2}$ Error history per iteration. Note that this $L_{2}$ Norm computation is based on the estimated $w$ vector vs the truth one. Now, this first plot is when $N$ is constant through every $X$ vector.

This next picture is when $N$ has a different value for every $X$ based on some Gaussian random sampling. Note, truth $\bar{N}$ in this case is the mixtureCenter variable (trying to estimate the mean $N$).

As you can see, estimating $N$ when it is constant results in some better estimations than when you have a slightly different $N$ for each $X$ sample, but it still gets a decent estimate in the harder situation. To further improve the convergence rates of the solution, one could do a batch processing via "sensor" fusion techniques. Below I show similar plots as above, but using a batch of 5 random $X$ vectors:

The first of these plots is for when the unknown $N$ is constant and the latter is when you are trying to estimate a mean $N$ based on random $N$ values for each $X$. Within a very small number of iterations, you can see great estimation performance in both examples.

My code can be found at the following github repository.

• thank you for your time and efforts! You made a pretty useful and educational answer, good job. The case I was interested in, is where for each observation of X, we have a random sample of N, from a Gaussian distribution. So N is not constant. At the end, I need parameters of mean and covariance of estimated N. May I ask to include the code you used to produce the results and plots in the answer, if it's possible? – PickleRick Jan 23 '16 at 22:55
• You should be able to use the estimate for $A$ and $\bar{N}$ to estimate the covariance matrix. All you end up doing is finding the random difference between unknown $N_{i}$ and $\bar{N}$ via: $Y_{i}-A*X_{i}-\bar{N} = \delta N_{i}$. Then the covariance, $Q$, is estimated via $Q = \frac{1}{n}\sum_{i}^{n} \delta N_{i}\delta N_{i}^{T}$. Also, the code isn't brief so I don't know if I will post it here. – spektr Jan 24 '16 at 0:10
• but wait...isn't Kalman filter for time series? my observations are not really time dependent and in $X_i, i$ is not time. they are only different observations from a distribution. – PickleRick Jan 24 '16 at 12:56
• @Anoosh No, you can modify Kalman filters, as I have, to estimate parameters. It's growing in usage in machine learning, as well. If you notice the dynamics equation with $w$, it has steady state dynamics (plus some noise). This is a typical formulation for estimating parameters. – spektr Jan 24 '16 at 13:51
• thank you! matlab has an implementation of kalman filter. I will try to use it for the example problem in my question. – PickleRick Jan 24 '16 at 13:52

Just to summarize the answers given for this question, here I provide a visualization of different methods carried out to learn the projection matrix $A$ given clean samples $X_i$ and noisy samples $Y_i$:

In the first row we can see 5 different projection matrices ($A$) chosen to project $X$ before adding Gaussian noise $N$.

In the second row, you can find the random projections initialized to be learned via different methods.

Third and fourth rows are related to a solution suggested in here.

In the third row, using the normal equation, we directly derive the projection matrix: $$\hat A = (X^T \times X)^{-1} \times X^T \times Y$$

In the fourth row, a linear regression using gradient descent is used to learn the projections via $X$ and $Y$.

And finally, the last row, using kalman filtering method suggested by choward.

As you can observe, the GD solution and the normal equation are very similar and since in this example, $(X^T \times X)^{-1}$ is feasible to calculate, it sounds a more reasonable solution and can be done in only one step. If the dimensionality of our features were too high, then we can use the GD since both $X$ and $Y$ are Gaussian.

The KF solution converges faster than GD, but can not achieve a good performance as GD, I believe.

After learning $A$, $N$ can be calculated via: $$\hat N=Y- \hat A \times X$$

and further $\hat {\mu_N}$ and $\hat {\Sigma_N}$ can be computed using $\hat N$.

• Cool visual! I will have to implement a gradient descent and see how it compares. I should say that if you just try to estimate just $A$ first, the KF does better (tried it this morning). I also believe you could tune the KF to be better than the one I implemented (I didn't really tune it). But nonetheless, cool visual and result! – spektr Jan 28 '16 at 16:02
• @choward thank you! you should try the regression solutions. I guess since my problem is much easier than a time-dependent problem with no Gaussian assumption, regression can solve it easily. Kalman should be more powerful in more difficult situation. e.g. a time serrie from an unknown distribution. Maybe you are right, maybe kalman also can be optimized. anyways, enjoyed a lot working on this question and learned about kalman filters power from you! – PickleRick Jan 28 '16 at 17:02
• This answer was copy-pasted on on CS.SE too. – D.W. Jan 28 '16 at 18:55
• @Anoosh Hey this is random at this point, but I have been doing some more advanced filtering stuff lately and improved working with Kalman Filters. I found the filter I gave you to use was not very well tuned and it can be improved. After retuning it, it performs VERY good. – spektr May 14 '16 at 2:04
• @choward very good! would be glad to give it a try tho! – PickleRick May 14 '16 at 13:11