[PYTHON MATRIX-SIG] max_element array method

[Konrad HINSEN]

> I don't know whether this has been discussed before because I have yet
> to find the matrix-sig archive. I'm currently implementing some clustering 
> algorithms and came across two features that would be very helpful to
> me:
> 
> 1. A space-efficient implementation of a symmetric matrix object, and

That is not so difficult. The most straightforward way is to define
symmetric matrices as a Python class and use a one-dimensional array
object as the data representation. The most reasonable arrangement is
the packed format used by LAPACK, and of course by using it you have
immediate access to all the linear algebra functions from that
library. All the elementwise operations can be mapped directly on the
internal representation. Indexing is straightforward. Matrix
multiplication has to be rewritten, and is a bit tricky, since the
product of two symmetric matrices is in general not a symmetric
matrix.

But if you want full compatibility with ordinary arrays, the
implementation gets complicated. Slices, for example, are messy.  So
are reduction operations, for example. But these operations don't make
much sense for symmetric matrices anyway. So the best solution is not
to pretend that symmetric matrices are a special case of general
arrays and implement only a meaningful subset of operations.

> 2. a max_elem() and min_elem() method that would return the tuple
>    (max_elem, index). For example:

Right, there is very little support now for efficient determination
of specific indices. I had a similar problem recently with my
class for function interpolation: finding the index of the largest
element smaller than a given number in a sorted array.

> Currently, I'm using a Python class that does just that. But a C-level 
> implementation would be probably faster and more space-efficient.

I don't know how you do it in Python, but if you are doing the whole
search in Python, consider the following alternative:

  def max_index(a):
    max = maximum.reduce(a)
    index = compress(equal(a, max), arange(a.shape[0]))
    return max, index

This version is only for rank-1 arrays; a general version is possible,
but messy. It is still inefficient compared to C, but there is no
loop in Python, which probably makes it faster for long arrays.
And if the maximum occurs more than once, you even get all the
indices.

-------------------------------------------------------------------------------
Konrad Hinsen                     | E-Mail: hinsenk@ere.umontreal.ca
Departement de chimie             | Tel.: +1-514-343-6111 ext. 3953
Universite de Montreal            | Fax:  +1-514-343-7586
C.P. 6128, succ. Centre-Ville     | Deutsch/Esperanto/English/Nederlands/
Montreal (QC) H3C 3J7             | Francais (phase experimentale)
-------------------------------------------------------------------------------

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

send messages to: matrix-sig@python.org
administrivia to: matrix-sig-request@python.org
=================