# The derivative of the Cholesky Factor

Let $A$ be a symmetric, positive definite $p\times p$ matrix, and let $f(A)$ be it's Cholesky factor. That is, $f(A)$ is a lower triangular $p\times p$ matrix such that $A = f(A) f(A)^{\top}$. I am wondering if the derivative $$\frac{\mathrm{d}\operatorname{vech}\left(f(A)\right)}{\mathrm{d}\operatorname{vech}\left(A\right)}$$ is known, where $\operatorname{vech}$ is the half-vectorization function.

(I am conjecturing that it is something like $L^{\top} \left(f(A)\otimes A^{-1}\right) L$, where $L$ is the elimination matrix, but I would need a reference or proof anyway.)

The derivative can be found via implicit differentiation. That is, $$\frac{\mathrm{d}\operatorname{vec}\left(Y\right)}{\mathrm{d}\operatorname{vec}\left(X\right)} = \left(\frac{\mathrm{d} \operatorname{vec}\left(X\right)}{\mathrm{d}\operatorname{vec}\left(Y\right)}\right)^{-1}.$$ It is relatively easy to compute the derivative of $A$ with respect to $f(A)$ since $A = f(A)f(A)^{\top}$. The only trick part is restricting $f(A)$ to be lower triangular.

For general $X$, we have $$\frac{\mathrm{d} \operatorname{vec}\left(XX^{\top}\right)}{\mathrm{d} \operatorname{vec}\left(X\right)} = \left(I + K\right)\left(X\otimes I\right),$$ where $K$ is the Commutation Matrix.

Now to get the derivative with respect to the $\operatorname{vech}$ requires use of the chain rule. This gives $$\frac{\mathrm{d} \operatorname{vech}\left(XX^{\top}\right)}{\mathrm{d} \operatorname{vech}_{\Delta}\left(X\right)} = L \left(I + K\right)\left(X\otimes I\right) D,$$ where here $L$ is the elimination matrix, and $D$ is the "lower triangular duplication matrix" which has the property that $D \operatorname{vech}\left(M\right) = \operatorname{vec}\left(M\right)$ for lower triangular matrices $M$. The sought derivative is the matrix inverse of the above expression.

## numerical confirmation:

Here is a numerical confirmation in R: (note that the chol function in R is an operator from upper triangular matrices to upper triangular matrices, thus some mucking about with transposes):

require(matrixcalc)
set.seed(2349024)
n <- 6
X <- cov(matrix(rnorm(1000*n),ncol=n))
fnc <- function(X) t(chol(X))

Y <- fnc(X)
d0 <- (diag(1,nrow=n^2) + commutation.matrix(r=n)) %*% (Y %x% diag(1,nrow=n))
L <- elimination.matrix(n)
d1 <- L %*% d0 %*% t(L)
dfin <- solve(d1)

# now compute the approximate derivative
apx.d <- matrix(rep(NA,length(dfin)),nrow=dim(dfin)[1])
my.eps <- 1e-6
low.idx <- which(lower.tri(diag(1,n),diag=TRUE))
for (iii in c(1:length(low.idx))) {
Xalt <- X
tweak <- low.idx[iii]
Xalt[tweak] <- Xalt[tweak] + my.eps
# "Note that only the upper triangular part of 'x' is used..."
Yalt <- fnc(t(Xalt))
dY <- (Yalt - Y) / my.eps
apx.d[,iii] <- dY[low.idx]
}
apx.error <- apx.d - dfin
max(abs(apx.error))
apx.error


The maximum absolute error I get is 5.606e-07, on the order of the delta in the input variable, 1e-06.

• I am lead to believe there is a form of this in a paper by Fujikoshi and Okamoto in J. Japan. Stat. Soc., but I can find no more than the TOC for this journal: jss.gr.jp/ja/journal/jjss1976.html . – Steven Pav Jun 24 '14 at 19:13
• I think the paper is here, but it appears to talk about the eigendecomposition, not the Cholesky. – pete Aug 30 '14 at 3:43
• good find, @pete ! I believe you are right, but wonder if the F&O derivative could be somehow used for Cholesky decomposition--perhaps by considering the matrix 'square root' function. – Steven Pav Aug 30 '14 at 5:44
• Theorem A.1 of Simo Särkkä's Bayesian Filtering and Smoothing gives a result for the scalar derivative $\frac{\partial A}{\partial\theta} = A\Phi\left(A^{-1}\frac{\partial P}{\partial\theta}{}A^{-\top}\right)$, where $P=AA^{\top}$, and $\Phi_{ij}(M) = M_{ij}$ where $i > j$; $\tfrac12M_{ij}$ where $i=j$; and $0$ where $i < j$. Using this I get the same thing as using your result above. – pete Aug 31 '14 at 8:25
• I can generalise that result to give me $\frac{d\operatorname{vech}A}{d\operatorname{vech}P} = L(I\otimes A)Q(A^{-1}\otimes A^{-1})D$, where $Q$ is the diagonal matrix that gives us $Q\operatorname{vec}M = \operatorname{vec}(\Phi(M))$, and $D$ is the usual duplication matrix for symmetric matrices. This is more direct, and inverts $A$ rather than the larger matrix, but needs us to introduce the function $\Phi$ and matrix $Q$. – pete Sep 1 '14 at 2:42

I've written a relevant note on arXiv: http://arxiv.org/abs/1602.07527

I included the neat closed form solution pete gives in a comment, and also a messy expression (converted to the notation f = chol(A)): $$\frac{\partial f_{ij}}{\partial A_{kl}} = \bigg(\sum_{m>j} f_{im}f_{mk}^{-1} + \tfrac{1}{2}f_{ij}f_{jk}^{-1}\bigg)f_{jl}^{-1} + (1-\delta_{kl})\bigg(\sum_{m>j} f_{im}f_{ml}^{-1} + \tfrac{1}{2}f_{ij}f_{jl}^{-1}\bigg)f_{jk}^{-1}.$$

However, if you're interested in differentiating a larger expression, you can do that in $O(N^3)$, without computing all $O(N^4)$ derivatives in $\frac{\partial \mathrm{vech}(f)}{\partial \mathrm{vech}A}$. The note explains different ways to do that.

(pete: if you tell me who you are, I'll add a proper acknowledgement to my note in any future revision.)

The following may be of help. If $A=T'T$ (where $T$ is upper triangular), then you can show that (see. Thm. 2.1.9 in Aspects of Multivariate Statistical Theory by R. J. Muirhead): \begin{equation*} (dA) = 2^p\prod_{i=1}^p t_{ii}^{p+1-i}(dT), \end{equation*} where using exterior products, we define $(dA) := \bigwedge_{i \le j}^p da_{ij}$; similarly, $(dT)=\bigwedge_{i\le j}^p dt_{ij}$.

Note: You are essentially asking for $dA$ not $(dA)$ --- if all elements of $A$ are distinct, then of course both coincide, otherwise not.

Let $T=\{M|M \text{ lower triangular }\},T^+=\{M\in T|\;m_{i,i}>0 \},S^+=\{M|M \text{ SPD }\}$ and $\phi:C=[c_{i,j}]\in T^+\rightarrow h(CC^T)\in h(S^+)\subset T$ where the "half" function $h$ sends to zero the strict upper part and keeps invariant the lower part of a symmetric matrix. $D\phi_C:K\in T\rightarrow H=h(KC^T+CK^T)\in T$. Then $\det(D\phi_C)=2^p\Pi_{j=1}^pc_{j,j}^{2p-2j+1}>0$ and $D\phi_C$ is an isomorphism. Here we are interested with $(D\phi_C)^{-1}:H\in T\rightarrow K\in T$. The linear equation, in the unknown $K$, $H=h(KC^T+CK^T)$ has a unique solution which is $D(vechf)_{vechA}(H)$ and we are done.