Transforming vector to upper triangular matrix

[Theodore Papamarkou]

Is there a neat and efficient way to solve the following problem in Julia? Assume a vector that contains n*(n+1)/2 elements. The aim is to pack these elements into the diagonal and upper right triangular portion of an n*n matrix. For example, the vector 1:10 would be transformed to the 4*4 matrix 

 1  2  3   4

 0  5  6   7

 0  0  8   9

 0  0  0  10

I am not so familiar with handling of symmetric and triangular matrices in Julia. Is there some built-in support for this?

[Stefan Karpinski]

I'm not at a computer right now, but this should be possible with a comprehension. Something like this:

[ i<=j ? v[f(i,j)] : 0 for i=1:n, j=1:n ]

Where f is an appropriate function of i and j for indexing into v.

[Theodore Papamarkou]

Thanks a lot Stefan, I haven't thought of combining a comprehension with the ternary operator, this is very helpful. I'm not in front of the computer either at the moment; after your tip I should be able to complete the rest of details (such as f) when I get back to the laptop. Thanks!

[Gunnar Farnebäck]

Straighforward and efficient but not so neat:

function f1{T}(x::Vector{T})

    n = int(sqrt(2 * length(x) + 0.25) - 0.5)

    y = zeros(T, n, n)

    k = 1

    for i = 1:n

        for j = i:n

            y[i,j] = x[k]

            k += 1

        end

    end

    y

end

More compact but not efficient:

function f2{T}(x::Vector{T})

    n = int(sqrt(2 * length(x) + 0.25) - 0.5)

    y = zeros(T, n, n)

    y[tril(trues(n, n))] = x

    y'

end

[Theodore Papamarkou]

Cheers Gunnar, that's time saving and helpful! I'll check it later on today from my computer.

[Stefan Karpinski]

I just worked out the correct expression:

[ i<=j ? v[j*(j-1)/2+i] : 0 for i=1:n, j=1:n ]

This is a one-liner and is very efficient, so it's what I would do. If you wanted to write this more generically, it would be something like this:

function vec2utri{T}(v::AbstractVector{T}, z::T=zero(T))

    n = length(v)

    s = iround((sqrt(8n+1)-1)/2)

    s*(s+1)/2 == n || error("vec2utri: length of vector is not triangular")

    [ i<=j ? v[j*(j-1)/2+i] : z for i=1:s, j=1:s ]

end

That's sufficiently complicated that maybe we want to provide it as a built-in, but hopefully with a better name. And of course, we'd want the opposite operation too: extracting the upper-triangular part of a matrix as a vector.

[Theodore Papamarkou]

Thank you Stefan, this saved me time and gave me a generic way of doing it. With this function I can now perform operations on vectors (rather than whole triangular matrices) and convert the final vector to a matrix which is very efficient. Yes, the converse (triangular matrix to vector) would also be helpful in various cases. May I also mention one more possibility as a built in, which is trivial, after having the function you wrote; vector to symmetric matrix, which means filling in the elements below the diagonal too as x_ij=x_ji. Thanks, your code helps get used to the Julia idiom.

[Diego Javier Zea]

What do you think about something like that ?

importall Base

type UpperTriDiag{T<:Number} <: AbstractMatrix{T}
  data::Vector{T}
  size::(Int,Int)
  UpperTriDiag(data::Vector{T}) = new( data, _get_uppertridiag_size( data ) )
end

size{T<:Number}(mat::UpperTriDiag{T}) = mat.size
size{T<:Number}(mat::UpperTriDiag{T},i::Int) = i<3 ? mat.size[i] : thorw("Only 2x2 matrix")

function _get_uppertridiag_size(vec::Vector)
  len = length(vec)
  numerator = sqrt( 8len + 1 ) -1
  if numerator % 2 != 0
    throw("Not valid vector length")
  end
  n = div( numerator , 2 )
  (n,n)
end

function getindex{T<:Number}(mat::UpperTriDiag{T},i::Int,j::Int)
  n=size(mat)[1]
  i > j ? zero(T) : mat.data[ i + (n*(j-1)) - ((n-(j-1))*(j-1)) - (0.5 * ( (j-1)*(j-2) ) ) ]
end

ndims{T<:Number}(mat::UpperTriDiag{T}) = 2

Example:

julia> UpperTriDiag{Int}([1,2,3])
2x2 Int64 UpperTriDiag:
 1  2
 0  3

julia> UpperTriDiag{Int}([1,2,3,4])
ERROR: "Not valid vector length"
 in _get_uppertridiag_size at none:5
 in UpperTriDiag at none:4

julia> UpperTriDiag{Int}([1,2,3,4,5,6])
3x3 Int64 UpperTriDiag:
 1  2  4
 0  3  5
 0  0  6

julia> a = UpperTriDiag{Int}([1,2,3,4,5,6])
3x3 Int64 UpperTriDiag:
 1  2  4
 0  3  5
 0  0  6

julia> a[1,1]
1

julia> a[1,2]
2

julia> a[3,1]
0

julia> dump(a)
UpperTriDiag{Int64}
  data: Array(Int64,(6,)) [1, 2, 3, 4, 5, 6]
  size: (Int64,Int64) (3,3)

[Theodore Papamarkou]

Thanks Diego. It is helpful to get so many responses with code and extra functionality. I will save your functions too in case i need them, for now Stefan's one liner comprehension does the trick :)

[Stefan Karpinski]

Just realized that I should have used div or >>> in this instead of /:

[ i<=j ? v[j*(j-1)>>>1+i] : 0 for i=1:n, j=1:n ]

[Gunnar Farnebäck]

Comprehension wins over first creating with zeros. But this is somewhat more efficient:

k=0; [ i<=j ? (k+=1; v[k]) : 0 for i=1:n, j=1:n ]

[Stefan Karpinski]

Yeah, that makes sense since it avoids a bunch of arithmetic in favor of an extra counter in the inner loop.

[Gunnar Farnebäck]

Actually with >>> in the arithmetic there's nearly no difference at all on my computer. Let's call it a draw. :-)

[Bradley Alpert]

Very nice: judicious use of side-effects in comprehensions.  Quibble, to get the correct answer:

k=0; [ j<=i ? (k+=1; v[k]) : 0 for i=1:n, j=1:n ]'

[Viral Shah]

On Thursday, April 25, 2013 8:16:43 PM UTC+5:30, Theodore Papamarkou wrote:

I am not so familiar with handling of symmetric and triangular matrices in Julia. Is there some built-in support for this?

See `base/linalg/triangular.jl` and `base/linalg/hermitian.jl` for this.

-viral

 

[Theodore Papamarkou]

Cheers Viral, thanks for the reference. From your experience, is it preferable to perform arithmetic on vectors (which is possible in the context of the problem I am coding) and simply convert the final vector to a hermitian matrix (using Stefan's-Gunnar's-Bradley's clever synthesis above), or would you recommend doing arithmetic on triangular matrices directly instead? My instinct would be to go for the former, I am asking cause you know better the (triangular and hermitian) matrix algebra performance of Julia.

[Viral Shah]

The way you use the symmetric matrices, is to specifically pass Hermitian(A) to linear algebra functions that can perform better on symmetric matrices.

My personal preference would be to work with matrices, unless you are working with vectors the way you are due to the rest of your program. Although the comprehensions are clever, the complexity is still O(n^2), and you iterate over n^2 elements. Working with matrices, I would typically write nested for loops that only touch the non-zero part (0.5*n^2), and that should end up being faster. This is what we do in functions like triu() etc. It would be good to verify these claims.

You should certainly try out the interfaces for the special matrix types, and we can certainly enhance them to make them more usable.

-viral

[Theodore Papamarkou]

Thanks Viral, I didn't know all this about the special matrix types in Julia. I am willing to try out the interfaces. The only difficulty I have at the moment is that I haven't managed to formulate the problem I am working on (computation of Hessian matrices via automatic differentiation) in a matrix fashion. Applying the chain rule twice to get the second order derivatives (i.e. the elements of the Hessian) is naturally done element-wise... the only other expression I can think of at the moment for the Hessian is to see it as the Jacobian of the involved gradients, which is not of much help though. If I manage to think of a matrix formulation of the problem, I will certainly follow your advice to try out the special matrix interfaces.

[Viral Shah]

Why not iterate over the 2d matrix with linear indexing?

[Theodore Papamarkou]

Yes, linear indexing is one more possibility... I just found the attached recent report which suggests a way of computing the Hessian by performing only o(n) calculations, which is impressive... I think I will finish first the naive implementation to have a working case, and then will try to speed it up incrementally by applying the more efficient approaches, as my head started spinning from reading :)

[Theodore Papamarkou]

P.S. The attachment :)

[Tom Short]

Interesting. Also, did you see that the Philip Caplan (the author) is a student of Steven G. Johnson, a significant Julia contributor? Maybe they have some code or examples they could share.

[Theodore Papamarkou]

I didn't notice this Tom, and it is potentially handy. Have you ever contacted Steven or Philip before? If so, would you like to initiate contact to see if there is some code available to share...?

[Theodore Papamarkou]

Thank you very much Philip!

On Friday, April 26, 2013 7:48:04 PM UTC+1, Philip Caplan wrote:

Hello, I've posted my code at: web.mit.edu/pcaplan/www/codes/SecondDerivOpt.zip, enjoy!

Philip

[Ahmed Mazari]

How to transform a lower triangular matrix to a vector ? l used the function vec( )  and it doesn't work l got this code error

ERROR: ArgumentError: Triangular matrix must have two dimensions
 in similar at linalg/triangular.jl:27
 in reshape at abstractarray.jl:213
 in vec at abstractarraymath.jl:14