矩阵取下三角部分

Question

我想知道在Julia中是否有一种命令或包可以直接提取矩阵的下三角部分，不包括对角线。我可以调用R命令来实现这一点（如gdata包的lowerTriangle），但我想知道Julia是否有类似的功能。例如，假设我有以下矩阵：

1.0    0.751   0.734    
0.751   1.0    0.948    
0.734  0.948    1.0

我不想创建一个类似于下三角矩阵的东西，就像这样：

NA     NA      NA     
0.751   NA      NA    
0.734  0.948    NA

但是将矩阵的下部分提取为数组：0.751 0.734 0.948

- coolsv

---

可能是Julia中的下三角矩阵的重复问题。 - HarmonicaMuse

不，我的问题不同。我会修改我的问题。 - coolsv

对于那些投票将此问题关闭为重复的人：这个问题最多只是另一个问题的相反，但实际上它们是具有不同要求的完全不同的问题。不要让标题中的一些共同单词欺骗你。 - Sundar R

抱歉…尝试寻找一些暗示性关键字和合适的标题。 - coolsv

Answer 1

如果您可以接受创建一个下三角矩阵作为中间步骤，那么您可以使用逻辑索引和 tril! 以及额外的参数来获得所需结果。

julia> M = [1.0 0.751 0.734
0.751 1.0 0.948
0.734 0.948 1.0];
julia> v = M[tril!(trues(size(M)), -1)]
3-element Array{Float64, 1}:
0.751
0.734
0.948

trues函数返回一个布尔值为true的、形状与M相同的数组。由于我们仅需要矩阵的一部分，tril!函数会将其进一步削减。这个函数的第二个参数指定了从哪个超对角线开始削减，我们用它来避免主对角线上的值。

我们使用上述方式得到一个索引数组来操作M，然后返回所需数值的数组。

- Sundar R

---

这只有在您不介意放弃出现在下三角部分的零时才有效。您可以使用逻辑数组来实现相同的想法，以使其更安全、更高效：M[tril!(trues(size(M)), -1)]。 - mbauman

那是个好点子。我会马上用你的解决方案更新答案，谢谢。 - Sundar R

尽管这种中间分配很好，但最好避免它。对称距离矩阵也经常出现这个问题，有很多人想要一个无分配的解决方案。 - juliohm

这个解决方案看起来不错，谢谢大家！它并不需要复杂的代码，很简单且快速。 - coolsv

Answer 2

使用推导式：

julia> [M[m, n] for m in 2:size(M, 1) for n in 1:m-1]
3-element Array{Float64,1}:
 0.751
 0.734
 0.948

但是它比 sundar/Matt B. 的解决方案 慢得多：

lower_triangular_1(M) = [M[m, n] for m in 2:size(M, 1) for n in 1:m-1]
lower_triangular_2(M) = [M[m, n] for n in 1:size(M, 2) for m in n+1:size(M, 1)]
lower_triangular_3(M) = M[tril!(trues(size(M)), -1)]
using BenchmarkTools
using LinearAlgebra  # avoid warning in 0.7

M=rand(100, 100)

使用Julia版本0.7.0-alpha.0进行测试：

julia> @btime lower_triangular_1(M);
  73.179 μs (10115 allocations: 444.34 KiB)

julia> @btime lower_triangular_2(M);
  71.157 μs (10117 allocations: 444.41 KiB)

julia> @btime lower_triangular_3(M);
  16.325 μs (6 allocations: 40.19 KiB)

虽然不够优雅，但使用@views可以更快：

function lower_triangular_4(M)
    # works only for square matrices
    res = similar(M, ((size(M, 1)-1) * size(M, 2)) ÷ 2)
    start_idx = 1
    for n = 1:size(M, 2)-1
        @views column = M[n+1:end, n]
        last_idx = start_idx -1 + length(column)
        @views res[start_idx:last_idx] = column[:]
        start_idx = last_idx + 1
    end
end

julia> @btime lower_triangular_4(M);
  4.272 μs (101 allocations: 44.95 KiB)

- ederag

---

一个不错的方法！顺便说一下，我正在尝试避免使用循环，但还是谢谢你！ - coolsv