I have implemented the forward-filtering-backwards-sampling (ffbs) algorithm. It consists of kalman filtering forward in time (to obtain mean and sigma). Then it uses these values and the Kalman smoother backwards in time to get mean and Sigma for all time points.

My problem is that the Kalman smoother seems to get negative variances. Can this be the case - and if so, how does it affect my results?

The code is made in R (and Sigma_predictions and Sigma_upates are obtained from the Kalman filter)


# Kalman filter iterations

# Kalman loop starts
for(i in 1:big_t){ # Loop through each time point

  # In the first iteration the inittial value is set
  thisIsTheFirstValue <- (i == 1)
    phi_upd <- mean_0
    sigma_upd <- sigma_0

  # Prediction step
  phi_pre <- (Gamma %*% phi_upd) + Theta
  sigma_pre <- (Gamma %*% sigma_upd %*% t(Gamma)) + Psi_cov

  # Update step
  v <- y[, i] - (alpha + B %*% phi_pre) # y_t - f_t
  big_F <- (B %*% sigma_pre %*% t(B)) + epsilon_cov # epsilon_cov = sigma * Idnetity matrix
  inv_F <- solve(big_F, tol = 1e-40)
  # Calculate updates
  phi_upd <- phi_pre + (sigma_pre %*% t(B) %*% inv_F %*% v)
  sigma_upd <- sigma_pre - (sigma_pre %*% t(B) %*% inv_F %*%  B %*% sigma_pre)

  # save values in the arrays (used later in backward sampling step)
  phis[, i] <- phi_upd
  phi_predictions[ , , i] <- phi_pre
  sigma_predictions[ , , i] <- sigma_pre
  phi_updates[ , , i] <- phi_upd
  sigma_updates[ , , i] <- sigma_upd

} # end for loop


# Find h_t and H_t by backwards sampling
loop_seq <- seq(from = (big_t - 1), to = 1, by = -1)
for(t in loop_seq){

  # In the first iteration the inittial value is set
  thisIsTheSecondLastTimePoint <- (t == (big_t-1))
    phi_next <- matrix(rmvnorm(1, mean = phi_updates[ , , big_t], sigma = (sigma_updates[ , ,  big_t])), nrow = (p+1), ncol = 1)
# Save phi drawing
  phi_array[ , , big_t] <- phi_next

  # Calculate mean value dep. on this value
  inv_sig_pred <- solve(sigma_predictions[ , , t+1]) 
  mean_t <- phi_updates[, , t] + (sigma_updates[, , t] %*% t(Gamma) %*% inv_sig_pred %*% (phi_array[ , , (t+1)] - phi_predictions[ , , t+1]))
  sigma_t <- sigma_updates[ , , t] - (sigma_updates[, , t] %*% t(Gamma) %*% inv_sig_pred %*% Gamma %*% sigma_updates[, , t])
  # Make the drawing
  phi_next <- t(rmvnorm(1, mean = mean_t, sigma = sigma_t))

  # Save phi drawings
  phi_array[, , t] <- phi_next

It is a state-space system, I am considering. lambda is set to 0.5. All is set due to an article, I am following. enter image description here

The mathematical calculations are (sorry for a changed notation): enter image description here enter image description here

The problems, I am facing are problems with diagonal elements being negative in the H_t-matrix. The diagonal of the first calculated H_t I get: enter image description here The diagonal of the second calculated H_t I get: enter image description here

  • $\begingroup$ Is t(Gamma) > 1? I don't know what your algorithm is supposed to do, so I don't know whether your problem is due to algorithm deficiency or implementation error. Is sigma_t actually a covariance matrix (not just 1 by 1 variance)? $\endgroup$ May 26, 2018 at 1:42
  • $\begingroup$ What is t(Gamma)? Where is the closing parenthesis for the first character in: (sigma_updates[, , t] %*% t(Gamma) %*% $\endgroup$ May 26, 2018 at 1:53
  • $\begingroup$ I think it is there (or am I wrong?). The line was just broken into two lines. It should be seen clearer now, as I set them to one. $\endgroup$
    – Dorthe
    May 26, 2018 at 2:06
  • $\begingroup$ Show us mathematically what the covariance update formulas are, and specifically what covariance is coming out with negative (variances?, meaning diagonal elements? or eigenvalues?). Are there some number of updates before the negative variances occur? $\endgroup$ May 26, 2018 at 2:09
  • $\begingroup$ Hope it makes sense. It is already in the first iteration, that negative variances occur. $\endgroup$
    – Dorthe
    May 26, 2018 at 2:53

1 Answer 1


You appear to be getting roundoff error level negative variances. Subtraction is a dangerous thing. There is a perhaps a better "root cause" fixing up which can be done, by modifying the updating equation.

But in the absence of that, one quick "fix", which I'm not saying is a good thing to do here, is to adjust the covariance matrix on each iteration to be symmetric positive semidefinite, with some minimum eigenvalue value, eigmin (a very small positive number). Do this by first symmetrizing, then eigenvalue eigenvector decomposiition, then replacing eigenvalues below eigmin with eigmin, then reconstituting the matrix. Symmetrize if needed at the end.

Symmetrization: $H = 0.5*(H + H^T)$

Let V = eigenvector matrix of H, D = diagonal; matrix of eigenvalues of H. Therefore, $H = VDV^{-1}$. Replace diagonal elements of D below eigmin by eigmin, resulting in $D_{adjusrted}$. Now form $VD_{adjusted}V^{-1}$, and symmetrize if necessary.

  • $\begingroup$ @Dorthe Does my answer help you or is it off the mark? $\endgroup$ May 30, 2018 at 0:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.