abs fails on union domain and PiecewiseSegment domain

[macd]

if we have f = abs(sin(5x))^3 we get a union of domains. If we have g = abs(sin(5x)^3), we get a PiecewiseSegment domain. abs(f) fails because "ERROR: ArgumentError: Domain cannot be empty" and abs(g) also fails (which should be no-op here). What I really want is to be able to calculate the total variation of f''' (which I define as tv(f) = sum(abs(f')) because ApproxFun's norm(f, 1) is still broken. If curious, see p. 51 of ATAP). I can, however, get tv(g''') to work, but that gets the same wrong answer as ATAP 16527.836634664207 (but that was corrected in the ATAP errata). So that looks like another problem where ApproxFun is missing the contributions from the Dirac delta functions. Here is a repro for the first problem. For reference for others, as I am sure you know, ATAP is Trefethen's "Approximation Theory and Approximation Practice"

julia v1.2> versioninfo()
Julia Version 1.2.0-rc2.0
Commit 9248bf7687 (2019-07-08 19:42 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-6770HQ CPU @ 2.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, skylake)
Environment:
  JULIA_DEPOT_PATH = /home/macd/julia-versions/dot_julia-1.2

julia v1.2> using ApproxFun

julia v1.2> x = Fun()
Fun(Chebyshev(),[0.0, 1.0])

julia v1.2> f = abs(sin(5x))^3;

julia v1.2> domain(f)
a union of 4 domains:
    1.  : the segment [-1.0,-0.6283185307179585]
    2.  : the segment [-0.6283185307179585,-5.070867905569396e-16]
    3.  : the segment [-5.070867905569396e-16,0.6283185307179586]
    4.  : the segment [0.6283185307179586,1.0]

julia v1.2> fa = abs(f)
ERROR: ArgumentError: Domain cannot be empty
Stacktrace:
 [1] Type at /home/macd/julia-versions/dot_julia-1.2/packages/ApproxFunOrthogonalPolynomials/OM7o1/src/Spaces/Chebyshev/Chebyshev.jl:16 [inlined]
 [2] Type at /home/macd/julia-versions/dot_julia-1.2/packages/ApproxFunOrthogonalPolynomials/OM7o1/src/Spaces/Chebyshev/Chebyshev.jl:22 [inlined]
 [3] canonicalspace at /home/macd/julia-versions/dot_julia-1.2/packages/ApproxFunOrthogonalPolynomials/OM7o1/src/Spaces/Ultraspherical/DirichletSpace.jl:30 [inlined]
...

[dlfivefifty]

yes delta functions are barely supported. I’m in the process of redesigning the backend of ApproxFun in ContinuumArrays.jl, and one of the targets is hp-FEM which should give more robust support for piecewise smooth and continuous functions. But it will take at least a year before this is realised

[macd]

Fair enough. Do you have a high level road map in mind? I know it depends on funding, post-docs, grad students, changing research agendas, but it would be nice to know your thinking.

[dlfivefifty]

I think the goal is to make it work robustly for multidimensional PDEs. Theres proofs of concepts for triangles and disks in MultivariateOrthogonalPolynomials.jl. To make this high performance there needs to be support for

[dlfivefifty]

Sorry finger slip

[dlfivefifty]

Broad road map: make ApproxFun applicable multidimensional high performance settings. Some initial exploration on disks/triangles/spheres is in MultivariateOrthogonalPolynomials.jl, though substantial redesign is needed to support distributed arrays, etc.

To tackle this I'm trying to separate the linear algebra from the complexity of ApproxFun. This is part of the goal of LazyArrays.jl+InfiniteArrays.jl+InfiniteBandedMatrices.jl: this will allow the creation of operators and infinite linear algebra to know nothing about Spaces or Funs. Since the initial representation is lazy realisation on distributed memory should be really easy: there's minimal data communication needed.

ContinuumArrays.jl is another important piece of the puzzle in that it cleans up a lot of the language in ApproxFun. Here both bases and operators become (quasi-)arrays which makes many of the complicated types like SumSpace, SubSpace, SubOperator, InterlaceOperator etc. into simple array operations, many of which already exist. This would make a simpler way of producing an infinite array that pipes into (possibly distributed) InfiniteBandedMatrices.jl solvers. This package already supports concepts that are not currently possible in ApproxFun: for example, FEM via splines/weak formulation can be implemented in a few lines of array-like code, this solves Poisson:

    L = LinearSpline(range(0,stop=1,length=10)) # a quasi-matrix with axes(L) == (0..1, OneTo(10))
    B = L[:,2:end-1] # Zero dirichlet by dropping first and last spline
    D = Derivative(axes(L,1)) # Derivative on 0..1, just a placeholder
    Δ = -((B'D')*(D*B)) # Weak Laplacian. Returns a `SymTridiagonal`
    f = L*exp.(L.points) # project exp(x) onto linear splines
    u = B * (Δ \ (B'f)) # solve Poisson
    @test u[0.1] ≈ -0.06612902692412974

ApproxFun.jl will then become a lightweight wrapper around ContinuumArrays.jl that translates continuous language (f(x)) to array language (f[x]). Operators will become functions of functions instead of weird pseudo-array like objects.

[dlfivefifty]

But note that multiple roots is a bit "ill-posed". It shouldn't error out but I'm not entirely sure how to avoid this.

[macd]

Thanks for your thoughts on this. I appreciate it. Let me know if / how I can help.

[dlfivefifty]

If you are serious and want to help with the roots you could look into what Chebfun does in this situation.

[macd]

Sure, I can take a look. Is there a canonical test case?

[dlfivefifty]

I meant for this particular issue. That is, how does chebfun do f = abs(sin(5x))^3?

[macd]

OK, will do.,

[macd]

So, I did get back to this and while it is an interesting problem, I'm not certain what's to be done either.

So, f = abs(sin(5x))^3 should have only 3 roots, but it turns out that both Chebfun and ApproxFun find more because they are doing root finding on the approximation to f, which you are well aware of. Chebfun finds that there are 7 roots while ApproxFun finds 4 roots. Subsequent behavior is very different.

In ApproxFun, if we try to perform abs(f) we get the error that is reported in this issue. With Chebfun, we get a new approximation that now has 8 smooth pieces, rather than the original 4 smooth pieces that both Chebfun and ApproxFun have.

chebfun column (8 smooth pieces) interval length endpoint values
[ -1, -0.63] 20 0.88 6.1e-16 [ -0.63, -0.63] 1 2.8e-16 2.8e-16 [ -0.63, 0] 25 3e-16 -3.6e-17 [ 0, 2.1e-06] 4 3.6e-16 7.7e-16 [ 2.1e-06, 0.63] 25 9.4e-16 7.8e-16 [ 0.63, 0.63] 4 6.3e-16 7.2e-16 [ 0.63, 0.63] 3 9e-17 5.4e-16 [ 0.63, 1] 20 8.1e-16 0.88 vertical scale = 1 Total length = 102

And, in fact, we can continue to take the abs of the results, accumulating more smooth pieces as we go, so that k = abs(abs(abs(abs(abs(f))))) will yield

chebfun column (16 smooth pieces) interval length endpoint values
[ -1, -0.63] 20 0.88 6.4e-16 [ -0.63, -0.63] 1 2.8e-16 2.8e-16 [ -0.63, -0.63] 3 3.5e-16 1.5e-15 [ -0.63,-2.3e-06] 25 1.4e-15 1.6e-15 [-2.3e-06,-6.7e-07] 3 1.5e-15 0 [-6.7e-07,-2.8e-07] 1 8.6e-18 8.6e-18 [-2.8e-07, 0] 1 3.8e-17 3.8e-17 [ 0, 2.1e-06] 4 3.6e-16 7.7e-16 [ 2.1e-06, 0.63] 25 8e-16 5.7e-16 [ 0.63, 0.63] 4 6.3e-16 7.2e-16 [ 0.63, 0.63] 1 3.6e-17 3.6e-17 [ 0.63, 0.63] 1 8.6e-19 8.6e-19 [ 0.63, 0.63] 1 3e-17 3e-17 [ 0.63, 0.63] 2 5.5e-27 4.8e-16 [ 0.63, 0.63] 1 1e-15 1e-15 [ 0.63, 1] 20 7.1e-16 0.88 vertical scale = 1 Total length = 113

But you can see that the original pieces are now sandwiched between very small intervals with very small expansions. But, that appears to have no bad side effects, beyond getting many more roots (there are now 15 of them). That is, sum(f - k) is 4.115362652361836e-17

So it looks like the Chebfun is much more robust, at least if you are not dependent upon the exact number of roots

So reflecting on this and other behavior, here is a kind of wacky idea that goes against the philosophy of Chebfun / ApproxFun. That is to encode some of the properties that a specific Fun would have, such as even, odd, allpositive, etc. Then even Funs would have only coefficients for even Tn, etc, and abs(f) for an allpositive Fun would be a no-op. But, having said that, it is probably hard to infer these properties and maybe not that generally useful.

But, getting back to the issue at hand, what does it take to replicate the Chebfun behavior? (I tried to look, but lost patience when I kept stepping past the error using Debugger)

Oh, BTW, in Chebfun, if you take the roots of a function that has a jump discontinuity that crosses zero, it finds that root too. This is a feature that can be turned off, but it is on by default. ApproxFun does not do this and misses these "roots".

[dlfivefifty]

I like your wacky idea: this would be pretty naturally expressed as a space Chebyshev()[1:2:end]. We could in theory check if the even coefficients are numerically zero and proceed.

we don’t necessarily want to replicate Chebfun exactly: the underlying implementation is the same but some of the discretisation sizes for the colleague matrices differ as Julia uses faster Hessenberg eigs which allows for larger matrices.

It’s probably just a simple bug that’s creating empty intervals, I’ll look into it.

Getting roots with functions with jumps is more involved. This is where the hp-FEM would be really helpful: we’d want an analogue of Heaviside space for piecewise polynomials, which would result from differentiating a ContinuousSpace. (The current behaviour of producing a PiecewiseSpace doesn’t scale well.) Do it’s something that should wait until ContinuumArrays becomes the backend and all of the flakey choices made in ApproxFun can be fixed

[macd]

So it looks like a major redo of the architecture of ApproxFun (and a boatload of work!). I will probably continue to try some of the problems in ATAP. Should I still file issues against ApproxFun? Also, are the Townsend and Trefethen, and the Olver and Townsend the best references on quasi-arrays? (The hp-FEM stuff looks scary to me.)

[dlfivefifty]

I think file in ApproxFun: things like roots will likely live here (though there is a quasi-vector analogue findall(iszero, v)...)