# R/eigen.r In shabbychef/madness: Automatic Differentiation of Multivariate Operations

# /usr/bin/r
#
# Author: Steven E. Pav
#
# This file is part of madness.
#
# madness is free software: you can redistribute it and/or modify
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# madness is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
#
# Created: 2016.01.09
# Copyright: Steven E. Pav, 2016
# Author: Steven E. Pav <[email protected]>

#' @include AllClass.r
#' @include utils.r
NULL

#' @title Spectral Decomposition of a Matrix
#'
#' @description
#'
#' Computes eigenvalues and eigenvectors of numeric (double, integer, logical) or
#'
#' @details
#'
#' The singular value decomposition of the matrix \eqn{X}{X} is
#' \deqn{X = U D V',}{X = U D V',}
#' where \eqn{U} and \eqn{V} are orthogonal, \eqn{V'} is \eqn{V}
#' transposed, and \eqn{D} is a diagonal matrix with the singular
#' values on the diagonal.
#'
#' @include AllClass.r
#' @inheritParams base::eigen
#' @param x \code{madness} object representing a numeric matrix
#' whose spectral decomposition is to be computed.
#' @return a list with components
#' \describe{
#' \item{values}{a \code{madness} object of a vector containing
#' the \eqn{p} eigenvalues of \code{x}, sorted in \emph{decreasing} order,
#' according to \code{Mod(value)} in the assymetric case when they might
#' be complex (even for real matrices). For real asymmetric matrices
#' the vector will be complex only if complex conjugate pairs of eigenvalues are
#' detected.}
#' \item{vectors}{either a \eqn{p \times p}{p * p} matrix whose columns contain the
#' eigenvectors of \code{x} or \code{NULL} if \code{only.values} is
#' \code{TRUE}. The vectors are normalized to unit length.
#'
#' Recall that the eigenvectors are only defined up to a constant:
#' even when the length is specified they are still only defined up to a
#' scalar of modulus one (the sign for real matrices).
#' If \code{r <- eigen(A)}, and \code{V <- r$vectors; lam <- r$values}, then
#' \deqn{A = V Lmbd V^{-1}}{A = V Lmbd V^(-1)}
#' (up to numerical fuzz), where \code{Lmbd =diag(lam)}.
#' }
#' }
#'
#' @references
#'
#' Izenman, Alan Julian. "Reduced-Rank Regression for the Multivariate Linear
#' Model." Journal of Multivariate Analysis 5, pp 248-264 (1975).
#' \url{http://www.sciencedirect.com/science/article/pii/0047259X75900421}
#'
#' Kato, Tosio. "Perturbation Theory for Linear Operators."
#' Springer (1995).
#' \url{http://www.maths.ed.ac.uk/~aar/papers/kato1.pdf}
#'
#' @name eigen
#' @template etc
NULL

## MM: do _NOT_  setGeneric() on existing functions! (==> conflict other pkg methods!)
#setGeneric('eigen', function(x,symmetric,only.values=FALSE,EISPACK=FALSE) standardGeneric('eigen'))

#' @name eigen
#' @rdname eigen
#' @exportMethod eigen
#' @rdname eigen
function(x,symmetric,only.values=FALSE,EISPACK=FALSE) {
if (!missing(symmetric) && symmetric) {
# if symmetric is set true, then we use only the lower
# triangular part? in which case we need to enforce symmetry...
vtag <- x@vtag
x <- vech(x,k=0)
x <- ivech(x,k=0,symmetric=TRUE)
x@vtag <- vtag
}

xtag <- x@xtag
if (missing(symmetric)) {
vvals <- eigen(x@val,only.values=only.values,EISPACK=EISPACK)
} else {
vvals <- eigen(x@val,symmetric=symmetric,only.values=only.values,EISPACK=EISPACK)
}

# eigen*values*
p <- length(vvals$values) rowr <- lapply(seq_len(p),function(jj) { (t(vvals$vectors[,jj,drop=FALSE]) %x% t(vvals$vectors[,jj,drop=FALSE])) %*% x@dvdx }) dvdx <- do.call(rbind,rowr) vtag <- paste0('eigen(',x@vtag,')$values')
varx <- x@varx
retv <- list()
retv$values <- new("madness", val=array(vvals$values,dim=c(1,p)), dvdx=dvdx, vtag=vtag, xtag=xtag, varx=varx)

# eigen*vectors* are trickier
rowr <- lapply(seq_len(p),
function(jj) {
scals <- 1 / (vvals$values[jj] - vvals$values[-jj])
krons <- t(vvals$vectors[,-jj,drop=FALSE]) %x% t(vvals$vectors[,jj,drop=FALSE]) %*% x@dvdx
rv <- vvals$vectors[,-jj,drop=FALSE] %*% diag(scals) %*% krons }) dvdx <- do.call(rbind,rowr) vtag <- paste0('eigen(',x@vtag,')$vectors')
varx <- x@varx
retv$vectors <- new("madness", val=vvals$vectors, dvdx=dvdx, vtag=vtag, xtag=xtag, varx=varx)
retv
})

#for vim modeline: (do not edit)
# vim:fdm=marker:fmr=FOLDUP,UNFOLD:cms=#%s:syn=r:ft=r

shabbychef/madness documentation built on May 29, 2019, 8:06 p.m.