Chapter 3
Week 3 lecture 3
Challenges I
- Confirm that the following function generates a symmetric matrix for any square matrix
X
.
<- function(X) {
make_sym 5 * (X + t(X))
. }
Solution
<- 3
n <- matrix(rnorm(n * n), n, n)
A <- make_sym(A)
A2 all.equal(A2, t(A2))
- In
R
, the functionstop()
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
<- function(X) {
make_sym_check # 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)
<- matrix(rnorm(2 * n), n, 2)
B make_sym_check(B)
- Confirm that
crossprod(A, B)
andt(A) %*% B
give the same result and then benchmark the two commands for \(1000 \times 200\) and \(1000 \times 100\) matricesA
andB
.
Solution
<- matrix(rnorm(1e3 * 2e2), 1e3)
A <- matrix(rnorm(1e3 * 1e2), 1e3)
B all.equal(crossprod(A, B), t(A) %*% B)
::microbenchmark(
microbenchmarkcrossprod(A, B),
t(A) %*% B)
- 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
<- rbind(c(0.91, 0.32, 0.62),
A c(0.32, 0.31, 0.15),
c(0.62, 0.15, 0.47))
<- c(-0.12, 0.51, 0.22)
x_1 <- c(-0.32, -1.07, -1.36)
x_2 <- c(-0.97, -0.16, -0.36)
x_3
t(x_1) %*% A %*% x_1 > 0
t(x_2) %*% A %*% x_2 > 0
t(x_3) %*% A %*% x_3 > 0
- If \(\mathbf{x}_i\) is stored in
R
asx_i
, how can all \(\mathbf{x}_i\) be simultaneously checked by formingX = cbind(x_1, x_2, x_3)
?
Solution
<- cbind(x_1, x_2, x_3)
X colSums(X * (A %*% X))
Week 4 lecture 1
Challenges I
- The function
dmvn1()
is a bit irritating because it returns a \(1 \times 1\) matrix and has alogarithm
attribute. Modify the function so that it returns a scalar, and has no attributes. [Consider usingas.vector()
.]
Solution
<- function(y, mu, Sigma, log = TRUE) {
dmvn1 # 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.
<- length(y)
p <- y - mu
res <- - 0.5 * determinant(Sigma)$modulus - 0.5 * p * log(2 * pi) -
out 0.5 * t(res) %*% solve(Sigma) %*% res
if (!log)
<- exp(out)
out attributes(out) <- NULL
out
}
<- c(.7, 1.3, 2.6)
y <- 1:3
mu <- matrix(c(4, 2, 1, 2, 3, 2, 1, 2, 2), 3, 3)
Sigma dmvn1(y, mu, Sigma)
Challenges II
- Using
solve(A, b)
inR
, whereA
\(= \mathbf{A}\) andb
\(= \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
<- cbind(
A c(2, 4, -6, 4),
c(3, 7, -10, 6),
c(0, 2, 0, 4),
c(0, 0, 1, 5)
)<- c(1, 2, 1, 0)
b <- c(23/2, -22/3, 11/3, -10/3)
x
all.equal(solve(A, b), x)
Week 4 lecture 2
Challenges I
- Using
backsolve()
andforwardsolve()
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
<- forwardsolve(L, z)
x1
x1solve(L, z)
<- backsolve(U, z)
x2
x2solve(U, z)
Challenges II
- Use the following function to generate a \(5 \times 5\) symmetric and positive definite matrix.
<- function(n) {
rsympdmat # function to generate symmetric positive definite matrix
# of dimension n x n
# n is an integer
<- matrix(rnorm(n * n), n)
A 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
<- rsympdmat(5)
A <- chol(A)
U all.equal(crossprod(U), A)
all(U[lower.tri(U)] == 0)
all(diag(U) > 0)
# Inverse of A
<- solve(A)
iA all.equal(chol2inv(U), iA)
# log determinant of A
determinant(A)
2 * sum(log(diag(U)))
- 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 howchol()
can check whether it’s positive definite.
<- function(X) {
make_sym # function to generate a symmetric matrix
# X is a square matrix
5 * (X + t(X))
. }
Solution
<- matrix(rnorm(16), 4, 4)
A
<- function(X) {
make_sym # function to generate a symmetric matrix
# X is a square matrix
5 * (X + t(X))
.
}
<- make_sym(A)
A_sym chol(A_sym)
Week 4 lecture 3
Challenges I
- 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
<- function(n) {
hilbert # Function to evaluate n by n Hilbert matrix.
# Returns n by n matrix.
<- 1:n
ind 1 / (outer(ind, ind, FUN = '+') - 1)
}
Solution
<- hilbert(4)
H4 <- c(.04, 3.19, 3.26, 1.93)
b
<- solve(H4, b)
x1 <- chol(H4)
U <- forwardsolve(U,
x2 backsolve(U, b, transpose = TRUE),
upper.tri = TRUE)
all.equal(x1, x2)
Week 5 lecture 1
Challenges I
- 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
<- matrix(c(4, 2, 1, 2, 3, 2, 1, 2, 2), 3, 3)
Sigma <- eigen(Sigma, symmetric = TRUE)
eS eS
- Show that its eigenvectors are orthogonal.
Solution
<- eS$vectors
U all.equal(crossprod(U), diag(3))
- Find \(|\boldsymbol{\Sigma}|\) and \(\log(|\boldsymbol{\Sigma}|)\), avoiding first calculating \(|\boldsymbol{\Sigma}|\) for the latter.
Solution
<- eS$values
lambda <- prod(lambda)
detS all.equal(detS, det(Sigma))
log(detS)
sum(log(lambda))
determinant(Sigma)
Challenges II
- Write two functions,
matpow1(A, k)
andmatpow2(A, k)
, that compute \(\mathbf{A}^k\) for a positive integer \(k\) and a symmetric matrix \(\mathbf{A}\), wherematpow1()
uses afor()
loop andmatpow2()
uses the eigen-decomposition of \(\mathbf{A}\).
<- function(A, k) {
matpow1 # 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 **
}
<- function(A, k) {
matpow2 # as above, but based on the eigen-decomposition of A
** insert code here **
}
Solution
# one loop-based approach
<- function(A, k) {
matpow1 # 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
<- A
Ak if (k > 1) {
for (i in 2:k) {
<- Ak %*% A
Ak
}
}
Ak
}
# an alternative loop-based approach
<- function(A, k) {
matpow1 # 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
<- diag(nrow(A))
Ak for (i in 1:k) {
<- Ak %*% A
Ak
}
Ak
}
<- hilbert(4)
H <- matpow1(H, 5)
Hpow5 all.equal(matpow1(H, 5), Hpow5)
<- function(A, k) {
matpow2 # as above, but based on the eigen-decomposition of A
<- eigen(A, symmetric = TRUE)
eA <- eA$vectors
U <- eA$values
lambda %*% tcrossprod(diag(lambda^k), U)
U
}
all.equal(matpow2(H, 5), Hpow5)
Week 5 lecture 2
Challenges I
- 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
<- c(.7, 1.3, 2.6)
y <- 1:3
mu <- matrix(c(4, 2, 1, 2, 3, 2, 1, 2, 2), 3, 3)
Sigma
<- svd(Sigma)
S.svd <- S.svd$d
S.d <- S.svd$v
S.V
<- crossprod(S.V, y - mu)
b2 <- b2 / S.d
x2 <- S.V %*% x2 z
Challenges II
- Use SVD to compute the rank-5 representation of \(\mathbf{H}_7\), the \(7 \times 7\) Hilbert matrix.
Solution
<- hilbert(7)
H <- svd(H)
svdH <- 5
r <- svdH$u[, 1:r]
U_r <- svdH$v[, 1:r]
V_r <- svdH$d[1:r]
D_r <- U_r %*% (D_r * t(V_r)) H_r
- 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
- 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
<- function(n) {
hilbert # Function to evaluate n by n Hilbert matrix.
# Returns n by n matrix.
<- 1:n
ind 1 / (outer(ind, ind, FUN = '+') - 1)
}<- hilbert(4)
H4
<- qr(H4)
qrH4 <- qr.Q(qrH4)
Q <- qr.R(qrH4)
R all.equal(Q %*% R, H4)
- Use the QR decomposition to find \(|\text{det}(\mathbf{H}_4)|\).
Solution
prod(abs(diag(R)))
det(H4)
- 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
<- c(.04, 3.19, 3.26, 1.93)
b backsolve(R, crossprod(Q, b))
qr.solve(H4, b)
qr.solve(qrH4, b)
Challenges II
- Generate an arbitrary \(m \times n\) matrix \(\mathbf{A}\) comprising positive values for \(m = 5\) and \(n = 3\). Then use function
qr()
inR
to find \(\mathbf{Q}\) and \(\mathbf{R}\) of its QR decomposition. Confirm that \(\mathbf{QR} = \mathbf{A}\).
Solution
<- 5
m <- 3
n <- matrix(rexp(m * n), m , n)
A <- qr(A)
qrA <- qr.Q(qrA)
Q <- qr.R(qrA)
R all.equal(Q %*% R, A)
- Investigate what
R
has actually given withqr()
, given Theorem 3.3.
Solution
# No sketch solution given
- 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
<- svd(A)
svdA all.equal(prod(svdA$d), abs(prod(diag(R))))
- 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
<- crossprod(A)
crossA all.equal(crossprod(R), crossA)
- 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
<- t(chol(crossA))
L all.equal(abs(R), t(L))
Week 6 lecture 1
Challenges I
- Generate an arbitrary \(m \times n\) matrix \(\mathbf{A}\) comprising positive values for \(m = 5\) and \(n = 3\). Then use function
qr()
inR
to find \(\mathbf{Q}\) and \(\mathbf{R}\) of its QR decomposition. Confirm that \(\mathbf{QR} = \mathbf{A}\).
Solution
<- 5
m <- 3
n <- matrix(rexp(m * n), m , n)
A <- qr(A)
qrA <- qr.Q(qrA)
Q <- qr.R(qrA)
R all.equal(Q %*% R, A)
- Investigate what
R
has actually given withqr()
, given Theorem 3.3.
Solution
# no sketch solution given
- 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
<- svd(A)
svdA all.equal(prod(svdA$d), abs(prod(diag(R))))
- 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
<- crossprod(A)
crossA all.equal(crossprod(R), crossA)
- 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
<- t(chol(crossA))
L all.equal(abs(R), t(L))
Challenges II
- 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
<- hilbert(8)
A
<- cbind(c(4.7, 3.2, 0.2, 1.1, 5, 8.3, 1.1, 5.6),
U c(10.6, 22.5, 14.1, 11.5, 3.3, 12.9, 9, 3.9))
<- function(iA, U, C = diag(ncol(U)), V) {
wsm # 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
<- solve(C)
iC <- iA %*% U
iAU - iAU %*% solve(iC + crossprod(V, iAU), crossprod(V, iA))
iA
}
wsm(solve(A), U, V = U)
- Confirm that Woodbury’s formula gives the same result as directly inverting \(\mathbf{B}\) with
solve()
(or equivalent).
Solution
<- A + tcrossprod(U)
B <- solve(B)
iB all.equal(iB, wsm(solve(A), U, V = U))
- 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
<- cbind(
C c(3.9, -1.7),
c(-1.7, 2.9)
)
wsm(solve(A), U, C, U)
- Confirm that the Sherman-Morrison-Woodbury formula gives the same result as directly inverting \(\mathbf{D}\)
Solution
<- A + U %*% tcrossprod(C, U)
D all.equal(solve(D), wsm(solve(A), U, C, U))
Challenges III
- 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}^2} = 2 \boldsymbol{\Sigma}^{-1}. \]