[Numpy-discussion] Large symmetrical matrix

[Anne Archibald]

2008/6/10 Simon Palmer <simon.palmer at gmail.com>:
> Hi I have a problem which involves the creation of a large square matrix
> which is zero across its diagonal and symmetrical about the diagonal i.e.
> m[i,j] = m[j,i] and m[i,i] = 0.  So, in fact, it is a large triangular
> matrix.  I was wondering whether there is any way of easily handling a
> matrix of this shape without either incurring a memory penalty or a whole
> whack of proprietary code?
>
> To get through this I have implemented a 1D array which has ((n-1)^2)/2
> elements inside a wrapper class which manpulates the arguments of array
> accessors with some arithmetic to return the approriate value.  To be honest
> I'd love to throw this away, but I haven't yet come across a feasible
> alternative.
>
> Any ideas?

Well, there's a risk of shooting yourself in the foot, but here's a
way to get at the elements quickly and conveniently: store the
coefficients in a 1D array (as you do already). Then if your matrix is
M by M, element (i,j) is at position (M-2)*j+i-1, assuming i<j. Of
course, you can simply use this expression every time you want to look
an element up, as you do now, but that's no fun. So now create a
"view" A of this array having offset -1, shape (M,M), and strides
(1,(M-2)). (This kind of weird view can be created with the ndarray
constructor.) Then A[i,j] addresses the right element. A[j,i] is of
course totally wrong. But as long as you avoid i<j, you can use all
numpy's magical indexing tricks. Unfortunately you can't really use
ufuncs, since they'll waste time operating on the useless lower half.

Personally, I would live with the factor-of-two memory hit.

Anne