LowerTriangularMatrices

LowerTriangularMatrices is a submodule that has been developed for SpeedyWeather.jl which is technically independent (SpeedyWeather.jl however imports it and so does SpeedyTransforms) and can also be used without running simulations. It is just not put into its own respective repository.

This module defines LowerTriangularArray, a lower triangular matrix format, which in contrast to LinearAlgebra.LowerTriangular does not store the entries above the diagonal. SpeedyWeather.jl uses LowerTriangularArray which is defined as a subtype of AbstractArray to store the spherical harmonic coefficients (see Spectral packing). For 2D LowerTriangularArray the alias LowerTriangularMatrix exists. Higher dimensional LowerTriangularArray are 'batches' of 2D LowerTriangularMatrix. So, for example a $(10\times 10\times 10)$ LowerTriangularArray holds 10 LowerTriangularMatrix of size $(10\times 10)$ in one array.

LowerTriangularMatrix is actually a vector

LowerTriangularMatrix and LowerTriangularArray can in many ways be used very much like a Matrix or Array, however, because they unravel the lower triangle into a vector their dimensionality is one less than their Array counterparts. A LowerTriangularMatrix should therefore be treated as a vector rather than a matrix with some (limited) added functionality to allow for matrix-indexing (vector or flat-indexing is the default though). More details below.

Creation of `LowerTriangularArray`

A LowerTriangularMatrix and LowerTriangularArray can be created using zeros, ones, rand, or randn

julia> using SpeedyWeather.LowerTriangularMatrices
julia> L = rand(LowerTriangularMatrix{Float32}, 5, 5)15-element, 5x5 LowerTriangularMatrix{Float32}
 0.0701276  0.0       0.0        0.0       0.0
 0.653854   0.469095  0.0        0.0       0.0
 0.199065   0.926987  0.110597   0.0       0.0
 0.996502   0.307519  0.107916   0.344857  0.0
 0.252573   0.692873  0.0690462  0.453723  0.0721205
julia> L2 = rand(LowerTriangularArray{Float32}, 5, 5, 5)15×5 LowerTriangularArray{Float32, 2, Matrix{Float32}}:
 0.828761  0.728882   0.450571   0.173966    0.213849
 0.514123  0.923521   0.593651   0.917904    0.330551
 0.608286  0.224022   0.785178   0.0742081   0.892341
 0.497111  0.460817   0.255194   0.147877    0.285532
 0.138869  0.995486   0.0180768  0.793455    0.36218
 0.383271  0.0522733  0.567525   0.214809    0.261065
 0.453477  0.557918   0.493369   0.378259    0.649281
 0.602534  0.246538   0.831705   0.507474    0.788921
 0.259997  0.755477   0.872156   0.608314    0.290026
 0.389943  0.023003   0.992606   0.871165    0.570309
 0.614961  0.832771   0.947476   0.00257832  0.799297
 0.808585  0.470364   0.339224   0.268653    0.764273
 0.457523  0.980233   0.324052   0.192175    0.593033
 0.429167  0.722978   0.377841   0.0877343   0.442746
 0.340129  0.959374   0.333665   0.032434    0.542169

or the undef initializer LowerTriangularMatrix{Float32}(undef, 3, 3). The element type is arbitrary though, you can use any type T too.

Note how for a matrix both the upper triangle and the lower triangle are shown in the terminal. The zeros are evident. However, for higher dimensional LowerTriangularArray we fall back to show the unravelled first two dimensions. Hence, here, the first column is the first matrix with 15 elements forming a 5x5 matrix, but the zeros are not shown.

Alternatively, it can be created through conversion from Array, which drops the upper triangle entries and sets them to zero (which are not stored however)

julia> M = rand(Float16, 3, 3)3×3 Matrix{Float16}:
 0.919   0.9116  0.1616
 0.3257  0.127   0.9497
 0.925   0.7837  0.5894
julia> L = LowerTriangularMatrix(M)6-element, 3x3 LowerTriangularMatrix{Float16}
 0.919   0.0     0.0
 0.3257  0.127   0.0
 0.925   0.7837  0.5894
julia> M2 = rand(Float16, 3, 3, 2)3×3×2 Array{Float16, 3}:
[:, :, 1] =
 0.787    0.581  0.78
 0.737    0.957  0.4272
 0.08936  0.824  0.9585

[:, :, 2] =
 0.01611   0.478   0.2627
 0.858     0.5645  0.1284
 0.007812  0.708   0.0293
julia> L2 = LowerTriangularArray(M2)6×2 LowerTriangularArray{Float16, 2, Matrix{Float16}}:
 0.787    0.01611
 0.737    0.858
 0.08936  0.007812
 0.957    0.5645
 0.824    0.708
 0.9585   0.0293

Size of `LowerTriangularArray`

There are three different ways to describe the size of a LowerTriangularArray. For example with L

julia> L = rand(LowerTriangularMatrix, 5, 5)15-element, 5x5 LowerTriangularMatrix{Float64}
 0.674048  0.0       0.0       0.0       0.0
 0.140671  0.178235  0.0       0.0       0.0
 0.498739  0.272608  0.164466  0.0       0.0
 0.876611  0.300971  0.765141  0.439821  0.0
 0.678391  0.102206  0.951725  0.838675  0.540317

we have (additional dimensions follow naturally thereafter)

1-based vector indexing (default)

julia> size(L)    # equivalently size(L, OneBased, as=Vector)(15,)

The lower triangle is unravelled hence the number of elements in the lower triangle is returned.

1-based matrix indexing

julia> size(L, as=Matrix)    # equivalently size(L, OneBased, as=Matrix)(5, 5)

If you think of a LowerTriangularMatrix as a matrix this is the most intuitive size of L, which, however, does not agree with the size of the underlying data array (hence it is not the default).

0-based matrix indexing

Because LowerTriangularArrays are used to represent the coefficients of spherical harmonics which are commonly indexed based on zero (i.e. starting with the zero mode representing the mean), we also add ZeroBased to get the corresponding size.

julia> size(L, ZeroBased, as=Matrix)(4, 4)

which is convenient if you want to know the maximum degree and order of the spherical harmonics in L. 0-based vector indexing is not implemented.

Indexing `LowerTriangularArray`

We illustrate the two types of indexing LowerTriangularArray supports.

Matrix indexing, by denoting two indices, column and row [l, m, ..]
Vector/flat indexing, by denoting a single index [lm, ..].

The matrix index works as expected

julia> L15-element, 5x5 LowerTriangularMatrix{Float64}
 0.674048  0.0       0.0       0.0       0.0
 0.140671  0.178235  0.0       0.0       0.0
 0.498739  0.272608  0.164466  0.0       0.0
 0.876611  0.300971  0.765141  0.439821  0.0
 0.678391  0.102206  0.951725  0.838675  0.540317
julia> L[2, 2]0.17823520845110719

But the single index skips the zero entries in the upper triangle, i.e. a 2, 2 index points to the same element as the index 6

julia> L[6]0.17823520845110719

which, important, is different from single indices of an AbstractMatrix

julia> Matrix(L)[6]0.0

which would point to the first element in the upper triangle (hence zero).

In performance-critical code a single index should be used, as this directly maps to the index of the underlying data vector. The matrix index is somewhat slower as it first has to be converted to the corresponding single index.

Consequently, many loops in SpeedyWeather.jl are build with the following structure

julia> n, m = size(L, as=Matrix)(5, 5)
julia> ij = 00
julia> for j in 1:m, i in j:n
           ij += 1
           L[ij] = i+j
       end

which loops over all lower triangle entries of L::LowerTriangularArray and the single index ij is simply counted up. However, one could also use [i, j] as indices in the loop body or to perform any calculation (i+j here).

`end` doesn't work for matrix indexing

Indexing LowerTriangularMatrix and LowerTriangularArray in matrix style ([i, j]) with end doesn't work. It either returns an error or wrong results as the end is lowered by Julia to the size of the underlying flat array dimension.

The setindex! functionality of matrixes will throw a BoundsError when trying to write into the upper triangle of a LowerTriangularArray, for example

julia> L[2, 1] = 0    # valid index0
julia> L[1, 2] = 0    # invalid index in the upper triangleERROR: BoundsError: attempt to access 15-element, 5x5 LowerTriangularMatrix{Float64} at index [1, 2]

But reading from it will just return a zero

julia> L[2, 3]     # in the upper triangle0.0

Higher dimensional LowerTriangularArray can be indexed with multidimensional array indices like most other arrays types. Both the vector index and the matrix index for the lower triangle work as shown here

julia> L = rand(LowerTriangularArray{Float32}, 3, 3, 5)6×5 LowerTriangularArray{Float32, 2, Matrix{Float32}}:
 0.715975  0.46843    0.931653  0.770093  0.32089
 0.254572  0.602575   0.625341  0.224724  0.56322
 0.539403  0.240248   0.959195  0.781257  0.830921
 0.602112  0.0259385  0.258088  0.988082  0.379072
 0.694874  0.35007    0.528572  0.761276  0.187882
 0.269102  0.356771   0.776419  0.994024  0.851801
julia> L[2, 1] # second lower triangle element of the first lower triangle matrix0.25457168f0
julia> L[2, 1, 1] # (2,1) element of the first lower triangle matrix0.25457168f0

The setindex! functionality follows accordingly.

Iterators

An iterator over all entries in the array can be created with eachindex

julia> L = rand(LowerTriangularArray, 5, 5, 5)15×5 LowerTriangularArray{Float64, 2, Matrix{Float64}}:
 0.659927   0.648042   0.646975   0.29414    0.715417
 0.243402   0.282938   0.275351   0.79192    0.106051
 0.649028   0.535636   0.700707   0.0840294  0.2754
 0.168419   0.359482   0.444256   0.597392   0.0054515
 0.372432   0.0719973  0.220474   0.185322   0.397055
 0.834506   0.500914   0.325992   0.768174   0.00536504
 0.901698   0.256945   0.560146   0.484899   0.499889
 0.0578256  0.106157   0.589795   0.323754   0.918656
 0.214793   0.669033   0.596973   0.507699   0.198438
 0.771488   0.889643   0.0718354  0.728853   0.824028
 0.343845   0.15741    0.645165   0.59735    0.519777
 0.694673   0.458893   0.96655    0.497754   0.761636
 0.940299   0.0260604  0.272538   0.301714   0.0207686
 0.609249   0.507024   0.209843   0.73768    0.473342
 0.963064   0.721696   0.693275   0.92303    0.363696
julia> for ij in eachindex(L)
           # do something
       end
julia> eachindex(L)Base.OneTo(75)

In order to only loop over the harmonics (essentially the horizontal, ignoring other dimensions) use eachharmonic

julia> eachharmonic(L)Base.OneTo(15)

If you only want to loop over the other dimensions use eachmatrix

julia> eachmatrix(L)CartesianIndices((5,))

together they can be used as

julia> for k in eachmatrix(L)
           for lm in eachharmonic(L)
               L[lm, k]
           end
       end

Note that k is a CartesianIndex that will loop over all other dimensions, whether there's only 1 (representing a 3D variable) or 5 (representing a 6D variable with the first two dimensions being a lower triangular matrix).

Linear algebra with `LowerTriangularArray`

The LowerTriangularMatrices module's main purpose is not linear algebra, and typical matrix operations will not work with LowerTriangularMatrix because it's treated as a vector not as a matrix, meaning that the following will not work as expected

julia> L = rand(LowerTriangularMatrix{Float32}, 3, 3)6-element, 3x3 LowerTriangularMatrix{Float32}
 0.506564  0.0        0.0
 0.526874  0.0118356  0.0
 0.854014  0.25879    0.310888
julia> L * LERROR: MethodError: no method matching *(::LowerTriangularMatrix{Float32}, ::LowerTriangularMatrix{Float32})
The function `*` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  *(::Any, ::Any, ::Any, ::Any...)
   @ Base operators.jl:596
  *(::ChainRulesCore.ZeroTangent, ::Any)
   @ ChainRulesCore ~/.julia/packages/ChainRulesCore/6Pucz/src/tangent_arithmetic.jl:104
  *(::Any, ::ChainRulesCore.NotImplemented)
   @ ChainRulesCore ~/.julia/packages/ChainRulesCore/6Pucz/src/tangent_arithmetic.jl:38
  ...
julia> inv(L)ERROR: MethodError: no method matching inv(::LowerTriangularMatrix{Float32})
The function `inv` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  inv(::CoordinateTransformations.SphericalFromCartesian)
   @ CoordinateTransformations ~/.julia/packages/CoordinateTransformations/jYKYg/src/coordinatesystems.jl:249
  inv(::Proj.Transformation; always_xy, area, ctx)
   @ Proj ~/.julia/packages/Proj/8LMij/src/coord.jl:196
  inv(::Geodesy.ECEFfromLLA)
   @ Geodesy ~/.julia/packages/Geodesy/otwWk/src/transformations.jl:191
  ...

And many other operations that require L to be a AbstractMatrix which it isn't. In contrast, typical vector operations like a scalar product between two "LowerTriangularMatrix" vectors does work

julia> L' * L1.4273062f0

Broadcasting with `LowerTriangularArray`

In contrast to linear algebra, many element-wise operations work as expected thanks to broadcasting, so operations that can be written in . notation whether implicit +, 2*, ... or explicitly written .+, .^, ... or via the @. macro

julia> L + L6-element, 3x3 LowerTriangularMatrix{Float32}
 1.01313  0.0        0.0
 1.05375  0.0236712  0.0
 1.70803  0.517579   0.621776
julia> 2L6-element, 3x3 LowerTriangularMatrix{Float32}
 1.01313  0.0        0.0
 1.05375  0.0236712  0.0
 1.70803  0.517579   0.621776
julia> @. L + 2L - 1.1*L / L^26-element, 3x3 LowerTriangularMatrix{Float64}
 -0.651799    0.0       0.0
 -0.507164  -92.9046    0.0
  1.27401    -3.47419  -2.60558

GPU

LowerTriangularArray{T, N, ArrayType} wraps around an array of type ArrayType. If this array is a GPU array (e.g. CuArray), all operations are performed on GPU as well (work in progress). The implementation was written so that scalar indexing is avoided in almost all cases, so that GPU operation should be performant. To use LowerTriangularArray on GPU you can e.g. just adapt an existing LowerTriangularArray.

using Adapt
L = rand(LowerTriangularArray{Float32}, 5, 5, 5)
L_gpu = adapt(CuArray, L)

Function and type index

SpeedyWeather.LowerTriangularMatrices.LowerTriangularArray — Type

A lower triangular array implementation that only stores the non-zero entries explicitly. L<:AbstractArray{T,N-1} although we do allow both "flat" N-1-dimensional indexing and additional N-dimensional or "matrix-style" indexing.

Supports n-dimensional lower triangular arrays, so that for all trailing dimensions L[:, :, ..] is a matrix in lower triangular form, e.g. a (5x5x3)-LowerTriangularArray would hold 3 lower triangular matrices.