RE: Octave 3.0

[Lippert, Ross A.]

(this post is long than I thought it would be.
contents:
  re: zero-copy reshape (and its myriad uses)
  re: sparse matrices
  re: matlab compatability
  re: MPI and PVM
)

re: zero-copy reshape
>reshape as currently implemented requires a matrix copy,
>but that would be easy to fix.
I agree that a copy is inconvenient.  Usually it is not an issue
since the other ops are more expensive than the copy.

But, one thing that I think could be really cool on this topic
is to abstract the block of doubles allocated to a matrix/vector
from the way the block is indexed, supporting multiple
views for the same data, e.g.:
 B = A'
or
 B = reshape(A,prod(size(A)),1)
would allocate no new memory for B, just create a new 'view'.  A new
block only gets allocated when one writes to an entry of B.  This
kind of view business should be compatable with the basic BLAS ops
since they allow transpose bits and strides to be set appropriately.

I noticed that GSL supports this trick, to the extent that
 v = A(1,:)
is just another view.  I'm not sure if they handle
 B = A([1 3 4],:)
as a view, and I don't see how one could do a gemm on such a B
without forcing a tmp copy.

BTW I am not an advocate of GSL-ing octave.  I've looked at some
of the GSL sources, and I'm not as impressed with the implementation
as I am with the interface.

>There is something to be said for convenient data
>structures, especially if it is convenient to operate
>with slices.  The results are going to be more readable
IMHO I'd change 'especially' to 'only'.  Sometimes the operations
we think we want aren't the right ones, and when we figure out what
ops we want, we figure out better representations for our data.

>than tricky indexing and reshaping operations.  It's
>not clear to me that Matlab does the right thing --- not
>while squeeze/reshape are required to transform a
>slice using a matrix operation.
Definitely.  I think it is very hard to work out what one
really should have basic ops on multidimensional data.  One
ugly, but sensible, thing one should have is something like
 C = CONTRACT(A,B,[2 3],[4 1])
where C(a,b,c) = SUM_{xy} A(a,x,y)*B(y,b,c,x)
but that's kind of nasty.  On the other hand one also can do
 C = reshape( reshape(A,na,nx*ny)*reshape( reshape(B,ny*nb*nc,nx)' 
,nx*ny,nb*nc) ,na,nb,nc)
which is also kind of nasty, but at least it introduces no new ops
to the language.

And yes, you have to have the cock-eyed view of 
matrices like I have , so you are right that this is for the
sophisticated user and not for everyone.  Maybe being forced
to do things like this makes one sophisticated.  Maybe I'd just
too proud of myself for tricks like this.

If zero-copy reshape/transpose were implemented, and some support
for inline functions were around, then perhaps one could just
write .m scripts to support N-D arrays and sparse matrices.

re: sparsity
For sparse matrices, I agree that idxop is definitely for sophisticates
and wrapping it up might require making some sort of special 'struct'
with some kind of overloading.  I'm kind of partial to some special
OO struct type which has a special field which says 'I have functions
in me', and then when someone does A*B and the interpreter sees that one
of them is a struct, and instead of the usual 'I don't know what to do with
structs' error message, it can check for the OO flag and look in the struct
fields for 'times' to pull out the appropriate functions.  This is a
nice way to take advantage of the fact that functions are first-class
(or nearly first class) objects in octave, and I think it is better than
MATLABs "greedy overloading by directory name" crap.

function s = sparse(A)
  [i,j,a] = find(A);
  s.OO = true;
  s.i = i;  s.j = j; s.a = a;
  s.times = sparse_times;
  s.plus  = sparse_plus;
  ...
endfunction

BTW just one last thing on idxop.  The thing that really necessitates it
is the fact that += doesn't do 'the right thing' (which may not be right
to some people) when the LHS is multiply indexed:
 v = zeros(2,1);
 v([1 2 2]) += 5
currently sets v = [5 5] instead of [5 10].  A sparse multiply could then
be
 [i,j,a] = find(A);
 y = zeros(size(A,1),1)
 y(i) += a.*x(j)
I had this argument before with some people on this list (and I'm content tha
my position may not be the right one), and I am aware that it would require
some serious fiddling with octave sources to do it my way.  Then there is
the issue of not having a max= or a min= operator for the other idxops,
though one might co-opt &= and |= for these services.

re: MATLAB compatability
The bit about rewriting code on the web is a really good one, for
which I don't have a good answer.  MATLAB has allowed people to
write some obnoxious things, and for octave to support that code
it needs to support some of these obnoxious things.  Perhaps, it
can be done with little change to octave other than some mods to
the basic operations then wrapped up in .m's.  That would be nice.

My point is that octave presents an opportunity not just to be
compatable with MATLAB, but to go that extra mile and do things
smarter than MATLAB.  (BTW another personal example of this is that
octave includes the LAPACK sylvester equation solver, which, to my
shock, MATLAB doesn't, and I was solving sylvester equations a lot
at one point and probably will need to again in the future).  I'd
like to see new features put in because it is a good idea, whether
from MATLAB or elsewhere, with a MATLAB-like interface when
applicable.

re: MPI and PVM
>Has anyone looked at the PVM project for scilab?
Three years ago I wrote some PVM DLD functions for octave.  The mkoctfile
script has become way more handy and useful for these purposes (either
that or I have learned how to use mkoctfile).  I had a lot of fun creating
matrices on one octave process and passing them to another octave process.
It was a fun toy, but it was a hack.  If someone would like to unhack it
it is at
http://www.eskimo.com/~ripper/research/pvmoct.html

I had thought of doing the same for MPI, but that whole ARGV thing made
me stop.


-r
-----Original Message-----
From: Paul Kienzle [mailto:address@hidden
Sent: Saturday, March 01, 2003 7:32 AM
To: Lippert, Ross A.
Cc: octave-maintainers mailing list
Subject: Re: Octave 3.0


Lippert, Ross A. wrote:

>I think it is important to separate those items which are on this list
>for the sake of MATLAB compatibility and those items on this list which
>are for the general improvement of the package.  For example, MATLAB's
>N-dimensional arrays are a convenient storage class, but they don't
>interface with the mathematical operations very well.  Many times people
>think they want an N-dimensional array because they are dealing with,
>say a 3D mesh, but I think one quickly finds that it is more worthwhile
>to keep these values in a 1D vector and the does 'reshape' as necessary.
>
reshape as currently implemented requires a matrix copy,
but that would be easy to fix.

There is something to be said for convenient data
structures, especially if it is convenient to operate
with slices.  The results are going to be more readable
than tricky indexing and reshaping operations.  It's
not clear to me that Matlab does the right thing --- not
while squeeze/reshape are required to transform a
slice using a matrix operation.

>* Maybe I am just shooting my mouth off, but is there anyone out there who
>makes any serious use of N-d arrays in their MATLAB work?
>
>** On the other hand, if this is a move to further increase octave's .m
>compatability, then it is necessary, however awful.
>
>As for sparse matrices, I have gotten a lot of mileage out of a DLD function
>Paul Kienzle and I came up with called idxop which does the following:
>   y = idxop(idx,x,['sum'|'prod'|'max'|'min'])
>where idx is an index vector whose length is equal to the number of rows of x
>and the resulting y is essentially given by
>   y = zeros(max(idx),size(x,2))+[0|1|-inf|inf]
>   for i=1:size(x,1),
>      y(idx(i),:) = [plus|times|max|min] ( y(idx(i),:), x(i,:) )
>   endfor
>This single function allows for a variety of sparse operations to occur on
>regular vectors/matrices without the need for a sparse data structure.  The
>advantage of doing things in terms of sparse ops and not sparse data is that
>quite often the operations one wishes to perform are a mix of sparse and
>dense ops (e.g. you have a graph laplacian on a disconnected graph and you
>wish to form the dense laplacian for each graph component and take its
>eigenvalues).  It is my belief that in practise, one spends quite a bit
>of time doing 'sparse' and 'full' when using MATLAB sparse matrix functions
>just so they can do a sparse op on their data at some point in the compute.
>
>Again * and ** apply.
>
* applies until you can provide recipes for
common 3-D and sparse operations that
unsophisticated users can apply --- yes you
can program everything in machine code but
macro assemblers sure help a lot.  The octave
wiki is a good forum for doing this:
    http://www.scarymath.org/octave

** applies while there exists interesting
code on the web that you want to use
without having to rewrite.

>Private functions would be a nice way to clean up the ever growing octave
>namespace, whether MATLAB allows this currently or not.
>
Not enough I expect.  To get that you need something
like namespace support.

>MPI support would raise a lot of eyebrows.  It is something my PhD advisor
>was once trying to shoe-horn into a customized version of MATLAB, but with
>little support from the mathworks.  Such a feature would really differentiate
>octave.
>
Has anyone looked at the PVM project for scilab?