VectorMatrixAddon         package:fMultivar         R Documentation

_V_e_c_t_o_r/_M_a_t_r_i_x _A_r_i_t_h_m_e_t_i_c_s _a_n_d _L_i_n_e_a_r _A_l_g_e_b_r_a

_D_e_s_c_r_i_p_t_i_o_n:

     A collection and description of functions  available for matrix
     arithmetics and linear  algebra.  

     Only functions which have been added by Rmetrics  are documented.
     These functions are often useful  for the manipulation of
     multivariate time series.  

     The origin of the functions is marked in the following way:

     R: part of R's base packages, 
      B: part of Rmetrics' fBasics package, 
      M: part of this Rmetrics package, fMultivar. 

     The functions are listed by topic. 

     General Matrix Functions:

       'matrix'     Creates a matrix from the given set of values,
       'diag'       R: Creates a diagonal matrix or extracts diagonals,
       'triang'     S: Extracs the lower tridiagonal part from a matrix,
       'Triang'     S: Extracs the upper tridiagonal part from a matrix,
       'pascal'     S: Creates a Pascal matrix,
       'colVec'     S: Creates a column vector from a vector,
       'rowVec'     S: Creates a row vector from a vector,
       'as.matrix'  R: Attempts to turn its argument into a matrix,
       'is.matrix'  R: Tests if its argument is a (strict) matrix,
       'dimnames'   R: Retrieves or sets the dimnames of an object,
       'colnames'   R: Retrieves or sets the column names,
       'colIds'     R: ... use alternatively,
       'rownames'   R: Retrieves or sets the row names,
       'rowIds'     R: ... use alternatively.

     Simple Matrix Operations:

       'dim'     R: Returns the dimension of a matrix object,
       'ncol'    R: Counts columns of a matrix object,
       'nrow'    R: Counts rows of a matrix object,
       'length'  R: Counts elements of a matrix object,
       '"["'     R: Subsets a matrix object,
       '"[["'    R: Subsets a matrix object,
       'cbind'   R: Augments a matrix object by columns,
       'rbind'   R: Augments a matrix object by rows.

     Basic Statistical Functions:

       'var'          R: Returns the variance matrix,
       'cov'          R: Returns the covariance matrix,
       'colStats'     B: Calculates column statistics,
       'rowStats'     B: Calculates row statistics,
       'colMeans'     B: Calculates column means,
       'rowMeans'     B: Calculates row means,
       'colAvgs'      B: Calculates column averages,
       'rowAvgs'      B: Calculates row averages,
       'colVars'      B: Calculates column variances,
       'rowVars'      B: Calculates row variances,
       'colStdevs'    B: Calculates column standard deviations,
       'rowStdevs'    B: Calculates row standard deviations,
       'colSkewness'  B: Calculates column skewness,
       'rowSkewness'  B: Calculates row skewness,
       'colKurtosis'  B: Calculates column kurtosis,
       'rowKurtosis'  B: Calculates row kurtosis,
       'colCumsums'   B: Calculates column cumulated sums,
       'rowCumsums'   B: Calculates row cumulated sums.

     Linear algebra:

       '%*%'          R: Returns the product of two matrices,
       '%x%', 'kron'  R: Returns the Kronecker product,
       'mexp'         M: Computes the exponential of a square matrix,
       'det'          R: Returns the determinante of a matrix,
       'inv'          S: Returns the inverse of a matrix,
       'norm'         S: Returns the norm of a matrix,
       'rk'           S: Returns the rank of a matrix,
       'tr'           S: Returns trace of a matrix,
       't'            R: Returns the transposed matrix,
       'vech'         M: Is the operator that stacks the lower triangle,
       'vec'          M: Is the operator that stacks a matrix.

     More linear algebra:

       'chol'          R: Returns the Cholesky factor matrix,
       'eigen'         R: Returns eigenvalues and eigenvectors,
       'svd'           R: Returns the singular value decomposition,
       'kappa'         R: Returns the condition number of a matrix,
       'qr'            R: Returns the QR decomposition of a matrix,
       'solve'         R: Solves a system of linear equations,
       'backsolve'     R: ... use when the matrix is upper triangular,
       'forwardsolve'  R: ... use when the matrix is lower triangular.

     Time Series Generation:

       'tslag'  R: Lagged or leading vector/matrix of selected order(s),
       'pdl'    R: Regressor matrix for polynomial distributed lags.

_U_s_a_g_e:

     triang(x)
     Triang(x)
     pascal(n)

     colVec(x)
     rowVec(x)
     colIds(x, ...)
     rowIds(x, ...)

     inv(x)
     norm(x, p = 2)
     rk(x, method = c("qr", "chol"))
     tr(x)
     kron(x, y)

     mexp(x, order = 8, method = c("pade", "taylor"))

     vec(x)
     vech(x)

     pdl(x, d = 2, q = 3, trim = FALSE)
     tslag(x, k = 1, trim = FALSE)

_A_r_g_u_m_e_n_t_s:

       d: [pdl] - 
           an integer specifying the order of the polynomial.  

       k: [tslag] - 
           an integer value, the number of positions the new series is 
          to lag or to lead the input series.  

  method: [mexp] - 
           Two methods are provided for the computation of the
          exponential of a matrix: Taylor series selected by
          'method="taylor"', and Pad\'e approximation with scaling and
          squaring to increase  precision selected by 'method="pade"'
          which is the default method. 
           [rk] - 
           a character value, the dimension of the square matrix. One
          can choose from two methods: For 'method = "qr"' the rank is
          computed as 'qr(x)\$rank', or alternatively for 
          'method="chol"' the rank is computed as 'attr(chol(x,
          pivot=TRUE), "rank")'. 

       n: [pascal] - 
           an integer value, the dimension of the square matrix. 

   order: [mexp] - 
           The order of approximation to be used, an integer value. The
          best value for this depends on machine precision. 

       p: [norm] - 
           an integer value, '1', '2' or 'Inf'. 'p=1' - The maximum
          absolute column sum norm which is defined  as the maximum of
          the sum of the absolute valued elements of columns  of the
          matrix. 'p=2' - The spectral norm is "the norm" of a matrix
          'X'.  This value is computed as the square root of the
          maximum eigenvalue  of 'CX' where 'C' is the conjugate
          transpose. 'p=Inf' - The maximum absolute row sum norm is
          defined  as the maximum of the sum of the absolute valued
          elements of rows of the matrix. 

       q: [pdl] - 
           an integer specifying the number of lags to use in  creating
          polynomial distributed lags. This must be  greater than d.  

    trim: [pdl] - 
           a logical flag; if TRUE, the missing values at  the
          beginning of the returned matrix will be trimmed.  
           [tslag] - 
           a logical flag, if TRUE, the missing values at the 
          beginning or end of the returned series will be trimmed.  The
          default value is FALSE.  

    x, y: a numeric matrix. 
           [tslag] - 
           a numeric vector, missing values are allowed. 
           [pdl] - 
           a numeric vector. 

     ...: arguments to be passed. 

_D_e_t_a_i_l_s:

     *Function from R's Base Package:* 

     Most of the functions are described in their R help pages which we
     recommend to consult for further information. For the additiotnal 
     functions added by Rmetrics we give a brief introduction. 

     *General Functions:* 

     Functions to generate matrices and related functions are described
      in the help page 'matrix'. To "decorate" these objects  several
     naming functions are available, a description can be found  on the
     help pages 'dimnames' and 'rownames'. 

     The function 'pascal' generates a Pascal matrix of order 'n' 
     which is a symmetric, positive, definite matrix with integer
     entries  made up from Pascal's triangle. The determinant of a
     Pascal matrix is 1.  The inverse of a Pascal matrix has integer
     entries. If 'lambda'  is an eigenvalue of a Pascal matrix, then
     '1/lambda' is also an  eigenvalue of the matrix. Pascal matrices
     are ill-conditioned.  

     The functions 'triang' and 'Triang' allow to transform a square
     matrix to a lower or upper triangular form.  A triangular matrix
     is either an upper triangular matrix or lower  triangular matrix.
     For the first case all matrix elements 'a[i,j]' of matrix 'A' are
     zero for 'i>j', whereas in the second case we have just the
     opposite situation. A lower triangular matrix is  sometimes also
     called left triangular. In fact, triangular matrices  are so
     useful that much computational linear algebra begins with 
     factoring or decomposing a general matrix or matrices into
     triangular  form. Some matrix factorization methods are the
     Cholesky factorization  and the LU-factorization. Even including
     the factorization step,  enough later operations are typically
     avoided to yield an overall  time savings. Triangular matrices
     have the following properties: the  inverse of a triangular matrix
     is a triangular matrix, the product of  two triangular matrices is
     a triangular matrix, the determinant of a  triangular matrix is
     the product of the diagonal elements, the  eigenvalues of a
     triangular matrix are the diagonal elements. 

     The functions 'colVec' and 'rowVec' transform a vector into  a
     column and row vector, respectively. A column vector is a matrix 
     object with one column, and a row vector is a matrix object with
     one row. 

     The functions 'dimnames', 'colname', and 'rowname' can be used to
     retrieve and set the names of matrices. The functions  'rowIds',
     'colIds', are S-Plus like synonyms. 

     *Simple Matrix Operations:* 

     The functions 'dim', 'nrow' and  'ncol' are functions to extract
     the dimension and   the number of rows and columns of a matrix.  

     The usual arithmetic operators, logical operators and mathematical
      functions like 'sqrt' or 'exp' and 'log' operate on matrices
     element by element. Note, that '"*"' is not matrix multiplication,
     instead we have to use '"%*%"'. 

     The methods '"["' and '"[["' are suited to extract subsets from a
     matrix, to delete rows and columns, or to permute rows and
     columns.  

     *Basic Statistical Functions:* 

     The functions 'var' and 'cov' compute the variance and covariance
     of a matrix. 

     Beside these functions 'colMeans' and 'rowMeans'  are R functions
     which compute the mean of columns and rows of a matrix.  Rmetrics
     has added further functions to compute column- or rowwise 
     variances, standard deviations, skewness, kurtosis and cumulated
     sums.  Two general  functions named 'rowStats' and  'colStats'
     allow to apply through the argument list any  function to compute
     row and column statistics from matrices. 

     'Linear Algebra:' 

     Matrix multiplication is done using the operator '%*%'.  

     The _Kronecker product_ can be computed using the operator  '%x%'
     or alternatively using the function 'kron'. 

     The function 'mexp' computes the exponential of a square  matrix
     x, defined as the sum from r=0 to infinity of x^r/r!. Two methods
     are provided Taylor series and Pad\'e approximation with scaling
     and squaring to increase precision. Also reported, as the
     '"accuracy"' attribute of the result, is an upper bound for the L2
     norm of the Cauchy error 'mexp(a, order + 10, method) - mexp(a,
     order, method)', and the used 'method' and 'order'. 

     The function 'det' computes the determinant of a matrix. 

     The inverse of a square matrix 'inv(X)' of dimension 'n' is
     defined so that  'X %*% inv(X) = inv(X) %*% X = diag(n)' where the
     matrix 'diag(n)' is the 'n'-dimensional identity matrix. A
     precondition for the existence of the matrix inverse is that the
     determinant of the matrix 'det(X)' is nonzero.  

     The function 't' computes the transposed of a square matrix. 

     The function 'vec' implements the operator that stacks a matrix as
     a column vector, to be more precise in a matrix with one column.
     vec(X) = (X11, X21, ..., XN1, X12, X22, ..., XNN). The function
     'vech' implements the operator that stacks the lower  triangle of
     a NxN matrix as an N(N+1)/2x1 vector: vech(X) =(X11, X21, X22,
     X31, ..., XNN), to be more precise in a  matrix with one row. 

     The function 'tr' computes the trace of a square matrix which is
     the sum of the diagonal elements of the matrix under
     consideration. 

     The function 'rk' computes the rank of a matrix which is  the
     dimension of the range of the matrix corresponding to the number 
     of linearly independent rows or columns of the matrix, or to the 
     number of nonzero singular values. The rank of a matrix is also
     named inear map.  

     The function 'norm' computes the norm of a matrix. Three choices 
     are possible:  'p=1' - The maximum absolute column sum norm which
     is defined  as the maximum of the sum of the absolute valued
     elements of columns  of the matrix. 'p=2' - The spectral norm is
     "the norm" of a matrix 'X'.  This value is computed as the square
     root of the maximum eigenvalue  of 'CX' where 'C' is the conjugate
     transpose. 'p=Inf' - The maximum absolute row sum norm is defined 
     as the maximum of the sum of the absolute valued elements of rows
     of the matrix. 

     'More Linear Algebra:' 

     The function 'chol' returns the Cholesky factor matrix, 'eigen'
     returns eigenvalues and eigenvectors, 'svd' returns the singular
     value decomposition, 'kappa' estimate the condition number of a
     matrix, 'qr' returns the QR decomposition of a matrix, 'ginv'
     returns the Moore-Penrose generalized inverse, 'solve' solves a
     system of linear equations, use 'backsolve' when the matrix is
     upper triangular, and 'forwardsolve' when the matrix is lower
     triangular. 

     'Time Series:' 
      The function 'pdl' returns a regressor matrix suitable  for
     polynomial distributed lags.  
      The function 'tslag' returns a lagged/led vector or matrix  for
     given time series data.

_A_u_t_h_o_r(_s):

     Marina Shapira and David Firth for the 'mexp' function, 
      Diethelm Wuertz for the Rmetrics R-port.

_R_e_f_e_r_e_n_c_e_s:

     Higham N.J., (2002); _Accuracy and Stability of Numerical
     Algorithms_,  2nd ed., SIAM.

     Golub, van Loan, (1996); _Matrix Computations_,  3rd edition.
     Johns Hopkins University Press. 

     Ward, R.C., (1977);   _Numerical computation of the matrix
     exponential with  accuracy estimate_, SIAM J. Num. Anal. 14,
     600-610.

     Moler C., van Loan C., (2003); _Nineteen dubious ways to compute
     the exponential of a matrix,  twenty-five years later_, SIAM
     Review 45, 3-49.

_E_x_a_m_p_l_e_s:

     ## SOURCE("fMultivar.6A-MatrixAddon")

     ## Create Pascal Matrix:
        P = pascal(3)
        P
        # Create lower triangle matrix
        L = triang(P)
        L
        # Extract diagonal part
        diag(P)
        
     ## Add/Subtract/Multiply/Divide:  
        X = P
        # Multiply matrix with a constant
        3 * X
        # Multiply two matrices elementwise
        X * P                     
        # Multiplies rows/columns of a matrix by a vector
        X %*% diag(P)            
        diag(P) %*% X           
            
     ## Operate on Subsets of a Matrix:
        n = 3; i = 2; j = 3
        D = diag(1:3)
        # Return the dimension of a matrix
        dim(P)                         
        # Get the last colum of a matrix
        P[, ncol(P)]                   
        # Delete a column of a matrix
        P[, -i]                      
        # Permute the columns of a matrix
        P[c(3, 1, 2), ]              
        # Augments matrix horizontally 
        cbind(P, D)                           
           
     ## Apply a function to all Elements of a Matrix: 
        # Return square root for each element
        sqrt(P)
        # Exponentiate the matrix elementwise
        exp(P)
        # Compute the median of each column
        apply(P, 2, "median") 
        # Test on all elements of a matrix       
        all( P > 2 )   
        # test on any element in a matrix                
        any( P > 2 )                  
          
     ## More Matrix Operations:
        # Return the product of two matrices
        P %*% D   
        # Return the Kronecker Product                     
        P %x% D                        
        # Return the transposed matrix
        t(P)                           
        # Return the inverse of a matrix
        inv(P)  
        # Return the norm of a matrix                      
        norm(P)                        
        # Return the determinante of a matrix
        det(P)                         
        # Return the rank of a matrix
        rk(P)                            
        # Return trace of a matrix
        tr(P)                          
        # Return the variance matrix
        var(P)     
        # Return the covariance matrix                   
        cov(P) 
        # Stack a matrix
        vec(P) 
        # Stack the lower triangle
        vech(P)
           
     ## Matrix Exponential:
        # Test case 1 from Ward (1977)
        test1 = t(matrix(c(
          4, 2, 0,
          1, 4, 1,
          1, 1, 4), 3, 3))
        mexp(test1)
        # Results on Power Mac G3 under Mac OS 10.2.8 
        #                    [,1]               [,2]               [,3]
        # [1,] 147.86662244637000 183.76513864636857  71.79703239999643
        # [2,] 127.78108552318250 183.76513864636877  91.88256932318409
        # [3,] 127.78108552318204 163.67960172318047 111.96810624637124
        # -- these agree with ward (1977, p608)
     ## Not run: 
        # A naive alternative to mexp, using spectral decomposition:
        mexp2 = function(matrix){
          z = eigen(matrix, sym = FALSE)
          Re(z$vectors %*% diag(exp(z$values)) %*% solve(z$vectors)) }
        mexp2(test1)
        # -- hopelessly inaccurate in all but the first column.
        # Test case 4 from Ward (1977)
        test4 <- structure(c(
          0, 0, 0, 0, 0, 0, 0, 0, 0, 1e-10,
          1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
          0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
          0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
          0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
          0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
          0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
          0, 0, 0, 0, 0, 0, 0, 0, 1, 0), .Dim = c(10, 10))
        attributes(mexp(test4))
        max(abs(mexp(test4) - mexp2(test4)))
        #[1] 8.746826694186494e-08
        # -- mexp2 is accurate to 7 d.p., whereas mexp to at least 14 d.p.              
     ## End(Not run)
        
     ## More Linear Algebra:
        X = P; b = c(1, 2, 3)
        # Return the Cholesky factor matrix
        chol(X)                        
        # Return eigenvalues and eigenvectors
        eigen(X)                       
        # Return the singular value decomposition
        svd(X)                         
        # Return the condition number of a matrix
        kappa(X)                       
        # Return the QR decomposition of a matrix
        qr(X)                          
        # Solve a system of linear equations
        # ... use backsolve when the matrix is upper triangular
        # ... use forwardsolve when the matrix is lower triangular
        solve(X, b)  
        backsolve(Triang(X), b)
        solve(Triang(X), b)                 
        forwardsolve(triang(X), b) 
        solve(triang(X), b)                                        

