#' 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)
}