Spherical Harmonic bugs/enhancements

[MikaelSlevinsky]

This issue is a placeholder for spherical harmonic bugs/enhancements that will be part of a patch:

• the way the butterfly algorithm is currently implemented, it only works for dimensions that are powers of two, though nothing stops this from working for general dimensions. This afflicts FastSphericalHarmonicPlan and ThinSphericalHarmonicPlan. For the time being, it's recommended to just pad with zeros to the nextpow2.

• the commands sph2fourier and fourier2sph should work on different band limits in l and m.

• the command plan_sph2fourier is not type-stable as it very crudely dispatches to a different SphericalHarmonicPlan based on the bandlimit. This is negligible compared to pre-computation costs for large degrees (and is somewhat likened to Base's factorize).

• the convenience constructors sphrand, sphrandn, etc... take m and n and create an array with dimensions m and 2n-1 which is arguably misleading.

• it would be great to have an allocation-free FFTW override for converting the bivariate Fourier array to function values on the sphere at tensor product grids equispaced in angle. (Added starting with 67a6b56)

• remove segmentation fault when using an even number of columns, thanks to @meggart for finding this. (Added starting with c0e4a54)

• add methods for complex data,

• and, NUFFT variants for nonuniform point.

[MikaelSlevinsky]

As this eventually plugs into ApproxFun's multivariate side, @dlfivefifty, do you have a recommendation regarding the dimensions of the convenience constructors?

[dlfivefifty]

What's the usage of the convenience constructors?

[marcusdavidwebb]

What separates a convenience constructor from a constructor? Marcus

…

On 31 May 2017, at 12:16, Sheehan Olver ***@***.***> wrote: What's the usage of the convenience constructors? — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#16 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ADbXaPhKnr1NRN_bqJY2csB9F3w_uejEks5r_T3mgaJpZM4NfY0a>.

[dlfivefifty]

I think what this package really needs is some inconvenience constructors 😜

I think he just means "convenience" as in they are not necessary for the package, but included for convenience.

[marcusdavidwebb]

Ah ok. I would have called them “developer constructors” as I imagine only developers of the software will actually use them. Marcus

…

On 31 May 2017, at 12:25, Sheehan Olver ***@***.***> wrote: I think what this package really needs is some inconvenience constructors 😜 I think he just means "convenience" as in they are not necessary for the package, but included for convenience. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#16 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ADbXaOkTTXlrv--ayZWH1wkOYZLxbiDHks5r_UAigaJpZM4NfY0a>.

[MikaelSlevinsky]

Well, no constructor is more inconvenient than a useless constructor! They stay as is for now.

[MikaelSlevinsky]

Another potential improvement to add to the list is to transpose the data. Givens rotations have better memory localization in row-major ordering, and the 1D row FFTs wouldn't need to transpose if they were 1D column FFTs.

[meggart]

Is it possible to allow Spherical harmonics analysis for arrays with an even number of columns? This currently kills my Julia session:

using FastTransforms
x=rand(100,100)

PA = FastTransforms.plan_analysis(x)
O=zeros(x)

A_mul_B!(O,PA,x)

while it works fine for an uneven number of columns, e.g. x=rand(100,99)

[MikaelSlevinsky]

In principle, yes, it just requires choosing a convention. The problem is equivalent to creating a Fourier interpolant through an even number of points. For three points, we return an expansion in the basis {1, sin(θ), cos(θ)}, but for four points should the expansion be in the basis {1, sin(θ), cos(θ), sin(2θ)}, {1, sin(θ), cos(θ), cos(2θ)} or a weighted average of the two? (In the code, the natural choice would be the first option.)

The current segfault behaviour is undeniably undesirable, apologies.

[meggart]

No problem. I think since there is no obvious reason for picking one or the other option, choosing the the one that is most natural to code should be ok? Thanks for implementing this in Julia, so far I had to call into SHTOOLS, and having a Julia-native solution helps a lot.

[MikaelSlevinsky]

The segfault issue is now fixed on master. The added test shows that using either odd or even equispaced longitudinal grids results in the resolution of a function on the sphere to approximately machine precision --- roughly the same spherical harmonic coefficients.

[meggart]

Thanks again, everything is working for me now.

[MikaelSlevinsky]

Super!

[meggart]

Sorry to bother again, but there seems to be a similar issue when doing the actual transform, which occurs only for larger arrays. Here is an MWE:

using FastTransforms
four2d = rand(720,1440)
coefs = FastTransforms.fourier2sph(four2d)
println("Works!")
four2d = rand(1440,2879)
coefs = FastTransforms.fourier2sph(four2d)
println("Works!")
four2d = rand(1440,2880)
coefs = FastTransforms.fourier2sph(four2d)
println("Segfault!")

The first two transforms work, while the third one does not, so to my guess is a combination of even number of columns + large array seems to cause this...

[MikaelSlevinsky]

This is now because you're using a larger number of columns than the dimension of the space of degree 1439 harmonics. I'll try to push a fix shortly, but for now this works:

julia> padfour2d = [four2d; zeros(1,2880)];

julia> coefs = fourier2sph(padfour2d; sketch = :none); # Pro-tip: the sketch keyword reduces the pre-computation by about half

[gradywright]

Amazing package @MikaelSlevinsky!

I'm wondering if it is possible to add a plan_fourier2sph function. Also it would be helpful if you had an example of how to use current plan_sph2fourier.

[MikaelSlevinsky]

Thanks! Actually, P = plan_sph2fourier(A) is all you need because P*A is equivalent to sph2fourier(A), but P\A is also the planned version of fourier2sph(A).

The best usage examples are in tests:

https://github.com/MikaelSlevinsky/FastTransforms.jl/blob/master/test/sphericalharmonics/synthesisanalysistests.jl

or there's a handful here:

https://github.com/JuliaApproximation/SpectralTimeStepping.jl/tree/master/NonlocalSphere

Those are all of the codes for the experiments and figures of https://arxiv.org/abs/1801.04902.

[gradywright]

Thanks for the quick reply and the links. New question:

F = sphrandn(Float64, 1024, 1024);
G = sph2fourier(F; sketch = :none);
P = plan_sph2fourier(F);
K = P*F;
norm(K-G)

Results in

8.017458048386342e-11

Why is the norm of the difference not zero?

[MikaelSlevinsky]

The sketch keyword gets passed directly to the interpolative decompositions in the butterfly algorithm. The planned functions can also take the same keywords. The default for sketch is a type of random column selection whereas :none uses RRQR to compute ID's. More info on the keywords for ID's can be found here https://github.com/JuliaMatrices/LowRankApprox.jl. For spherical harmonic transforms, I always use sketch = :none for high degrees: it gives better accuracy and it's roughly twice as fast. For accuracy, you can always compare with the SlowSphericalHarmonicPlan(A).

Note that the reason for the API inconsistency (leg2cheb and cheb2leg both have their own plans), is that the spherical harmonic connection problem is not 1-1. Rather than having errors in assuming that fourier2sph works like the inverse of a square matrix, I stubbornly want to enforce that it's a weighted least squares ;).

[MikaelSlevinsky]

The latest release (0.3.0) has multithreading support for the O(n³) algorithms, so the degree n = 1024 cutoff is somewhat obsolete. Just fire up Julia with export JULIA_NUM_THREADS=8 (or your maximum) in terminal.

[MikaelSlevinsky]

Looking at the tests, it appears that julia 0.7 has introduced a performance degradation that must be fixed. Compare v0.6 https://travis-ci.org/MikaelSlevinsky/FastTransforms.jl/jobs/414184117 with v0.7 https://travis-ci.org/MikaelSlevinsky/FastTransforms.jl/jobs/414184119 in particular the fast plan and thin plan tests. Also, I don't know why the in-place A_mul_B! (mul!) methods are allocating on the first and second calls, but if the thinplantests.jl are run again, then the allocations disappear.

[dlfivefifty]

It’s mostly a type inference change, so checking @code_warntype should put that out.

[MikaelSlevinsky]

This is a v0.7 MWE setup that's destroying the FMM-based transforms:

import Base: size, getindex, setindex!
import LinearAlgebra: rank
struct ThreadSafeVector{T} <: AbstractVector{T}
    V::Matrix{T}
    function ThreadSafeVector{T}(V::Matrix{T}) where T
        @assert size(V, 2) == Threads.nthreads()
        new{T}(V)
    end
end

ThreadSafeVector(V::Matrix) = ThreadSafeVector{eltype(V)}(V)
ThreadSafeVector(v::Vector) = ThreadSafeVector(repmat(v, 1, Threads.nthreads()))

@inline size(v::ThreadSafeVector) = (size(v.V, 1), )
@inline getindex(v::ThreadSafeVector{T}, i::Integer) where T = v.V[i, Threads.threadid()]
@inline setindex!(v::ThreadSafeVector{T}, x, i::Integer) where T = v.V[i, Threads.threadid()] = x

threadsafezeros(::Type{T}, n::Integer) where T = ThreadSafeVector{T}(zeros(T, n, Threads.nthreads()))
threadsafeones(::Type{T}, n::Integer) where T = ThreadSafeVector{T}(ones(T, n, Threads.nthreads()))

abstract type AbstractLowRankMatrix{T} <: AbstractMatrix{T} end

"""
Store the singular value decomposition of a matrix:

    A = UΣV'

"""
struct LowRankMatrix{T} <: AbstractLowRankMatrix{T}
    U::Matrix{T}
    Σ::Diagonal{T}
    V::Matrix{T}
    temp::ThreadSafeVector{T}
end

LowRankMatrix(U::Matrix{T}, Σ::Diagonal{T}, V::Matrix{T}) where T = LowRankMatrix(U, Σ, V, threadsafezeros(T, length(Σ.diag)))

size(L::LowRankMatrix) = size(L.U, 1), size(L.V, 1)
rank(L::LowRankMatrix) = length(L.Σ.diag)

function getindex(L::LowRankMatrix{T},i::Integer,j::Integer) where T
    ret = zero(T)
    U, Σ, V, r = L.U, L.Σ, L.V, rank(L)
    for k = r:-1:1
        ret += U[i,k]*Σ[k,k]*V[j,k]
    end

    ret
end


function test_mul!(y::AbstractVecOrMat{T}, L::LowRankMatrix{T}, x::AbstractVecOrMat{T}, istart::Int, jstart::Int, INCX::Int, INCY::Int) where T
    m, n = size(L)
    ishift, jshift = istart-INCY, jstart-INCX
    temp = L.temp

    @inbounds for k = 1:rank(L)
        temp[k] = zero(T)
        for j = 1:n
            temp[k] += L.V[j,k]*x[jshift+j*INCX]
        end
        temp[k] *= L.Σ[k,k]
    end

    @inbounds for k = 1:rank(L)
        tempk = temp[k]
        for i = 1:m
            y[ishift+i*INCY] += L.U[i,k]*tempk
        end
    end

    y
end

After running that, this shows a significant degradation (in v0.7)

x = randn(10000)
y = zero(x)
L = LowRankMatrix(rand(1000, 20), Diagonal(rand(20)), rand(1000, 20))

@time test_mul!(y, L, x, 1, 1, 2, 2);
@time test_mul!(y, L, x, 1, 1, 2, 2);

[MikaelSlevinsky]

Actually, it has nothing to do with the ThreadSafeVector type. This is all pure Julia code... surprising! Even smaller example:

import Base: size, getindex, setindex!
import LinearAlgebra: rank
abstract type AbstractLowRankMatrix{T} <: AbstractMatrix{T} end

"""
Store the singular value decomposition of a matrix:

    A = UΣV'

"""
struct LowRankMatrix{T} <: AbstractLowRankMatrix{T}
    U::Matrix{T}
    Σ::Diagonal{T}
    V::Matrix{T}
    temp::Vector{T}
end

LowRankMatrix(U::Matrix{T}, Σ::Diagonal{T}, V::Matrix{T}) where T = LowRankMatrix(U, Σ, V, zeros(T, length(Σ.diag)))

size(L::LowRankMatrix) = size(L.U, 1), size(L.V, 1)
rank(L::LowRankMatrix) = length(L.Σ.diag)

function getindex(L::LowRankMatrix{T},i::Integer,j::Integer) where T
    ret = zero(T)
    U, Σ, V, r = L.U, L.Σ, L.V, rank(L)
    for k = r:-1:1
        ret += U[i,k]*Σ[k,k]*V[j,k]
    end

    ret
end


function test_mul!(y::AbstractVecOrMat{T}, L::LowRankMatrix{T}, x::AbstractVecOrMat{T}, istart::Int, jstart::Int, INCX::Int, INCY::Int) where T
    m, n = size(L)
    ishift, jshift = istart-INCY, jstart-INCX
    temp = L.temp

    @inbounds for k = 1:rank(L)
        temp[k] = zero(T)
        for j = 1:n
            temp[k] += L.V[j,k]*x[jshift+j*INCX]
        end
        temp[k] *= L.Σ[k,k]
    end

    @inbounds for k = 1:rank(L)
        tempk = temp[k]
        for i = 1:m
            y[ishift+i*INCY] += L.U[i,k]*tempk
        end
    end

    y
end

x = randn(10000)
y = zero(x)
L = LowRankMatrix(rand(1000, 20), Diagonal(rand(20)), rand(1000, 20))

@time test_mul!(y, L, x, 1, 1, 2, 2);
@time test_mul!(y, L, x, 1, 1, 2, 2);

[MikaelSlevinsky]

Ah, perhaps it's because Diagonal is now typed by the underlying array.

[MikaelSlevinsky]

In less recent versions, there was only one A_mul_B! (before the name change). Now, HierarchicalMatrices and FastTransforms have their own mul! that is sometimes separate from LinearAlgebra.

This change sabotages some attempts to use mul! for allocation-free transforms since that calls the abstract dense fallbacks. Some transforms need an explicit FastTransforms.mul!, which is unacceptable: this package has every right to override mul! when using different types. @dlfivefifty

[dlfivefifty]

Not sure why there’s a separate mul! here. But mul! in Base is a bit broken, because it should use traits. I made a PR but we couldn’t agree to names, so that PR now lives in LazyArrays.jl

[MikaelSlevinsky]

I think part of the problem may be that mul! (A_mul_B!) used to reside in Base, but now it's in LinearAlgebra. Perhaps HierarchicalMatrices needs more than just the Compat patch, i.e. using LinearAlgebra.

[dlfivefifty]

It's more than that. Adjoint, Transpose, UpperTriangular, etc. overrides of mul! have a cascading effect where to add your own overrides requiring more and more ambiguity overrides. The solution to this is to use traits.

My feeling is forget about LinearAlgebra.mul! and switch over to LazyArrays, which also has a more descriptive syntax:

c .= α .* Mul(A,b) .+ β .* c

This will support all the Base types that mul! does and more (e.g., it knows that view(Adjoint(view(A,1:2,1:3)), 1:5, 2:3) is BLAS row major and lowers to the correct BLAS routines.).

The plan is that a trait-like mul! or LazyArray's approach will be incorporated into Base/LinearAlgebra in the future.

[AshtonSBradley]

how can I set up allocation free spherical harmonic transforms? I may have missed something, but I can't track down a version of

function *(P::SphericalHarmonicPlan, X::AbstractMatrix)
    mul!(zero(X), P, X)
end

that lets me pass my own field to write to instead of zero(X). Should I just be writing my own method, or is that somehow dangerous?

Also, is there any fix for the performance regression?

[MikaelSlevinsky]

As long as X and Y have the same sizes and P = plan_sph2fourier(X), then I think mul!(Y, P, X) is what you're looking for. (There is no support for lmul!(P, X) at the moment.)

[maximilian-gelbrecht]

The segfault issue is now fixed on master. The added test shows that using either odd or even equispaced longitudinal grids results in the resolution of a function on the sphere to approximately machine precision --- roughly the same spherical harmonic coefficients.

I was wondering, support for even numbered polar angles/longitudes does not seem to be implemented anymore. Could it be reimplemented? I'd argue that most grids used by models/data (e.g. in meteorology) that I have encountered, have an even number of angles/longitudes, so that it would be extremely convenient to have this working.

Thanks.

[MikaelSlevinsky]

Yes, that feature (annoying to the developer) didn't make it to the new world. MikaelSlevinsky/FastTransforms#23