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

18.S096 PSET 1 Solutions IAP 2023 February 3, 2023 Problem 0 (4+4+4+4 points) The hyperbolic Corgi notebook may be found at https://mit-c25.netlify.app/notebooks/1_hyperbolic_corgi. Compute the 2 × 2 Jacobian matrix for each of the following image transformations from that notebook: (a) rotate(θ): (x, y) → (cos(θ)x + sin(θ)y, − sin(θ)x + cos(θ)y) Solution: This is simply a linear function from R2 → R2 ! ! ! x cos θ sin θ x → y − sin θ cos θ y | {z } {z } | ⃗ x R(θ) By the same reasoning as in problem 1, the derivative (Jacobian) is simply the rotation operator R(θ): ⃗ and hence the Jacobian is R(θ) . d(R⃗x) = Rdx, (b) hyperbolic_rotate(θ): (x, y) → (cosh(θ)x + sinh(θ)y, sinh(θ)x + cosh(θ)y) Solution: This is another linear transformation: ! ! ! x cosh θ sinh θ x → y sinh θ cosh θ y | {z } | {z } ⃗ x H(θ) with Jacobian H(θ) . (c) nonlin_shear(θ): (x, y) → (x, y + θx2 ) Solution: The differential is: x d y + θx2 ! = ! dx = dy + 2θx dx ! 1 0 2θx 1 dx dy ! so the Jacobian is the boxed matrix. p (d) warp(θ): (x, y) → rotate(θ x2 + y 2 )(x, y) Solution: This is the function ⃗x → R(θ∥⃗x∥)⃗x in terms of the rotation matrix R(θ) from part (a), so we we can use the product rule: ⃗ d (R(θ∥⃗x∥)⃗x) = dR⃗x + Rdx 1 where by the chain rule: dR = R′ (θ∥⃗x∥)d(θ∥⃗x∥) = θR′ (θ∥⃗x∥)d(∥⃗x∥) with R′ (ϕ) = − sin ϕ cos ϕ − cos ϕ − sin ϕ ! by familiar 18.01 derivatives of each component—which follows from the definition dR = R(ϕ + dϕ) − R(ϕ) = R′ (ϕ)dϕ, since the scalar dϕ multiplies R′ elementwise. To get d(∥⃗x∥) we can apply the chain rule again: d(∥⃗x∥) = d((⃗xT ⃗x)1/2 ) = T ⃗ d(⃗xT ⃗x)
2⃗x dx = , 2(⃗xT ⃗x)1/2 2
∥⃗x∥ noting that familiar 18.01 calculus rules work fine when applying the chain rule to scalar terms.1 Hence, putting it all together and rearranging scalar terms (which we can move freely), we have: d(warp ⃗x) = = θ ′ ⃗ + Rdx ⃗ R (θ∥⃗x∥)⃗x⃗xT dx ∥⃗x∥ ⃗x⃗xT θ∥⃗x∥ R (θ∥⃗x∥) T + R(θ∥⃗x∥) ⃗x ⃗x ′ ! ⃗ dx in terms of R and R′ defined above, with the boxed term being the Jacobian, and we have re-arranged terms T ⃗ x⃗ xT to “beautify” the expression by making it clear that ⃗⃗xx⃗xT ⃗x = ∥⃗ x∥2 is an orthogonal projection operator. Problem 1 (5+4 points) (a) Suppose that L[x] is a linear operation (for x in some vector space V , with outputs L[x] in some other vector space W ). If f (x) = L[x] + y for a constant y ∈ W , what is f ′ (x) (in terms of L and/or y)? Solution: This problem is mainly about knowing the definitions of linear operators and derivatives. If f (x) = L[x] + y, then  − y) = L[dx] L[x] − ( df = f (x + dx) − f (x) = ( L[x + dx] +y)
| {z }
 L[x]+L[dx] = so we have f ′ (x)[dx] = L[dx] or equivalently f ′ (x) = L . For affine functions, the derivative is just the linear part. (b) Give the derivatives of f (A) = AT (transpose) and g(A) = 1 + tr A (trace) as special cases of the rule you derived in the previous part. Solution: Again, the key is simply to understand linearity. In both of these examples, we have a linear operator that you cannot easily write as a matrix × vector product (unless you "vectorize" the inputs and/or outputs). (i) f (A) = AT is a linear operator because transposition is linear : (A + B)T = AT + B T and (αA)T = αAT . So, in the notation of part (a), L[x] = AT and y = 0, so f ′ (A)[dA] = (dA)T . Equivalently, d(AT ) = (dA)T . T 1 We can alternatively let r = ∥x∥ =⇒ r 2 = xT x =⇒ 2rdr = d(xT x) = 2xT dx =⇒ dr = 2x dx . But this is basically re-deriving r a rule from first-year calculus. Once we hit a scalar term we needn’t be shy about applying 18.01 rules. 2 (ii) Here, the key is that trace is linear : tr(A + B) = tr A + tr B and tr(αA) = α tr A by inspection of the tr A is an affine function with definition of the trace. So, in the notation of part (a), g(x) = |{z} 1 + |{z} y L[A] ′ g (A)[dA] = tr(dA) , or equivalently d(1 + tr A) = tr(dA) . Problem 2 (5+6+5+5 points) Calculate derivatives of each of the following functions in the requested forms—as a linear operator f ′ (x)[dx], a Jacobian matrix, or a gradient ∇f —as specified in each part. T (a) f (x) = xT (A + diagm(x))2 x, where the inputs x ∈ Rn are vectors, the outputs are scalars,  A = A isa x1   n×n  . x2 constant symmetric n × n matrix ∈ R , and diagm(x) denotes the n × n diagonal matrix   .. . Give the gradient ∇f , such that f ′ (x)dx = (∇f )T dx. Solution: Applying the product rule, we have df = dxT (A + diagm(x))2 x + xT (A + diagm(x))2 dx + xT d(diagm x)(A + diagm(x))x + xT (A + diagm(x)) diagm(dx)x | {z } =diagm(dx) where d(A + diagm(x)) = d(diagm x) since A is a constant, and because diagm is linear (as in problem 1) we have d(diagm x) = diagm(dx). Now, in order to get this in the form ∇f · dx, we neee to move all of our dx factors to the right. The first trick is one we showed in class for a very similar problem: every scalar equals the transpose of itself, giving dxT (A + diagm(x))2 x = [dxT (A + diagm(x))2 x]T = xT (A + diagm(x))2 dx using the fact that A + diagm(x) is symmetric (A = AT was given and diagm x is diagonal). Similarly combining the other pair of terms in df , we get: df = 2xT (A + diagm(x))2 dx + 2xT (A + diagm(x)) diagm(dx)x . The second trick is more subtle: if you think carefully about diagm(dx)x, you will realize that it is simply an elementwise product (denoted by .∗ in Julia), so: diagm(dx)x = dx .∗ x = x .∗ dx = diagm(x)dx Hence   df = 2xT (A + diagm(x))2 + 2xT (A + diagm(x)) diagm(x) dx and ∇f = [· · · ]T therefore gives   ∇f = 2 (A + diagm(x))2 + diagm(x)(A + diagm(x)) x = 2(A + 2 diagm(x))(A + diagm(x))x . (b) f (x) = (A + yxT )−1 b, where the inputs x and outputs f (x) are n-component (column) vectors in Rn , y and b are constant vectors ∈ Rn , and A is a constant n × n matrix ∈ Rn×n . 3 (i) Give f ′ (x) as a Jacobian matrix. Solution: The key here is the formula derived in class for the derivative of a matrix inverse: d(B −1 ) = −B −1 dB B −1 . Applying this to B = A + yxT and dB = y(dx)T , and hence to f (x) via the product rule, gives: df = −(A + yxT )−1 y(dx)T (A + yxT )−1 b {z } | f (x) = −(A + yxT )−1 yf (x)T dx , where we have again used (dx)T f (x) = f (x)T dx to move dx to the right. By inspection, our Jacobian matrix is then the rank-1 matrix: f ′ (x) = −(A + yxT )−1 yf (x)T . (ii) If you are given A−1 , then you can compute (A + yxT )−1 and hence f (x) for any x in ∼ n2 scalararithmetic operations (i.e., roughly proportional to n2 , or in computer-science terms Θ(n2 ) “complexity”), using the “Sherman–Morrison” formula (Google it). Explain how your Jacobian matrix can therefore also be computed in ∼ n2 operations for any x given A−1 (i.e. give a sequence of computational steps, each of which costs no more than ∼ n2 arithmetic). Solution: Since we have (A + yxT )−1 in ∼ n2 operations for any x, we can also use it to compute c = (A + yxT )−1 y by an additional matrix–vector multiplication (∼ n2 scalar arithmetic operations). Our Jacobian is then the outer product (column × row) f ′ (x) = −cf (x)T which requires an additional n2 multiplications (and n negations of c) to yield an n × n matrix. Hence, overall, the whole process requires an operation count that scales proportional to n2 . Note that the order in which we do the operations matters! If we computed it in the order f ′ (x) = −(A + yxT )−1 yf (x)T
we would have had a matrix–matrix multiplication costing ∼ n3 operations, even if the matrix inversion had a cost ∼ n2 . T (c) f (x) = xxxT x , with vector inputs x ∈ Rn and matrix outputs f ∈ Rn×n . Give f ′ (x) as a linear operator, i.e. a linear formula for f ′ (x)[dx].
Solution: We mainly just apply the product rule here, noting that d (xT x)−1 simplifes to the ordinary 4 quotient rule because xT x is a scalar:
d(xxT ) + xxT d (xT x)−1 T x x dx xT + x dxT xxT d(xT x) = − xT x (xT x)2 df = = dx xT + x dxT xxT (xT dx) = f ′ (x)[dx] − 2 xT x (xT x)2 which could be simplified in various ways, but we cannot simply ut all of the dx factors on the right since dxxT ̸= xdxT (very different from the scalar dxT x = xT dx). T (d) g(x) = xxxT x b, with vector inputs x ∈ Rn and vector outputs f ∈ Rn , where b ∈ Rn is a constant vector. Give g ′ (x) as a Jacobian matrix. Solution: We can use the solution from in the previous part since g(x) = f (x)b, but we can simplify it further because dxT b = bT dx, and xT b is a scalar that can be commuted freely, allowing us to move all of the dx factors to the right: dx xT b + x dxT b xxT b(2xT dx) − xT x (xT x)2   xxT bxT 1 T T (x b)I + xb − 2 dx , = T x x xT x | {z } dg = df b = g ′ (x) where I is the n × n identity matrix (since dx(xT b) = (xT b)Idx). This again could be simplified in various ways. Problem 3 (5+5+5 points) (a) Argue briefly that linear functions that map n × n matrices to n × n matrices themselves form a vector space V . What is the dimension of this vector space? Solution: Suppose L1 , L2 ∈ V are two such linear functions. Then this is a vector space if we let L = αL1 + βL2 be the linear map L[X] = αL1 [X] + βL2 [X] for some scalars α, β—it is clear by inspection that L satisfies the axioms of linearity if L1 , L2 do, so this is a vector space (we can add, subtract, and scale). How many parameters does such a map have? It has n2 inputs and n2 outputs, so a linear function has n4 parameters—we could equivalently write an L ∈ V in “vectorized” form as an n2 × n2 matrix multiplying vec(X) to produce vec(L[X]). (b) Argue briefly that linear functions of n × n matrices of the form X → AX, where A is n × n, form a vector space. What is the dimension of this vector space? Solution: This is clearly a subspace of V : if we let LA [X] = AX, then by inspection LA1 ± LA2 = LA1 ±A2 and αLA = LαA using the definitions above. But it is of dimension n2 , the number of parameters in the 5 n × n matrix A. (c) Argue briefly why it follows that there must be infinitely many linear functions ∈ V that are not of the form X → AX. Solution: Since the X → AX functions are an n2 -dimensional subspace of the n4 -dimensional V , it clearly cannot be all of V unless n = 1. Indeed, simply counting dimensions we know that there are n4 − n2 = n2 (n2 − 1) dimensions left. 6