[Numpy-discussion] Question: LAPACK-LITE consequences for Numpy speed (win32 platfor m)
Shuetrim, Geoff
gshuetrim at kpmg.com.au
Tue Jan 2 00:56:21 CST 2001
Apologies for asking an overly vague question like this but:
on an Intel/win32 platform (I only have a windows version of
Gauss), I am comparing Numpy matrix inversion with that
of Gauss (very much the same type of software as Matlab
which at least some of you used to work with). As the size
of the matrix to be inverted increases, the speed of Numpy
appears to asymptote (on my machine) to about half that of
Gauss. For small matrices, it is much worse than that
because of the various overheads that are dealt with by Numpy.
Would this speed differential be largely eliminated if I was not
using LAPACK-LITE? If so I will try to figure my way through
hooking into Intel's MKL - anyone got hints on doing this - I saw
mention of it in the mailing list archives. Would I be better off,
speed-wise, eschewing win32 altogether and using native LAPACK
and BLAS libraries on my Linux box?
This is relevant to me in the context of a multivariate Kalman filtering
module that I am working up to replace one I have been using on
the Gauss platform for years. The Numpy version of my filter has a
very similar logic structure to that of Gauss but is drastically slower.
I have only been working with Numpy for a month or so which may
mean that my code is relatively innefficient. I have been assuming
that Numpy - as an interpreted language is mainly slowed by
looping structures.
Thanks in advance,
Geoff Shuetrim
____________________________________________________________________________
________
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
_as = _a
_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]
_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. "
**********************************************************************...
More information about the Numpy-discussion
mailing list