Julia中的下三角矩阵

Posted 2023-02-16

【中文标题】Julia中的下三角矩阵

【英文标题】：Lower triangular matrix in julia

【发布时间】：2016-08-19 12:52:47

【问题描述】：

我的列数等于行数。并且对角线等于零。如何构建这个矩阵？

#mat
#     [,1] [,2] [,3] [,4]
#[1,]   0   NA   NA   NA
#[2,]    1   0   NA   NA
#[3,]    2    4   0   NA
#[4,]    3    5    6   0

我试过了

x=rand(4,4)
4x4 ArrayFloat64,2:
 0.60064   0.917443  0.561744   0.135717 
 0.106728  0.72391   0.0894174  0.0656103
 0.410262  0.953857  0.844697   0.0375045
 0.476771  0.778106  0.469514   0.398846 

c=LowerTriangular(x)

4x4 LowerTriangularFloat64,ArrayFloat64,2:
 0.60064   0.0       0.0       0.0     
 0.106728  0.72391   0.0       0.0     
 0.410262  0.953857  0.844697  0.0     
 0.476771  0.778106  0.469514  0.398846

但我正在寻找这样的东西

c=LowerTriangular(x)
4x4 LowerTriangularFloat64,ArrayFloat64,2:
 0.0   NA      NA       NA    
 0.106728  0.0    NA      NA    
 0.410262  0.953857  0.0 NA     
 0.476771  0.778106  0.469514  0

对角线应该为零。

【问题讨论】：

仅供参考，julia 不使用 NA 值。我们有 NaN

外部参考：groups.google.com/forum/#!topic/julia-users/VGbdlfbRfbc

另外，请考虑使用零而不是 NA（或 NaN）。为什么？ 1. 线性代数的含义在使用 NA 时大部分会丢失，这将阻止对现有线性代数代码的任何使用。 2. 已知 NA 具有传染性并会产生错误。如果你再描述一下用例，这个选择可能会更清楚。

【参考方案1】：

您可能想要使用列表推导。但是，如果您可以在问题中提供有关您正在尝试做什么的更多信息，那就太好了。

numrows =4
numcols = 4
[ x>y ? 1 : (x == y ? 0 : NaN) for x in 1:numrows, y in 1:numcols]

这将给出：

 0  NaN  NaN  NaN
 1    0  NaN  NaN
 1    1    0  NaN
 1    1    1    0

对于任意数量的行和列。然后你就可以从那里开始工作了。

请参阅文档了解列表推导和条件：

http://docs.julialang.org/en/release-0.4/manual/arrays/#comprehensions

http://docs.julialang.org/en/release-0.4/manual/control-flow/#man-conditional-evaluation

【讨论】：

现在，我想在 foor 循环中探索下三角矩阵。我如何有效地做到这一点，形成一个快速优化的观点，因为我有一个 5000 x 5000 的矩阵。谢谢

【参考方案2】：

这里有一些灵感来自 Stefan Karpinski 在 Julia 用户list 上的代码：

function vec2ltri_altT(v::AbstractVectorT, z::T=zero(T))
    n = length(v)
    v1 = vcat(0,v)
    s = round(Int,(sqrt(8n+1)-1)/2)
    s*(s+1)/2 == n || error("vec2utri: length of vector is not triangular")
    s+=1
    [ i>j ? v1[round(Int, j*(j-1)/2+i)] : (i == j ? z : NaN) for i=1:s, j=1:s ]
end

julia> vec2ltri_alt(collect(1:6))
4x4 ArrayAny,2:
 0  NaN  NaN  NaN
 1    0  NaN  NaN
 2    3    0  NaN
 3    4    6    0

注意：如果需要，请查看三元运算符上的 official documentation，以更清楚地了解此处的 ? ... : 语法发生了什么。

对于那些寻找更“标准”的对角矩阵解决方案的人：

这是一个创建更标准解决方案的版本：

function vec2ltriT(v::AbstractVectorT, z::T=zero(T))
    n = length(v)
    s = round(Int,(sqrt(8n+1)-1)/2)
    s*(s+1)/2 == n || error("vec2utri: length of vector is not triangular")
    [ i>=j ? v[round(Int, j*(j-1)/2+i)] : z for i=1:s, j=1:s ]
end

a = vec2ltri(collect(1:6))

julia> a = vec2ltri(collect(1:6))
3x3 ArrayInt64,2:
 1  0  0
 2  3  0
 3  4  6

julia> istril(a)  ## verify matrix is lower triangular
true

如果你想要上三角：而不是下三角，只需将i<=j 更改为i>=j。

其他随机工具注意还有像tril!(a)这样的函数，它将给定矩阵转换为下三角形，用零替换主对角线以上的所有内容。有关此功能以及各种其他相关工具的更多信息，请参阅 Julia documentation。

【参考方案3】：

接受的解决方案没有按顺序索引向量的所有元素，并且输出矩阵有重复的元素。公式是错误的。 这是我的建议，灵感来自以前的答案：

对于下三角矩阵：

function vec2ltriT(v::VectorT)
    d = length(v)
    n = Int((sqrt(8d+1)+1)/2)
    n*(n-1)/2 == d || error("vec2ltri: length of vector is not triangular")
    [ i>j ? v[Int((2n-j)*(j-1)/2)+i-j] : 0 for i=1:n, j=1:n ]
end

将输出：

julia> vec2ltri(collect(1:6))
4×4 ArrayInt64,2:
 0  0  0  0
 1  0  0  0
 2  4  0  0
 3  5  6  0

对于上三角矩阵：

function vec2utriT(v::VectorT)
    d = length(v)
    n = Int((sqrt(8d+1)+1)/2)
    n*(n-1)/2 == d || error("vec2utri: length of vector is not triangular")
    [ i<j ? v[Int((j-1)*(j-2)/2)+i] : 0 for i=1:n, j=1:n ]
end

将输出：

julia> vec2utri(collect(1:6))
4×4 ArrayInt64,2:
 0  1  2  4
 0  0  3  5
 0  0  0  6
 0  0  0  0