5.6 Parallel computing
The scale and magnitude of the computational effort in estimating parameters are extensive. The combination of an optimization of several optimizations generates a large need for computational power. Adding to this workload the program language Matlab results in longer execution times compared to other languages. An implementation in another language would increase the speed significantly but it would still be a large computational task.
The preferred solution and what seems to be the trend is to use parallel com-puting. There are several areas in the PSM prototype where parallel computing would be beneficial. It has been chosen to implement parallelization in the prototype in the split from overall minimization to minimization pr. individ-ual. The parallelization results inN separate calculations which is seen clearly by all the arrows in Figure 5.2. Furthermore, but not implemented it would be advantageous to introduce parallel execution of the numerically calculated gradients.
Matlab does not have the option of parallel computing by default. MatMPI7 enables parallel computing in Matlab by creating a set of scripts that is executed in separate processes. MatMPI uses message parsing but it was found faster to parse all data and parameters through files. The individual calculation extracts its unique part of data by using its identifier number. The individual log-likelihood result is parsed back into the leader thread by proper message parsing to avoid deadlocks or race conditions. A shared memory environment is required as the message parsing is implemented through shared files.
MatMPI works by issuing commands in blocks. This means that all working processors are started at the same time and next block of data is not started until all have finished. Some consideration should be put into matching the number of CPUs with the data or simply just choosing the maximum number of CPUs available.
In Figure5.3 a simulated data set consisting of 20 series of objects in free fall was generated. A single APL was calculated for a different number of working processors. The server handled a large number of jobs simultaneously for other users which have influenced the computation times. It can be seen that the computation time is reduced to around a fifth of the original. In this example it appears that the overhead is relatively large, however for more computational intensive models the benefit of adding more CPUs may be substantial.
7MatMPI is a package to Matlab - developed at MIT Lincoln Laboratory by Dr. Jeremy Kepner. - http://www.ll.mit.edu/MatlabMPI/
0 2 4 6 8 10 12 14 16 18 20 0
50 100 150 200 250
CPUs
Secs
Data: 20 subjects, 26 samplepoints, η∈ R2
Figure 5.3: Computation Times using Parallel Computing
PSM summary
A prototype has been developed in Matlab implementing a non-linear mixed-effect model with SDEs. The numerical implementation is based on an ODE-solver, minimizers and own code.
Functions for estimating the parameter uncertainty have been implemented.
The Hessian calculation has been analysed and optimized to avoid unnecessary calculations of the costly APL function.
The computational work in solving parameter estimation in the prototype is substantial but parallelization has been shown to improve execution time. A more thorough validation of the implemented prototype is performed in the next chapter.
Chapter 6
Validation of PSM
The main purpose of the validation of PSM is to ensure that the results produced can be trusted. This will be done by comparing with existing PK/PD software programs. However, no single program can handle non-linear mixed effects models with SDEs and for this reason the prototype is validated against two different programs. The first is CTSM (Continuous Time Stochastic Modelling) is able to handle stochastic state space models for a single subject. The second program NONMEM is able to handle non-linear mixed effects model but only in ODE context.
The implemented prototype consists of several functions which together are able to handle both the stochastic state space model and the non-linear mixed effect part. Many of the functions are interlinked making it possible to test several functions at the same time as it is unlikely that erroneous subfunctions will produce a correct result in the overall function.
The validation will be done in two steps. Firstly, the Kalman filter will be verified by comparing results with CTSM using a model based on SDEs on a single individual. Secondly, the APL function will be compared with NONMEM.
This can only be done with ODEs in NONMEM, and thus the validation is based upon that the SDE extension is validated against CTSM.
CTSM 259.028450 PSM EKF 259.027492 PSM KF 259.027543
Table 6.1: Comparison of CTSM and PSM log-likelihood functions.
It is also needed to validate that the parameter estimation yields correct popula-tion parameters as well as individual parameters. Finally, the Kalman smoother will be validated.
All of the above must further more be carried out for both the linear and non-linear versions of the implementation.
6.1 Kalman filter
The Kalman filter is the basis for the implemented population model since the likelihood function is based on the one-step predictions and variances. The Kalman filter in PSM is validated against CTSM by comparing the values of the likelihood functions for one individual, see Eq. (3.28). This is done by comparing with a two-compartment C-Peptide model which will be described in further detail in Section7.2. The data are taken from subject 1.
The model parameters to be estimated are the initial states,C10, C20, ISR0, mea-surement noiseS, and Wiener magnitude for the secretion rateσISR. Normally, the initial state is constrained to steady state giving just one parameter to estimate, but it is not possible to set this constraint in CTSM. The model is es-timated using CTSM, and the resulting estimates are ˆθ= (C1, C2, ISR, S, σISR)
= (906.62,−8.7908,104.58,13763,7.1796). It is noted thatC2is estimated to a negative initial concentration, which is of course not physically possible. How-ever, the estimated standard deviation in CTSM is 29.8 and it can therefore not be assumed significantly different from 0.
The estimated parameter values from CTSM are used with both the non-linear and linear Kalman filters in PSM and the resulting likelihood values are shown in Table 6.1. It is seen that the values deviate but they are still acceptably close to being equal. The deviation can be caused by numerical differences in implementation or that the parameters from CTSM were only available with 5 significant digits.
As mentioned the likelihood value is only based on the predictions from the Kalman filter. To verify the filtered values from the Kalman filter, it is compared
6.1 Kalman filter 43
to CTSM by plotting the secretion rate (3rd state) and by visual inspection it is seen that both the non-linear and linear Kalman filters in PSM yields filtered estimates approximately equal to CTSM.
6.1.1 Kalman smoother
The Kalman smoother is validated against CTSM using the same example as in the previous section. Smoothed estimates are generated in CTSM and in Fig-ure 6.1 the smoothed ISR-state is compared to the non-linear PSM smoother estimates. By looking at the figure it is seen that they are practically identical in estimate as well as standard deviation. CTSM only returns values for obser-vation time points whereas PSM has inserted subpoints between obserobser-vations.
The comparison should only be based on the value at observation time points.
0 120 360 600 840 1140 1440
50 100 150 200 250
PSM EKF CTSM
0 120 360 600 840 1140 1440
0 10 20 30 40 50
PSM EKF CTSM Smoothed ISR st.dev.
Smoothed ISR state
Figure 6.1: Comparison of EKF smoothed secretion rate values.
The linear version of the Kalman Smoother based on the Bryson-Frazier (B-F) formulas is also tested against CTSM and the output is shown in Figure 6.2.
At a first glance this also corresponds with the CTSM smoother. However, there seems to be some small numerical problems with the method just before T = 600s for this model, see Figure 6.3. It is expected that the smoothed estimate of the state and standard deviation are continuous, and not exhibit the discontinuous behavior seen in the figure. However, right after T = 600 the
B-F method returns to the CTSM estimates. A thorough investigation of the problem should be carried out, but it has not been done as it is considered of minor interest.
0 120 360 600 840 1140 1440
50 100 150 200 250
PSM KF(BF) CTSM
0 120 360 600 840 1140 1440
0 10 20 30 40 50
PSM KF(BF) CTSM Smoothed ISR st.dev.
Smoothed ISR state
Figure 6.2: Comparison of KF smoothed secretion rate values.
The Kalman filter and smoother have been successfully validated with CTSM as reference. This is considered sufficient for the present need and in line with focus of the thesis.