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

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

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

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

其中 0 是与 A 大小相同且填充零的矩阵。有什么方便的方法吗？

换句话说，是否有等效于 Python 的 linalg.block_diag（来自 scipy）？

解决方法

您可以使用 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

,

如果您不想为矩阵和 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 提供的全部功能。

,

您可以使用 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