# Discrepancy in Cholesky decomposition matrix from variance covariance matrix obtained in Stata and R using regression

This is a question to someone who knows both R and Stata. For a simple linear regression model, I was trying to estimate the Cholesky decomposition matrix from the variance-covariance matrix of the parameters obtained. The obtained Cholesky decomposition matrix varies between R and Stata. Does anyone know the reason?

## The code for R:

example <-foreign::read.dta("https://dss.princeton.edu/training/states.dta")
lmfit <- lm(csat ~ expense, data = example)
vcov_matrix <- vcov(lmfit)
chol_decomp_matrix <- chol(vcov_matrix)


chol_decomp_matrix provides

            (Intercept)      expense
(Intercept)     32.7009 -0.005835735
expense          0.0000  0.001546270


## The code for Stata:

use https://dss.princeton.edu/training/states.dta
regress csat expense
matrix list e(V)
matrix vcov_mat  = e(V)
matrix chol_mat =  cholesky(vcov_mat)
matrix list chol_mat


chol_mat provides

            expense      _cons
expense     0.006037    0
-cons      -31.610095   8.375592


Can anyone please explain the difference?

From the Stata documentation, the cholesky function returns a lower triangular matrix $$G$$ such that $$A = GG'$$ given input matrix $$A$$. Meanwhile, from the R documentation, chol returns an upper triangular matrix $$M$$ such that $$X = M'M$$ for input matrix $$X$$. From this we would expect them to be transposes of each other. However, the outputs differ even more significantly because R and Stata have flipped the variable ordering -- the R function lists intercept first and expense second, while Stata lists expense first and intercept second.

The second difference (flipping the variable ordering) is what yields the very different-looking Cholesky decompositions. To see why, note that the following are the Cholesky decompositions of a matrix and the same matrix with its rows and columns flipped:

$$\left[\begin{matrix} a & b \\ b & c\\ \end{matrix}\right] = \left[\begin{matrix} \sqrt{a} & 0 \\ b/\sqrt{a} & \sqrt{c-b^2/a}\end{matrix}\right]\left[\begin{matrix} \sqrt{a} & b/\sqrt{a} \\ 0 & \sqrt{c-b^2/a}\end{matrix}\right] \\ \left[\begin{matrix} c & b \\ b & a\\ \end{matrix}\right] = \left[\begin{matrix} \sqrt{c} & 0 \\ b/\sqrt{c} & \sqrt{a-b^2/c}\end{matrix}\right]\left[\begin{matrix} \sqrt{c} & b/\sqrt{c} \\ 0 & \sqrt{a-b^2/c}\end{matrix}\right]$$

Note that none of the corresponding non-zero entries in the Cholesky decomposition of a 2x2 symmetric positive definite matrix will match after flipping the rows and columns unless the two elements on the main diagonal ($$a$$ and $$c$$ above) are equal.

Though the Cholesky decompositions are different, they both encode the exact same matrix (up to row and column reordering), and there is no real advantage to using one versus the other (for instance, all entries are real-valued whenever the input matrix is positive semidefinite in either case). You can check that indeed the two decompositions in your example encode an identical variance covariance matrix (up to numerical precision of the output in your question):

chol.r <- rbind(c(32.7009, -0.005835735), c(0.0000, 0.001546270))
t(chol.r) %*% chol.r
#              [,1]          [,2]
# [1,] 1069.3488608 -1.908338e-01
# [2,]   -0.1908338  3.644675e-05
chol.stata <- rbind(c(0.006037, 0), c(-31.610095, 8.375592))
(chol.stata %*% t(chol.stata))[2:1,2:1]
#              [,1]          [,2]
# [1,] 1069.3486473 -1.908301e-01
# [2,]   -0.1908301  3.644537e-05

• Thanks for the response @josliber. But does this order matter? if yes, how is it explained? any idea? Mar 2, 2020 at 16:03
• @SheejaKrishnan the order does not matter, but it will lead to a different Cholesky decomposition. Give me a few minutes and I can edit in an example. Mar 2, 2020 at 16:03
• @ josliber, For regression with many covariates and other types of regression, this is difficult to spot. I am using the Cholesky decomposition matrix in my model to generate the regression parameters. This will have then a significant effect that depends on R or STATA. Can you suggest which is the correct method, then? Mar 2, 2020 at 16:08
• @SheejaKrishnan this should have no effect -- the Cholesky decompositions are encoding the same matrix, up to row and column ordering. The only thing you would have to do is to make sure your code that uses the Cholesky decomposition is using the same variable ordering that is being returned by your software (R or Stata). Mar 2, 2020 at 16:22
• +1 This is an impressive answer because it quickly establishes its scope and the relevant definitions, is clearly argued, has a nice (and explicit) theoretical answer, and is supplemented by a computed example.
– whuber
Mar 2, 2020 at 20:22