18.S096 (IAP 2023) Matrix Calculus for Machine Learning and Beyond: Problem Set 2 Solution

18.S096 Pset 2 SOLUTIONS, IAP 2023 Problem 1 (5+5+5 points) Suppose that A(p) takes a vector p ∈ Rn−1 and returns  a1 p1   p1 a2   A(p) =  p2     the n × n tridiagonal real-symmetric matrix  p2 .. . .. . .. . an−1 pn−1 pn−1 an     ,     where a ∈ Rn−1 is some constant vector. Now, define a scalar-valued function f (p) by
2 f (p) = cT A(p)−1 b for some constant vectors b, c ∈ Rn (assuming we choose p and a so that A is invertible). Note that, in practice A(p)−1 b is not computed by explicitly inverting the matrix A—instead, it can be computed in Θ(n) (i.e., roughly proportional to n) arithmetic operations using Gaussian elimination that takes advantage of the “sparsity” of A (the pattern of zero entries), a “tridiagonal solve”. (a) Write down a formula for computing ∂f /∂p1 (in terms of matrix–vector products and matrix inverses). (Hint: once you know df in terms of dA, you can get ∂f /∂p1 by “dividing” both sides by ∂p1 , so that dA becomes ∂A/∂p1 .) Solution: From the chain rule and the formula for the differential of a matrix inverse, we have df = −2(cT A−1 b)cT A−1 dA A−1 b (noting that cT A−1 b is a scalar so we can commute it as needed). Hence ∂f ∂A −1 = −2(cT A−1 b)cT A−1 A b | {z } ∂p1 | {z } ∂p1 x vT   0 1    1 0 0      . . T  . .  x = v1 x2 + v2 x1 , . . =v  0    ..  . 0 0    0 0 {z } | ∂A ∂p1 where we have simplified the result in terms of x and v for the next part. 1 (b) Outline a sequence of steps to compute both f and ∇f (with respect to p) using only two tridiagonal solves x = A−1 b and an “adjoint” solve v = A−1 (something), plus Θ(n) (i.e., roughly proportional to n) additional arithmetic operations. Solution: Using the notation from the previous part, exploiting the fact that AT = A, we can choose v = A−1 [−2(cT x)c] , which is a single tridiagonal solve. Given x and v, the results of our two Θ(n) tridiagonal solves, we can compute each component of the gradient similar to above by ∂f /∂pk = vk xk+1 + vk+1 xk for k = 1, . . . , n − 1, which costs Θ(1) arithmetic per k and hence Θ(n) arithmetic to obtain all of ∇f . (c) Write a program implementing your ∇f procedure (in Julia, Python, Matlab, or any language you want) from the previous part. (You don’t need to use a fancy tridiagonal solve if you don’t know how to do this in your language; you can solve A−1 (vector) inefficiently if needed using your favorite matrix libraries.) Implement a finite-difference test: Choose a, b, c, p at random, and check that ∇f · δp ≈ f (p + δp) − f (p) (to a few digits) for a randomly chosen small δp. Solution: See accompanying Julia notebook Problem 2 (5+5 points) Suppose that we have a two-argument function f (x, y), where x, y and f may belong to arbitrary vector (Banach) ∂f spaces. Let’s define “partial” derivatives fx and fy (also denoted ∂f ∂x and ∂y ) by the linearization: df = f (x + dx, y + dy) − f (x, y) = fx (x, y)[dx] + fy (x, y)[dy], implicitly dropping higher-order terms as usual. Compute the partial derivatives of the following functions: (a) f (A, x) = A−1 x for n × n matrices A ∈ Rn×n and vectors x ∈ Rn : give fA as a linear operator, and fx as a Jacobian matrix. Solution: By the product rule: −1 df = −A−1 dA A−1 x + A dx , | {z } | {z } fA [dA] fx [dx] so fA [dA] = −A−1 dA A−1 x is a linear operator (input = dA, output = vector) and fx = A−1 is the Jacobian matrix with respect to x. (b) f (A, B) = tr(AT BA), for matrices A, B ∈ Rn×n : give the gradients ∇A f and ∇B f such that fA [dA] = ∇A f · dA and fB [dB] = ∇B f · dB under the Frobenius inner product X · Y = tr(X T Y ) = tr(Y T X). Solution: By the product rule and the usual trace properties (tr XY = tr Y X, tr X = tr X T , tr(X + Y ) = tr X + tr Y ): df = tr(dAT BA) + tr(AT dB A) + tr(AT B dA) = tr(AT B T dA) + tr(AAT dB) + tr(AT B dA) = tr(AT (B + B T ) dA) + tr(AAT dB) so we have ∇A f = (B + B T )A and ∇B f = AAT . 2 Problem 3 (5+5 points) If S is an m×m real-symmetric matrix with a “simple” (multiplicity = 1) eigenvalue λ and corresponding eigenvector q (Sq = λq), normalized to q T q = 1, then the “Hellman–Feynman theorem” states that dλ = q T dS q for a change dS in the matrix S. (a) Derive the Hellman–Feynman theorem by considering the differentials of both sides of the equations d(λ = q T Sq) and d(q T q = 1). Solution: By the product rule, and the eigen-equation Sq = λq, we get dλ = dq Sq + q T dS q + q T S dq = λ (dq q + q T dq) +q T dS q | {z } =d(q T q)=d(1)=0 T = q dS q . Q.E.D. (b) What is the gradient ∇λ with respect to S, for the usual Frobenius inner product ∇λ · dS = tr((∇λ)T dS) Solution: We use the fact that dλ = tr(dλ) since it is a scalar, combined with the cyclic property of the trace, to obtain: dλ = tr(dλ) = tr(q T dS q) = tr(qq T dS) and hence ∇λ = (qq T )T = qq T . Problem 4 (6+6 points) The Jacobian determinant (sometimes called simply “the Jacobian,” clashing with the concept of the Jacobian ∂fi )1≤i,j≤n matrix) is the determinant of the Jacobian matrix. Specifically if f (x) is a function from Rn to Rn and ( ∂x j ′ ′ is the Jacobian matrix f (x), then its determinant det f (x) is the Jacobian determinant. Sometimes we take the absoute value and not worry too much about the sign. (a) The Jacobian determinant represents the local scaling of volume. Compute the Jacobian determinant of the hyperbolic rotation defined in Pset 1, problem 1b, in simplest form. Use this to describe how a little square around a point generally transforms with a hyperbolic rotation. Solution: Recall that “hyperbolic rotation” from pset 1 was defined by the linear transformation ! ! ! x cosh θ sinh θ x → y y sinh θ cosh θ | {z } | {z } ⃗ x H(θ) with Jacobian H(θ), so its Jacobian determinant is simply det H(θ) = cosh2 (θ) − sinh2 (θ) = 1 . This means that the transformation preserves area, i.e. an infinitesimal square around a point is transformed to a rhombus with the same area. (Why a rhombus? Because the columns of H, corresponding to the edges 3 of the transformed square, have equal length but are not orthogonal.) (b) There are many ways to equivalently take a scalar function f (α) and extend it to a matrix function F (M ), which takes in a square matrix and returns a square matrix of the same size. The simplest is to define f (M ) = Xf (Λ)X −1 , where M = XΛX −1 is an eigen-decomposition of M (and use continuity to include non-diagonalizable matrices). Here, f (Λ) denotes the application of a scalar function f (λ) to the eigenvalues λ (on the diagonal of Λ). (e.g., you’ve probably seen eM defined in terms of eλ .) One could then write f ′ (M ) as an explicit n2 ×n2 Jacobian matrix (e.g. via vec(dM ) and Kronecker products), and could then compute its determinant. (i) Write a computer program (in any language) to find the 9 × 9 Jacobian matrix of f (M ) and then the Jacobian determinant by either finite differences or by using automatic differentiation, for f (λ) being eλ , λ2 , and sin(λ) on the 3 × 3 matrix M = [0 1 4; 1 0 1; 4 1 0] with entries Mi,j = (i − j)2 . Solution: The Jacobian determinants should be about 939.059, 4096, and −8.41346 × 10−6 , respectively. See accompanying Julia notebook. (ii) Compare with the following known theoretical formula for the Jacobian determinant for a scalar function f (λ) applied to a diagonalizable matrix M , in terms of M ’s eigenvalues λ: Q 2 Y i