Thanks George for the update. I agree with you that the next step is finer profiling of the code. You may want to start by profiling the Matlab code. It will give you an idea of which steps takes the longer and compare with the time it takes in C++ to do the same thing.
I'm a bit puzzled by your comparison of time it takes to do a matrix inversion. In your previous message, you mentioned that most of C++ time was spend in matrix inversion (on the small example). Could you please check that it is also the case with the Matlab code? I always thought that the longest operation was the update of the variance matrix (P). Maybe it depends upon the size of the problem.
Best
Michel
G. Perendia wrote:
Dear Michel
The tortoise and the hare race seem to be closing-up ... but slowly:
- Using the larger sw_euro 3 countries model , with 57x57 system matrix
T (though, this time on my desktop machine since the laptop "melted down" due to overheating), Matlab Dynare KF was, this time, leading "only" around 45% over the execution time needed for the C++DLLs, one called repeatedly externally, the other within its inner loop:
1000 calls to matlab KF time = 97.5000 1000 calls to dll time = 143.8590 1000 innrloop_dll time = 139.4370
(i.e. this is "catching-up" in comparison with the 8x8 matrix test done earlier which provided the Matlab (the "tortoise"?) with the lead of around 250% )
- The C++DLL delay does not appear to be due to the slower matrix
inversion: a 57x57 Matrix inversion in 100000 loops showed that a C++ dll, using GeneralMatrix.InvLeft (which calls ztrsv), running in its inner loop, was marginally quicker:
- Inner loop dll_mx_inv =etime(clock, t)= 117.8690
- Matlab_mx_inv=etime(clock, t) = 119.6220
(....though a bit less precise too!)
- Also, the delay is not due to the shortcut he Matlab KF does with cov.
conversion since the steady state is not reached for the large matrix. However, this feature would be a later stage, further improvement to C++ KF which, as additional tests show, appear not to be currently taking advantage of that feature.
... and I am now in doubt who is the hare and who the tortoise and whether if the "hare" can really overtake at any reasonably sized matrix.
- However, the C++ KF is highly modular and rich in functionality, so that
its performance could probably benefit from some refactoring (I have already identified some spots), and, of course, of having Andrea's code integrated.
I think I should focus on profiling and analysing its inner-workings to figure out why is C++ KF so much slower before I make any substantial changes.
Best regards George
----- Original Message ----- From: "Michel Juillard" michel.juillard@ens.fr To: "List for Dynare developers" dev@dynare.org Sent: Thursday, June 04, 2009 3:09 PM Subject: Re: [DynareDev] Kalman Filter
Dear George
1)My first idea is that calls to BLAS and dgemm are too heavy for small matrices.
- Can you compare the time spent inversing the matrix in C++ and in
Matlab and make sure that the two matrices have indeed the same size?
- I'm attaching an example made of Smets and Wouters piled up for 3
countries with artificial data. It doesn't make and economic sense but represents a medium to large model, a bit smaller than what I had in mind.
Best,
Michel
G. Perendia wrote:
Hi
- Profiling results:
Initial profiling of both the stand alone Executable test and DLL run in Matlab thread indicates the 10,000 loop is spending far the most of time
in
ztrsv (~40%) (and, in some instances, its peer dtrsv ?) BLAS triangular matrices equation solver(s), that appear to be extensively used to invert the KF matrices (although direct trace to the caller functions was unclear and it is unknown most of the time but it appears to be dtrsv but also dgemm(?)), followed by around 7.5% equally in:
Wcstombs lsame freebuffers modfthat appear to be invoked by ztrsv or similar BLAS routines such as caxpy and dtrsv.
- Performance tests:
Despite a search, I still could not find a large, around 100 variable
model
suitable for performance testing KF on large models so please let me know
if
you have anything suitable.
PS: Some more performance tests:
test 6) a stand alone executable (.exe) test version of C++ in inner
10,000
loop (with simulated data) took nearly the same as running the DLL from Matlab ~ 2’ 40” or ~ 160” (i.e. in a repeated test on a refreshed
machine,
Matlab DLL had a bit shorter a but still high execution time: innrloop_dll_time = 154.9330")
test 7) Overall Estimation of ls2003.mod took about 17" or 18% longer with DLL KF than with Matlab KF: 108" vs 91"
Best regards
George
----- Original Message ----- From: "G. Perendia" george@perendia.orangehome.co.uk To: "List for Dynare developers" dev@dynare.org Sent: Monday, June 01, 2009 6:19 PM Subject: Re: [DynareDev] Kalman Filter
I used a rather small model with 4 observables, T with 8 state out of
total
9 endogenous variables only.
I will try a larger one once I find a suitable one - do you have any suggestions or examples of 100 state var model I could use?
Best regards
George
----- Original Message ----- From: "Michel Juillard" michel.juillard@ens.fr To: "G. Perendia" george@perendia.orangehome.co.uk; "List for Dynare developers" dev@dynare.org; "'Andrea Pagano'" andrea.pagano@jrc.it Sent: Monday, June 01, 2009 6:00 PM Subject: Re: [DynareDev] Kalman Filter
One more thing: what is the size of your testing model. How many variables total ? How many state variables? How many observed variables?
Can you redo the experiments 3 and 4 for a model with more than 100 state variables? Probably 100 calls will be enough
Best
Michel
G. Perendia wrote:
Dear Michel
Results of run lapse time using etime(clock,t) are rather surprising
and
very disappointing (I had to repeat some in a sheer disbelief) - far
the
best is Matlab Dynare kalman_filter.m (even after applying some
refactoring
of the kalman.cpp's deep inner loop, not applied to all DLL tests
yet)!!!.
Timing tests on Pentium 2.4 GHz
1)calling Kalman_filters.m in the 10000 loop, itself calling the
optimised
version of DLL with preparation of H and Z matrices in each loop: total_dll_time = 202.7320
2)calling core optimised version of Kalman_filter_dll in the 10000
loop
after (i.e., loop without) preparation of H&Z matrices in
kalman_filters.m:
core_dll_time= = 1st run: 172.5890 2nd run: 195.0500 at ~93% of CPU
time
3)calling the Kalman_filter_loop.dll with its own inner 10000 loop
from
kalman_filters.m (including the one-off call time) innrloop_dll_time = 1st run 197.8640 2nd run: 192.5970 3rd run (after some refactoring of the kalman.cpp filter deep inner
loop
not yet used in other DLL tests): 184.9760
4)calling Dynare SVN kalman_filter.m in 10000 loop matlabKF_time = 1st run: 48.9600 2nd run: 48.4600
- for a comparison and sanity check: calling Kalman_filters running
the
debug version of DLL in the 10000 loop dbug total_dll_time=etime(clock, t) = 642.4140 running at ~93% of CPU
time
It is possible that using proper, optimised lapack and BLAS may give a better performance for C++. (NOTE: the DLLs used Matlab lapack and
blas
libraries.)
Shall I nevertheless continue with the C++ version of DsgeLikelihood
as
planned or try to revise the performance tests first?
Best regards
George artilogica@btconnect.com
----- Original Message ----- From: "Michel Juillard" michel.juillard@ens.fr To: "G. Perendia" george@perendia.orangehome.co.uk Sent: Monday, June 01, 2009 10:32 AM Subject: Re: [DynareDev] Kalman Filter
Thanks George, now it is clear. I'm also interested by a timing comparison. Could you time a) 10'000 calls to the Matlab filter routine b) 10'000 calls to the C++ filter routine c) add the 10'000 loop inside the DLL and time one call to this
modified
function.
I guess after that we should pass to the coding in C++ of
DsgeLikelihood.m
All the best,
Michel
G. Perendia wrote:
> Dear Michel > > No, I do not think there is additional problem (other than with > >
Diffuse
> Likelihood*.m files when if used) or that correction was not > >
sufficient.
> Over the weekend I was just trying full 1st stage estimation with a > > >
model we
> use for Part Info, with and without the corrections, which just > >
shows
the
> scale of the effect of the log likelihood calculation error (now > > >
properly
> corrected for the SVN likelihood files) . As you can see in the .rtf > > >
file,
> in the few of the tests with SVN version (using > > >
/likelihood/kalman.filter.m)
> I used your new code too, but for comparing results in older > >
versions
of
> Dynare- v4.0.2 and 4.0.3 - (using Diffuse Likelihood1.m), I used my > > >
quick
> fix applied to the 4.02 DiffuseLikelihod*.m that I also enclosed. > > Also, some additional tests I did over the weekend show that, after > >
the
> correction, Dynare now reports estimates more comparable to the new > >
Part
> Info method than before the correction, also proving that the > >
correction
> works ! > > As you suggested, I run a test of your new correction that can be > > >
compared
> to C++ and, for the test model I used, F converges at t=3 (start is > >
41),
and
> your new kalman_filter.m code gives total > LIK=1640935.5855619563 > which is virtually the same as C++ and the quick fix for Diffuse > > >
Likelihood
> I tested last week (requoted from below): > > > > > >>>>> the reported likelihood for presample=40 from >>>>> >>>>>
(quick
> fix) Matlab KF > > > > >>>>> 1640935.5855267849 >>>>> which is nearly the same as that from C++ KF >>>>> >>>>>
below:
>>>>> 1640935.5854489324 >>>>> >>>>> >>>>> >>>>> > Let me know if you want me to run more tests or work on correcting > >
any
other
> Dynare Matlab code! > > Best regards > > George > > ----- Original Message ----- > From: "Michel Juillard" michel.juillard@ens.fr > To: "G. Perendia" george@perendia.orangehome.co.uk > Sent: Monday, June 01, 2009 9:08 AM > Subject: Re: [DynareDev] Kalman Filter > > > > > > >> Dear George, >> >> I'm a bit confused. We agreed that all the versions of Kalman >> likelihood, both in ./matlab and in ./matlab/kalman/likelihood were >> wrong when presample was larger than the period where the filter >> >> >>
reaches
>> its stationary state. So I don't see the point in still using this >> >>
code
>> for testing. >> >> the code in matlab/kalman/likelihood was corrected last Wednesday >> >>
on
SVN
>> r2007 >> >> Do you have reasons to believe that the correction wasn't >> >>
sufficient?
Or
>> is there another source of discrepancy between the Matlab code and >> >>
the
>> C++ code? >> Before comparing results of entire estimation procedures, we have >> >>
to
>> make sure that all the outputs from the Matlab and C++ functions >> >>
with
>> the same inputs are close enough. >> >> Note that we still have to correct DsgeLikelihood_hh.m that calls >> >>
the
>> obsolete DiffuseLikelihood*.m routines instead of the ones in >> kalman/likelihood >> >> Then we should just remove the DiffuseLikelihood*.m functions >> >> >> Best >> >> Michel >> >> >> G. Perendia wrote: >> >> >> >> >>> Dear Michel >>> >>> Problem and discrepancies with log likelihood calculation with >>> >>> >>> >>> > presampling > > > > >>> (start>1) is something I have been encountering for some time >>> >>>
(i.e.
few
>>> weeks) already , whilst doing some tests for the new partial >>> >>> >>>
information
>>> kalman estimation, but, expecting that the problem is in the new >>> >>> >>>
Partial
>>> info code, I did not get to the bottom of the cause until I >>> >>> >>>
encountered
>>> discrepancy between the C++ KF (modified to handle presampling >>> >>> >>>
properly)
> and > > > > >>> the Dynare KF too. It was that second discrepancy that led me to >>> >>> >>>
analyse
>>> Dynare KF handling of presampling and to find the problem there. >>> >>> I did some more preliminary tests yesterday with the older and the >>> >>> >>> >>> > latest > > > > >>> SVN versions of Dynare, without and with the correction (including >>> correction using the fix based around the min() function I >>> >>>
initially
>>> suggested and now applied as a quick fix to the older >>> >>> >>> >>> > DiffuseLikelihood1.m) > > > > >>> and the estimation results with and without correction are >>> >>> >>>
significantly
>>> different as can be seen from the enclosed rtf document which uses >>> >>>
two
> very > > > > >>> similar models, (i.e. one with gamma>0 the other with gamma=0), to >>> >>> >>> >>> > compare > > > > >>> the results (also very similar to the one I sent to Stephane). >>> >>> Although the initial calculations of log-likelihood ("Initial >>> >>>
value
of
> the > > > > >>> log posterior (or likelihood): ") are same, (i.e. when reste is 0 >>> >>>
or
> very > > > > >>> small), it appears that the differences builds up during the >>> >>> >>>
estimation
> as > > > > >>> the kalman filter cov matrix start to converges earlier and >>> >>>
earlier,
>>> eventually before the given presample start is reached, as the >>> >>>
later
in
> the > > > > >>> ML and estimation error minimisation we are (and probably >>> >>>
throughout
the
>>> MCMC estimation process), when the parameters get close to the >>> >>>
mode
>>> (especially if the parameter priors were taken from previously >>> >>> >>>
estimated
>>> posteriors as is the case with the test models I used). >>> >>> (I did not yet run MCMC due to estimation failing with the chol() >>> non-positive definite error so I would need to change starting >>> >>> >>>
values.)
>>> Let me know if you would like me to develop and perform more >>> >>> >>>
structured
>>> tests on a large range of models. >>> >>> >>> Best regards >>> >>> George >>> >>> ----- Original Message ----- >>> From: "Michel Juillard" michel.juillard@ens.fr >>> To: "List for Dynare developers" dev@dynare.org >>> Cc: "G. Perendia" george@perendia.orangehome.co.uk >>> Sent: Wednesday, May 27, 2009 12:16 PM >>> Subject: Re: [DynareDev] Kalman Filter >>> >>> >>> >>> >>> >>> >>> >>>> In fact the discrepancy found by George occurs only for start ~= >>>> >>>>
1
>>>> So it is clearly related to the handling of the constants. >>>> >>>> Best >>>> >>>> Michel >>>> >>>> >>>> Stéphane Adjemian wrote: >>>> >>>> >>>> >>>> >>>> >>>>> George, Can you share your example. In my experience (with >>>>> >>>>>
medium
>>>>> scaled models) it takes much more iterations to get to the >>>>> >>>>>
steady
>>>>> state kalman filter... >>>>> >>>>> I didn't look at the cc codes, but one possible source of >>>>> >>>>> >>>>>
discrepancy
>>>>> between matlab and cc may be the steady state KF criterion. I >>>>> >>>>>
test
for
>>>>> the stationarity of KF considering the changes in the biggest >>>>> >>>>> >>>>>
elements
>>>>> of the Kalman Gain matrix, not the changes in the determinant of >>>>> >>>>>
the
>>>>> forecast errors covariance matrix... >>>>> >>>>> Best, >>>>> Stéphane. >>>>> >>>>> >>>>> 2009/5/27 G. Perendia <george@perendia.orangehome.co.uk >>>>> mailto:george@perendia.orangehome.co.uk> >>>>> >>>>> >>>>> 1. In my experience, the filter cov. determinants dF becomes >>>>> stationary after several steps (5-10), e.g. 7 >>>>> and number of periods the filter is stationary is >>>>> reste=smpl-7 >>>>> whilst, if start =41 (presample=40) than reste is then >>>>> >>>>>
larger
than
>>>>> smpl -start+1=smpl-40. >>>>> >>>>> 2. it is not problematic to add all the determinants at >>>>> >>>>>
once
in
>>>>> the last period of the stationary filter >>>>> as long as we subtract from the total sum >>>>> >>>>> start-1-(smpl-reste) = 41 -1 -(smpl-smpl+7)=33 >>>>> determinants, or, add only >>>>> min(reste,(smpl-start+1)) >>>>> Best regards >>>>> >>>>> George >>>>> >>>>> ----- Original Message ----- >>>>> *From:* Stéphane Adjemian >>>>> >>>>> >>>>>
mailto:stephane.adjemian@gmail.com
>>>>> *To:* List for Dynare developers mailto:dev@dynare.org >>>>> *Cc:* G. Perendia >>>>> >>>>>
mailto:george@perendia.orangehome.co.uk
>>>>> *Sent:* Wednesday, May 27, 2009 11:04 AM >>>>> *Subject:* Re: [DynareDev] Kalman Filter >>>>> >>>>> Hi all, >>>>> >>>>> I agree, the matlab code is very unclear (even if I had >>>>> >>>>>
fun
>>>>> writting it this way ;-) and prone to errors if one uses >>>>> >>>>>
the
>>>>> vector lik (Marco is using it). I would rather prefer to >>>>> >>>>>
add
>>>>> the constants outside of the loop with a (sub)vector >>>>> operation, this should be more efficient. I will do it >>>>> >>>>>
today
>>>>> or tomorrow. >>>>> >>>>> Best, >>>>> Stéphane. >>>>> >>>>> 2009/5/27 Michel Juillard <michel.juillard@ens.fr >>>>> mailto:michel.juillard@ens.fr> >>>>> >>>>> On closer inspection, I don't think that the >>>>> >>>>>
expression
>>>>> pointed by George in kalman_filter.m is wrong: >>>>> >>>>> 1. reste = smpl-t or the number of periods during >>>>> >>>>>
which
>>>>> the filter is stationary. This shouldn't be larger >>>>> >>>>>
than
>>>>> T-start+1 >>>>> >>>>> 2. it is problematic (see below) but not wrong to >>>>> >>>>>
add
all
>>>>> the determinants at once in the last period of the >>>>> stationary filter >>>>> >>>>> 3. I don't think this explains the difference with >>>>> >>>>>
the
C++
>>>>> version of the filter and we still have to look for >>>>> >>>>>
it.
>>>>> 4. it remains that the current code is very unclear >>>>> >>>>>
and
>>>>> that if LIK is correct the vector lik doesn't have >>>>> >>>>>
the
>>>>> correct constants on each elements. >>>>> >>>>> 5. I would like to simplify the code and add the >>>>> >>>>>
correct
>>>>> constant to each element of the lik vector. It would >>>>> >>>>>
be
a
>>>>> little bit less efficient in Matlab than the current >>>>> >>>>> >>>>>
code,
>>>>> but I doubt it would be noticeable. >>>>> Stephane, what do you think? >>>>> >>>>> >>>>> Best >>>>> >>>>> Michel >>>>> >>>>> G. Perendia wrote: >>>>> >>>>> Dear Michel >>>>> >>>>> I think I found an error in Dynare Matlab >>>>> kalman_filter. suite of utilities >>>>> which affects the likelihood LIK results with >>>>> >>>>> >>>>>
start>1
>>>>> (i.e. presampling>0): >>>>> >>>>> the calculation speed-up construct which relies >>>>> >>>>>
on
>>>>> converged covariance >>>>> matrix >>>>> >>>>> lik(t) = lik(t) + reste*log(dF); >>>>> >>>>> adds reste * log(dF) to the last-1 (i.e. the >>>>> >>>>>
smpl)
>>>>> member of lik >>>>> (the last, the lik(smpl+1) one contains >>>>> >>>>> >>>>> >>>>> >>>>> >>> smpl*pp*log(2*pi)) >>> >>> >>> >>> >>> >>>>> but reste is usually larger than T-start+1 so >>>>> >>>>>
that
>>>>> LIK = >>>>> >>>>> >>>>> >>>>> >>>>> >>> .5*(sum(lik(start:end))-(start-1)*lik(smpl+1)/smpl) >>> >>> >>> >>> >>> >>>>> has much more log(dF)s added than required since >>>>> >>>>> >>>>>
they
>>>>> are all concentrated >>>>> in the last-1 (the T) member >>>>> >>>>> For example, if I change the above construct to >>>>> lik(t) = lik(t) + >>>>> >>>>>
min(reste,(smpl-start+1))*log(dF);
>>>>> the reported likelihood for presample=40 from >>>>> >>>>>
Matlab
> KF > > > > >>> is >>> >>> >>> >>> >>> >>>>> 1640935.5855267849 >>>>> which is nearly the same as that from C++ KF >>>>> >>>>>
below:
>>>>> 1640935.5854489324 >>>>> >>>>> Shall I make changes to kalman/likelihood/ KFs >>>>> >>>>>
and
>>>>> upload the .m files? >>>>> This problem affects also the older versions of >>>>> DiffuseLikelihood**.m too. >>>>> >>>>> Best regards >>>>> >>>>> George >>>>> artilogica@btconnect.com >>>>> >>>>> >>>>> >>>>> >>>>> >>> mailto:artilogica@btconnect.com >>> >>> >>> >>> >>> >>>>> ----- Original Message ----- From: "Michel >>>>> >>>>>
Juillard"
>>>>> <michel.juillard@ens.fr >>>>> >>>>> >>>>> >>>>> > mailto:michel.juillard@ens.fr> > > > > >>>>> To: "G. Perendia" >>>>> >>>>>
<george@perendia.orangehome.co.uk
>>>>> mailto:george@perendia.orangehome.co.uk> >>>>> Sent: Tuesday, May 26, 2009 10:32 AM >>>>> Subject: Re: Kalman Filter+PS >>>>> >>>>> >>>>> >>>>> >>>>> Hi George, >>>>> >>>>> >>>>> >>>>> >>>>> Re 1) below: >>>>> I modified C++ KF so that it reports >>>>> log-likelihood for given >>>>> start/preampling in same/similar manner >>>>> >>>>>
as
the
>>>>> Matlab KFs do and I am >>>>> getting approximately close results, >>>>> >>>>>
e.g.
>>>>> ll= -1640935.5854489324 for C++ and >>>>> (-) 1640482.4179242959 for Matlab KF >>>>> >>>>>
(for
>>>>> start=41, i.e. >>>>> >>>>> >>>>> presample=40). >>>>> >>>>> >>>>> whilst they appear same for presample=0 >>>>> (e.g.2.5906e+006), i.e. >>>>> -2590556.989730841 vs >>>>> 2590556.989778722 >>>>> >>>>> Are those results acceptably close or >>>>> >>>>>
should
I
>>>>> investigate further where >>>>> >>>>> >>>>> the >>>>> >>>>> >>>>> above difference may come form? >>>>> >>>>> >>>>> >>>>> >>>>> This indicates a problem . The difference >>>>> >>>>>
should
>>>>> be the same with and >>>>> without presample. It may come from the >>>>> computation of the likelihood >>>>> constant. This is done in a very obscure >>>>> >>>>>
manner
in
>>>>> Dynare Matlab. >>>>> >>>>> >>>>> >>>>> >>>>> >>>>> >>>>> >>>>> >>>>> >>>>> _______________________________________________ >>>>> Dev mailing list >>>>> Dev@dynare.org mailto:Dev@dynare.org >>>>> http://www.dynare.org/cgi-bin/mailman/listinfo/dev >>>>> >>>>> >>>>> >>>>> >>>>> -- >>>>> Stéphane Adjemian >>>>> CEPREMAP & Université du Maine >>>>> >>>>> Tel: (33)1-43-13-62-39 >>>>> >>>>> >>>>> >>>>> >>>>> >>>>>
>>>>> _______________________________________________ >>>>> Dev mailing list >>>>> Dev@dynare.org mailto:Dev@dynare.org >>>>> http://www.dynare.org/cgi-bin/mailman/listinfo/dev >>>>> >>>>> >>>>> _______________________________________________ >>>>> Dev mailing list >>>>> Dev@dynare.org mailto:Dev@dynare.org >>>>> http://www.dynare.org/cgi-bin/mailman/listinfo/dev >>>>> >>>>> >>>>> >>>>> >>>>> -- >>>>> Stéphane Adjemian >>>>> CEPREMAP & Université du Maine >>>>> >>>>> Tel: (33)1-43-13-62-39 >>>>> >>>>> >>>>> >>>>> > --------------------------------------------------------------------- >
>>>>> _______________________________________________ >>>>> Dev mailing list >>>>> Dev@dynare.org >>>>> http://www.dynare.org/cgi-bin/mailman/listinfo/dev >>>>> >>>>> >>>>> >>>>> >>>>> >>>>> >>>> _______________________________________________ >>>> Dev mailing list >>>> Dev@dynare.org >>>> http://www.dynare.org/cgi-bin/mailman/listinfo/dev >>>> >>>> >>>> >>>> >>>> >>>> >
Dev mailing list Dev@dynare.org http://www.dynare.org/cgi-bin/mailman/listinfo/dev
Dev mailing list Dev@dynare.org http://www.dynare.org/cgi-bin/mailman/listinfo/dev
Dev mailing list Dev@dynare.org http://www.dynare.org/cgi-bin/mailman/listinfo/dev
Dev mailing list Dev@dynare.org http://www.dynare.org/cgi-bin/mailman/listinfo/dev
Dev mailing list Dev@dynare.org http://www.dynare.org/cgi-bin/mailman/listinfo/dev