[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[0] <= 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[0]
    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
administrivia to: matrix-sig-request@python.org
=================