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?
 A: 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

