5

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.

nbro
  • 15,395
  • 32
  • 113
  • 196
vincet
  • 917
  • 3
  • 13
  • 26
  • 1
    Just FYI julia does not use NA values. We have NaN – Alexander Morley Aug 19 '16 at 14:03
  • x-ref: https://groups.google.com/forum/#!topic/julia-users/VGbdlfbRfbc – StefanKarpinski Aug 19 '16 at 15:01
  • 1
    Also, consider still using zeros instead of NA (or NaN). Why? 1. The meaning of linear algebra is mostly lost using NAs and this would prevent any use of existing linear algebra code. 2. NAs are known to be infectious and bug generating. If you describe the use-case some more, this choice could be clearer. – Dan Getz Aug 19 '16 at 16:51

3 Answers3

6

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

Alexander Morley
  • 4,111
  • 12
  • 26
  • now , l want to explore in a foor loop just the lower triangular matrix. How can l do that efficiently form a fast optimization point of view because l have a matrix of 5000 by 5000 . thanks – vincet Aug 22 '16 at 14:25
4

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.

Graham
  • 7,431
  • 18
  • 59
  • 84
Michael Ohlrogge
  • 10,559
  • 5
  • 48
  • 76
1

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 matrix:

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 matrix:

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
Zranz
  • 57
  • 9