Nothing
#' Calculate the (truncated) higher-order SVD (HOSVD).
#'
#' Calculates the left singular vectors of each matrix unfolding of an array,
#' then calculates the core array. The resulting output is a Tucker
#' decomposition.
#'
#' If \code{r} is equal to the rank of \code{Y}, then \code{Y} is equal to
#' \code{atrans(S, U)}, up to numerical accuracy.
#'
#' More details on the HOSVD can be found in
#' \href{https://doi.org/10.1137/S0895479896305696}{ De Lathauwer
#' et. al. (2000)}.
#'
#' @param Y An array of numerics.
#' @param r A vector of integers. The rank of the truncated HOSVD.
#'
#' @return \code{U} A list of matrices with orthonormal columns. Each matrix
#' contains the mode-specific singular vectors of its mode.
#'
#' \code{S} An all-orthogonal array. This is the core array from the HOSVD.
#'
#' @export
#'
#' @author Peter Hoff.
#'
#' @references De Lathauwer, L., De Moor, B., & Vandewalle, J. (2000).
#' \href{http://doi.org/10.1137/S0895479896305696}{A
#' multilinear singular value decomposition}. \emph{SIAM journal on Matrix
#' Analysis and Applications}, 21(4), 1253-1278.
#'
#' @keywords decompositions
#'
#' @examples
#' #Generate random data.
#' p <- c(2, 3, 4)
#' X <- array(stats::rnorm(prod(p)), dim = p)
#'
#' #Calculate HOSVD.
#' hosvd_x <- hosvd(X)
#' S <- hosvd_x$S
#' U <- hosvd_x$U
#'
#' #Recover X.
#' trim(X - atrans(S, U))
#'
#' #S is all-orthogonal.
#' trim(mat(S, 1) %*% t(mat(S, 1)))
#' trim(mat(S, 2) %*% t(mat(S, 2)))
#' trim(mat(S, 3) %*% t(mat(S, 3)))
hosvd <- function(Y, r = NULL) {
## higher order svd
m <- dim(Y)
K <- length(m)
if (is.null(r)) {
r <- pmin(m, prod(m) / m)
}
## get standard tsvd
U <- list()
for (k in 1:K) {
U[[k]] <- svd(mat(Y, k))$u[, 1:r[k], drop = FALSE]
}
S <- atrans(Y, lapply(U, t))
## now put into identifiable form
if (all(r > 1)) {
G <- sign(S)
a <- list()
for (k in 1:K) {
a[[k]] <- rep(1, r[k])
if (k == 1) {
a[[k]][1] <- G[1]
}
idx <- matrix(1, nrow = r[k] - 1, ncol = K)
idx[, k] <- 2:r[k]
a[[k]][2:r[k]] <- a[[k]][1] * G[idx] / G[1]
}
aprod <- a[[1]]
for (k in 2:K) {
aprod <- outer(aprod, a[[k]])
}
U <- mapply(function(X, b) {
X %*% diag(b)
}, U, a, SIMPLIFY = FALSE)
S <- S * aprod
}
list(U = U, S = S)
}
#' Calculate the incredible higher-order LQ decomposition (HOLQ).
#'
#' This function calculates a generalization of the LQ decomposition to tensors.
#' This decomposition has a close connection to likelihood inference in
#' Kronecker structured covariande models.
#'
#' Given an array \code{X}, the default version of this function will calculate
#' (1) \code{L} a list of lower triangular matricies with positive diagonal
#' elements and unit determinant, \code{Z} an array of the same dimensions as
#' \code{X} that has special orthogonal properties, and (3) \code{sig} a numeric
#' such that \code{X} is the same as \code{sig * atrans(Z,L)} up to numeric
#' precision.
#'
#' This output (1) can be considered a generalization of the LQ decomposition to
#' tensors, (2) solves an optimization problem which the matrix LQ decomposition
#' solves, and (3) has a special connection to likelihood inference in the array
#' normal model.
#'
#' There are options to constrain the matrices in \code{L} to either be
#' diagonal, lower triangular with unit diagonal, or the identity matrix. Each
#' of these correspond to submodels in Kronecker structured covariance models.
#' The core array corresponding to each of these options has different
#' properities (see \href{https://doi.org/10.1016/j.laa.2016.04.033}{Gerard and Hoff
#' (2016)}). These more constrained tensor decompositions are called HOLQ
#' juniors.
#'
#' The MLE of the \eqn{i}th component covariance matrix under \emph{any}
#' elliptically contoured Kronecker structured covariance model is given by
#' \code{L[[i]] \%*\% t(L[[i]])}. The MLE for the total variation pamarameter
#' will be different depending on the distribution of the array, but for the
#' array normal it is \code{sig ^ 2 / prod(p)} (the "variance" form for the
#' total variation parameter).
#'
#' The likelihood ratio test statistic depends only on \code{sig} and can be
#' implemented in \code{\link{lrt_stat}}.
#'
#' The algorithm used to fit the HOLQ iteratively repeats the LQ decomposition
#' along each mode.
#'
#' For more details on the incredible HOLQ, see
#' \href{https://doi.org/10.1016/j.laa.2016.04.033}{Gerard and Hoff (2016)}.
#'
#' @param X An array of numerics.
#' @param tol A numeric. The maximum difference in frobenius norm between two
#' successive core arrays before stopping. Or maximum difference of the ratio
#' of sigs from 1 before stopping (which depends on the value of
#' \code{use_sig}).
#' @param itermax An integer. The maximum number of iterations of the LQ
#' decomposition to do before stopping.
#' @param mode_rep A vector of integers. The optional mode(s) to be considered
#' identity matrices.
#' @param mode_diag A vector of integers. The optional mode(s) to be considered
#' as independent but heteroscedastic.
#' @param mode_ldu A vector of integers. The optional modes(s) to be considered
#' as having unit diagonal.
#' @param print_diff A logical. Should the updates be printed to the screen each
#' iteration?
#' @param start_vals Determines how to start the optimization algorithm. Either
#' 'identity' (default), or 'residuals', which results in using the cholesky
#' square roots of the sample covariance matrices along each mode scaled to
#' have unit determinant. You can also use your own start values.
#' @param use_sig A logical. What stopping criterion should we use? Frobenius
#' norm of difference of cores (FALSE) or absolute value of difference of
#' ratio of \code{sig} from 1 (TRUE).
#'
#' @return \code{Z} The core array with scaled all-orthonormality property.
#'
#' \code{A} A list of matrices.
#'
#' \code{sig} A numeric. The total variation parameter. This is the "standard
#' deviation" form.
#'
#' @seealso \code{\link{array_bic_aic}} for using the output of \code{holq} to
#' calculate AIC and BIC,
#'
#' \code{\link{get_isvd}} for using the output of \code{holq} to calculate a
#' tensor generalization of the singular value decomposition.
#'
#' \code{\link{lq}} for the matrix LQ decomposition.
#'
#' \code{\link{lrt_stat}} for using the output of \code{holq} to calculate
#' likelihood ratio test statistics.
#'
#' \code{\link{mle_from_holq}} for using the output of \code{holq} to
#' calculate the maximum likelihood estimates of the component covariance
#' matrices under the array normal model.
#'
#' @export
#'
#' @references Gerard, D., & Hoff, P. (2016). A higher-order LQ
#' decomposition for separable covariance models.
#' \emph{Linear Algebra and its Applications}, 505, 57-84.
#' \url{https://doi.org/10.1016/j.laa.2016.04.033}
#' \url{http://arxiv.org/pdf/1410.1094v1.pdf}
#'
#' @author David Gerard.
#'
#' @keywords decompositions likelihood
#'
#' @examples
#' #Genrate random data.
#' p <- c(2, 3, 4)
#' X <- array(stats::rnorm(prod(p)), dim = p)
#'
#' #Calculate HOLQ with unit diagonal on 2nd mode,
#' # and diagonal along 3rd mode.
#' holq_x <- holq(X, mode_ldu = 2, mode_diag = 3)
#' Z <- holq_x$Z
#' A <- holq_x$A
#' sig <- holq_x$sig
#'
#' #Reconstruct X
#' trim(X - sig * atrans(Z, A), 10^-5)
#'
#' #Properties of core
#' #First mode has orthonormal rows.
#' trim(mat(Z, 1) %*% t(mat(Z, 1)), 10^-5)
#'
#' #Second mode has orthogonal rows.
#' trim(mat(Z, 2) %*% t(mat(Z, 2)), 10^-5)
#'
#' #Third mode has unit diagonal (up to scale).
#' diag(mat(Z, 3) %*% t(mat(Z, 3)))
holq <- function(X, tol = 10 ^ -9, itermax = 1000, mode_rep = NULL, mode_diag = NULL, mode_ldu = NULL, print_diff = TRUE, start_vals = "identity",
use_sig = TRUE) {
p <- dim(X) # the dimension of X
n <- length(p) # the number of modes (order)
if (is.null(mode_rep)) {
mode_rep <- n + 1
}
if (is.null(mode_diag)) {
mode_diag <- n + 1
}
if (is.null(mode_ldu)) {
mode_ldu <- n + 1
}
if (start_vals[1] == "identity") {
## start values are the identity
A <- list()
for (index in 1:n) {
A[[index]] <- diag(p[index])
}
X_fnorm <- fnorm(X)
Z_new <- X / X_fnorm
sig <- X_fnorm
} else if (start_vals[1] == "residuals") {
## start values are the empirical residuals
resid_X <- start_resids(X, mode_rep)
A <- resid_X$Sig
Z_new <- atrans(X, resid_X$Sig_inv)
Z_fnorm <- fnorm(Z_new)
Z_new <- Z_new / Z_fnorm
sig <- Z_fnorm
} else {
## if supplied own starting values
stopifnot(p == sapply(start_vals, nrow))
A_det1 <- list()
A_inv <- list()
for (k in 1:n) {
A_det1[[k]] <- start_vals[[k]] / det(start_vals[[k]]) ^ (1 / p[k])
A_inv[[k]] <- backsolve(A_det1[[k]], diag(p[k]), upper.tri = FALSE)
}
Z_new <- atrans(X, A_inv)
Z_fnorm <- fnorm(Z_new)
Z_new <- Z_new / Z_fnorm
sig <- Z_fnorm
A <- A_det1
}
crit_vec <- c(tol + 1, tol + 1)
iter_index <- 0
while (crit_vec[use_sig + 1] > tol & iter_index <= itermax) {
Z_old <- Z_new
sig_old <- sig
sig_old <- sig
for (mat_index in (1:n)[-mode_rep]) {
if (any(mat_index == mode_diag)) {
## assume unstructured along a mode multiply in sig first so that don't have to worry about it later
Z_mat_temp <- sig * diag(A[[mat_index]]) * mat(Z_new, mat_index)
L_k <- diag(sqrt(p[mat_index] * rowSums(Z_mat_temp ^ 2) / prod(p)))
Z_mat_temp <- diag(1 / diag(L_k)) %*% Z_mat_temp
Z_arr <- array(Z_mat_temp, dim = c(p[mat_index], p[-mat_index]))
Z_new <- aperm(Z_arr, match(1:n, c(mat_index, (1:n)[-mat_index])))
L_k_mult <- prod(diag(L_k)) ^ (1 / p[mat_index])
L_k <- L_k / L_k_mult
Z_new_fnorm <- fnorm(Z_new)
Z_new <- Z_new / Z_new_fnorm
sig <- L_k_mult * Z_new_fnorm
A[[mat_index]] <- L_k
} else if (any(mat_index == mode_ldu)) {
## assume unit lower triangular along a mode
lq_z <- lq(mat(Z_new, mat_index))
L_k_temp <- lq_z$L
L_k_diag <- diag(L_k_temp)
L_k <- L_k_temp %*% diag(1 / L_k_diag)
Z_arr <- array(diag(L_k_diag) %*% t(lq_z$Q), dim = c(p[mat_index], p[-mat_index]))
Z_new <- aperm(Z_arr, match(1:n, c(mat_index, (1:n)[-mat_index])))
Z_new_fnorm <- fnorm(Z_new)
Z_new <- Z_new / Z_new_fnorm
sig <- sig * Z_new_fnorm
A[[mat_index]] <- A[[mat_index]] %*% L_k
} else {
## assume diagonal along a mode
lq_z <- lq(mat(Z_new, mat_index))
L_k <- lq_z$L
Z_arr <- array(t(lq_z$Q), dim = c(p[mat_index], p[-mat_index]))
Z_new <- aperm(Z_arr, match(1:n, c(mat_index, (1:n)[-mat_index])))
L_k_mult <- prod(diag(L_k)) ^ (1 / p[mat_index])
L_k <- L_k / L_k_mult
Z_new_fnorm <- fnorm(Z_new)
Z_new <- Z_new / Z_new_fnorm
sig <- sig * L_k_mult * Z_new_fnorm
A[[mat_index]] <- A[[mat_index]] %*% L_k
}
}
z_diff <- fnorm(Z_old - Z_new)
sig_diff <- abs(sig_old / sig - 1)
crit_vec <- c(z_diff, sig_diff)
if (print_diff) {
cat("Scale Diff = ", sig_diff, "\n")
cat(" Scale = ", sig, "\n\n")
}
iter_index <- iter_index + 1
}
return(list(Z = Z_new, A = A, sig = sig))
}
#' The incredible higher-order polar decomposition (IHOP).
#'
#' Mmm, pancakes.
#'
#' This function will calculate the higher-order polar decomposition, a
#' generalization of the polar decomposition to tensors. It generalizes a
#' minimization formulation of the polar decomposition.
#'
#' Given an array \code{X}, \code{ihop} will output \code{L} a list of lower
#' triangular matrices with positive diagonal elements and unit Frobenius norm,
#' \code{R} a core array with certain orthogonality properties, and \code{sig} a
#' total variation parameter. We have that \code{X} is equal to \code{sig *
#' atrans(R, L)} up to numerical precision.
#'
#' \code{t(solve(L[[i]])) \%*\% mat(R, i)} will have orthonormal rows for all
#' \code{i}.
#'
#' For more details on the IHOP, see
#' \href{https://doi.org/10.1016/j.laa.2016.04.033}{Gerard and Hoff (2016)}.
#'
#' @param X An array of numerics.
#' @param itermax An integer. The maximum number of iterations to perform during
#' the optimization procedure.
#' @param tol A numeric. The algorithm will stop when the Frobenius norm of the
#' difference of core arrays between subsequent iterations is below \code{tol}
#' (for \code{use_sig = FALSE}) or when the absolute difference between the
#' ratio of subsequent values of \code{sig} and 1 is less than \code{tol} (for
#' \code{use_sig = TRUE}).
#' @param print_diff A logical. Should we print the updates of the algorithm?
#' @param mode_rep A vector. Which component matrices should be set to be the
#' identity?
#' @param use_sig A logical. See \code{tol}.
#' @return \code{R} A core array which, in combination with \code{L}, has
#' certain orthogonality properties.
#'
#' \code{L} A list of lower triangular matrices with unit Frobenius norm.
#'
#' \code{sig} A numeric.
#'
#' @export
#'
#' @keywords decompositions
#'
#' @references Gerard, D., & Hoff, P. (2016). A higher-order LQ
#' decomposition for separable covariance models.
#' \emph{Linear Algebra and its Applications}, 505, 57-84.
#' \url{https://doi.org/10.1016/j.laa.2016.04.033}
#' \url{http://arxiv.org/pdf/1410.1094v1.pdf}
#'
#' @author David Gerard.
#'
#' @examples
#' #Generate random data.
#' p <- c(2, 3, 4)
#' X <- array(stats::rnorm(prod(p)), dim = p)
#'
#' #Calculate IHOP.
#' ihop_x <- ihop(X)
#' R <- ihop_x$R
#' L <- ihop_x$L
#' sig <- ihop_x$sig
#'
#' #Reconstruct X
#' trim(X - sig * atrans(R, L))
#'
#' #Orthogonality properties
#' ortho_1 <- t(solve(L[[1]])) %*% mat(R, 1)
#' trim(ortho_1 %*% t(ortho_1))
#'
#' ortho_2 <- t(solve(L[[2]])) %*% mat(R, 2)
#' trim(ortho_2 %*% t(ortho_2))
#'
#' ortho_3 <- t(solve(L[[3]])) %*% mat(R, 3)
#' trim(ortho_3 %*% t(ortho_3))
ihop <- function(X, itermax = 100, tol = 10 ^ -9, print_diff = TRUE, mode_rep = NULL, use_sig = TRUE) {
p <- dim(X)
n <- length(p)
### set starting values to identity
A <- list()
for (index in 1:n) {
A[[index]] <- diag(p[index]) / sqrt(p[index])
}
X_fnorm <- fnorm(X)
Z_new <- X / X_fnorm
sig <- X_fnorm * prod(sqrt(p))
iter_index <- 1
if (is.null(mode_rep)) {
mode_rep <- n + 1
}
crit_vec <- c(tol + 1, tol + 1) ### use z_diff or sig_diff?
while (iter_index <= itermax & crit_vec[use_sig + 1] > tol) {
Z_old <- Z_new
sig_old <- sig
for (mode_index in 1:n) {
if (any(mode_index == mode_rep)) {
## Do nothing
} else {
polar_z <- polar(A[[mode_index]] %*% mat(Z_new, mode_index))
chol_p <- t(chol(polar_z$P))
Z_arr <- array(t(chol_p) %*% polar_z$Z, dim = c(p[mode_index], p[-mode_index]))
Z_new <- aperm(Z_arr, match(1:n, c(mode_index, (1:n)[-mode_index])))
L_k_mult <- fnorm(chol_p)
L_k <- chol_p / L_k_mult
Z_new_fnorm <- fnorm(Z_new)
Z_new <- Z_new / Z_new_fnorm
sig <- sig * L_k_mult * Z_new_fnorm
A[[mode_index]] <- L_k
iter_index <- iter_index + 1
}
}
z_diff <- fnorm(Z_old - Z_new)
sig_diff <- abs(sig_old / sig - 1)
crit_vec <- c(z_diff, sig_diff)
if (print_diff) {
cat("Scale Diff = ", sig_diff, "\n")
cat(" Scale = ", sig, "\n\n")
}
}
return(list(R = Z_new, L = A, sig = sig))
}
#' Calculate the incredible SVD (ISVD).
#'
#' The ISVD is a generalization of the SVD to tensors. It is derived from the
#' incredible HOLQ.
#'
#' Let \code{sig * atrans(Z, L)} be the HOLQ of \code{X}. Then the ISVD
#' calculates the SVD of each \code{L[[i]]}, call it \code{U[[i]] \%*\% D[[i]]
#' \%*\% t(W[[i]])}. It then returns \code{l = sig}, \code{U}, \code{D}, and
#' \code{V = atrans(Z, W)}. These values have the property that \code{X} is
#' equal to \code{l * atrans(atrans(V, D), U)}, up to numerical precision.
#' \code{V} is also scaled all-orthonormal.
#'
#' For more details on the ISVD, see
#' \href{https://doi.org/10.1016/j.laa.2016.04.033}{Gerard and Hoff (2016)}.
#'
#' @param x_holq The output from \code{\link{holq}}.
#'
#' @return l A numeric.
#' @return U A list of orthogonal matrices.
#' @return D A list of diagonal matrices with positive diagonal entries and unit
#' determinant. The diagonal entries are in descending order.
#' @return V A scaled all-orthonormal array.
#'
#' @export
#'
#' @author David Gerard.
#'
#' @keywords decompositions
#'
#' @references Gerard, D., & Hoff, P. (2016). A higher-order LQ
#' decomposition for separable covariance models.
#' \emph{Linear Algebra and its Applications}, 505, 57-84.
#' \url{https://doi.org/10.1016/j.laa.2016.04.033}
#' \url{http://arxiv.org/pdf/1410.1094v1.pdf}
#'
#' @examples
#' #Generate random data.
#' p <- c(4,4,4)
#' X <- array(stats::rnorm(prod(p)), dim = p)
#'
#' #Calculate HOLQ, then ISVD
#' holq_x <- holq(X)
#' isvd_x <- get_isvd(holq_x)
#' l <- isvd_x$l
#' U <- isvd_x$U
#' D <- isvd_x$D
#' V <- isvd_x$V
#'
#' #Recover X
#' trim(X - l * atrans(atrans(V, D), U))
#'
#' #V is scaled all-orthonormal
#' trim(mat(V, 1) %*% t(mat(V, 1)), epsilon = 10^-5)
#'
#' trim(mat(V, 2) %*% t(mat(V, 2)), epsilon = 10^-5)
#'
#' trim(mat(V, 3) %*% t(mat(V, 3)), epsilon = 10^-5)
get_isvd <- function(x_holq) {
## x_holq is the output of the holq
p <- dim(x_holq$Z)
n <- length(p)
A_svd <- lapply(x_holq$A, svd)
D <- vector(length = n, mode = "list")
U <- vector(length = n, mode = "list")
V_k <- vector(length = n, mode = "list")
for (mode_index in 1:n) {
U[[mode_index]] <- A_svd[[mode_index]]$u
D[[mode_index]] <- diag(A_svd[[mode_index]]$d)
V_k[[mode_index]] <- A_svd[[mode_index]]$v
}
V <- atrans(x_holq$Z, lapply(V_k, t))
l <- x_holq$sig
return(list(l = l, U = U, D = D, V = V))
}
#' Calculate the higher-order orthogonal iteration (HOOI).
#'
#' This function will calculate the best rank \code{r} (where \code{r} is a
#' vector) approximation (in terms of sum of squared differences) to a given
#' data array.
#'
#' Given an array \code{X}, this code will find a core array \code{G} and a list
#' of matrices with orthonormal columns \code{U} that minimizes \code{fnorm(X -
#' atrans(G, U))}. If \code{r} is equal to the dimension of \code{X}, then it
#' returns the HOSVD (see \code{\link{hosvd}}).
#'
#' For details on the HOOI see
#' \href{https://doi.org/10.1137/S0895479898346995}{Lathauwer et
#' al (2000)}.
#'
#' @param X An array of numerics.
#' @param r A vector of integers. This is the given low multilinear rank of the
#' approximation.
#' @param tol A numeric. Stopping criterion.
#' @param print_fnorm Should updates of the optimization procedure be printed?
#' This number should get larger during the optimizaton procedure.
#' @param itermax The maximum number of iterations to run the optimization procedure.
#'
#' @return \code{G} An all-orthogonal core array.
#'
#' \code{U} A vector of matrices with orthonormal columns.
#'
#' @export
#'
#' @author David Gerard.
#'
#' @keywords decompositions
#'
#' @references De Lathauwer, L., De Moor, B., & Vandewalle, J. (2000).
#' \href{http://doi.org/10.1137/S0895479898346995}{On the best
#' rank-1 and rank-(\eqn{r_1, r_2,..., r_n}) approximation of higher-order tensors}.
#' \emph{SIAM Journal on Matrix Analysis and Applications}, 21(4), 1324-1342.
#'
#' @examples
#' ## Generate random data.
#' p <- c(2, 3, 4)
#' X <- array(stats::rnorm(prod(p)), dim = p)
#'
#' ## Calculate HOOI
#' r <- c(2, 2, 2)
#' hooi_x <- hooi(X, r = r)
#' G <- hooi_x$G
#' U <- hooi_x$U
#'
#' ## Reconstruct the hooi approximation.
#' X_approx <- atrans(G, U)
#' fnorm(X - X_approx)
hooi <- function(X, r, tol = 10 ^ -6, print_fnorm = FALSE, itermax = 500) {
p <- dim(X)
n <- length(p)
## Check input -----------------------------------------------------
assertthat::are_equal(length(r), n)
assertthat::assert_that(all(p >= r))
assertthat::assert_that(all(r >= 0))
assertthat::assert_that(tol > 0)
assertthat::assert_that(itermax > 0)
assertthat::assert_that(is.logical(print_fnorm))
## Corner case of Rank = 0 -----------------------------------------
if (any(r == 0)) {
G <- array(0, dim = p)
U_new <- list()
for (mode_index in 1:n) {
U_new[[mode_index]] <- matrix(0, nrow = p[mode_index], ncol = p[mode_index])
}
return(list(G = G, U = U_new))
}
U_new <- hosvd(X, r = r)$U
epsilon <- tol + 1 ## how far U_old is from U_new
iter_index <- 1
fnorm_val <- -Inf ## the fnorm
while (epsilon > tol & iter_index < itermax) {
U_old <- U_new
## do svd along each mode
for (mode_index in 1:n) {
U_mult <- U_new
U_mult[[mode_index]] <- diag(p[mode_index])
Y <- atrans(X, lapply(U_mult, t))
U_new[[mode_index]] <- svd(mat(Y, mode_index), nu = r[mode_index], nv = 0)$u
}
## calculate how far the matrices are from one another
dist_current <- rep(NA, length = n)
for (mode_index in 1:n) {
dist_current[mode_index] <- fnorm(U_new[[mode_index]] - U_old[[mode_index]])
}
epsilon <- max(dist_current)
if (print_fnorm == TRUE) {
fnorm_old <- fnorm_val
fnorm_val <- fnorm(atrans(X, lapply(U_new, t)))
assertthat::assert_that(fnorm_val >= fnorm_old)
cat(fnorm_val, "\n")
}
iter_index <- iter_index + 1
}
G <- atrans(X, lapply(U_new, t))
return(list(G = G, U = U_new))
}
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.