Let $A$ be a symmetric, positive definite $p\times p$ matrix, and let $f(A)$ be its 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.)


4 Answers 4


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):

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

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

  • $\begingroup$ 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 . $\endgroup$
    – Steven Pav
    Jun 24, 2014 at 19:13
  • 1
    $\begingroup$ I think the paper is here, but it appears to talk about the eigendecomposition, not the Cholesky. $\endgroup$
    – pete
    Aug 30, 2014 at 3:43
  • $\begingroup$ 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. $\endgroup$
    – Steven Pav
    Aug 30, 2014 at 5:44
  • 1
    $\begingroup$ 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. $\endgroup$
    – pete
    Aug 31, 2014 at 8:25
  • 2
    $\begingroup$ 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$. $\endgroup$
    – pete
    Sep 1, 2014 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.


Your Answer

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

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