Browse files

Add toeplitz.propack decomposition method

  • Loading branch information...
1 parent e31abb0 commit 280ad23af313b65953dd3a142aae851e97e857d2 @eodus committed Apr 26, 2012
Showing with 47 additions and 2 deletions.
  1. +47 −2 R/toeplitz.R
View
49 R/toeplitz.R
@@ -155,6 +155,51 @@ tmatmul <- function(tmat, v, transposed = FALSE) {
stop("'SVD' method is not applicable to toeplitz SSA");
}
-"decompose.toeplitz-ssa.propack" <- function(x, ...) {
- stop("'PROPACK' method is not applicable to toeplitz SSA");
+"decompose.toeplitz-ssa.propack" <- function(x,
+ neig = min(50, L, K),
+ ...,
+ force.continue = FALSE) {
+ N <- x$length; L <- x$window; K <- N - L + 1;
+
+ # Check, whether continuation of decomposition is requested
+ if (!force.continue && nlambda(x) > 0)
+ stop("Continuation of decompostion is not yet implemented for this method.");
+
+ h <- .get(x, "hmat", allow.null = TRUE);
+ if (is.null(h)) {
+ F <- .get(x, "F");
+ h <- new.hmat(F, L = L);
+ }
+
+ olambda <- .get(x, "olambda", allow.null = TRUE);
+ U <- .get(x, "U", allow.null = TRUE);
+
+ T <- .get(x, "tmat", allow.null = TRUE);
+ if (is.null(T)) {
+ T <- new.tmat(F, L = L);
+ }
+
+ S <- propack.svd(T, neig = neig, ...);
+
+ # Save results
+ .set(x, "hmat", h);
+ .set(x, "tmat", T);
+ .set(x, "olambda", S$d);
+ if (!is.null(S$u))
+ .set(x, "U", S$u);
+
+ num <- length(S$d);
+ lambda <- numeric(num);
+ V <- matrix(nrow = K, ncol = num);
+ for (i in 1:num) {
+ Z <- hmatmul(h, S$u[, i], transposed = TRUE);
+ lambda[i] <- sqrt(sum(Z^2));
+ V[, i] <- Z / lambda[i];
+ }
+
+ # Save results
+ .set(x, "lambda", lambda);
+ .set(x, "V", V);
+
+ x;
}

0 comments on commit 280ad23

Please sign in to comment.