vincet
vincet

Reputation: 977

Lower triangular matrix in julia

I have the number of columns equals the number of rows. And the the diagonal is is equal to zero. How can I build this matrix?

#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

I tried this

x=rand(4,4)
4x4 Array{Float64,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 LowerTriangular{Float64,Array{Float64,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

but l'm looking for something like this

c=LowerTriangular(x)
4x4 LowerTriangular{Float64,Array{Float64,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

The diagonal should be equal to zero.

Upvotes: 5

Views: 7649

Answers (3)

Zranz
Zranz

Reputation: 67

The accepted solution does not index all elements of the vector in order, and the output matrix has repeated elements. The formula is wrong. Here is my proposal, inspired by previous answers:

For lower triangular entries:

function vec2ltri{T}(v::Vector{T})
    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

which will output:

julia> vec2ltri(collect(1:6))
4×4 Array{Int64,2}:
 0  0  0  0
 1  0  0  0
 2  4  0  0
 3  5  6  0

For upper triangular entries:

function vec2utri{T}(v::Vector{T})
    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

which will output:

julia> vec2utri(collect(1:6))
4×4 Array{Int64,2}:
 0  1  2  4
 0  0  3  5
 0  0  0  6
 0  0  0  0

Upvotes: 1

Michael Ohlrogge
Michael Ohlrogge

Reputation: 10980

Here is something taking inspiration from Code by Stefan Karpinski on the Julia User's list:

function vec2ltri_alt{T}(v::AbstractVector{T}, 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 Array{Any,2}:
 0  NaN  NaN  NaN
 1    0  NaN  NaN
 2    3    0  NaN
 3    4    6    0

Note: If desired, check out the official documentation on the ternary operator for a bit more clarity on what is going on with the ? ... : syntax here.

For those looking for a more "standard" diagonal matrix solution:

Here is a version that creates a more standard solution:

function vec2ltri{T}(v::AbstractVector{T}, 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 Array{Int64,2}:
 1  0  0
 2  3  0
 3  4  6

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

If you want upper triangular: instead of lower, just change the i<=j to i>=j.

Other random tools Note also functions like tril!(a) which will convert in place a given matrix to lower triangular, replacing everything above the main diagonal with zeros. See the Julia documentation for more info on this function, as well as various other related tools.

Upvotes: 4

Alexander Morley
Alexander Morley

Reputation: 4181

You might want to use a list comprehension. But it would be good if you could give more information in the question about what you are trying to do.

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

which will give:

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

For any number of rows and columns. And then you can work from there.

See docs for list comprehensions & conditionals:

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

Upvotes: 6

Related Questions