triangular matrix fill

classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

triangular matrix fill

Tom Waite
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

_______________________________________________
Numpy-discussion mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: triangular matrix fill

Charles R Harris


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



_______________________________________________
Numpy-discussion mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: triangular matrix fill

Robert Kern-2
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
_______________________________________________
Numpy-discussion mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: 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 <[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
_______________________________________________
Numpy-discussion mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/numpy-discussion


_______________________________________________
Numpy-discussion mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: triangular matrix fill

Charles R Harris


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 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


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



_______________________________________________
Numpy-discussion mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/numpy-discussion