Matrix calculus

Introduction

Estimators are usually written down as the solution to an optimization problem in many variables at once. The least squares coefficient minimizes a sum of squared residuals over a \(k\)-vector, a maximum likelihood estimator maximizes a log-likelihood over a parameter vector, and a method of moments estimator minimizes a quadratic form in a vector of sample moments. In each case the object being differentiated is a scalar and the object being differentiated with respect to is a vector or a matrix, so the derivative is a collection of partial derivatives arranged in the shape of the argument. Matrix calculus is the set of conventions that makes that arrangement systematic, together with a handful of results that are used often enough to be worth memorizing.

This post derives the results that cover most of what applied work requires. Four of them concern a scalar differentiated with respect to a vector or a matrix: the linear form \(a^{T}x\), the quadratic form \(x^{T}Ax\), the least squares objective \(\lVert Ax - b \rVert_{2}^{2}\), and \(x^{T}Ax\) differentiated with respect to \(A\) rather than \(x\). Two further sections add what estimation beyond linear least squares needs: the Jacobian of a vector-valued function together with the chain rule, which delivers the first-order conditions of nonlinear least squares and the generalized method of moments as well as the delta method, and the derivatives of the log determinant and the matrix inverse, which are what maximum likelihood for a covariance matrix requires. Every result is obtained the same way, by writing the scalar in summation notation, differentiating with respect to a single component, and recognizing the resulting expression as a matrix or vector product. The matrix identities are bookkeeping on top of that procedure. A final section repeats the derivations with differentials, which is faster once the notation is familiar and serves here as a check.

Conventions

Two conventions have to be fixed before anything can be derived, because the literature contains both and they differ by a transpose.

Let \(f : \mathbb{R}^{n} \to \mathbb{R}\) be a scalar function of a column vector \(x\) of size \(n \times 1\). The derivative \(\partial f / \partial x\) is defined here as the \(n \times 1\) column vector whose \(i\)-th entry is \(\partial f / \partial x_{i}\):

\[ \frac{\partial f}{\partial x} = \begin{bmatrix} \partial f / \partial x_{1} \\ \partial f / \partial x_{2} \\ \vdots \\ \partial f / \partial x_{n} \end{bmatrix}. \]

This is the denominator layout, also called the gradient convention, and it has the property that the derivative has the same shape as the argument. The alternative, numerator layout, stores the same numbers in a \(1 \times n\) row. Nothing substantive depends on the choice, but signs and transposes in intermediate steps do, so it is stated once and used throughout.

The same rule extends to a scalar function of a matrix. If \(f : \mathbb{R}^{m \times n} \to \mathbb{R}\) is a scalar function of a matrix \(A\), then \(\partial f / \partial A\) is the \(m \times n\) matrix whose entry in row \(k\) and column \(l\) is \(\partial f / \partial A_{kl}\). Again the derivative has the shape of the argument.

One piece of notation recurs below. The Kronecker delta \(\delta_{ik}\) equals \(1\) when \(i = k\) and \(0\) otherwise, so that \(\partial x_{i} / \partial x_{k} = \delta_{ik}\). It is the componentwise statement that the entries of \(x\) are separate variables, and it is what allows a double sum to collapse to a single sum after differentiation.

The linear form

Let \(x\) and \(a\) both be \(n \times 1\) column vectors, with \(a\) constant, and consider

\[ f(x) = a^{T} x . \]

The product \(a^{T}x\) is \(1 \times n\) times \(n \times 1\), so it is a scalar, and in summation notation it is

\[ f(x) = \sum_{i=1}^{n} a_{i} x_{i} . \]

Differentiate with respect to one component \(x_{k}\). Only the term with \(i = k\) involves \(x_{k}\), and the coefficient \(a_{k}\) is constant, so

\[ \frac{\partial f}{\partial x_{k}} = \sum_{i=1}^{n} a_{i} \frac{\partial x_{i}}{\partial x_{k}} = \sum_{i=1}^{n} a_{i} \delta_{ik} = a_{k} . \]

The \(k\)-th partial derivative is the \(k\)-th entry of \(a\), for every \(k\), so stacking the partials into a column gives

\[ \frac{\partial}{\partial x} \left( a^{T} x \right) = a . \]

Since \(a^{T}x\) is a scalar it equals its own transpose, \(a^{T}x = x^{T}a\), and the same derivative applies to that form as well. The result is the matrix analogue of \(\mathrm{d}(ax)/\mathrm{d}x = a\) for scalars, which is the reason the denominator layout is convenient: the shape of the answer matches the shape of \(x\) with no transpose inserted.

The quadratic form

Let \(x\) be \(n \times 1\) and let \(A\) be a constant \(n \times n\) matrix that is not assumed to be symmetric. Consider

\[ f(x) = x^{T} A x . \]

Checking dimensions, \(x^{T}\) is \(1 \times n\), \(A\) is \(n \times n\) and \(x\) is \(n \times 1\), so the product is \(1 \times 1\) and \(f\) is a scalar.

The case \(n = 2\)

Before doing the general case it helps to write out a small one, because the pattern is visible there. With

\[ x = \begin{bmatrix} x_{1} \\ x_{2} \end{bmatrix}, \qquad A = \begin{bmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{bmatrix}, \]

the inner product \(Ax\) is

\[ Ax = \begin{bmatrix} A_{11}x_{1} + A_{12}x_{2} \\ A_{21}x_{1} + A_{22}x_{2} \end{bmatrix}, \]

and multiplying on the left by \(x^{T}\) gives

\[ f(x) = x_{1}\left( A_{11}x_{1} + A_{12}x_{2} \right) + x_{2}\left( A_{21}x_{1} + A_{22}x_{2} \right) = A_{11}x_{1}^{2} + A_{12}x_{1}x_{2} + A_{21}x_{1}x_{2} + A_{22}x_{2}^{2} . \]

This is an ordinary quadratic polynomial in two variables, so it can be differentiated by the usual rules:

\[ \frac{\partial f}{\partial x_{1}} = 2A_{11}x_{1} + \left( A_{12} + A_{21} \right) x_{2}, \qquad \frac{\partial f}{\partial x_{2}} = \left( A_{12} + A_{21} \right) x_{1} + 2A_{22}x_{2} . \]

Stacking the two partials into a column and factoring out \(x\),

\[ \frac{\partial f}{\partial x} = \begin{bmatrix} 2A_{11} & A_{12} + A_{21} \\ A_{12} + A_{21} & 2A_{22} \end{bmatrix} \begin{bmatrix} x_{1} \\ x_{2} \end{bmatrix} . \]

The matrix in front is \(A + A^{T}\), entry by entry, so the conjecture is \(\partial f / \partial x = (A + A^{T})x\). Two features of the two-variable case are worth carrying forward. The off-diagonal coefficients appear only through the sum \(A_{12} + A_{21}\), so the cross term cannot distinguish \(A\) from any other matrix with the same off-diagonal sum, and the diagonal terms pick up a factor of two, which is the derivative of \(x_{k}^{2}\).

The general case

In summation notation the quadratic form is a double sum,

\[ f(x) = x^{T} A x = \sum_{i=1}^{n} \sum_{j=1}^{n} A_{ij} x_{i} x_{j} , \]

which can be verified by expanding \(Ax\) into \(\sum_{j} A_{ij}x_{j}\) for row \(i\) and then taking the inner product with \(x\). Differentiating with respect to \(x_{k}\) requires the product rule, because \(x_{k}\) can appear in either factor:

\[ \frac{\partial f}{\partial x_{k}} = \sum_{i=1}^{n} \sum_{j=1}^{n} A_{ij} \frac{\partial (x_{i}x_{j})}{\partial x_{k}} = \sum_{i=1}^{n} \sum_{j=1}^{n} A_{ij} \left( \delta_{ik} x_{j} + x_{i} \delta_{jk} \right) . \]

Split the expression into two double sums. In the first, \(\delta_{ik}\) kills every term except \(i = k\), leaving a single sum over \(j\); in the second, \(\delta_{jk}\) kills every term except \(j = k\), leaving a single sum over \(i\):

\[ \frac{\partial f}{\partial x_{k}} = \sum_{j=1}^{n} A_{kj} x_{j} + \sum_{i=1}^{n} A_{ik} x_{i} . \]

Both sums are entries of matrix products. The first is row \(k\) of \(A\) times \(x\), that is \((Ax)_{k}\). The second runs over the row index of \(A\) with the column index fixed at \(k\), so it is column \(k\) of \(A\) times \(x\), which is row \(k\) of \(A^{T}\) times \(x\), that is \((A^{T}x)_{k}\). Hence

\[ \frac{\partial f}{\partial x_{k}} = (Ax)_{k} + (A^{T}x)_{k} = \left( \left( A + A^{T} \right) x \right)_{k} , \]

and since this holds for every \(k\), stacking gives

\[ \frac{\partial}{\partial x} \left( x^{T} A x \right) = \left( A + A^{T} \right) x . \]

When \(A\) is symmetric, \(A^{T} = A\) and the result simplifies to \(2Ax\), which is the analogue of \(\mathrm{d}(ax^{2})/\mathrm{d}x = 2ax\). The symmetric case is the one that occurs in practice, since a quadratic form only ever sees the symmetric part of its matrix: because \(x^{T}Ax\) is a scalar it equals its transpose \(x^{T}A^{T}x\), so

\[ x^{T} A x = \tfrac{1}{2} \left( x^{T} A x + x^{T} A^{T} x \right) = x^{T} \left[ \tfrac{1}{2} \left( A + A^{T} \right) \right] x . \]

Replacing \(A\) by its symmetric part \(\tfrac{1}{2}(A + A^{T})\) changes neither the value of the quadratic form nor its gradient, which is why the antisymmetric part of \(A\) is unidentified from the function \(f\) alone. In the \(n = 2\) case above, this is the observation that only \(A_{12} + A_{21}\) entered.

The least squares objective

Let \(A\) be \(m \times n\), \(x\) be \(n \times 1\) and \(b\) be \(m \times 1\), and consider the squared Euclidean distance between \(Ax\) and \(b\),

\[ f(x) = \lVert Ax - b \rVert_{2}^{2} . \]

The residual \(Ax - b\) is \(m \times 1\), and the squared norm of a vector \(v\) is \(v^{T}v\), so \(f\) is a scalar. This is the sum of squared residuals of a linear regression with design matrix \(A\), outcome \(b\) and coefficient vector \(x\), written with the symbols of linear algebra rather than those of econometrics.

Componentwise

Write the \(i\)-th residual as \(r_{i} = \sum_{j=1}^{n} A_{ij} x_{j} - b_{i}\), so that

\[ f(x) = \sum_{i=1}^{m} r_{i}^{2} = \sum_{i=1}^{m} \left( \sum_{j=1}^{n} A_{ij} x_{j} - b_{i} \right)^{2} . \]

The chain rule applied to the square, together with \(\partial r_{i} / \partial x_{k} = A_{ik}\), gives

\[ \frac{\partial f}{\partial x_{k}} = \sum_{i=1}^{m} 2 r_{i} \frac{\partial r_{i}}{\partial x_{k}} = 2 \sum_{i=1}^{m} A_{ik} r_{i} . \]

The sum runs over the row index of \(A\) with the column index fixed at \(k\), which as before is row \(k\) of \(A^{T}\) applied to \(r\). Therefore

\[ \frac{\partial f}{\partial x} = 2 A^{T} r = 2 A^{T} \left( Ax - b \right) . \]

By expansion

The same result follows from the two previous sections without any indices. Expanding the squared norm,

\[ f(x) = (Ax - b)^{T} (Ax - b) = x^{T} A^{T} A x - x^{T} A^{T} b - b^{T} A x + b^{T} b . \]

The two middle terms are scalars and are transposes of each other, since \((x^{T}A^{T}b)^{T} = b^{T}Ax\), so they are equal and can be combined:

\[ f(x) = x^{T} \left( A^{T} A \right) x - 2 \left( A^{T} b \right)^{T} x + b^{T} b . \]

The three terms are now a quadratic form, a linear form and a constant. The matrix \(A^{T}A\) is symmetric, because \((A^{T}A)^{T} = A^{T}A\), so the quadratic form contributes \(2A^{T}Ax\). The linear form is \(a^{T}x\) with \(a = A^{T}b\) and contributes \(-2A^{T}b\). The constant contributes nothing. Adding them,

\[ \frac{\partial f}{\partial x} = 2 A^{T} A x - 2 A^{T} b = 2 A^{T} \left( Ax - b \right) , \]

which agrees with the componentwise derivation. The dimensions check: \(A^{T}\) is \(n \times m\) and the residual is \(m \times 1\), so the gradient is \(n \times 1\), the shape of \(x\).

What the first-order condition says

Setting the gradient to zero gives the normal equations,

\[ A^{T} A x = A^{T} b , \]

and if \(A\) has full column rank then \(A^{T}A\) is invertible and the stationary point is unique:

\[ \hat{x} = \left( A^{T} A \right)^{-1} A^{T} b . \]

That the stationary point is a minimum rather than a maximum or a saddle follows from the second derivative. Differentiating the gradient once more gives the Hessian

\[ \frac{\partial^{2} f}{\partial x \, \partial x^{T}} = 2 A^{T} A , \]

and for any \(v \neq 0\), \(v^{T}(A^{T}A)v = \lVert Av \rVert_{2}^{2} \geq 0\), so the Hessian is positive semidefinite everywhere and \(f\) is convex. Under full column rank, \(Av = 0\) only for \(v = 0\), so the Hessian is positive definite and the minimum is strict.

The first-order condition also has a geometric reading. Written as \(A^{T}(Ax - b) = 0\), it says that the residual at the optimum is orthogonal to every column of \(A\), since row \(k\) of \(A^{T}\) is column \(k\) of \(A\). The fitted vector \(A\hat{x}\) is therefore the projection of \(b\) onto the column space of \(A\), and the orthogonality that is usually taken as the defining property of a projection appears here as the first-order condition of a minimization problem rather than as an assumption.

Vector-valued functions

Every function differentiated so far has been scalar-valued, and the estimators that occur in practice are rarely presented that way. Nonlinear least squares squares a vector of residuals, the generalized method of moments weights a vector of sample moments, and the delta method transforms an estimator through a vector of smooth functions. Each of these is a scalar criterion composed with a vector-valued function of the parameters, so what is needed is the derivative of a vector with respect to a vector together with a rule for composing it with the scalar results above.

Let \(g : \mathbb{R}^{n} \to \mathbb{R}^{m}\) have components \(g_{1}, \dots, g_{m}\), each a scalar function of the \(n \times 1\) vector \(x\). The derivative of \(g\) with respect to \(x\) is the \(m \times n\) Jacobian matrix

\[ J(x) = \frac{\partial g}{\partial x^{T}} = \begin{bmatrix} \partial g_{1} / \partial x_{1} & \cdots & \partial g_{1} / \partial x_{n} \\ \vdots & & \vdots \\ \partial g_{m} / \partial x_{1} & \cdots & \partial g_{m} / \partial x_{n} \end{bmatrix}, \]

whose entry in row \(i\) and column \(j\) is \(\partial g_{i} / \partial x_{j}\). Row \(i\) is therefore the transpose of the gradient of the \(i\)-th component. The notation \(\partial g / \partial x^{T}\) records the shape: a column of outputs differentiated by a row of inputs produces a matrix with one row per output and one column per input.

The chain rule

Let \(f : \mathbb{R}^{m} \to \mathbb{R}\) be a scalar function and consider the composition \(h(x) = f(g(x))\). The ordinary chain rule for functions of several variables says that a change in \(x_{k}\) reaches \(h\) through each of the \(m\) intermediate variables in turn,

\[ \frac{\partial h}{\partial x_{k}} = \sum_{i=1}^{m} \left. \frac{\partial f}{\partial u_{i}} \right|_{u = g(x)} \frac{\partial g_{i}}{\partial x_{k}} = \sum_{i=1}^{m} J_{ik} \left( \frac{\partial f}{\partial u} \right)_{i} , \]

and the sum over the row index of \(J\) with the column index fixed at \(k\) is entry \(k\) of \(J^{T}\) applied to the gradient of \(f\). Stacking over \(k\),

\[ \frac{\partial h}{\partial x} = J(x)^{T} \left. \frac{\partial f}{\partial u} \right|_{u = g(x)} . \]

The transpose is forced by the shapes: \(J^{T}\) is \(n \times m\) and the gradient of \(f\) is \(m \times 1\), so the product is \(n \times 1\), which is the shape of \(x\).

The least squares gradient is the first check. With \(g(x) = Ax - b\), each component is \(g_{i}(x) = \sum_{j} A_{ij}x_{j} - b_{i}\), so \(\partial g_{i} / \partial x_{j} = A_{ij}\) and the Jacobian is \(A\) itself. With \(f(u) = u^{T}u\), the quadratic form result with identity matrix gives \(\partial f / \partial u = 2u\). The chain rule then returns \(2A^{T}(Ax - b)\), which is what the two earlier derivations produced.

Weighted criteria

The criterion that appears in minimum distance estimation, in the generalized method of moments and in generalized least squares is a quadratic form in a vector-valued function,

\[ Q(\theta) = g(\theta)^{T} W g(\theta) , \]

where \(g(\theta)\) is \(m \times 1\), the weight matrix \(W\) is \(m \times m\), symmetric, positive definite and does not depend on \(\theta\), and \(\theta\) is \(n \times 1\) with \(m \geq n\). Since \(W\) is symmetric, the quadratic form result gives \(\partial (u^{T}Wu) / \partial u = 2Wu\), and the chain rule gives

\[ \frac{\partial Q}{\partial \theta} = 2 J(\theta)^{T} W g(\theta) , \qquad J(\theta) = \frac{\partial g}{\partial \theta^{T}} . \]

The first-order condition \(J(\theta)^{T} W g(\theta) = 0\) is a system of \(n\) equations in \(n\) unknowns regardless of how many moments there are. When \(m = n\) and \(J\) is invertible it reduces to \(g(\hat{\theta}) = 0\), so the moment conditions hold exactly in the sample, which is the just identified case. When \(m > n\) they cannot all be set to zero, and the condition instead requires the weighted moment vector \(Wg(\hat{\theta})\) to be orthogonal to the columns of \(J(\hat{\theta})\). Setting \(g(\theta) = A\theta - b\) and \(W = I\) recovers the normal equations of the previous section, so weighted nonlinear estimation and least squares share one first-order condition with different ingredients.

The delta method

The Jacobian also governs how sampling variance is transported through a transformation. Suppose \(\sqrt{N}(\hat{\theta} - \theta_{0})\) converges in distribution to \(N(0, V)\) and let \(h : \mathbb{R}^{n} \to \mathbb{R}^{m}\) be continuously differentiable at \(\theta_{0}\) with Jacobian \(H = \partial h / \partial \theta^{T}\) evaluated at \(\theta_{0}\). A first-order Taylor expansion around \(\theta_{0}\) gives \(h(\hat{\theta}) \approx h(\theta_{0}) + H(\hat{\theta} - \theta_{0})\), and a linear map applied to a normal limit produces a normal limit with a sandwiched variance,

\[ \sqrt{N} \left( h(\hat{\theta}) - h(\theta_{0}) \right) \xrightarrow{d} N \left( 0, \, H V H^{T} \right) . \]

For a scalar transformation, \(m = 1\), the Jacobian is the transpose of the gradient and the asymptotic variance is the quadratic form \((\partial h / \partial \theta)^{T} V (\partial h / \partial \theta)\). This is the standard error reported for a ratio of coefficients, an elasticity or a marginal effect. Taking \(h(\beta) = \beta_{1} / \beta_{2}\) as an example, the gradient is \((1/\beta_{2}, \, -\beta_{1}/\beta_{2}^{2})^{T}\), so the asymptotic variance is

\[ \frac{V_{11}}{\beta_{2}^{2}} - \frac{2 \beta_{1} V_{12}}{\beta_{2}^{3}} + \frac{\beta_{1}^{2} V_{22}}{\beta_{2}^{4}} , \]

which makes visible why such standard errors become unreliable as \(\beta_{2}\) approaches zero.

Differentiating with respect to a matrix

The final case reverses the roles. Let \(x\) be a constant \(n \times 1\) vector and let the \(n \times n\) matrix \(A\) be the variable, and consider

\[ f(A) = x^{T} A x . \]

The function is the same quadratic form as before, but it is now read as a function of \(n^{2}\) variables, the entries of \(A\), with \(x\) held fixed. By the convention stated above, \(\partial f / \partial A\) is the \(n \times n\) matrix whose \((k,l)\) entry is \(\partial f / \partial A_{kl}\).

Start from the summation form,

\[ f(A) = \sum_{i=1}^{n} \sum_{j=1}^{n} A_{ij} x_{i} x_{j} . \]

Viewed as a function of the entries of \(A\), this is linear: each entry \(A_{ij}\) appears exactly once, multiplied by the constant \(x_{i}x_{j}\). Differentiating with respect to a single entry \(A_{kl}\) therefore returns that constant,

\[ \frac{\partial f}{\partial A_{kl}} = \sum_{i=1}^{n} \sum_{j=1}^{n} \frac{\partial A_{ij}}{\partial A_{kl}} x_{i} x_{j} = \sum_{i=1}^{n} \sum_{j=1}^{n} \delta_{ik} \delta_{jl} x_{i} x_{j} = x_{k} x_{l} , \]

where the two deltas select \(i = k\) and \(j = l\). The entry in row \(k\) and column \(l\) of the answer is the product of the \(k\)-th and \(l\)-th entries of \(x\), and the matrix with that entry pattern is the outer product \(x x^{T}\):

\[ \frac{\partial}{\partial A} \left( x^{T} A x \right) = x x^{T} . \]

For \(n = 2\) the claim can be read off the expansion used earlier, \(f = A_{11}x_{1}^{2} + A_{12}x_{1}x_{2} + A_{21}x_{2}x_{1} + A_{22}x_{2}^{2}\), whose four partial derivatives with respect to \(A_{11}, A_{12}, A_{21}, A_{22}\) are \(x_{1}^{2}, x_{1}x_{2}, x_{2}x_{1}, x_{2}^{2}\), arranged as

\[ \begin{bmatrix} x_{1}^{2} & x_{1}x_{2} \\ x_{2}x_{1} & x_{2}^{2} \end{bmatrix} = \begin{bmatrix} x_{1} \\ x_{2} \end{bmatrix} \begin{bmatrix} x_{1} & x_{2} \end{bmatrix} = x x^{T} . \]

Three remarks follow. The derivative does not depend on \(A\), which is the matrix version of the fact that the derivative of a linear function is constant. The result is symmetric and of rank one, since it is an outer product of a vector with itself. And the placement of the transpose matters in the same way as in the projection formula: \(x^{T}x\) is a scalar while \(x x^{T}\) is an \(n \times n\) matrix, so the two expressions built from the same symbols are objects of different types.

A caveat applies when \(A\) is constrained to be symmetric. The derivation above treats the \(n^{2}\) entries as free variables, which is the correct reading of the question as posed. If instead \(A_{kl}\) and \(A_{lk}\) are required to move together, the off-diagonal partials collect both contributions and the derivative becomes \(2xx^{T} - \mathrm{diag}(x x^{T})\). The two answers disagree about which variables are being varied, not about the calculus, so it is worth stating which convention is in force whenever a derivative with respect to a symmetric matrix is reported.

The log determinant and the inverse

Likelihood functions for multivariate models contain two objects that the rules so far do not cover. The normalizing constant of a multivariate normal density contributes a determinant, and the exponent contributes a quadratic form in an inverse covariance matrix, so maximizing the likelihood over that covariance matrix requires derivatives of \(\log \det A\) and of \(A^{-1}\) with respect to the entries of \(A\). The same two objects appear in the Jacobian term \(\log \lvert I - \rho W \rvert\) of a spatial autoregressive model and in the concentrated likelihood of a seemingly unrelated regressions system.

The determinant

Let \(A\) be \(n \times n\) and invertible, and expand the determinant in cofactors along row \(k\),

\[ \det A = \sum_{l=1}^{n} A_{kl} C_{kl} , \qquad C_{kl} = (-1)^{k+l} M_{kl} , \]

where the minor \(M_{kl}\) is the determinant of the matrix obtained by deleting row \(k\) and column \(l\) from \(A\). The point of the expansion is that none of the cofactors \(C_{k1}, \dots, C_{kn}\) contains any entry of row \(k\), since row \(k\) is exactly what was deleted. The determinant is therefore a linear function of the entries of row \(k\) with the cofactors as coefficients, and differentiating a linear function returns its coefficient:

\[ \frac{\partial \det A}{\partial A_{kl}} = C_{kl} . \]

The cofactors are related to the inverse through the adjugate. By definition \(\operatorname{adj}(A)_{lk} = C_{kl}\), and for invertible \(A\) the adjugate satisfies \(A^{-1} = \operatorname{adj}(A) / \det A\), so \(C_{kl} = \det(A) \, (A^{-1})_{lk}\). The entry in row \(k\) and column \(l\) of the derivative is thus the entry in row \(l\) and column \(k\) of \(A^{-1}\), scaled by the determinant, and the matrix with that entry pattern is the transposed inverse:

\[ \frac{\partial \det A}{\partial A} = \det(A) \left( A^{-1} \right)^{T} . \]

This is Jacobi’s formula. Applying the chain rule for the scalar logarithm removes the leading factor,

\[ \frac{\partial \log \det A}{\partial A} = \frac{1}{\det A} \frac{\partial \det A}{\partial A} = \left( A^{-1} \right)^{T} , \]

which for symmetric \(A\), the case of interest for covariance matrices, is simply \(A^{-1}\). The result is the matrix analogue of \(\mathrm{d} \log a / \mathrm{d}a = 1/a\), with the reciprocal replaced by the inverse and a transpose inserted by the layout convention. In differential form the same statement reads \(\mathrm{d} \log \det A = \operatorname{tr}(A^{-1} \mathrm{d}A)\).

The inverse

The derivative of the inverse follows from differentiating an identity rather than from an expansion. Since \(A A^{-1} = I\) holds for every invertible \(A\), and the right-hand side is constant, the product rule gives

\[ 0 = \mathrm{d} \left( A A^{-1} \right) = \left( \mathrm{d}A \right) A^{-1} + A \, \mathrm{d} \left( A^{-1} \right) , \]

and solving for the second term by multiplying on the left by \(A^{-1}\),

\[ \mathrm{d} \left( A^{-1} \right) = - A^{-1} \left( \mathrm{d}A \right) A^{-1} . \]

This is the analogue of \(\mathrm{d}(1/a)/\mathrm{d}a = -1/a^{2}\), with the square replaced by a product that keeps the two inverses on either side of the perturbation, an order that cannot be rearranged because matrices do not commute. Componentwise, taking \(\mathrm{d}A\) to be the matrix with a one in position \((k,l)\) and zeros elsewhere,

\[ \frac{\partial (A^{-1})_{ij}}{\partial A_{kl}} = - (A^{-1})_{ik} (A^{-1})_{lj} . \]

The combination needed for likelihoods is the trace of an inverse against a constant matrix, \(f(A) = \operatorname{tr}(A^{-1}S)\). Using linearity of the trace, the result above, and the cyclic property \(\operatorname{tr}(XYZ) = \operatorname{tr}(ZXY)\),

\[ \mathrm{d}f = \operatorname{tr}\left( \mathrm{d}(A^{-1}) \, S \right) = - \operatorname{tr}\left( A^{-1} \left( \mathrm{d}A \right) A^{-1} S \right) = - \operatorname{tr}\left( A^{-1} S A^{-1} \, \mathrm{d}A \right) . \]

The identification rule \(\mathrm{d}f = \operatorname{tr}((\partial f / \partial A)^{T} \mathrm{d}A)\), stated in the next section, then gives

\[ \frac{\partial}{\partial A} \operatorname{tr}\left( A^{-1} S \right) = - \left( A^{-1} S A^{-1} \right)^{T} , \]

which reduces to \(-A^{-1}SA^{-1}\) when \(A\) and \(S\) are both symmetric.

The Gaussian likelihood

The two results combine in the derivation that motivates them. Let \(y_{1}, \dots, y_{N}\) be independent draws from a \(p\)-variate normal distribution with mean \(\mu\) and covariance matrix \(\Sigma\). The log-likelihood is

\[ \ell(\mu, \Sigma) = - \frac{Np}{2} \log (2\pi) - \frac{N}{2} \log \det \Sigma - \frac{1}{2} \sum_{i=1}^{N} (y_{i} - \mu)^{T} \Sigma^{-1} (y_{i} - \mu) . \]

The sum in the last term can be rewritten as a trace, which is the manipulation that makes the matrix derivative available. Each summand is a scalar and therefore equals its own trace, so the cyclic property moves the two vectors together:

\[ \sum_{i=1}^{N} (y_{i} - \mu)^{T} \Sigma^{-1} (y_{i} - \mu) = \sum_{i=1}^{N} \operatorname{tr}\left( \Sigma^{-1} (y_{i} - \mu)(y_{i} - \mu)^{T} \right) = N \operatorname{tr}\left( \Sigma^{-1} S(\mu) \right) , \]

where \(S(\mu) = N^{-1} \sum_{i} (y_{i} - \mu)(y_{i} - \mu)^{T}\). The log-likelihood is then

\[ \ell(\mu, \Sigma) = - \frac{Np}{2} \log (2\pi) - \frac{N}{2} \log \det \Sigma - \frac{N}{2} \operatorname{tr}\left( \Sigma^{-1} S(\mu) \right) , \]

and the two terms that depend on \(\Sigma\) are exactly the two derivatives just obtained. Since \(\Sigma\) and \(S\) are symmetric,

\[ \frac{\partial \ell}{\partial \Sigma} = - \frac{N}{2} \Sigma^{-1} + \frac{N}{2} \Sigma^{-1} S(\mu) \Sigma^{-1} = \frac{N}{2} \Sigma^{-1} \left( S(\mu) - \Sigma \right) \Sigma^{-1} . \]

Setting this to zero and multiplying on both sides by \(\Sigma\) leaves \(\hat{\Sigma} = S(\hat{\mu})\). The mean is obtained from the other score, using the Jacobian chain rule with \(g(\mu) = y_{i} - \mu\), whose Jacobian is \(-I_{p}\), and \(\partial (u^{T} \Sigma^{-1} u) / \partial u = 2 \Sigma^{-1} u\):

\[ \frac{\partial \ell}{\partial \mu} = \Sigma^{-1} \sum_{i=1}^{N} (y_{i} - \mu) = 0 \quad \Longrightarrow \quad \hat{\mu} = \bar{y} , \]

since \(\Sigma^{-1}\) is invertible. The maximum likelihood estimator of the covariance matrix is therefore the sample second moment about the sample mean, \(\hat{\Sigma} = N^{-1} \sum_{i} (y_{i} - \bar{y})(y_{i} - \bar{y})^{T}\), with divisor \(N\) rather than \(N - 1\), which is the source of its downward bias.

The symmetry caveat of the previous section applies to the score for \(\Sigma\) and turns out to be harmless. Imposing \(\Sigma_{kl} = \Sigma_{lk}\) doubles each off-diagonal partial derivative and leaves the diagonal ones unchanged, so it rescales the entries of the gradient by positive constants. The set of points at which the gradient vanishes is unaffected, and \(\hat{\Sigma} = S(\bar{y})\) stands under either convention.

The same results with differentials

Once the componentwise derivations are understood, a faster route is available that avoids indices. The first-order change in a scalar \(f\) can be written for a vector argument as

\[ \mathrm{d}f = \left( \frac{\partial f}{\partial x} \right)^{T} \mathrm{d}x , \]

and for a matrix argument as

\[ \mathrm{d}f = \sum_{k} \sum_{l} \frac{\partial f}{\partial A_{kl}} \, \mathrm{d}A_{kl} = \mathrm{tr}\left( \left( \frac{\partial f}{\partial A} \right)^{T} \mathrm{d}A \right) . \]

Both statements are identification rules: compute \(\mathrm{d}f\) by ordinary differentiation, write it in one of the two forms above, and read off the derivative as whatever multiplies \(\mathrm{d}x\) or \(\mathrm{d}A\).

For the quadratic form in \(x\), the product rule gives \(\mathrm{d}f = \mathrm{d}x^{T} A x + x^{T} A \, \mathrm{d}x\). The first term is a scalar and equals its transpose \(x^{T}A^{T}\mathrm{d}x\), so

\[ \mathrm{d}f = \left( x^{T} A^{T} + x^{T} A \right) \mathrm{d}x = \left( \left( A + A^{T} \right) x \right)^{T} \mathrm{d}x , \]

and the gradient is \((A + A^{T})x\). For the least squares objective, with \(r = Ax - b\) and \(\mathrm{d}r = A \, \mathrm{d}x\),

\[ \mathrm{d}f = \mathrm{d}\left( r^{T} r \right) = 2 r^{T} \mathrm{d}r = 2 r^{T} A \, \mathrm{d}x = \left( 2 A^{T} r \right)^{T} \mathrm{d}x , \]

so the gradient is \(2A^{T}(Ax - b)\). For the matrix case, \(x\) is constant and \(\mathrm{d}f = x^{T} \mathrm{d}A \, x\), which is a scalar and therefore equal to its own trace. Using the cyclic property of the trace,

\[ \mathrm{d}f = \mathrm{tr}\left( x^{T} \mathrm{d}A \, x \right) = \mathrm{tr}\left( x x^{T} \mathrm{d}A \right) , \]

and since \(xx^{T}\) is symmetric this matches \(\mathrm{tr}((\partial f / \partial A)^{T} \mathrm{d}A)\) with \(\partial f / \partial A = x x^{T}\). All three agree with the componentwise results.

The determinant is the case in which the differential route is clearly the shorter one. Jacobi’s formula \(\mathrm{d} \det A = \det(A) \operatorname{tr}(A^{-1} \mathrm{d}A)\) divided by \(\det A\) gives \(\mathrm{d} \log \det A = \operatorname{tr}(A^{-1} \mathrm{d}A)\), and the identification rule reads off \(\partial \log \det A / \partial A = (A^{-1})^{T}\) in one line, against a cofactor expansion in several. The derivative of the inverse and the trace result of the previous section were obtained this way already, since no componentwise route to them is any shorter.

Conclusion

The results can be collected as follows, with \(x\) an \(n \times 1\) vector, \(a\) an \(n \times 1\) constant vector, \(b\) an \(m \times 1\) constant vector, \(A\) and \(S\) matrices of the indicated size, \(g\) a vector-valued function with Jacobian \(J = \partial g / \partial x^{T}\), and \(W\) a constant symmetric weight matrix.

Function	Differentiated with respect to	Derivative
\(a^{T}x\)	\(x\)	\(a\)
\(x^{T}Ax\), \(A\) of size \(n \times n\)	\(x\)	\((A + A^{T})x\), or \(2Ax\) if \(A = A^{T}\)
\(\lVert Ax - b \rVert_{2}^{2}\), \(A\) of size \(m \times n\)	\(x\)	\(2A^{T}(Ax - b)\)
\(f(g(x))\), \(f\) scalar	\(x\)	\(J(x)^{T} \, \partial f / \partial u\) at \(u = g(x)\)
\(g(x)^{T} W g(x)\)	\(x\)	\(2 J(x)^{T} W g(x)\)
\(x^{T}Ax\), \(A\) of size \(n \times n\)	\(A\)	\(xx^{T}\)
\(\log \det A\), \(A\) invertible	\(A\)	\((A^{-1})^{T}\), or \(A^{-1}\) if \(A = A^{T}\)
\(\operatorname{tr}(A^{-1}S)\), \(A\) invertible	\(A\)	\(-(A^{-1} S A^{-1})^{T}\)

Almost none of these required a rule that does not already hold for functions of one variable. The linear form differentiates to its coefficient, the quadratic form to twice its matrix once the matrix has been symmetrized, the composition of a square with a linear map picks up the transpose of that map by the chain rule, and the logarithm of a determinant differentiates to an inverse in the same way that \(\log a\) differentiates to \(1/a\). The two places where genuinely matrix-specific reasoning entered are the cofactor expansion behind Jacobi’s formula and the differentiation of the identity \(AA^{-1} = I\), and even the second of those is the familiar trick of differentiating both sides of a relation that holds identically.

What matrix notation contributes is a way of arranging \(n\) or \(n^{2}\) partial derivatives so that the arrangement is itself a matrix expression. That is what makes a first-order condition such as \(A^{T}(Ax - b) = 0\) readable as a statement about orthogonality rather than as a list of equations, and what allows the maximum likelihood estimator of a covariance matrix to be read off from \(\Sigma^{-1}(S - \Sigma)\Sigma^{-1} = 0\) in a single step. When a result is in doubt, the componentwise derivation settles it, and the \(n = 2\) case is usually enough to see which of two candidate expressions is correct.