[Numpy-discussion] triangular matrix fill

[Tom Waite]

Thanks Robert and Chuck.

The matrix calculation is a bottleneck computation and that is my interest
in only doing the lower triangular part of the outer product. This is part
of a nonlinear brainwarp for the NIPY project. For a typical anatomical MRI
volume size this is a 768 x 768 matrix that is inside three nested loops
that walk the volume. This matrix is actually the curvature matrix
computation and is "inspired" by the Marquardt routine (mrqcof) in Numerical
Recipes. Having this lower triangular outer product in numpy could be of
value for scipy.optimize down the road if Levenberg-Marquardt is added. Is
this something I should post to be added to numpy or write my own extension
code as this might be too specialized?

Tom


On Thu, May 22, 2008 at 7:13 PM, Robert Kern <robert.kern at gmail.com> wrote:

> On Thu, May 22, 2008 at 9:07 PM, Charles R Harris
> <charlesr.harris at gmail.com> wrote:
> >
> > On Thu, May 22, 2008 at 7:19 PM, Tom Waite <twaite at berkeley.edu> wrote:
> >>
> >> I have a question on filling a lower triangular matrix using numpy. This
> >> is essentially having two loops and the inner loop upper limit is the
> >> outer loop current index. In the inner loop I have a vector being
> >> multiplied by a constant set in the outer loop. For a matrix N*N in
> size,
> >> the C the code is:
> >>
> >> for(i = 0; i < N; ++i){
> >>     for(j = 0; j < i; ++j){
> >>         Matrix[i*N + j] = V1[i] * V2[j];
> >>     }
> >> }
> >>
> >
> > You can use numpy.outer(V1,V2) and just ignore everything  on and above
> the
> > diagonal.
> >
> > In [1]: x = arange(3)
> >
> > In [2]: y = arange(3,6)
> >
> > In [3]: outer(x,y)
> > Out[3]:
> > array([[ 0,  0,  0],
> >        [ 3,  4,  5],
> >        [ 6,  8, 10]])
> >
> > You can mask the upper part if you want:
> >
> > In [16]: outer(x,y)*fromfunction(lambda i,j: i>j, (3,3))
> > Out[16]:
> > array([[0, 0, 0],
> >        [3, 0, 0],
> >        [6, 8, 0]])
> >
> >  Or you could use fromfunction directly.
>
> Or numpy.tril().
>
> --
> Robert Kern
>
> "I have come to believe that the whole world is an enigma, a harmless
> enigma that is made terrible by our own mad attempt to interpret it as
> though it had an underlying truth."
>  -- Umberto Eco
> _______________________________________________
> Numpy-discussion mailing list
> Numpy-discussion at scipy.org
> http://projects.scipy.org/mailman/listinfo/numpy-discussion
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20080526/4a198496/attachment.html>