矩阵块索引

Question

我使用的是Julia 1.1版。我正在研究可以由更小的矩阵构造的矩阵，例如泡利矩阵。我不清楚如何通过在Julia中使用一组较小的矩阵来有效地构造大矩阵，即直接将较小的矩阵写入某个索引位置。

Julias kron并不令人满意，因为我需要生成几个“大矩阵”才能得到最终结果。例如，我想创建类似这样的东西(这只是一个很小的例子)

sy = [[0 -im]; [im 0]]
M = [[0 sy adjoint(sy)]; 
     [adjoint(sy) 0 sy];
     [sy adjoint(sy) 0]]

可以通过执行两个kronecker产品将两个结果相加来实现此目的。然而，这将是一个巨大的浪费，特别是如果矩阵变得更大的话。

我也已经尝试过使用BlockArrays.jl包，但意识到它不能完全满足我的需求。

最后，我希望能够处理我的大矩阵的“矩阵块”，这样我就可以直接将构造矩阵分配到正确的位置，对于上面的示例，这将如下所示(我在这里没有使用循环来更清楚地说明我的观点)：

M[1, 2] = sy
M[1, 3] = adjoint(sy)
M[2, 1] = adjoint(sy)
M[2, 3] = sy
M[3, 1] = sy
M[3, 2] = adjoint(sy)

我意识到这意味着将我原来的大数组索引减少到类似数组“块索引”的东西。

我考虑过用视图来做这件事，我创建了一个SubArrays矩阵，然后我可以用矩阵块索引符号来处理这个矩阵。

S0 = view(M, 1:2, 1:2)
S1 = view(M, 1:2, 2:4)
S2 = view(M, 1:2, 4:6)
...

Viewmatrix = [[S0 S1 S2]; [S3 S4 S5]; [S6 S7 S8]]
Viewmatrix[1, 2] .= sy
Viewmatrix[1, 3] .= adjoint(sy)
...

现在我还不清楚如何去写这样一个视图矩阵，或者这是不是一个可行的解决问题的方法。如果有更好的方法来解决这个问题，我很想知道。

Answer 1

BlockArrays.jl不仅支持2×2分块数组，尽管他们过去在文档中只使用它们。您可以轻松创建所需的3×3分块6×6数组，如下所示：

M = BlockArray(fill(0im, 6, 6), [2, 2, 2], [2, 2, 2])
M[Block(1, 2)] = sy
M[Block(1, 3)] = adjoint(sy)
M[Block(2, 1)] = adjoint(sy)
M[Block(2, 3)] = sy
M[Block(3, 1)] = sy
M[Block(3, 2)] = adjoint(sy)

julia> M
3×3-blocked 6×6 BlockArray{Complex{Int64},2}:
 0+0im  0+0im  │  0+0im  0-1im  │  0+0im  0-1im
 0+0im  0+0im  │  0+1im  0+0im  │  0+1im  0+0im
 ──────────────┼────────────────┼──────────────
 0+0im  0-1im  │  0+0im  0+0im  │  0+0im  0-1im
 0+1im  0+0im  │  0+0im  0+0im  │  0+1im  0+0im
 ──────────────┼────────────────┼──────────────
 0+0im  0-1im  │  0+0im  0-1im  │  0+0im  0+0im
 0+1im  0+0im  │  0+1im  0+0im  │  0+0im  0+0im

但要小心:这些块是按引用存储的。因此，如果您后来修改了sy，则包含它的所有块也将更改，反之亦然。如果你想避免这种情况，可以使用广播赋值(用.=代替=)。

如果您的问题实际上与示例一样简单，并且更密集，那么使用mortar函数将可用块“粘合在一起”可能会更简单：

julia> mortar(reshape([z, sy, sy', sy', z, sy, sy, sy', z], (3, 3)))
3×3-blocked 6×6 BlockArray{Complex{Int64},2,Array{AbstractArray{Complex{Int64},2},2},BlockArrays.BlockSizes{2,Array{Int64,1}}}:
 0+0im  0+0im  │  0+0im  0-1im  │  0+0im  0-1im
 0+0im  0+0im  │  0+1im  0+0im  │  0+1im  0+0im
 ──────────────┼────────────────┼──────────────
 0+0im  0-1im  │  0+0im  0+0im  │  0+0im  0-1im
 0+1im  0+0im  │  0+0im  0+0im  │  0+1im  0+0im
 ──────────────┼────────────────┼──────────────
 0+0im  0-1im  │  0+0im  0-1im  │  0+0im  0+0im
 0+1im  0+0im  │  0+1im  0+0im  │  0+0im  0+0im

尽管这在内部使用抽象类型，而不是提升分配的数组。