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)"?

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.

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]])
a * a.transpose() # [[1, 6], [6, 16]]

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())".

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.

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!

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)".

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

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.

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.

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