Even faster code... but potentially less precise

config-i1 · Oct 24, 2024 · ee76e1e · ee76e1e
1 parent c28d0dd
commit ee76e1e
Show file tree

Hide file tree

Showing 11 changed files with 80 additions and 29 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: legion
 Type: Package
 Title: Forecasting Using Multivariate Models
-Version: 0.2.1.41002
-Date: 2024-10-22
+Version: 0.2.1.41003
+Date: 2024-10-24
 Authors@R: c(person("Ivan", "Svetunkov", email = "ivan@svetunkov.ru", role = c("aut", "cre"),
                   comment="Lecturer at Centre for Marketing Analytics and Forecasting, Lancaster University, UK"),
              person("Kandrika Fadhlan", "Pritularga", email = "k.pritularga@lancaster.ac.uk", role = c("aut"),

diff --git a/NEWS b/NEWS
@@ -1,8 +1,9 @@
-legion v0.2.1 (Release data: 2024-10-22)
+legion v0.2.1 (Release data: 2024-10-24)
 =======
 
 Changes:
 * We now use Matrix package for sparse matrix computations. Very useful for data with lots of series.
+* We now calculate eigen values inside the C++ using the eigens_gen() function from the Armadillo library. This works faster but might produce less accurate estimates of parameters (need to experiment).
 
 Bugfixes:
 * Corrected the author details in the documentation.

diff --git a/R/RcppExports.R b/R/RcppExports.R
@@ -5,6 +5,10 @@ vSimulatorWrap <- function(arrayStates, arrayErrors, arrayF, arrayW, arrayG, mod
     .Call('_legion_vSimulatorWrap', PACKAGE = 'legion', arrayStates, arrayErrors, arrayF, arrayW, arrayG, modelLags)
 }
 
+discounter <- function(matrixF, matrixW, matrixG, k) {
+    .Call('_legion_discounter', PACKAGE = 'legion', matrixF, matrixW, matrixG, k)
+}
+
 vFitterWrap <- function(matrixY, matrixV, matrixF, matrixW, matrixG, lags, E, T, S, matrixO) {
     .Call('_legion_vFitterWrap', PACKAGE = 'legion', matrixY, matrixV, matrixF, matrixW, matrixG, lags, E, T, S, matrixO)
 }

diff --git a/R/ves.R b/R/ves.R
@@ -282,7 +282,13 @@ ves <- function(data, model="PPP", lags=c(frequency(data)),
         # Check the bounds
         # Edit KFP: change symmetric to FALSE
         if(bounds=="a"){
-            eigenValues <- eigen(elements$matF - elements$matG %*% elements$matW, only.values=TRUE, symmetric=FALSE)$values;
+            # eigenValues <- eigen(elements$matF - elements$matG %*% elements$matW, only.values=TRUE, symmetric=FALSE)$values;
+            if(nComponentsAll>10){
+                eigenValues <- discounter(elements$matF, elements$matW, elements$matG, 5);
+            }
+            else{
+                eigenValues <- eigen(elements$matF - elements$matG %*% elements$matW, only.values=TRUE)$values;
+            }
             if(max(abs(eigenValues)>(1 + 1E-50))){
                 return(max(abs(eigenValues))*1E+100);
             }
@@ -621,10 +627,10 @@ ves <- function(data, model="PPP", lags=c(frequency(data)),
 
         ### Persistence matrix
         matG <- switch(seasonalType,
-                       "i" =  matrix(0,nSeries*nComponentsAll,nSeries),
-                       "c" = matrix(0,nSeries*nComponentsNonSeasonal+1,nSeries));
+                       "i" =  Matrix(0, nSeries*nComponentsAll, nSeries, sparse=TRUE),
+                       "c" = Matrix(0, nSeries*nComponentsNonSeasonal+1, nSeries, sparse=TRUE));
         if(!persistenceEstimate){
-            matG <- persistenceValue;
+            matG <- Matrix(persistenceValue, sparse=TRUE);
         }
 
         ### Damping parameter
@@ -665,17 +671,17 @@ ves <- function(data, model="PPP", lags=c(frequency(data)),
                                 setdiff(c(1:nSeries*nComponentsAll),c(1:nComponentsAll)+nComponentsAll*(i-1))] <- 0.1;
             }
         }
-        matF <- switch(seasonalType,
-                       "i"=transitionValue,
-                       "c"=rbind(cbind(transitionValue[-(c(1:nSeries)*nComponentsAll),
-                                                       -(c(1:nSeries)*nComponentsAll)],
-                                       0),
-                                 c(transitionValue[nComponentsAll*nSeries,
-                                                   -(c(1:nSeries)*nComponentsAll)],1)));
+        matF <- Matrix(switch(seasonalType,
+                              "i"=transitionValue,
+                              "c"=rbind(cbind(transitionValue[-(c(1:nSeries)*nComponentsAll),
+                                                              -(c(1:nSeries)*nComponentsAll)],
+                                              0),
+                                        c(transitionValue[nComponentsAll*nSeries,
+                                                          -(c(1:nSeries)*nComponentsAll)],1))), sparse=TRUE);
 
         ### Measurement matrix
         if(seasonalType=="i"){
-            matW <- matrix(0,nSeries,nSeries*nComponentsAll);
+            matW <- Matrix(0, nSeries, nSeries*nComponentsAll, sparse=TRUE);
             for(i in 1:nSeries){
                 matW[i,c(1:nComponentsAll)+nComponentsAll*(i-1)] <- 1;
             }
@@ -686,7 +692,7 @@ ves <- function(data, model="PPP", lags=c(frequency(data)),
             }
         }
         else{
-            matW <- matrix(0,nSeries,nSeries*nComponentsNonSeasonal+1);
+            matW <- Matrix(0, nSeries, nSeries*nComponentsNonSeasonal+1, sparse=TRUE);
             for(i in 1:nSeries){
                 matW[i,c(1:nComponentsNonSeasonal)+nComponentsNonSeasonal*(i-1)] <- 1;
             }
@@ -1669,8 +1675,8 @@ ves <- function(data, model="PPP", lags=c(frequency(data)),
 
     ##### Return values #####
     model <- list(model=modelname,timeElapsed=Sys.time()-startTime,
-                  states=NA,persistence=persistenceValue,transition=transitionValue,
-                  measurement=matW, phi=dampedValue, B=B,
+                  states=NA,persistence=as.matrix(persistenceValue),transition=as.matrix(transitionValue),
+                  measurement=as.matrix(matW), phi=dampedValue, B=B,
                   lagsAll=lagsModel,
                   initialType=initialType,initial=initialValue,initialSeason=initialSeasonValue,
                   nParam=parametersNumber, imodel=ovesModel,

diff --git a/R/vets.R b/R/vets.R
@@ -337,7 +337,7 @@ vets <- function(data, model="PPP", lags=c(frequency(data)),
         else{
             matGBlockGamma <- NULL;
         }
-        matG <- rbind(matGBlockAlpha, matGBlockBeta, matGBlockGamma);
+        matG <- Matrix(rbind(matGBlockAlpha, matGBlockBeta, matGBlockGamma), sparse=TRUE);
 
         if(damped){
             dampedValue <- 0.95;
@@ -440,9 +440,9 @@ vets <- function(data, model="PPP", lags=c(frequency(data)),
             matFBlock23 <- NULL;
             matFBlock32 <- NULL;
         }
-        matF <- rbind(cbind(matFBlock11,matFBlock12,matFBlock13,deparse.level=0),
-                      cbind(matFBlock21,matFBlock22,matFBlock23,deparse.level=0),
-                      cbind(matFBlock31,matFBlock32,matFBlock33,deparse.level=0));
+        matF <- Matrix(rbind(cbind(matFBlock11,matFBlock12,matFBlock13,deparse.level=0),
+                             cbind(matFBlock21,matFBlock22,matFBlock23,deparse.level=0),
+                             cbind(matFBlock31,matFBlock32,matFBlock33,deparse.level=0)), sparse=TRUE);
 
         #### Blocks for measurement matrix ####
         if(componentsCommonLevel){
@@ -473,7 +473,7 @@ vets <- function(data, model="PPP", lags=c(frequency(data)),
         else{
             matWBlock3 <- NULL;
         }
-        matW <- cbind(matWBlock1,matWBlock2,matWBlock3,deparse.level=0);
+        matW <- Matrix(cbind(matWBlock1,matWBlock2,matWBlock3,deparse.level=0), sparse=TRUE);
 
         #### Vector of states ####
         matVt <- matrix(NA, nComponentsAll, obsStates);
@@ -791,16 +791,28 @@ vets <- function(data, model="PPP", lags=c(frequency(data)),
 
         # Check the bounds
         if(bounds=="a"){
-            eigenValues <- eigen(elements$matF - elements$matG %*% elements$matW, only.values=TRUE)$values;
+            # eigenValues <- eigen(elements$matF - elements$matG %*% elements$matW, only.values=TRUE)$values;
+            # eigenValues <- eigen(discounter(Matrix(elements$matF, sparse=TRUE),
+            #                                 Matrix(elements$matW, sparse=TRUE),
+            #                                 Matrix(elements$matG, sparse=TRUE)),
+            #                      only.values=TRUE)$values;
+            # eigenValues <- eigen(discounter(elements$matF, elements$matW, elements$matG),
+            #                      only.values=TRUE)$values;
+            if(nComponentsAll>10){
+                eigenValues <- discounter(elements$matF, elements$matW, elements$matG, 5);
+            }
+            else{
+                eigenValues <- eigen(elements$matF - elements$matG %*% elements$matW, only.values=TRUE)$values;
+            }
             if(max(abs(eigenValues)>(1 + 1E-50))){
                 return(max(abs(eigenValues))*1E+100);
             }
         }
 
         # Fit the model
         fitting <- vFitterWrap(switch(Etype, "M"=log(yInSample), yInSample),
-                               elements$matVt, Matrix(elements$matF, sparse=TRUE),
-                               Matrix(elements$matW, sparse=TRUE), Matrix(elements$matG, sparse=TRUE),
+                               elements$matVt, elements$matF,
+                               elements$matW, elements$matG,
                                lagsModel, Etype, Ttype, Stype, ot);
 
         # Calculate the loss
@@ -1462,7 +1474,7 @@ vets <- function(data, model="PPP", lags=c(frequency(data)),
 
     ##### Return values #####
     model <- list(model=modelName,timeElapsed=Sys.time()-startTime,
-                  states=NA, persistence=matG, transition=matF, measurement=matW, B=B,
+                  states=NA, persistence=as.matrix(matG), transition=as.matrix(matF), measurement=as.matrix(matW), B=B,
                   lagsAll=lagsModel,
                   # initialType=initialType,initial=initialValue,initialSeason=initialSeasonValue,
                   nParam=parametersNumber, occurrence=ovesModel,

diff --git a/R/vssFunctions.R b/R/vssFunctions.R
@@ -1193,7 +1193,7 @@ vssFitter <- function(...){
     ParentEnvironment <- ellipsis[['ParentEnvironment']];
 
     fitting <- vFitterWrap(switch(Etype, "M"=log(yInSample), yInSample),
-                           matVt, Matrix(matF, sparse=TRUE), Matrix(matW, sparse=TRUE), Matrix(matG, sparse=TRUE),
+                           matVt, matF, matW, matG,
                            lagsModel, Etype, Ttype, Stype, ot);
     matVt[] <- fitting$matVt;
     yFitted[] <- fitting$yfit;
@@ -1250,7 +1250,7 @@ vssForecaster <- function(...){
     if(h>0){
         yForecast[] <- vForecasterWrap(matrix(matVt[,(obsInSample+1):(obsInSample+lagsModelMax)],
                                               ncol=lagsModelMax),
-                                     Matrix(matF, sparse=TRUE), Matrix(matW, sparse=TRUE),
+                                     matF, matW,
                                      nSeries, h, Etype, Ttype, Stype, lagsModel);
     }
     else{

diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
@@ -27,6 +27,20 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// discounter
+arma::cx_vec discounter(arma::sp_mat const& matrixF, arma::sp_mat& matrixW, arma::sp_mat const& matrixG, int const& k);
+RcppExport SEXP _legion_discounter(SEXP matrixFSEXP, SEXP matrixWSEXP, SEXP matrixGSEXP, SEXP kSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< arma::sp_mat const& >::type matrixF(matrixFSEXP);
+    Rcpp::traits::input_parameter< arma::sp_mat& >::type matrixW(matrixWSEXP);
+    Rcpp::traits::input_parameter< arma::sp_mat const& >::type matrixG(matrixGSEXP);
+    Rcpp::traits::input_parameter< int const& >::type k(kSEXP);
+    rcpp_result_gen = Rcpp::wrap(discounter(matrixF, matrixW, matrixG, k));
+    return rcpp_result_gen;
+END_RCPP
+}
 // vFitterWrap
 RcppExport SEXP vFitterWrap(arma::mat const& matrixY, arma::mat matrixV, arma::sp_mat& matrixF, arma::sp_mat& matrixW, arma::sp_mat& matrixG, arma::uvec& lags, char const& E, char const& T, char const& S, arma::sp_mat& matrixO);
 RcppExport SEXP _legion_vFitterWrap(SEXP matrixYSEXP, SEXP matrixVSEXP, SEXP matrixFSEXP, SEXP matrixWSEXP, SEXP matrixGSEXP, SEXP lagsSEXP, SEXP ESEXP, SEXP TSEXP, SEXP SSEXP, SEXP matrixOSEXP) {
@@ -69,6 +83,7 @@ END_RCPP
 
 static const R_CallMethodDef CallEntries[] = {
     {"_legion_vSimulatorWrap", (DL_FUNC) &_legion_vSimulatorWrap, 6},
+    {"_legion_discounter", (DL_FUNC) &_legion_discounter, 4},
     {"_legion_vFitterWrap", (DL_FUNC) &_legion_vFitterWrap, 10},
     {"_legion_vForecasterWrap", (DL_FUNC) &_legion_vForecasterWrap, 9},
     {NULL, NULL, 0}

diff --git a/src/RcppExports.o b/src/RcppExports.o
diff --git a/src/legion.so b/src/legion.so
diff --git a/src/vssGeneral.cpp b/src/vssGeneral.cpp
@@ -33,6 +33,19 @@ arma::vec vErrorValue(arma::vec const &vectorY, arma::vec const &vectorYFit, cha
     return returnedValue;
 }
 
+/* # Wrapper for fitter */
+// [[Rcpp::export]]
+arma::cx_vec discounter(arma::sp_mat const &matrixF, arma::sp_mat &matrixW, arma::sp_mat const &matrixG, int const &k){
+    arma::cx_vec eigval;
+    // If eigen decomposition works, return the values
+    // Otherwise return a large number
+    if(!arma::eigs_gen(eigval, matrixF - matrixG * matrixW, k)){
+        eigval.fill(1e+300);
+    };
+    return eigval;
+    // return arma::conv_to<arma::mat>::from(matrixF - matrixG * matrixW);
+}
+
 // Fitter for vector models
 List vFitter(arma::mat const &matrixY, arma::mat &matrixV, arma::sp_mat const &matrixF,
              arma::sp_mat &matrixW, arma::sp_mat const &matrixG,

diff --git a/src/vssGeneral.o b/src/vssGeneral.o