# Half Vectorization

**URL:** <https://discourse.julialang.org/t/half-vectorization/7399>\
**Category:** New to Julia\
**Tags:** question\
**Created:** [November 30, 2017, 3:12am UTC](https://discourse.julialang.org/t/half-vectorization/7399 "2017-11-30T03:12:21Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![slakeon](https://avatars.discourse-cdn.com/v4/letter/s/e274bd/32.png) [@slakeon](https://discourse.julialang.org/u/slakeon)\
**Post date:** [November 30, 2017, 3:12am UTC](https://discourse.julialang.org/t/half-vectorization/7399/1 "2017-11-30T03:12:21Z")

</div>

Is there an algorithm that can do half vectorization, the vech ([here](https://en.wikipedia.org/wiki/Vectorization_(mathematics)#Half-vectorization)) or generates elimination of duplication matrices ([here](https://en.wikipedia.org/wiki/Duplication_and_elimination_matrices))? I’m trying to recover only the unique elements of a covariance matrix in vector form.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [November 30, 2017, 4:22am UTC](https://discourse.julialang.org/t/half-vectorization/7399/2 "2017-11-30T04:22:06Z")

</div>

Sure, the algorithm is: write a loop. (Or rather, two nested loops are probably the easiest way to implement vech.)

---

<div class="post-metadata">

**Author:** ![GunnarFarneback](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunnarfarneback/32/1827_2.png) [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)\
**Post date:** [November 30, 2017, 10:13am UTC](https://discourse.julialang.org/t/half-vectorization/7399/3 "2017-11-30T10:13:47Z")

</div>

If compact code is more important than performance you can do `A[tril(trues(A))]`. Otherwise the loop solution can easily be implemented by a comprehension:  
`[A[i, j] for i = 1:size(A, 1), j = 1:size(A, 2) if i >= j]`

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [November 30, 2017, 1:14pm UTC](https://discourse.julialang.org/t/half-vectorization/7399/4 "2017-11-30T13:14:04Z")

</div>

> [@GunnarFarneback](#):
>
> If compact code is more important than performance you can do A[tril(trues(A))]. Otherwise the loop solution can easily be implemented by a comprehension: `[A[i, j] for i = 1:size(A, 1), j = 1:size(A, 2) if i >= j]`

Actually, `A[tril(trues(A))]` is significantly _faster_ than the comprehension on my machine.

Explicit loops are even better, about 50x faster than this comprehension on my machine:

```julia
vech0(A) = A[tril(trues(A))]
vech1(A) = [A[i, j] for i = 1:size(A, 1), j = 1:size(A, 2) if i >= j]

function vech(A::AbstractMatrix{T}) where T
    m = LinAlg.checksquare(A)
    v = Vector{T}((m*(m+1))>>1)
    k = 0
    for j = 1:m, i = j:m
        @inbounds v[k += 1] = A[i,j]
    end
    return v
end

B = rand(1000,1000); B = B + B';
using BenchmarkTools, Compat
@btime vech0($B); @btime vech1($B); @btime vech($B);

```

gives

```julia
  893.148 μs (16 allocations: 4.06 MiB)
  21.599 ms (24 allocations: 5.00 MiB)
  385.942 μs (2 allocations: 3.82 MiB)

```

on my machine with Julia 0.6.

I knew loops would be faster than comprehensions, but I’m honestly surprised that the difference is so large. The comprehension version does multiple allocations (comprehensions with a filter, like `i >= j`, have to grow the array as they go along), loops over the entire array, and is not type-stable according to `@code_warntype`. Still, it seems like filtered comprehensions could faster…

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [November 30, 2017, 1:20pm UTC](https://discourse.julialang.org/t/half-vectorization/7399/5 "2017-11-30T13:20:46Z")

</div>

> [@stevengj](#):
>
> Still, it seems like filtered comprehensions could faster…

In Julia 0.7 on my machine, things are somewhat better:

```julia
  vech0: 1.200 ms (9 allocations: 4.06 MiB)
  vech1: 5.514 ms (23 allocations: 5.00 MiB)
  vech: 377.758 μs (2 allocations: 3.82 MiB)

```

Even so, the factor of 14 difference is larger than I would have initially thought.

---

<div class="post-metadata">

**Author:** ![GunnarFarneback](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunnarfarneback/32/1827_2.png) [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)\
**Post date:** [November 30, 2017, 1:22pm UTC](https://discourse.julialang.org/t/half-vectorization/7399/6 "2017-11-30T13:22:47Z")

</div>

For completeness,

```julia
i = ceil(Int, sqrt(2 * n + 0.25) - 0.5)
j = n - i * (i - 1) ÷ 2

```

gives an enumeration of the lower triangle indices.

---

<div class="post-metadata">

**Author:** ![GunnarFarneback](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunnarfarneback/32/1827_2.png) [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)\
**Post date:** [November 30, 2017, 1:56pm UTC](https://discourse.julialang.org/t/half-vectorization/7399/7 "2017-11-30T13:56:59Z")

</div>

Right, that’s what you get for guessing at performance with insufficient insight. It should have been obvious that it would require pretty advanced compiler optimizations to correctly predict the size of the filtered comprehension and avoid growing and other overhead.

---

<div class="post-metadata">

**Author:** ![slakeon](https://avatars.discourse-cdn.com/v4/letter/s/e274bd/32.png) [@slakeon](https://discourse.julialang.org/u/slakeon)\
**Post date:** [November 30, 2017, 6:46pm UTC](https://discourse.julialang.org/t/half-vectorization/7399/8 "2017-11-30T18:46:42Z")

</div>

Thank you all! These work great!

---

<div class="post-metadata">

**Author:** ![reinhardhansen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/reinhardhansen/32/29383_2.png) [@reinhardhansen](https://discourse.julialang.org/u/reinhardhansen)\
**Post date:** [September 18, 2020, 6:55pm UTC](https://discourse.julialang.org/t/half-vectorization/7399/9 "2020-09-18T18:55:33Z")

</div>

This code will not run with Julia 1.5.  
`vech0` will run with the modification: `trues(A)` → `trues(size(A))`  
`vech` will not run. (I did make the change `LinAlg.checksquare` → `LinearAlgebra.checksquare)`

```julia
ERROR: LoadError: MethodError: no method matching Array{Float64,1}(::Int64)
Closest candidates are:
  Array{Float64,1}() where T at boot.jl:425
  Array{Float64,1}(::UndefInitializer, ::Int64) where T at boot.jl:406
  Array{Float64,1}(::UndefInitializer, ::Int64...) where {T, N} at boot.jl:412

```

What else must be modified?

```julia
vech0(A) = A[tril(trues(size(A)))]
vech1(A) = [A[i, j] for i = 1:size(A, 1), j = 1:size(A, 2) if i >= j]

function vech(A::AbstractMatrix{T}) where T
    m = LinearAlgebra.checksquare(A)
    v = Vector{T}((m*(m+1))>>1)
    k = 0
    for j = 1:m, i = j:m
        @inbounds v[k += 1] = A[i,j]
    end
    return v
end

B = rand(1000,1000); B = B + B';
using BenchmarkTools, Compat
@btime vech0($B); @btime vech1($B); @btime vech($B);

```

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [September 18, 2020, 7:01pm UTC](https://discourse.julialang.org/t/half-vectorization/7399/10 "2020-09-18T19:01:35Z")

</div>

> [@reinhardhansen](#):
>
> `Vector{T}((m*(m+1))>>1)`

it’s erroring about this line, I think it should be

```julia
Vector{T}(undef, (m*(m+1))>>1)

```

? to initialize

---

<div class="post-metadata">

**Author:** ![reinhardhansen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/reinhardhansen/32/29383_2.png) [@reinhardhansen](https://discourse.julialang.org/u/reinhardhansen)\
**Post date:** [September 18, 2020, 7:07pm UTC](https://discourse.julialang.org/t/half-vectorization/7399/11 "2020-09-18T19:07:50Z")

</div>

> [@jling](#):
>
> `Vector{T}(undef, (m*(m+1))>>1)`

Thank you! That solved it. Below are the times in Julia 1.5 on my laptop, in case someone is interested.

```julia
537.001 μs (8 allocations: 4.06 MiB)
  4.250 ms (20 allocations: 5.00 MiB)
  304.750 μs (2 allocations: 3.82 MiB)

```