# [PYTHON MATRIX-SIG] Speed...cg...and more!

Thomas Schwaller et@appl-math.tu-muenchen.de
Mon, 13 Nov 1995 12:57:32 +0100

```Hi all. Unfortunately I missed the first 5 weeks or so :-( of the
very interesting but sometimes quite heavy discussion!
So let me begin slow.... ;-)

1) The ofunc object is very interesting, but I miss a feature which I always
missed also in the mathmodule.c of the Python distribution. When I want to
use e.g. sin(a) of an not predefined type (the new math functions in
Matrix just extend the standard ones by comples numbers, matrices,..)
I would like that sin (and all the other ones) test at last if
the object has a sin attribute and execute it. In Mathmodule.c this was
just a minor hack but here I don't now. Why do I want that behaviour.
I write a lot of math stuff (ince contributes a 1D automatic differentiation
module) and I also want to be able to call sin(a) of these types.
In case of the Matrix objects one problem we will ultimately be faced with
is the problem of using e.g. exp(A) but not in the pointwise sense
but in the mathematical sense (defined by power series for example)
So what to do then. Extending exp in python is done very eay:

from Matrix import *

Matrix_sin=sin

def sin(x):
"""uses builtin sin function or attribute sin"""
try:
return Matrix_sin(x)
except AttributeError:
if hasattr(x, 'sin'): return x.sin()
else: raise AttributeError, 'sin not available for this type'

but I think this should be in the python code, so I just can do something like

class dummy:
def __init__(self, x)
self.x=x

def sin(self):
return sin(x), sin(x-1)

a=dummy()
print sin(a)

This is a thing which should be there once and for all.

2) When writing mathematical stuff with matrices one uses very often
the multiplication (in the Matrix sense) with transposed matrices.
Something like:
s  = Multiply(A.transpose(), r)

Is it possible to write that directly with the available methods.
When doing this kind of operation you don't really need to trnaspose
the matrix, just switch the order of the the loops. When you do a lot of
mulitplications with a lot of transposed matrices this is quite important.
Please tell me when I'm missing something here. As far as I understand, at
the moment you can influence the indices, but can you also influence the
order of execution whitout producing a new matrix on the fly?
As a goodie here's a simple conjugate gradient method for testing:
(Hope it' correct!)

from Matrix import *
import sys

def cg(A, b, x=None, eps=1e-12, k=None, printer=None):
m, n = A.shape
if k == None: k = m*2   		# indeed much too high!
if x == None: x, r = zeros(n,1), b
else : r = b - Multiply(A, x)
p = Multiply(A.transpose(), r)
gamma = dot(p, p)
for j in range(k):
q = Multiply(A, p)
alpha = gamma / dot(q, q)
x  = x + alpha*p
r  = r - alpha*q
s  = Multiply(A.transpose(), r)
gamma_new = dot(s, s)
beta = gamma_new / gamma
gamma = gamma_new
if gamma <= eps*eps: return x
p = s + beta*p
if printer: printer(x, r, s, p, q, beta, gamma)
raise error,'cg-method diverges with %d iterations and eps = %.0e' %
(k,eps)

def cg_test():
def pr(x, r, s, p, q, beta, gamma):
print gamma
n =20
b= ones(n, 1)
for i in range(100):
a = rand(n, n)
a=0.5*(a+a.transpose())
x = cg(a, b)
sys.stdout.write('.')
sys.stdout.flush()

cg_test()

3) As a matter of style, perhaps we should discuss how to write such things.
I guess as soon as the matrixmodule will be in the standard distribution
there will be a lot of matlab style writing. We should try to keep the
procedures, classes ... as uniform as possible, so that everybody can
rely on a certain behaviour. We want to be better than Matlab with such a
powerfull Language as Python, don't we?

4) Last but not least I need a simple example how to use matrixobjects
writing own extension.

Suppose I have a method for comuting the first eigenvalue of a sqaure
matrix:

int n=1000;
double **a=....;
double lambda1=first_eigenvalue(a, n);

What is the standard way of using that in python?
a=Matrix_(n,n)
lambda1=first_eigenvalue(a)

I mean how should I produce (access) a double** having a matrixobject?
(this is not a question about how to write an extenssion, just about
how to use the matrixobject :-))

5) Are the handouts of the mentionned talk about the matrix module available
for    all? Would be nice!

Tom

=================
MATRIX-SIG  - SIG on Matrix Math for Python

send messages to: matrix-sig@python.org