Broadcast operations with LowerTriangularMatrix

[milankl]

Broadcasting with LowerTriangularMatrix still returns unexpected results. The biggest problem I found was dotted operations like .*= , see how suddenly a zero appears where it shouldn’t

julia> L = rand(LowerTriangularMatrix,3,3)
3×3 LowerTriangularMatrix{Float64}:
 0.998746  0.0       0.0
 0.185934  0.783976  0.0
 0.181445  0.950264  0.970376

julia> L .*= 2
3×3 LowerTriangularMatrix{Float64}:
 1.99749   0.0  0.0
 0.371868  0.0  0.0
 0.0       0.0  1.94075

This is (I think) because it loops over each index (i,j) with an @inbounds such that the @boundscheck in setindex! is disabled and then converts i,j to a single index k which isn’t safe if j>i… The problem is apparently that the broadcasting of * creates a Matrix which is then fused with .= and writes the entries back into a LowerTriangularMatrix

julia> L = rand(LowerTriangularMatrix,3,3)
3×3 LowerTriangularMatrix{Float64}:
 0.0948286  0.0        0.0
 0.10467    0.821121   0.0
 0.971535   0.0221127  0.744497

julia> L .* 2
3×3 Matrix{Float64}:
 0.189657  0.0        0.0
 0.20934   1.64224    0.0
 1.94307   0.0442254  1.48899

[milankl]

@giordano just raised another related issue

julia> A = LowerTriangularMatrix(rand(3,3))
3×3 LowerTriangularMatrix{Float64}:
 0.382579  0.0       0.0
 0.554634  0.779331  0.0
 0.453919  0.236713  0.559563

julia> sum(A)
2.966738421294002

julia> prod(A)
0.009942566052321158

The sum is correct, but the product should obviously be zero.

[maximilian-gelbrecht]

@giordano just raised another related issue

julia> A = LowerTriangularMatrix(rand(3,3))
3×3 LowerTriangularMatrix{Float64}:
 0.382579  0.0       0.0
 0.554634  0.779331  0.0
 0.453919  0.236713  0.559563

julia> sum(A)
2.966738421294002

julia> prod(A)
0.009942566052321158

The sum is correct, but the product should obviously be zero.

Is it a bug or a feature?

It's currently implemented solely with spherical harmonics coefficients in mind and for those prod(A) then returns the product of the coefficients, I actually like that. Or in more general terms: it does return the product of the lower triangular of the matrix. This could also be intended.

If the product of the full matrix is meant, yes of course that's always zero, but then one could also just add

Base.prod(A::LowerTriangularMatrix{T}) where T = zero(T)

[maximilian-gelbrecht]

To original problem: I never customised broadcast before. But something like broadcast!(*,A,....) = broadcast!(*,A.data,..) could maybe work. So "forwarding" all critical broadcast operations to the actual vector with the data. I had no time to test this (yet)

[maximilian-gelbrecht]

There are also BroadcastStyles and more specifically, ArrayStyle and DefaultArrayStyle{N}. Because essentially, we have a 2D array that should broadcast like a 1D Array.

Here's the manual on it

PS: It's still fascinating to me how customizable Julia is

I made a quick attempt and just added

Base.BroadcastStyle(::Type{<:LowerTriangularMatrix}) = Broadcast.DefaultArrayStyle{1}

but that's not working, we'll have to look a bit deeper into that.

Your suspicion about the created Matrix should be correct though, as the manual states, this is indeed the default fallback. But the Base.similar(bc::Broadcasted{DestStyle}, ::Type{ElType}) is also customizable.

[giordano]

It's currently implemented solely with spherical harmonics coefficients in mind and for those prod(A) then returns the product of the coefficients, I actually like that. Or in more general terms: it does return the product of the lower triangular of the matrix. This could also be intended.

But that's a different operation than the product of all the elements of an array. Changing the semantic of a function can cause surprising and confusing results.

[maximilian-gelbrecht]

Rename LowerTriangularMatrix to SphericalHarmonicsCoefficients and the current product would be correct ;)

But whatever, I have no strong feelings about this. Add

Base.prod(A::LowerTriangularMatrix{T}) where T = zero(T)

and you have the correct answer for the full 2D matrix that's consistent with using them for general purposes.

I'll have a bit further look in the broadcasting issue in the next days.

[maximilian-gelbrecht]

So, a hacky way to fix it, is just to make LowerTriangularMatrix a subtype of AbstractVector instead of AbstractMatrix. Then, the MWE in the OP works. The 2D indexing also still works, obviously it is not a subtype of AbstractMatrix in this case, so functions that assume that would need to be changed.

For a way how to fix this and keep it a subtype of AbstractMatrix, one would need to specialise

Base.BroadcastStyle(::Type{SrcType}) = SrcStyle()
and
Base.similar(bc::Broadcasted{DestStyle}, ::Type{ElType})

The manual gives one concrete example for this, but especially for the Base.similar I don't quite get it how to generalise this.

[milankl]

Thanks for looking into this!

So, a hacky way to fix it, is just to make LowerTriangularMatrix a subtype of AbstractVector instead of AbstractMatrix.

I tried that earlier, and thought it was a good idea as 2D indexing then comes as a bonus rather than assumed to exist. However, I ran into even more problems down the line, which gave me the feeling that <:AbstractMatrix is better. I don't have hard feelings on this and believe both is possible; both will come with advantages and disadvantages, but it's hard to foresee the consequences.

I'd suggest to try first changing the broadcasting and if that imposes even bigger obstacles then we can try <:AbstractVector again?!

[maximilian-gelbrecht]

I give up. No helpful replies on Slack and I tried a lot of different variations of what's given in the manual. Doesn't seem to work.

So I'd change it to AbstractVector

@milankl The printing / show will be no problem to adjust. There are one or two routines that need slight adjustment that assume LowerTriangularMatrix <: AbstractMatrix, but not many. Are there any other things that might need adjustment as soon as all tests pass?

Overall there a few routines that allow spectral input to be an AbstractMatrix. Personally, I would just require that all spectral input needs to be LowerTriangularMatrix. What do you think?

[milankl]

I've added a method like

function gridded(   alms::AbstractMatrix,       # spectral coefficients
                    S::SpectralTransform{NF},   # struct for spectral transform parameters
                    ) where NF                  # number format NF
 
    map = zeros(S.Grid{NF},S.nresolution)               # preallocate output
    almsᴸ = LowerTriangularMatrix{Complex{NF}}(alms)    # drop the upper triangle and convert to NF
    gridded!(map,almsᴸ,S)                               # now execute the in-place version
    return map
end

So that someone can also come with a Matrix of coefficients without the need to do LowerTriangularMatrix(::Matrix) first. So the ::AbstractMatrix here isn't an issue and we'd only need to do vec(alms) inside the function body.

[milankl]

Thanks for trying!! Then let's do <:AbstractVector and proceed. I feel like we've already spent too much time on this ;)

[milankl]

The only request I have it that we retain show similar to this

julia> rand(LowerTriangularMatrix,4,4)
4×4 LowerTriangularMatrix{Float64}:
 0.0025516  0.0        0.0        0.0
 0.2128     0.906056   0.0        0.0
 0.589329   0.0822436  0.0122374  0.0
 0.865687   0.838233   0.741264   0.685711

As otherwise troubleshooting with spectral fields isn't easy...

[maximilian-gelbrecht]

In case we ever come back to this, to remind myself and archive it, this was my best effort so far:

# Broadcasting 

import Base.BroadcastStyle, Base.similar

struct LowerTriangularStyle <: Broadcast.AbstractArrayStyle{1} end

Base.BroadcastStyle(::Type{<:LowerTriangularMatrix}) = LowerTriangularStyle()

Base.similar(bc::Broadcast.Broadcasted{LowerTriangularStyle}, ::Type{ElType}) where ElType = similar(LowerTriangularMatrix{Eltype}, axes(bc))

LowerTriangularStyle(::Val{0}) = LowerTriangularStyle()
LowerTriangularStyle(::Val{1}) = LowerTriangularStyle()
LowerTriangularStyle(::Val{2}) = LowerTriangularStyle()

No errors, just the same wrong results