如何保留带状对角矩阵并在 julia 的大矩阵中用 0 替换其他元素

Question

我想保留对角矩阵，并在 julia 的大矩阵中用 0 替换其他元素。例如，A 是我的矩阵，我只想保留 A 中的 2 x 2 对角线元素，并将所有其他元素替换为 0。B 矩阵是我想要的。我只是想知道有没有一种优雅的方法来做到这一点。

A = [1 2 3 4 5 6 7 8; 
     1 2 3 4 5 6 7 8; 
     1 2 3 4 5 6 7 8; 
     1 2 3 4 5 6 7 8; 
     1 2 3 4 5 6 7 8;
     1 2 3 4 5 6 7 8; 
     1 2 3 4 5 6 7 8; 
     1 2 3 4 5 6 7 8]

B = [1 2 0 0 0 0 0 0; 
     1 2 0 0 0 0 0 0; 
     0 0 3 4 0 0 0 0; 
     0 0 3 4 0 0 0 0; 
     0 0 0 0 5 6 0 0; 
     0 0 0 0 5 6 0 0; 
     0 0 0 0 0 0 7 8; 
     0 0 0 0 0 0 7 8]

Answer 1

可能某处有一个高级别 API，但是，编写一个 for 循环应该可行。

function change_zero!(a)
 lo = 1
 for j in 1:size(a, 2)
   if isodd(j)
     lo += 2
   end
   for i in 1:lo-3
     a[i,j]=0
   end
   for i in lo:size(a,1)
      a[i,j]=0
   end
 end
 a
end

change_zero!(A)

Answer 2

为了完整起见，以下是对您（原始）标题中问题的惯用回答：

julia> function zeronondiag!(A)
           di = diagind(A)
           for i in eachindex(A)
               i ∉ di && (A[i] = zero(A[i]))
           end
           return A
       end
zeronondiag! (generic function with 1 method)

julia> zeronondiag!(copy(A))
8×8 Matrix{Int64}:
 1  0  0  0  0  0  0  0
 0  2  0  0  0  0  0  0
 0  0  3  0  0  0  0  0
 0  0  0  4  0  0  0  0
 0  0  0  0  5  0  0  0
 0  0  0  0  0  6  0  0
 0  0  0  0  0  0  7  0
 0  0  0  0  0  0  0  8

请注意 diagind returns 线性索引的范围，因此 ∉ 检查相当有效。

julia> diagind(A)
1:9:64

您应该能够使用与 BlockArrays.jl 非常相似的逻辑来获得块对角线形式。

Answer 3

实现此目的的最短代码是 using BlockBandedMatrices，如下所示：

julia> BlockBandedMatrix(A,repeat([2],4),repeat([2],4),(0,0))
4×4-blocked 8×8 BlockBandedMatrix{Int64}:
 1  2  │  ⋅  ⋅   │  ⋅  ⋅  │  ⋅  ⋅
 1  2  │  ⋅  ⋅   │  ⋅  ⋅  │  ⋅  ⋅
 ──────┼────────┼───────┼──────
  ⋅  ⋅  │  3  4  │  ⋅  ⋅  │  ⋅  ⋅
  ⋅  ⋅  │  3  4  │  ⋅  ⋅  │  ⋅  ⋅
 ──────┼────────┼───────┼──────
  ⋅  ⋅  │  ⋅  ⋅   │  5  6 │  ⋅  ⋅
  ⋅  ⋅  │  ⋅  ⋅   │  5  6 │  ⋅  ⋅
 ──────┼────────┼───────┼──────
  ⋅  ⋅  │  ⋅  ⋅   │  ⋅  ⋅  │  7  8
  ⋅  ⋅  │  ⋅  ⋅   │  ⋅  ⋅  │  7  8

另一件值得关注的事情是 BandedMatrices 包，它提供了此类功能以及一组专用的线性代数函数，用于高效处理此类数据结构。

julia> using BandedMatrices

julia> BandedMatrix(A, (1,0))
8×8 BandedMatrix{Int64} with bandwidths (1, 0):
 1  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 1  2  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  2  3  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  3  4  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  4  5  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  5  6  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  6  7  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  7  8

Answer 4

方法一：
这是使用 CartesianIndices:

的一种有趣方式

julia> B = zero(A);
julia> blocksize = 2; 
julia> d = diag(CartesianIndices(A))
8-element Vector{CartesianIndex{2}}:
 CartesianIndex(1, 1)
 CartesianIndex(2, 2)
 CartesianIndex(3, 3)
 CartesianIndex(4, 4)
 CartesianIndex(5, 5)
 CartesianIndex(6, 6)
 CartesianIndex(7, 7)
 CartesianIndex(8, 8)

julia> for p in Iterators.partition(d, blocksize)
         block = first(p):last(p)
         B[block] .= @view A[block]
       end

在每次迭代中，Iterators.partition returns blocksize 个对角线元素，因此属于一个块的所有对角线元素。

关于 CartesianIndices 的一个有用的事情是范围已经按块运行：CartesianIndex(1,1):CartesianIndex(2,2) returns CartesianIndex (1,1),(2,1) 的值， (1,2) 和 (2,2) 自动。所以 first(p):last(p) returns 在每次迭代中我们想要的块中所有元素的索引。

方法二：
在这种情况下，因为事物是对称的，非 CartesianIndices 方式也非常简洁：

julia> B = zero(A);
julia> for b in Iterators.partition(1:size(A, 1), blocksize)
         B[b,b] .= @view A[b,b]
       end
julia> B
8×8 Matrix{Int64}:
 1  2  0  0  0  0  0  0
 1  2  0  0  0  0  0  0
 0  0  3  4  0  0  0  0
 0  0  3  4  0  0  0  0
 0  0  0  0  5  6  0  0
 0  0  0  0  5  6  0  0
 0  0  0  0  0  0  7  8
 0  0  0  0  0  0  7  8

在第一次迭代中（作为例子），partition returns 1:2（假设blocksize = 2），所以我们分配给B[1:2, 1:2]这是我们想要的街区。

将其概括为允许 non-standard 索引（例如 OffsetArrays）：

julia> for (r, c) in zip(Iterators.partition.(axes(A), blocksize)...)
         B[r, c] .= @view A[r, c] 
       end

（感谢@phipsgabler 的 .= @view 建议，避免不必要的分配，以及 axes(A) 方法。）

Answer 5

就像其中一个答案一样，我更喜欢为这种操作编写一个内部带有循环的简单函数。一个稍微更通用的函数，允许您指定 off-diagonal 元素的值和块对角线的大小：

function setoffdiag!(A::AbstractMatrix{T}, value::T = zero(T); block_size::Integer = 1) where {T}

    m, n = size(A)
    k = 1
    r = 1

    @inbounds for j = 1:n

        @inbounds for i = 1:(k - 1)
            A[i, j] = value
        end

        @inbounds for i = (k + block_size):m
            A[i, j] = value
        end

        k += (r == block_size) * block_size
        r += 1 - (r == block_size) * block_size
    end

    return A
end

Answer 6

这个快速修复可能会成功：（假设输入是一个平方矩阵）

function two_by_two_diag(A)
    B = zeros(Int64,size(A))
    for i in 1:2:size(A,1)
        B[i,i] = A[i,i]
        if i != size(A,1)
            B[i,i+1] = A[i,i+1]
            B[i+1,i] = A[i+1,i]
            B[i+1, i+1] = A[i+1, i+1]
        end
    end
    return B
end

如何保留带状对角矩阵并在 julia 的大矩阵中用 0 替换其他元素

How to keep a banded diagonal matrix and replace other elements by 0 in a large matrix for julia

julia