

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]; } } Thanks Tom
_______________________________________________
Numpydiscussion mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/numpydiscussion


On Thu, May 22, 2008 at 7:19 PM, Tom Waite < [hidden email]> 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.
Chuck
_______________________________________________
Numpydiscussion mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/numpydiscussion


On Thu, May 22, 2008 at 9:07 PM, Charles R Harris
< [hidden email]> wrote:
>
> On Thu, May 22, 2008 at 7:19 PM, Tom Waite < [hidden email]> 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
_______________________________________________
Numpydiscussion mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/numpydiscussion


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 LevenbergMarquardt 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 < [hidden email]> wrote:
On Thu, May 22, 2008 at 9:07 PM, Charles R Harris
< [hidden email]> wrote:
>
> On Thu, May 22, 2008 at 7:19 PM, Tom Waite < [hidden email]> 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
_______________________________________________
Numpydiscussion mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/numpydiscussion
_______________________________________________
Numpydiscussion mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/numpydiscussion


On Mon, May 26, 2008 at 11:06 AM, Tom Waite < [hidden email]> wrote:
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 LevenbergMarquardt 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 I think your own extension code is the place to start. If something like this goes into scipy, then it should probably be part of a larger package containing other relevant optimizations.
Chuck
_______________________________________________
Numpydiscussion mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/numpydiscussion

