A gentle introduction to matrix calculus

Jan R. Magnus

Abstract

Matrix calculus is an important tool when we wish to optimize functions involving matrices or perform sensitivity analyses. This tutorial is designed to make matrix calculus more accessible to graduate students and young researchers. It contains the theory that would suffice in most applications, many fully worked-out exercises and examples, and presents some of the ‘tacit knowledge’ that is prevalent in this field.

1. Introduction

Section Summary: This tutorial aims to make matrix calculus accessible to graduate students and young researchers by distilling its essential ideas into a practical, self-contained guide rather than a lengthy theoretical text. It emphasizes hands-on techniques, worked solutions to exercises, useful shortcuts, and common pitfalls, while building on two core ideas: a clear definition of the matrix derivative and the advantages of working with differentials, which preserve the original dimensions of the functions involved. The paper then outlines its structure, moving from basic tools like traces and Kronecker products through optimization examples to more advanced applications.

The purpose of this paper is to make matrix calculus more accessible to graduate students and young researchers. It is written as a stand-alone tutorial with minimal references to underlying results, so that it can be used in an advanced undergraduate course and as an aid to graduate and Ph.D. students in one of their courses or projects.

Matrix calculus is much used in economics, statistics, mathematics, psychology, control theory, engineering, and elsewhere, and it is required whenever we wish to optimize functions involving matrices or perform sensitivity analyses. Heinz Neudecker and I started to work on our monograph Matrix Differential Calculus with Applications in Statistics and Econometrics in the early 1980s, and the book came out in 1988. It went through several revisions and new editions, and I regard the 2019 third edition, which I wrote after Heinz’ death, as the definitive text.

Not everyone is willing to struggle through a 450-page monograph, and it is for those people that the current tutorial may be of some use. It contains the essence of matrix calculus, leaving out the subtleties that are not used in most applications, such as how to differentiate eigenvalues or Moore–Penrose inverses. Obviously, the tutorial provides less information than the monograph, but in two directions it provides more. First, the monograph contains many exercises, but there exists no answer book, so there are no model answers. The current tutorial also contains many exercises, but they all come with fully worked-out solutions. Second, while the emphasis in the monograph is primarily on the theory, the emphasis in this tutorial is primarily on the practice, and I have included suggestions on shortcuts, warnings for pitfalls, and certain ‘tricks’ that are useful to know when doing matrix calculus in practice.

Matrix calculus rests on two pillars and it requires six tools. The tools are presented in Section 2 and consist of some basic results on the trace and linear and quadratic functions, some less basic results on the Kronecker product and the vec operator, and finally some more advanced results on the commutation and duplication matrix. Almost all these results will be proved using only elementary mathematics.

The two pillars are discussed in Sections 3, 4. The first pillar is the definition of a matrix derivative, which I present in Section 3. Suppose that we have a matrix function, say $F (X) = X^{- 1}$. The derivative $D F (X)$ will be a matrix, but how are the elements in this matrix to be organized? This seemingly simple question was not fully resolved until the mid-1980s. I will show that there exists only one correct definition of a matrix derivative.

Section 4 introduces the concept of a differential, the second pillar on which matrix calculus rests. The key advantage of the differential over the more common derivative is the following. Consider an $m \times 1$ vector function $f$, for example $f (x) = A x$ where $A$ is an $m \times n$ matrix of constants. Then the derivative $D f (x)$ is an $m \times n$ matrix, in this case the matrix $A$. But the differential $d f (x)$ remains an $m \times 1$ vector. The differential $d f (x)$ has the same dimension as $f$, irrespective of the dimension of the vector $x$. The advantage is even larger for matrices. The differential $d F (X)$ of a matrix function $F (X)$ has the same dimension as $F$, irrespective of the dimension of the matrix $X$. The practical importance of working with differentials will be demonstrated through many examples.

I discuss optimization in Section 5 and this leads to our first example (least squares) in Section 6, both with and without constraints. In Section 7, I introduce matrix calculus, showing how the previous results can be straightforwardly generalized from vector calculus to matrix calculus, at least if one employs the correct definition of matrix derivative. Section 8 contains exercises, and Sections 9–11 are devoted to the second differential with associated exercises in Section 12. Sections 13–15 contain three further (more advanced) examples on how to apply matrix calculus. In the concluding Section 16, I provide some hints from my own experience in an attempt to transform tacit (unwritten) knowledge into explicit (written) knowledge.

The following notation is used. Lower-case symbols ($a$, $x$) denote scalars or vectors, upper-case symbols ($A$, $X$) denote matrices. Thus, $f$ denotes a scalar or vector function, and $F$ a matrix function. I write $A^{\prime}$, $A^{- 1}$, $\operatorname{tr} A$, $| A |$ for the transpose, inverse, trace, and determinant of $A$. All functions and variables are real. Parentheses are used sparingly. I write $d X$, $\operatorname{tr} X$, and $\operatorname{vec} X$ without parentheses, and also $d X Y$, $\operatorname{tr} X Y$, and $\operatorname{vec} X Y$ instead of $d (X Y)$, $\operatorname{tr} (X Y)$, and $\operatorname{vec} (X Y)$. However, I write $\operatorname{vech} (X)$ with parentheses for historical reasons.

2. The tools

Section Summary: The section introduces three fundamental tools for matrix calculus. The trace sums the diagonal entries of a square matrix and obeys useful identities such as invariance under transposition and cyclic permutation of products. Linear forms Ax and quadratic forms x′Ax are examined to determine when they vanish for all vectors, revealing that this occurs precisely when A is the zero matrix or, for quadratics, when A is skew-symmetric; the Kronecker product is then defined as a block-wise scaling operation that multiplies matrices of any compatible dimensions and satisfies associative, distributive, and mixed-product rules.

We present six tools that are indispensable in matrix calculus.

2.1. The trace

An important function of a square matrix $A = (a_{i j})$ is the trace, defined as the sum of its diagonal elements: $\operatorname{tr} A = \sum_{i} a_{i i}$. Since a matrix and its transpose contain the same diagonal elements, we have

$ \operatorname{tr} A^{\prime} = \operatorname{tr} A. \tag{1} $

Somewhat less trivial is

Proposition 1

For any two matrices $A$ and $B$ of the same order,

$ \operatorname{tr} A^{\prime} B = \operatorname{tr} B A^{\prime}. $

Proof

This follows because

$ \begin{aligned} \operatorname{tr} A^{\prime} B &= \sum_{j} (A^{\prime} B){j j} = \sum{j} \sum_{i} a_{i j} b_{i j} \ &= \sum_{i} \sum_{j} b_{i j} a_{i j} = \sum_{i} (B A^{\prime})_{i i} = \operatorname{tr} B A^{\prime}. \quad \square \end{aligned} $

The proposition implies that $\operatorname{tr} A B C = \operatorname{tr} C A B = \operatorname{tr} B C A$, because the proposition allows cyclical permutations of the matrices. But the proposition does not imply that $\operatorname{tr} A B C = \operatorname{tr} A C B$, which is not cyclical.

2.2. Linear and quadratic forms

A linear form is an expression such as $A x$. When $A x = 0$, this does not imply that either $A = 0$ or $x = 0$. For example, if

$ A = \left( \begin{array}{cc} \phantom{-} 1 & - 1 \ - 2 & \phantom{-} 2 \ \phantom{-} 3 & - 3 \end{array} \right), \quad x = \left( \begin{array}{c} 1 \ 1 \end{array} \right), $

then $A x = 0$, but neither $A = 0$ nor $x = 0$. However, when $A x = 0$ for every $x$, then $A$ must be zero.

Things are different with a quadratic form, that is, an expression such as $x^{\prime} A x$. When $x^{\prime} A x = 0$, this does not imply that $A = 0$ or $x = 0$ or $A x = 0$. Even when $x^{\prime} A x = 0$ for every $x$, it does not follow that $A = 0$, as the example

$ A = \left( \begin{array}{cc} \phantom{-} 0 & 1 \ - 1 & 0 \end{array} \right) $

demonstrates. This matrix is skew-symmetric, that is, it satisfies $A^{\prime} = - A$. In fact, when $x^{\prime} A x = 0$ for every $x$ then it follows that $A$ must be skew-symmetric.

Proposition 2

We have

• (i) $A x = 0$ for every $x$ $: ⟺ :$ $A = 0$,
• (ii) $x^{\prime} A x = 0$ for every $x$ $: ⟺ :$ $A$ is skew-symmetric,
• (iii) $x^{\prime} A x = 0$ for every $x$ and $A = A^{\prime}$ $: ⟺ :$ $A = 0$.

Proof

We shall prove only one direction of these equivalences, the other being trivial. Let $e_{i}$ denote the $i$th elementary vector, that is the vector with one in the $i$th position and zeros elsewhere. To prove (i), let $x = e_{i}$. Then, $A e_{i} = 0$ for every $i$, that is, every column of $A$ is zero, and hence $A = 0$. To prove (ii), let $x = e_{i} + e_{j}$. Then,

$ a_{i i} + a_{i j} + a_{j i} + a_{j j} = 0, $

implying in particular, by setting $i = j$, that $a_{i i} = a_{j j} = 0$. Hence, $a_{i j} + a_{j i} = 0$, as requested. To prove (iii), note that if $A$ is symmetric and skew-symmetric, then $A = - A$, implying that $A$ must be the null matrix. □

Proposition 2(i) has important practical consequences. Suppose we wish to prove that $A = B$. Then we could try and prove that $a_{i j} = b_{i j}$ for all $i$ and $j$. But it is often easier to show that $A x = B x$ for arbitrary $x$. Instead of showing directly that $A = B$, we show that $A$ and $B$ have the same effect when working on an arbitrary vector $x$.

Proposition 2(ii) also has important consequences. It tells us that if $x^{\prime} A x = x^{\prime} B x$ for every $x$, then it is not necessarily true that $A = B$, but it is true that $A + A^{\prime} = B + B^{\prime}$. If $A$ is symmetric and $x^{\prime} A x = x^{\prime} B x$ for every $x$, then it still does not follow that $A = B$, but it does follow that $A = (B + B^{\prime}) / 2$. This will be important when we discuss the second identification theorem in Proposition 13 and Eq. (42).

2.3. The Kronecker product

Let $A$ be an $m \times n$ matrix and $B$ a $p \times q$ matrix. The $m p \times n q$ matrix defined by

$ \left( \begin{array}{ccc} a_{11} B & \ldots & a_{1 n} B \ \vdots & & \vdots \ a_{m 1} B & \ldots & a_{m n} B \end{array} \right) \tag{2} $

is called the Kronecker product of $A$ and $B$ and is written as $A \otimes B$. The Kronecker product $A \otimes B$ is thus defined for any pair of matrices $A$ and $B$, unlike the matrix product $A B$ which exists only if the number of columns in $A$ equals the number of rows in $B$ or if either $A$ or $B$ is a scalar.

Proposition 3

The following three properties show that the Kronecker product is in fact a product:

$ \begin{aligned} (A \otimes B) \otimes C &= A \otimes (B \otimes C), \ &(A + B) \otimes (C + D) = A \otimes C + A \otimes D + B \otimes C + B \otimes D, \ &(A \otimes B) (C \otimes D) = A C \otimes B D, \end{aligned} $

where it is assumed that $A + B$ and $C + D$ are defined in the second equality, and that $A C$ and $B D$ are defined in the third equality.

Proof

These equalities are easiest proved by direct application of the definition, which takes quite a bit of space and is therefore omitted. □

Exercise 1

For any two column vectors $a$ and $b$ (not necessarily of the same order), we have

$ a \otimes b^{\prime} = a b^{\prime} = b^{\prime} \otimes a. $

Solution. This follows, in essence, from the fact that a vector $x$ can be written as $x \otimes 1$ and also as $1 \otimes x$. Define elementary vectors $e_{i}$ ($m \times 1$) and $u_{j}$ ($n \times 1$), where $m$ and $n$ are the orders of $a$ and $b$, respectively. Then,

$ \begin{aligned} e_{i}^{\prime} (a \otimes b^{\prime}) u_{j} = (e_{i}^{\prime} \otimes 1) (a \otimes b^{\prime}) (1 \otimes u_{j}) = (e_{i}^{\prime} a) \otimes (b^{\prime} u_{j}) = a_{i} b_{j} = (a b^{\prime}){i j}, \ e{i}^{\prime} (b^{\prime} \otimes a) u_{j} = (1 \otimes e_{i}^{\prime}) (b^{\prime} \otimes a) (u_{j} \otimes 1) = (b^{\prime} u_{j}) \otimes (e_{i}^{\prime} a) = b_{j} a_{i} = (a b^{\prime})_{i j}, \end{aligned} $

using Proposition 3. $; \parallel$

The transpose, trace, and inverse of the Kronecker product are given in the next proposition.

Proposition 4

We have

$ \begin{aligned} \begin{aligned} (A \otimes B)^{\prime} &= A^{\prime} \otimes B^{\prime}, \ &\operatorname{tr} (A \otimes B) = (\operatorname{tr} A) (\operatorname{tr} B), \end{aligned} \ (A \otimes B)^{- 1} = A^{- 1} \otimes B^{- 1}, \end{aligned} $

where $A$ and $B$ are assumed to be square in the second equality and nonsingular in the third equality.

Proof

The first two results follow again by direct application of the definition. The third result follows from Proposition 3 by writing

$ (A \otimes B) (A^{- 1} \otimes B^{- 1}) = (A A^{- 1}) \otimes (B B^{- 1}) = I_{m} \otimes I_{n} = I_{m n}, $

where $A$ has order $m \times m$, and $B$ has order $n \times n$. □

2.4. The vec operator

Consider an $m \times n$ matrix $A$. This matrix has $n$ columns, say $a_{1}, \ldots, a_{n}$. We can transform the matrix into a vector by defining the $m n \times 1$ vector $\operatorname{vec} A$ as the vector which stacks the columns of $A$ one underneath the other:

$ \operatorname{vec} A = \left( \begin{array}{c} a_{1} \ a_{2} \ \vdots \ a_{n} \end{array} \right). \tag{3} $

The $\operatorname{vec}$ operator is not the only operator that transforms a matrix into a vector. We could also take the rows of $A$ and stack them one underneath the other, but the $\operatorname{vec}$ transformation is now the most common.

Exercise 2

If

$ A = \left( \begin{array}{ccc} 1 & 2 & 3 \ 4 & 5 & 6 \end{array} \right), $

then what is $\operatorname{vec} A$? And what is $\operatorname{vec} A^{\prime}$?

Solution. We have $\operatorname{vec} A = (1, 4, 2, 5, 3, 6)^{\prime}$ and $\operatorname{vec} A^{\prime} = (1, 2, 3, 4, 5, 6)^{\prime}$. $; \parallel$

The $\operatorname{vec}$ operator and the trace are related through the following result.

Proposition 5

If $A$ and $B$ are matrices of the same order, then

$ (\operatorname{vec} A)^{\prime} (\operatorname{vec} B) = \operatorname{tr} A^{\prime} B. $

Proof

The product $(\operatorname{vec} A)^{\prime} (\operatorname{vec} B)$ multiplies the corresponding elements in $A$ and $B$ and adds them together, so that $(\operatorname{vec} A)^{\prime} (\operatorname{vec} B) = \sum_{i j} a_{i j} b_{i j}$. But this is precisely the expression for $\operatorname{tr} A^{\prime} B$ by the proof of Proposition 1. □

Exercise 3

If $a$ and $b$ are vectors of arbitrary order, show that

$ \operatorname{vec} a b^{\prime} = b \otimes a. $

Solution. The $j$th column of $a b^{\prime}$ is given by $b_{j} a$. $; \parallel$

The previous exercise relates the Kronecker product and the vec operator. Its generalization is an important and frequently used tool.

Proposition 6

For any matrices $A$, $B$, and $C$ for which the product $A B C$ is defined, we have

$ \operatorname{vec} A B C = (C^{\prime} \otimes A) \operatorname{vec} B. $

Proof

Using Exercise 3, we see that

$ \begin{aligned} \operatorname{vec} (A b e^{\prime} C) &= \operatorname{vec} \left( (A b) (C^{\prime} e)^{\prime} \right) = (C^{\prime} e) \otimes (A b) \ &= (C^{\prime} \otimes A) (e \otimes b) = (C^{\prime} \otimes A) \operatorname{vec} (b e^{\prime}) \end{aligned} $

for any vectors $b$ and $e$. Then, writing $B = \sum_{j} b_{j}^{} e_{j}^{\prime}$ where $b_{j}$ and $e_{j}$ denote the $j$th column of $B$ and $I$, respectively, the result follows. □

2.5. The commutation matrix

Let $A$ be an $m \times n$ matrix. The vectors $\operatorname{vec} A$ and $\operatorname{vec} A^{\prime}$ contain the same $m n$ elements, but in a different order. Hence, there exists a unique $m n \times m n$ matrix, which transforms $\operatorname{vec} A$ into $\operatorname{vec} A^{\prime}$. This matrix contains $m n$ ones and $m n (m n - 1)$ zeros and is called the commutation matrix, denoted by $K_{m n}$. Thus,

$ K_{m n} \operatorname{vec} A = \operatorname{vec} A^{\prime}. \tag{4} $

If $m = n$, we write $K_{n}$ instead of $K_{n n}$. The commutation matrix is a square matrix with exactly one entry of $1$ in each row and each column with all other entries $0$. (Such matrices are called permutation matrices.)

Exercise 4

Show that $\operatorname{tr} K_{n} = n$.

Solution. If $A$ is a square matrix of order $n \times n$, then its $i$th diagonal element is placed in the same position in $\operatorname{vec} A$ as in $\operatorname{vec} A^{\prime}$. Hence, the $i$th diagonal block of $K_{n}$ has 1 in the $i i$th position and zeros elsewhere. Summing up the diagonal elements of $K_{n}$ thus adds 1 for each block, that is, $n$ in total. $; \parallel$

Proposition 7

The commutation matrix satisfies

$ K_{m n}^{\prime} = K_{m n}^{- 1} = K_{n m}. $

Proof

All permutation matrices are orthogonal, hence $K_{m n}^{\prime} = K_{m n}^{- 1}$. Also, premultiplying (4) by $K_{n m}$ gives $K_{n m} K_{m n} \operatorname{vec} A = \operatorname{vec} A$, which shows that $K_{n m} K_{m n} = I_{m n}$. □

The key property of the commutation matrix enables us to interchange (commute) the two matrices of a Kronecker product.

Proposition 8

$ K_{p m} (A \otimes B) = (B \otimes A) K_{q n} $

for any $m \times n$ matrix $A$ and $p \times q$ matrix $B$.

Proof

This is easiest shown, not by proving a matrix identity but by proving that the effect of the two matrices on an arbitrary vector is the same, in the spirit of Proposition 2(i). Thus, let $X$ be an arbitrary $q \times n$ matrix. Then, by repeated application of (4) and Proposition 6,

$ \begin{aligned} K_{p m} (A \otimes B) \operatorname{vec} , X &= K_{p m} \operatorname{vec} , B X A^{\prime} = \operatorname{vec} , A X^{\prime} B^{\prime} \ &= (B \otimes A) \operatorname{vec} , X^{\prime} = (B \otimes A) K_{q n} \operatorname{vec} , X. \end{aligned} $

Since $X$ is arbitrary, the result follows. □

Exercise 5

Let $N_{n} = (I_{n^{2}} + K_{n}) / 2$. Show that $N_{n}$ is symmetric idempotent and that $N_{n} (A \otimes A) = (A \otimes A) N_{n}$ for every $n \times n$ matrix $A$.

Solution. Since $K_{n}$ is symmetric and orthogonal by Proposition 7, we have $K_{n}^{2} = I_{n^{2}}$ and hence $N_{n}^{} = N_{n}^{\prime} = N_{n}^{2}$. Also, using Proposition 8, $K_{n} (A \otimes A) = (A \otimes A) K_{n}$ and hence $N_{n} (A \otimes A) = (A \otimes A) N_{n}$. $; \parallel$

2.6. The duplication matrix

Many matrices in statistics and econometrics are symmetric. When we differentiate with respect to symmetric matrices, we must take the symmetry into account and we need the duplication matrix.

Let $A$ be a square $n \times n$ matrix. Then $\operatorname{vech} (A)$ (the ‘vec-half’ operator) denotes the $\frac{1}{2} n (n + 1) \times 1$ vector that is obtained from $\operatorname{vec} A$ by eliminating all elements of $A$ above the diagonal. For example, when $n = 3$, crossing out the elements above the diagonal,

$ \left( \begin{array}{ccc} a_{11} & \cancel{a_{12}} & \cancel{a_{13}} \ a_{21} & a_{22} & \cancel{a_{23}} \ a_{31} & a_{32} & a_{33} \end{array} \right), $

and vectorizing the remaining elements, we obtain

$ \operatorname{vech} (A) = (a_{11}, a_{21}, a_{31}, a_{22}, a_{32}, a_{33})^{\prime}. $

In this way, for symmetric $A, \operatorname{vech} (A)$ contains only the generically distinct elements of $A$. Since the elements of $\operatorname{vec} A$ are those of $\operatorname{vech} (A)$ with some repetitions, there exists a unique $n^{2} \times \frac{1}{2} n (n + 1)$ matrix which transforms, for symmetric $A, \operatorname{vech} (A)$ into $\operatorname{vec} A$. This matrix is called the duplication matrix and is denoted by $D_{n}$. Thus,

$ D_{n} \operatorname{vech} (A) = \operatorname{vec} A \quad (A = A^{\prime}). \tag{5} $

The matrix $D_{n}$ has full column rank $\frac{1}{2} n (n + 1)$, so that $D_{n}^{\prime} D_{n}^{}$ is nonsingular. This implies that $\operatorname{vech} (A)$ can be uniquely solved from (5), and we have

$ \operatorname{vech} (A) = (D_{n}^{\prime} D_{n}^{})^{- 1} D_{n}^{\prime} \operatorname{vec} A \quad (A = A^{\prime}). \tag{6} $

Proposition 9

The duplication matrix is connected to the commutation matrix by

$ K_{n} D_{n} = D_{n}, \quad D_{n}^{} (D_{n}^{\prime} D_{n}^{})^{- 1} D_{n}^{\prime} = \frac{1}{2} (I_{n^{2}} + K_{n}). $

Proof

Let $X$ be a symmetric $n \times n$ matrix. Then,

$ K_{n} D_{n} \operatorname{vech} (X) = K_{n} \operatorname{vec} X = \operatorname{vec} X = D_{n} \operatorname{vech} (X). $

The symmetry of $X$ does not restrict $\operatorname{vech} (X)$, which is therefore arbitrary. Hence, the first result follows. To prove the second result, let

$ M_{n} = D_{n}^{} (D_{n}^{\prime} D_{n}^{})^{- 1} D_{n}^{\prime}, \quad N_{n} = \frac{1}{2} (I_{n^{2}} + K_{n}), \quad \Delta_{n} = M_{n} - N_{n}. $

Both $M_{n}$ and $N_{n}$ are symmetric idempotent, the latter by Exercise 5. Since $K_{n} D_{n} = D_{n}$, we have $M_{n} N_{n} = N_{n} M_{n} = M_{n}$, so that $\Delta_{n}$ is also symmetric idempotent. Now,

$ \operatorname{tr} M_{n} = \operatorname{tr} I_{n (n + 1) / 2} = n (n + 1) / 2 = (n^{2} + n) / 2 = \operatorname{tr} N_{n}, $

since $\operatorname{tr} K_{n} = n$ (Exercise 4). This gives $r (\Delta_{n}) = \operatorname{tr} \Delta_{n} = 0$, and hence $\Delta_{n} = 0$. □

Much of the interest in the duplication matrix is due to the importance of the matrix $D_{n}^{\prime} (A \otimes A) D_{n}^{}$, caused by the fact that, for symmetric $X$, the scalar expression $\operatorname{tr} A X A^{\prime} X$ occurs frequently in statistics and econometrics; see Exercise 28 and Section 13.

Proposition 10

Let $A$ be an $n \times n$ matrix, not necessarily symmetric. The matrix $D_{n}^{\prime} (A \otimes A) D_{n}^{}$ satisfies the following properties:

$ \begin{aligned} D_{n}^{} (D_{n}^{\prime} D_{n}^{})^{- 1} D_{n}^{\prime} (A \otimes A) D_{n}^{} = (A \otimes A) D_{n}^{}, \ \begin{aligned} (D_{n}^{\prime} (A \otimes A) D_{n}^{})^{- 1} = (D_{n}^{\prime} D_{n}^{})^{- 1} D_{n}^{\prime} (A^{- 1} \otimes A^{- 1}) D_{n}^{} (D_{n}^{\prime} D_{n}^{})^{- 1}, \ | D_{n}^{\prime} (A \otimes A) D_{n}^{} | = 2^{\frac{1}{2} n (n - 1)} | A |^{n + 1}, \end{aligned} \end{aligned} $

where $A$ is assumed to be nonsingular in the second equality.

Proof

The first result follows from Exercise 5 and $N_{n} D_{n} = D_{n}$. To prove the second result we write

$ \begin{aligned} &D_{n}^{\prime} (A \otimes A) D_{n}^{} (D_{n}^{\prime} D_{n}^{})^{- 1} D_{n}^{\prime} (A^{- 1} \otimes A^{- 1}) D_{n}^{} (D_{n}^{\prime} D_{n}^{})^{- 1} \ &= D_{n}^{\prime} (A \otimes A) (A^{- 1} \otimes A^{- 1}) D_{n}^{} (D_{n}^{\prime} D_{n}^{})^{- 1} \ &= D_{n}^{\prime} D_{n}^{} (D_{n}^{\prime} D_{n}^{})^{- 1} = I_{n (n + 1) / 2}. \end{aligned} $

The third result is more difficult to prove. It follows from

$ | D_{n}^{\prime} (A \otimes A) D_{n}^{} | = | D_{n}^{\prime} D_{n}^{} | | (D_{n}^{\prime} D_{n}^{})^{- 1} D_{n}^{\prime} (A \otimes A) D_{n}^{} |, $

and the fact that $| D_{n}^{\prime} D_{n}^{} | = 2^{\frac{1}{2} n (n - 1)}$ and the eigenvalues of the matrix $(D_{n}^{\prime} D_{n}^{})^{- 1} D_{n}^{\prime} (A \otimes A) D_{n}^{}$ are $\lambda_{i} \lambda_{j}$ ($1 \leq j \leq i \leq n$) where ${\lambda_{i} }$ denote the eigenvalues of $A$; see Magnus (1988, Theorems 4.4 and 4.10). □

3. The definition of matrix derivative

Section Summary: The section explains that the matrix derivative is defined by vectorizing both the output matrix F and the input matrix X, then arranging all partial derivatives into a single matrix where each row shows how one element of F changes with respect to every element of X, and each column shows how every element of F changes with respect to one element of X. This specific layout, rather than preserving the original matrix shapes, is required to keep the essential idea of a derivative intact when moving from vectors to matrices. Any other arrangement of the partial derivatives fails to qualify as a true derivative and should be avoided.

Matrix calculus uses six tools (discussed in the previous section) and it rests on two pillars: the correct definition of a matrix derivative and the concept of a differential. In the current section I discuss how to define a matrix derivative; in the next section I introduce differentials.

Let $x$ ($n \times 1$) and $y$ ($m \times 1$) be two vectors and let $y$ be a function of $x$, say $y = f (x)$. What is the derivative of $y$ with respect to $x$? Let us first consider the linear equation $y = f (x) = A x$, where $A$ is an $m \times n$ matrix of constants. The derivative is $A$ and we write

$ \frac{\partial f (x)}{\partial x^{\prime}} = A. \tag{7} $

The notation $\partial f (x) / \partial x^{\prime}$ is just notation, nothing else. We sometimes write the derivative as $D f (x)$ or as $f^{\prime} (x)$, but we only use the latter notation if it does not cause confusion with the transpose. The proposed notation emphasizes that we differentiate an $m \times 1$ column vector $f$ with respect to a $1 \times n$ row vector $x^{\prime}$, resulting in an $m \times n$ derivative matrix.

More generally, the derivative of $f (x)$ is an $m \times n$ matrix containing all partial derivatives $\partial f_{i} (x) / \partial x_{j}$, but in a specific ordering, namely

$ \frac{\partial f (x)}{\partial x^{\prime}} = \left( \begin{array}{cccc} \partial f_{1} (x) / \partial x_{1} & \partial f_{1} (x) / \partial x_{2} & \ldots & \partial f_{1} (x) / \partial x_{n} \ \partial f_{2} (x) / \partial x_{1} & \partial f_{2} (x) / \partial x_{2} & \ldots & \partial f_{2} (x) / \partial x_{n} \ \vdots & \vdots & & \vdots \ \partial f_{m} (x) / \partial x_{1} & \partial f_{m} (x) / \partial x_{2} & \ldots & \partial f_{m} (x) / \partial x_{n} \end{array} \right). \tag{8} $

There is only one definition of a vector derivative, and this is it. Of course, one can organize the $m n$ partial derivatives in different ways, but these other combinations of the partial derivatives are not derivatives, have no practical use, and should be avoided.

Notice that each row of the derivative in (8) contains the partial derivatives of one element of $f$ with respect to all elements of $x$, and that each column contains the partial derivatives of all elements of $f$ with respect to one element of $x$. This is an essential characteristic of the derivative. As a consequence, the derivative of a scalar function, such as $y = a^{\prime} x$ (where $a$ is a vector of constants), is a row vector; in this case, $a^{\prime}$. So the derivative of $a^{\prime} x$ is $a^{\prime}$, not $a$.

The essence of the concept derivative must be preserved when we generalize from vector functions $f (x)$ to matrix functions $F (X)$. While it is tempting to keep the structure of the matrices $X$ and $F$ intact — and this line of thought was the original idea about matrix derivatives at the end of the 1960s — it was not until the early 1980s that we realized that this is the wrong generalization, since it does not maintain the concept of derivative; see Magnus (2010) for further discussion on the concept of matrix derivative and Magnus (2024) for a sketch of the historical development.

We now understand that when we have a collection of functions and a collection of variables, which may be organized in matrices (such as $F$ and $X$) or otherwise, then what is required is a one-dimensional ordering of the functions and of the variables, for which we choose $\operatorname{vec} F$ and $\operatorname{vec} X$, because these are now the most common vectorizations. In fact, which vectorization we choose is irrelevant, we may even use a different vectorization for $F$ and for $X$, although this would be confusing and is not recommended.

However, once we have chosen our vectorization, we must stick to it, so that the derivative

$ \frac{\partial \operatorname{vec} F (X)}{\partial (\operatorname{vec} X)^{\prime}} \tag{9} $

has the essential property that each row contains the partial derivatives of one element of $F$ with respect to all elements of $X$, and that each column contains the partial derivatives of all elements of $F$ with respect to one element of $X$. This is the only correct definition of a matrix derivative.

4. The differential

Section Summary: The differential of a function associates to each point a linear approximation that captures how small changes in the input affect the output, and it is obtained simply by multiplying the derivative by an increment dx. This object obeys familiar algebraic rules for constants, sums, and products, and it remains unchanged in form when functions are composed, which yields Cauchy's rule of invariance and lets one differentiate complicated expressions step by step. For vector-valued functions the same construction preserves the dimension of the original mapping, offering a practical advantage over the matrix of partial derivatives.

4.1. Definition, notation, and rules

If $f$ is a scalar function of one variable, such as $f (x) = x^{2}$, then we define the differential $d f (x)$ at a point $x$ as the linear function (of $u$)

$ d f (x) (u) = f^{\prime} (x) u, \tag{10} $

where $x$ is a point where the derivative $f^{\prime} (x)$ exists and $u$ is an arbitrary point in $R$. The differential is closely related to the derivative, but the two concepts are not the same: the differential is a geometric concept, while the derivative is an algebraic concept, representing the differential by some number (the slope or the Jacobian).

Exercise 6

Show that

$ \begin{aligned} f (x) = x ⟹ d f (x) (u) = u, \ f (x) = x^{2} ⟹ d f (x) (u) = 2 x u, \ f (x) = \sin (x) ⟹ d f (x) (u) = \cos (x) u, \ f (x) = \sin^{2} (x) ⟹ d f (x) (u) = 2 \sin (x) \cos (x) u. \end{aligned} $

Solution. These results follow directly from the definition in (10), where in the last case we apply the chain rule:

$ \frac{d \sin^{2} (x)}{d x} = \frac{d \sin^{2} (x)}{d \sin (x)} \frac{d \sin (x)}{d x} = 2 \sin (x) \cos (x). \quad \parallel $

Exercise 6 shows that the differential associated with the identity function $f (x) = x$ is given by $d x (u) = u$. Hence, $d f (x) (u) = f^{\prime} (x) , d x (u)$, or for short,

$ d f (x) = f^{\prime} (x) , d x, \tag{11} $

which is the notation we shall use hereafter.

Exercise 7

Show that

$ \begin{aligned} d x^{2} = 2 x , d x, \ d \sin (x) = \cos (x) , d x, \ d \sin^{2} (x) = 2 \sin (x) \cos (x) , d x. \end{aligned} $

Solution. This follows immediately from the previous exercise by replacing $u$ with $d x$. $; \parallel$

The differential is an operator, in fact a linear operator, and it obeys the following simple rules:

$ \begin{aligned} d a = 0, \quad (a \text{constant}), \ d (a x) = a , d x \quad (a \text{constant}), \ d (f (x) + g (x)) = d f (x) + d g (x), \ d (f (x) g (x)) = (d f (x)) , g (x) + f (x) , d g (x). \end{aligned} \tag{12} $

4.2. Cauchy’s rule of invariance

The chain rule, well-known for derivatives, also applies to differentials and is then called Cauchy’s rule of invariance. Cauchy’s rule states that taking differentials of functions preserves composition. Formally, if $z = h (x) = f (g (x))$ is the composition of the functions $y = g (x)$ and $z = f (y)$, then $d y = g^{\prime} (x) , d x$ and $d z = f^{\prime} (y) , d y$, so that

$ d z = f^{\prime} (y) , d y = f^{\prime} (y) g^{\prime} (x) , d x = f^{\prime} (g (x)) , g^{\prime} (x) , d x. \tag{13} $

Exercise 8

Use Cauchy’s rule of invariance to find the differentials of the functions $\sin^{2} (x)$, $e^{x^{2}}$, and $e^{\sin (x^{2})}$.

Solution. We have

$ \begin{aligned} d \sin^{2} (x) = 2 \sin (x) , d \sin (x) = 2 \sin (x) \cos (x) , d x, \ d e^{x^{2}} = e^{x^{2}} d x^{2} = 2 e^{x^{2}} x , d x, \ \begin{aligned} d e^{\sin (x^{2})} &= e^{\sin (x^{2})} , d \sin (x^{2}) = e^{\sin (x^{2})} , \cos (x^{2}) , d x^{2} \ &= 2 x , e^{\sin (x^{2})} , \cos (x^{2}) , d x. \end{aligned} \end{aligned} $

Notice the subtle difference between the derivation of $d \sin^{2} (x)$ here, compared with Exercise 6, Exercise 7. In Exercise 6 we used the chain rule, while here we use Cauchy’s rule of invariance. The latter is simpler. $; \parallel$

Cauchy’s rule, like the chain rule, is a key instrument in differential calculus. Suppose we realize that $x$ in the function $\sin^{2} (x)$ depends on $t$, say $x = t^{2}$. Then, we do not need to compute the differential of $\sin^{2} (t^{2})$ all over again. Instead, we write simply

$ d \sin^{2} (t^{2}) = 2 \sin (t^{2}) \cos (t^{2}) , d t^{2} = 4 t \sin (t^{2}) \cos (t^{2}) , d t. $

Cauchy’s rule thus allows us to apply the rules of calculus sequentially, one after another.

If we know the derivative then we can compute the differential. For our purposes we need the opposite: to compute the derivative from the differential. This is possible, because

$ d f (x) = \alpha (x) , d x : ⟺ : f^{\prime} (x) = \alpha (x), \tag{14} $

where $\alpha$ may depend on $x$, but not on $d x$. Eq. (14) is a special case of the first identification theorem (Proposition 11). It shows that we can identify the derivative from the differential (and vice versa), and it shows that the concept differential is equivalent to the familiar concept derivative.

4.3. Vector functions

So far we have only concerned ourselves with scalar functions of one variable, and the reader may wonder why we bother to introduce differentials. They do not seem to have a great advantage over the more familiar derivatives. This may be true for scalar functions of one variable, but it is not true for vector functions of several variables, to which we now turn.

Thus, let $f$ now denote an $m \times 1$ vector function of an $n \times 1$ vector $x$. The derivative $D f (x)$ is an $m \times n$ matrix, but the differential $d f (x)$ has the same dimension as $f$, namely $m \times 1$. The practical advantage of working with differentials is therefore that it preserves the order of the function. This simple fact will be of great practical importance, especially when we generalize from vectors to matrices.

The differential of a vector function $f$ is defined as

$ d f (x) (u) = \left( D f (x) \right) , u, \tag{15} $

generalizing (10), where, as in the one-dimensional case, the identity function $f (x) = x$ gives $d x (u) = u$, which allows us to write

$ d f (x) = \left( D f (x) \right) , d x, \tag{16} $

generalizing (11).

Exercise 9

Find the differential of the function $f (x) = x_{1}^{2} x_{2}$.

Solution. We have

$ d f (x) = \frac{\partial f (x)}{\partial x_{1}} , d x_{1} + \frac{\partial f (x)}{\partial x_{2}} , d x_{2} = 2 x_{1} x_{2} , d x_{1} + x_{1}^{2} , d x_{2}, $

where we note that both $f (x)$ and $d f (x)$ are scalars. $; \parallel$

The operating rules (12) concerning scalar functions of one variable generalize straightforwardly to vector functions of several variables, and we have

$ \begin{aligned} d a = 0, \quad d (x^{\prime}) = (d x)^{\prime}, \quad d (a^{\prime} x) = a^{\prime} d x, \ d (x + y) = d x + d y, \quad d (x^{\prime} y) = (d x)^{\prime} y + x^{\prime} d y, \end{aligned} \tag{17} $

where $x$ and $y$ are vectors and $a$ is a vector of real constants, all of the same order.

In Exercise 9 we found the differential by first calculating the derivative, but our purpose is to find the derivative by first calculating the differential. This would work like this.

Exercise 10

Find again the differential of the function $f (x) = x_{1}^{2} x_{2}$ without first calculating the partial derivatives.

Solution. Using the operating rules, we write

$ \begin{aligned} d f (x) &= d \left( x_{1}^{2} x_{2} \right) = (d x_{1}^{2}) x_{2} + x_{1}^{2} (d x_{2}) = 2 x_{1} x_{2} , d x_{1} + x_{1}^{2} , d x_{2} \ &= (2 x_{1} x_{2}, : x_{1}^{2}) \left( \begin{array}{c} d x_{1} \ d x_{2} \end{array} \right), \end{aligned} $

from which we find the derivative as $D f (x) = (2 x_{1} x_{2}, : x_{1}^{2})$. $; \parallel$

4.4. First identification theorem

By identifying the derivative from the differential, we have generalized the identification result (14) from scalar functions of one variable to vector functions of several variables. This one-to-one relationship is formally stated as follows.

Proposition 11

$ d f (x) = A (x) , d x : ⟺ : D f (x) = A (x). $

Let us apply Proposition 11 to linear and quadratic functions, the most frequently used functions in statistics and econometrics.

Exercise 11

Find the derivative of the linear function $a^{\prime} x$, where $a$ is a vector of constants, and of the quadratic function $x^{\prime} A x$, where $A$ is a square matrix of constants.

Solution. Writing out all the steps (which one would not normally do), we have

$ \begin{aligned} d (a^{\prime} x) = (d a)^{\prime} x + a^{\prime} (d x) = a^{\prime} , d x, \ \begin{aligned} d (x^{\prime} A x) &= (d x)^{\prime} A x + x^{\prime} (d A) x + x^{\prime} A (d x) \ &= x^{\prime} A^{\prime} , d x + x^{\prime} A , d x = x^{\prime} (A + A^{\prime}) , d x, \end{aligned} \end{aligned} $

where we used the fact that the differential of a constant is 0, and also the fact that the transpose of a scalar is the same scalar, so that $(d x)^{\prime} A x = x^{\prime} A^{\prime} , d x$. Hence, the derivatives are

$ \frac{\partial a^{\prime} x}{\partial x^{\prime}} = a^{\prime}, \quad \frac{\partial x^{\prime} A x}{\partial x^{\prime}} = x^{\prime} (A + A^{\prime}), $

and in the special case where $A$ is symmetric, the derivative of $x^{\prime} A x$ is $2 x^{\prime} A$. In either case, the derivative is a row vector, not a column vector. $; \parallel$

4.5. Cauchy invariance for vector functions

Now suppose that $z = f (y)$ and that $y = g (x)$, so that $z = f (g (x))$. Then,

$ \frac{\partial z}{\partial x^{\prime}} = \frac{\partial z}{\partial y^{\prime}} , \frac{\partial y}{\partial x^{\prime}}. \tag{18} $

This is the chain rule for vector functions. The corresponding result for differentials is the following.

Proposition 12

Let $z = f (y)$ and $y = g (x)$, so that $z = f (g (x))$. Then,

$ d z = A (y) B (x) , d x, $

where $A (y)$ and $B (x)$ are defined through

$ d z = A (y) , d y, \quad d y = B (x) , d x. $

In practice, we avoid introducing a new variable $z$ or a new function $h (x) = f (g (x))$, and we write $f (y)$ when we think of $z$ as a function of $y$ and we write $f (x)$ when we think of $z$ as a function of $x$. With this convention, Proposition 12 can be reformulated as

$ d f (x) = A (y) B (x) , d x, \tag{19} $

where $d f (y) = A (y) , d y$ and $d y = d g (x) = B (x) , d x$. Such abuses of notation are very common in mathematics when one has to balance pedantic correctness and readability. In fact, the chain rule (18) is typically written as $\partial f (x) / \partial x^{\prime} = (\partial f (y) / \partial y^{\prime}) , (\partial y / \partial x^{\prime})$, without anybody complaining about it.

Exercise 12

Find the derivative of

$ f (x) = \left( \begin{array}{c} x_{1}^{2} - x_{2}^{2} \ x_{1} x_{2} x_{3} \end{array} \right), \quad x = (x_{1}, x_{2}, x_{3})^{\prime}. $

Solution. The differential is

$ \begin{aligned} d f (x) &= \left( \begin{array}{c} d (x_{1}^{2}) - d (x_{2}^{2}) \ d (x_{1} x_{2} x_{3}) \end{array} \right) = \left( \begin{array}{c} 2 x_{1} d x_{1} - 2 x_{2} , d x_{2} \ (d x_{1}) x_{2} x_{3} + x_{1} (d x_{2}) x_{3} + x_{1} x_{2} , d x_{3} \end{array} \right) \ &= \left( \begin{array}{ccc} 2 x_{1} & - 2 x_{2} & 0 \ x_{2} x_{3} & x_{1} x_{3} & x_{1} x_{2} \end{array} \right) \left( \begin{array}{c} d x_{1} \ d x_{2} \ d x_{3} \end{array} \right), \end{aligned} $

which identifies the derivative as

$ \frac{\partial f (x)}{\partial x^{\prime}} = \left( \begin{array}{ccc} 2 x_{1} & - 2 x_{2} & 0 \ x_{2} x_{3} & x_{1} x_{3} & x_{1} x_{2} \end{array} \right). \quad \parallel $

Exercise 13

Let

$ f (y) = e^{y_{1}} \sin (y_{2}), \quad y_{1} = x_{1}^{} x_{2}^{2}, \quad y_{2} = x_{1}^{2} x_{2}^{}. $

Find the derivative of $f$ with respect to $x = (x_{1}, x_{2})^{\prime}$.

Solution. Using the notational convention (19), we write

$ \begin{aligned} d f (y) &= (d e^{y_{1}}) \sin (y_{2}) + e^{y_{1}} , d \sin (y_{2}) \ &= e^{y_{1}} \sin (y_{2}) , d y_{1} + e^{y_{1}} \cos (y_{2}) , d y_{2} = a (y)^{\prime} d y, \end{aligned} $

where

$ a (y) = e^{y_{1}} \left( \begin{array}{c} \sin (y_{2}) \ \cos (y_{2}) \end{array} \right), \quad d y = \left( \begin{array}{c} d y_{1} \ d y_{2} \end{array} \right). $

Also,

$ d y = \left( \begin{array}{cc} x_{2}^{2} & 2 x_{1} x_{2} \ 2 x_{1} x_{2} & x_{1}^{2} \end{array} \right) \left( \begin{array}{c} d x_{1} \ d x_{2} \end{array} \right) = B (x) , d x. $

Hence,

$ d f (x) = a (y)^{\prime} , d y = a (y)^{\prime} B (x) , d x = c_{1} , d x_{1} + c_{2} , d x_{2}, $

where

$ \begin{aligned} c_{1} = x_{2} e^{y_{1}} \left( x_{2} \sin (y_{2}) + 2 x_{1} \cos (y_{2}) \right), \ c_{2} = x_{1} e^{y_{1}} \left( x_{1} \cos (y_{2}) + 2 x_{2} \sin (y_{2}) \right), \end{aligned} $

so that the derivative is $\partial f (x) / \partial x^{\prime} = (c_{1}, c_{2})$. $; \parallel$

5. Optimization

Section Summary: To optimize a scalar function of a vector variable, one computes its differential and sets the resulting linear expression to zero, which yields the critical point; for a positive definite quadratic this produces an explicit global minimum. When equality constraints are present, the approach is extended by forming a Lagrangian that incorporates one or more undetermined multipliers, taking the differential of the Lagrangian, and solving the resulting system of first-order conditions together with the original constraints. The procedure supplies necessary conditions for a local extremum, and these conditions identify an absolute minimum whenever the objective is strictly convex.

Let $f (x)$ be a scalar function that we wish to optimize with respect to an $n \times 1$ vector $x$. We achieve this by first computing the differential $d f (x) = a (x)^{\prime} , d x$, and then setting $a (x) = 0$.

Exercise 14

Minimize the function

$ f (x) = \frac{1}{2} x^{\prime} A x - b^{\prime} x, $

where the matrix $A$ is positive definite.

Solution. The differential is

$ d f (x) = x^{\prime} A , d x - b^{\prime} , d x = (A x - b)^{\prime} , d x, $

since a positive definite matrix is symmetric, by definition.[^1] The solution $\hat{x}$ needs to satisfy $A \hat{x} - b = 0$, and hence $\hat{x} = A^{- 1} b$. The function $f$ has an absolute minimum at $\hat{x}$, which can be seen by defining $y = x - \hat{x}$ and writing

$ y^{\prime} A y = (x - A^{- 1} b)^{\prime} A (x - A^{- 1} b) = 2 f (x) + b^{\prime} A^{- 1} b. $

Since $A$ is positive definite, $y^{\prime} A y$ has a minimum at $y = 0$ and hence $f (x)$ has a minimum at $x = \hat{x}$. Alternatively, we note that $f (x)$ is strictly convex and use the fact that any (strictly) convex function attains a (strict) absolute minimum. $; \parallel$

Next suppose there is a constraint, say $g (x) = 0$. Then we need to optimize subject to the constraint, and we need Lagrangian theory. This works as follows. First define the Lagrangian function (the Lagrangian)

$ L (x) = f (x) - \lambda g (x), \tag{20} $

where $\lambda$ is the Lagrange multiplier. Then we obtain the differential of $L (x)$ with respect to $x$,

$ d L (x) = d f (x) - \lambda d g (x), \tag{21} $

and set it equal to zero. The equations

$ \frac{\partial f (x)}{\partial x^{\prime}} = \lambda \frac{\partial g (x)}{\partial x^{\prime}}, \quad g (x) = 0 \tag{22} $

are the first-order conditions. From these $n + 1$ equations in $n + 1$ unknowns ($x$ and $\lambda$), we solve $x$ and $\lambda$.

Exercise 15

Let $A$ be positive definite. Minimize the function

$ f (x) = \frac{1}{2} x^{\prime} A x - b^{\prime} x, $

subject to the constraint $c^{\prime} x = 0$.

Solution. We write the Lagrangian as

$ L (x) = \frac{1}{2} x^{\prime} A x - b^{\prime} x - \lambda c^{\prime} x. $

The differential of $L (x)$ is

$ d L (x) = (A x - b)^{\prime} , d x - \lambda c^{\prime} , d x = (A x - b - \lambda c)^{\prime} , d x, $

from which we obtain the first-order conditions $A x - b - \lambda c = 0$ and $c^{\prime} x = 0$. This gives $x = A^{- 1} (b + \lambda c)$ and hence $0 = c^{\prime} x = c^{\prime} A^{- 1} (b + \lambda c)$. We then solve $\tilde{\lambda} = - c^{\prime} A^{- 1} b / c^{\prime} A^{- 1} c$, so that the constrained optimum is achieved at

$ \tilde{x} = A^{- 1} (b + \tilde{\lambda} c) = A^{- 1} b - \frac{c^{\prime} A^{- 1} b}{c^{\prime} A^{- 1} c} , A^{- 1} c. $

Alternatively, we can write the first-order conditions as

$ \left( \begin{array}{cc} A & c \ c^{\prime} & 0 \end{array} \right) \left( \begin{array}{c} x \ - \lambda \end{array} \right) = \left( \begin{array}{c} b \ 0 \end{array} \right) $

and solve for $x$ and $- \lambda$ by inverting the matrix on the left-hand side. Notice that by expressing the equation in $x$ and $- \lambda$ (rather than $x$ and $\lambda$) I achieve symmetry of the matrix on the left-hand side. Symmetric matrices are easier, both theoretically and computationally, than non-symmetric matrices, so, if possible, we work with symmetric matrices.

Either way, since $f (x)$ is linear-quadratic (hence strictly convex) and the constraint is linear, $f (x)$ attains an absolute minimum at $\tilde{x}$ under the constraint. $; \parallel$

If the constraint $g$ is a vector rather than a scalar, then we have not one but several (say, $m$) constraints. In that case we need $m$ multipliers and it works like this. First, define the Lagrangian

$ L (x) = f (x) - l^{\prime} g (x), \tag{23} $

where $l = (\lambda_{1}, \lambda_{2}, \ldots, \lambda_{m})^{\prime}$ is a vector of Lagrange multipliers. Then, we obtain the differential of $L (x)$ with respect to $x$:

$ d L (x) = d f (x) - l^{\prime} d g (x) \tag{24} $

and set it equal to zero. The equations

$ \frac{\partial f (x)}{\partial x^{\prime}} = l^{\prime} \frac{\partial g (x)}{\partial x^{\prime}}, \quad g (x) = 0 \tag{25} $

constitute $n + m$ equations (the first-order conditions). If we can solve these equations, then we obtain the solutions, say $\hat{x}$ and $\hat{l}$.

Exercise 16

Show that optimizing a function subject to constraints is not equivalent to optimizing the Lagrangian function.

Solution. A simple counterexample is provided by the optimization problem: $\max (x y)$ under the constraint $x + y = 2$. The solution is $x = y = 1$, but the Lagrangian $L (x, y) = x y - \lambda (x + y - 2)$ has a saddle-point at $(1, 1)$. $; \parallel$

The Lagrangian method gives necessary conditions for a local constrained extremum to occur at a given point $\hat{x}$. But how do we know that this point is in fact a maximum or a minimum? Sufficient conditions are available but they may be difficult to verify. However, in the often occurring situation where $L (x)$ is (strictly) convex, as in Exercise 14, Exercise 15, $f (x)$ attains a (strict) absolute minimum at the solution $\hat{x}$ under the constraint $g (x) = 0$. Notice that the Lagrangian $L (x, y)$ in Exercise 16 is not convex, because $\partial L (x, y) / \partial x = y - \lambda$ and $\partial L (x, y) / \partial y = x - \lambda$, so that the Hessian takes the form

$ H (x, y) = \left( \begin{array}{cc} 0 & 1 \ 1 & 0 \end{array} \right), $

which is not positive semidefinite.

[^1]: A (semi)definite matrix must be symmetric. For example, $A = \left( \begin{array}{cc} \phantom{-} 1 & 3 \ - 6 & 9 \end{array} \right)$ satisfies $x^{\prime} A x \gt 0$ for all nonzero $x$, but its eigenvalues are not real, let alone positive. This matrix is not positive definite.

6. Example 1: Least squares

Section Summary: The section explains how to find the vector x that makes the product Ax as close as possible to a given vector b, in the sense of minimizing the sum of squared differences. Differentiating the squared-error expression and setting the result to zero produces the familiar least-squares formula x-hat equals (A transpose A) inverse times A transpose b. When x must also obey linear constraints, the same approach is extended with Lagrange multipliers, yielding an adjusted solution that satisfies the constraints while still minimizing the error.

Suppose we are given an $n \times k$ matrix $A$ with linearly independent columns, so that $r (A) = k$, and an $n \times 1$ vector $b$. We wish to find a $k \times 1$ vector $x$, such that $A x$ is ‘as close as possible’ to $b$ in the sense that the ‘error’ vector $e = b - A x$ is minimized. A convenient scalar measure of the ‘error’ would be $e^{\prime} e$ and our objective is to minimize

$ f (x) = \frac{e^{\prime} e}{2} = \frac{(b - A x)^{\prime} (b - A x)}{2}, \tag{26} $

where we note that we write $e^{\prime} e / 2$ rather than $e^{\prime} e$. This makes no difference, since any $x$ which minimizes $e^{\prime} e$ will also minimize $e^{\prime} e / 2$, but it is a common trick, useful because we know that we are minimizing a quadratic function, so that a 2 will appear in the derivative. The $1 / 2$ neutralizes this $2$.

Differentiating $f (x)$ in (26) gives

$ d f (x) = e^{\prime} , d e = e^{\prime} , d (b - A x) = - e^{\prime} A , d x. $

Hence, the optimum is obtained when $A^{\prime} e = 0$, that is, when $A^{\prime} A x = A^{\prime} b$, from which we obtain

$ \hat{x} = (A^{\prime} A)^{- 1} A^{\prime} b, \tag{27} $

the least-squares solution.

If there are constraints on $x$, say $R x = r$, then we need to solve

$ \begin{aligned} \text{minimize} f (x) \ \text{subject to} R x = r. \end{aligned} \tag{28} $

We assume that the $m$ rows of $R$ are linearly independent, and define the Lagrangian

$ L (x) = (b - A x)^{\prime} (b - A x) / 2 - l^{\prime} (R x - r), $

where $l$ is a vector of Lagrange multipliers.

Exercise 17

Show that the matrix $V = R (A^{\prime} A)^{- 1} R^{\prime}$ is nonsingular if and only if the $m$ rows of $R$ are linearly independent.

Solution. Let $B = R (A^{\prime} A)^{- 1 / 2}$, so that $r (B) = r (R)$. Since $V = B B^{\prime}$, it follows that $r (V) = r (B)$, and hence that $r (V) = r (R)$. Since $V$ is an $m \times m$ matrix, it is nonsingular if and only if $r (V) = m$, that is, if and only if $r (R) = m$. $; \parallel$

We now write the differential of $L (x)$ as

$ \begin{aligned} d L (x) &= d (b - A x)^{\prime} (b - A x) / 2 - l^{\prime} , d (R x - r) \ &= (b - A x)^{\prime} , d (b - A x) - l^{\prime} R , d x \ &= - (b - A x)^{\prime} A , d x - l^{\prime} R , d x. \end{aligned} $

Setting the differential equal to zero and denoting the solutions by $\tilde{x}$ and $\tilde{l}$, we obtain the first-order conditions

$ (b - A \tilde{x})^{\prime} A + \tilde{l}^{\prime} R = 0, \quad R \tilde{x} = r, $

or, written differently,

$ A^{\prime} A \tilde{x} - A^{\prime} b = R^{\prime} \tilde{l}, \quad R \tilde{x} = r. $

We do not know $\tilde{x}$ but we know $R \tilde{x}$. Hence, we premultiply by $R (A^{\prime} A)^{- 1}$. Letting $\hat{x} = (A^{\prime} A)^{- 1} A^{\prime} b$ as in (27), this gives

$ r - R \hat{x} = R (A^{\prime} A)^{- 1} R^{\prime} \tilde{l}. $

Since we have assumed that $R$ has full row rank, we can solve for $l$:

$ \tilde{l} = \left( R (A^{\prime} A)^{- 1} R^{\prime} \right)^{- 1} (r - R \hat{x}), $

and hence for $x$:

$ \tilde{x} = \hat{x} + (A^{\prime} A)^{- 1} R^{\prime} \tilde{l} = \hat{x} + (A^{\prime} A)^{- 1} R^{\prime} \left( R (A^{\prime} A)^{- 1} R^{\prime} \right)^{- 1} (r - R \hat{x}). \tag{29} $

Since the constraint is linear and the function $f (x)$ is linear-quadratic, it follows that the solution $\tilde{x}$ indeed minimizes $f (x) = e^{\prime} e / 2$ under the constraint $R x = r$.

7. Matrix calculus

Section Summary: This section extends the rules of differentials from vectors to matrices, covering operations such as transposition, traces, sums, products, determinants, logarithms of determinants, and inverses for nonsingular cases. It shows that defining matrix derivatives through vectorization lets standard results like the chain rule carry over directly. The material also adapts constrained optimization to matrix-valued functions and constraints by using a matrix of Lagrange multipliers in the Lagrangian, yielding first-order conditions via the trace operator.

We have moved from scalar calculus to vector calculus, now we move from vector calculus to matrix calculus. The rules for vector differentials in Section 4 carry over to matrix differentials. Let $A$ be a matrix of constants and let $\alpha$ be a scalar. Then, for any $X$,

$ d A = 0, \quad d (\alpha X) = \alpha , d X, \quad d (X^{\prime}) = (d X)^{\prime}, $

and, for square $X$,

$ d , \operatorname{tr} X = \operatorname{tr} , d X. $

If $X$ and $Y$ are of the same order, then

$ d (X + Y) = d X + d Y, $

and, if the matrix product $X Y$ is defined,

$ d (X Y) = (d X) Y + X d Y. $

Two less trivial differentials are the determinant and the inverse. For nonsingular $X$ we have

$ d | X | = | X | \operatorname{tr} X^{- 1} , d X, \tag{30} $

and in particular, when $| X | \gt 0$,

$ d \log | X | = \frac{d | X |}{| X |} = \operatorname{tr} X^{- 1} , d X. \tag{31} $

The proof of (30) is a little tricky and can be found in Matrix Differential Calculus, Section 8.3, where two proofs are provided. The differential of the inverse is, for nonsingular $X$,

$ d X^{- 1} = - X^{- 1} (d X) X^{- 1}. \tag{32} $

This we can prove easily by considering the equation $X^{- 1} X = I$. Differentiating both sides gives

$ (d X^{- 1}) X + X^{- 1} d X = 0 $

and the result then follows by postmultiplying with $X^{- 1}$.

If we use the correct definition of matrix derivative (but only then), all results from vector calculus carry over to matrix calculus. In particular, the first equivalence theorem for vector functions (Proposition 11) now becomes

$ d \operatorname{vec} F (X) = A (X) d \operatorname{vec} X : ⟺ : D F (X) = A (X). \tag{33} $

The same holds for the chain rule. More precisely, if $Z = F (Y)$ and $Y = G (X)$, so that $Z = F (G (X))$, then

$ d Z = A (Y) B (X) , d X, $

where $A (Y)$ and $B (X)$ are defined through

$ d Z = A (Y) , d Y, \quad d Y = B (X) , d X, $

as in Proposition 12.

Constrained optimization, treated for vector functions in Section 5, can easily and elegantly be extended to matrix constraints. If we have a matrix $G$ (rather than a vector $g$) of constraints and a matrix $X$ (rather than a vector $x$) of variables, then we define a matrix of multipliers $L = (\lambda_{i j})$ of the same dimension as $G = (g_{i j})$. The Lagrangian then becomes

$ L (X) = f (X) - \operatorname{tr} L^{\prime} G (X), \tag{34} $

where we have used the fact, as in the proof of Proposition 1, that

$ \operatorname{tr} L^{\prime} G = \sum_{i} \sum_{j} \lambda_{i j} g_{i j}. $

If the constraint matrix $G$ is symmetric, we may take the matrix of Lagrange multipliers to be symmetric as well.

Exercise 18

Consider the optimization problem

$ \begin{aligned} \text{minimize} f (x) \ \text{subject to} G (X) = 0. \end{aligned} $

If $G (X)$ is symmetric for all $X$, then the Lagrangian function is

$ L (X) = f (X) - \operatorname{tr} L G (X), $

where $L$ may be assumed to be symmetric.

Solution. Since $G$ is symmetric, we have

$ \operatorname{tr} L^{\prime} G = \operatorname{tr} L^{\prime} G^{\prime} = \operatorname{tr} (G L)^{\prime} = \operatorname{tr} G L = \operatorname{tr} L G, $

so that $\operatorname{tr} L^{\prime} G = \operatorname{tr} L_{∗} G$, where $L_{∗} = (L + L^{\prime}) / 2$ is symmetric. $; \parallel$

From the Lagrangian (34) we obtain the differential

$ d L (X) = d f (X) - \operatorname{tr} L^{\prime} , d G (X), \tag{35} $

and, setting the differential equal to 0, we obtain the first-order conditions

$ \frac{\partial f (X)}{\partial (\operatorname{vec} X)^{\prime}} = (\operatorname{vec} L)^{\prime} : \frac{\partial \operatorname{vec} G (X)}{\partial (\operatorname{vec} X)^{\prime}}, \quad G (X) = 0, \tag{36} $

using Proposition 5 and the fact that

$ \operatorname{tr} L^{\prime} d G (X) = (\operatorname{vec} L)^{\prime} , d \operatorname{vec} G (X) = (\operatorname{vec} L)^{\prime} \frac{\partial \operatorname{vec} G (X)}{(\partial \operatorname{vec} X)^{\prime}} , d \operatorname{vec} X. $

8. Four exercises: first derivative

Section Summary: This section presents four exercises that show how to find first derivatives of scalar and matrix functions built from traces, determinants, and inverses, all by first working out the differential and then converting it to a derivative via the vec operator. The initial three exercises treat scalar cases, such as the trace of X'AX, the log-determinant of X'X, and the trace of powers of X, each yielding a compact row-vector result after simplification with cyclic properties of the trace. The last exercise derives the full matrix derivative of AX^{-1}B for nonsingular X by expressing the differential in Kronecker-product form.

Exercise 19

Obtain the derivative of the scalar function $f (X) = \operatorname{tr} X^{\prime} A X$.

Solution. The differential is

$ \begin{aligned} d f (X) &= d (\operatorname{tr} X^{\prime} A X) = \operatorname{tr} d (X^{\prime} A X) \ &= \operatorname{tr} (d X)^{\prime} A X + \operatorname{tr} X^{\prime} A , d X = \operatorname{tr} X^{\prime} (A + A^{\prime}) , d X \ &= \operatorname{tr} C^{\prime} d X = (\operatorname{vec} C)^{\prime} d \operatorname{vec} X, \end{aligned} $

using Proposition 5 and letting $C = (A + A^{\prime}) X$. Hence the derivative is $D f (X) = (\operatorname{vec} C)^{\prime}$.

Exercise 20

Obtain the derivative of the scalar function $f (X) = \log | X^{\prime} X |$, where $X$ has full column rank.

Solution. From the differential

$ \begin{aligned} d f (X) &= d \log | X^{\prime} X | = \operatorname{tr} (X^{\prime} X)^{- 1} d (X^{\prime} X) \ &= \operatorname{tr} (X^{\prime} X)^{- 1} (d X)^{\prime} X + \operatorname{tr} (X^{\prime} X)^{- 1} X^{\prime} d X = 2 \operatorname{tr} (X^{\prime} X)^{- 1} X^{\prime} d X \ &= 2 \operatorname{tr} C^{\prime} d X = 2 (\operatorname{vec} C)^{\prime} d \operatorname{vec} X, \end{aligned} $

where $C = X (X^{\prime} X)^{- 1}$, we obtain $D f (X) = 2 (\operatorname{vec} C)^{\prime}$.

Exercise 21

Let $f_{k} (X) = \operatorname{tr} X^{k}$ ($k = 1, 2, \ldots ,$). Find the derivative.

Solution. We have

$ \begin{aligned} d f_{k} (X) &= \operatorname{tr} (d X) X^{k - 1} + \operatorname{tr} X (d X) X^{k - 2} + \cdots + \operatorname{tr} X^{k - 1} d X \ &= k \operatorname{tr} X^{k - 1} d X = k (\operatorname{vec} {X^{\prime}}^{k - 1})^{\prime} d \operatorname{vec} X. \end{aligned} $

This gives $D f_{k} (X) = k (\operatorname{vec} {X^{\prime}}^{k - 1})^{\prime}$. In particular,

$ D f_{1} (X) = D \operatorname{tr} X = (\operatorname{vec} I)^{\prime}, \quad D f_{2} (X) = D \operatorname{tr} X^{2} = 2 (\operatorname{vec} X^{\prime})^{\prime}. $

Exercise 22

Find the derivative of the matrix equation $F (X) = A X^{- 1} B$, where $X$ is nonsingular.

Solution. Since

$ d F (X) = A (d X^{- 1}) B = - A X^{- 1} (d X) X^{- 1} B, $

we find

$ d \operatorname{vec} F (X) = - \left( (X^{- 1} B)^{\prime} \otimes (A X^{- 1}) \right) d \operatorname{vec} X, $

using Proposition 6, so that the derivative is given by

$ D F (X) = \frac{\partial \operatorname{vec} F (X)}{\partial (\operatorname{vec} X)^{\prime}} = - (X^{- 1} B)^{\prime} \otimes (A X^{- 1}), $

by the first identification theorem for matrices (33).

9. The second differential

Section Summary: The second differential is simply the result of differentiating the first differential once more, which for a scalar function produces a quadratic expression in the increment dx that encodes all second-order partial derivatives. This expression corresponds to the Hessian matrix through the identity d²f(x) = (dx)' H f(x) dx, with the understanding that any matrix appearing in the quadratic form must be replaced by its symmetric part to recover the Hessian uniquely. The approach yields the same second-derivative information as traditional partial differentiation but does so more directly and with less algebraic effort, as illustrated by the quadratic example f(x) = x' A x.

The second differential is simply the differential of the first differential, that is, $d^{2} f = d (d f)$. Higher-order differentials are similarly defined, but they are seldom needed.

Exercise 23

Let $f (x) = x^{\prime} A x$. Show that $d^{2} f (x) = (d x)^{\prime} (A + A^{\prime}) , d x$.

Solution. We know from Exercise 11 that $d f (x) = x^{\prime} (A + A^{\prime}) , d x$. Then,

$ \begin{aligned} d^{2} f (x) &= d \left( x^{\prime} (A + A^{\prime}) , d x \right) = (d x)^{\prime} (A + A^{\prime}) , d x + x^{\prime} (A + A^{\prime}) , d^{2} x \ &= (d x)^{\prime} (A + A^{\prime}) , d x, \end{aligned} $

since $d^{2} x = 0$. $; \parallel$

Exercise 24

Why is $d^{2} x = 0$?

Solution. If $f$ is a function of $x$, then $d x$ is short-hand for the differential $d x (u) = u$ associated with the identity function, see Section 4.1. The first derivative of the identity function is 1 (if $f$ is a scalar function of one variable), and the second derivative is 0. Hence $d^{2} x = 0$. In simpler words, if $f$ is a function of $x$ and $x$ is the ‘endpoint’, then $d^{2} x = 0$. But if $x$ is only an intermediary variable and in fact $x = x (t)$, then $d^{2} x$ is not 0 (unless $x$ is a linear function of $t$). This relationship is further developed and made more precise in Section 10. $; \parallel$

The first differential leads to the first derivative (sometimes called the Jacobian matrix) and the second differential leads to the second derivative (called the Hessian matrix). We emphasize that the concept of Hessian matrix is only useful for scalar functions, not for vector or matrix functions. When we have a vector function $f$ we shall consider the Hessian matrix of each component of $f$ separately, and when we have a matrix function $F$ we shall consider the Hessian matrix of each element of $F$ separately.

Thus, let $f$ be a scalar function and let

$ d f (x) = a (x)^{\prime} d x, \quad d a (x) = (H f (x)) , d x, \tag{37} $

where

$ a (x)^{\prime} = \frac{\partial f (x)}{\partial x^{\prime}}, \quad H f (x) = \frac{\partial a (x)}{\partial x^{\prime}} = \frac{\partial}{\partial x^{\prime}} \left( \frac{\partial f (x)}{\partial x^{\prime}} \right)^{\prime}. $

The $i j$th element of the Hessian matrix $H f (x)$ is thus obtained by first calculating $a_{j} (x) = \partial f (x) / \partial x_{j}$ and then $(H f (x)){i j} = \partial a{j} (x) / \partial x_{i}$. The Hessian matrix contains all second-order partial derivatives $\partial^{2} f (x) / \partial x_{i} , \partial x_{j}$, and it is symmetric if $f$ is twice differentiable.

The Hessian matrix is often written as

$ H f (x) = \frac{\partial^{2} f (x)}{\partial x , \partial x^{\prime}}, \tag{38} $

where the expression on the right-hand side is a notation, the precise meaning of which is given by

$ \frac{\partial^{2} f (x)}{\partial x , \partial x^{\prime}} = \frac{\partial}{\partial x^{\prime}} \left( \frac{\partial f (x)}{\partial x^{\prime}} \right)^{\prime}. \tag{39} $

Given (37) and using the symmetry of $H f (x)$, we obtain the second differential as

$ d^{2} f (x) = (d a)^{\prime} , d x = (d x)^{\prime} (H f (x)) , d x, \tag{40} $

which shows that the second differential of $f$ is a quadratic form in $d x$.

Now, suppose that we have obtained $d^{2} f (x) = (d x)^{\prime} B (x) , d x$. We also know that $d^{2} f (x) = (d x)^{\prime} H f (x) , d x$ by the definition of the Hessian. Hence,

$ (d x)^{\prime} (H f (x) - B (x)) , d x = 0 $

for all $d x$, where $H f (x)$ is symmetric, but $B (x)$ is not necessarily symmetric. Does this imply that $H f (x) = B (x)$? No, it does not, as we have seen in Section 2.2. It does, however, imply that

$ (H f (x) - B (x))^{\prime} + (H f (x) - B (x)) = 0, $

and hence that $H f (x) = (B (x) + B (x)^{\prime}) / 2$, using the symmetry of $H f (x)$. This proves the following result.

Proposition 13

$ d^{2} f (x) = (d x)^{\prime} B (x) , d x : ⟺ : H f (x) = \frac{B (x) + B (x)^{\prime}}{2}. $

The second identification theorem shows that there is a one-to-one correspondence between second-order differentials and second-order derivatives, but only if we make the matrix $B (x)$ in the quadratic form symmetric. Hence, the second differential identifies the second derivative.

Exercise 25

Consider again the quadratic function $f (x) = x^{\prime} A x$. Find the second differential without writing out the first differential in its final form.

Solution. We can start with $d f (x) = x^{\prime} (A + A^{\prime}) , d x$, as in Exercise 23, and obtain $d^{2} f (x) = (d x)^{\prime} (A + A^{\prime}) , d x$. The matrix in the quadratic form is already symmetric, so we obtain directly $H f (x) = A + A^{\prime}$. Alternatively — and this is often quicker — we differentiate $f$ twice without writing out the first differential in its final form, as follows. From

$ d f (x) = (d x)^{\prime} A x + x^{\prime} A , d x, $

we obtain

$ d^{2} f (x) = 2 (d x)^{\prime} A , d x, $

which immediately identifies the Hessian matrix as $H f (x) = A + A^{\prime}$. $; \parallel$

Even with such a simple function as $f (x) = x^{\prime} A x$, the advantage and elegance of using differentials is clear. Without differentials we would need to prove first that $\partial a^{\prime} x / \partial x^{\prime} = a^{\prime}$ and $\partial x^{\prime} A x / \partial x^{\prime} = x^{\prime} (A + A^{\prime})$, and then use (39) to obtain

$ \frac{\partial^{2} x^{\prime} A x}{\partial x , \partial x^{\prime}} = \frac{\partial (x^{\prime} (A + A^{\prime}))^{\prime}}{\partial x^{\prime}} = \frac{\partial (A + A^{\prime}) x}{\partial x^{\prime}} = A + A^{\prime}, $

which is cumbersome in this simple case and not practical in more complex situations.

10. Chain rule for second differentials

Section Summary: The section explains why the second differential of a variable may be zero or nonzero depending on whether that variable is treated as independent or as a function of something else. It derives the general chain rule for the second differential of a composite function z = f(g(x)), expressing d²z first in terms of the intermediate variable y and then fully in terms of dx. Two worked examples then illustrate computing the Hessian matrix of a composition, once using the expanded chain-rule formula and once by retaining the intermediate differentials throughout.

Let us now further analyze the question why sometimes $d^{2} x = 0$ and sometimes $d^{2} x \neq 0$, first considered in Exercise 24. If $f$ is a function of $x$, and $x$ is the argument of interest, then $d^{2} x = 0$. But if $f$ is a function of $x$, which in turn is a function of $t$, then it is no longer true that $d^{2} x$ equals zero.

More generally, suppose that $z = f (y)$ and that $y = g (x)$, so that $z = f (g (x))$. Then,

$ d z = A (y) , d y $

and

$ d^{2} z = (d A) , d y + A (y) , d^{2} y. \tag{41} $

This is true whether or not $y$ depends on some other variables. If we think of $z$ as a function of $y$, then $d^{2} y = 0$, but if $y$ depends on $x$ then $d^{2} y$ is not zero; in fact,

$ d y = B (x) , d x, \quad d^{2} y = (d B) , d x. $

This leads to the following result.

Proposition 14

Let $z = f (y)$ and $y = g (x)$, so that $z = f (g (x))$. Then,

$ d^{2} z = (d A) B (x) , d x + A (y) (d B) , d x, $

where $A (y)$ and $B (x)$ are defined through

$ d z = A (y) , d y, \quad d y = B (x) , d x. $

Personally, I usually avoid Proposition 14. Instead, I express $d^{2} z$ in terms of $d y$ and $d^{2} y$, as in (41), and proceed from there. Let me give two examples, one using Proposition 14, the other not using the proposition.

Exercise 26

Let

$ f (y_{1}, y_{2}) = e^{y_{1}} \sin (y_{2}), \quad y_{1} = x_{1}^{} x_{2}^{2}, \quad y_{2} = x_{1}^{2} x_{2}^{}. $

Find the Hessian matrix by using Proposition 14.

Solution. By Proposition 14,

$ d^{2} f (x) = (d a)^{\prime} B (x) , d x + a (y)^{\prime} (d B) , d x, $

where

$ a (y) = e^{y_{1}} \left( \begin{array}{c} \sin (y_{2}) \ \cos (y_{2}) \end{array} \right), \quad B (x) = \left( \begin{array}{cc} x_{2}^{2} & 2 x_{1} x_{2} \ 2 x_{1} x_{2} & x_{1}^{2} \end{array} \right). $

Now, letting

$ C (y) = e^{y_{1}} \left( \begin{array}{cc} \sin (y_{2}) & \phantom{-} \cos (y_{2}) \ \cos (y_{2}) & - \sin (y_{2}) \end{array} \right) $

and

$ D_{1} (x) = 2 \left( \begin{array}{cc} 0 & x_{2} \ x_{2} & x_{1} \end{array} \right), \quad D_{2} (x) = 2 \left( \begin{array}{cc} x_{2} & x_{1} \ x_{1} & 0 \end{array} \right), $

we obtain

$ d a = C (y) , d y = C (y) B (x) , d x $

and

$ d B = (d x_{1}) D_{1} (x) + (d x_{2}) D_{2} (x). $

Let us write $d x_{1}$ and $d x_{2}$ in terms of $d x$, which can be done by defining $e_{1} = (1, 0)^{\prime}$ and $e_{2} = (0, 1)^{\prime}$. Then, $d x_{1} = e_{1}^{\prime} d x$ and $d x_{2} = e_{2}^{\prime} d x$, and hence

$ \begin{aligned} d^{2} f (x) &= (d a)^{\prime} B (x) , d x + a (y)^{\prime} (d B) , d x \ &= (d x)^{\prime} B (x) C (y) B (x) , d x + a (y)^{\prime} \left( (d x_{1}) D_{1} (x) + (d x_{2}) D_{2} (x) \right) , d x \ &= (d x)^{\prime} B (x) C (y) B (x) , d x + (d x)^{\prime} e_{1} a (y)^{\prime} D_{1} (x) , d x + (d x)^{\prime} e_{2} a (y)^{\prime} D_{2} (x) , d x \ &= (d x)^{\prime} \left[ B (x) C (y) B (x) + e_{1} a (y)^{\prime} D_{1} (x) + e_{2} a (y)^{\prime} D_{2} (x) \right] , d x. \end{aligned} $

Some care is required where to position the scalars $e_{1}^{\prime} d x$ and $e_{2}^{\prime} d x$ in the matrix product. A scalar can be positioned anywhere in a matrix product, but we wish to position the two scalars in such a way that the usual matrix multiplication rules still apply.

Having obtained the second differential in the desired form, Proposition 13 implies that the Hessian is equal to

$ \begin{aligned} H f (x) &= B (x) C (y) B (x) + \frac{1}{2} \left( e_{1} a (y)^{\prime} D_{1} (x) + D_{1} (x) a (y) e_{1}^{\prime} \right) \ &\quad + \frac{1}{2} \left( e_{2} a (y)^{\prime} D_{2} (x) + D_{2} (x) a (y) e_{2}^{\prime} \right). \quad \parallel \end{aligned} $

Exercise 27

Now find the Hessian matrix, not using Proposition 14.

Solution. We have

$ d f (y) = (d e^{y_{1}}) \sin (y_{2}) + e^{y_{1}} d \sin (y_{2}) = e^{y_{1}} \sin (y_{2}) , d y_{1} + e^{y_{1}} \cos (y_{2}) , d y_{2}, $

and hence

$ \begin{aligned} d^{2} f (y) &= [ d e^{y_{1}} ] \sin (y_{2}) , d y_{1} + e^{y_{1}} [ d \sin (y_{2}) ] , d y_{1} + e^{y_{1}} \sin (y_{2}) , d^{2} y_{1} \ &\quad + [ d e^{y_{1}} ] \cos (y_{2}) , d y_{2} + e^{y_{1}} [ d \cos (y_{2}) ] , d y_{2} + e^{y_{1}} \cos (y_{2}) , d^{2} y_{2} \ &= e^{y_{1}} \sin (y_{2}) , (d y_{1})^{2} + e^{y_{1}} \cos (y_{2}) , d y_{1} , d y_{2} + e^{y_{1}} \sin (y_{2}) , d^{2} y_{1} \ &\quad + e^{y_{1}} \cos (y_{2}) , d y_{1} , d y_{2} - e^{y_{1}} \sin (y_{2}) , (d y_{2})^{2} + e^{y_{1}} \cos (y_{2}) , d^{2} y_{2}, \end{aligned} $

where we emphasize that $d^{2} y_{1}$ and $d^{2} y_{2}$ are not 0, because they depend on $x_{1}$ and $x_{2}$. In fact,

$ d y_{1} = x_{2}^{2} , d x_{1} + 2 x_{1} x_{2} , d x_{2}, \quad d^{2} y_{1} = 4 x_{2} , d x_{1} , d x_{2} + 2 x_{1} (d x_{2})^{2}, $

and

$ d y_{2} = 2 x_{1} x_{2} , d x_{1} + x_{1}^{2} , d x_{2}, \quad d^{2} y_{2} = 2 x_{2} (d x_{1})^{2} + 4 x_{1} , d x_{1} , d x_{2}. $

Inserting these expressions into $d^{2} f (y)$ gives the required result. $; \parallel$

11. The Hessian matrix

Section Summary: In moving from vector to matrix calculus, the second differential of a scalar function is converted into its Hessian matrix by symmetrizing the coefficient matrix that multiplies the vectorized differentials. Trace expressions commonly arising in this setting are rewritten using Kronecker products and commutation matrices, which automatically enforce the correct ordering of variables and yield explicit formulas for the Hessian. The same technique extends without change to symmetric-matrix cases by inserting duplication matrices that respect the symmetry constraints.

When we move from vector calculus to matrix calculus, we need an ordering of the functions and of the variables. As motivated in Section 3, we shall view the matrix function $F (X)$ as a vector function $f (x)$, where $f = \operatorname{vec} F$ and $x = \operatorname{vec} X$. We already obtained in (33) the extension of the first identification theorem,

$ d \operatorname{vec} F (X) = A (X) , d \operatorname{vec} X : ⟺ : \frac{\partial \operatorname{vec} F (X)}{\partial (\operatorname{vec} X)^{\prime}} = A (X). $

For the second identification theorem, we obtain

$ d^{2} f (X) = (d \operatorname{vec} X)^{\prime} B (X) , d \operatorname{vec} X : ⟺ : H f (X) = \frac{B (X) + B (X)^{\prime}}{2}. \tag{42} $

Notice that we only provide the Hessian matrix for scalar functions, not for vector or matrix functions, as explained in Section 9.

The commutation matrix, introduced in Section 2.5, has many applications in matrix theory, and it is essential in identifying the Hessian matrix from the second differential. The second differential of a scalar function often takes the form of a trace, either $\operatorname{tr} A (d X)^{\prime} B d X$ or $\operatorname{tr} A (d X) B d X$. The following result is then of importance.

Proposition 15

Let $f$ be a twice differentiable real-valued function of an $n \times q$ matrix $X$. Then,

$ d^{2} f (X) = \operatorname{tr} A (d X)^{\prime} B d X : ⟺ : H f (X) = \frac{1}{2} (A^{\prime} \otimes B + A \otimes B^{\prime}) $

and

$ d^{2} f (X) = \operatorname{tr} A (d X) B d X : ⟺ : H f (X) = \frac{1}{2} K_{q n} (A^{\prime} \otimes B + B^{\prime} \otimes A). $

Proof

We write

$ \begin{aligned} \operatorname{tr} A (d X)^{\prime} B d X &= \operatorname{tr} (d X)^{\prime} B (d X) A = (\operatorname{vec} d X)^{\prime} \operatorname{vec} B (d X) A \ &= (\operatorname{vec} d X)^{\prime} (A^{\prime} \otimes B) \operatorname{vec} d X = (d \operatorname{vec} X)^{\prime} (A^{\prime} \otimes B) d \operatorname{vec} X \end{aligned} $

and

$ \begin{aligned} \operatorname{tr} A (d X) B d X &= \operatorname{tr} (d X)^{\prime} B^{\prime} (d X)^{\prime} A^{\prime} = (\operatorname{vec} d X)^{\prime} \operatorname{vec} B^{\prime} (d X)^{\prime} A^{\prime} \ &= (\operatorname{vec} d X)^{\prime} (A \otimes B^{\prime}) \operatorname{vec} d X^{\prime} = (d \operatorname{vec} X)^{\prime} (A \otimes B^{\prime}) K_{n q} d \operatorname{vec} X, \end{aligned} $

using Proposition 5, Proposition 6, and properties of the commutation matrix. The result now follows from Proposition 13, using the fact that

$ (A \otimes B^{\prime}) K_{n q} + K_{q n} (A^{\prime} \otimes B) = K_{q n} (A^{\prime} \otimes B + B^{\prime} \otimes A). \quad \square $

Exercise 28

Suppose that $d^{2} f (X) = 2 \operatorname{tr} A (d X) A , d X$ and that both $A$ and $X$ are known to be symmetric $n \times n$ matrices. Show that $H f (X) = 2 D_{n}^{\prime} (A \otimes A) D_{n}^{}$.

Solution. This follows from

$ \begin{aligned} \operatorname{tr} A (d X)^{\prime} A , d X &= (d \operatorname{vec} X)^{\prime} (A \otimes A) d \operatorname{vec} X \ &= (d \operatorname{vech} (X))^{\prime} D_{n}^{\prime} (A \otimes A) D_{n}^{} , d \operatorname{vech} (X). \quad \parallel \end{aligned} $

12. Four exercises: second derivative

Section Summary: This section derives the Hessians, or second derivatives, for four matrix functions previously studied through first differentials. For each case it starts from the known first differential, computes the second differential using rules for traces, inverses, and products, and converts the result into an explicit Hessian expression that involves Kronecker products, commutation matrices, and projection terms. The functions treated are the quadratic trace tr(X′AX), the log-determinant log|X′X|, the power traces tr(X^k), and the linear fractional map AX^{-1}B.

In this section we consider the same functions as in Section 8, but now we obtain the second derivatives.

Exercise 29

Obtain the Hessian of $f (X) = \operatorname{tr} X^{\prime} A X$.

Solution. We know that $d f (X) = \operatorname{tr} C^{\prime} d X$, where $C = (A + A^{\prime}) X$. The second differential is

$ d^{2} f (X) = d \operatorname{tr} X^{\prime} (A + A^{\prime}) , d X = \operatorname{tr} (d X)^{\prime} (A + A^{\prime}) , d X. $

Hence, the Hessian is $H f (X) = I_{q} \otimes (A + A^{\prime})$.

Exercise 30

Obtain the Hessian of $f (X) = \log | X^{\prime} X |$, where $X$ has full column rank.

Solution. Since $d f (X) = 2 \operatorname{tr} C^{\prime} d X$, where $C = X (X^{\prime} X)^{- 1}$, the second differential is

$ \begin{aligned} d^{2} f (X) &= 2 , d \left( \operatorname{tr} (X^{\prime} X)^{- 1} X^{\prime} d X \right) \ &= 2 \operatorname{tr} (d (X^{\prime} X)^{- 1}) X^{\prime} d X + 2 \operatorname{tr} (X^{\prime} X)^{- 1} (d X)^{\prime} d X \ &= - 2 \operatorname{tr} (X^{\prime} X)^{- 1} (d X^{\prime} X) (X^{\prime} X)^{- 1} X^{\prime} d X + 2 \operatorname{tr} (X^{\prime} X)^{- 1} (d X)^{\prime} d X \ &= - 2 \operatorname{tr} (X^{\prime} X)^{- 1} (d X)^{\prime} X (X^{\prime} X)^{- 1} X^{\prime} d X \ &\quad - 2 \operatorname{tr} (X^{\prime} X)^{- 1} X^{\prime} (d X) (X^{\prime} X)^{- 1} X^{\prime} d X + 2 \operatorname{tr} (X^{\prime} X)^{- 1} (d X)^{\prime} d X \ &= 2 \operatorname{tr} (X^{\prime} X)^{- 1} (d X)^{\prime} M d X - 2 \operatorname{tr} (X^{\prime} X)^{- 1} X^{\prime} (d X) (X^{\prime} X)^{- 1} X^{\prime} d X \ &= 2 \operatorname{tr} (X^{\prime} X)^{- 1} (d X)^{\prime} M d X - 2 \operatorname{tr} C^{\prime} (d X) C^{\prime} d X, \end{aligned} $

where $M = I_{n} - X (X^{\prime} X)^{- 1} X^{\prime}$. The second equality in this derivation follows from considering $(X^{\prime} X)^{- 1} X^{\prime} d X$ as a product of three matrices: $(X^{\prime} X)^{- 1}$, $X^{\prime}$, and $d X$ (a matrix of constants), the third equality uses the differential of the inverse in (32), and the fourth equality separates $d X^{\prime} X$ into $(d X)^{\prime} X + X^{\prime} d X$. Hence, we obtain from Proposition 15,

$ H f (X) = 2 (X^{\prime} X)^{- 1} \otimes M - 2 K_{q n} (C \otimes C^{\prime}). $

Exercise 31

Let $f_{k} (X) = \operatorname{tr} X^{k}$ ($k = 1, 2, \ldots ,$). Find the Hessian.

Solution. We have $d f_{k} (X) = k \operatorname{tr} X^{k - 1} d X$, and, in particular, $d f_{1} (X) = \operatorname{tr} d X$ and $d f_{2} (X) = 2 \operatorname{tr} X , d X$. Hence, $d^{2} f_{1} (X) = 0$ and, for $k \geq 2$,

$ d^{2} f_{k} (X) = k , \operatorname{tr} , (d X^{k - 1}) , d X = k , \sum_{j = 0}^{k - 2} \operatorname{tr} X^{j} (d X) X^{k - 2 - j} d X. $

This gives $H f_{1} (X) = 0$, and, for $k \geq 2$,

$ H f_{k} (X) = (k / 2) \sum_{j = 0}^{k - 2} K_{n} ({X^{\prime}}^{j} \otimes X^{k - 2 - j} + {X^{\prime}}^{k - 2 - j} \otimes X^{j}). $

Exercise 32

Find the Hessian of the matrix equation $F (X) = A X^{- 1} B$, where $X$ is nonsingular.

Solution. We know that $d F (X) = - A X^{- 1} (d X) X^{- 1} B$. The second differential is

$ \begin{aligned} d^{2} F (X) &= - A (d X^{- 1}) (d X) X^{- 1} B - A X^{- 1} (d X) (d X^{- 1}) B \ &= 2 A X^{- 1} (d X) X^{- 1} (d X) X^{- 1} B. \end{aligned} $

To obtain the Hessian matrix of the $s t$th element of $F$, we let

$ C_{t s} = X^{- 1} B e_{t}^{} e_{s}^{\prime} A X^{- 1}, $

where $e_{s}$ and $e_{t}$ are elementary vectors with 1 in the $s$th (respectively, $t$th) position and zeros elsewhere. Then,

$ d^{2} F_{s t} (X) = 2 e_{s}^{\prime} A X^{- 1} (d X) X^{- 1} (d X) X^{- 1} B e_{t}^{} = 2 \operatorname{tr} C_{t s} (d X) X^{- 1} (d X), $

and hence

$ H F_{s t} (X) = \frac{\partial^{2} F_{s t}}{(\partial \operatorname{vec} X) (\partial \operatorname{vec} X)^{\prime}} = K_{n}^{} (C_{t s}^{\prime} \otimes X^{- 1} + {X^{\prime}}^{- 1} \otimes C_{t s}). $

13. Example 2: Maximum likelihood

Section Summary: This section derives the maximum likelihood estimators for the mean vector and covariance matrix of a multivariate normal distribution from a sample of independent observations. It begins with the log-likelihood function expressed in terms of the parameters, differentiates it using matrix calculus to locate the maximum, and obtains the familiar sample mean and (biased) sample covariance formulas. The section then computes the observed information matrix from the second differential, employing the duplication matrix to respect symmetry of the covariance estimator, and reports its explicit inverse and determinant.

Consider a sample of $m \times 1$ vectors $y_{1}, y_{2}, \ldots, y_{n}$ from the multivariate normal distribution with mean $\mu$ and variance $\Omega$, where $\Omega$ is positive definite and $n \geq m + 1$. The density of $y_{i}$ is

$ f (y_{i}) = (2 \pi)^{- m / 2} | \Omega |^{- 1 / 2} \exp \left( - \frac{1}{2} (y_{i} - \mu)^{\prime} \Omega^{- 1} (y_{i} - \mu) \right), $

and since the $y_{i}$ are independent and identically distributed, the joint density of $(y_{1}, \ldots, y_{n})$ is given by $\prod_{i} f (y_{i})$. The likelihood is equal to the joint density, but now thought of as a function of the parameters $\mu$ and $\Omega$, rather than of the observations. Its logarithm is the loglikelihood, which here takes the form

$ \Lambda (\mu, \Omega) = - \frac{m n}{2} \log 2 \pi - \frac{n}{2} \log | \Omega | - \frac{1}{2} \sum_{i = 1}^{n} (y_{i} - \mu)^{\prime} \Omega^{- 1} (y_{i} - \mu). \tag{43} $

The maximum likelihood (ML) estimators are obtained by maximizing the loglikelihood (which is the same, but usually easier, as maximizing the likelihood). Thus, we differentiate $\Lambda$ and obtain

$ \begin{aligned} d \Lambda &= - \frac{n}{2} , d \log | \Omega | + \frac{1}{2} \sum_{i = 1}^{n} (d \mu)^{\prime} \Omega^{- 1} (y_{i} - \mu) \ &\quad - \frac{1}{2} \sum_{i = 1}^{n} (y_{i} - \mu)^{\prime} (d \Omega^{- 1}) (y_{i} - \mu) + \frac{1}{2} \sum_{i = 1}^{n} (y_{i} - \mu)^{\prime} \Omega^{- 1} d \mu \ &= - \frac{n}{2} , d \log | \Omega | - \frac{1}{2} \sum_{i = 1}^{n} (y_{i} - \mu)^{\prime} (d \Omega^{- 1}) (y_{i} - \mu) + \sum_{i = 1}^{n} (y_{i} - \mu)^{\prime} \Omega^{- 1} d \mu \ &= - \frac{n}{2} \operatorname{tr} (\Omega^{- 1} d \Omega + S d \Omega^{- 1}) + \sum_{i = 1}^{n} (y_{i} - \mu)^{\prime} \Omega^{- 1} d \mu, \end{aligned} \tag{44} $

where

$ S = S (\mu) = \frac{1}{n} \sum_{i = 1}^{n} (y_{i} - \mu) (y_{i} - \mu)^{\prime}. $

Denoting the ML estimators by $\hat{\mu}$ and $\hat{\Omega}$, letting $\hat{S} = S (\hat{\mu})$, and setting $d \Lambda = 0$ then implies that

$ \operatorname{tr} \left( \hat{\Omega}^{- 1} - \hat{\Omega}^{- 1} \hat{S} \hat{\Omega}^{- 1} \right) d \Omega = 0, \quad \sum_{i = 1}^{n} (y_{i} - \hat{\mu})^{\prime} \hat{\Omega}^{- 1} d \mu = 0, $

for all $d \Omega$ and all $d \mu$. This, in turn, implies that

$ \hat{\Omega}^{- 1} = \hat{\Omega}^{- 1} \hat{S} \hat{\Omega}^{- 1}, \quad \sum_{i = 1}^{n} (y_{i} - \hat{\mu}) = 0. $

Hence, the ML estimators are given by

$ \hat{\mu} = \frac{1}{n} \sum_{i = 1}^{n} y_{i} = \bar{y}, \quad \hat{\Omega} = \frac{1}{n} \sum_{i = 1}^{n} (y_{i} - \bar{y}) (y_{i} - \bar{y})^{\prime}. \tag{45} $

We note that the condition that $\Omega$ is symmetric has not been imposed. But since the solution $\hat{\Omega}$ is symmetric, imposing the condition would have made no difference.

The second differential is obtained by differentiating (44) again. This gives

$ \begin{aligned} d^{2} \Lambda &= - \frac{n}{2} \operatorname{tr} \left( (d \Omega^{- 1}) d \Omega + (d S) d \Omega^{- 1} + S d^{2} \Omega^{- 1} \right) - n (d \mu)^{\prime} \Omega^{- 1} d \mu \ &\quad + \sum_{i = 1}^{n} (y_{i} - \mu)^{\prime} (d \Omega^{- 1}) d \mu. \end{aligned} \tag{46} $

From the second differential (46) we can obtain the Hessian matrix by differentiating further, but if one is only interested in the information matrix (minus the expectation of the Hessian matrix), then an important shortcut is possible by taking expectations at this stage. Since $E (S) = \Omega$ and $E (d S) = 0$, we obtain

$ \begin{aligned} E d^{2} \Lambda &= - \frac{n}{2} \operatorname{tr} \left( (d \Omega^{- 1}) d \Omega + \Omega d^{2} \Omega^{- 1} \right) - n (d \mu)^{\prime} \Omega^{- 1} d \mu \ &= \frac{n}{2} \operatorname{tr} \Omega^{- 1} (d \Omega) \Omega^{- 1} d \Omega - n \operatorname{tr} (d \Omega) \Omega^{- 1} (d \Omega) \Omega^{- 1} - n (d \mu)^{\prime} \Omega^{- 1} d \mu \ &= - \frac{n}{2} \operatorname{tr} \Omega^{- 1} (d \Omega) \Omega^{- 1} d \Omega - n (d \mu)^{\prime} \Omega^{- 1} d \mu, \end{aligned} \tag{47} $

using the facts that $d \Omega^{- 1} = - \Omega^{- 1} (d \Omega) \Omega^{- 1}$ and

$ \begin{aligned} d^{2} \Omega^{- 1} &= - (d \Omega^{- 1}) (d \Omega) \Omega^{- 1} - \Omega^{- 1} (d \Omega) d \Omega^{- 1} \ &= 2 \Omega^{- 1} (d \Omega) \Omega^{- 1} (d \Omega) \Omega^{- 1}. \end{aligned} $

To obtain the information matrix we need to take the symmetry of $\Omega$ into account and this is where the duplication matrix appears. So far, we have avoided the vec operator and in practical situations one should work with differentials (rather than with derivatives) as long as possible. But we cannot go further than (47) without use of the vec operator. Thus, from (47),

$ \begin{aligned}

• E d^{2} \Lambda &= \frac{n}{2} \operatorname{tr} \Omega^{- 1} (d \Omega) \Omega^{- 1} d \Omega + n (d \mu)^{\prime} \Omega^{- 1} d \mu \ &= \frac{n}{2} (d \operatorname{vec} \Omega)^{\prime} (\Omega^{- 1} \otimes \Omega^{- 1}) , d \operatorname{vec} \Omega + n (d \mu)^{\prime} \Omega^{- 1} d \mu \ &= \frac{n}{2} (d \operatorname{vech} (\Omega))^{\prime} D_{m}^{\prime} (\Omega^{- 1} \otimes \Omega^{- 1}) D_{m} , d \operatorname{vech} (\Omega) + n (d \mu)^{\prime} \Omega^{- 1} d \mu, \end{aligned} $

which implies that the information matrix for $\mu$ and $\operatorname{vech} (\Omega)$ takes the form

$ F = n \left( \begin{array}{cc} \Omega^{- 1} & 0 \ 0 & \frac{1}{2} D_{m}^{\prime} (\Omega^{- 1} \otimes \Omega^{- 1}) D_{m} \end{array} \right). \tag{48} $

The results on the duplication matrix in Section 2.6 allow us to obtain the inverse:

$ (F / n)^{- 1} = \left( \begin{array}{cc} \Omega & 0 \ 0 & 2 (D_{m}^{\prime} D_{m})^{- 1} D_{m}^{\prime} (\Omega \otimes \Omega) D_{m} (D_{m}^{\prime} D_{m})^{- 1} \end{array} \right) $

and the determinant:

$ | F / n | = | \Omega | \cdot | 2 (D_{m}^{\prime} D_{m}^{})^{- 1} D_{m}^{\prime} (\Omega \otimes \Omega) D_{m}^{} (D_{m}^{\prime} D_{m}^{})^{- 1} | = 2^{m} | \Omega |^{m + 2}. $

14. Example 3: Maximum likelihood with parameters in the design matrix

Section Summary: This section extends the standard linear regression model by allowing the design matrix X itself to depend on unknown parameters ψ, in addition to the usual regression coefficients β and variance parameters θ in Ω. Under normality the authors derive the log-likelihood, concentrate it analytically with respect to β to obtain a simple closed-form estimator, and then maximize the reduced likelihood over ψ and θ (via grid search if needed). They also compute the information matrix and show that it is block-diagonal, so uncertainty about the variance parameters does not affect the precision of the estimates of β and ψ.

Let us consider the linear model $y = X \beta + u$, where $E (u) = 0$ and $\operatorname{var} (u) = \Omega (\theta)$. In contrast to the standard linear model we assume, in addition, that $X$ depends on a vector of parameters $ψ$. We wish to find the ML estimators of $\beta$, $ψ$, and $\theta$, and the information matrix, whose inverse approximates the variance of the ML estimators. This model was considered in Ikefuji et al. (2022).

Under normality, the loglikelihood takes the form

$ \Lambda (\beta, ψ, \theta) = \text{constant} - (1 / 2) \log | \Omega | - (1 / 2) (y - X \beta)^{\prime} \Omega^{- 1} (y - X \beta). \tag{49} $

Maximizing $\Lambda$ with respect to $\beta$ and $\theta$ is (relatively) easy, while maximization with respect to $ψ$ is more difficult. Upon differentiating $X \beta$ we obtain

$ d (X \beta) = X d \beta + (d X) \beta = X d \beta + (\beta^{\prime} \otimes I_{n}) Z d ψ, $

where $Z = \partial \operatorname{vec} X / \partial ψ^{\prime}$. Differentiating the loglikelihood then gives

$ \begin{aligned} d \Lambda &= - (1 / 2) \operatorname{tr} (\Omega^{- 1} d \Omega) + (1 / 2) (y - X \beta)^{\prime} \Omega^{- 1} (d \Omega) \Omega^{- 1} (y - X \beta) \ &\quad + (y - X \beta)^{\prime} \Omega^{- 1} X d \beta + (y - X \beta)^{\prime} \Omega^{- 1} (d X) \beta, \end{aligned} \tag{50} $

from which we obtain the first-order conditions

$ \begin{aligned} (y - X \beta)^{\prime} \Omega^{- 1} X d \beta = 0, \ (y - X \beta)^{\prime} \Omega^{- 1} (d \Omega) \Omega^{- 1} (y - X \beta) = \operatorname{tr} (\Omega^{- 1} d \Omega), \ (y - X \beta)^{\prime} \Omega^{- 1} (d X) \beta = 0, \end{aligned} $

for $\beta$, $\theta$, and $ψ$, respectively. This implies that $\hat{\beta}$ takes the simple form

$ \hat{\beta} (ψ, \theta) = (X^{\prime} \Omega^{- 1} X)^{- 1} X^{\prime} \Omega^{- 1} y, \tag{51} $

so we can concentrate the likelihood with respect to $\beta$. The concentrated loglikelihood is

$ \Lambda_{c} (ψ, \theta) = \text{constant} - (1 / 2) \log | \Omega | - (1 / 2) \hat{u}^{\prime} \Omega^{- 1} \hat{u}, $

where $\hat{u} = y - X \hat{\beta} = y - X (X^{\prime} \Omega^{- 1} X)^{- 1} X^{\prime} \Omega^{- 1} y$. We obtain $\hat{ψ}$ and $\hat{\theta}$ by maximizing $\Lambda_{c}$, and then $\hat{\beta}$ through (51).

One way to maximize $\Lambda_{c}$ is to perform a grid search over $ψ$. Fixing $ψ = ψ_{0}$, we have $X = X (ψ_{0})$ and $d ψ = 0$. Then,

$ \begin{aligned} d \hat{\beta} &= [ d (X^{\prime} \Omega^{- 1} X)^{- 1} ] X^{\prime} \Omega^{- 1} y + (X^{\prime} \Omega^{- 1} X)^{- 1} d (X^{\prime} \Omega^{- 1} y) \ &= - (X^{\prime} \Omega^{- 1} X)^{- 1} X^{\prime} \Omega^{- 1} (d \Omega) \Omega^{- 1} \hat{u}, \end{aligned} $

and hence

$ \begin{aligned} d \Lambda_{c} (ψ_{0}, \theta) &= - (1 / 2) \operatorname{tr} (\Omega^{- 1} d \Omega) + (1 / 2) \hat{u}^{\prime} \Omega^{- 1} (d \Omega) \Omega^{- 1} \hat{u} \ &\quad - \hat{u}^{\prime} \Omega^{- 1} X (X^{\prime} \Omega^{- 1} X)^{- 1} X^{\prime} \Omega^{- 1} (d \Omega) \Omega^{- 1} \hat{u}. \end{aligned} $

Solving $\theta$ from $d \Lambda_{c} (ψ_{0}, \theta) = 0$ gives $\hat{\theta} (ψ_{0})$ and hence $\Lambda_{c} (ψ_{0}, \hat{\theta} (ψ_{0}))$. Performing a grid search over $ψ$ we find $\hat{ψ}$ where $\Lambda_{c} (ψ, \hat{\theta} (ψ))$ is maximized. Given $\hat{ψ}$ and $\hat{\theta}$ we then compute $\hat{\beta}$ from (51).

To obtain the information matrix we take the differential of $d \Lambda$ in (50). This gives, letting $\mu = X \beta$,

$ \begin{aligned} d^{2} \Lambda &= (1 / 2) \operatorname{tr} (\Omega^{- 1} d \Omega)^{2} - (y - X \beta)^{\prime} \Omega^{- 1} (d \Omega) \Omega^{- 1} (d \Omega) \Omega^{- 1} (y - X \beta) \ &\quad - (d \mu)^{\prime} \Omega^{- 1} (d \mu) - 2 (y - X \beta)^{\prime} \Omega^{- 1} (d \Omega) \Omega^{- 1} (d \mu) + (y - X \beta)^{\prime} \Omega^{- 1} (d^{2} \mu) \ &\quad - (1 / 2) \operatorname{tr} (\Omega^{- 1} d^{2} \Omega) + (1 / 2) (y - X \beta)^{\prime} \Omega^{- 1} (d^{2} \Omega) \Omega^{- 1} (y - X \beta). \end{aligned} \tag{52} $

Minus the expectation of the second differential takes the simple form

$

• E (d^{2} \Lambda) = (1 / 2) \operatorname{tr} (\Omega^{- 1} d \Omega)^{2} + (d \mu)^{\prime} \Omega^{- 1} (d \mu), $

which implies that the information matrix is block-diagonal in $(\beta, ψ)$ and $\theta$. Therefore we do not have to take the variance of the ML estimator $\hat{\theta}$ into account when calculating the variance of the ML estimators $(\hat{\beta}, \hat{ψ})$. Now writing

$ \operatorname{tr} (\Omega^{- 1} d^{2} \Omega) = (d \operatorname{vec} \Omega)^{\prime} (\Omega^{- 1} \otimes \Omega^{- 1}) (d \operatorname{vec} \Omega) = (d \theta)^{\prime} F_{\theta \theta} , d \theta, $

with

$ F_{\theta \theta} = \left( \frac{\partial \operatorname{vec} \Omega}{\partial \theta^{\prime}} \right)^{\prime} (\Omega^{- 1} \otimes \Omega^{- 1}) \left( \frac{\partial \operatorname{vec} \Omega}{\partial \theta^{\prime}} \right), $

and

$ (d \mu)^{\prime} \Omega^{- 1} (d \mu) = \left( \begin{array}{c} d \beta \ d ψ \end{array} \right)^{\prime} \left( \begin{array}{cc} F_{\beta \beta} & F_{\beta ψ} \ F_{ψ \beta} & F_{ψ ψ} \end{array} \right) \left( \begin{array}{c} d \beta \ d ψ \end{array} \right), $

with

$ F_{\beta \beta} = X^{\prime} \Omega^{- 1} X, \quad F_{\beta ψ}^{} = F_{ψ \beta}^{\prime} = X^{\prime} \left( \beta^{\prime} \otimes \Omega^{- 1} \right) Z, $

and

$ F_{ψ ψ} = Z^{\prime} \left( \beta \beta^{\prime} \otimes \Omega^{- 1} \right) Z, $

we obtain the information matrix

$ F = \left( \begin{array}{ccc} F_{\beta \beta} & F_{\beta ψ} & 0 \ F_{ψ \beta} & F_{ψ ψ} & 0 \ 0 & 0 & F_{\theta \theta} \end{array} \right), $

whose inverse $F^{- 1}$ approximates the variance matrix of the ML estimators.

15. Example 4: The Eckart–Young theorem

Section Summary: The section presents the Eckart–Young theorem as a method for finding the best low-rank approximation to a given matrix A, by minimizing the total squared difference between A and a product XZ' where the factors have a smaller inner dimension r. Using Lagrange multipliers to enforce the required orthogonality constraint on Z, the derivation shows that the solution is obtained by taking Z from the eigenvectors linked to the r largest eigenvalues of A'A, yielding the approximation AZZ' and a minimal error equal to the sum of the remaining smaller eigenvalues. This approach thus extracts the dominant structure in A while quantifying exactly how much information is lost.

My next example is somewhat related to restricted least squares (Section 6), but instead of approximating a vector $A b$, we wish to approximate $A$ itself. More specifically, I wish to approximate a given $m \times n$ matrix $A$ by a matrix $X Z^{\prime}$, such that

$ f (X, Z) = \operatorname{tr} (A - X Z^{\prime}) (A - X Z^{\prime})^{\prime} \tag{53} $

is minimized. Let $r$ be a given (small) integer such that $X$ has dimension $m \times r$ and $Z$ has dimension $n \times r$, where the latter is normalized by the semi-orthogonality condition $Z^{\prime} Z = I_{r}$.

The Eckart–Young theorem (Eckart and Young, 1936, Theorem 17.7 Magnus and Neudecker, 2019) tells us that the minimum of $f$ is obtained when $Z$ contains the eigenvectors associated with the $r$ largest eigenvalues of $A^{\prime} A$ and $X = A Z$. The ‘best’ approximation $\tilde{A}$ (of rank $r$) to $A$ is then $\tilde{A} = A Z Z^{\prime}$, and the constrained minimum of $f$ is the sum of the $n - r$ smallest eigenvalues of $A^{\prime} A$.

To prove the Eckart–Young theorem, define the Lagrangian function

$ L (X, Z) = \frac{1}{2} \operatorname{tr} (A - X Z^{\prime}) (A - X Z^{\prime})^{\prime} - \frac{1}{2} \operatorname{tr} L (Z^{\prime} Z - I_{r}), \tag{54} $

where $L$ is a symmetric $r \times r$ matrix of Lagrange multipliers (see Exercise 18). Differentiating $L$, we obtain

$ \begin{aligned} d L (X, Z) &= \operatorname{tr} (A - X Z^{\prime}) d (A - X Z^{\prime})^{\prime} - \frac{1}{2} \operatorname{tr} L \left( (d Z)^{\prime} Z + Z^{\prime} d Z \right) \ &= - \operatorname{tr} (A - X Z^{\prime}) Z (d X)^{\prime} - \operatorname{tr} (A - X Z^{\prime}) (d Z) X^{\prime} - \operatorname{tr} L Z^{\prime} d Z \ &= - \operatorname{tr} (A - X Z^{\prime}) Z (d X)^{\prime} - \operatorname{tr} (X^{\prime} A - X^{\prime} X Z^{\prime} + L Z^{\prime}) d Z. \end{aligned} $

The first-order conditions are therefore

$ \begin{aligned} (A - X Z^{\prime}) Z = 0, \ X^{\prime} A - X^{\prime} X Z^{\prime} + L Z^{\prime} = 0, \ Z^{\prime} Z = I_{r}. \end{aligned} \tag{55} $

These conditions give $X = A Z$ and hence $L = X^{\prime} X - X^{\prime} A Z = 0$, so that we arrive at the equation

$ (A^{\prime} A) Z = Z (Z^{\prime} A^{\prime} A Z). \tag{56} $

Now, let $P$ be an orthogonal $r \times r$ matrix such that

$ P^{\prime} (Z^{\prime} A^{\prime} A Z) P = \Lambda_{1}, $

where $\Lambda_{1}$ is a diagonal $r \times r$ matrix containing the eigenvalues of $Z^{\prime} A^{\prime} A Z$ on its diagonal. Let $T_{1} = Z P$. Then (56) can be written as $A^{\prime} A T_{1} = T_{1} \Lambda_{1}$, where $T_{1}$ is a semi-orthogonal $n \times r$ matrix (that is, it satisfies $T_{1}^{\prime} T_{1}^{} = I_{r}$) that diagonalizes $A^{\prime} A$, and the $r$ diagonal elements in $\Lambda_{1}$ are eigenvalues of $A^{\prime} A$. (The matrix $A^{\prime} A$ has $n$ eigenvalues and $r$ of these are contained in $\Lambda_{1}$.)

Given $X = A Z$, we have

$ (A - X Z^{\prime}) (A - X Z^{\prime})^{\prime} = A (I - Z Z^{\prime}) A^{\prime}, $

and thus

$ \operatorname{tr} (A - X Z^{\prime}) (A - X Z^{\prime})^{\prime} = \operatorname{tr} A^{\prime} A - \operatorname{tr} \Lambda_{1}. \tag{57} $

To minimize (57), we must maximize $\operatorname{tr} \Lambda_{1}$. Hence $\Lambda_{1}$ must contain the $r$ largest eigenvalues of $A^{\prime} A$, and $T_{1}$ contains eigenvectors associated with these $r$ eigenvalues. The ‘best’ approximation to $A$ is then

$ X Z^{\prime} = A Z Z^{\prime} = A T_{1}^{} T_{1}^{\prime}, $

so that an optimal choice is $Z = T_{1}, : X = A T_{1}$. From (57), it is clear that the value of the constrained minimum is the sum of the $n - r$ smallest eigenvalues of $A^{\prime} A$.

16. Tacit knowledge

Section Summary: Tacit knowledge refers to the practical insights and skills gained through hands-on experience in a field, rather than from formal courses or textbooks. In this section, the author shares a handful of such tips accumulated over years of working with matrices, aimed at making proofs and calculations more efficient. These cover general strategies like testing ideas on small or simplified cases before generalizing, as well as specific techniques for derivatives, differentials, and Hessians in matrix calculus.

When you want to become a carpenter or a doctor or a chef, you follow courses, read books, and pass exams, and at some point you become a trainee in a studio, clinic, or restaurant. There you learn important things that are not available in books and were not taught to you in courses. Such knowledge is called ‘tacit knowledge’. Let me try, in this final section, to write down a few things that I have learned over the years in the hope that they may be of use to the reader. The following suggestions refer to matrix algebra in general.

• • Think of matrices as units of a higher-order algebra. Do not think in terms of the elements of a matrix.
• • When proving a theorem about $n \times n$ matrices, try to prove it first for $n = 2$ and $n = 3$. If it works for $n = 2$, then this is good news, but it is no guarantee. But if it works for $n = 3$, then it probably works in general.
• • When a matrix is symmetric, set it equal to a diagonal matrix, and see if the theorem works in that case.
• • If you need to prove that $A = B$, it is almost always easier to prove that $C = A - B = 0$. There are many ways to prove that $C = 0$, but the method of showing that $c_{i j} = 0$ for all $i$ and $j$ is probably not the most efficient. Maybe you can prove that $\operatorname{tr} C^{\prime} C = 0$, which is equivalent.
• • If we have a matrix $A$ and a matrix $B$ which is the same as $A$ except that (for example) its last column is missing, then do not write $B = \tilde{A}$, but write $B$ explicitly in terms of $A$, in the present case $B = A E$, where $E$ is a matrix of zeros and ones.

The next suggestions refer specifically to matrix calculus.

• • Always use the correct definition of matrix derivative. The derivative of $\operatorname{tr} X$ is not the identity matrix but $(\operatorname{vec} I)^{\prime}$, a row vector.
• • When you need the Hessian, do not start with the first derivative, but with the first differential. Then obtain the second differential, then the Hessian (see the exercises in Section 12).
• • Remember that Cauchy invariance does not hold for second differentials. I usually avoid the chain rule for second differentials, and write the first differential in terms of the eventual matrices, unless there is linearity (symmetry, diagonality), which can be safely added at the end.
• • If $d^{2} f (x) = (d x)^{\prime} B (x) , d x$, remember to ‘symmetrize’ the matrix $B$, because $B$ will not be symmetric in general. The Hessian is then $H f (x) = (B (x) + B (x)^{\prime}) / 2$.

Declaration of competing interest {-}

The author declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments {-}

Section Summary: The author thanks Serena Ng for pointing out the how-to paper series and the editors for backing the tutorial concept on matrix calculus. Two referees and Yannick van Etten receive credit for their helpful suggestions, while Henk Pijls is acknowledged for providing steady mathematical support. The publisher John Wiley is also thanked for permitting reuse of material from the book Matrix Differential Calculus.

I thank Serena Ng for alerting me to the ‘how-to’ paper series, the editors for supporting the idea of a tutorial on matrix calculus, and two unusually constructive referees and Yannick van Etten for valuable comments. As always, Henk Pijls has been my unfailing ‘helpdesk’ for mathematical questions of all types. My publisher John Wiley gracefully permitted me to freely use material from Matrix Differential Calculus, 3rd edition.

References {-}

Section Summary: This section compiles a short bibliography of scholarly works focused on matrix algebra, approximation techniques, and differential calculus, along with their uses in statistics and econometrics. Several entries highlight the contributions of J.R. Magnus, including his books and articles on linear structures and matrix derivatives, while one paper examines real-world data on earthquake risks reflected in Japanese property prices. The references also credit earlier foundational research from the 1930s alongside more recent publications.

• Eckart C., Young G. The approximation of one matrix by another of lower rank. Psychometrika, 1 (1936), pp. 211-218 View in Scopus · Google Scholar

• Ikefuji M., Laeven R.J.A., Magnus J.R., Yue Y. Earthquake risk embedded in property prices: Evidence from five Japanese cities. J. Amer. Statist. Assoc., 117 (2022), pp. 82-93 Crossref · View in Scopus · Google Scholar

• Magnus J.R. Linear Structures. Oxford University Press, New York (1988) Google Scholar

• Magnus J.R. On the concept of matrix derivative. J. Multivariate Anal., 101 (2010), pp. 2200-2206 View PDF · View article · View in Scopus · Google Scholar

• Magnus J.R. Matrix derivatives: Why and where did it go wrong? IMAGE. Bull. Int. Linear Algebra Soc., 72 (Spring) (2024), pp. 3-8 Google Scholar

• Magnus J.R., Neudecker H. Matrix Differential Calculus with Applications in Statistics and Econometrics. (third ed.), John Wiley, Chichester/New York (2019) Google Scholar