[Numpy-discussion] Question: LAPACK-LITE consequences for Numpy speed (win32 platfor m)
Janko Hauser
jhauser at ifm.uni-kiel.de
Tue Jan 2 09:17:48 CST 2001
First, have you tested the execution time of inverse() alone in NumPy
and Gauss, or only with the appended function. There are some ways to
speed up the function without using a specialized library. The
function also has a possible error. (see comments in the function).
In general the penality for array creation (memory allocation) is
rather big in Numpy, especially if you use long expressions, as the
result of every term is stored in a newly created array. To speed up
these kinds of operations use the inplace operations with some
workspace arrays, like:
>>> # Translating _ap[0] = _T[0] * _a0 + _c[0]
>>> _ap[0]=add(multiply(_T[0],_a0,_ap[0]),_c[0],_ap[0])
Then it should be possible to code this without a loop at all. This
would also speed up the function.
HTH,
__Janko
Shuetrim, Geoff writes:
> ________
>
> A simple version of the filter is given below:
> (Note that I have modified Matrix.py in my installation to include a
> transpose method for the Matrix class, T()).
>
> # ***************************************************************
> # kalman.py module by Geoff Shuetrim
> #
> # Please note - this code is thoroughly untested at this stage
> #
> # You may copy and use this module as you see fit with no
> # guarantee implied provided you keep this notice in all copies.
> # ***************************************************************
>
> # Minimization routines
> """kalman.py
>
> Routines to implement Kalman filtering and smoothing
> for multivariate linear state space representations
> of time-series models.
>
> Notation follows that contained in Harvey, A.C. (1989)
> "Forecasting Structural Time Series Models and the Kalman Filter".
>
> Filter --- Filter - condition on data to date
>
> """
>
> # Import the necessary modules to use NumPy
> import math
> from Numeric import *
> from LinearAlgebra import *
> from RandomArray import *
> import MLab
> from Matrix import *
>
> # Initialise module constants
> Num = Numeric
> max = MLab.max
> min = MLab.min
> abs = Num.absolute
> __version__="0.0.1"
>
> # filtration constants
> _obs = 100
> _k = 1
> _n = 1
> _s = 1
>
> # Filtration global arrays
> _y = Matrix(cumsum(standard_normal((_obs,1,_n))))
> _v = Matrix(zeros((_obs,1,_n),Float64))
> _Z = Matrix(ones((_obs,_k,_n),Float64)) + 1.0
> _d = Matrix(zeros((_obs,1,_n),Float64))
> _H = Matrix(zeros((_obs,_n,_n),Float64)) + 1.0
> _T = Matrix(zeros((_obs,_k,_k),Float64)) + 1.0
> _c = Matrix(zeros((_obs,1,_k),Float64))
> _R = Matrix(zeros((_obs,_k,_s),Float64)) + 1.0
> _Q = Matrix(zeros((_obs,_s,_s),Float64)) + 1.0
> _a = Matrix(zeros((_obs,1,_k),Float64))
> _a0 = Matrix(zeros((_k,1),Float64))
> _ap = _a
!!! Are you sure? This does not copy, but only makes a new reference
> _as = _a
!!! Same here
> _P = Matrix(zeros((_obs,_k,_k),Float64))
> _P0 = Matrix(zeros((_k,_k),Float64))
> _Pp = _P
> _Ps = _P
> _LL = Matrix(zeros((_obs,1,1),Float64))
>
> def Filter(): # Kalman filtering routine
>
> _ap[0] = _T[0] * _a0 + _c[0]
> _Pp[0] = _T[0] * _P0 * _T[0].T() + _R[0] * _Q[0] * _R[0].T()
>
> for t in range(1,_obs-1):
>
> _ap[t] = _T[t] * _a[t-1] + _c[t]
!!! You are changing _a and _as and _ap at the same time
> _Pp[t] = _T[t] * _P0 * _T[t].T() + _R[t] * _Q[t] * _R[t].T()
>
> Ft = _Z[t] * _Pp[t] * _Z[t].T() + _H[t]
> Ft_inverse = inverse(Ft)
> _v[t] = _y[t] - _Z[t] * _ap[t] - _d[t]
>
> _a[t] = _ap[t] + _Pp[t] * _Z[t].T() * Ft_inverse * _v[t].T()
> _P[t] = _Pp[t] - _Pp[t].T() * _Z[t].T() * Ft_inverse * _Z[t] *
> _Pp[t]
> _LL[t] = -0.5 * (log(2*pi) + log(determinant(Ft)) + _v[t] *
> Ft_inverse * _v[t].T())
>
> Filter()
> ____________________________________________________________________________
> ________
>
>
>
> **********************************************************************
> " This email is intended only for the use of the individual or entity
> named above and may contain information that is confidential and
> privileged. If you are not the intended recipient, you are hereby
> notified that any dissemination, distribution or copying of this
> Email is strictly prohibited. When addressed to our clients, any
> opinions or advice contained in this Email are subject to the
> terms and conditions expressed in the governing KPMG client
> engagement letter. If you have received this Email in error, please
> notify us immediately by return email or telephone +61 2 93357000
> and destroy the original message. Thank You. "
> **********************************************************************...
>
> _______________________________________________
> Numpy-discussion mailing list
> Numpy-discussion at lists.sourceforge.net
> http://lists.sourceforge.net/mailman/listinfo/numpy-discussion
More information about the Numpy-discussion
mailing list