www.digitalmars.com         C & C++   DMDScript  

digitalmars.D.learn - Multiplying transposed matrices in mir

reply p.shkadzko <p.shkadzko gmail.com> writes:
I'd like to calculate XX^T where X is some [m x n] matrix.

// create a 3 x 3 matrix
Slice!(double*, 2LU) a = [2.1, 1.0, 3.2, 4.5, 2.4, 3.3, 1.5, 0, 
2.1].sliced(3, 3);
auto b = a * a.transposed; // error

Looks like it is not possible due to "incompatible types for (a) 
* (transposed(a)): Slice!(double*, 2LU, cast(mir_slice_kind)2) 
and Slice!(double*, 2LU, cast(mir_slice_kind)0)"

I'd like to understand why and how should this operation be 
performed in mir.
Also, what does the last number "0" or "2" means in the type 
definition "Slice!(double*, 2LU, cast(mir_slice_kind)0)"?
Apr 19 2020
parent reply jmh530 <john.michael.hall gmail.com> writes:
On Sunday, 19 April 2020 at 17:07:36 UTC, p.shkadzko wrote:

 I'd like to calculate XX^T where X is some [m x n] matrix.

 // create a 3 x 3 matrix
 Slice!(double*, 2LU) a = [2.1, 1.0, 3.2, 4.5, 2.4, 3.3, 1.5, 0, 
 2.1].sliced(3, 3);
 auto b = a * a.transposed; // error

 Looks like it is not possible due to "incompatible types for 
 (a) * (transposed(a)): Slice!(double*, 2LU, 
 cast(mir_slice_kind)2) and Slice!(double*, 2LU, 
 cast(mir_slice_kind)0)"

 I'd like to understand why and how should this operation be 
 performed in mir.
 Also, what does the last number "0" or "2" means in the type 
 definition "Slice!(double*, 2LU, cast(mir_slice_kind)0)"?


2 is Contiguous, 0 is Universal, 1 is Canonical. To this day, I 
don’t have the greatest understanding of the difference.

Try the mtimes function in lubeck.
Apr 19 2020
parent reply p.shkadzko <p.shkadzko gmail.com> writes:
On Sunday, 19 April 2020 at 17:22:12 UTC, jmh530 wrote:

 On Sunday, 19 April 2020 at 17:07:36 UTC, p.shkadzko wrote:

 I'd like to calculate XX^T where X is some [m x n] matrix.

 // create a 3 x 3 matrix
 Slice!(double*, 2LU) a = [2.1, 1.0, 3.2, 4.5, 2.4, 3.3, 1.5, 
 0, 2.1].sliced(3, 3);
 auto b = a * a.transposed; // error

 Looks like it is not possible due to "incompatible types for 
 (a) * (transposed(a)): Slice!(double*, 2LU, 
 cast(mir_slice_kind)2) and Slice!(double*, 2LU, 
 cast(mir_slice_kind)0)"

 I'd like to understand why and how should this operation be 
 performed in mir.
 Also, what does the last number "0" or "2" means in the type 
 definition "Slice!(double*, 2LU, cast(mir_slice_kind)0)"?


 2 is Contiguous, 0 is Universal, 1 is Canonical. To this day, I 
 don’t have the greatest understanding of the difference.

 Try the mtimes function in lubeck.


Ah, I see. There are docs on internal representations of Slices 
but nothing about the rationale. It would be nice to have them 
since it is pretty much the core of Slice.

"a.mtimes(a.transposed);" works but the results are different 
from what NumPy gives.

For example:

a = np.array([[1, 2], [3, 4]])


Slice!(int*, 2LU) a = [1, 2, 3, 4].sliced(2,2);
writeln(a.mtimes(a.transposed)); // [[5, 11], [11, 25]]


So, lubeck mtimes is equivalent to NumPy "a.dot(a.transpose())".
Apr 19 2020
parent reply jmh530 <john.michael.hall gmail.com> writes:
On Sunday, 19 April 2020 at 17:55:06 UTC, p.shkadzko wrote:

 snip

 So, lubeck mtimes is equivalent to NumPy "a.dot(a.transpose())".


There are elementwise operation on two matrices of the same size 
and then there is matrix multiplication. Two different things. 
You had initially said using an mxn matrix to do the calculation. 
Elementwise multiplication only works for matrices of the same 
size, which is only true in your transpose case when they are 
square. The mtimes function is like dot or   in python and does 
real matrix multiplication, which works for generic mxn matrices. 
If you want elementwise multiplication of a square matrix and 
it’s transpose in mir, then I believe you need to call 
assumeContiguous after transposed.
Apr 19 2020
parent reply p.shkadzko <p.shkadzko gmail.com> writes:
On Sunday, 19 April 2020 at 18:59:00 UTC, jmh530 wrote:

 On Sunday, 19 April 2020 at 17:55:06 UTC, p.shkadzko wrote:

 snip

 So, lubeck mtimes is equivalent to NumPy 
 "a.dot(a.transpose())".


 There are elementwise operation on two matrices of the same 
 size and then there is matrix multiplication. Two different 
 things. You had initially said using an mxn matrix to do the 
 calculation. Elementwise multiplication only works for matrices 
 of the same size, which is only true in your transpose case 
 when they are square. The mtimes function is like dot or   in 
 python and does real matrix multiplication, which works for 
 generic mxn matrices. If you want elementwise multiplication of 
 a square matrix and it’s transpose in mir, then I believe you 
 need to call assumeContiguous after transposed.


"assumeContiguous" that's what I was looking for. Thanks!
Apr 19 2020
parent reply p.shkadzko <p.shkadzko gmail.com> writes:
On Sunday, 19 April 2020 at 19:13:14 UTC, p.shkadzko wrote:

 On Sunday, 19 April 2020 at 18:59:00 UTC, jmh530 wrote:

 On Sunday, 19 April 2020 at 17:55:06 UTC, p.shkadzko wrote:

 snip

 So, lubeck mtimes is equivalent to NumPy 
 "a.dot(a.transpose())".


 There are elementwise operation on two matrices of the same 
 size and then there is matrix multiplication. Two different 
 things. You had initially said using an mxn matrix to do the 
 calculation. Elementwise multiplication only works for 
 matrices of the same size, which is only true in your 
 transpose case when they are square. The mtimes function is 
 like dot or   in python and does real matrix multiplication, 
 which works for generic mxn matrices. If you want elementwise 
 multiplication of a square matrix and it’s transpose in mir, 
 then I believe you need to call assumeContiguous after 
 transposed.


 "assumeContiguous" that's what I was looking for. Thanks!


well no, "assumeContiguous" reverts the results of the 
"transposed" and it's "a * a".
I would expect it to stay transposed as NumPy does "assert 
np.all(np.ascontiguous(a.T) == a.T)".
Apr 19 2020
parent reply jmh530 <john.michael.hall gmail.com> writes:
On Sunday, 19 April 2020 at 19:20:28 UTC, p.shkadzko wrote:

 [snip]
 well no, "assumeContiguous" reverts the results of the 
 "transposed" and it's "a * a".
 I would expect it to stay transposed as NumPy does "assert 
 np.all(np.ascontiguous(a.T) == a.T)".


Ah, you're right. I use it in other places where it hasn't been 
an issue.

I can do it with an allocation (below) using the built-in syntax, 
but not sure how do-able it is without an allocation (Ilya would 
know better than me).

/+dub.sdl:
dependency "lubeck" version="~>1.1.7"
dependency "mir-algorithm" version="~>3.7.28"
+/
import mir.ndslice;
import lubeck;

void main() {
     auto a = [2.1, 1.0, 3.2, 4.5, 2.4, 3.3, 1.5, 0, 
2.1].sliced(3, 3);
     auto b = a * a.transposed.slice;
}
Apr 19 2020
parent reply p.shkadzko <p.shkadzko gmail.com> writes:
On Sunday, 19 April 2020 at 20:06:23 UTC, jmh530 wrote:

 On Sunday, 19 April 2020 at 19:20:28 UTC, p.shkadzko wrote:

 [snip]
 well no, "assumeContiguous" reverts the results of the 
 "transposed" and it's "a * a".
 I would expect it to stay transposed as NumPy does "assert 
 np.all(np.ascontiguous(a.T) == a.T)".


 Ah, you're right. I use it in other places where it hasn't been 
 an issue.

 I can do it with an allocation (below) using the built-in 
 syntax, but not sure how do-able it is without an allocation 
 (Ilya would know better than me).

 /+dub.sdl:
 dependency "lubeck" version="~>1.1.7"
 dependency "mir-algorithm" version="~>3.7.28"
 +/
 import mir.ndslice;
 import lubeck;

 void main() {
     auto a = [2.1, 1.0, 3.2, 4.5, 2.4, 3.3, 1.5, 0, 
 2.1].sliced(3, 3);
     auto b = a * a.transposed.slice;
 }


Thanks. I somehow missed the whole point of "a * a.transposed" 
not working because "a.transposed" is not allocated.
Apr 19 2020
next sibling parent reply jmh530 <john.michael.hall gmail.com> writes:
On Sunday, 19 April 2020 at 20:29:54 UTC, p.shkadzko wrote:

 [snip]

 Thanks. I somehow missed the whole point of "a * a.transposed" 
 not working because "a.transposed" is not allocated.


a.transposed is just a view of the original matrix. Even when I 
tried to do a raw for loop I ran into issues because modifying 
the original a in any way caused all the calculations to be wrong.

Honestly, it's kind of rare that I would do an element-wise 
multiplication of a matrix and its transpose.
Apr 19 2020
parent reply p.shkadzko <p.shkadzko gmail.com> writes:
On Sunday, 19 April 2020 at 21:27:43 UTC, jmh530 wrote:

 On Sunday, 19 April 2020 at 20:29:54 UTC, p.shkadzko wrote:

 [snip]

 Thanks. I somehow missed the whole point of "a * a.transposed" 
 not working because "a.transposed" is not allocated.


 a.transposed is just a view of the original matrix. Even when I 
 tried to do a raw for loop I ran into issues because modifying 
 the original a in any way caused all the calculations to be 
 wrong.

 Honestly, it's kind of rare that I would do an element-wise 
 multiplication of a matrix and its transpose.


It is. I was trying to calculate the covariance matrix of some 
dataset X which would be XX^T.
Apr 20 2020
parent jmh530 <john.michael.hall gmail.com> writes:
On Monday, 20 April 2020 at 19:06:53 UTC, p.shkadzko wrote:

 [snip]
 It is. I was trying to calculate the covariance matrix of some 
 dataset X which would be XX^T.


Incorrect. The covariance matrix is calculated with matrix 
multiplication, not element-wise  multiplication. For instance, I 
often work with time series data that is TXN where T > N. 
Couldn't do a calculation with element-wise multiplication in 
that case.

Try using Lubeck's covariance function or checking your results 
with the covariance function in other languages.
Apr 20 2020
prev sibling parent reply 9il <ilyayaroshenko gmail.com> writes:
On Sunday, 19 April 2020 at 20:29:54 UTC, p.shkadzko wrote:

 On Sunday, 19 April 2020 at 20:06:23 UTC, jmh530 wrote:

 On Sunday, 19 April 2020 at 19:20:28 UTC, p.shkadzko wrote:

 [...]


 Ah, you're right. I use it in other places where it hasn't 
 been an issue.

 I can do it with an allocation (below) using the built-in 
 syntax, but not sure how do-able it is without an allocation 
 (Ilya would know better than me).

 /+dub.sdl:
 dependency "lubeck" version="~>1.1.7"
 dependency "mir-algorithm" version="~>3.7.28"
 +/
 import mir.ndslice;
 import lubeck;

 void main() {
     auto a = [2.1, 1.0, 3.2, 4.5, 2.4, 3.3, 1.5, 0, 
 2.1].sliced(3, 3);
     auto b = a * a.transposed.slice;
 }


 Thanks. I somehow missed the whole point of "a * a.transposed" 
 not working because "a.transposed" is not allocated.


In the same time, the SliceKind isn't matter for assignment 
operations:

auto b = a.slice; // copy a to b
b[] *= a.transposed; // works well
Apr 19 2020
parent reply 9il <ilyayaroshenko gmail.com> writes:
On Monday, 20 April 2020 at 02:42:33 UTC, 9il wrote:

 On Sunday, 19 April 2020 at 20:29:54 UTC, p.shkadzko wrote:

 On Sunday, 19 April 2020 at 20:06:23 UTC, jmh530 wrote:

 [...]


 Thanks. I somehow missed the whole point of "a * a.transposed" 
 not working because "a.transposed" is not allocated.


 In the same time, the SliceKind isn't matter for assignment 
 operations:

 auto b = a.slice; // copy a to b
 b[] *= a.transposed; // works well


BTW for the following operation

auto b = a * a.transposed.slice;

`b` isn't allocated as well because `*` is lazy.


 auto b = a.slice; // copy a to b
 b[] *= a.transposed; // works well


So, the assignment operations are preferable anyway.
Apr 19 2020
parent p.shkadzko <p.shkadzko gmail.com> writes:
On Monday, 20 April 2020 at 02:50:29 UTC, 9il wrote:

 On Monday, 20 April 2020 at 02:42:33 UTC, 9il wrote:

 On Sunday, 19 April 2020 at 20:29:54 UTC, p.shkadzko wrote:

 On Sunday, 19 April 2020 at 20:06:23 UTC, jmh530 wrote:

 [...]


 Thanks. I somehow missed the whole point of "a * 
 a.transposed" not working because "a.transposed" is not 
 allocated.


 In the same time, the SliceKind isn't matter for assignment 
 operations:

 auto b = a.slice; // copy a to b
 b[] *= a.transposed; // works well


 BTW for the following operation

 auto b = a * a.transposed.slice;

 `b` isn't allocated as well because `*` is lazy.


 auto b = a.slice; // copy a to b
 b[] *= a.transposed; // works well


 So, the assignment operations are preferable anyway.


Interesting, thanks for the examples.
Apr 20 2020