NumPy and Octave (qestion and discussion)

[Huaiyu Zhu]

Question 1 (call octave from python):

How to call octave scripts from Python?  Since both support C interface, is
this the way to go?  Is there anywhere where I can find specific examples?
What about default arguments and optional arguments?


Question 2 (modify NumPy):

I do like octave / matlab syntax for numerical computations, which is
simple, elegant, and resembles the formulas written on paper. For example:

x = [1,2,3,4]       # row vector (1x4 matrix)
y = x'+2            # col vector (4x1 matrix)
x*y                 # scalar (1x1 matrix)
y*x                 # 4x4 matrix
a = x(1:2)          # row vector (1x2 matrix)
b = y(3:4)          # col vector (2x1 matrix)
a*b                 # scalar
b*a                 # 2x2 matrix

>From what I know, NumPy has at least the following defect compared with octave:

(1) Don't distinguish row and col for  1-dim vector.
(2) Can't identify 1xn and nx1 matrices with vectors.
    Can't identify 1x1 matrices with scalars.
(3) Default multiplication is element-wise, instead of linear algebra.
(4) mystical expansion of 1xn * nx1 into matrix (against rule (3) above).
(5) No \ operator (for implicit solve of linear equations).

Can NumPy be modified to look similar to octave?  I bet there would be a lot
of people wanting similar feature.  Right now the above example is so much
more complicated in Python: 

from Numeric import *
x = [1,2,3,4]       # list
x = array(x)        # array (neither row nor col)
y = transpose(x)    # same array
print x*y           # array = x**2
print y*x           # array = x**2

x = x[NewAxis,:]    # row vector (1x4 matrix)
y = transpose(x)+2  # col vector (4x1 matrix)
print x*y           # 4x4 matrix
print y*x           # 4x4 matrix
print matrixmultiply(x, y)  # 1x1 matrix (not scalar)
print matrixmultiply(y, x)  # 4x4 matrix
print matrixmultiply(x, y)[0][0]    # scalar

a = x[0:2]          # 1x4 matrix (=x)
b = y[2:4]          # 2x1 vector (as expected)
a = transpose(transpose(x)[0:2]) # 1x2 matrix (as expected)
print a*b           # 2x2 matrix
print b*a           # 2x2 matrix
print matrixmultiply(a, b)  # 1x1 matrix (not scalar)
print matrixmultiply(b, a)  # 2x2 matrix
print matrixmultiply(a, b)[0][0]    # scalar

The above are only basic operations.  You can easily trip over things like
a=x[0,0:2]; b=y[2:4,0]; print a*b, and many more.

For a typical use of numerical computation in octave:

X = rand(5,3)            # 5x3 random matrix
y = rand(3,1)            # 3x1 vector
b = X\y                  # LMS solution of a linear equation y = X*b
b = (X'*X)\(X'*y)        # or written out in more details
b = inv(X'*X)*(X'*y)     # or in a less efficient form.

the corresponding Python notation is horrendous:
b = inv(matrixmultiply(transpose(X), X) * \
(matrixmultiply(transpose(X)*y[:,NewAxis])

Are there better ways to do numerical computations in python?
Thanks.

-- 
Huaiyu Zhu                               hzhu at knowledgetrack.com