R/compact_svd.R

In kschmaus/gammmbest: Method Of Moments Based Estimation for Generalized Additive Mixed Models

#' Compact Singular Value Decomposition of Design Matrix
#'
#' Singular value decomposition is performed, where only the first \code{r}
#' singular values and (left/right) singular vectors are retained, where
#' \code{r} is the matrix rank. The design matrix is assumed to be the
#' concatenation between fixed effects and random effects design matrices.
#'
#' @param x,x_list A matrix or list of matrices.
#' @param fixef_index An integer vector specifying which columns of \code{x}
#'        have fixed effects coefficients.
#' @param ranef_index An integer vector specifying which columns of \code{x}
#'        have random effects coefficients.
#' @return An object of class \code{gammmbest} containing
#'   \item{rank}{The rank of the design matrix.}
#'   \item{D}{A diagonal matrix of nonzero singular values.}
#'   \item{U}{A matrix whose columns contain left singular vectors.}
#'   \item{V}{A matrix whose columns contain right singular vectors.}
#'   \item{V1}{The right singular vectors associated with fixed effects}
#'   \item{V2}{The right singular vectors associated with random effects}
#'
#' @examples
#' x <- matrix(1:9, 3, 3)
#' xx <- cbind(x, x)
#' compact_svd(xx, fixef_index = 1:3, ranef_index = 4:6)
compact_svd <- function(x, fixef_index, ranef_index) {
	x <- x[, c(fixef_index, ranef_index), drop = FALSE]
	svd_x <- svd(x)

	rank <- Matrix::rankMatrix(x, sval = svd_x[["d"]])
	attributes(rank) <- NULL
	seq_len_rank <- seq_len(rank)

	D <- diag(svd_x[["d"]][seq_len_rank], rank, rank)
	U <- svd_x[["u"]][, seq_len_rank, drop = FALSE]
	V <- svd_x[["v"]][, seq_len_rank, drop = FALSE]
	V1 <- V[seq_along(fixef_index), , drop = FALSE]
	V2 <- V[seq_along(ranef_index) + length(fixef_index), , drop = FALSE]

    list(rank = rank, D = D, U = U, V = V, V1 = V1, V2 = V2)
}

#' @describeIn compact_svd Calls \code{compact_svd} for a list of matrices.
compact_svd_list <- function(x_list, fixef_index, ranef_index) {
    Map(compact_svd, x_list, list(fixef_index), list(ranef_index))
}