40 lines
1.3 KiB
R
40 lines
1.3 KiB
R
#' Squared Frobenius norm of difference of Kronecker product matrices
|
|
#'
|
|
#' \|(A1 %x% ... %x% Ar - B1 %x% ... %x% Br\|_F
|
|
#'
|
|
#' This is equivalent to the expression
|
|
#' \code{norm(Reduce(kronecker, As) - Reduce(kronecker, Bs), "F")}, but faster.
|
|
#'
|
|
#' @examples
|
|
#' A1 <- matrix(rnorm(5^2), 5)
|
|
#' A2 <- matrix(rnorm(7^2), 7)
|
|
#' B1 <- matrix(rnorm(5^2), 5)
|
|
#' B2 <- matrix(rnorm(7^2), 7)
|
|
#' stopifnot(all.equal(
|
|
#' dist.kron.norm(list(A1, A2), list(B1, B2)),
|
|
#' norm(kronecker(A1, A2) - kronecker(B1, B2), "F")
|
|
#' ))
|
|
#'
|
|
#' p <- c(3, 7, 5, 2)
|
|
#' As <- Map(function(pj) matrix(rnorm(pj^2), pj), p)
|
|
#' Bs <- Map(function(pj) matrix(rnorm(pj^2), pj), p)
|
|
#' stopifnot(all.equal(
|
|
#' dist.kron.norm(As, Bs),
|
|
#' norm(Reduce(kronecker, As) - Reduce(kronecker, Bs), "F")
|
|
#' ))
|
|
#'
|
|
#' @export
|
|
dist.kron.norm <- function(As, Bs, eps = .Machine$double.eps) {
|
|
if (is.list(As) && is.list(Bs)) {
|
|
norm2 <- prod(unlist(Map(function(x) sum(x^2), As))) -
|
|
2 * prod(unlist(Map(function(a, b) sum(a * b), As, Bs))) +
|
|
prod(unlist(Map(function(x) sum(x^2), Bs)))
|
|
} else if (is.matrix(As) && is.matrix(Bs)) {
|
|
norm2 <- sum((As - Bs)^2)
|
|
} else {
|
|
stop("Unexpected input")
|
|
}
|
|
|
|
if (abs(norm2) < .Machine$double.eps) 0 else sqrt(norm2)
|
|
}
|