spam (version 0.13-3)

ordering: Extract the permutation

Description

Extract the (inverse) permutation used by the Cholesky decomposition

Usage

ordering( x, inv=FALSE)

Arguments

x

object of class spam.chol.method returned by the function chol.

inv

Return the permutation (default) or inverse thereof.

Details

Recall that calculating a Cholesky factor from a sparse matrix consists of finding a permutation first, then calculating the factors of the permuted matrix. The ordering is important when working with the factors themselves. The ordering from a full/regular matrix is 1:n. Note that there exists many different algorithms to find orderings. See the examples, they speak more than 10 lines.

References

Ng, E. G. and B. W. Peyton (1993), "Block sparse Cholesky algorithms on advanced uniprocessor computers", SIAM J. Sci. Comput., 14, pp. 1034-1056.

See Also

chol, solve.

Examples

# Construct a pd matrix S to work with (size n)
n <- 100    # dimension
S <- .25^abs(outer(1:n,1:n,"-"))
S <- as.spam( S, eps=1e-4)
I <- diag(n)  # Identity matrix

cholS <- chol( S)
ord <- ordering(cholS)
iord <- ordering(cholS, inv=TRUE)

R <- as.spam( cholS ) # R'R = P S P', with P=I[ord,],
  # a permutation matrix (rows permuted).
RtR <- t(R) %*% R

# the following are equivalent:
as.spam( RtR -            S[ord,ord] )
as.spam( RtR[iord,iord] - S )
as.spam( t(R[,iord]) %*% R[,iord] - S )

# trivially:
as.spam( t(I[iord,]) - I[ord,])  # (P^-1)' = P  
as.spam( t(I[ord,]) - I[,ord])  # 
as.spam( I[iord,] - I[,ord])
as.spam( I[ord,]%*%S%*%I[,ord] - S[ord,ord] )
   # pre and post multiplication with P and P' is ordering