如何在 Julia 中以对角线形式连接矩阵

Question

我有一个矩阵A，我想构建一个块对角矩阵M as

A 0 ⋯ 0
0 A ⋯ 0
⋮ ⋮ ⋱ ⋮
0 0 ⋯ A

其中0s 是与A 大小相同且填充零的矩阵。有没有方便的方法来做到这一点？

换句话说，有没有相当于 Python 的 linalg.block_diag（来自 scipy）？

Answer 1

您可以使用 SparseArrays.jl 标准库中的 blockdiag。例如：

julia> using SparseArrays

julia> A = sparse(rand(3,2)) # A must be sparse for blockdiag
3×2 SparseMatrixCSC{Float64, Int64} with 6 stored entries:
 0.465226  0.473656
 0.230678  0.391924
 0.928409  0.772551

julia> M = blockdiag(A, A, A)
9×6 SparseMatrixCSC{Float64, Int64} with 18 stored entries:
 0.465226  0.473656   ⋅         ⋅         ⋅         ⋅
 0.230678  0.391924   ⋅         ⋅         ⋅         ⋅
 0.928409  0.772551   ⋅         ⋅         ⋅         ⋅
  ⋅         ⋅        0.465226  0.473656   ⋅         ⋅
  ⋅         ⋅        0.230678  0.391924   ⋅         ⋅
  ⋅         ⋅        0.928409  0.772551   ⋅         ⋅
  ⋅         ⋅         ⋅         ⋅        0.465226  0.473656
  ⋅         ⋅         ⋅         ⋅        0.230678  0.391924
  ⋅         ⋅         ⋅         ⋅        0.928409  0.772551

如果你真的需要 M 密集，你可以把它转换回来

julia> Matrix(M)
9×6 Matrix{Float64}:
 0.465226  0.473656  0.0       0.0       0.0       0.0
 0.230678  0.391924  0.0       0.0       0.0       0.0
 0.928409  0.772551  0.0       0.0       0.0       0.0
 0.0       0.0       0.465226  0.473656  0.0       0.0
 0.0       0.0       0.230678  0.391924  0.0       0.0
 0.0       0.0       0.928409  0.772551  0.0       0.0
 0.0       0.0       0.0       0.0       0.465226  0.473656
 0.0       0.0       0.0       0.0       0.230678  0.391924
 0.0       0.0       0.0       0.0       0.928409  0.772551

Answer 2

您可以使用cat 执行此操作，给它两个维度：

julia> A = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> cat(A, A', reverse(A); dims=(1,2))
6×6 Matrix{Int64}:
 1  2  0  0  0  0
 3  4  0  0  0  0
 0  0  1  3  0  0
 0  0  2  4  0  0
 0  0  0  0  4  3
 0  0  0  0  2  1

Answer 3

如果您不想为矩阵和 0 的副本分配额外空间，可以使用我编写的以下结构，它模拟 block-diagonal matrix。

struct BlockDiagonalMatrix{T} <: AbstractMatrix{T}
    matrix::AbstractMatrix{T}
    repeats::Integer

    BlockDiagonalMatrix(mat::AbstractMatrix{T}, repeats) where {T} = new{T}(mat, repeats)
end

Base.IndexStyle(::Type{BlockDiagonalMatrix}) = IndexLinear()
Base.size(bdm::BlockDiagonalMatrix) = size(bdm.matrix) .* bdm.repeats
function Base.getindex(bdm::BlockDiagonalMatrix{T}, i) where {T}
    bdm_height = size(bdm, 1)

    i -= 1  # Math is slightly simpler when indexing is 0-based
    i_col, i_row = divrem(i, bdm_height)

    block_height, block_width = size(bdm.matrix)
    block_col = fld(i_col, block_width)
    block_row = fld(i_row, block_height)

    if block_row == block_col
        inner_mat_col = 1 + (i_col % block_width)  # convert back to 1-based
        inner_mat_row = 1 + (i_row % block_height)
        return bdm.matrix[inner_mat_row, inner_mat_col]
    else
        return zero(T)
    end
end

现在，这模拟了完整的块对角矩阵，同时仅使用了一个数组的存储。事实上，即使是漂亮的印刷品也是免费的。

julia> bdm = BlockDiagonalMatrix([1 2; 3 4; 5 6], 4)
12×8 BlockDiagonalMatrix{Int64}:
 1  2  0  0  0  0  0  0
 3  4  0  0  0  0  0  0
 5  6  0  0  0  0  0  0
 0  0  1  2  0  0  0  0
 0  0  3  4  0  0  0  0
 0  0  5  6  0  0  0  0
 0  0  0  0  1  2  0  0
 0  0  0  0  3  4  0  0
 0  0  0  0  5  6  0  0
 0  0  0  0  0  0  1  2
 0  0  0  0  0  0  3  4
 0  0  0  0  0  0  5  6

julia> bdm = BlockDiagonalMatrix([1.0 2.0], 4)
4×8 BlockDiagonalMatrix{Float64}:
 1.0  2.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  2.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0  2.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  1.0  2.0

您可以了解有关如何实现自定义数组类型 here 的更多信息 - 实现内置 Array 提供的全部功能只需很少的工作。