[Numpy-discussion] Question: LAPACK-LITE consequences for Numpy speed (win32 platfor m)

[Shuetrim, Geoff]

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. "
**********************************************************************...