Changeset 9a6c6f7 in tol


Ignore:
Timestamp:
Mar 9, 2009, 8:58:41 AM (18 years ago)
Author:
Víctor de Buen <vdebuen@…>
Branches:
1941-docker-image-to-test, devel, master
Children:
c6ede5b
Parents:
7168a54
Message:

Drop close to zero cells of all VMatrix to force maximum performance in large periodic ARMA process

git-svn-id: ​https://www.tol-project.org/svn/tolp/trunk@707 be4846ae-f051-483a-831f-b64e7d08527f

File:
1 edited

Legend:

Unmodified
Added
Removed

  • tol/stdlib/math/stat/models/bayesian/arima/_arma_process.tol

    r7168a54 r9a6c6f7  
    1010[[
    1111  Text _.autodoc.description = "Methods related with ARMA process";
     12
     13  Real drop = 1.E-6;
    1214
    1315  ////////////////////////////////////////////////////////////////////////////
    … …  
    5052
    5153  ////////////////////////////////////////////////////////////////////////////
     54  VMatrix DoDrop(VMatrix X)
     55  ////////////////////////////////////////////////////////////////////////////
     56  {
     57    Drop(X,drop)
     58  };
     59
     60  ////////////////////////////////////////////////////////////////////////////
     61  VMatrix CholeskiFactorCov(VMatrix cov)
     62  ////////////////////////////////////////////////////////////////////////////
     63  {
     64  //WriteLn("TRACE CholeskiFactorCov("+FullName(cov)+") Rows(cov)="<<VRows(cov)+" Sparsity(cov)="<<Real(VNonNullCells(cov)/VRows(cov)^2));
     65    VMatrix c = ARMAProcess::DoDrop(cov);
     66  //WriteLn("TRACE CholeskiFactorCov Sparsity(c)="<<Real(VNonNullCells(c)/VRows(c)^2));
     67  //Real t0 = Copy(Time);
     68  //VMatrix L = CholeskiFactor(c,"X",True);
     69  //Real tm = Time-t0;
     70  //WriteLn("TRACE CholeskiFactorCov time="<<tm);
     71  //L
     72    Pack(CholeskiFactor(c,"X",True))
     73  };
     74
     75  ////////////////////////////////////////////////////////////////////////////
    5276    VMatrix AvoidSymRoundErr(VMatrix S)
    5377  ////////////////////////////////////////////////////////////////////////////
    … …  
    6185  {
    6286    Real n = VRows(S);
    63     AvoidSymRoundErr(Case(
     87    ARMAProcess::DoDrop(ARMAProcess::AvoidSymRoundErr(Case(
    6488      m==n, S,
    6589      m<n,  Sub(S,1,1,m,m),
    … …  
    76100        )+Pol2VMat(pol, m, m);
    77101        aux
    78       }))
     102      })))
    79103  };
    80104
    … …  
    84108  {
    85109    Real n = VRows(L);
    86     Case(
     110    ARMAProcess::DoDrop(Case(
    87111      m==n, L,
    88112      m<n,  Sub(L,1,1,m,m),
    … …  
    95119        )+Pol2VMat(pol, m, m);
    96120        aux
    97       })
     121      }))
    98122  };
    99123
    … …  
    101125  NameBlock FastCholeskiCovFactor(Polyn ar_, Polyn ma_, Real m_)
    102126  ////////////////////////////////////////////////////////////////////////////
    103   {[[
     127  {
     128  //WriteLn("TRACE FastCholeskiCovFactor Begin "+Time);
     129    [[
    104130    Text _.autodoc.description =
    105131    "Computes a Fast Choleski decomposition of the covariance S=L*L' of an "
    … …  
    136162    //Stores arguments as members
    137163    Polyn _.ar = ar_;
    138   //Real _trc.1 = { WriteLn("TRACE _.ar=" << _.ar); 0 };
     164  //Real _trc.1 = { WriteLn("TRACE FastCholeskiCovFactor _.ar=" << _.ar); 0 };
    139165    Polyn _.ma = ma_;
    140   //Real _trc.2 = { WriteLn("TRACE _.ma=" << _.ma); 0 };
     166  //Real _trc.2 = { WriteLn("TRACE FastCholeskiCovFactor _.ma=" << _.ma); 0 };
    141167    Real _.m  = m_;
    142   //Real _trc.3 = { WriteLn("TRACE _.m=" << _.m); 0 };
     168  //Real _trc.3 = { WriteLn("TRACE FastCholeskiCovFactor _.m=" << _.m); 0 };
    143169    //Stores degrees of AR and MA polynomials
    144170    Real _.p = Degree(_.ar);
    145   //Real _trc.4 = { WriteLn("TRACE _.p=" << _.p); 0 };
     171  //Real _trc.4 = { WriteLn("TRACE FastCholeskiCovFactor _.p=" << _.p); 0 };
    146172    Real _.q = Degree(_.ma);
    147   //Real _trc.5 = { WriteLn("TRACE _.q=" << _.q); 0 };
     173  //Real _trc.5 = { WriteLn("TRACE FastCholeskiCovFactor _.q=" << _.q); 0 };
    148174    //Calculates de minimum length of non Toeplitz main minor of factors
    149175    Real _.n = Max(_.p,_.q)+_.p+1;
    150   //Real _trc.6 = { WriteLn("TRACE _.n=" << _.n); 0 };
     176  //Real _trc.6 = { WriteLn("TRACE FastCholeskiCovFactor _.n=" << _.n); 0 };
    151177    //Calculates de main minor of length _.n of joint ARMA covariance
    152178    VMatrix _.cov_arma_n =
    153       Pol2VMat(ARMATACov(_.ar,_.ma,_.n), _.n, _.n);
     179      ARMAProcess::DoDrop(Pol2VMat(ARMATACov(_.ar,_.ma,_.n), _.n, _.n));
    154180    //Calculates de main minor of length _.n of just pure AR covariance
    155181    VMatrix _.cov_ar_n =
    156       Pol2VMat(ARMATACov(_.ar,1,_.n), _.n, _.n);
     182      ARMAProcess::DoDrop(Pol2VMat(ARMATACov(_.ar,1,_.n), _.n, _.n));
    157183    //Does Choleski factorization of pure AR
    158   //Real _trc.7 = { WriteLn("TRACE _.cov_ar_n=" << (Matrix VMat2Mat(_.cov_ar_n))); 0 };
    159     VMatrix _.L_ar_n = CholeskiFactor(_.cov_ar_n,"X", True);
    160   //Real _trc.8 = { WriteLn("TRACE _.L_ar_n=" << (Matrix VMat2Mat(_.L_ar_n))); 0 };
     184  //Real _trc.7 = { WriteLn("TRACE FastCholeskiCovFactor _.cov_ar_n=" << (Matrix VMat2Mat(_.cov_ar_n))); 0 };
     185    VMatrix _.L_ar_n = ARMAProcess::CholeskiFactorCov(_.cov_ar_n);
     186  //Real _trc.8 = { WriteLn("TRACE FastCholeskiCovFactor _.L_ar_n=" << (Matrix VMat2Mat(_.L_ar_n))); 0 };
    161187    //Calculates the inverse of Choleski factorization of pure AR
    162     VMatrix _.Li_ar_n = CholeskiSolve(_.L_ar_n,Eye(_.n),"PtL");
    163   //Real _trc.9 = { WriteLn("TRACE _.Li_ar_n=" << (Matrix VMat2Mat(_.Li_ar_n))); 0 };
     188    VMatrix _.Li_ar_n = ARMAProcess::DoDrop(CholeskiSolve(_.L_ar_n,Eye(_.n),"PtL"));
     189  //Real _trc.9 = { WriteLn("TRACE FastCholeskiCovFactor _.Li_ar_n=" << (Matrix VMat2Mat(_.Li_ar_n))); 0 };
    164190    //Extends main minor of length _.n to full data length _.m as Toplitz from 
    165191    //last row. Terms _.n^2 is added and substracted to clean rounding error.
    166192    VMatrix _.Li_ar=
    167193      (ARMAProcess::ExtendBandedLowTriangToeplitz(_.Li_ar_n,_.m)+_.n^2)-_.n^2;
    168   //Real _trc.10 = { WriteLn("TRACE _.Li_ar=" << (Matrix VMat2Mat(_.Li_ar))); 0 };
     194  //Real _trc.10 = { WriteLn("TRACE FastCholeskiCovFactor _.Li_ar=" << (Matrix VMat2Mat(_.Li_ar))); 0 };
    169195    //Filters full ARMA covariance of AR factor
    170196    VMatrix _.cov_ma_n = (_.Li_ar_n*_.cov_arma_n*Tra(_.Li_ar_n)+_.n^2)-_.n^2;
    171197    //Extends main minor of length _.n to full data length _.m as Toplitz from
    172198    //last row and column getting a symmetric and sparse banded matrix
    173   //Real _trc.11 = { WriteLn("TRACE _.cov_ma_n=" << (Matrix VMat2Mat(_.cov_ma_n))); 0 };
     199  //Real _trc.11 = { WriteLn("TRACE FastCholeskiCovFactor _.cov_ma_n=" << (Matrix VMat2Mat(_.cov_ma_n))); 0 };
    174200    VMatrix _.cov_ma_m=
    175201      (ARMAProcess::ExtendBandedSymToeplitz(_.cov_ma_n,_.m)+_.n^2)-_.n^2;
    176   //Real _trc.12 = { WriteLn("TRACE _.cov_ma_m=" << (Matrix VMat2Mat(_.cov_ma_m))); 0 };
     202  //Real _trc.12 = { WriteLn("TRACE FastCholeskiCovFactor _.cov_ma_m=" << (Matrix VMat2Mat(_.cov_ma_m))); 0 };
    177203    //Does Choleski factorization of pseudo MA factor
    178     VMatrix _.L_ma = CholeskiFactor(_.cov_ma_m,"X",True);
     204    VMatrix _.L_ma = ARMAProcess::CholeskiFactorCov(_.cov_ma_m);
    179205    //Converts to Cholmod.R.Sparse
    180     VMatrix _.L_ma_sp = Convert(_.L_ma,"Cholmod.R.Sparse");
     206    VMatrix _.L_ma_sp = ARMAProcess::DoDrop(Convert(_.L_ma,"Cholmod.R.Sparse"));
    181207    //Calculates de logarithm of determinant of _.m ARMA covarianze
    182208    Real _.logDetCov = 2*(-VMatSum(Log(SubBand(_.Li_ar,0)))+
    … …  
    227253    Real    p     = Degree(ar);
    228254    Real    q     = Degree(ma);
     255  //WriteLn("TRACE Eval.Almagro Begin "+Time);
     256  //WriteLn("TRACE Eval.Almagro(p"<<p+",q="<<q+",m="<<m);
    229257    If(!p & !q, [[
    230258    Polyn   _.ar    = ar;
    … …  
    259287
    260288    ////////////////////////////////////////////////////////////////////////////
    261     VMatrix factor(VMatrix cov)
    262     ////////////////////////////////////////////////////////////////////////////
    263     {
    264       CholeskiFactor(cov,"X",True)
    265     };
    266     ////////////////////////////////////////////////////////////////////////////
    267289    VMatrix solve.LLt(VMatrix L, VMatrix b)
    268290    ////////////////////////////////////////////////////////////////////////////
    … …  
    283305    };
    284306
     307  //Real _trc.1 = { WriteLn("TRACE Eval.Almagro Creating _.cov_u "+Time); 0 };
    285308    VMatrix _.cov_u = Case(
    286309    _.p & _.q,
    287310    {
    288       VMatrix _.cov_zz = Pol2VMat(ARMATACov(_.ar,_.ma,_.p), _.p, _.p);
     311      VMatrix _.cov_zz = ARMAProcess::DoDrop(Pol2VMat(ARMATACov(_.ar,_.ma,_.p), _.p, _.p));
    289312    //WriteLn("_.cov_zz = "<<_.cov_zz);
    290313      VMatrix _.cov_aa = Eye(_.q);
    291       VMatrix _.cov_za = Pol2VMat(Expand(_.ma/_.ar,_.mxpq), _.p, _.q);
     314      VMatrix _.cov_za = ARMAProcess::DoDrop(Pol2VMat(Expand(_.ma/_.ar,_.mxpq), _.p, _.q));
    292315  //WriteLn("_.cov_za = "<<_.cov_za);
    293316
    … …  
    297320    _.p & !_.q,
    298321    {
    299       Pol2VMat(ARMATACov(_.ar,_.ma,_.p), _.p, _.p)
     322      ARMAProcess::DoDrop(Pol2VMat(ARMATACov(_.ar,_.ma,_.p), _.p, _.p))
    300323    //WriteLn("_.cov_zz = "<<_.cov_zz);
    301324    },
    … …  
    306329    });
    307330
     331  //Real _trc.2 = { WriteLn("TRACE Eval.Almagro Creating _.L_u "+Time); 0 };
    308332  //WriteLn("_.cov_u = "<<_.cov_u);
    309     VMatrix _.L_u  = factor(_.cov_u);
     333    VMatrix _.L_u  = ARMAProcess::CholeskiFactorCov(_.cov_u);
     334  //Real _trc.3 = { WriteLn("TRACE Eval.Almagro Creating _.covi_u "+Time); 0 };
    310335    VMatrix _.covi_u =
    311336    {
    312337      aux = solve.LLt(_.L_u,Eye(_.q_p));
    313       ARMAProcess::AvoidSymRoundErr(aux)
    314     };
    315     VMatrix _.L_u_sp = Convert(_.L_u,"Cholmod.R.Sparse");
     338      ARMAProcess::DoDrop(ARMAProcess::AvoidSymRoundErr(aux))
     339    };
     340  //Real _trc.4 = { WriteLn("TRACE Eval.Almagro Creating _.L_u_sp "+Time); 0 };
     341    VMatrix _.L_u_sp = ARMAProcess::DoDrop(Convert(_.L_u,"Cholmod.R.Sparse"));
     342  //Real _trc.5 = { WriteLn("TRACE Eval.Almagro Creating cov_G "+Time); 0 };
    316343    Real _.det_cov_u = 2*VMatSum(Log(SubBand(_.L_u_sp,0)));
    317     VMatrix _.H1 = Pol2VMat(F^_.p*(1-_.ar), _.m, _.p);
    318     VMatrix _.H2 = Pol2VMat(F^_.q*(_.ma-1), _.m, _.q);
     344    VMatrix _.H1 = ARMAProcess::DoDrop(Pol2VMat(F^_.p*(1-_.ar), _.m, _.p));
     345    VMatrix _.H2 = ARMAProcess::DoDrop(Pol2VMat(F^_.q*(_.ma-1), _.m, _.q));
    319346    VMatrix _.H = _.H1 | _.H2;
    320     VMatrix _.G =
    321       DifEq(1/_.ma, _.H, Rand(0,_.q_p,0,0), Rand(_.q,_.q_p,0,0), 0.8, 1.E-17);
    322     VMatrix _.GtG = MtMSqr(_.G);
     347    VMatrix _.G = ARMAProcess::DoDrop(
     348      DifEq(1/_.ma, _.H, Rand(0,_.q_p,0,0), Rand(_.q,_.q_p,0,0), 0.8, 1.E-17));
     349    VMatrix _.GtG = ARMAProcess::DoDrop(MtMSqr(_.G));
    323350    VMatrix _.cov_G = _.covi_u + _.GtG;
    324     VMatrix _.L_G   = factor(_.cov_G);
    325     VMatrix _.L_G_sp = Convert(_.L_G,"Cholmod.R.Sparse");
     351  //Real _trc.6 = { WriteLn("TRACE Eval.Almagro Creating _.L_G "+Time); 0 };
     352    VMatrix _.L_G   = ARMAProcess::CholeskiFactorCov(_.cov_G);
     353  //Real _trc.7 = { WriteLn("TRACE Eval.Almagro Creating _.L_G_sp "+Time); 0 };
     354    VMatrix _.L_G_sp = ARMAProcess::DoDrop(Convert(_.L_G,"Cholmod.R.Sparse"));
    326355    Real _.det_cov_G = 2*VMatSum(Log(SubBand(_.L_G_sp,0)));
    327356    Real _.det_cov_z = _.det_cov_u + _.det_cov_G;
    328357    VMatrix _.pi_z = DifEq(_.ar/_.ma,_.z);
    329358    VMatrix _.G_pi_z = Tra(Tra(_.pi_z)*_.G);
     359  //Real _trc.8 = { WriteLn("TRACE Eval.Almagro Creating various "+Time); 0 };
    330360    VMatrix _.u = solve.LLt(_.L_G, _.G_pi_z);
    331361
    … …  
    355385      _.det_cov_z+(aa+uu)/_.s2
    356386    );
     387  //Real _trc.9 = { WriteLn("TRACE Eval.Almagro End "+Time); 0 };
    357388
    358389    ////////////////////////////////////////////////////////////////////////////
Note: See TracChangeset for help on using the changeset viewer.