[help] How does the transpose-function work and can this be optimized somehow?

[ClmnsRck]

This is my code for Multiple Linear Regression:

template <typename T, long M, long N>
void blitzMLR(){
    blitz::Array<T, 2> x(randVec<T>(M*(N+1)).data(), blitz::shape(M,N+1), blitz::duplicateData);
    blitz::Array<T, 1> y(randVec<T>(M).data(), blitz::shape(M), blitz::duplicateData);
    blitz::Array<T, 1> w(randVec<T>(N+1).data(), blitz::shape(N+1), blitz::duplicateData);
    for(long i = 0; i < M; ++i){
        x(0,i) = static_cast<T>(1);
    }
    {
        blitz::firstIndex a;
        blitz::secondIndex c;
        blitz::thirdIndex b;
        blitz::Array<T, 2> tmp1 = sum(x.transpose(1,0)(a,b)*x(b,c),b);
        blitz::Array<T, 2> tmp2; //TODO: Inverse of tmp1 **If somebody knows a way to do this, please let me know**
        blitz::Array<T, 2> tmp3 = sum(tmp2(a,b)*x.transpose(1,0)(b,c),b);
        w = sum(tmp3(a,b) * y(b), b);
    }
}

How do i do those to transposes? I would like to switch the first with the second dimension.

EDIT: IF you have anything to add on how to improve performance of this peace of code, anything inside the second scope-brackets is to be optimized ...

[lutorm]

If you're already working with index placeholders, just transpose the indices. The transpose of A(firstIndex, secondIndex) is A(secondIndex, firstIndex). But I don't understand why you have three names for firstIndex in your code, is that just a typo? If both a and b are firstIndex, then A(a,b) is the same as A(a,a), ie the diagonal vector.

On Mon, Apr 15, 2019 at 9:02 AM ClmnsRck notifications@github.com wrote:

This is my code for Multiple Linear Regression:

template <typename T, long M, long N> void blitzMLR(){ blitz::Array<T, 2> x(randVec(M(N+1)).data(), blitz::shape(M,N+1), blitz::duplicateData); blitz::Array<T, 1> y(randVec(M).data(), blitz::shape(M), blitz::duplicateData); blitz::Array<T, 1> w(randVec(N+1).data(), blitz::shape(N+1), blitz::duplicateData); for(long i = 0; i < M; ++i){ x(0,i) = static_cast(1); } { blitz::firstIndex a; blitz::firstIndex c; blitz::firstIndex b; blitz::Array<T, 2> tmp1 = sum(x.transpose(1,0)(a,b)x(b,c),b); blitz::Array<T, 2> tmp2; //TODO: Inverse of tmp1 If somebody knows a way to do this, please let me know blitz::Array<T, 2> tmp3 = sum(tmp2(a,b)x.transpose(1,0)(b,c),b); w = sum(tmp3(a,b) y(b), b); } }

How do i do those to transposes? I would like to switch the first with the second dimension.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/blitzpp/blitz/issues/114, or mute the thread https://github.com/notifications/unsubscribe-auth/AFK8GNr0YWe3OlTxA3-KuUL9dEgQFJM9ks5vhMy1gaJpZM4cwrD- .

[ClmnsRck]

Nope, that is a artifact from copy and pasting, i will fix that.

But about the inverse, i will probably have to implement that by myself, won't i?

[ClmnsRck]

And could i just say

sum(A(b,a)*B(b,c), b) //matrix-product of A^T * B

?

[lutorm]

There's a limitation on the partial reductions that say they can only reduce the last dimension, see 3.14 in the blitz user's guide. So while you'd want to do

C(a,c) = sum(A(a,b)*B(b,c),b)

to sum over the intermediate dimension, you actually have to swap the last two dimensions and do

C(a,b) = sum(A(a,c)*B(c,b),c)

(basically, the way to think about this is that the "output" dimensions of the expression has to map onto the dimensions of the thing you are assigning it to. Since the target array has the first two indices, those are the ones that have to remain in the final expression, after reductions, too.)

Hope that helps,

On Mon, Apr 15, 2019 at 11:53 AM ClmnsRck notifications@github.com wrote:

And could i just say

sum(A(b,a)B(b,c), b) //matrix-product of A^T B

?

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/blitzpp/blitz/issues/114#issuecomment-483432908, or mute the thread https://github.com/notifications/unsubscribe-auth/AFK8GPCWS7GkJ665D5kq53GPM1yjyqleks5vhPTrgaJpZM4cwrD- .

[lutorm]

Oh I just remembered that you want to use the transpose of one matrix in the product. Then yeah, you just transpose the indices for that Array.

On Mon, Apr 15, 2019 at 1:54 PM Patrik Jonsson code@familjenjonsson.org wrote:

There's a limitation on the partial reductions that say they can only reduce the last dimension, see 3.14 in the blitz user's guide. So while you'd want to do

C(a,c) = sum(A(a,b)*B(b,c),b)

to sum over the intermediate dimension, you actually have to swap the last two dimensions and do

C(a,b) = sum(A(a,c)*B(c,b),c)

(basically, the way to think about this is that the "output" dimensions of the expression has to map onto the dimensions of the thing you are assigning it to. Since the target array has the first two indices, those are the ones that have to remain in the final expression, after reductions, too.)

Hope that helps,

On Mon, Apr 15, 2019 at 11:53 AM ClmnsRck notifications@github.com wrote:

And could i just say

sum(A(b,a)B(b,c), b) //matrix-product of A^T B

?

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/blitzpp/blitz/issues/114#issuecomment-483432908, or mute the thread https://github.com/notifications/unsubscribe-auth/AFK8GPCWS7GkJ665D5kq53GPM1yjyqleks5vhPTrgaJpZM4cwrD- .

[ClmnsRck]

So which of my matrix-products do i have to change then? The one where the transpose is on the first? because this one cant be transposed just by changing the indices?

And the other thing is, do i have to write it like this: C(a,c) = sum(A(a,b)*B(b,c),b); or in other words, is the index on the lhs needed?