我想保留julia的对角矩阵，并在大矩阵中用0替换其他元素。例如，A是我的矩阵，我只想在A中保留2乘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]
方法1:
这里有一个使用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返回blocksize个对角元素，因此属于一个块的所有对角元素。
CartesianIndices的一个有用之处在于，范围已经按块操作：CartesianIndex(1,1):CartesianIndex(2,2)自动返回(1,1(、(2,1(、(1,2(和(2,2(的CartesianIndex值。所以first(p):last(p)在每次迭代中返回我们想要的块中所有元素的索引。
方法2:
在这种情况下，因为事物是对称的，所以非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返回1:2(假设blocksize = 2(，因此我们将其分配给B[1:2, 1:2]，这是我们想要的块。
将其推广到允许非标准索引(例如OffsetArrays(：
julia> for (r, c) in zip(Iterators.partition.(axes(A), blocksize)...)
B[r, c] .= @view A[r, c]
end
(感谢@phipsgabler提出的避免不必要分配的.= @view建议以及axes(A)方法。(
可能在某个地方有一个高级别的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)
实现这一点的最短代码是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
为了完整起见，以下是您(原始(标题中问题的惯用答案：
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返回一系列线性索引，所以∉检查是相当有效的。
julia> diagind(A)
1:9:64
您应该能够使用与BlockArrays.jsl非常相似的逻辑来获得块对角线形式。
就像其中一个答案一样，我更喜欢编写一个内部有循环的简单函数来进行这种操作。一个稍微更通用的函数，允许您指定非对角线元素的值和块对角线的大小：
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
这个快速修复可能会奏效：(假设输入是一个平方矩阵(
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
using LinearAlgebra
d = diag(CartesianIndices(A))
blocksize = 2;
mapreduce(p->A[first(p):last(p)],(x,y)->cat(x,y,dims=(1,2)), Iterators.partition(d, blocksize))
低于性能高得多的版本(>10X(
function blockmat(A,bs)
n=size(A, 1)
B = zero(A);
doff=1
soff=1
for g in 1:div(n,bs)
for i in 1:bs
copyto!(B,doff,A,soff,bs)
doff += n
soff += n
end
doff+=bs
soff+=bs
end
B
end
julia> @btime begin
for b in Iterators.partition(1:size(A, 1), blocksize)
B[b,b] .= @view A[b,b]
end
B
end
2.156 μs (37 allocations: 2.41 KiB)
15×15 Matrix{Int64}:
1 16 31 0 0 0 0 0 0 0 0 0 0 0 0
2 17 32 0 0 0 0 0 0 0 0 0 0 0 0
3 18 33 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 49 64 79 0 0 0 0 0 0 0 0 0
0 0 0 50 65 80 0 0 0 0 0 0 0 0 0
0 0 0 51 66 81 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 97 112 127 0 0 0 0 0 0
0 0 0 0 0 0 98 113 128 0 0 0 0 0 0
0 0 0 0 0 0 99 114 129 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 145 160 175 0 0 0
0 0 0 0 0 0 0 0 0 146 161 176 0 0 0
0 0 0 0 0 0 0 0 0 147 162 177 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 193 208 223
0 0 0 0 0 0 0 0 0 0 0 0 194 209 224
0 0 0 0 0 0 0 0 0 0 0 0 195 210 225
julia> @btime blockmat(A,3)
207.407 ns (1 allocation: 1.98 KiB)
15×15 Matrix{Int64}:
1 16 31 0 0 0 0 0 0 0 0 0 0 0 0
2 17 32 0 0 0 0 0 0 0 0 0 0 0 0
3 18 33 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 49 64 79 0 0 0 0 0 0 0 0 0
0 0 0 50 65 80 0 0 0 0 0 0 0 0 0
0 0 0 51 66 81 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 97 112 127 0 0 0 0 0 0
0 0 0 0 0 0 98 113 128 0 0 0 0 0 0
0 0 0 0 0 0 99 114 129 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 145 160 175 0 0 0
0 0 0 0 0 0 0 0 0 146 161 176 0 0 0
0 0 0 0 0 0 0 0 0 147 162 177 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 193 208 223
0 0 0 0 0 0 0 0 0 0 0 0 194 209 224
0 0 0 0 0 0 0 0 0 0 0 0 195 210 225
这里有任何残留物的版本
function blockmatwr(A,bs)
n=size(A, 1)
B = zero(A);
doff=1
soff=1
foreach(div(n,bs)) do
for i in 1:bs
copyto!(B,doff,A,soff,bs)
doff += n
soff += n
end
doff+=bs
soff+=bs
end
r=%(n,bs)
foreach(r) do
copyto!(B,doff,A,soff,r)
doff += n
soff += n
end
B
end
julia> A=reshape(1:225,15,15);
julia> @btime blockmatwr(A,4)
201.155 ns (1 allocation: 1.98 KiB)
15×15 Matrix{Int64}:
1 16 31 46 0 0 0 0 0 0 0 0 0 0 0
2 17 32 47 0 0 0 0 0 0 0 0 0 0 0
3 18 33 48 0 0 0 0 0 0 0 0 0 0 0
4 19 34 49 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 65 80 95 110 0 0 0 0 0 0 0
0 0 0 0 66 81 96 111 0 0 0 0 0 0 0
0 0 0 0 67 82 97 112 0 0 0 0 0 0 0
0 0 0 0 68 83 98 113 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 129 144 159 174 0 0 0
0 0 0 0 0 0 0 0 130 145 160 175 0 0 0
0 0 0 0 0 0 0 0 131 146 161 176 0 0 0
0 0 0 0 0 0 0 0 132 147 162 177 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 193 208 223
0 0 0 0 0 0 0 0 0 0 0 0 194 209 224
0 0 0 0 0 0 0 0 0 0 0 0 195 210 225