# Triangular Arrays in julia?

**URL:** <https://discourse.julialang.org/t/triangular-arrays-in-julia/930>\
**Category:** Offtopic\
**Created:** [December 13, 2016, 8:11pm UTC](https://discourse.julialang.org/t/triangular-arrays-in-julia/930 "2016-12-13T20:11:00Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ward9250](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ward9250/32/42768_2.png) [@Ward9250](https://discourse.julialang.org/u/Ward9250)\
**Post date:** [December 13, 2016, 8:11pm UTC](https://discourse.julialang.org/t/triangular-arrays-in-julia/930/1 "2016-12-13T20:11:00Z")

</div>

Does anyone know if triangular arrays have been implemented in julia or in a package yet?

> **[Triangular array](https://en.wikipedia.org/wiki/Triangular_array)**
>
> In mathematics and computing, a triangular array of numbers, polynomials, or the like, is a doubly indexed sequence in which each row is only as long as the row's own index. That is, the ith row contains only i elements.
> Notable particular examples include these:
> Triangular arrays of integers in which each row is symmetric and begins and ends with 1 are sometimes called generalized Pascal triangles; examples include Pascal's triangle, the Narayana numbers, and the triangle of Eulerian numbers.
> ...

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 13, 2016, 8:18pm UTC](https://discourse.julialang.org/t/triangular-arrays-in-julia/930/2 "2016-12-13T20:18:22Z")

</div>

Not that I know of. It would be really nice to have a package which makes it easy to do parallel computations on triangular arrays, since that’s very common in numerical analysis.

---

<div class="post-metadata">

**Author:** ![Ward9250](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ward9250/32/42768_2.png) [@Ward9250](https://discourse.julialang.org/u/Ward9250)\
**Post date:** [December 13, 2016, 11:33pm UTC](https://discourse.julialang.org/t/triangular-arrays-in-julia/930/3 "2016-12-13T23:33:24Z")

</div>

I was going to write a small type wrapping a vector, since what I really want is the indexing behaviour and memory savings, and element-wise operations, as opposed to any more mathematical operations. I could try start a dedicated package, but I’m not a mathematician or expert enough to know what people would want.

---

<div class="post-metadata">

**Author:** ![Keno](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/keno/32/285_2.png) [@Keno](https://discourse.julialang.org/u/Keno)\
**Post date:** [December 13, 2016, 11:36pm UTC](https://discourse.julialang.org/t/triangular-arrays-in-julia/930/4 "2016-12-13T23:36:47Z")

</div>

Base has `UpperTriangular` and `LowerTriangular` with dispatch to various BLAS functions. I do believe they still store the entire matrix through.

---

<div class="post-metadata">

**Author:** ![Ward9250](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ward9250/32/42768_2.png) [@Ward9250](https://discourse.julialang.org/u/Ward9250)\
**Post date:** [December 13, 2016, 11:39pm UTC](https://discourse.julialang.org/t/triangular-arrays-in-julia/930/5 "2016-12-13T23:39:28Z")

</div>

Yes, I think that’s true, the upper triangular and lower triangular refers to which values are not zero.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 13, 2016, 11:43pm UTC](https://discourse.julialang.org/t/triangular-arrays-in-julia/930/6 "2016-12-13T23:43:45Z")

</div>

They do store the whole matrix. Also, they are restricted to being square.

---

<div class="post-metadata">

**Author:** ![andreasnoack](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasnoack/32/27_2.png) [@andreasnoack](https://discourse.julialang.org/u/andreasnoack)\
**Post date:** [December 14, 2016, 1:07am UTC](https://discourse.julialang.org/t/triangular-arrays-in-julia/930/7 "2016-12-14T01:07:16Z")

</div>

The wiki article is quite brief. Could you explain a little more about what kind of operations you’d like a triangular array to support. In particular, what is “the indexing behavior” and “element-wise” operations in the context of triangular arrays?

---

<div class="post-metadata">

**Author:** ![Ward9250](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ward9250/32/42768_2.png) [@Ward9250](https://discourse.julialang.org/u/Ward9250)\
**Post date:** [December 14, 2016, 11:37am UTC](https://discourse.julialang.org/t/triangular-arrays-in-julia/930/8 "2016-12-14T11:37:50Z")

</div>

So if you’re familiar with the `dist` class from R, it is very similar.

It represents a matrix, or rather the lower half of a matrix underneath the diagonal (as the matrix is assumed to be symmetrical with `0` values as the diagonal, there is no point in storing both halfs), and so each row has a different number of cells.

Row one is one cell, row two is two cells, three has three and so on:

```julia
as.dist(matrix(rexp(100, rate=.1), ncol=10))
             1 2 3 4 5 6
2 30.96809401                                                            
3 4.29517367 25.35477096                                                
4 11.05042388 47.12602515 9.05635823                                    
5 24.80794737 15.96819158 5.57638651 10.38488874                        
6 3.73206725 26.75391206 8.32492987 8.59412784 0.35015252            
7 33.00415806 0.20160700 2.03088283 0.45323016 13.87089156 22.28731310
8 9.10296273 7.98441481 1.90013604 8.93347893 7.18713469 7.02904816
9 3.74577306 3.19449860 5.11027106 1.63343495 14.50446960 1.38334442
10 24.93290636 2.64181538 1.29009316 0.04882485 3.75045984 19.01684707

```

Looking at the structure in R, you can see the data is laid out linearly:

```julia
> str(as.dist(matrix(rexp(100, rate=.1), ncol=10)))
Class 'dist' atomic [1:45] 0.441 4.219 7.834 0.657 0.112 ...
  ..- attr(*, "Size")= int 10
  ..- attr(*, "call")= language as.dist.default(m = matrix(rexp(100, rate = 0.1), ncol = 10))
  ..- attr(*, "Diag")= logi FALSE
  ..- attr(*, "Upper")= logi FALSE

```

And the array support indexing linearly i.e. `m[43]` but also by row/column `m[4, 5]`.

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [December 14, 2016, 2:59pm UTC](https://discourse.julialang.org/t/triangular-arrays-in-julia/930/9 "2016-12-14T14:59:21Z")

</div>

Question: is it possible to store upper/lower triangular matrices in compact form and still pass that representation efficiently to BLAS operations? I.e. can that storage format be represented with the kinds of strides that BLAS understands?

---

<div class="post-metadata">

**Author:** ![nalimilan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nalimilan/32/147_2.png) [@nalimilan](https://discourse.julialang.org/u/nalimilan)\
**Post date:** [December 14, 2016, 3:01pm UTC](https://discourse.julialang.org/t/triangular-arrays-in-julia/930/10 "2016-12-14T15:01:08Z")

</div>

Sounds like `Symmetric`, except that it saves some space by not storing upper entries (as our implementation currently does) nor diagonal entries (which are assumed to be zero). I’d start by making `Symmetric` more efficient as regards memory layout in Julia Base. Not storing diagonal elements sounds like a secondary issue since for large arrays the cost will be relatively low.

LAPACK supports a [packed storage format](http://www.netlib.org/lapack/lug/node123.html) for symmetric matrices, so you could still get high performance for common linear algebra operations (e.g. [matrix \* vector product](http://www.netlib.org/lapack/explore-html/d6/d30/group __single__ blas__level2_gad1af0d0777da05d1c27ea99a69c8017c.html#gad1af0d0777da05d1c27ea99a69c8017c)). (These optimized operations can be added progressively.)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 14, 2016, 5:10pm UTC](https://discourse.julialang.org/t/triangular-arrays-in-julia/930/11 "2016-12-14T17:10:35Z")

</div>

> [@StefanKarpinski](#):
>
> Question: is it possible to store upper/lower triangular matrices in compact form and still pass that representation efficiently to BLAS operations? I.e. can that storage format be represented with the kinds of strides that BLAS understands?

Is it still going to be much better to use BLAS? For square lower triangular matrices you have less than half of the operations to do, and with `@polly` being added, could the difference between BLAS and not-BLAS be around 2x?

(Is there a good way to switch between the BLAS and generic operations to do a quick test what the difference is like?)

---

<div class="post-metadata">

**Author:** ![Ward9250](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ward9250/32/42768_2.png) [@Ward9250](https://discourse.julialang.org/u/Ward9250)\
**Post date:** [December 14, 2016, 5:30pm UTC](https://discourse.julialang.org/t/triangular-arrays-in-julia/930/12 "2016-12-14T17:30:14Z")

</div>

I did not know about Symmetric, I can live with storing the 0 values in the diagonal!

I’ve no idea when it comes to BLAS I’m afraid.

---

<div class="post-metadata">

**Author:** ![nalimilan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nalimilan/32/147_2.png) [@nalimilan](https://discourse.julialang.org/u/nalimilan)\
**Post date:** [December 14, 2016, 5:36pm UTC](https://discourse.julialang.org/t/triangular-arrays-in-julia/930/13 "2016-12-14T17:36:08Z")

</div>

> [@ChrisRackauckas](#):
>
> Is it still going to be much better to use BLAS? For square lower triangular matrices you have less than half of the operations to do, and with @polly being added, could the difference between BLAS and not-BLAS be around 2x?

Even if there are only about half operations to do, that 1/2 is a constant factor: the cost is still increasing with the size of the matrix. So that doesn’t really make a difference from standard arrays, except for changing the threshold above which calling BLAS is beneficial.

> [@ChrisRackauckas](#):
>
> (Is there a good way to switch between the BLAS and generic operations to do a quick test what the difference is like?)

For matrix multiplication, you can call `Base.LinAlg.generic_matmatmul` directly. I’m not sure for what operations and types BLAS is called, and for which pure Julia code is used instead.

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [December 14, 2016, 6:25pm UTC](https://discourse.julialang.org/t/triangular-arrays-in-julia/930/16 "2016-12-14T18:25:30Z")

</div>

You could just wrap a float in a new type and define the needed operations. This should dispatch to the generic versions of the linear algebra library.

---

<div class="post-metadata">

**Author:** ![andreasnoack](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasnoack/32/27_2.png) [@andreasnoack](https://discourse.julialang.org/u/andreasnoack)\
**Post date:** [December 14, 2016, 8:20pm UTC](https://discourse.julialang.org/t/triangular-arrays-in-julia/930/17 "2016-12-14T20:20:13Z")

</div>

For linear algebra, there is a significant overhead when using the packed format. This is partly because the indexing operations are slightly more complicated but more importantly, it is because the compact storage doesn’t allow blocked (BLAS-3) linear algebra operations. This was the main reason I decided not to use the compact storage scheme when I implemented `XTriangular` and `Symmetric` since I thought that 2x memory saving was not sufficient to justify slower linear algebra. However, for other applications, this might, of course, be completely different.

An interesting alternative is the [Rectangular Full Packaged storage scheme](https://github.com/andreasnoack/LinearAlgebra.jl/blob/master/src/rectfullpacked.jl) which maintains BLAS-3 performance without wasting half the memory. However, I never really needed the memory saving so the development stalled.

---

<div class="post-metadata">

**Author:** ![andreasnoack](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasnoack/32/27_2.png) [@andreasnoack](https://discourse.julialang.org/u/andreasnoack)\
**Post date:** [December 14, 2016, 8:24pm UTC](https://discourse.julialang.org/t/triangular-arrays-in-julia/930/18 "2016-12-14T20:24:02Z")

</div>

Thanks for the example. Is the linear index similar to the indexing into the underlying vector or is it similar to a full matrix representation, i.e. can you `m[n^2]` or is the last index `m[div(n*(n-1),2)]`?

Can you tell a little more about what kind of operations you’d like to perform on the triangular/symmetric array?

---

<div class="post-metadata">

**Author:** ![mkborregaard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkborregaard/32/556_2.png) [@mkborregaard](https://discourse.julialang.org/u/mkborregaard)\
**Post date:** [December 16, 2016, 10:19am UTC](https://discourse.julialang.org/t/triangular-arrays-in-julia/930/19 "2016-12-16T10:19:32Z")

</div>

It’d be cool to be able to do stuff like

```julia
a = UpperTriangular(ones(4,4))
a[:] = 1:10

```

perhaps one solution would be to make it possible to make UpperTriangular views on the square matrix that just stored the indices of the upper triangle and allowed vector operations? That should not interfere with BLAS.

---

<div class="post-metadata">

**Author:** ![andreasnoack](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasnoack/32/27_2.png) [@andreasnoack](https://discourse.julialang.org/u/andreasnoack)\
**Post date:** [December 16, 2016, 11:27am UTC](https://discourse.julialang.org/t/triangular-arrays-in-julia/930/20 "2016-12-16T11:27:24Z")

</div>

There have been discussions about adding functionality to only reference the non-zero part of a matrix more generally (i.e. also for `SparseMatrixCSC` and `SymTiridagonal`) and that could possibly solve this issue.

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [December 19, 2016, 11:26am UTC](https://discourse.julialang.org/t/triangular-arrays-in-julia/930/21 "2016-12-19T11:26:21Z")

</div>

That would clash with generic linear indexing. Of course, we’re looking to move away from that since it’s not an efficient mechanism for writing generic algorithms on all array types, but giving it a _different_ meaning seems like asking for trouble.

---

<div class="post-metadata">

**Author:** ![mkborregaard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkborregaard/32/556_2.png) [@mkborregaard](https://discourse.julialang.org/u/mkborregaard)\
**Post date:** [December 19, 2016, 3:07pm UTC](https://discourse.julialang.org/t/triangular-arrays-in-julia/930/22 "2016-12-19T15:07:30Z")

</div>

Ok, I see. Well, looking forward to see what the new alternatives will be 🙂 It’d be cool if they also made it easy to work with only the defined elements of Triangular matrices.

[Next page](https://discourse.julialang.org/t/triangular-arrays-in-julia/930.md?page=2)