Iterators, reshape and the mathematical consistency (Kronecker)

[iagobaapellaniz]

I'll try to explain why I don't like the way reshape(), eachindex() and other iterative functions iterate a multidimensional array

With Kronecker product:

kron([1, 0],[0, 1, 0])
6-element Array{Int64,1}:
 0
 1
 0
 0
 0
 0

So to say, the 1st element of first matrix times the 2nd element of the second is the only non-zero element. See how reshape changes all this:

v = fill(0,2,3)
v[1,2] = 1

Again, it seems that the 1st element of 2 dimensional representation and the 2nd element of a 3 dimensional representation es the only element different from zero. But,

reshape(v,6)
6-element Array{Int64,1}:
 0
 0
 1
 0
 0
 0

is not consistent with how the Kronecker product does.

Probably, there are some reasons to do it in that way. Meanwhile, I've to redefine a function similar to reshape() to work with.

[iagobaapellaniz]

Just to sketch all this...

[johnmyleswhite]

I feel like you'd benefit from reading about column-major and row-major arrays. After reading about the topic more, I suspect you'll feel that this issue shouldn't be kept open.

[timholy]

"Preferred for you" has one big problem: it's not cache-efficient. Julia (and other languages that use LAPACK for their matrix operations) stores matrices in column-major order. Traversing in the order you propose would be a tenfold performance hit, unless we switched the internal representation to row-major order (which is how C and other languages do it).

I'd say the borked one is kron: it violates the "first dimension is fastest" maxim. Is matlab the source of this horrible function? (Spoken as a former matlab-enthusiast, but somehow I always used repmat and bsxfun instead of kron.)

[jiahao]

I'd say the borked one is kron: it violates the "first dimension is fastest" maxim. Is matlab the source of this horrible function?

Yes

[iagobaapellaniz]

Thanks! I lead to the same conclusion that first index should be first on iterate. The issue is now closed. I'm using the software for quantum systems and I will readecuate it for my goals.

But now I think this is no longer an open issue. For now obvious reasons.

[iagobaapellaniz]

I want to mention that Kronecker Product definition itself was what inspired the Matlab function kron. Cheers

[timholy]

Yeah, good point, it's implicit in the definition. It's interesting that mathematics has been so inconsistent. Numbering for matrices is (i,j), which is different from Cartesian (x,y) because we write i as vertical. Likewise, given how matrices are numbered, it seems that the Kronecker product should be defined with A and B swapped, because in A⊗B, A comes first and so its rows should take precedence.

Or better, the abstract properties of the Kronecker product suggest that its definition is correct. So it's really just evidence that the numbering of matrices opposite to that of Cartesian order was a big mistake. One that was made in probably the 17th century (plus or minus a few centuries), but we're stuck with it today.

[El-Zurdo-Sordo]

@timholy Given that, would it be so awful to have array types in Julia that are optionally row-major? (and optionally 0-based)? (keeping the column-major / 1 based defaults). For many data structures, say vectors of vectors, or vectors of vectors of vectors, they are naturally "row-major" as far as cache-efficiency, etc. It would also allow for easier dealing with C/C++/Java/Python etc. libraries.

[timholy]

It's been discussed many times before (search old issues and the mailing lists). Best place to start is trying it out by creating a package (and search first, there may already be one).

[eschnett]

See https://github.com/eschnett/FlexibleArrays.jl for arrays with arbitrary lower bounds. You would write

typealias Arr3d_lb0 FlexArray(0, 0, 0)
a3 = Arr3d_lb0{Float64}(9, 9, 9)

to define a type and create a 10x10x10 array with indices starting at 0.

Adding row-major arrays to this implementation would be straightforward.

[iagobaapellaniz]

Thank you so much. It's very inspiring all this.

[mcabbott]

Here is perhaps a clearer example. The following is is, I think, the natural definition of an outer product of two real vectors. The result is a 2-array with both of the given vectors' indices, in the same order:

julia> a = Real[1, 0.01]; b = 1:3;

julia> a * b'
2×3 Array{Real,2}:
 1     2     3
 0.01  0.02  0.03

julia> ans[2,3] == a[2] * b[3] # etc
true

But kron wants to return a 1-array, and instead of returning a reshape of this, it returns a reshape of the transpose of this:

julia> vec(a * b')
6-element Array{Real,1}:
 1
 0.01
 2
 0.02
 3
 0.03

julia> kron(b, a) == ans
true

julia> kron(a, b)
6-element Array{Any,1}:
 1
 2
 3
 0.01
 0.02
 0.03

julia> reshape(kron(a, b), 3,2)
3×2 Array{Any,2}:
 1  0.01
 2  0.02
 3  0.03

I suppose this would be the natural thing to do in a row-major language. It also matches Matlab's convention, despite Matlab being column-major. And it matches numpy's convention.