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

Question

我想保留对角矩阵并将其他元素替换为 0 在朱莉娅的大矩阵中。例如，A我有一个矩阵，我只想保留 2 x 2 对角线元素A，并将所有其他元素替换为 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

方法 1：
\n这里有一个有趣的方法，使用CartesianIndices：

\n

julia> B = zero(A);\njulia> blocksize = 2; \njulia> d = diag(CartesianIndices(A))\n8-element Vector{CartesianIndex{2}}:\n CartesianIndex(1, 1)\n CartesianIndex(2, 2)\n CartesianIndex(3, 3)\n CartesianIndex(4, 4)\n CartesianIndex(5, 5)\n CartesianIndex(6, 6)\n CartesianIndex(7, 7)\n CartesianIndex(8, 8)\n\njulia> for p in Iterators.partition(d, blocksize)\n         block = first(p):last(p)\n         B[block] .= @view A[block]\n       end\n\n

\n

在每次迭代中，Iterators.partition返回blocksize对角元素的数量，即属于一个块的所有对角元素。

\n

一个有用的事情CartesianIndices是范围已经按块操作：自动CartesianIndex(1,1):CartesianIndex(2,2)返回CartesianIndex(1,1)、(2,1)、(1,2) 和 (2,2) 的值。因此，first(p):last(p)在每次迭代中返回我们想要的块中所有元素的索引。

\n\n

方法 2：
\n在这种情况下，由于事物是对称的，因此非CartesianIndices方法也非常整洁和简单：

\n

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

\n

在第一次迭代中（作为示例），partition返回1:2（假设blocksize = 2），因此我们分配给B[1:2, 1:2]我们想要的块。

\n

概括地说，允许非标准索引（例如 OffsetArrays）：

\n

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

\n\n

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

\n