Rev 4513 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
\documentclass{article}\usepackage{myVignette}\usepackage[authoryear,round]{natbib}\bibliographystyle{plainnat}%%\VignetteIndexEntry{Introduction to the Matrix Package}%%\VignetteDepends{Matrix}\SweaveOpts{engine=R,eps=FALSE,pdf=TRUE,width=5,height=3,strip.white=TRUE,keep.source=TRUE}\title{Introduction to the Matrix package --- as of Feb.~2005\footnote{There's an unfinished ``2nd Introduction to the Matrix package''which contains partly newer information, but is not at all self-contained.Eventually that will replace this one.}}\author{Douglas Bates\\R Core Development Group\\\email{bates@r-project.org}}\date{\today}\begin{document}\maketitle\begin{abstract}Linear algebra is at the core of many areas of statistical computingand from its inception the \Slang{} has supported numerical linearalgebra via a matrix data type and several functions and operators,such as \code{\%*\%}, \code{qr}, \code{chol}, and \code{solve}.However, these data types and functions do not provide direct accessto all of the facilities for efficient manipulation of densematrices, as provided by the Lapack subroutines, and they do notprovide for manipulation of sparse matrices.The \code{Matrix} package provides a set of S4 classes for dense andsparse matrices that extend the basic matrix data type. Methods fora wide variety of functions and operators applied to objects fromthese classes provide efficient access to BLAS (Basic Linear AlgebraSubroutines), Lapack (dense matrix), TAUCS (sparse matrix) andUMFPACK (sparse matrix) routines. One notable characteristic of thepackage is that whenever a matrix is factored, the factorization isstored as part of the original matrix so that further operations onthe matrix can reuse this factorization.\end{abstract}<<preliminaries, echo=FALSE>>=options(width=75)@\section{Introduction}\label{sec:Intro}Linear algebra is at the core of many statistical computing techniquesand, from its inception, the \Slang{} has supported numerical linearalgebra via a matrix data type and several functions and operators,such as \code{\%*\%}, \code{qr}, \code{chol}, and \code{solve}.Initially the numerical linear algebra functions in \RR{} calledunderlying Fortran routines from the Linpack~\citep{Linpack} andEispack~\cite{Eispack} libraries but over the years most of thesefunctions have been switched to use routines from theLapack~\cite{Lapack} library. Furthermore, \RR{} can be configured touse accelerated BLAS (Basic Linear Algebra Subroutines), such as thosefrom the Atlas~\cite{Atlas} project or Goto's BLAS~\cite{GotosBLAS}.Lapack provides routines for operating on several special forms ofmatrices, such as triangular matrices and symmetric matrices.Furthermore,matrix decompositions like the QR decompositions producemultiple output components that should be regarded as parts of asingle object. There is some support in R for operations on specialforms of matrices (e.g. the \code{backsolve}, \code{forwardsolve} and\code{chol2inv} functions) and for special structures (e.g. a QRstructure is implicitly defined as a list by the \code{qr},\code{qr.qy}, \code{qr.qty}, and related functions) but it is not asfully developed as it could be.Also there is no direct support for sparse matrices in R although\citet{koen:ng:2003} have developed a contributed package for sparsematrices based on SparseKit.The \code{Matrix} package provides S4 classes and methods for denseand sparse matrices. The methods for dense matrices use Lapack andBLAS. The sparse matrix methods use TAUCS~\citep{Taucs},UMFPACK~\citep{Umfpack}, and Metis~\citep{Metis}.\section{Classes for dense matrices}\label{sec:DenseClasses}The \code{Matrix} package will provide classes for real (stored asdouble precision) and complex (stored as double precision complex)dense matrices. At present only the real classes have beenimplemented. These classes are\begin{description}\item[dgeMatrix] Real matrices in general storage mode\item[dsyMatrix] Symmetric real matrices in non-packed storage\item[dspMatrix] Symmetric real matrices in packed storage (one triangle only)\item[dtrMatrix] Triangular real matrices in non-packed storage\item[dtpMatrix] Triangular real matrices in packed storage (triangle only)\item[dpoMatrix] Positive semi-definite symmetric real matrices innon-packed storage\item[dppMatrix] \ \ ditto \ \ in packed storage\end{description}Methods for these classes include coercion between these classes, whenappropriate, and coercion to the \code{matrix} class; methods formatrix multiplication (\code{\%*\%}); cross products(\code{crossprod}), matrix norm (\code{norm}); reciprocal conditionnumber (\code{rcond}); LU factorization (\code{lu}) or, for the\code{poMatrix} class, the Cholesky decomposition (\code{chol}); andsolutions of linear systems of equations (\code{solve}).Further, group methods have been defined for the \code{Arith} (basicarithmetic, including with scalar numbers) and the \code{Math} (basicmathematical functions) group..Whenever a factorization or a decomposition is calculated it ispreserved as a (list) element in the \code{factors} slot of theoriginal object. In this way a sequence of operations, such asdetermining the condition number of a matrix then solving a linearsystem based on the matrix, do not require multiple factorizations ofthe same matrix nor do they require the user to store the intermediateresults.\section{Classes for sparse matrices}\label{sec:SparseClasses}\subsection{Representations of sparse matrices}\label{ssec:SparseReps}Conceptually, the simplest representation of a sparse matrix is as atriplet of an integer vector \code{i} giving the row numbers, aninteger vector \code{j} giving the column numbers, and a numericvector \code{x} giving the non-zero values in the matrix. An S4 classdefinition might be\begin{Schunk}\begin{Sinput}setClass("dgTMatrix",representation(i = "integer", j = "integer", x = "numeric",Dim = "integer"))\end{Sinput}\end{Schunk}The triplet representation is row-oriented if elements in the same rowwere adjacent and column-oriented if elements in the same column wereadjacent. The compressed sparse row (csr) (or compressed sparsecolumn - csc) representation is similar to row-oriented triplet(column-oriented triplet) except that \code{i} (\code{j}) just storesthe index of the first element in the row (column). (There are acouple of other details but that is the gist of it.) These compressedrepresentations remove the redundant row (column) indices and providefaster access to a given location in the matrix because you only needto check one row (column).The preferred representation of sparse matrices in the SparseM packageis csr. Matlab uses csc. We hope that Octave will also use thisrepresentation. There are certain advantages to csc in systems like Rand Matlab where dense matrices are stored in column-major order. Forexample, Sivan Toledo's TAUCS~\cite{Taucs} library and Tim Davis'sUMFPACK~\cite{Umfpack} library are both based on csc and can both uselevel-3 BLAS in certain sparse matrix computations.The Matrix package provides the following classes for sparse matrices\begin{description}\item[dgTMatrix] general, numeric, sparse matrices in (a possiblyredundant) triplet form. This can be a convenient form in which toconstruct sparse matrices.\item[dgCMatrix] general, numeric, sparse matrices in the (sorted) compressedsparse column format.\item[dsCMatrix] symmetric, real, sparse matrices in the (sorted)compressed sparse column format. Only the upper or the lower triangle isstored. Although there is provision for both forms, the lowertriangle form works best with TAUCS.\item[dtCMatrix] triangular, real, sparse matrices in the (sorted)compressed sparse column format.\end{description}\bibliography{Matrix}\end{document}