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)))
  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 <- hlibert(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.3&-0.9 \\-1.7&2.9&1.0&1.2 \\-1.3&1.0&9.4&1.3 \\-0.9&1.2&1.3&0.7 \\\end{pmatrix} \]
Solution 
C <- cbind(
  c(3.9, -1.7, -1.3, -.9),
  c(-1.7, 2.9, 1, 1.2),
  c(-1.3, 1, 9.4, 1.3),
  c(-.9, 1.2, 1.3, .7)
)

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 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}^2} = 2 \boldsymbol{\Sigma}^{-1}. \]