Change default parameters

1. For lbfgs, the curvature could be very small
2. When iterate between A and Z, their line-search step sizes could be
difficult to identify. Now using a large maxiter to ensure descend and
adaptive init_stepsize to make it faster.

Author: junfengwen

RarmaSolvers.m

@@ -3,10 +3,10 @@
 % regularized ARMA
 
     properties(Constant)
-        SUCCESS = 0; 
-        ERROR_MAXITER = 10; 
-        ERROR_BACKTRACK = 20; 
-        ERROR_BACKTRACK_MAXITER = 21        
+        SUCCESS = 0;
+        ERROR_MAXITER = 10;
+        ERROR_BACKTRACK = 20;
+        ERROR_BACKTRACK_MAXITER = 21;
     end
 
     methods(Static)
@@ -16,7 +16,7 @@
         % x0 can be a matrix
 
             backtrack_backoff = 0.5;
-            backtrack_maxiter = 10;
+            backtrack_maxiter = 100; % make this large to ensure descend
 
             for iter = 1:backtrack_maxiter
 
@@ -38,8 +38,8 @@
 
             if ~isa(fun,'function_handle'), error('fmin_LBFGS -> improper function handle'); end
 
-            DEFAULTS.curvTol = 1e-3;            
-            DEFAULTS.funTol = 1e-6;   
+            DEFAULTS.curvTol = 1e-6;  % for curvature condition
+            DEFAULTS.funTol = 1e-6;
             DEFAULTS.m = 50;          % number gradients in bundle
             DEFAULTS.maxiter = 1000;   
             DEFAULTS.verbose = 0;
@@ -58,7 +58,7 @@
                 opts = RarmaUtilities.getOptions(opts, DEFAULTS);
             end
 
-            x = x0;   
+            x = x0;
             flag = RarmaSolvers.SUCCESS;
 
             t = length(x0);
@@ -70,7 +70,7 @@
             slope = Inf;
 
             % damped limited memory BFGS method
-            [f,g] = fun(x);  
+            [f,g] = fun(x);
             for iter = 1:opts.maxiter
                 % compute search direction
                 dir = RarmaSolvers.invhessmult(-g,Y,S,Rho,H0,inds,opts.m);
@@ -85,7 +85,7 @@
                 s = xnew - x;
                 y = gnew - g;
                 curvature = y'*s;
-                if curvature > opts.curvTol	
+                if curvature > opts.curvTol
                     rho = 1/curvature;
                     if length(inds) < opts.m
                         i = length(inds)+1;

demo.m

@@ -1,10 +1,11 @@
 clear
 clc
 
-rng(100);
+rng(10);
 num_reps = 10;
 Models = cell(num_reps,1);
 Xpredictall = cell(num_reps,1);
+isStable = zeros(num_reps,1);
 Err = zeros(num_reps,1);
 
 %% Generate data
@@ -14,15 +15,15 @@
 
 %% Learn RARMA
 % In practice, the following paramters should be 
-% cross-validated before applied
+% cross-validated for EACH dataset before applied
 opts = [];
 opts.ardim = 2;
 opts.madim = 2;
-opts.reg_wgt_ar = 1e-2;
-opts.reg_wgt_ma = 1e-1;
+opts.reg_wgt_ar = 0.07; % stronger regularization on A -> more stable
+opts.reg_wgt_ma = 0.01;
 for ii = 1:num_reps
     Models{ii} = rarma(Xtrainall{ii},opts);
-    [isStable, eigs] = RarmaUtilities.checkStable(Models{ii}.A)
+    isStable(ii) = RarmaUtilities.checkStable(Models{ii}.A);
 end
 
 %% Prediction and Evaluation

rarma.m

@@ -27,7 +27,7 @@
 DEFAULTS.ardim = 5;
 DEFAULTS.init_stepsize = 10;
 DEFAULTS.Loss = @RarmaFcns.euclidean_rarma; 
-DEFAULTS.maxiter = 100;
+DEFAULTS.maxiter = 1000;
 DEFAULTS.madim = 5;
 DEFAULTS.recover = 1;  % Recover B and Epsilon from learned Z
 DEFAULTS.reg_ar = @RarmaFcns.frob_norm_sq;
@@ -107,14 +107,14 @@
     Z = zeros(sizeZ);
     [A, obj, iter, msg] =opts.optimizer(@(Avec)(objA(Avec, X, Z)), Ainit(:));
     A = reshape(A, sizeA);
-    [Z, prev_obj] = iterateZ(Z,A);
+    [Z, prev_obj] = iterateZ(Z,A,opts.init_stepsize);
 
   for i = 1:opts.maxiter
     % Do A first since it returns the incorrect obj
-    A = iterateA(A,Z);
-    [Z, obj] = iterateZ(Z,A);
+    A = iterateA(A,Z,opts.init_stepsize/i); % adaptive stepsize
+    [Z, obj] = iterateZ(Z,A,opts.init_stepsize/i);
 
-    if abs(obj-prev_obj) < opts.TOL
+    if abs(prev_obj-obj) < opts.TOL % doing minimization
       break; 
     end
     prev_obj = obj;
@@ -133,15 +133,15 @@
   model.predict = @(Xstart, horizon, opts)(RarmaFcns.iterate_predict_ar(Xstart, model, horizon, opts));
 end
 
-function [A, f] = iterateA(A, Z)
+function [A, f] = iterateA(A, Z, init_stepsize)
   [f,g] = objA(A, X, Z);
-  stepsize = RarmaSolvers.line_search(A, f, g, @(A)(objA(A, X, Z)),opts.init_stepsize);
+  stepsize = RarmaSolvers.line_search(A, f, g, @(A)(objA(A, X, Z)),init_stepsize);
   A = A - stepsize*(g);
 end
 
-function [Z, f] = iterateZ(Z, A)
+function [Z, f] = iterateZ(Z, A, init_stepsize)
   [f,g] = objZ(Z, A);
-  stepsize = RarmaSolvers.line_search(Z, f, g, @(Z)(objZ(Z, A)), opts.init_stepsize/10); 
+  stepsize = RarmaSolvers.line_search(Z, f, g, @(Z)(objZ(Z, A)), init_stepsize); 
   Z = Z - stepsize*g;
 end