Chapter 3

Week 3 lecture 3

Challenges I

  1. Confirm that the following function generates a symmetric matrix for any square matrix X.
make_sym <- function(X) {
  .5 * (X + t(X))
}
Solution 
n <- 3
A <- matrix(rnorm(n * n), n, n)
A2 <- make_sym(A)
all.equal(A2, t(A2))
  1. In R, the function stop() can be used to stop a function when an error occurs. So
if (condition)
  stop('some error message')

Modify make_sym() so that it tests whether X is square, and returns an appropriate error if it isn’t.

Solution 
make_sym_check <- function(X) {
  # function to generate a symmetric matrix
  # X is a square matrix
  if (nrow(X) != ncol(X))
    stop('Supplied matrix X is not square')
  .5 * (X + t(X))
}

make_sym_check(A)
B <- matrix(rnorm(2 * n), n, 2)
make_sym_check(B)
  1. Confirm that crossprod(A, B) and t(A) %*% B give the same result and then benchmark the two commands for \(1000 \times 200\) and \(1000 \times 100\) matrices A and B.
Solution 
A <- matrix(rnorm(1e3 * 2e2), 1e3)
B <- matrix(rnorm(1e3 * 1e2), 1e3)
all.equal(crossprod(A, B), t(A) %*% B)

microbenchmark::microbenchmark(
  crossprod(A, B), 
  t(A) %*% B)
  1. Confirm that \(\mathbf{x}_i^{\text{T}} \mathbf{Ax}_i > 0\) for

\[ \mathbf{A} = \left(\begin{array}{ccc} 0.91 & 0.32 & 0.62\\ 0.32 & 0.31 & 0.15\\ 0.62 & 0.15 & 0.47 \end{array}\right),~ \mathbf{x}_1 = \left(\begin{array}{r} -0.12\\ 0.51\\ 0.22 \end{array}\right),~ \mathbf{x}_2 = \left(\begin{array}{r} -0.32\\ -1.07\\ -1.36 \end{array}\right),~ \mathbf{x}_3 = \left(\begin{array}{r} -0.97\\ -0.16\\ -0.36 \end{array}\right). \]

Solution 
A <- rbind(c(0.91, 0.32, 0.62),
           c(0.32, 0.31, 0.15),
           c(0.62, 0.15, 0.47))
x_1 <- c(-0.12, 0.51, 0.22)
x_2 <- c(-0.32, -1.07, -1.36)
x_3 <- c(-0.97, -0.16, -0.36)

t(x_1) %*% A %*% x_1 > 0
t(x_2) %*% A %*% x_2 > 0
t(x_3) %*% A %*% x_3 > 0
  1. If \(\mathbf{x}_i\) is stored in R as x_i, how can all \(\mathbf{x}_i\) be simultaneously checked by forming X = cbind(x_1, x_2, x_3)?
Solution 
X <- cbind(x_1, x_2, x_3)
colSums(X * (A %*% X))

Week 4 lecture 1

Challenges I

  1. The function dmvn1() is a bit irritating because it returns a \(1 \times 1\) matrix and has a logarithm attribute. Modify the function so that it returns a scalar, and has no attributes. [Consider using as.vector().]
Solution 
dmvn1 <- function(y, mu, Sigma, log = TRUE) {
  # Function to evaluate multivariate Normal pdf
  # at y given mean mu and variance-covariance
  # matrix Sigma.
  # Returns 1x1 matrix, on log scale, if log == TRUE.
  p <- length(y)
  res <- y - mu
  out <- - 0.5 * determinant(Sigma)$modulus - 0.5 * p * log(2 * pi) -
             0.5 * t(res) %*% solve(Sigma) %*% res
  if (!log) 
    out <- exp(out)
  attributes(out) <- NULL
  out
}

y <- c(.7, 1.3, 2.6)
mu <- 1:3
Sigma <- matrix(c(4, 2, 1, 2, 3, 2, 1, 2, 2), 3, 3)
dmvn1(y, mu, Sigma)

Challenges II

  1. Using solve(A, b) in R, where A \(= \mathbf{A}\) and b \(= \mathbf{b}\), confirm that \(\mathbf{x} = (23/2, -22/3, 11/3, -10/3)^\text{T}\) is the solution to \(\mathbf{Ax} = \mathbf{b}\) for \[ \mathbf{A} = \left(\begin{array}{rrrr} 2 & 3 & 0 & 0\\ 4 & 7 & 2 & 0\\ -6 & -10 & 0 & 1\\ 4 & 6 & 4 & 5\\ \end{array}\right) \text{ and } \mathbf{b} = \left(\begin{array}{r} 1\\ 2\\ 1\\ 0\\ \end{array}\right). \]
Solution 
A <- cbind(
  c(2, 4, -6, 4),
  c(3, 7, -10, 6),
  c(0, 2, 0, 4),
  c(0, 0, 1, 5)
)
b <- c(1, 2, 1, 0)
x <- c(23/2, -22/3, 11/3, -10/3)

all.equal(solve(A, b), x)

Week 4 lecture 2

Challenges I

  1. Using backsolve() and forwardsolve() accordingly solve
  • \(\mathbf{Lx}_1 = \mathbf{z}\)
  • \(\mathbf{Ux}_2 = \mathbf{z}\)

for \(\mathbf{x}_1\) and \(\mathbf{x}_2\), where

\[ \hspace{-.5cm}\mathbf{L} = \left( \begin{array}{cccc} 2.0 & 0.0 & 0.0 & 0.0\\ 2.3 & 1.2 & 0.0 & 0.0\\ 0.8 & 2.1 & 1.8 & 0.0\\ 0.3 & 1.4 & 1.1 & 0.5 \end{array} \right) ~,\mathbf{U} = \left(\begin{array}{cccc} 0.3 & 1.9 & 1.9 & 1.4\\ 0.0 & 2.5 & 0.6 & 0.1\\ 0.0 & 0.0 & 2.5 & 0.5\\ 0.0 & 0.0 & 0.0 & 2.1 \end{array}\right) ~\mathbf{z} = \left(\begin{array}{r} -0.47\\ 0.07\\ -0.48 \\1.45 \end{array}\right). \]

and confirm your result with solve()

Solution 
x1 <- forwardsolve(L, z)
x1
solve(L, z)

x2 <- backsolve(U, z)
x2
solve(U, z)

Challenges II

  1. Use the following function to generate a \(5 \times 5\) symmetric and positive definite matrix.
rsympdmat <- function(n) {
# function to generate symmetric positive definite matrix
# of dimension n x n
# n is an integer
A <- matrix(rnorm(n * n), n)
crossprod(A, diag(abs(rnorm(n))) %*% A)
}
  • Calculate its Cholesky decomposition and confirm that it is a Cholesky decomposition, i.e. that its cross product or transposed cross product give the required result
  • Invert it through its Cholesky decomposition, and confirm your result with solve()
  • Compute the logarithm of its determinant, and confirm your result with determinant()
Solution 
A <- rsympdmat(5)
U <- chol(A)
all.equal(crossprod(U), A)
all(U[lower.tri(U)] == 0)
all(diag(U) > 0)

# Inverse of A
iA <- solve(A)
all.equal(chol2inv(U), iA)

# log determinant of A
determinant(A)
2 * sum(log(diag(U)))

See here for an RMarkdown example embedding some of the above code.

  1. The Cholesky decomposition exists for any symmetric positive definite matrix. Generate an arbitrary \(4 \times 4\) matrix, use make_sym() to make it symmetric and then consider how chol() can check whether it’s positive definite.
make_sym <- function(X) {
  # function to generate a symmetric matrix
  # X is a square matrix
  .5 * (X + t(X))
}
Solution 
A <- matrix(rnorm(16), 4, 4)

make_sym <- function(X) {
  # function to generate a symmetric matrix
  # X is a square matrix
  .5 * (X + t(X))
}

A_sym <- make_sym(A)
chol(A_sym)

Week 4 lecture 3

Challenges I

  1. Compute the Cholesky decomposition of \(\mathbf{H}_4\), i.e. the \(4 \times 4\) Hilbert matrix, and then use this to solve \(\mathbf{H}_4 \mathbf{x} = \mathbf{b}\) for \(\mathbf{x}\) where \[ \mathbf{b} = \left(\begin{array}{c} 0.04\\ 3.19\\ 3.26\\ 1.93 \end{array}\right) \] and the following function can be used to form the Hilbert matrix
hilbert <- function(n) {
  # Function to evaluate n by n Hilbert matrix.
  # Returns n by n matrix.
  ind <- 1:n
  1 / (outer(ind, ind, FUN = '+') - 1)
}
Solution 
H4 <- hilbert(4)
b <- c(.04, 3.19, 3.26, 1.93)

x1 <- solve(H4, b)
U <- chol(H4)
x2 <- forwardsolve(U, 
  backsolve(U, b, transpose = TRUE), 
  upper.tri = TRUE)

all.equal(x1, x2)

Week 5 lecture 1

Challenges I

  1. Compute the eigen-decomposition of \(\boldsymbol{\Sigma}\) from Example 3.2, i.e.

\[ \boldsymbol{\Sigma} = \begin{pmatrix}4&2&1 \\2&3&2 \\1&2&2 \\\end{pmatrix}. \]

Solution 
Sigma <- matrix(c(4, 2, 1, 2, 3, 2, 1, 2, 2), 3, 3)
eS <- eigen(Sigma, symmetric = TRUE)
eS
  1. Show that its eigenvectors are orthogonal.
Solution 
U <- eS$vectors
all.equal(crossprod(U), diag(3))
  1. Find \(|\boldsymbol{\Sigma}|\) and \(\log(|\boldsymbol{\Sigma}|)\), avoiding first calculating \(|\boldsymbol{\Sigma}|\) for the latter.
Solution 
lambda <- eS$values
detS <- prod(lambda)
all.equal(detS, det(Sigma))
log(detS)
sum(log(lambda))
determinant(Sigma)

Challenges II

  1. Write two functions, matpow1(A, k) and matpow2(A, k), that compute \(\mathbf{A}^k\) for a positive integer \(k\) and a symmetric matrix \(\mathbf{A}\), where matpow1() uses a for() loop and matpow2() uses the eigen-decomposition of \(\mathbf{A}\).
matpow1 <- function(A, k) {
  # function to calculate matrix power of symmetric matrix A
  # based on a for loop
  # k is an integer
  # return matrix of same dimension as A
  ** insert code here **
}

matpow2 <- function(A, k) {
  # as above, but based on the eigen-decomposition of A
  ** insert code here **
}
Solution 
# one loop-based approach
matpow1 <- function(A, k) {
  # function to calculate matrix power of symmetric matrix A
  # based on a for loop
  # k is an integer
  # return matrix of same dimension as A
  Ak <- A
  if (k > 1) {
    for (i in 2:k) {
      Ak <- Ak %*% A
    }
  }
  Ak
}

# an alternative loop-based approach
matpow1 <- function(A, k) {
  # function to calculate matrix power of symmetric matrix A
  # based on a for loop
  # k is an integer
  # return matrix of same dimension as A
  Ak <- diag(nrow(A))
  for (i in 1:k) {
      Ak <- Ak %*% A
  }
  Ak
}

H <- hilbert(4)
Hpow5 <- matpow1(H, 5)
all.equal(matpow1(H, 5), Hpow5)

matpow2 <- function(A, k) {
  # as above, but based on the eigen-decomposition of A
  eA <- eigen(A, symmetric = TRUE)
  U <- eA$vectors
  lambda <- eA$values
  U %*% tcrossprod(diag(lambda^k), U)
}

all.equal(matpow2(H, 5), Hpow5)

Week 5 lecture 2

Challenges I

  1. Consider Example 3.22 from the lecture notes. Use the SVD of \(\boldsymbol{\Sigma}\) to solve \(\boldsymbol{\Sigma} \mathbf{z} = \mathbf{y} - \boldsymbol{\mu}\) for \(\mathbf{z}\), where

\[ {\bf y} = \left( \begin{array}{c} 0.7\\ 1.3\\ 2.6 \end{array} \right),~~ \boldsymbol{\mu} = \left( \begin{array}{c} 1\\ 2\\ 3 \end{array} \right)~~\text{and}~~ \boldsymbol{\Sigma} = \left( \begin{array}{rrr} 4 & 2 & 1\\ 2 & 3 & 2\\ 1 & 2 & 2\\ \end{array} \right). \]

Solution 
y <- c(.7, 1.3, 2.6)
mu <- 1:3
Sigma <- matrix(c(4, 2, 1, 2, 3, 2, 1, 2, 2), 3, 3)

S.svd <- svd(Sigma)
S.d <- S.svd$d
S.V <- S.svd$v

b2 <- crossprod(S.V, y - mu)
x2 <- b2 / S.d
z <- S.V %*% x2

Challenges II

  1. Use SVD to compute the rank-5 representation of \(\mathbf{H}_7\), the \(7 \times 7\) Hilbert matrix.
Solution 
H <- hilbert(7)
svdH <- svd(H)
r <- 5
U_r <- svdH$u[, 1:r]
V_r <- svdH$v[, 1:r]
D_r <- svdH$d[1:r]
H_r <- U_r %*% (D_r * t(V_r))
  1. Compare the rank-5 representation to \(\mathbf{H}_7\).
Solution 
range(H - H_r)
sum(diag(solve(H, H_r))

Week 5 lecture 3

Challenges I

  1. Compute the QR decomposition of \(\mathbf{H}_4\), the \(4 \times 4\) Hilbert matrix, in the form \(\mathbf{QR} = \mathbf{H}_4\) by finding \(\mathbf{Q}\) and \(\mathbf{R}\).
Solution 
hilbert <- function(n) {
  # Function to evaluate n by n Hilbert matrix.
  # Returns n by n matrix.
  ind <- 1:n
  1 / (outer(ind, ind, FUN = '+') - 1)
}
H4 <- hilbert(4)

qrH4 <- qr(H4)
Q <- qr.Q(qrH4)
R <- qr.R(qrH4)
all.equal(Q %*% R, H4)
  1. Use the QR decomposition to find \(|\text{det}(\mathbf{H}_4)|\).
Solution 
prod(abs(diag(R)))
det(H4)
  1. Use the QR decomposition to solve \(\mathbf{H}_4 \mathbf{x} = \mathbf{b}\) for \(\mathbf{x}\) where \[ \mathbf{b} = \begin{pmatrix}0.04 \\3.19 \\3.26 \\1.93 \\\end{pmatrix}. \]
Solution 
b <- c(.04, 3.19, 3.26, 1.93)
backsolve(R, crossprod(Q, b))
qr.solve(H4, b)
qr.solve(qrH4, b)

Challenges II

  1. Generate an arbitrary \(m \times n\) matrix \(\mathbf{A}\) comprising positive values for \(m = 5\) and \(n = 3\). Then use function qr() in R to find \(\mathbf{Q}\) and \(\mathbf{R}\) of its QR decomposition. Confirm that \(\mathbf{QR} = \mathbf{A}\).
Solution 
m <- 5
n <- 3
A <- matrix(rexp(m * n), m , n)
qrA <- qr(A)
Q <- qr.Q(qrA)
R <- qr.R(qrA)
all.equal(Q %*% R, A)
  1. Investigate what R has actually given with qr(), given Theorem 3.3.
Solution 
# No sketch solution given
  1. Calculate the SVD of your matrix \(\mathbf{A}\) from Question 1, and confirm that the absolute value of the product of the diagonal entries of \(\mathbf{R}\) matches the product of the singular values of its SVD.
Solution 
svdA <- svd(A)
all.equal(prod(svdA$d), abs(prod(diag(R))))
  1. Compute \(\mathbf{A}^{\text{T}}\mathbf{A}\) and confirm that \(\mathbf{R}^{\text{T}}\mathbf{R} = \mathbf{A}^{\text{T}}\mathbf{A}\), so that \(\mathbf{R}\) is a upper-triangular Cholesky decomposition of \(\mathbf{A}^{\text{T}}\mathbf{A}\).
Solution 
crossA <- crossprod(A)
all.equal(crossprod(R), crossA)
  1. Compute the Cholesky decomposition \(\mathbf{LL}^{\text{T}} = \mathbf{A}^{\text{T}}\mathbf{A}\). What do you notice about the values of \(\mathbf{R}\) and \(\mathbf{L}^{\text{T}}\)?
Solution 
L <- t(chol(crossA))
all.equal(abs(R), t(L))

Week 6 lecture 1

Challenges I

  1. Generate an arbitrary \(m \times n\) matrix \(\mathbf{A}\) comprising positive values for \(m = 5\) and \(n = 3\). Then use function qr() in R to find \(\mathbf{Q}\) and \(\mathbf{R}\) of its QR decomposition. Confirm that \(\mathbf{QR} = \mathbf{A}\).
Solution 
m <- 5
n <- 3
A <- matrix(rexp(m * n), m , n)
qrA <- qr(A)
Q <- qr.Q(qrA)
R <- qr.R(qrA)
all.equal(Q %*% R, A)
  1. Investigate what R has actually given with qr(), given Theorem 3.3.
Solution 
# no sketch solution given
  1. Calculate the SVD of your matrix \(\mathbf{A}\) from Question 1, and confirm that the absolute value of the product of the diagonal entries of \(\mathbf{R}\) matches the product of the singular values of its SVD.
Solution 
svdA <- svd(A)
all.equal(prod(svdA$d), abs(prod(diag(R))))
  1. Compute \(\mathbf{A}^{\text{T}}\mathbf{A}\) and confirm that \(\mathbf{R}^{\text{T}}\mathbf{R} = \mathbf{A}^{\text{T}}\mathbf{A}\), so that \(\mathbf{R}\) is a upper-triangular Cholesky decomposition of \(\mathbf{A}^{\text{T}}\mathbf{A}\).
Solution 
crossA <- crossprod(A)
all.equal(crossprod(R), crossA)
  1. Compute the Cholesky decomposition \(\mathbf{LL}^{\text{T}} = \mathbf{A}^{\text{T}}\mathbf{A}\). What do you notice about the values of \(\mathbf{R}\) and \(\mathbf{L}^{\text{T}}\)?
Solution 
L <- t(chol(crossA))
all.equal(abs(R), t(L))

Challenges II

  1. Use Woodbury’s formula to compute \(\mathbf{B}^{-1}\) where \(\mathbf{B} = \mathbf{H}_8 + \mathbf{UU}^\text{T}\) where \[ \mathbf{U} = \begin{pmatrix}4.7&10.6 \\3.2&22.5 \\0.2&14.1 \\1.1&11.5 \\5.0&3.3 \\8.3&12.9 \\1.1&9.0 \\5.6&3.9 \\\end{pmatrix} \] and \(\mathbf{H}_8\) is the \(8 \times 8\) Hilbert matrix.
Solution 
A <- hilbert(8)

U <- cbind(c(4.7, 3.2, 0.2, 1.1, 5, 8.3, 1.1, 5.6), 
           c(10.6, 22.5, 14.1, 11.5, 3.3, 12.9, 9, 3.9))

wsm <- function(iA, U, C = diag(ncol(U)), V) {
  # function to compute Woodbury-Sherman-Morrison formula
  # iA is a m x m matrix
  # U and V are m x n matrices
  # C is a n x n matrix, and defaults to nrow(U) identity matrix
  # returns m x m matrix
  iC <- solve(C)
  iAU <- iA %*% U
  iA - iAU  %*% solve(iC + crossprod(V, iAU), crossprod(V, iA))
}

wsm(solve(A), U, V = U)
  1. Confirm that Woodbury’s formula gives the same result as directly inverting \(\mathbf{B}\) with solve() (or equivalent).
Solution 
B <- A + tcrossprod(U)
iB <- solve(B)
all.equal(iB, wsm(solve(A), U, V = U))
  1. Use the Sherman-Morrison-Woodbury formula to compute \(\mathbf{D}^{-1}\) where \(\mathbf{D} = \mathbf{A} + \mathbf{UCU}^\text{T}\) with \(\mathbf{A}\) and \(\mathbf{U}\) as in Question 1 and \[ \mathbf{C} = \begin{pmatrix}3.9&-1.7 \\-1.7&2.9 \\\end{pmatrix} \]
Solution 
C <- cbind(
  c(3.9, -1.7),
  c(-1.7, 2.9)
)

wsm(solve(A), U, C, U)
  1. Confirm that the Sherman-Morrison-Woodbury formula gives the same result as directly inverting \(\mathbf{D}\)
Solution 
D <- A + U %*% tcrossprod(C, U)
all.equal(solve(D), wsm(solve(A), U, C, U))

Challenges III

  1. Consider the \(MVN_p(\boldsymbol{\mu}, \boldsymbol{\Sigma})\) likelihood’s squared Mahalanobis distance, i.e. \[ g(\mathbf{y}) = (\mathbf{y} - \boldsymbol{\mu})^{\text{T}}\boldsymbol{\Sigma}^{-1}(\mathbf{y} - \boldsymbol{\mu}). \] Find \(\dfrac{\partial g(\mathbf{y}) }{\partial \mathbf{y}}\) and \(\dfrac{\partial^2 g(\mathbf{y})}{\partial \mathbf{y} \partial \mathbf{y}^{\text{T}}}\).
Solution

\[ \dfrac{\partial g(\mathbf{y}) }{\partial \mathbf{y}} = 2 \boldsymbol{\Sigma}^{-1}(\mathbf{y} - \boldsymbol{\mu}),~~\dfrac{\partial^2 g(\mathbf{y}) }{\partial \mathbf{y} \mathbf{y}^\top} = 2 \boldsymbol{\Sigma}^{-1}. \]