Nothing
context("Solve Quadratic Positive Definite")
test_that("can solve Ax = b", {
# http://papers.nips.cc/paper/5322-a-differential-equation-for-modeling-nesterovs-accelerated-gradient-method-theory-and-insights.pdf
# A is a 500×500 random
# positive definite matrix and b a random vector.
# The eigenvalues of A are between 0.001 and 1.
# The vector b is generated as i. i. d. Gaussian random variables with mean 0 and variance 25.
posdef_quad_fg <- function(n = 500, min_eig = 0.001, max_eig = 1,
bmean = 0, bvar = 25) {
b <- stats::rnorm(n = n, mean = bmean, sd = sqrt(bvar))
A <- rposdef(n = n, ev = stats::runif(n = n, min = min_eig, max = max_eig))
quad_matrix_fg(A, b)
}
# Generating a random positive-definite matrix with user-specified positive
# eigenvalues
# If eigenvalues are not specified, they are generated from a uniform
# distribution
# From https://stat.ethz.ch/pipermail/r-help/2008-February/153708.html
rposdef <- function(n, ev = stats::runif(n, 0, 10)) {
Z <- matrix(ncol = n, stats::rnorm(n^2))
decomp <- qr(Z)
Q <- qr.Q(decomp)
R <- qr.R(decomp)
d <- diag(R)
ph <- d / abs(d)
O <- Q %*% diag(ph)
t(O) %*% diag(ev) %*% O
}
# Returns function and gradient for the strongly convex function
# 0.5 xA'x - b't
# gradient is 0.5A'x + 0.5Ax - b
# http://math.stackexchange.com/questions/222894/how-to-take-the-gradient-of-the-quadratic-form
# Strongly convex
quad_matrix_fg <- function(A, b) {
f <- function(par) {
drop(0.5 * t(par) %*% A %*% par - t(b) %*% par)
}
f0 <- f(solve(A, b))
res <-
list(
fn = function(par) {
f(par) - f0
},
gr = function(par) {
as.vector(0.5 * (t(A) %*% par + A %*% par) - b)
},
hs = function(par) {
A
},
A = A,
b = b,
f0 = f
)
res
}
n <- 250
pd_fg <- posdef_quad_fg(
n = n, min_eig = 0.001, max_eig = 1,
bmean = 0, bvar = 25
)
pd_res <- mize(rep(0, n), fg = pd_fg, max_iter = 1000)
expect_equal(pd_res$par, solve(pd_fg$A, pd_fg$b), tol = 1e-5)
})
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.